# Stochastic Multi Configuration Time-Dependent Hartree for Dissipative
Quantum Dynamics with Strong Intramolecular Coupling

Souvik Mandal∗<EMAIL_ADDRESS>Laboratoire de Physique et Chimie
Théoriques, UMR 7019 CNRS/Université de Lorraine, 1 Bd Arago, 57070 Metz,
France  Fabien Gatti  Institut de Sciences Moléculaires d’Orsay, UMR 8214
CNRS/Université Paris-Saclay, Bât 520, rue André Rivière, 91405 Orsay CEDEX,
France  Oussama Bindech  Roberto Marquardt  Laboratoire de Chimie Quantique -
Institut de Chimie - UMR 7177 CNRS/Unistra
Université de Strasbourg
4, rue Blaise Pascal - CS 90032 - 67081 STRASBOURG CEDEX - France  Jean
Christophe Tremblay∗<EMAIL_ADDRESS>Laboratoire de
Physique et Chimie Théoriques, UMR 7019 CNRS/Université de Lorraine, 1 Bd
Arago, 57070 Metz, France

###### Abstract

In this article, we explore the dissipation dynamics of a strongly coupled
multidimensional system in contact with a Markovian bath following a system-
bath approach. We use in this endeavour the recently developed stochastic
Multi-Configuration Time-Dependent Hartree approach within the Monte Carlo
wave packet formalism [J. Chem. Phys. 156, 094109 (2022)]. The method proved
to yield thermalized ensembles of wave packets when intramolecular coupling is
weak. To treat strongly coupled systems, new Lindblad dissipative operators
are constructed as linear combinations of the system coordinates and
associated momenta. These are obtained by an unitary transformation to a
normal mode representation, which reduces intermode coupling up to second
order. Additionally, we use combinations of generalized raising/lowering
operators to enforce the Boltzmann distribution in the dissipation operators,
which yield perfect thermalization in the harmonic limit. The two ansatz are
tested using a model two-dimensional hamiltonian parameterized to disentangle
the effects of intramolecular potential coupling, of strong mode mixing
observed in Fermi resonances, and of anharmonicity.

##  I Introduction

The study of the quantum dynamics of multidimensional systems is an extremely
demanding task, both experimentally and theoretically. From an experimental
point of view, the isolation of features that are characteristic of the
quantum mechanics of the underlying processes is a crucial aspect of
measurements that often suffer from perturbing incoherent processes. The main
theoretical difficulty originates from the nature of the wavefunctions which
are typically spread over a large number of degrees of freedom and therefore
require a large numerical basis to be characterized. If the number of strongly
interacting particles is very large, the standard solutions of the quantum
mechanical equations of motion become prohibitively expensive with presently
available computational facilities because of the highly increased density of
states. However, often we are interested in understanding the dynamics of a
certain subsystem of strongly interacting particles, which interact weakly
with a large environment. This situation then becomes more manageable within a
so-called system-bath Breuer and Petruccione (2006); Nitzan (2006); Blum
(2012) approach, where the system is tackled full quantum mechanically, while
the bath acts as a reservoir for energy transfer.

In a previous work, Mandal  _et al._ (2022) we examined the possibility of
describing an open quantum system with energy dissipation and population
relaxation in connection to a thermostat within the Lindblad formalism. This
new method complements other work in high dimensional system-bath dynamics.
Gao (1997, 1998); Nest and Meyer (2003); Nest and Saalfrank (2000); Nest and
Meyer (2002); Pesce and Saalfrank (1997, 1998); Pesce  _et al._ (1998);
Saalfrank (2006); Tremblay, Beyvers, and Saalfrank (2008); Tremblay and
Saalfrank (2009); Tremblay, Monturet, and Saalfrank (2010, 2011); Tremblay
(2011); Tremblay, Füchsel, and Saalfrank (2012); Füchsel  _et al._ (2012);
Tremblay (2013); Serwatka and Tremblay (2019) In our approach, a Markovian
master equation Gao (1998); Ford and O’Connell (1999) is derived for a
Lindblad operator in sum-of-product form. The master equation is solved by
propagation of stochastic wave packets Breuer and Petruccione (2007) using the
Monte Carlo Wave Packet (MCWP) methodology,Dalibard, Castin, and Mølmer
(1992); Dum, Zoller, and Ritsch (1992); Carmichael (1993); Mølmer, Castin, and
Dalibard (1993) which is now implemented in the Heidelberg Multi-Configuration
Time-Dependent Hartree (MCTDH) code.Worth  _et al._ The method was benchmarked
for a two-dimensional vibrational model system. The population of different
states was found to behave thermally within their statistical fluctuations in
the weak intramolecular coupling regime. For strong coupling, a hyperthermal
behaviour was observed, which originates from both intramolecular coupling and
anharmonicity, and which lacks a sound physical interpretation.

In this paper, we present solutions to overcome the unphysical hyperthermal
behaviour. First, we introduce a linear coordinate transformation to reduce
the potential coupling upto second order, similarly to a normal mode picture.
Secondly, we introduce generalized raising/lowering operators to enforce
thermal behaviour more pertinently. The paper is organized as follows. In
Section II, the dissipative operators and their rotation to normal coordinates
are discussed. Section III describes the different model systems used to
benchmark the new approach. In Section IV, the thermalization dynamics of the
different model systems is studied. Section V presents our conclusions.

##  II Theory

We have recently implemented the Monte Carlo Wave Packet (MCWP) method
Dalibard, Castin, and Mølmer (1992); Dum, Zoller, and Ritsch (1992);
Carmichael (1993); Mølmer, Castin, and Dalibard (1993) in the Heidelberg
Multi-Configuration Time-Dependent Hartree (MCTDH) code Worth  _et al._ to
study dissipative quantum dynamics. It is a stochastic unraveling of the
reduced density matrix propagation, which we dubbed sMCTDH for “stochastic
MCTDH”. In a nutshell the quantum dynamics of a system in contact with a
thermal bath is represented by a set of solutions of a stochastic time
dependent Schrödinger equation. These solutions are called realizations. Each
stochastic wavepacket is propagated piecewise deterministicallyBreuer and
Petruccione (2006) under the operation of a nonlinear hamiltonian, which
contains a dissipative term. The environment induces fluctuations in the
system, which we account for by sampling the distribution of waiting times. On
average over many realizations of the stochastic process, we recover the
ensemble dynamics of the reduced density matrix. The details of the MCWP
method and its algorithm can be found elsewhere. Serwatka and Tremblay (2019);
Mandal  _et al._ (2022)

Here we only focus on the specific form of the nonlinear hamiltonian chosen
for the sMCTDH method

$\displaystyle{\hat{H}}$ $\displaystyle=$ $\displaystyle{\hat{H}}_{\rm
sys}-\frac{{\rm i}\hbar}{2}\sum\limits_{k=1}^{f}\sum\limits_{j=1}^{M_{\rm
op}}\hat{A}_{kj}^{\dagger}(T)\hat{A}_{kj}(T)$ (1)

where $\hat{H}_{\rm sys}$ is the system hamiltonian, the specific form of
which will be discussed in Section III. $\sum\limits_{j=1}^{M_{\rm
op}}\hat{A}_{kj}^{\dagger}(T)\hat{A}_{kj}(T)$ is a temperature-dependent
dissipation operator with bath temperature $T$, and $\hbar$ is the reduced
Planck constant. $f$ is the total number of degrees of freedom and $M_{\rm
op}$ is the number of operators for each mode $k$. This factorizable form of
the dissipation operator is also consistent with the form of the MCTDH
wavefunction. We chose $M_{\rm op}$ = 3 for each mode $k$ and the
corresponding three different dissipation channels can be defined for each
individual mode of the system of interest by the following operators Gao
(1997); Mandal  _et al._ (2022); Ford and O’Connell (1999)

$\hat{A}_{k\,1}(T)=\mu_{k}(T)\,\hat{w}_{k},\qquad\hat{A}_{k\,2}(T)=\nu_{k}(T)\,\hat{p}_{k},\qquad\hat{A}_{k\,3}=\kappa_{k}\,\left(\hat{w}_{k}+\frac{\rm
i}{m_{k}\,\omega_{k}}\,\hat{p}_{k}\right)\,.$ (2)

Here, $\hat{w}_{k}$ is the multiplicative operator of a generalized coordinate
for mode $k$ and $\hat{p}_{k}$ is the associated momentum operator. As
$\hat{w}$ is multiplicative operator, we may leave out the circumflex over the
symbol, when the symbol represents the coordinate. $\omega_{k}$ is the
harmonic angular frequency of the potential considered and $m_{k}$ is the
corresponding reduced mass. The first two operators $\hat{A}_{k1}(T)$ and
$\hat{A}_{k2}(T)$ induce fluctuations in the system. The operator
$\hat{A}_{k3}$ is the generalized lowering operator for mode $k$. With the
above definition, Eq. (2), the dynamics of a two-dimensional model system
consisting of two coupled vibrational modes was investigated in our previous
work. Mandal  _et al._ (2022) In the strong intermode coupling regime a
hyperthermal behaviour of the high-energy mode was observed at asymptotic
times. This hyperthermal behaviour was particularly strong for the case where
the high-energy mode relaxes more slowly than the low-energy mode. In our the
previous work Mandal  _et al._ (2022), we gave convincing evidence that the
hyperthermal behaviour arises from intermode coupling. In order to reduce the
effect of intermode coupling, we advocate here a linear transformation of the
generalized set of coordinates $w_{k}$ to a normal mode representation denoted
${Q}_{k}$ defined by the following reciprocal relations:

$\displaystyle Q_{k}=\sum_{l=1}^{f}c_{kl}\sqrt{m_{l}}w_{l}$ (3) $\displaystyle
w_{l}=\sum_{k=1}^{f}\frac{c_{lk}}{\sqrt{m_{l}}}Q_{k}$ (4)

For the sake of clarity, we call in this paper the coordinates $w_{l}$ local
mode coordinates. The coefficients $c_{kl}$ compose the orthonormal
eigenvectors of the hessian matrix

$\displaystyle
H_{f}=\begin{pmatrix}\frac{1}{m_{1}}\frac{\partial^{2}V}{\partial
w_{1}^{2}}&\frac{1}{\sqrt{m_{1}m_{2}}}\frac{\partial^{2}V}{\partial
w_{1}\partial
w_{2}}&\ldots&\frac{1}{\sqrt{m_{1}m_{f}}}\frac{\partial^{2}V}{\partial
w_{1}\partial w_{f}}\\\
\frac{1}{\sqrt{m_{2}m_{1}}}\frac{\partial^{2}V}{\partial w_{2}\partial
w_{1}}&\frac{1}{m_{2}}\frac{\partial^{2}V}{\partial
w_{2}^{2}}&\ldots&\frac{1}{\sqrt{m_{2}m_{f}}}\frac{\partial^{2}V}{\partial
w_{2}\partial w_{f}}\\\ \vdots&\vdots&\ddots&\vdots\\\
\frac{1}{\sqrt{m_{f}m_{1}}}\frac{\partial^{2}V}{\partial w_{f}\partial
w_{1}}&\frac{1}{\sqrt{m_{f}m_{2}}}\frac{\partial^{2}V}{\partial w_{f}\partial
w_{2}}&\ldots&\frac{1}{m_{f}}\frac{\partial^{2}V}{\partial
w_{f}^{2}}\end{pmatrix}$ (5)

where $V$ is the potential energy hypersurface. That allows us to decouple the
modes exactly upto second order in the potential expansion. All higher-order
couplings will still remain after the transformation. The angular frequency
$\omega_{k}^{\prime}$ for the corresponding normal mode is calculated from the
eigenvalue ($\lambda_{k}$) of $H_{f}$.

$\displaystyle\omega^{\prime}_{k}=\sqrt{\lambda_{k}}$ (6)

The conjugate momentum of $\hat{Q}_{k}$ can be defined as

$\hat{p}_{Q_{k}}=-{\rm i}\hbar\sum_{l=1}^{f}\frac{\partial w_{l}}{\partial
Q_{k}}\frac{\partial}{\partial w_{l}}$ (7)

The quantity $\frac{\partial w_{l}}{\partial
Q_{k}}=\frac{c_{lk}}{\sqrt{m_{l}}}$ is just a number, which can be
straightforwardly calculated from the linear transformation, Eq. (4). Using
the normal mode coordinate $Q_{k}$ with frequency $\omega^{\prime}_{k}$ and
the corresponding conjugate momentum $\hat{p}_{Q_{k}}$ as defined in Eq. (3),
(6) and (7), the dissipation operators can be defined for individual normal
modes following Eq. (2)

$\hat{A}_{k\,1}(T)=\mu_{k}(T)\,Q_{k},\qquad\hat{A}_{k\,2}(T)=\nu_{k}(T)\,\hat{p}_{Q_{k}},\qquad\hat{A}_{k\,3}=\kappa_{k}\,\left(Q_{k}+\frac{\rm
i}{\omega^{\prime}_{k}}\,\hat{p}_{Q_{k}}\right)$ (8)

Since the normal modes are linear combination of mass-weighted generalized
coordinates, the mass factor is included in the mass-weighted momentum
$\hat{p}_{Q_{k}}$. Substituting Eq. (3) and Eq. (7) in Eq. (8) yields a
redefinition of the dissipation operators in terms of local mode coordinates

$\displaystyle\hat{A}_{k\,1}(T)=\mu_{k}(T)\,\sum_{l=1}^{f}c_{kl}\sqrt{m_{l}}w_{l},\qquad\hat{A}_{k\,2}(T)=\,-{\rm
i}\hbar\nu_{k}(T)\,\sum_{l=1}^{f}\,\frac{c_{lk}}{\sqrt{m_{l}}}\frac{\partial}{\partial
w_{l}},$
$\displaystyle\hat{A}_{k\,3}=\kappa_{k}\,\left(\sum_{l=1}^{f}\left(c_{kl}\sqrt{m_{l}}w_{l}+\frac{{\hbar}}{\omega^{\prime}_{k}}\,\frac{c_{lk}}{\sqrt{m_{l}}}\frac{\partial}{\partial
w_{l}}\right)\right)$ (9)

We recognize that the dissipation operators in the normal mode representation
are simple linear combinations of the same operators in local mode
coordinates. As such, they retain a simple sum-of-product form compatible with
the sMCTDH algorithm. For each mode 3 different mechanisms for energy
dissipation can be identified. Altogether, there will be $3\times f$
dissipation channels that will give different significance to the overall,
effectively temperature-dependent dissipation rate. The coefficients
$\mu_{k}(T)$, $\nu_{k}(T)$, and $\kappa_{k}$ are chosen to enforce thermal
behaviour in an uncoupled one-dimensional harmonic mode Gao (1998); Ford and
O’Connell (1999), and they take the form

$\displaystyle\mu_{k}^{2}(T)$ $\displaystyle=$
$\displaystyle\frac{\gamma_{k}\omega^{\prime}_{k}}{2\hbar}\left(\coth\Big{(}\frac{\hbar\omega^{\prime}_{k}}{2k_{\rm
B}\,T}\Big{)}-1\right)$ (10) $\displaystyle\nu_{k}^{2}(T)$ $\displaystyle=$
$\displaystyle\frac{\gamma_{k}}{2\hbar\omega^{\prime}_{k}}\left(\coth\Big{(}\frac{\hbar\omega^{\prime}_{k}}{2k_{\rm
B}\,T}\Big{)}-1\right)$ (11) $\displaystyle\kappa_{k}^{2}$ $\displaystyle=$
$\displaystyle\frac{\gamma_{k}\omega^{\prime}_{k}}{2\hbar}$ (12)

where $k_{\rm B}$ is the Boltzmann constant and the quantity $\gamma_{k}$ is a
rate parameter. The frequency $\omega_{k}^{\prime}$ is obtained from Eq. (6)
and thus it contains the information on the intermode coupling upto second
order. For a derivation of Eqs. (10)-(12), the reader is referred to the
literature Mandal  _et al._ (2022).

In order to achieve thermal distribution exactly in the limit of the
bilinearly coupled harmonic oscillators, the temperature-dependent
coefficients $\mu_{k}(T)$ and $\nu_{k}(T)$ can be modified as

$\displaystyle\tilde{\mu}k(T)=\kappa_{k}e^{\hbar\omega_{k}^{\prime}/2k_{\rm
B}T}$ (13) $\displaystyle\tilde{\nu}_{k}(T)=-\frac{\rm
i}{\omega_{k}^{\prime}}\kappa_{k}e^{\hbar\omega_{k}^{\prime}/2{k_{\rm B}T}}$
(14)

With this new definition of the coefficients, the dissipation operators in Eq.
(8) can be rewritten as

$\hat{A}_{k\,1}(T)+\hat{A}_{k\,2}(T)=\kappa_{k}\,e^{\hbar\omega_{k}^{\prime}/2k_{\rm
B}T}\left(\hat{Q}_{k}-\frac{\rm
i}{\omega^{\prime}_{k}}\,\hat{p}_{Q_{k}}\right)\qquad\hat{A}_{k\,3}=\kappa_{k}\,\left(\hat{Q}_{k}+\frac{\rm
i}{\omega^{\prime}_{k}}\,\hat{p}_{Q_{k}}\right)$ (15)

It can be recognized that the sum of the first two dissipation operators,
$\hat{A}_{k\,1}(T)+\hat{A}_{k\,2}$(T), becomes a generalized raising operator.
The rates associated with the upward and downward transitions are now related
via the constant $\kappa_{k}$, and they are constrained to impose thermal
detailed balance exactly. Note that the factor of 2 in the exponent is
necessary because the rates appear squared in the Lindblad equation. The
physical interpretation of the new set of operators as temperature-dependent
fluctuations and dissipation terms is blurred, with the expected advantage of
reaching thermal equilibrium at asymptotic times in exactly separable cases.

##  III Model

To study the dissipation dynamics and its thermalization behaviour, we chose a
two-dimensional (2D) model consisting of two coupled oscillators as in our
previous work Mandal  _et al._ (2022). The model is parameterized such as to
represent a dioxygen molecule adsorbed on a Pt(111) surface Gland, Sexton, and
Fisher (1980); Gao (1998), where the O-O and the Pt-O2 bonds are the two
coupled oscillators. The O-O vibration ($R$-mode) will be considered to be a
harmonic mode, and the O2-Pt vibration ($Z$-mode) will be either a harmonic
mode or an anharmonic mode. The $R$ and $Z$-modes are represented by the
coordinates $y_{R}$ and $y_{Z}$, respectively. We investigate three different
types of model potentials for this O2-Pt system below.

###  III.1 Anharmonic strong coupling potential

We call anharmonic strong coupling potential a potential in which the $Z$-mode
is anharmonic, the $R$-mode is harmonic and the modes are bilinearly coupled.
The hamiltonian of such a system can be written as

$\displaystyle\hat{H}_{\rm
sys}=\left(-\frac{\hbar^{2}}{2m_{R}}\frac{\partial^{2}}{\partial
R^{2}}\right)+\left(-\frac{\hbar^{2}}{2m_{Z}}\frac{\partial^{2}}{\partial
Z^{2}}\right)+D_{\rm
e}\bigg{(}a_{Z}^{2}y_{Z}^{2}\,+\,a_{R}^{2}\left(y_{R}\;+\;Cy_{Z}\right)^{2}\bigg{)}$
(16)

where the last term is the potential $V(R,Z)$. $y_{Z}$ is a Morse coordinate
and $y_{R}$ is a harmonic coordinate:

$\displaystyle y_{Z}=\frac{1-{\rm e}^{-a_{Z}\,(Z-Z_{\rm e})}}{a_{Z}},\qquad
y_{R}=R-R_{\rm e}$ (17)

The quantities $a_{Z}$ and $a_{R}$ are defined as

$\displaystyle a_{Z}=\sqrt{\frac{m_{Z}\,\omega_{Z}^{2}}{2\,D_{\rm
e}}}\;\;\;\;\;\;{\rm
and}\;\;\;\;\;\;a_{R}=\sqrt{\frac{m_{R}\,\omega_{R}^{2}}{2\,D_{\rm e}}}$ (18)

where $D_{\rm e}$ is the depth of the potential well. Here, $R$ is the
distance between the two oxygen atoms, with $R_{\rm e}$ being its value at
equilibrium. $Z$ is the distance from the center of mass of the oxygen
molecule to a platinum atom, with $Z_{\rm e}$ being the equilibrium distance.
$\omega_{Z/R}$ and $m_{Z/R}$ are the corresponding frequency and mass of the
respective modes. The quantity $C$ is a dimensionless coupling constant. As we
focus on the thermalization behaviour in the strong coupling regime the value
of the constant $C$ is chosen as 0.5. Mandal  _et al._ (2022) The reader is
referred to ref. 4 for results with moderate and weaker couplings. For this
anharmonic strong coupling potential, the local mode operators used to define
the dissipative operators in Eq. (II) with $k=\\{1,2\\}$ take the form

$w_{1}\equiv
w_{Z}=y_{Z}-\left(\frac{\hbar}{4m_{Z}\omega_{Z}\Lambda_{Z}}\right)^{1/2};\qquad
w_{2}\equiv w_{R}=y_{R}\\\ $ (19)

The quantity $\Lambda_{k}$ is defined by $\Lambda_{k}=2D_{\rm
e}\,/\,(\hbar\,\omega_{k})$. The eigenstates of $\hat{H}_{\rm sys}$ in Eq.
(16) will be denoted by the ket $\lvert n_{v_{Z}\,v_{R}}\rangle$, where
$n=0,1,2,\ldots$ is a cardinal number specifying the energetic order of the
states, while $v_{Z}$ and $v_{R}$ signify the quantum numbers of those zeroth-
order ($C$ = 0 case) states which have the strongest participation in the
given eigenstate. The set of the first twelve zeroth-order states and
eigenstates is given in Table 1. They are computed using the Heidelbebrg MCTDH
package. A basis of 5 single particle functions is used for each mode. The
integration is performed on a primitive grid composed of 64 discrete variable
functions (sine-DVR) for each one of the two modes.

State | Relative energies in units of $hc$ cm-1
---|---
| zeroth-order | eigenstates
| $C=0$ | $C=0.5$
$\lvert 0_{00}\rangle$ | 000.0 | 000.0
$\lvert 1_{10}\rangle$ | 399.5 | 384.4
$\lvert 2_{20}\rangle$ | 770.7 | 739.6
$\lvert 3_{01}\rangle$ | 870.0 | 910.5
$\lvert 4_{30}\rangle$ | 1113.5 | 1071.9
$\lvert 5_{11}\rangle$ | 1269.5 | 1289.6
$\lvert 6_{40}\rangle$ | 1427.0 | 1385.4
$\lvert 7_{21}\rangle$ | 1640.7 | 1617.1
$\lvert 8_{50}\rangle$ | 1713.9 | 1690.3
$\lvert 9_{02}\rangle$ | 1740.0 | 1819.2
$\lvert 10^{\ast}\rangle$ | 1971.5 | 1894.8
$\lvert 11^{\ast}\rangle$ | 1983.5 | 1995.2
Table 1:  Ordering and energies of zeroth-order states and eigenstates of the
anharmonic strong coupling potential system (see Eq. (16)). All reported
energies are relative to the respective ground vibrational state in units of
$hc$ cm-1. The states 10 and 11 with the ∗ symbol constitute a Fermi resonance
pair between the two zeroth-order states $\lvert v_{Z}=6,v_{R}=0\rangle$ and
$\lvert v_{Z}=3,v_{R}=1\rangle$. The parameters used in Eq. (16) are $m_{Z}$ =
27.48 u, $m_{R}$ = 7.9995 u, $\hbar\omega_{Z}$ = 53 meV, $\hbar\omega_{R}$ =
108 meV, $Z_{e}$ = 2.11 Å, $R_{e}$ = 1.37 Åand $D_{\rm e}$ = 0.4 eV. The
parameters are taken from refs. 6, 30, 31 and 32.

###  III.2 Fermi resonance potential

Relative energies in units of $hc$ cm-1
---
zeroth-order |  | eigenstates
State | $C=0$ | State | $C=0.1$ | $C=0.5$
$\lvert{0,0}\rangle$ | 000.0 | $\lvert 0_{1}\rangle$ | 000.0 | 000.0
$\lvert{1,0}\rangle$ | 400.0 | $\lvert(1/2)_{1}\rangle$ | 400.0 | 400.0
$\lvert{2,0}\rangle$ | 800.0 | $\lvert 1_{2}\rangle$ | 799.2 | 796.0
$\lvert{0,1}\rangle$ | 800.0 | $\lvert 1_{1}\rangle$ | 800.8 | 804.0
$\lvert{3,0}\rangle$ | 1200.0 | $\lvert(3/2)_{2}\rangle$ | 1198.6 | 1193.1
$\lvert{1,1}\rangle$ | 1200.0 | $\lvert(3/2)_{1}\rangle$ | 1201.4 | 1207.0
$\lvert{4,0}\rangle$ | 1600.0 | $\lvert 2_{3}\rangle$ | 1597.7 | 1588.8
$\lvert{2,1}\rangle$ | 1600.0 | $\lvert 2_{2}\rangle$ | 1600.0 | 1600.0
$\lvert{0,2}\rangle$ | 1600.0 | $\lvert 2_{1}\rangle$ | 1602.3 | 1611.4
$\lvert{5,0}\rangle$ | 2000.0 | $\lvert(5/2)_{3}\rangle$ | 1996.8 | 1984.1
$\lvert{3,1}\rangle$ | 2000.0 | $\lvert(5/2)_{2}\rangle$ | 2000.0 | 2000.1
$\lvert{1,2}\rangle$ | 2000.0 | $\lvert(5/2)_{1}\rangle$ | 2003.2 | 2016.2
Table 2:  Ordering and energies of zeroth-order states and eigenstates of the
harmonic Fermi resonance potential system (see Eq. (20)). All reported
energies are relative to the respective ground vibrational state in units of
$hc$ cm-1. The parameters used in Eq. (20) are $m_{Z}$ = 27.48 u, $m_{R}$ =
7.9995 u, $\hbar\omega_{Z}$ = 49.59 meV and $\hbar\omega_{R}$ = 99.18 meV.

Here we investigate another type of strong coupling, i.e., due to a Fermi
resonance. To this end, we use coupled harmonic oscillators both for the $R$\-
and $Z$-modes with coordinates $y_{R}=(R-R_{\rm e})$ and $y_{Z}=(Z-Z_{\rm
e})$, respectively. The $Z$-mode is quadratically coupled with the $R$-mode by
a dimensionless coupling constant $C$. The hamiltonian of the system can be
written as

$\displaystyle\hat{H}_{\rm
sys}=\left(-\frac{\hbar^{2}}{2m_{R}}\frac{\partial^{2}}{\partial
R^{2}}\right)+\left(-\frac{\hbar^{2}}{2m_{Z}}\frac{\partial^{2}}{\partial
Z^{2}}\right)+\left(\frac{1}{2}\,m_{Z}\omega_{Z}^{2}y_{Z}^{2}+\frac{1}{2}m_{R}\omega_{R}^{2}\left(y_{R}+\frac{C}{Z_{\rm
e}}y_{Z}^{2}\right)^{2}\right)$ (20)

where the last term is the potential $V(R,Z)$. Here we set $\omega_{R}$ =
2$\omega_{Z}$ with the remaining parameters as in the previous example. Since
there is a 2:1 ratio between the frequencies of the two modes, due to the
coupling, the creation/annihilation of one vibrational quanta in the $R$-mode
will indirectly infer the creation/annihilation of two vibrational quanta in
the $Z$-mode. To leading order, the coupling is linear in the $R$-mode, but
quadratic in the $Z$-mode. Because of the energetic resonance the two modes
are strongly coupled even for small value of $C$. Table 2 shows the energies
of the first twelve zeroth-order states and eigenstates ($C$ = 0.1 and 0.5) of
$\hat{H}_{\rm sys}$ defined in Eq. (20). The zeroth-order states are defined
by $\lvert v_{Z},v_{R}\rangle$, where $v_{Z}$ and $v_{R}$ are the quantum
numbers for the $Z$\- and $R$-mode respectively. The eigenstates are labeled
using the polyad notation, where the $j^{\rm th}$ state of the $N$ polyad is
denoted as $\mid N_{j}\rangle$. Here, $N=\max({v_{R}})+\max({v_{Z}})/2$, where
$\max({v_{k}})$ refers to the maximal quantum number of the zeroth-order
state, $j=1,\ldots,N+1$, if $N$ is integer, or $j=1,\ldots,N+1/2$, if $N$ is
half-integer, and $j=1$ designates the eigenstate with the highest energy
within the polyad. This notation was used in refs. 33; 34. For the Fermi
resonance, the local mode operators in Eq. (II) take the form (with
$k\in\\{1,2\\}$)

$\displaystyle\hat{w}_{1}\equiv\hat{w}_{Z}=y_{Z};\qquad\hat{w}_{2}\equiv\hat{w}_{R}=y_{R}$
(21)

###  III.3 Bilinearly coupled harmonic potential

Finally, we explore the bilinear coupled harmonic potential, where both $R$\-
and $Z$-modes are harmonic coordinates $y_{R}=(R-R_{\rm e})$ and
$y_{Z}=(Z-Z_{\rm e})$, respectively. The model hamiltonian follows from that
of the strongly coupled anharmonic case in Eq. (16) in the limit $D_{\rm
e}\to\infty$:

$\displaystyle\hat{H}_{\rm
sys}=\left(-\frac{\hbar^{2}}{2m_{R}}\frac{\partial^{2}}{\partial
R^{2}}\right)+\left(-\frac{\hbar^{2}}{2m_{Z}}\frac{\partial^{2}}{\partial
Z^{2}}\right)+\left(\frac{1}{2}\,m_{Z}\omega_{Z}^{2}y_{Z}^{2}+\frac{1}{2}m_{R}\omega_{R}^{2}\left(y_{R}+Cy_{Z}\right)^{2}\right)$
(22)

where the last term is the potential $V(R,Z)$. The $R$\- and $Z$-modes are
bilinearly coupled by the dimensionless coupling constant $C$. The bilinearly
coupled harmonic system is useful because it can be perfectly decoupled using
the normal mode approach.

State | Relative energies in units of $hc$ cm-1
---|---
| zeroth-order | eigenstates
| $C=0$ | $C=0.5$
$\lvert 0_{00}\rangle$ | 000.0 | 000.0
$\lvert 1_{10}\rangle$ | 395.1 | 378.4
$\lvert 2_{20}\rangle$ | 790.1 | 756.8
$\lvert 3_{01}\rangle$ | 869.9 | 908.2
$\lvert 4_{30}\rangle$ | 1185.2 | 1135.3
$\lvert 5_{11}\rangle$ | 1265.1 | 1286.6
$\lvert 6_{40}\rangle$ | 1580.2 | 1513.7
$\lvert 7_{21}\rangle$ | 1660.1 | 1665.1
$\lvert 8_{02}\rangle$ | 1739.9 | 1816.4
Table 3:  Ordering and energies of zeroth-order states and eigenstates of the
bilinearly coupled harmonic potential system (see Eq. (22)). All reported
energies are relative to the respective vibrational ground state in units of
$hc$ cm-1. The parameters used in Eq. (22) are $m_{Z}$ = 27.48 u, $m_{R}$ =
7.9995 u, $\hbar\omega_{Z}$ = 53 meV and $\hbar\omega_{R}$ = 108 meV.

Therefore, it reaches well-defined limits that can be used to benchmark the
dynamics of the system. The value of $C$ is 0.5, as above, for the
investigation in the strong coupling regime. The eigenstates of $\hat{H}_{\rm
sys}$ in Eq. (22) are denoted by $\lvert n_{v_{Z}\,v_{R}}\rangle$ as for the
case of an anharmonic bilinear strong coupling potential. The energy values of
the first nine zeroth-order states and eigenstates are presented in Table 3.
In this case, the local mode operators take the same form as in Eq. (21).

##  IV Results and Discussions

We perform sMCTDH simulations with the model potentials described above. The
population dynamics and their thermalization behaviours at different bath
temperatures are investigated in the following subsections.

###  IV.1 Anharmonic strong coupling potential

Figure 1 presents the population dynamics of the eigensates for the anharmonic
strong coupling potential ($C$ = 0.5) described by the system hamiltonian in
Eq. (16) at a bath temperature of 200 K. The initial wavefunction is chosen as
a product of a Gaussian distribution shifted in $Z$ and a Gaussian
distribution in $R$ centered at $R_{\rm e}$. We perform the simulation using
two different types of dissipation coordinates. Figure 1(a) displays the
population dynamics using the dissipation operators as described in Eq. (2),
where the dissipation operators are represented in terms of the local mode
coordinates ${w}_{k}$ defined in Eq. (19). On the other hand, Figure 1(b)
shows the population dynamics using the dissipation operators as described in
Eq. (II).

As the initial wavefunction is a displaced wave packet along the $Z$
coordinate the initial populations are mostly dominated by the states excited
in the $Z$-mode (e.g. states $\lvert 10_{60}\rangle$, $\lvert 8_{50}\rangle$
and $\lvert 6_{40}\rangle$). A few states excited along both the $Z$\- and
$R$-modes (e.g. states $\lvert 11_{31}\rangle$ and $\lvert 7_{21}\rangle$) are
also contributing to the initial population due to strong intermode coupling.
With time the population of the excited states transfers to the lower states,
mostly via the $Z$-channel because the dissipation rate along the $Z$-mode is
faster ($\gamma_{Z}^{-1}$ = 0.5 ps) compared to that along the $R$-mode
($\gamma_{R}^{-1}$ = 2 ps ). Finally, the populations go to the ground state,
which takes the energy from the bath and re-excites state $\lvert
3_{01}\rangle$ along the $R$-mode. The early time dynamics for both the local
mode and normal coordinates are mostly similar in nature. Their behaviour is
different at asymptotic times, when the equilibrium is reached. As our
interest is on the thermalization of the system, this long-time behaviour will
be the focus of the analysis.

When dissipation is mediated by operators in the local mode representation,
Figure 1(a) reveals that the first excited states along each mode reach about
the same population. This behaviour is in principle inline with a Boltzmann
distribution since state $\lvert 3_{01}\rangle$ has a higher energy than state
$|1_{10}\rangle$. The overtone of the $Z$-mode, state $\lvert 2_{20}\rangle$,
is found at an energy similar to state $\lvert 3_{01}\rangle$, with a
separation of about 170 cm-1 (see Table 1). Consequently, their character
become highly mixed by the strong potential coupling. This apparently leads to
an unphysically high population of the $R$-mode (state $\lvert 3_{01}\rangle$)
at asymptotic times, while also suppressing the population of the overtone,
state $\lvert 2_{20}\rangle$.

Figure 1: Evolution of population of selected eigenstates within the two-
dimensional model for O2/Pt(111) with the anharmonic strong coupling potential
(see Eq. (16)) at 200 K bath temperature for (a) local mode approach (Eq. (2))
and (b) normal mode approach (Eq. (II)). The initial wavefunction is chosen as
a product of a shifted Gaussian in $Z$ (center = 2.40 Å, width = 0.039 Å) and
a Gaussian in $R$ (center = 1.37 Å, width = 0.049 Å). The relaxation parameter
are $\gamma_{Z}^{-1}$ = 0.5 ps and $\gamma_{R}^{-1}$ = 2 ps and the coupling
strength is $C$ = 0.5. The solid black line indicates the sum of populations
of the states which are presented.

This problem is resolved by using the normal mode approach for the
dissipation, i.e. when using Eq. (II) instead of Eq. (2), as shown in Figure
1(a). The asymptotic population of state $\lvert 3_{01}\rangle$ which
corresponds to equilibrium is found to be strongly reduced. Interestingly, the
population of state $\lvert 1_{10}\rangle$ is almost unchanged when moving
from the local mode (Eq. (2)) to the normal mode representation, Eq. (II), of
the dissipation operators. This could be explained by the fact that state
$\lvert 1_{10}\rangle$ ($Z$-mode) exhibits less mixing character even for
strong intermode coupling. In Table 1, the energetic effect of this strong
coupling is nonetheless clearly observed. De facto, at 200 K, the system at
equilibrium can be reduced to a two-level system. The temperature associated
with this equilibrium can be estimated from the population of the ground state
and of state $\lvert 1_{10}\rangle$ at asymptotic times by inverting the
Boltzmann distribution and using the energies in Table 1. From Figure 1(a) we
find (202 $\pm$ 4) K by averaging the populations over the last two
picoseconds, which matches well with the simulation temperature (200 K).

The population dynamics using the normal modes at a bath temperature of 400 K
is shown in Figure 2 with the same initial condition as in Figure 1. The
populations for different states reach almost constant values at asymptotic
times. Here the population of states $\lvert 2_{20}\rangle$ and $\lvert
3_{01}\rangle$ are almost the same. Yet, state $\lvert 2_{20}\rangle$ has
lower energy than state $\lvert 3_{01}\rangle$ (see Table 1), such that the
latter should have a population at equilibrium that is smaller by a factor of
about 1.8 from that of the former if the system is to obey thermal detailed
balance. This is an indication that, despite the improvement in the
thermalization behaviour resulting from the normal mode representation of the
dissipation operators, population dissipation from the states $\lvert
2_{20}\rangle$ and $\lvert 3_{01}\rangle$ cannot be fully decoupled. As a
consequence, state $\lvert 2_{20}\rangle$ and $\lvert 3_{01}\rangle$
apparently do not lead to a proper thermal behaviour at higher temperature.

Figure 2: Evolution of population of selected eigenstates basis within the
two-dimensional model for O2/Pt(111) with anharmonic strong coupling potential
(see Eq. (16)) using normal mode approach (Eq. (II)) at 400 K bath
temperature. The initial conditions, the relaxation parameters and the
coupling constant are same as Figure 1.

In order to better understand the thermalized nature of the equilibrated
dynamics from the above populations, $\ln(P_{n}/P_{0}$) is calculated as a
function of the energy difference to the ground state, $E_{n}-E_{0}$. For this
study, dissipation dynamics is calculated with Eq. (II), exclusively. The
results are shown in Figures 3 (a) and (b) for 200 and 400 K, respectively. To
evaluate the convergence of the dynamics with respect to the number of
realizations of the stochastic process, the population $P_{n}$ of the $n^{\rm
th}$ eigenstate is first averaged over a particular time interval at
asymptotic times. Figure 3 depicts the linear regression of $\ln(P_{n}/P_{0}$)
as a function of $E_{n}-E_{0}$ obtained for four different time intervals
([10,11], [11,12], [12,13] and [13,14] ps). All eigenstates with populations
below 10-6 are neglected for the fitting. While the average slopes are drawn
as solid lines, the maximum and minimum range among the statistical errors
obtained as the uncertainty on the slopes of all fits is shown by a shaded
grey cone. The dashed red line is the theoretically calculated slope
($1/k_{\rm B}T$) from the Boltzmann distribution at 200 K (panel (a)) and 400
K (panel (b)). In both cases, the fits are well grouped and they reveal
statistically significant trends. Figure 3(a) shows that all fitted curves
including their statistical error cone pass above the theoretical expected
dashed red line. The system thus remains hyperthermal at 200 K for the
anharmonic strong coupling potential despite the normal mode approach of Eq.
(II). This hyperthermal behaviour arises because the dissipation operators are
only decoupled to second order. On the other hand, a slightly hypothermal
behaviour is seen in Figure 3(b) at a bath temperature of 400 K. This
different behaviour compared to a system at 200 K is possibly due to the
choice of the temperature-dependence for the coefficients $\mu(T)$ and
$\nu(T)$ (see Eq. (10) and (11)), which is physically well-motivated Gao
(1998) but not strictly valid at all temperatures. Ford and O’Connell (1999)
We shall review this point later on in this work. In general, although the
general trends in population for strongly coupled systems are in much better
agreement with the ensemble temperatures, the detailed behaviour remains
somewhat difficult to predict.

Figure 3: Fit (solid line) of $\ln(P_{n}/P_{0}$) as a function of
$E_{n}-E_{0}$ at bath temperatures (a) 200 K and (b) 400 K evaluated from
Figures 1(b) and 2 at different time interval as indicated by the colour code
key (valid for both panels). The gray shaded cone shows the maximum and
minimum range among the errors related to all fits. The dashed red line is the
theoretically expected behavior from the Boltzmann distribution at the
respective bath temperatures.

###  IV.2 Fermi resonance potential

Figure 4 displays the population dynamics for the harmonic Fermi resonance
potential (see Eq. (20)) with a bath temperature of 200 K for weak intermode
coupling ($C$ = 0.1). The left and right panels show the population dynamics
of the eigenstates (see Figure 4(a)) and zeroth-order states (see Figure 4(b))
of the system hamiltonian (see Eq. (20)), respectively. The dissipation
operator is represented by Eq. (II) in both cases. The initial wavefunction
for both cases is chosen as eigenstate $\lvert(5/2)_{2}\rangle$ of the
dissipation-free hamiltonian $\hat{H}_{\rm sys}$ (see Eq. (20)).

Figure 4: Evolution of population of selected (a) eigenstates and (b) zeroth-
order states for a two-dimensional harmonic Fermi resonance potential (see Eq.
(20)) coupled to a bath at 200 K. The initial wavefunction is chosen as an
eigenstate $\lvert(5/2)_{2}\rangle$ of the system hamiltonian (see Eq. (20)).
The relaxation parameters are $\gamma_{Z}^{-1}$ = 0.5 ps and $\gamma_{R}^{-1}$
= 2 ps and the coupling strength is $C$ = 0.1. Within a same polyad some
eigenstates have almost equal population throughout the simulation period.
Those states are {$\lvert 1_{2}\rangle$, $\lvert 1_{1}\rangle$},
{$\lvert(3/2)_{2}\rangle$, $\lvert(3/2)_{1}\rangle$}, {$\lvert 2_{3}\rangle$,
$\lvert 2_{1}\rangle$}, and {$\lvert(5/2)_{3}\rangle$,
$\lvert(5/2)_{1}\rangle$}.

As a result the initial population is found completely in eigenstate
$\lvert(5/2)_{2}\rangle$ when reporting the population dynamics in the
eigenstates basis (see Figure 4(a)). On the other hand, states $\lvert
1,2\rangle$ (major with 63% population) and $\lvert 5,0\rangle$ (minor with
37% population) dominate the initial population in the zeroth-order basis (see
Figure 4(b)). This shows that the initial state $\lvert(5/2)_{2}\rangle$ is a
strong mixture of $\lvert 1,2\rangle$ and $\lvert 5,0\rangle$, with almost no
contribution of $\lvert 3,1\rangle$. At early times (before 2ps), the dynamics
in the eigenstate basis follows two paths. One path is the sequential decay
from one polyad to the next. The population dynamics (see Fig. 4(a)) suggests
that this sequential path alternates between fast and slow population
transfers. We attribute this behaviour to the strong character mixing in the
Fermi resonance. At each step, the relaxation is dominated either by its fast
component ($Z$-mode, $\gamma_{Z}^{-1}=500$ fs) or its slow component
($R$-mode, $\gamma_{Z}^{-1}=2$ ps).

The second path involves redistribution of energy within the $(5/2)$ polyad
via a slow transition. Given that these states are eigenstates of the system,
this slow transition must involve the bath. It is followed by a fast
transition to the $\\{|2_{3}\rangle,|2_{1}\rangle\\}$ states, and a slow
transition to the $(3/2)$ polyad. These two paths can be depicted as (scheme
I)

$\lvert(5/2)_{2}\rangle$$\lvert
2_{2}\rangle$$\lvert(3/2)_{2}\rangle$$\lvert(3/2)_{1}\rangle$$\lvert
1_{2}\rangle$$\lvert 1_{1}\rangle$$\lvert(1/2)_{1}\rangle$$\lvert
0_{1}\rangle$$\lvert(5/2)_{3}\rangle$$\lvert(5/2)_{1}\rangle$$\lvert
2_{3}\rangle$$\lvert 2_{1}\rangle$FSFSFSFSF: FastS: Slow

The energy relaxation behaviour can be alternatively understood with the help
of zeroth-order states. Two relaxation paths can be identified from the two
dominantly populated zeroth-order states (scheme II)

minor (37%)$\lvert 5,0\rangle$$\lvert 4,0\rangle$$\lvert 3,0\rangle$$\lvert
2,0\rangle$$\lvert 1,0\rangle$F (0.1 ps)F (0.125 ps)F (0.167 ps)F (0.25
ps)major (63%)$\lvert 1,2\rangle$$\lvert 0,2\rangle$$\lvert 1,1\rangle$$\lvert
0,1\rangle$$\lvert 0,0\rangle$F (0.5 ps)S (1 ps)S (1 ps)F (0.5 ps)S (2 ps)F
(0.5 ps)IC (5.2 ps)IC (7.2 ps)IC (11.9 ps)IC (20.8 ps)F: FastS: SlowIC:
Internal Conversion

The timescales given in parenthesis are estimated from the rate parameters
$\gamma_{Z/R}$ in Eq. (10) using the harmonic scaling law. The internal
conversion (IC) timescales refer to population transfer within a polyad and
are calculated from the energetic splitting within that polyad. The early-time
dynamics in Figure 4(b) reveals that the minor component, state $\lvert
5,0\rangle$, which is a highly excited state along the $Z$-mode relaxes
sequentially and rapidly since the associate rate is faster. The major channel
starts with the fast $\lvert 1,2\rangle\to\lvert 0,2\rangle$ transition
mediated by the loss of one quantum along the $Z$-mode due to relaxation via
the bath. Internal conversion between states $\lvert 0,2\rangle$ and $\lvert
4,0\rangle$ is possible because the two states are part of the same polyad.
Consequently, the relaxation would continue along the fast sequential path
along the $Z$-mode (upper row in scheme II). The alternative path involves
relaxation of a quantum in the $R$-mode to reach the $\lvert 0,1\rangle$ state
(lower row in scheme II). It can also be reached via the $\lvert 1,1\rangle$
state in two steps, first by loosing a quantum in the slow $R$-mode and then
in the fast $Z$-mode. The population dynamics reveals that these two
sequential paths with alternating slow-fast relaxation towards state $\lvert
0,1\rangle$ are the dominant ones.

In the zeroth-order representation, population accumulates in the lowest-lying
excited state along the $R$-mode. This is because relaxation along $R$ is
slower. In the eigenstate basis, population accumulates in the polyad, $1$. It
is the lowest-lying polyad that contains strong character of the $R$-mode. The
energy of the $Z$-mode is half that of the $R$-mode. The state
$\lvert(1/2)_{1}\rangle$ is dominated by the $Z$-character, which undergoes
faster relaxation than the $R$-mode. Thus, population transfer from the $1$
polyad to state $|(1/2)_{1}\rangle$ is the rate limiting step.

By construction, populations for both eigenstates and zeroth-order states (see
Figure 4(a) and 4(b)) reach constant equilibrium values at asymptotic times.
The coefficients Eqs. (10) to (12) used in the dissipation operators are
chosen to aim for a thermal equilibrium, at least in the case of uncoupled
harmonic oscillators. The decoupling approach proposed in this work yields
thermal equilibrium also in the case of resonant coupling. To verify this, the
population $P_{N_{j}}$ of the eigenstates $\lvert N_{j}\rangle$ in Figure 4(a)
is used to fit $\ln(P_{N_{j}}/P_{0_{1}})$ as a function of energy $E_{N_{j}}$
-$E_{0_{1}}$. Figure 5 depicts the result of the linear fits (shown by the
colour coded solid lines) when the asymptotic populations are averaged over
time intervals [20,21] and [21,22] ps. The gray shaded cone shows the maximum
and minimum range among the associated statistical error related to the
fitting. The colour coded lines essentially overlap, which indicates that the
calculations are statistically converged. The theoretically expected dashed
red line from a Boltzmann distribution at 200 K fits very well within the
error of the fits. A similarly good result can be obtained for strong
intermode coupling ($C=0.5$). This implies that the Fermi resonance model
system can be decoupled and thermalized very effectively using the approach
based on normal modes in the Lindblad operator.

Figure 5: Fit (solid line) of $\ln(P_{N_{j}}/P_{0_{1}})$ as a function of
$E_{N_{j}}$ -$E_{0_{1}}$ for a bath temperature of 200 K evaluated from Figure
4(a) at two different time intervals as indicated by the colour code. The gray
shaded cone shows the maximum and minimum range among the errors related to
all fits. The dashed red line is the theoretically expected behavior from the
Boltzmann distribution at the bath temperature.

Despite the weak intermode potential coupling in the example discussed, the
presence of the Fermi resonance leads to a strong degree of mode mixing in the
model system. This implies that both degrees of freedom should reach the same
temperature, while this is not necessarily the case for coupled modes in
general (see next section and Appendix). Figure 6 depicts the time evolution
of temperatures calculated from the reduced populations of zeroth-order states
along each mode. These quantities are calculated from the dynamics reported in
Figure 4(b) as partial traces:

$\displaystyle P_{v_{Z}}^{({\rm r})}(t)$ $\displaystyle=$
$\displaystyle\sum_{v_{R}}P_{v_{Z}\;v_{R}}(t)$ (23) $\displaystyle
P_{v_{R}}^{({\rm r})}(t)$ $\displaystyle=$
$\displaystyle\sum_{v_{Z}}P_{v_{Z}\;v_{R}}(t)\,,$ (24)

where $P_{v_{Z}\;v_{R}}(t)$ is the population for state $\lvert
v_{Z},v_{R}\rangle$ at time $t$. The state temperature $T_{v_{Z/R}}(t)$ is
defined by inverting the Boltzmann distribution

$\displaystyle T_{v_{Z/R}}(t)=\frac{\epsilon_{v_{Z/R}}-\epsilon_{0}}{k_{\rm
B}}\;\;\Bigg{\\{}\ln\Bigg{(}\frac{P_{0}^{({\rm r})}(t)}{P_{v_{Z/R}}^{({\rm
r})}(t)}\Bigg{)}\Bigg{\\}}^{-1}\,.$ (25)

Here, $\epsilon_{v_{Z/R}}$ is the energy of the $v_{Z/R}$th vibrational state
of mode $Z/R$, respectively, defined from the 1D potentials given in Eq. (20)
by putting $C=0$. The temperatures calculated from Eq. (25) for the reduced
populations of states $v_{Z}=\\{1,2,3,4,5\\}$ (blue lines) and
$v_{R}=\\{1,2,3\\}$ (red lines) are presented as a function of time in Figure
6 (a) for $C=0.1$ and in Figure 6 (b) for $C=0.5$.

Figure 6: Temperature $T_{v_{Z/R}}$ for different $v_{Z}$ and $v_{R}$
evaluated using Eq. (25) for Fermi resonance potential as a function of time
for a bath at a temperature of 200 K. Panel (a) weak potential coupling, $C$ =
0.1. Panel (b) strong potential coupling, $C$ = 0.5.

All states in both the $Z$\- and $R$-modes start with different initial
temperatures – here depicted from 2 ps onward – which are generally high
compared to the bath temperature (200 K). Eventually, they all reach the bath
temperature. The $Z$-mode has the fast relaxation rate in the model and all
associated excited states rapidly reach the same temperature, following the
same type of exponential decay. The slower $R$-mode behaves differently. Since
relaxation from a higher excited state is faster, the $v_{R}=3$ population
reaches first the bath temperature, followed by the $v_{R}=2$ population. The
temperature associated with the $v_{R}=1$ state remains hyperthermal for a
longer time before also reaching the bath temperature. The fact that all
reduced populations of a given mode reach the same temperature is an
indication that the system is found at thermal equilibrium (see also Figure
5).

Despite the weak intermode coupling stemming from the potential, the character
of the individual modes is completely lost in the eigenstates because of the
exact resonant mixing of the zeroth-order states, independently of the size of
the coupling strength as shown by the result for stronger potential coupling
($C$ = 0.5), in Figure 6 (b). Here again, all temperatures stabilise at the
bath temperature 200 K, but the equilibration occurs even more rapidly due to
stronger intermode coupling. Interestingly, the stronger intermode coupling
leads to higher initial temperatures in the $Z$-mode while the $R$-mode is
found at lower temperature compared to the weak coupling case. Stronger
intermode coupling also induces faster relaxation and faster intramolecular
vibrational energy redistribution. While the Fermi resonance condition implies
that coupled modes loose their character in the composition of the system
eigenstates, the individual modes reach the same thermal equilibrium state as
in the later case.

The fact that the temperature of the two coupled modes is the same in Figure 6
is due to the Fermi resonance. It is related to the equivalence of time and
state average (ergodic theorem) as a result of an entirely coherent quantum
dynamics, was shown previously at the example of the wavepacket dynamics of
the infrared CH chromophore in $\rm CHX_{3}$ compounds Marquardt  _et al._
(1986). The Fermi resonance is mediated by the combination of two factors.
First, the zeroth order frequencies for the $Z$\- and $R$-mode are
commensurable, i.e., $\omega_{R}/\omega_{Z}$ = 2:1 in the present case.
Second, the specific form of the coupling chosen in the model hamiltonian is
$\sim y_{R}y_{Z}^{2}$ (see Eq. (20)). Such an operator can be expressed in
terms of creation and annihilation operators, where one quantum of vibration
in the $R$-mode is created while two quanta of vibration in the $Z$-mode are
destroyed, and vise versa. These two conditions being met, strong intermode
mixing becomes possible.

###  IV.3 Bilinearly coupled harmonic potentials

Figure 7: Evolution of population of selected eigenstates within the two-
dimensional model for O2/Pt(111) at (a) 200 K and (b) 400 K for the bilinearly
coupled harmonic potentials (see Eq. (22)). The initial wavefunction is chosen
as a product of a shifted Gaussian in $Z$ (center = 2.20 Å, width = 0.039 Å)
and a Gaussian in $R$ (center = 1.37 Å, width = 0.049 Å). The relaxation
parameters are $\gamma_{Z}^{-1}$ = 0.5 ps and $\gamma_{R}^{-1}$ = 2 ps and the
coupling strength is $C$ = 0.5.

Figure 7 displays the population dynamics of eigenstates for the bilinearly
coupled harmonic system (see Eq. (22)) with strong coupling ($C$ = 0.5) using
dissipation operators described by Eq. (II). The initial wavefunction is taken
as a product of a shifted Gaussian in $Z$ and a Gaussian in $R$ centered at
the equilibrium distance $R_{\rm e}$. The simulations are carried out at bath
temperatures of 200 K and 400 K, and shown in panels (a) and (b),
respectively. The initial state predominantly populates the ground state
($\lvert 0_{00}\rangle$), as well as excited states along the $Z$-mode
($\lvert 1_{10}\rangle$, $\lvert 2_{20}\rangle$ and $\lvert 4_{30}\rangle$)
due to initial displacement in $Z$. A few excited states in the $R$-mode
($\lvert 3_{01}\rangle$, $\lvert 5_{11}\rangle$ and $\lvert 7_{21}\rangle$)
are also weakly populated due to the strong intermode coupling. The figure
shows the usual sequential exchanges of population among different states and
modes, before finally reaching an equilibrium at asymptotic times. The
asymptotic populations of the excited states are higher for 400 K compared to
200 K, and the qualitative aspect of the equilibrium for the dominantly
populated states is physically correct.

Figure 8: Fit (solid line) of $\ln(P_{n}/P_{0}$) as a function of
$E_{n}-E_{0}$ for bath temperatres (a) 200 K and (b) 400 K evaluated from
Figure 7 at two different time intervals as indicated by the colour code. The
gray shaded cone shows the maximum and minimum range among the errors related
to all fits. The dashed red line is the theoretically expected behavior from
the Boltzmann distribution at the respective bath temperatures.

In order to understand the nature of the equilibrium, the populations of the
eigenstates from Figure 7 are used to plot $\ln(P_{\rm n}/P_{0})$ as a
function of energy $E_{\rm n}$ -$E_{\rm 0}$ using the method already discussed
in Figures 3 and 5. The results are plotted in Figure 8(a) and (b) for 200 K
and 400 K, respectively. Although the intermode coupling is strong the system
reaches thermal equilibrium at 200 K within the error bars. Despite the
improvement proposed by decoupling via normal mode analysis Eq. (II), the
population dynamics exhibits a hypothermal behaviour for a bath at a
temperature of 400 K. Since the normal mode transformation should be exact for
a bilinearly coupled harmonic system, this result is cause of concerns. The
error introduced in the model at higher temperature likely stems from the
definition of the coefficients used in the Lindblad operators. It is thus
interesting to investigate the behaviour of the other type of generalized
raising/lowering operators, defined by Eq. (15).

Figure 9: Evolution of population of selected eigenstates basis within the
two-dimensional model for O2/Pt(111) at (a) 200 K and (b) 400 K for the
bilinearly coupled harmonic potential (see Eq. (22)) and modified dissipation
operators based on Eqs. (13) to (15). The initial conditions, relaxation
parameters and the coupling constant are same as Figuer 7.

The population dynamics of the eigenstates for the bilinearly coupled harmonic
system is displayed in Figure 9. The initial conditions are the same as
described above. Panels (a) and (b) show the simulations at a bath temperature
of 200 K and 400 K, respectively. At 200 K, the population dynamics is
quantitatively the same as the one shown in Figure 7(a) using the original
definition of the constants $\mu_{k}(T)$ and $\nu_{k}(T)$. The population of
state $\lvert 3_{01}\rangle$ is visibly reduced at 400 K compared to Figure 7
(b), with all remaining states bearing little to no population. To provide a
more quantitative picture of the system temperature, averaged populations are
fitted according to the Boltzmann law as in the previous case studies. Figure
10, panels (a) and (b) display the fit of $\ln(P_{n}/P_{0}$) as a function of
$E_{n}-E_{0}$ at 200 K and 400 K, respectively. For the fit, the individual
state populations are averaged over four different time intervals: [10,11],
[11,12], [12,13] and [13,14] ps. All eigenstates with population below 10-6
are neglected for the fit. The maximum and minimum range among the errors
related to the fits (solid line) is shown by a shaded gray cone. The dashed
red line is the theoretical line calculated from the Boltzmann distribution
for 200 K and 400 K, respectively. In both cases, the theoretical dashed red
line matches perfectly within the fitting error.

Figure 10: Fit (solid line) of $\ln(P_{n}/P_{0}$) as a function of
$E_{n}-E_{0}$ evaluated for the bilinearly coupled harmonic potentials as
evaluated from Figure 9(a) and 9 (b) at bath temperatures 200 K (a) and 400 K
(b) and different time intervals as indicated by the colour code key (valid
for both panels). The gray shaded cone shows the maximum and minimum range
among the errors related to all fits. The dashed red line is the theoretically
expected behavior from the Boltzmann distribution at the respective bath
temperatures. The insets show the equivalent results for the anharmonic strong
coupling potential (see Eq. (16)).

This confirms that the ensemble behaves as a thermal distribution at both the
lower and higher temperature investigated here. At the lower temperature, the
population of the higher excited states is much smaller and the slope is
steeper. Further, the numerical uncertainty in their determination is larger,
such that the statistical error on the slope is larger than at the higher
temperature. Finally, the earlier time interval, [10,11] ps (green line in
Figure (10) (a)), yields a higher temperature than that of the bath, which is
indicative that the system is not completely equilibrated at those early time
intervals. Toward the end of the simulation (magenta and dark blue lines), the
fitted temperature is slightly lower than the bath temperature. It is worth
studying the effect of the new dissipation operators on the previously
investigated system of coupled anharmonic potentials and to compare the
results with those obtained in section IV.1, in particular in Figure 3. The
plots obtained for the anharmonic strong coupling potential with the
dissipation operators defined in Eq. ((15) are displayed in the insets of
Figure (10). The hyperthermal behaviour is still recognized at 200 K, but it
is much smaller at 400 K. We attribute this to a more significant population
of the higher excited states, which leads to smaller numerical errors in the
populations and in the fits.

The hyperthermal behaviour arises because the dissipation operators defined in
Eq. (15) do not depend on energy. The operators lead to the
creation/annihilation of a well defined energy quantum. While such operators
fit well the case of harmonic oscillators, they fail in the case of anharmonic
potentials. Morse-type anharmonic oscillators for instance have levels with
decreasing energetic spacing. Consequently, the raising/lowering dissipation
operators fulfill the detailed balance condition only approximately at higher
energies. This is likely the origin of the skewed population distribution.

In contrast to the Fermi resonance case discussed in the section IV.2, the two
modes are not found at the same temperature for the bilinearly coupled
harmonic potential even when the system has reached thermal equilibrium.
Figure 11 depicts the time dependent mode temperatures calculated from Eq.
(25) for the reduced populations of the lowest-lying states along the $Z$-mode
($v_{Z}=\\{1,2,3,4,5\\}$, blue lines) and along the $R$-mode
($v_{R}=\\{1,2,3\\}$, red lines) similarly to Figure 6. While the temperatures
for a given mode converge to the same value, the temperatures of the two modes
are different. This shows that the modes belong to different thermal
ensembles: the temperature of the low-frequency $Z$-mode is close to the bath
temperature of 200 K, but the high-frequency $R$-mode has a temperature of
about 370 K.

A rationale for the asymptotic stabilisation of the modes at a temperature
different from the bath and from the other modes is given in the Appendix,
where it is shown that a thermalized ensemble of multi-dimensional eigenstates
does not necessarily imply thermalization of individual modes of the zeroth-
order basis. This finding is in stark contrast with the Fermi resonance, where
strong mode mixing and thermalization at a single temperature was observed
even for weak intermode coupling. In the present case, the strong intermode
coupling does not lead to strong mode mixing. Rather, rapid coherent
population transfer between the two modes takes place at early times, as
observed from the fast oscillations of the individual state temperatures in
Figure 11. At longer times, these coherent oscillations disappear and
differential thermalization occurs. This is indicative of a very different
physical origin of the intermode coupling. The Fermi resonance redistributes
population among essentially isoenergetic states of different modes.

Figure 11: Temperature $T_{v_{Z/R}}$ for different $v_{Z}$ and $v_{R}$
evaluated using Eq. (25) for the bilinearly coupled harmonic potential using
modified dissipation operators based on Eqs. (13) to (15) as a function of
time for a bath at a temperature of 200 K.

##  V Conclusion

In conclusion, we have used the Monte Carlo wave packet method implemented in
the Heidelberg MCTDH package to study dissipative dynamics in 2D model systems
subject to strong intermode coupling. In our previous work, Mandal  _et al._
(2022) we have shown that under certain circumstances the dissipation
operators constructed from the system coordinates can lead to a hyperthermal
behaviour of modes due to intermode coupling. To alleviate this issue, we have
investigated the potential of a linear transformation of the dissipation
operators from the generalized coordinates describing the system to their
normal mode representation, which reduces intermode coupling up to the second
order. We have first studied the equilibration dynamics for three different
model potentials in the strong coupling regime. For the anharmonic strong
coupling potential, although the use of normal modes improves the thermal
behaviour at the asymptotic times compared to generalized coordinates, the
overall trend from the Boltzmann distribution of states remains hyperthermal
for a bath at 200 K temperature. This is because the transformation to the
normal modes can not fully eliminate the intermode coupling due to the
anharmonicity. On the other hand, the population dynamics at 400 K shows a
hypothermal behaviour. This is partly due to the temperature-dependent
constants related to dissipation operators which do not strictly enforce
thermal detailed balance at all temperatures.

In the case of the coupled harmonic Fermi resonance potential, the system
achieves the expected thermal distribution. The Fermi resonance also mixes
both modes strongly, even for weak coupling, and the individual mode
characteristics are completely lost. As a consequence, both modes are observed
to reach the same temperature at equilibrium.

Finally, we have studied the thermalization dynamics for a bimodal harmonic
potential with strong bilinear coupling, where the modes are perfectly
decoupled in the normal mode representation. Although the population dynamics
at 200 K reaches the expected thermal distribution, a hypothermal behaviour is
observed at 400 K, which stems in part from the definition of the temperature-
dependent constants in the dissipation operator.

An alternative form of generalized raising/lowering operators allowed us to
enforce exactly the correct thermal distribution in the case of the bilinearly
coupled harmonic potential. Despite the strong intramolecular coupling for the
harmonic case, the two modes reach different temperatures at equilibrium
because of the specific form of the intermode coupling. This is in contrast
with the harmonic Fermi resonance case, where strong mode mixing dominates the
thermalization dynamics of the individual modes due to exact commensurable
frequencies and the specific form of the leading coupling term. These new
dissipation operators also perform well for anharmonic potentials, in
particular at higher bath temperatures where higher excited states are more
significantly populated. The small hyperthermal behaviour for the anharmonic
case is related to the form of dissipation operators adopted in this work,
which are not energy dependent. Finding energy dependent dissipation operators
to get exact thermalization even for anharmonic cases is of future research
interest.

###### Acknowledgements.

We acknowledge financial support from ANR through the QDDA project 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

  * Breuer and Petruccione (2006) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2006).
  * Nitzan (2006) A. Nitzan,  _Chemical dynamics in condensed phases: relaxation, transfer and reactions in condensed molecular systems_ (Oxford university press, 2006).
  * Blum (2012) K. Blum,  _Density matrix theory and applications_ , Vol. 64 (Springer Science & Business Media, 2012).
  * Mandal  _et al._ (2022) S. Mandal, F. Gatti, O. Bindech, R. Marquardt, and J.-C. Tremblay, “Multidimensional stochastic dissipative quantum dynamics using a Lindblad operator,” J. Chem. Phys. 156, 094109 (2022).
  * Gao (1997) S. Gao, “Dissipative Quantum Dynamics with a Lindblad Functional,” Phys. Rev. Lett. 79, 3101–3104 (1997).
  * Gao (1998) S. Gao, “Lindblad approach to quantum dynamics of open systems,” Phys. Rev. B 57, 4509–4517 (1998).
  * Nest and Meyer (2003) M. Nest and H.-D. Meyer, “Dissipative quantum dynamics of anharmonic oscillators with the Multi-Configuration Time-Dependent Hartree (MCTDH) Method,” J. Chem. Phys. 119, 24 (2003).
  * Nest and Saalfrank (2000) M. Nest and P. Saalfrank, “Open-system quantum dynamics for gas-surface scattering: Nonlinear dissipation and mapped fourier grid methods,” J. Chem. Phys. 113, 8753–8761 (2000).
  * Nest and Meyer (2002) M. Nest and H.-D. Meyer, “Improving the mapping mechanism of the mapped Fourier method,” Chem. Phys. Lett. 352, 486–490 (2002).
  * Pesce and Saalfrank (1997) L. Pesce and P. Saalfrank, “‘Free’ nuclear density propagation in two dimensions the coupled-channel density matrix method and its application to inelastic molecule-surface scattering,” Chem. Phys. 219, 43–55 (1997).
  * Pesce and Saalfrank (1998) L. Pesce and P. Saalfrank, “The coupled channel density matrix method for open quantum systems: Formulation and application to the vibrational relaxation of molecules scattering from nonrigid surfaces,” J. Chem. Phys. 108, 3045–3056 (1998).
  * Pesce  _et al._ (1998) L. Pesce, T. Gerdts, U. Manthe, and P. Saalfrank, “Variational wave packet method for dissipative photodesorption problems,” Chem. Phys. Lett. 288, 383–390 (1998).
  * Saalfrank (2006) P. Saalfrank, “Quantum dynamical approach to ultrafast molecular desorption from surfaces,” Chem. Rev. 106, 4116 (2006).
  * Tremblay, Beyvers, and Saalfrank (2008) J. C. Tremblay, S. Beyvers, and P. Saalfrank, “Selective excitation of coupled co vibrations on a dissipative Cu (100) surface by shaped infrared laser pulses,” J. Chem. Phys. 128, 194709 (2008).
  * Tremblay and Saalfrank (2009) J. C. Tremblay and P. Saalfrank, “Selective subsurface absorption of hydrogen in palladium using laser distillation,” J. Chem. Phys. 131, 084716 (2009).
  * Tremblay, Monturet, and Saalfrank (2010) J. C. Tremblay, S. Monturet, and P. Saalfrank, “Electronic damping of adsorbate vibrations at metallic surfaces,” Phys. Rev. B 81, 125408 (2010).
  * Tremblay, Monturet, and Saalfrank (2011) J. C. Tremblay, S. Monturet, and P. Saalfrank, “The effects of electron–hole pair coupling on the infrared laser-controlled vibrational excitation of NO on Au(111),” J. Phys. Chem. A 115, 10698–10707 (2011).
  * Tremblay (2011) J. C. Tremblay, “Laser control of molecular excitations in stochastic dissipative media,” J. Chem. Phys. 134, 174111 (2011).
  * Tremblay, Füchsel, and Saalfrank (2012) J. C. Tremblay, G. Füchsel, and P. Saalfrank, “Excitation, relaxation, and quantum diffusion of CO on copper,” Phys. Rev. B 86, 045438 (2012).
  * Füchsel  _et al._ (2012) G. Füchsel, J. C. Tremblay, T. Klamroth, and P. Saalfrank, “Selective Excitation of Molecule-Surface Vibrations in H2 and D2 Dissociatively Adsorbed on Ru (0001),” Isr. J. Chem. 52, 438–451 (2012).
  * Tremblay (2013) J. C. Tremblay, “A unifying model for non-adiabatic coupling at metallic surfaces beyond the local harmonic approximation: From vibrational relaxation to scanning tunneling microscopy,” J. Chem. Phys. 138, 244106 (2013).
  * Serwatka and Tremblay (2019) T. Serwatka and J. C. Tremblay, “Stochastic wave packet approach to nonadiabatic scattering of diatomic molecules from metals,” J. Chem. Phys. 150, 184105 (2019).
  * Ford and O’Connell (1999) G. W. Ford and R. F. O’Connell, “Comment on “Dissipative Quantum Dynamics with a Lindblad Functional”,” Phys. Rev. Lett. 82, 3376–3376 (1999).
  * Breuer and Petruccione (2007) H. P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, UK, 2007).
  * Dalibard, Castin, and Mølmer (1992) J. Dalibard, Y. Castin, and K. Mølmer, “Wave-function approach to dissipative processes in quantum optics,” Phys. Rev. Lett. 68, 580 (1992).
  * Dum, Zoller, and Ritsch (1992) R. Dum, P. Zoller, and H. Ritsch, “Monte carlo simulation of the atomic master equation for spontaneous emission,” Phys. Rev. A 45, 4879 (1992).
  * Carmichael (1993) H. Carmichael,  _An Open Systems Approach to Quantum Optics_ , edited by H. Araki (Springer-Verlag, Berlin, 1993).
  * Mølmer, Castin, and Dalibard (1993) K. Mølmer, Y. Castin, and J. Dalibard, “Monte Carlo wave-function method in quantum optics,” J. Opt. Soc. Am. B 10, 524–538 (1993).
  * (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). Current version: 8.4.18 (2019). http://mctdh.uni-hd.de.
  * Gland, Sexton, and Fisher (1980) J. L. Gland, B. A. Sexton, and G. B. Fisher, “Oxygen interactions with the Pt(111) surface,” Surface Science 95, 587–602 (1980).
  * Wurth  _et al._ (1990) W. Wurth, J. Stöhr, P. Feulner, X. Pan, K. R. Bauchspiess, Y. Baba, E. Hudel, G. Rocker, and D. Menzel, “Bonding, structure, and magnetism of physisorbed and chemisorbed O2 on Pt(111),” Phys. Rev. Lett. 65, 2426–2429 (1990).
  * Steininger, Lehwald, and Ibach (1982) H. Steininger, S. Lehwald, and H. Ibach, “Adsorption of oxygen on Pt(111),” Surface Science 123, 1–17 (1982).
  * Dübal and Quack (1984) H.-R. Dübal and M. Quack, “Tridiagonal Fermi resonance structure in the IR spectrum of the excited CH chromophore in CF3H,” J. Chem. Phys. 81, 3779–3791 (1984).
  * Marquardt  _et al._ (1986) R. Marquardt, M. Quack, J. Stohner, and E. Sutcliffe, “Quantum-mechanical Wavepacket Dynamics of the CH Group in the Symmetric Top $\rm X_{3}CH$ Compounds using Effective Hamiltonians from High-resolution Spectroscopy,” J. Chem. Soc., Faraday Trans. 2 82, 1173–1187 (1986).

## Appendix:

The aim of this derivation is to demonstrate that a thermalized ensemble of
multi-dimensional eigenstate does not necessarily imply thermalization of the
individual modes at the same temperature. Our starting point is an incoherent
ensemble of eigenstates ($\lvert i\rangle$) described by the density matrix

$\displaystyle\rho=\sum_{i}p_{i}\lvert i\rangle\langle i\rvert$ (A1)

with populations distributed according to the Boltzmann law

$\displaystyle p_{i}=\frac{e^{-E_{i}/k_{\rm B}T}}{Z}$ (A2)

Here, $E_{i}$ is the energy of state $i$ and $Z$ is the partition function. In
general, we can write

$\displaystyle H\,\lvert i\rangle=E_{i}\,\lvert i\rangle$ (A3)

As a minimal model, we will study a system with two degrees of freedom. The
hamiltonian thus takes the form

$\displaystyle H=H_{1}+H_{2}+V$ (A4)

where $H_{1}$ and $H_{2}$ are the zeroth-order hamiltonian of modes 1 and 2,
respectively, and $V$ is the potential coupling between the two modes. The
individual mode hamiltonians can be used to defined a zeroth-order basis from
the stationary Schrödinger equation

$\displaystyle
H_{1}\lvert\alpha\rangle=\epsilon_{\alpha}^{(1)}\lvert\alpha\rangle~{}~{};~{}~{}H_{2}\lvert\beta\rangle=\epsilon_{\beta}^{(2)}\lvert\beta\rangle$
(A5)

where $\alpha$ and $\beta$ are zeroth order states with energy
$\epsilon_{\alpha}^{(1)}$ and $\epsilon_{\beta}^{(2)}$, respectively. The
coupled eigenstates $\lvert i\rangle$ of the full hamiltonian can be written
as linear combinations of Hartree products

$\displaystyle\lvert
i\rangle=\sum_{\alpha\beta}c_{i,\alpha\beta}\,\lvert\alpha\rangle\otimes\lvert\beta\rangle$
(A6)

This structure can be exploited to evaluate the reduced populations of each
mode by taking the partial trace of the density matrix $\rho$ over the
remaining degree of freedom (DOF).

### Case 1

The simplest case is that of uncoupled modes, which implies that the
eigenstates of the total hamiltonian can be written as a single Hartree
product with an energy defined from the single mode energies

$\displaystyle\lvert
i\rangle=\lvert\alpha\rangle\otimes\lvert\beta\rangle~{}~{};~{}~{}E_{i}=\epsilon_{\alpha}^{(1)}+\epsilon_{\beta}^{(2)}$
(A7)

In such a case the partition function takes a factorizable form

$\displaystyle
Z=\sum_{\alpha\beta}e^{-(\epsilon_{\alpha}^{(1)}+\epsilon_{\beta}^{(2)})/k_{\rm
B}T}=\Bigg{(}\sum_{\alpha}e^{-\epsilon_{\alpha}^{(1)}/k_{\rm
B}T}\Bigg{)}\Bigg{(}\sum_{\beta}e^{-(\epsilon_{\beta}^{(2)}/k_{\rm
B}T}\Bigg{)}=Z_{1}Z_{2}$ (A8)

The associated population computed from the Boltzmann distribution becomes
equally factorizable

$\displaystyle
p_{i}=\frac{e^{-(\epsilon_{\alpha}^{(1)}+\epsilon_{\beta}^{(2)})/k_{\rm
B}T}}{Z}=\Bigg{(}\frac{e^{-\epsilon_{\alpha}^{(1)}/k_{\rm
B}T}}{Z_{1}}\Bigg{)}\,\Bigg{(}\frac{e^{-\epsilon_{\beta}^{(2)}/k_{\rm
B}T}}{Z_{2}}\Bigg{)}$ (A9)

This simplification should lead to independent thermalization in each mode.

The reduced density matrix of a given mode can be obtained by tracing out one
DOF, e.g.,

$\displaystyle\rho^{(1)}=Tr_{(2)}\left(\rho\right)$ (A10)

In the factorizable case, we can write

$\displaystyle\rho=\sum_{\alpha\beta}\rho_{\alpha\beta}\lvert\alpha\rangle\langle\alpha\rvert\otimes\lvert\beta\rangle\langle\beta\rvert$
(A11)

where the element of the density matrix are the populations defined above,
$\rho_{\alpha\beta}=p_{i}$. Taking the partial trace over the second DOF leads
to

$\displaystyle\rho^{(1)}$ $\displaystyle=$ $\displaystyle
Tr_{(2)}\left(\sum_{\alpha\beta}\rho_{\alpha\beta}\lvert\alpha\rangle\langle\alpha\rvert\otimes\lvert\beta\rangle\langle\beta\rvert\right)$
(A12) $\displaystyle=$
$\displaystyle\sum_{\gamma}\sum_{\alpha\beta}\rho_{\alpha\beta}\lvert\alpha\rangle\langle\alpha\rvert\otimes\langle\gamma\lvert\beta\rangle\langle\beta\lvert\gamma\rangle$
$\displaystyle=$
$\displaystyle\sum_{\alpha}\left(\sum_{\beta}\rho_{\alpha\beta}\right)\lvert\alpha\rangle\langle\alpha\rvert$

Since the populations $\rho_{\alpha\beta}$ are factorizable, we can rewrite

$\displaystyle\sum_{\beta}\rho_{\alpha\beta}$ $\displaystyle=$
$\displaystyle\sum_{\beta}\frac{e^{-\epsilon_{\alpha}^{(1)}/k_{B}T}}{Z_{1}}\,\,\frac{e^{-\epsilon_{\beta}^{(1)}/k_{\rm
B}T}}{Z_{2}}$ (A13) $\displaystyle=$
$\displaystyle\frac{e^{-\epsilon_{\alpha}^{(1)}/k_{\rm
B}T}}{Z_{1}}\left(\frac{\sum_{\beta}e^{-\epsilon_{\beta}^{(2)}/k_{\rm
B}T}}{Z_{2}}\right)$ $\displaystyle=$
$\displaystyle\frac{e^{-\epsilon_{\alpha}^{(1)}/k_{\rm
B}T}}{Z_{1}}=\rho_{\alpha}^{(1)}$

This leads to the following reduced density matrix

$\displaystyle\rho^{(1)}=\sum_{\alpha}p_{\alpha}^{(1)}\lvert\alpha\rangle\langle\alpha\rvert$
(A14)

for which the populations follow a Boltzmann distribution,
$p_{\alpha}^{(1)}=\frac{e^{-\epsilon_{\alpha}^{(1)}/k_{\rm B}T}}{Z_{1}}$.
Since the two DOF are thermalized independently, we can conclude that, when
the hamiltonian is separable, i.e. the coupling vanishes, $V\rightarrow 0$,
the reduced probability are distributed thermally, at the same temperature.

### Case 2

In the more general case where the coupling does not vanish, thermalization of
the individual modes is not necessarily observed at equilibrium. Let us first
define the coupled density matrix from the coupled eigenstates in Eq. (A7),
that is

$\displaystyle\rho=\sum_{i}p_{i}\sum_{\alpha\beta\alpha^{\prime}\beta^{\prime}}c_{i,\alpha\beta}c_{i,\alpha^{\prime}\beta^{\prime}}^{\star}\lvert\alpha\rangle\langle\alpha^{\prime}\rvert\otimes\lvert\beta\rangle\langle\beta^{\prime}\lvert$
(A15)

Taking the partial trace of the coupled density matrix yields

$\displaystyle\rho^{(1)}=Tr_{(2)}\left(\rho\right)$ $\displaystyle=$
$\displaystyle\sum_{i}p_{i}\sum_{\gamma}\sum_{\alpha\beta\alpha^{\prime}\beta^{\prime}}c_{i,\alpha\beta}c_{i,\alpha^{\prime}\beta^{\prime}}^{\star}\lvert\alpha\rangle\langle\alpha^{\prime}\rvert\otimes\langle\gamma\lvert\beta\rangle\langle\beta^{\prime}\lvert\gamma\rangle$
(A16) $\displaystyle=$
$\displaystyle\sum_{i}p_{i}\sum_{\alpha\alpha^{\prime}\beta}c_{i,\alpha\beta}c^{\star}_{i,\alpha^{\prime}\beta}\lvert\alpha\rangle\langle\alpha^{\prime}\lvert$
$\displaystyle=$
$\displaystyle\sum_{\alpha\alpha^{\prime}}\left(\sum_{i\beta}p_{i}c_{i,\alpha\beta}c^{\star}_{i,\alpha^{\prime}\beta}\right)\lvert\alpha\rangle\langle\alpha^{\prime}\lvert$
$\displaystyle=$
$\displaystyle\sum_{\alpha\alpha^{\prime}}\rho_{\alpha\alpha^{\prime}}^{(1)}\lvert\alpha\rangle\langle\alpha^{\prime}\lvert$

This represents a fully coherent reduced density matrix, in which all states
$\lvert\alpha\rangle$ are coherently interacting with all other states
$\lvert\alpha^{\prime}\rangle$ of a given mode via the other DOF. The diagonal
elements of the reduced density matrix can be recast as

$\displaystyle\rho_{\alpha\alpha}=\sum_{i\beta}p_{i}c_{i,\alpha\beta}c^{\star}_{i,\alpha\beta}=\sum_{i}p_{i}\sum_{\beta}\lvert
c_{i,\alpha\beta}\rvert^{2}=\sum_{i}p_{i}w_{i,\alpha}$ (A17)

with weights $w_{i,\alpha}=\sum_{\beta}\lvert c_{i,\alpha\beta}\rvert^{2}$.
These diagonal elements are populations of the different states
$\lvert\alpha\rangle$ In general, the population do not have to obey a
Boltzmann distribution because the weight $w_{i,\alpha}$ are different for
different values of $\\{i,\alpha\\}$, and they depend on the coupling strength
$V$ and the other mode indirectly via the expansion coefficients
$c_{i,\alpha\beta}$.