# Quasi-Diabatic Propagation Scheme for Simulating Polariton Chemistry

Deping Hu<EMAIL_ADDRESS>Department of Chemistry, University of
Rochester, 120 Trustee Road, Rochester, NY 14627, U.S.A.  Arkajit Mandal
Department of Chemistry, University of Rochester, 120 Trustee Road, Rochester,
NY 14627, U.S.A.  Department of Chemistry, Columbia University, New York, NY
10027, U.S.A.  Braden M. Weight  Department of Physics and Astronomy,
University of Rochester, Rochester, NY 14627, U.S.A.  Pengfei Huo
<EMAIL_ADDRESS>Department of Chemistry, University of Rochester,
120 Trustee Road, Rochester, NY 14627, U.S.A.  The Institute of Optics, Hajim
School of Engineering, University of Rochester, Rochester, New York, 14627,
U.S.A.

###### Abstract

We generalize the quasi-diabatic (QD) propagation scheme to simulate the non-
adiabatic polariton dynamics in molecule-cavity hybrid systems. The adiabatic-
Fock states, which are the tensor product states of the adiabatic electronic
states of the molecule and photon Fock states, are used as the locally well-
defined diabatic states for the dynamics propagation. These locally well-
defined diabatic states allow using any diabatic quantum dynamics methods for
dynamics propagation, and the definition of these states will be updated at
every nuclear time step. We use several recently developed non-adiabatic
mapping approaches as the diabatic dynamics methods to simulate polariton
quantum dynamics in a Shin-Metiu model coupled to an optical cavity. The
results obtained from the mapping approaches provide very accurate population
dynamics compared to the numerically exact method and outperform the widely
used mixed quantum-classical approaches, such as the Ehrenfest dynamics and
the fewest switches surface hopping approach. We envision that the generalized
QD scheme developed in this work will provide a powerful tool to perform the
non-adiabatic polariton simulations by allowing a direct interface between the
diabatic dynamics methods and ab initio polariton information.

##  I Introduction

Coupling molecules to the quantized radiation field inside an optical cavity
creates a set of new photon-matter hybrid states, which are commonly referred
to as polaritons,[1, 2, 3, 4, 5, 6] which have been shown to facilitate new
chemical reactivities.[7, 1, 8, 6, 9] Theoretical investigations play a
crucial role in understanding the fundamental limit and basic principles in
this emerging field,[6, 5, 10, 11, 12, 12, 13] as these polariton chemical
reactions often involve a rich dynamical interplay among the electronic,
nuclear, and photonic degrees of freedom (DOFs). Accurately simulating
polaritonic quantum dynamics remains a challenging task and is beyond the
scope of photochemistry or quantum optics.[2]

The trajectory-based non-adiabatic dynamics approaches[14, 15, 16] play an
important role in simulating the non-adiabatic dynamics of the coupled
electronic-nuclear DOFs. Two of the most commonly used mixed quantum-classical
(MQC) methods are the Ehrenfest and fewest switches surface hopping (FSSH)
approaches.[17, 18] Both approaches describe the electronic subsystem quantum
mechanically, and treat the nuclear DOFs classically. It is thus a natural
idea for the theoretical chemistry community to extend these two approaches to
investigate polariton chemistry by treating the electronic-photonic DOFs (or
so-called polariton subsystem) quantum mechanically and the nuclear DOFs
classically. Incorporating the description of the photon field into the MQC
methods has become a basic strategy to simulate polariton chemistry.[10, 11,
12, 13, 19, 20, 21, 22] The key ingredient in the MQC simulations of polariton
dynamics is the expression of the nuclear gradient. Recently, we derived a
rigorous expression of the nuclear gradient using the quantum electrodynamics
(QED) Hamiltonian without making the usual approximations,[23] such as the
rotating wave approximation. These gradient expressions, together with the
corresponding MQC approaches (Ehrenfest and FSSH approaches), is valid for any
any number of electronic states or Fock states at any light-matter coupling
strength. However, the inherent semi-classical approximation in these
approaches can lead to the break-down of detailed balance[24] (incorrect long
time population) in Ehrenfest dynamics and the creation of artificial
electronic coherence[25] or incorrect chemical kinetics[25] for the FSSH
dynamics without invoking ad hoc decoherence corrections.

In response to these theoretical challenges, a wide range of non-adiabatic
dynamics approaches have been developed in the diabatic representation. Many
of them belong to the family of non-adiabatic mapping dynamics which are based
on the Meyer-Miller-Stock-Thoss (MMST) mapping formalism.[26, 27, 28] These
methods include partial linearized density matrix[29, 30] (PLDM), symmetrical
quasi-classical[31, 32] (SQC), the quantum-classical Liouville equation (QCLE)
dynamics.[33, 34] In particular, the recently developed $\gamma$-SQC has been
shown[35] to provide impressively accurate non-adiabatic photo-dissociation
quantum dynamics with coupled Morse potentials through the adjusted zero-point
energy (ZPE) parameter of the mapping variables, thus appearing to be a
promising method to simulate on-the-fly quantum dynamics of complex molecular
systems. In addition, the spin-mapping Linearized Semi-Classical approach
(spin-LSC)[36, 37, 38], which uses generalized spin mapping representation[37]
for the electronic DOF as well as the Linearization approximation[39, 40] for
the nuclear DOF, has also shown a significant improvement of the population
dynamics in the system-bath model problems (such as in spin-boson systems[36]
and many-state exciton Hamiltonians of light-harvesting complexes.[37]) The
$\gamma$-SQC approach has already demonstrated[41] its ability to outperform
MQC approaches (Ehrenfest and FSSH) in describing the electronic non-adiabatic
dynamics for ab initio on-the-fly simulations. These new mapping approaches
should, in principle, also outperform the MQC methods in simulating the
polaritonic non-adiabatic dynamics that happens in the electron-photon
subspace coupling to the motion of the nuclei. Unfortunately, to the best of
our knowledge, there are only limited studies of using mapping dynamics to
investigate polariton chemistry for model systems with strict diabatic
states.[6, 42, 9]

Recently, we have developed the quasi-diabatic (QD) propagation scheme[43, 44,
45, 46, 47, 41] as a general framework to seamlessly combine a diabatic
quantum dynamics approach, such as the mapping based methods,[35, 37] with the
adiabatic outputs of an electronic structure method. The QD propagation scheme
uses the adiabatic states at a reference nuclear geometry (the so-called
“crude adiabatic” states) as the locally well-defined diabatic states during a
short-time propagation and then dynamically updates the QD basis at each
consecutive nuclear propagation step. In this propagation scheme, one does not
construct a global diabatic representation but instead, uses a sequence of
locally diabatic representations (one for each short-time segment) to
propagate the dynamics. We have both analytically[43] and numerically[45, 46]
demonstrated that the QD scheme provides exactly the same results compared to
the direct diabatic quantum dynamics at the single trajectory level.

In this work, we generalize the QD propagation scheme to simulate polariton
non-adiabatic dynamics in a molecule-cavity hybrid system. In particular, we
use the adiabatic-Fock state at a reference nuclear geometry as the locally
well-defined diabatic basis to propagate the polariton dynamics, and
dynamically update the definition of these local diabatic states between two
consecutive propagation steps. These adiabatic-Fock states are tensor products
of the electronic adiabatic states for the molecular system and the Fock
states of the photon field inside an optical cavity. We use the Shin-Metiu
(SM) model[48, 49] as the “ab initio” model molecular system to investigate
strong and ultra-strong light-matter interactions between a molecule and an
optical cavity. Through numerical simulations, we demonstrate the accuracy of
using both $\gamma$-SQC[35] and spin-LSC[36, 37] to obtain non-adiabatic
polariton dynamics, which outperforms widely used MQC approaches.

##  II Theory and Methods

###  II.1 The Pauli-Fierz QED Hamiltonian

The Pauli-Fierz (PF) QED Hamiltonian for one molecule coupled to quantized
radiation field inside an optical cavity can be written as

$\hat{H}=\hat{T}_{\mathrm{n}}+\hat{H}_{\mathrm{en}}+\hat{H}_{\mathrm{p}}+\hat{H}_{\mathrm{enp}}+\hat{H}_{\mathrm{d}},$
(1)

where $\hat{T}_{\mathrm{n}}$ represents the nuclear kinetic energy operator,
$\hat{H}_{\mathrm{en}}$ is the electronic Hamiltonian that describes electron-
nucleus interactions. Further, $\hat{H}_{\mathrm{p}}$,
$\hat{H}_{\mathrm{enp}}$, and $\hat{H}_{\mathrm{d}}$ represent the photonic
Hamiltonian, electronic-nuclear-photonic interactions, and the dipole self-
energy (DSE) term, respectively. A full derivation of this Hamiltonian, as
well as its connection with the various atomic cavity QED models can be found
in the Appendix of Ref. 42.

The electronic-nuclear potential $\hat{H}_{\mathrm{en}}$, which describes the
common molecular Hamiltonian excluding the nuclear kinetic energy is described
as follows

$\hat{H}_{\mathrm{en}}=\hat{T}_{\mathrm{e}}+\hat{V}_{\mathrm{ee}}+\hat{V}_{\mathrm{en}}+\hat{V}_{\mathrm{nn}}.$
(2)

The above expression includes electronic kinetic energy, electron-electron
interaction, electron-nucleus interaction and nucleus-nucleus interaction. The
expressions of these four terms can be found in previous work.[50, 51, 52]
Modern electronic structure theory have been developed around solving the
eigenvalue problem of $\hat{H}_{\mathrm{en}}$, providing the following
electronically adiabatic energy and its corresponding state

$\hat{H}_{\mathrm{en}}|\phi_{\alpha}(\mathbf{R})\rangle=E_{\alpha}({\mathbf{R}})|\phi_{\alpha}(\mathbf{R})\rangle.$
(3)

Here, $|\phi_{\alpha}(\mathbf{R})\rangle$ represents the
$\alpha_{\mathrm{th}}$ many-electron adiabatic state for a given molecular
system, with the adiabatic energy $E_{\alpha}({\mathbf{R}})$.

For clarity, we restrict our discussions to the cavity with only one photonic
mode, and all the formulas presented here can be easily generalized into a
more realistic, many-mode cavity. The photonic Hamiltonian is written as

$\hat{H}_{\mathrm{p}}=\frac{1}{2}\left(\hat{p}_{\mathrm{c}}^{2}+\omega_{\mathrm{c}}^{2}\hat{q}_{\mathrm{c}}^{2}\right)=\hbar\omega_{\mathrm{c}}\Big{(}\hat{a}^{\dagger}\hat{a}+\frac{1}{2}\Big{)},$
(4)

where
$\hat{q}_{\mathrm{c}}=\sqrt{\hbar/2\omega_{\mathrm{c}}}(\hat{a}^{\dagger}+\hat{a})$
and
$\hat{p}_{\mathrm{c}}=i\sqrt{\hbar\omega_{\mathrm{c}}/2}(\hat{a}^{\dagger}-\hat{a})$
are photon field operators, $\hat{a}^{\dagger}$ and $\hat{a}$ are the photonic
creation and annihilation operators, respectively and $\omega_{\mathrm{c}}$ is
the photon frequency.

The light-matter coupling term (electronic-nuclear-photonic interactions)
under the dipole gauge is expressed as

$\hat{H}_{\mathrm{enp}}=\omega_{\mathrm{c}}\hat{q}_{\mathrm{c}}(\bm{\lambda}\cdot{\hat{\bm{\mu}}})=g_{\mathrm{c}}\bm{\epsilon}\cdot{\hat{\bm{\mu}}}(\hat{a}^{\dagger}+\hat{a}).$
(5)

where $\bm{\lambda}=\lambda\cdot{\bm{\epsilon}}$ characterizes the cavity
photon field strength, ${\bm{\epsilon}}$ is the direction of the field
polarization. The photon field strength is determined by the volume of the
cavity as $\lambda=\sqrt{1/\epsilon_{0}\mathcal{V}_{0}}$, where $\epsilon_{0}$
is the permittivity inside the cavity and $\mathcal{V}_{0}$ is the effective
quantization volume inside the cavity. Another way to characterize the light-
matter coupling strength is using
$g_{\mathrm{c}}=\sqrt{\hbar\omega_{\mathrm{c}}/2}\lambda$. Note that the
common notation used in the literature,[53, 6] the definition of
$g_{\mathrm{c}}$ also includes $\bm{\lambda}\cdot{\hat{\bm{\mu}}}$. Further,
the total dipole operator of both electrons and nuclei is defined as

${\hat{\bm{\mu}}}=-\sum_{i}e\hat{{\bf r}}_{i}+\sum_{j}Z_{j}e\hat{\bf{R}}_{j},$
(6)

where $-e$ is the charge of the electron and $Z_{j}e$ is the charge of the
$j_{\mathrm{th}}$ nucleus.

Finally, the DSE term is expressed as

$\hat{H}_{\mathrm{d}}=\frac{1}{2}(\bm{\lambda}\cdot\hat{\bm{\mu}})^{2}=\frac{g_{\mathrm{c}}^{2}}{\hbar\omega_{\mathrm{c}}}(\bm{\epsilon}\cdot\hat{\bm{\mu}})^{2}.$
(7)

This is a necessary term in the PF Hamiltonian and ensures both gauge
invariance of the Hamiltonian[54, 9] and a bounded ground state.[55, 9, 56] In
this work, we do not consider the cavity loss. The cavity loss can be
effectively incorporated by using Lindblad dynamics approaches with the MQC
simulations.[57]

For the molecule-cavity hybrid system, a convenient basis for quantum dynamics
simulations could be the photon-dressed electronic adiabatic states

$|\psi_{i}({\mathbf{R}})\rangle=|\phi_{\alpha}({\mathbf{R}})\rangle\otimes|n\rangle\equiv|\phi_{\alpha}({\bf
R}),n\rangle,$ (8)

where quantum number $i\equiv\\{\alpha,n\\}$ indicates both the adiabatic
electronic state of the molecule and the Fock state. Note that we have
introduced a shorthand notation in Eq. 8, which will be used throughout the
rest of this paper. This is one of the most straightforward choices of basis
for the hybrid system because of the readily available adiabatic electronic
information (e.g., wavefunctions, energies, and the dipole matrix) from
electronic structure calculations that we need to construct the elements of
Hamiltonian.

In the MQC simulation, such as the Ehrenfest or FSSH approach, or the recently
developed mapping non-adiabatic approaches, the total molecular Hamiltonian is
expressed as

$\hat{H}=\hat{T}_{\mathrm{n}}+\hat{V},$ (9)

where $\hat{T}_{\mathrm{n}}$ represents the nuclear kinetic energy operator,
and $\hat{V}$ represents the rest of the Hamiltonian. For a bare molecular
system, $\hat{V}=\hat{H}_{\mathrm{en}}$ expressed in Eq. 2. For a molecule-
cavity hybrid system,

$\hat{V}=\hat{H}_{\mathrm{en}}+\hat{H}_{\mathrm{p}}+\hat{H}_{\mathrm{enp}}+\hat{H}_{\mathrm{d}}\equiv\hat{H}_{\mathrm{pl}},$
(10)

which is commonly referred to as the polariton Hamiltonian,[58, 3] also
denoted as $\hat{H}_{\mathrm{pl}}$. In a similar way that electronic adiabatic
states are defined in Eq. 3, one can further define the polaritonic state[58,
3] as the eigenstate of $\hat{V}=\hat{H}_{\mathrm{pl}}$ (see definition in Eq.
10) through the following eigenequation

$\hat{H}_{\mathrm{pl}}|\mathcal{E}_{J}({\bf R})\rangle=\mathcal{E}_{J}({\bf
R})|\mathcal{E}_{J}({\bf R})\rangle,$ (11)

where $|\mathcal{E}_{J}({\bf R})\rangle$ is the polariton state with polariton
energy $\mathcal{E}_{J}({\bf R})$. The polariton eigenstate can be expressed
as

$|\mathcal{E}_{J}({\bf
R})\rangle=\sum_{\alpha,n}c^{J}_{\alpha,n}({\mathbf{R}})|\phi_{\alpha}({\bf
R}),n\rangle,$ (12)

where $c^{J}_{\alpha,n}({\mathbf{R}})=\langle\phi_{\alpha}({\bf
R}),n|\mathcal{E}_{J}({\bf R})\rangle$ and $\mathcal{E}_{J}({\bf R})$ can be
obtained by diagonalizing the matrix of $\hat{V}=\hat{H}_{\mathrm{pl}}$
(constructed from the adiabatic-Fock state basis in Eq. 8) as

$\mathbf{U}^{\dagger}[V({\mathbf{R}})]\mathbf{U}=[\mathcal{E}({\mathbf{R}})],$
(13)

where

$[V({\mathbf{R}})]_{ij}=\langle\psi_{i}({\mathbf{R}})|\hat{V}|\psi_{j}({\mathbf{R}})\rangle.$
(14)

Note that the ${\bf R}$-dependence of $|\mathcal{E}_{J}({\bf R})\rangle$ is
entirely coming from the ${\bf R}$-dependence of the adiabatic states
$|\phi_{\alpha}({\mathbf{R}})\rangle$, and the Fock state $|n\rangle$ is
completely ${\bf R}$-independent. Meanwhile, the $|\mathcal{E}_{J}({\bf
R})\rangle$ is the eigenstate of $\hat{V}$, whereas the adiabatic state
$|\phi_{\alpha}({\mathbf{R}})\rangle$ is only the eigenstate of
$\hat{H}_{\mathrm{en}}$, and not for $\hat{V}$.

###  II.2 Quasi-Diabatic Propagation Scheme for Molecular Cavity QED

The QD propagation scheme explicitly addresses the discrepancy between
accurate quantum dynamics methods in the diabatic representation and the
electronic structure methods in the adiabatic representation. The essential
idea of the QD scheme is to use the electronic adiabatic states associated
with a reference geometry as the local diabatic states during a short-time
quantum propagation and dynamically updates the definition of the QD states
along the time-dependent nuclear trajectory.[43, 44, 45, 46, 47, 41]

In this work, we apply the QD propagation scheme to the case of molecular
cavity QED. This requires the use of a convenient basis with a reference
nuclear geometry as the $\it locally$ well-defined diabatic basis, in the
sense that its character is fixed (which is automatically guaranteed because
of the fixed reference geometry by construction) as well as it is a complete
basis (which is only true when the geometry is close to this reference
geometry). The potential candidate for this basis is the adiabatic-Fock state
$|\psi_{i}({\mathbf{R}})\rangle=|\phi_{\alpha}({\mathbf{R}}),n\rangle$ (Eq.
8), which is not the same as the polariton states $|\mathcal{E}_{J}({\bf
R})\rangle$ (Eq. 12) except for the zero-coupling limit. In this work, we use
the adiabatic-Fock state as the convenient choice due to its simplicity in
terms of the polariton coupling and nuclear gradient expressions in the QD
propagation framework.

Consider a short-time propagation of the nuclear DOFs during
$t\in[t_{0},t_{1}]$, where the nuclear positions evolve from ${\bf R}(t_{0})$
to ${\bf R}(t_{1})$, and the corresponding adiabatic-Fock basis (defined in
Eq. 8) are $\\{|\psi_{i}({\bf R}(t_{0}))\rangle\\}$ and $\\{|\psi_{j}{\bf
R}(t_{1}))\rangle\\}$. We uses the basis $\\{|\psi_{i}({\bf
R}_{0})\rangle\equiv|\phi_{\alpha}({\mathbf{R}_{0}}),n\rangle\\}$ at the
reference nuclear geometry ${\bf R}(t_{0})$ as the diabatic basis during this
short-time propagation, such that

$|\psi_{i}({\bf R}_{0})\rangle\equiv|\psi_{i}({\bf
R}(t_{0}))\rangle,\leavevmode\nobreak\ \leavevmode\nobreak\
\mathrm{for}\leavevmode\nobreak\ t\in[t_{0},t_{1}].$ (15)

With the above QD basis defined independently of ${\bf R}(t)$ within each
propagation segment, the electronic derivative couplings vanish while
$\hat{V}({\bf R}(t))$ in the QD basis becomes off-diagonal. With this local
diabatic basis, all of the necessary diabatic quantities can be evaluated and
used to propagate quantum dynamics during $t\in[t_{0},t_{1}]$.

During this propagation step, the matrix element of $\hat{V}$ in the QD basis
is evaluated as

$V_{\alpha\beta,mn}({\bf R}(t))=\langle\phi_{\alpha}({\bf
R}_{0}),m|\hat{V}({\bf R}(t))|\phi_{\beta}({\bf R}_{0}),n\rangle.$ (16)

For on-the-fly simulations, this quantity is obtained from a linear
interpolation[59] between $V_{{\alpha}{\beta},mn}({\bf R}_{0})$ and
$V_{\alpha\beta,mn}({\bf R}(t_{1}))$ as follows

$\displaystyle V_{\alpha\beta,mn}({\bf R}(t))=V_{\alpha\beta,mn}({\bf R}_{0})$
(17) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+\frac{(t-t_{0})}{(t_{1}-t_{0})}\big{[}V_{\alpha\beta,mn}({\bf
R}(t_{1}))-V_{\alpha\beta,mn}({\bf R}_{0})\big{]}.$

The above linear interpolation scheme can be further improved in future work
and one potential choice is the recently developed norm-preserving
interpolation scheme.[60, 61]

It is straightforward to evaluate $V_{{\alpha}{\beta},mn}({\bf R}_{0})$ and
$V_{\alpha\beta,mn}({\bf R}(t_{1}))$ separately for the molecule-cavity hybrid
system, as discussed below. Using electronic ab initio calculation, as well as
the properties of $\hat{a}^{\dagger}$ and $\hat{a}$ for the photonic DOF, we
can explicitly evaluate each term of $V_{\alpha\beta,mn}({\bf R}_{0})$ (see
Eq. 10) as follows

$\displaystyle H^{\mathrm{en}}_{\alpha\beta,mn}({\bf
R}_{0})=\langle\phi_{\alpha}({\bf R}_{0}),m|\hat{H}_{\mathrm{en}}({\bf
R}_{0})|\phi_{\beta}({\bf R}_{0}),n\rangle$ $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ =E_{\alpha}({\bf
R}_{0})\delta_{{\alpha},{\beta}}\delta_{m,n}$ (18a) $\displaystyle
H^{\mathrm{p}}_{\alpha\beta,mn}({\bf R}_{0})=\langle\phi_{\alpha}({\bf
R}_{0}),m|\hat{H}_{\mathrm{p}}|\phi_{\beta}({\bf R}_{0}),n\rangle$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
=\hbar\omega_{\mathrm{c}}(n+\frac{1}{2})\delta_{\alpha,\beta}\delta_{m,n}$
(18b) $\displaystyle H^{\mathrm{enp}}_{\alpha\beta,mn}({\bf
R}_{0})=\langle\phi_{\alpha}({\bf R}_{0}),m|\hat{H}_{\mathrm{enp}}({\bf
R}_{0})|\phi_{\beta}({\bf R_{0}}),n\rangle$ $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
=g_{\mathrm{c}}\bm{\epsilon}\cdot{\bm{\mu}}_{\alpha\beta}({\bf
R}_{0})\big{(}\sqrt{n}\delta_{m,n-1}+\sqrt{n+1}\delta_{m,n+1}\big{)}$ (18c)
$\displaystyle H^{\mathrm{d}}_{\alpha\beta,mn}({\bf
R}_{0})=\langle\phi_{\alpha}({\bf R}_{0}),m|\hat{H}_{\mathrm{d}}({\bf
R}_{0})|\phi_{\beta}({\bf R}_{0}),n\rangle$ $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
=\frac{g_{\mathrm{c}}^{2}}{\hbar\omega_{\mathrm{c}}}\sum_{\gamma}(\bm{\epsilon}\cdot{\bm{\mu}}_{\alpha\gamma}({\bf
R}_{0}))(\bm{\epsilon}\cdot{\bm{\mu}}_{\gamma\beta}({\bf R}_{0}))\delta_{m,n}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\equiv D^{2}_{\alpha\beta}({\bf R}_{0})\delta_{m,n},$ (18d)

where $\hat{H}_{\mathrm{p}}$ (see its definition in Eq. 4) is an ${\bf
R}$-independent operator, the sum $\sum_{\gamma}$ in the matrix element of
$\hat{H}_{\mathrm{d}}$ runs over the diabatic states, and
$D^{2}_{\alpha\beta}$ denotes the elements of DSE. Further, the matrix element
of the dipole operator under the diabatic representation is expressed as

${\bm{\mu}}_{\alpha\beta}({\bf R}_{0})\equiv\langle\phi_{\alpha}({\bf
R}_{0})|\hat{\bm{\mu}}({\bf R}_{0})|\phi_{\beta}({\bf R}_{0})\rangle.$ (19)

Similarly, at time $t_{1}$, the matrix element $V_{\alpha\beta,mn}({\bf
R}(t_{1}))=\langle\phi_{\alpha}({\bf R}_{0}),m|\hat{V}({\bf
R}(t_{1}))|\phi_{\beta}({\bf R}_{0}),n\rangle$ can also be written explicitly,
with each term expressed as follows

$\displaystyle H^{\mathrm{en}}_{\alpha\beta,mn}({\bf
R}(t_{1}))=\langle\phi_{\alpha}({\bf R}_{0}),m|\hat{H}_{\mathrm{en}}({\bf
R}(t_{1}))|\phi_{\beta}({\bf R}_{0}),n\rangle$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
=H^{\mathrm{en}}_{\alpha\beta}({\bf R}(t_{1}))\delta_{m,n}$ (20a)
$\displaystyle H^{\mathrm{p}}_{\alpha\beta,mn}({\bf
R}(t_{1}))=\langle\phi_{\alpha}({\bf
R}_{0}),m|\hat{H}_{\mathrm{p}}|\phi_{\beta}({\bf R}_{0}),n\rangle$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
=\hbar\omega_{\mathrm{c}}(n+\frac{1}{2})\delta_{\alpha,\beta}\delta_{m,n}$
(20b) $\displaystyle H^{\mathrm{enp}}_{\alpha\beta,mn}({\bf
R}(t_{1}))=\langle\phi_{\alpha}({\bf R}_{0}),m|\hat{H}_{\mathrm{enp}}({\bf
R}(t_{1}))|\phi_{\beta}({\bf R}_{0}),n\rangle$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
=g_{\mathrm{c}}\bm{\epsilon}\cdot{\bm{\mu}}_{\alpha\beta}({\bf
R}(t_{1}))\big{(}\sqrt{n}\delta_{m,n-1}+\sqrt{n+1}\delta_{m,n+1}\big{)}$ (20c)
$\displaystyle H^{\mathrm{d}}_{\alpha\beta,mn}({\bf
R}(t_{1}))=\langle\phi_{\alpha}({\bf R}_{0}),m|\hat{H}_{\mathrm{d}}({\bf
R}(t_{1}))|\phi_{\beta}({\bf R_{0}}),n\rangle$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
=\frac{g_{\mathrm{c}}^{2}}{\hbar\omega}\sum_{\gamma}(\bm{\epsilon}\cdot{\bm{\mu}}_{\alpha\gamma}({\bf
R}(t_{1})))(\bm{\epsilon}\cdot{\bm{\mu}}_{\gamma\beta}({\bf
R}(t_{1})))\delta_{m,n}$ $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \equiv D^{2}_{\alpha\beta}({\bf
R}(t_{1}))\delta_{m,n},$ (20d)

where $H^{\mathrm{en}}_{\alpha\beta}({\bf
R}(t_{1}))\equiv\langle\phi_{\alpha}({\bf R}_{0})|\hat{H}_{\mathrm{en}}({\bf
R}(t_{1}))|\phi_{\beta}({\bf R}_{0})\rangle$, and
${\bm{\mu}}_{\alpha\beta}({\bf R}(t_{1}))\equiv\langle\phi_{\alpha}({\bf
R}_{0})|\hat{\bm{\mu}}({\bf R}(t_{1}))|\phi_{\beta}({\bf R}_{0})\rangle$.

To conveniently calculate $H^{\mathrm{en}}_{\alpha\beta}({\bf R}(t_{1}))$ and
${\bm{\mu}}_{\alpha\beta}({\bf R}(t_{1}))$, we use the following relations

$\displaystyle H^{\mathrm{en}}_{\alpha\beta}({\bf R}(t_{1}))$
$\displaystyle=\sum_{\lambda\nu}S_{\alpha\lambda}\tilde{H}^{\mathrm{en}}_{\lambda\nu}({\bf
R}(t_{1}))S^{\dagger}_{\beta\nu}$ (21a)
$\displaystyle{\bm{\mu}}_{\alpha\beta}({\bf R}(t_{1}))$
$\displaystyle=\sum_{\lambda\nu}S_{\alpha\lambda}\tilde{{\bm{\mu}}}_{\lambda\nu}({\bf
R}(t_{1}))S^{\dagger}_{\beta\nu},$ (21b)

where the matrix elements at ${\bf R}(t_{1})$ are expressed as

$\displaystyle\tilde{H}^{\mathrm{en}}_{\lambda\nu}({\bf R}(t_{1}))$
$\displaystyle=\langle\phi_{\lambda}({\bf
R}(t_{1}))|\hat{H}_{\mathrm{en}}({\bf R}(t_{1}))|\phi_{\nu}({\bf
R}(t_{1}))\rangle$ (22a) $\displaystyle=E_{\lambda}({\bf
R}(t_{1}))\delta_{\lambda\nu}$
$\displaystyle\tilde{\bm{\mu}}_{\lambda\nu}({\bf R}(t_{1}))$
$\displaystyle=\langle\phi_{\lambda}({\bf R}(t_{1}))|\hat{\bm{\mu}}({\bf
R}(t_{1}))|\phi_{\nu}({\bf R}(t_{1}))\rangle,$ (22b)

and the overlap matrix between two electronic adiabatic states (with two
different nuclear geometries) are

$\displaystyle S_{\alpha\lambda}$ $\displaystyle=\langle\phi_{\alpha}({\bf
R}_{0})|\phi_{\lambda}({\bf R}(t_{1}))\rangle$ (23a) $\displaystyle
S^{\dagger}_{\beta\nu}$ $\displaystyle=\langle\phi_{\nu}({\bf
R}(t_{1}))|\phi_{\beta}({\bf R}_{0})\rangle.$ (23b)

Using the above information, as well as Eq. 17, we can obtain each tern of
$V_{\alpha\beta,mn}({\bf R}(t))$ for propagating the dynamics of the quantum
subsystem that contains both electronic and photonic DOFs.

Next, we need to evaluate the nuclear gradients to propagate the dynamics of
the classical subsystem, which contains the nuclear DOFs. In particular, we
need to evaluate the nuclear gradients on each term of $\nabla
V_{\alpha\beta,mn}({\bf R}(t_{1}))$. First, let us focus on the gradient term
from the electron-nuclear coupling term[47] as follows

$\displaystyle{\nabla}H^{\mathrm{en}}_{\alpha\beta,mn}({\bf R}(t_{1}))$ (24)
$\displaystyle={\nabla}\langle\phi_{\alpha}({\bf
R}(t_{0})),m|\hat{H}_{\mathrm{en}}({\bf R}(t_{1}))|\phi_{\beta}({\bf
R}(t_{0})),n\rangle$ $\displaystyle=\langle\phi_{\alpha}({\bf
R}(t_{0}))|{\nabla}\hat{H}_{\mathrm{en}}({\bf R}(t_{1}))|\phi_{\beta}({\bf
R}(t_{0}))\rangle\cdot\langle m|n\rangle$
$\displaystyle=\sum_{\lambda\nu}S_{\alpha\lambda}\langle\phi_{\lambda}({{\bf
R}({t_{1}})})|\nabla{\hat{H}}_{\mathrm{en}}({\bf R}(t_{1}))|\phi_{\nu}({{\bf
R}(t_{1})})\rangle S^{\dagger}_{\beta\nu}\cdot\delta_{mn}.$

Here, from the first to the second line, we have used the fact that neither
$\langle\phi_{\alpha}({\bf R}_{0})|$ nor $\langle m|$ are ${\bf R}$-dependent,
which allows ${\nabla}$ to bypass both and directly act on $\hat{V}({\bf
R}(t_{1}))$. We have also used the fact that $\hat{H}_{\mathrm{en}}$ does not
contain any photonic operators. The gradient term $\langle\phi_{\lambda}({{\bf
R}({t_{1}})})|\nabla{\hat{H}}_{\mathrm{en}}({\bf R}(t_{1}))|\phi_{\nu}({{\bf
R}(t_{1})})\rangle$ can be evaluated using the following well-known
equality[18]

$\langle\phi_{\lambda}({\bf R})|\nabla\hat{H}_{\mathrm{en}}({\bf
R})|\phi_{\nu}({\bf R})\rangle=\begin{cases}\nabla
E_{\lambda}\quad&(\lambda=\nu)\\\ {\bf
d}_{\lambda\nu}\left(E_{\nu}-E_{\lambda}\right)\quad&(\lambda\neq\nu),\end{cases}$
(25)

where the non-adiabatic coupling (NAC) vector (or so-called derivative
coupling) is

${\bf d}_{\lambda\nu}=\langle\phi_{\lambda}({\bf R})|\nabla|\phi_{\nu}({\bf
R})\rangle.$ (26)

For the gradient on the matrix $H^{\mathrm{p}}_{\alpha\beta,mn}$, because
there is no nuclear DOF in $\hat{H}_{\mathrm{p}}$, thus

$\nabla H^{\mathrm{p}}_{\alpha\beta,mn}({\bf
R}({t_{1}}))=\nabla\big{[}\hbar\omega_{\mathrm{c}}(n+\frac{1}{2})\delta_{\alpha,\beta}\delta_{m,n}\big{]}=0.$
(27)

For the gradient on the light-matter interaction term
$H^{\mathrm{enp}}_{\alpha\beta,mn}$, we have

$\displaystyle\nabla H^{\mathrm{enp}}_{\alpha\beta,mn}({\bf
R}(t_{1}))={\nabla}\langle\phi_{\alpha}({\bf
R}(t_{0})),m|\hat{H}_{\mathrm{enp}}({\bf R}(t_{1}))|\phi_{\beta}({\bf
R}(t_{0})),n\rangle$ (28) $\displaystyle=\langle\phi_{\alpha}({\bf
R}(t_{0}))|{\nabla}\hat{\bm{\mu}}({\bf R}(t_{1}))|\phi_{\beta}({\bf
R}(t_{0}))\rangle\cdot\bm{\epsilon}\cdot
g_{\mathrm{c}}\big{(}\sqrt{n}\delta_{m,n-1}+\sqrt{n+1}\delta_{m,n+1}\big{)}$
$\displaystyle=\sum_{\lambda\nu}S_{\alpha\lambda}\langle\phi_{\lambda}({{\bf
R}({t_{1}})})|\nabla\hat{\bm{\mu}}({\bf R}(t_{1}))|\phi_{\nu}({{\bf
R}(t_{1})})\rangle S^{\dagger}_{\beta\nu}\cdot\bm{\epsilon}\cdot
g_{\mathrm{c}}\big{(}\sqrt{n}\delta_{m,n-1}+\sqrt{n+1}\delta_{m,n+1}\big{)},$

where $S_{\alpha\lambda}$ and $S^{\dagger}_{\beta\nu}$ are defined in Eq. 23a
and Eq. 23b, respectively. To evaluate the term $\langle\phi_{\lambda}({\bf
R})|\nabla\hat{\bm{\mu}}({\bf R})|\phi_{\nu}({\bf R)}\rangle$ that appears in
Eq. 28, we use a simple relation based on the chain rule as follows

$\displaystyle\langle\phi_{\lambda}({\bf R})|\nabla\hat{\bm{\mu}}({\bf
R})|\phi_{\nu}({\bf R})\rangle=\nabla\langle\phi_{\lambda}({\bf
R})|\hat{\bm{\mu}}({\bf R})|\phi_{\nu}({\bf
R})\rangle-\langle\nabla\phi_{\lambda}({\bf R})|\hat{\bm{\mu}}({\bf
R})|\phi_{\nu}({\bf R})\rangle-\langle\phi_{\lambda}({\bf
R})|\hat{\bm{\mu}}({\bf R})|\nabla\phi_{\nu}({\bf R})\rangle$
$\displaystyle=\nabla{\bm{\mu}}_{\lambda\nu}({\bf
R})-\sum_{\gamma}\langle\nabla\phi_{\lambda}({\bf R})|\phi_{\gamma}({\bf
R})\rangle{\bm{\mu}}_{\gamma\nu}({\bf
R})-\sum_{\gamma}{\bm{\mu}}_{\lambda\gamma}({\bf R})\langle\phi_{\gamma}({\bf
R})|\nabla\phi_{\nu}({\bf R})\rangle$
$\displaystyle=\nabla{\bm{\mu}}_{\lambda\nu}({\bf R})+\sum_{\gamma}\big{[}{\bf
d}_{\lambda\gamma}({\bf R}){\bm{\mu}}_{\gamma\nu}({\bf
R})-{\bm{\mu}}_{\lambda\gamma}({\bf R}){\bf d}_{\gamma\nu}({\bf R})\big{]},$
(29)

where ${\bf d}_{\lambda\gamma}({\bf R})=\langle\phi_{\lambda}({\bf
R})|\nabla\phi_{\gamma}({\bf R})\rangle=-\langle\nabla\phi_{\lambda}({\bf
R})|\phi_{\gamma}({\bf R})\rangle$ and ${\bf d}_{\gamma\nu}({\bf
R})=\langle\phi_{\gamma}({\bf R})|\nabla\phi_{\nu}({\bf
R})\rangle=-\langle\nabla\phi_{\gamma}({\bf R})|\phi_{\nu}({\bf R})\rangle$
are the electronic derivative couplings (defined in Eq. 26). Note that from
the first line to the second line, we have inserted
$\hat{\mathcal{P}}=\sum_{\gamma}|\phi_{\gamma}({\bf
R})\rangle\langle\phi_{\gamma}({\bf R})|$, which is the resolution of identity
in the electronic subspace. As one can clearly see, this term requires the
evaluation of derivative coupling ${\bf d}_{\lambda\gamma}({\bf R})$ and ${\bf
d}_{\gamma\nu}({\bf R})$, and the derivative on dipole matrix element
$\nabla{\bm{\mu}}_{\lambda\nu}({\bf R})$. For most of the electronic structure
methods, the derivatives on dipole matrix elements
$\nabla{\bm{\mu}}_{\lambda\nu}({\bf R})$ are not implemented. Nevertheless,
recent theoretical development has made these quantities available.[22]

Finally, the gradient from the DSE term is expressed as follows

$\displaystyle\nabla H^{\mathrm{d}}_{\alpha\beta,mn}({\bf
R}(t_{1}))={\nabla}\langle\phi_{\alpha}({\bf
R}_{0}),m|\hat{H}_{\mathrm{d}}({\bf R}(t_{1}))|\phi_{\beta}({\bf
R}_{0}),n\rangle$ (30) $\displaystyle=$
$\displaystyle{\nabla}\Big{[}\sum_{\gamma}(\bm{\epsilon}\cdot\langle\phi_{\alpha}({\bf
R}_{0})|\hat{\bm{\mu}}({\bf R}(t_{1}))|\phi_{\gamma}({\bf
R}_{0})\rangle)(\bm{\epsilon}\cdot\langle\phi_{\gamma}({\bf
R}_{0})|\hat{\bm{\mu}}({\bf R}(t_{1}))|\phi_{\beta}({\bf
R}_{0})\rangle)\Big{]}\frac{g_{\mathrm{c}}^{2}}{\hbar\omega}\delta_{m,n}$
$\displaystyle=$
$\displaystyle\Big{[}\sum_{\gamma}(\bm{\epsilon}\cdot\langle\phi_{\alpha}({\bf
R}_{0})|{\nabla}\hat{\bm{\mu}}({\bf R}(t_{1}))|\phi_{\gamma}({\bf
R}_{0})\rangle)(\bm{\epsilon}\cdot\langle\phi_{\gamma}({\bf
R}_{0})|\hat{\bm{\mu}}({\bf R}(t_{1}))|\phi_{\beta}({\bf R}_{0})\rangle)$
$\displaystyle+(\bm{\epsilon}\cdot\langle\phi_{\alpha}({\bf
R}_{0})|\hat{\bm{\mu}}({\bf R}(t_{1}))|\phi_{\gamma}({\bf
R}_{0})\rangle)(\bm{\epsilon}\cdot\langle\phi_{\gamma}({\bf
R}_{0})|{\nabla}\hat{\bm{\mu}}({\bf R}(t_{1}))|\phi_{\beta}({\bf
R}_{0})\rangle)\Big{]}\frac{g_{\mathrm{c}}^{2}}{\hbar\omega}\delta_{m,n},$

where the term $\langle\phi_{\alpha}({\bf R}(t_{0}))|\hat{\bm{\mu}}({\bf
R}(t_{1}))|\phi_{\gamma}({\bf R}(t_{0}))\rangle$ can be evaluated using Eq.
21b, and the $\langle\phi_{\alpha}({\bf R}_{0})|{\nabla}\hat{\bm{\mu}}({\bf
R}(t_{1}))|\phi_{\beta}({\bf R}_{0})\rangle$ type of derivative can be
computed in the same fashion as elaborated in Eq. 28.

Using the matrix elements $V_{\alpha\beta,mn}({\bf R}(t))$ (Eq. 17) and the
nuclear gradient $\nabla V_{\alpha\beta,mn}({\bf R})$ (as outlined in Eqs. 24
\- 30), one can in principle use any trajectory-based approaches or wavepacket
approaches with guiding trajectories[62, 63, 64] to propagate the quantum
dynamics in the time step of $t\in[t_{0},t_{1}]$. During the next short-time
propagation segment $t\in[t_{1},t_{2}]$, the QD scheme adopts a new reference
geometry ${{\bf R}^{\prime}_{0}}\equiv{\bf R}(t_{1})$ and new diabatic basis
$|\psi_{j}({{\bf R}^{\prime}_{0}})\rangle\equiv|\psi_{j}({\bf
R}(t_{1}))\rangle=|\phi_{\beta}({\bf R}(t_{1})),m\rangle$. Between the
$t\in[t_{0},t_{1}]$ propagation and the $t\in[t_{1},t_{2}]$ propagation
segments, all of the necessary quantities will be transformed from
$\\{|\psi_{i}({\bf R}_{0})\rangle\\}$ to the $\\{|\psi_{j}({\bf
R}^{\prime}_{0})\rangle\\}$ basis, using the relation

$|\psi_{j}({\bf R}(t_{1}))\rangle=\sum_{i}\langle\psi_{i}({\bf
R}(t_{0}))|\psi_{j}({\bf R}(t_{1}))\rangle|\psi_{i}({\bf R}(t_{0}))\rangle.$
(31)

For the model calculations in this work, these overlap integrals are evaluated
as

$\langle\psi_{i}({\bf R}(t_{0}))|\psi_{j}({\bf
R}(t_{1}))\rangle=\langle\phi_{\alpha}({\bf R}(t_{0}))|\phi_{\beta}({\bf
R}(t_{1}))\rangle\cdot\langle n|m\rangle,$ (32)

where the electronic adiabatic state overlaps $\langle\phi_{\alpha}({\bf
R}(t_{0}))|\phi_{\beta}({\bf R}(t_{1}))\rangle$ are directly calculated using
the discrete variable representation (DVR) basis, and the Fock states are
orthonormal to each other $\langle n|m\rangle=\delta_{n,m}$. When performing
the transformation in Eq. 31 (as well as in Eq. 43 for the non-adiabatic
mapping methods), the eigenvectors maintain their mutual orthogonality subject
to a very small error when they are expressed in terms of the previous basis
due to the incompleteness of the basis.[65, 66] Nevertheless, the
orthogonality remains to be well satisfied among $\\{|\psi_{i}({\bf
R}(t_{0}))\rangle\\}$ or $\\{|\psi_{j}({\bf R}(t_{1}))\rangle\\}$. This small
numerical error generated from each step can, however, accumulate over many
steps and cause a significant error at longer times, leading to non-unitary
dynamics.[65, 66] This potential issue can be easily resolved by using
orthonormalization procedure among the vectors of the overlap matrix composed
by $\langle\psi_{i}({\bf R}(t_{0}))|\psi_{j}({\bf R}(t_{1}))\rangle$, as been
done in our previous work[44] for simulating photo-induced charge transfer
dynamics. Here, we perform the Löwdin orthogonalization procedure[67] as
commonly used in the local diabatization approach[65] to ensure unitary
propagation.

As the nuclear geometry closely follows the reference geometry throughout the
propagation, the QD basis forms a convenient and compact basis. Note that, in
principle, one needs an infinite set of crude adiabatic states
$\\{|\psi_{i}({\bf R}_{0})\rangle\\}$ to represent the time-dependent
electronic wavefunction because the electronic wavefunction could change
rapidly with the motion of the nuclei, and the crude adiabatic basis is only
convenient when the reference geometry ${\bf R}_{0}$ is close to the nuclear
geometry ${\bf R}$. By dynamically updating the basis in the QD scheme, the
time-dependent electronic wavefunction is expanded with the “moving crude
adiabatic basis”[62] that explores the most relevant and important parts of
the Hilbert space, thus requiring only a few states for quantum dynamics
propagation.

###  II.3 Non-adiabatic Mapping Dynamics Methods

The Meyer-Miller-Stock-Thoss (MMST) formalism[26, 27, 28] maps the discrete
electronic DOFs onto continuous phase space variables. In the strict diabatic
basis $\\{|i\rangle\\}$ (in the sense that $\langle i|\nabla|j\rangle=0$ for
all $|i\rangle$ and $|j\rangle$), the total Hamiltonian in Eq. 9 is expressed
as

$\hat{H}={\frac{{\bf P}^{2}}{2{\bf M}}}+\sum_{i}V_{ii}(\hat{\bf
R})|i\rangle\langle i|+\sum_{i\neq j}V_{ij}(\hat{\bf R})|i\rangle\langle j|,$
(33)

where $V_{ij}(\hat{\bf R})=\langle i|\hat{V}(\hat{\bf r},\hat{\bf
R})|j\rangle$ are the matrix elements of the electronic Hamiltonian. Note that
here $|i\rangle$ is used to represent the strict diabatic basis, and not to be
confused with the adiabatic-Fock state
$|\psi_{i}({\mathbf{R}})\rangle=|\phi_{\alpha}({\bf R}),n\rangle$ introduced
in Eq. 8. Nevertheless, based on the QD scheme, these adiabatic-Fock states
with a reference geometry ${\bf R}_{0}$ will be used as the diabatic state in
the neighborhood of the reference geometries, as indicated in Eq. 15.

In the non-adiabatic mapping approach, the Hamiltonian operator in Eq. 33 is
transformed into the following MMST Hamiltonian

$\mathcal{H}_{\mathrm{m}}={\frac{{\bf P}^{2}}{2{\bf
M}}}+{\frac{1}{2}}\sum_{ij}V_{ij}({\bf
R})\left(p_{i}p_{j}+q_{i}q_{j}-2\gamma_{j}\delta_{ij}\right),$ (34)

where $2\gamma_{j}$ is viewed as a parameter[68] which specifies the ZPE of
the mapping oscillators.[69, 68, 36, 37] In principle, $2\gamma_{j}$ is state-
specific and trajectory-specific.[35] The MMST mapping Hamiltonian has been
historically justified by Stock and Thoss using harmonic oscillator’s raising
and lowering operators as the mapping operator.[27, 28] Recently, it has been
derived using the $SU(N)$ Lie group theory or so-called generalized spin
mapping approach.[37]

Classical trajectories are generated based on Hamilton’s equations of motion
(EOM)

$\displaystyle\dot{q}_{j}$ $\displaystyle=$
$\displaystyle\partial\mathcal{H}_{\mathrm{m}}/\partial
p_{j};\leavevmode\nobreak\ \leavevmode\nobreak\
\dot{p}_{i}=-\partial\mathcal{H}_{\mathrm{m}}/\partial q_{i}$ (35a)
$\displaystyle\dot{\bf R}$ $\displaystyle=$
$\displaystyle\partial\mathcal{H}_{\mathrm{m}}/\partial{\bf
P};\leavevmode\nobreak\ \leavevmode\nobreak\ \dot{\bf
P}=-\partial\mathcal{H}_{\mathrm{m}}/\partial{\bf R}={\bf F},$ (35b)

with the nuclear force expressed as

${\bf F}=-{\frac{1}{2}}\sum_{ij}\nabla V_{ij}({\bf
R})\big{(}p_{i}p_{j}+q_{i}q_{j}-2\gamma_{j}\delta_{ij}\big{)}.$ (36)

Overall, the MMST mapping provides a consistent classical footing for both
electronic and nuclear DOFs, and the non-adiabatic transitions between
electronic states are captured through the classical motion of the fictitious
harmonic oscillators. The non-adiabatic dynamics obtained from this formalism
have shown good performance in the ab initio on-the-fly dynamics.[47, 70, 41]

To sample the initial electronic condition and estimate the population, it is
also convenient to use the action-angle variables,
$\\{\varepsilon_{j},\theta_{j}\\}$, which are related to the canonical mapping
variables $\\{p_{j},q_{j}\\}$ through

$\varepsilon_{j}=\frac{1}{2}\left(p_{j}^{2}+q_{j}^{2}\right);\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\theta_{j}=-\tan^{-1}\left(\frac{p_{j}}{q_{j}}\right),$ (37)

and the inverse relations

$q_{j}=\sqrt{2\varepsilon_{j}}\cos(\theta_{j});\leavevmode\nobreak\
\leavevmode\nobreak\ p_{j}=-\sqrt{2\varepsilon_{j}}\sin(\theta_{j}),$ (38)

where $\varepsilon_{j}$ is a positive-definite action variable that is
directly proportional to the mapping variables’ radius in action-space.[35]

The SQC approach calculates the population of electronic state $|j\rangle$,
which is to be evaluated as[68]

$\displaystyle\rho_{jj}(t)$ $\displaystyle=\mathrm{Tr}_{\bf
R}\left[\hat{\rho}(0)e^{i\hat{H}t/\hbar}|j\rangle\langle
j|e^{-i\hat{H}t/\hbar}\right]$ (39) $\displaystyle\approx\int
d\bm{\tau}\rho_{\mathrm{W}}({\bf P},{\bf
R})W_{i}({\bm{\varepsilon}}(0))W_{j}({\bm{\varepsilon}}(t)),$

where $\hat{\rho}(0)=\hat{\rho}_{\bf R}\otimes|i\rangle\langle i|$ is the
initial density operator, $\rho_{\mathrm{W}}({\bf P},{\bf R})$ is the Wigner
transform of $\hat{\rho}_{\bf R}$ operator for the nuclear DOFs,
${\bm{\varepsilon}}=\\{\varepsilon_{1},\varepsilon_{2},...,\varepsilon_{\mathcal{N}}\\}$
is the positive-definite action variable vector for $\mathcal{N}$ electronic
states,[35]
$W_{i}({\bm{\varepsilon}})=\delta(\varepsilon_{i}-(1+\gamma_{i}))\prod_{i\neq
j}\delta(\varepsilon_{j}-\gamma_{j})$ is the Wigner transformed action
variables,[71] and $d\bm{\tau}\equiv d{\bf P}\cdot d{\bf R}\cdot
d{\bm{\varepsilon}}\cdot d\bm{\theta}$. For practical reasons, the above delta
functions in $W_{i}({\bm{\varepsilon}})$ are broadened using a distribution
function (so-called window function) that can be used to bin the resulting
electronic action variables in action-space.[68] Further, we use the
$\gamma$-SQC approach,[35] which uses a state-specific and trajectory-specific
$\gamma_{j}$ parameter in Eq. 34 to correct the initial force according to the
initially populated state. This method has been proven to provide very
accurate non-adiabatic dynamics in model photo-dissociation problems (coupled
Morse potential), as well as outperform FSSH (with decoherence correction) in
ab initio on-the-fly simulations.[70, 41] The details of $\gamma$-SQC are
provided in Appendix A.

For the spin-LSC approach,[36, 37] one chooses a universal ZPE parameter
$2\gamma_{j}=\Gamma$ for all states and trajectories. The spin-LSC population
dynamics is calculated as

$\displaystyle\rho_{jj}(t)$ $\displaystyle=\mathrm{Tr}_{\bf
R}\left[\hat{\rho}_{R}\otimes|i\rangle\langle
i|e^{i\hat{H}t/\hbar}|j\rangle\langle j|e^{-i\hat{H}t/\hbar}\right]$ (40)
$\displaystyle\approx\int d\bm{\tau}\rho_{\mathrm{W}}({\bf P},{\bf
R})[|i\rangle\langle i|]_{\mathrm{s}}(0)\cdot[|j\rangle\langle
j|]_{\mathrm{\bar{s}}}(t),$

where the population estimators are obtained from the Stratonovich-Weyl
transformed electronic projection operators, with the expressions as
follows[37]

$\displaystyle\big{[}|i\rangle\langle
i|\big{]}_{\mathrm{s}}=\frac{1}{2}(q^{2}_{i}+p^{2}_{i}-\Gamma)$ (41a)
$\displaystyle\big{[}|j\rangle\langle
j|\big{]}_{\mathrm{\bar{s}}}=\frac{\mathcal{N}+1}{2(1+\frac{\mathcal{N}\Gamma}{2})^{2}}\cdot(q^{2}_{j}+p^{2}_{j})-\frac{1-\frac{\Gamma}{2}}{1+\frac{\mathcal{N}\Gamma}{2}}.$
(41b)

The parameter $\Gamma$ is related to the radius of the generalized Bloch
sphere $r_{\mathrm{s}}$ through
$\Gamma=\frac{2}{\mathcal{N}}(r_{\mathrm{s}}-1)$, where $\mathrm{s}$ and
$\mathrm{\bar{s}}$ are complementary indices in the Stratonovich-Weyl
transform. Among the vast parameter space, one of the best performing
choices[36, 37] is when
$r_{\mathrm{s}}=r_{\mathrm{\bar{s}}}=\sqrt{\mathcal{N}+1}$, which is referred
to as $\mathrm{s}=\mathrm{W}$, leading to a ZPE parameter

$\Gamma=\frac{2}{\mathcal{N}}(\sqrt{\mathcal{N}+1}-1),$ (42)

as well as the identical expression of $[|i\rangle\langle i|]_{\mathrm{s}}$
and $[|j\rangle\langle j|]_{\mathrm{\bar{s}}}$ in Eq. 41. We further use the
focused initial condition[36, 37] that replaces the sampling of the mapping
variables in the $d\bm{\tau}$ integral of Eq. 40 with specific values of the
mapping variables, such that $\frac{1}{2}(q^{2}_{i}+p^{2}_{i}-\Gamma)=1$ for
initially occupied state $|i\rangle$ and
$\frac{1}{2}(q^{2}_{j}+p^{2}_{j}-\Gamma)=0$ for the initially unoccupied
states $|j\rangle$. The angle variables $\\{\theta_{j}\\}$ (Eq. 37) are
randomly sampled[37] in the range of $[0,2\pi)$.

Using the QD propagation scheme, one can directly perform non-adiabatic using
both $\gamma$-SQC and spin-LSC in their original diabatic formalism, with the
information from the “ab initio” polaritonic calculations of the molecule-
cavity hybrid system. Using the schemes outlined in Sec. II.2, one can obtain
the polariton coupling $\langle\psi_{i}({\bf R}_{0})|\hat{V}({\bf
R})|\psi_{j}({\bf R}_{0})\rangle$ (see Eq. 16 and Eq. 17) and nuclear gradient
$\nabla\langle\psi_{i}({\bf R}_{0})|\hat{V}({\bf R})|\psi_{j}({\bf
R}_{0})\rangle$ (see Eq. 24-Eq. 30), which are the necessary ingredients to
solve the MMST mapping EOMs in Eq. 35 and Eq. 36. Between two propagation
steps, the QD basis is transformed from $\\{|\psi_{i}({\bf
R}(t_{0}))\rangle\equiv|\phi_{\alpha}({\mathbf{R}}(t_{0})),n\rangle\\}$ to
$\\{|\psi_{j}({\bf
R}(t_{1}))\rangle\equiv|\phi_{\beta}({\mathbf{R}}(t_{1})),m\rangle\\}$. This
leads to the corresponding transform of mapping variables between the two
consecutive QD bases as follows[43, 47]

$\displaystyle\sum_{i}q_{i}\langle\psi_{i}({\bf R}(t_{0}))|\psi_{j}({\bf
R}(t_{1}))\rangle\rightarrow q_{j}$ (43a)
$\displaystyle\sum_{i}p_{i}\langle\psi_{i}({\bf R}(t_{0}))|\psi_{j}({\bf
R}(t_{1}))\rangle\rightarrow p_{j},$ (43b)

where the overlaps between the two steps are evaluated using Eq. 32 (see
discussion under that equation). More computational details for the
$\gamma$-SQC and spin-LSC are provided in section III.2.

##  III Computational Details

###  III.1 The Model System

In this work, we use the asymmetrical Shin-Metiu model[48, 49] as the “ab
initio” model molecular system to investigate strong light-matter interactions
between a molecule and an optical cavity. The model contains a transferring
proton (nucleus) and an electron, as well as two fixed ions labeled as donor
(D) and acceptor (A), as shown in Fig. 1a. This model is usually used to
describe the proton-coupled electron transfer (PCET) reaction and has been
studied recently using the exact factorization approach to investigate how
cavity can influence chemical reactivities.[72, 49, 73] The electron-nuclear
interaction potential operator $\hat{H}_{\mathrm{en}}$ (c.f. Eq. 2) is
expressed as

$\hat{H}_{\mathrm{en}}=\sum_{\sigma=\pm 1}\left(\frac{1}{|R+\frac{\sigma
L}{2}|}-\frac{\mathrm{erf}\Big{(}\frac{|r+\frac{\sigma
L}{2}|}{a_{\sigma}}\Big{)}}{|r+\frac{\sigma
L}{2}|}\right)-\frac{\mathrm{erf}\Big{(}\frac{|R-r|}{a_{f}}\Big{)}}{|R-r|},$
(44)

where the first term represents the potential of the transferring proton, the
second term represents the potential of the transferring electron, and the
third term represents the electron-proton coupling. We choose the same
parameters used in Ref. 49, which is $L=19$ a.u., $a_{+}=3.1$ a.u.,
$a_{-}=4.0$ a.u., $a_{f}=5.0$ a.u. and the proton mass is $M=1836$ a.u. To
calculate the electronic properties of the SM model, we use the Sinc discrete
variable representation (DVR) basis[74] to represent the electronic adiabatic
states. These adiabatic states $|\phi_{\alpha}(R)\rangle$ are computed on-the-
fly for a given nuclear configuration $R$ by solving Eq. 3. The details are
provided in Appendix B.

Figure 1: (a) The schematic illustration of the asymmetrical SM model, where
one proton and one electron can transfer between two fixed ions (donor and
acceptor). The distance between the donor and acceptor is 19 a.u. Here $V_{p}$
and $V_{e}$ are the potential of the transferring proton and electron,
respectively. (b) The first two adiabatic electronic states, (c) NAC between
these two electronic states and (d) their permanent and transition dipole
moments of the SM model. The PESs of the polaritonic states inside the cavity
are obtained with the light-matter coupling strength (e)
$g_{\mathrm{c}}=0.001$ and (f) $g_{\mathrm{c}}=0.005$. The color used in (e)
and (f) is coded according to $\langle\hat{a}^{\dagger}\hat{a}\rangle$, as
shown in the upper position of panel (e).

Fig. 1a also depicts the model potential in Eq. 44, with the black curve
depicting the proton potential (the first term in Eq. 44) and the green curve
depicting the electron potential (the second term in Eq. 44). Fig. 1b presents
the two lowest adiabatic electronic states of the SM model (red and blue
curves). There is an avoided crossing between the ground and the first excited
state potential energy surfaces (PESs) near $R=2.0$ a.u.. Fig. 1c presents the
NAC between them (the green curve), which shows a strong coupling near the
avoided crossing region. The matrix elements of the dipole operator under the
adiabatic representation (Eq. 19) of the SM model are presented in Fig. 1d.

When coupling the SM molecular model with the cavity, the photon frequency of
the cavity mode is chosen as $\hbar\omega_{\mathrm{c}}=2.721$ eV
($\approx\,0.1$ a.u.). Further, we assume that the cavity field polarization
direction ${\bm{\epsilon}}$ is always aligned with the direction of the dipole
operator $\hat{\bm{\mu}}$, such that
${\bm{\epsilon}}\cdot{\bm{\mu}}_{\gamma\nu}={\mu}_{\gamma\nu}$ (for
$\\{\nu,\gamma\\}=\\{e,g\\}$) where ${\mu}_{\gamma\nu}$ is the magnitude of
$\hat{\bm{\mu}}$. Explicitly considering the angle between ${\bm{\epsilon}}$
and $\hat{\bm{\mu}}$ will generate a polariton induced conical intersection
(even for a diatomic molecule), which will induce geometric phase effects.[75]
We consider two different light-matter coupling strengths
$g_{\mathrm{c}}=0.001$ a.u. and $g_{\mathrm{c}}=0.005$ a.u. in this work. The
normalized coupling strength is often defined as[76] $\eta\equiv
g_{\mathrm{c}}\cdot|\bm{\epsilon}\cdot{{\bm{\mu}}}_{eg}|/\omega_{\mathrm{c}}$
where $|\bm{\epsilon}\cdot{{\bm{\mu}}}_{eg}|$ is the typical magnitude of the
transition dipole projected along the field polarization direction. For the
coupling strength considered above (and taking
$\hbar\omega_{\mathrm{c}}\approx 2.721$ eV for the model calculation), the
normalized coupling strength is $\eta=0.06$ (for $g_{\mathrm{c}}=0.001$ a.u.)
and $\eta=0.3$ (for $g_{\mathrm{c}}=0.005$ a.u.). When $0.1<\eta<1$, the light
and matter interaction achieves the ultra-strong coupling regime,[76, 77]
which is difficult to achieve but still within the reach of the current
experimental setup.[78, 79] Thus, besides the pure theoretical value to derive
the exact nuclear gradient expression, our computational results are also
within the reach of the near future experimental setup.

The polaritonic PESs $\mathcal{E}_{J}(R)$ associated with polariton states
$|\mathcal{E}_{J}(R)\rangle$ (see their definition in Eq. 11) are presented in
Fig. 1e with the light-matter coupling strength $g_{\mathrm{c}}=0.001$ a.u.
and in Fig. 1f with the light-matter coupling strength $g_{\mathrm{c}}=0.005$.
These polariton potentials are color coded (as shown in the inset of panel
(e)) based on the expectation value of
$\langle\hat{a}^{\dagger}\hat{a}\rangle$ indicated on top of this panel. Note
that this should not be viewed as a “photon number” operator under the dipole
gauge used in the PF Hamiltonian[80, 56] because the rigorous photon number
operator should be obtained by applying the Power-Zienau-Woolley (PZW) Gauge
transformation[81, 82, 54] on the photon number operator
$\hat{a}^{\dagger}\hat{a}$. Nevertheless, it can be viewed as an approximate
estimation of the photon number when the light-matter couplings are not in the
ultra-strong coupling regime. [83]

The initial state (for $t=0$) of the molecule-cavity hybrid system is

$|\Phi(t=0)\rangle=|e,0\rangle\otimes|\chi\rangle,$ (45)

which corresponds to a Franck-Condon excitation of the hybrid system to the
$|e,0\rangle$ state, with $|\chi\rangle$ as the initial nuclear wavefunction.
For the SM model in this work, we use $\chi(R)=\langle
R|\chi\rangle\sim\exp[-M\omega_{0}(R-R_{0})^{2}/2\hbar]$, where $M$ is the
mass of the proton (nucleus in the SM model), $R_{0}$ is the position with a
minimum potential energy of the ground electronic state. Here, $\chi(R)$ is
the vibrational ground state wavefunction on the ground electronic states,
centered at $R_{0}$ under the harmonic approximation, with the harmonic
oscillation frequency being $\omega_{0}$. We use the parameters in the
original reference[49] for $R_{0}=-4$ and $\omega_{0}=0.000382$ a.u. To solve
the exact quantum dynamics, we use the DVR basis for the nuclear DOF and the
adiabatic-Fock state for the electronic-photonic subsystem. The details of the
exact quantum dynamics are provided in Appendix B.

###  III.2 Details of $\gamma$-SQC and spin-LSC Dynamics

To perform the $\gamma$-SQC dynamics, we need to sample the initial condition
for the quantum subsystem. In this work, we first sample the action-angle
variables $\\{\varepsilon_{j},\theta_{j}\\}$ then transform them to the
mapping variables$\\{p_{j},q_{j}\\}$ using Eq. 38. Among them, the action
variables $\\{\varepsilon_{j}\\}$ are sampled according to the window function
in Eq. 47, and the angle variables $\\{\theta_{j}\\}$ are randomly sampled
from $[0,2\pi)$. The triangle window is used in this work, although the square
window generates similar results.

For the spin-LSC dynamics, we use the focused initial conditions[37] as
described in section II.3, where the action variable $\varepsilon_{i}$ is set
to be $1+\Gamma/2$ for the initially occupied state and $\Gamma/2$ for the
initially unoccupied state, with $\Gamma$ expressed in Eq. 42. The angle
variables $\\{\theta_{j}\\}$ are randomly generated between $[0,2\pi)$ as in
the $\gamma$-SQC method. The canonical mapping variables are obtained from Eq.
38.

The initial nuclear distribution of all trajectory-based simulations
(Ehrenfest, FSSH, $\gamma$-SQC and spin-LSC) are generated by sampling the
Wigner density

$[\langle
R|\chi\rangle]_{\mathrm{w}}=\frac{1}{\hbar\pi}e^{-M(P^{2}+\omega^{2}(R-R_{0})^{2})/\omega\hbar},$
(46)

which is the Wigner transformation of the nuclear wavefunction
$\chi(R)=\langle R|\chi\rangle$ in the initial state (see Eq. 45). Here, $R$
and $P$ are the nuclear coordinate and momentum, respectively. The initial
state for the electronic-photonic subsystem is set to $|e0\rangle$. The
nuclear time step used in the QD-$\gamma$-SQC and QD-spin-LSC is $dt=0.1$ fs,
with $100$ equally spaced electronic time steps for the mapping variables’
integration during each nuclear time step. The equation of motion in Eq.
35-Eq. 36 are integrated using a second-order symplectic integrator for the
MMST variables[84, 85] for a given nuclear time step, and these mapping
variables are transformed based on Eq. 43 between two adjacent nuclear time
steps due to the change of the QD basis. The population dynamics using all MQC
and mapping methods were averaged over $5000$ trajectories, although $3000$
trajectories were enough to produce the basic trend of the polariton dynamics,
see Fig. S3 in the Supplementary Material. The light-matter coupling strength
$g_{\mathrm{c}}$ was chosen to be $0.001$ and $0.005$, according to our
previous work.[23]

We also benchmark the results of non-adiabatic mapping dynamics approaches
with commonly used MQC approaches, including the Ehrenfest dynamics and the
FSSH method. The details of these two MQC approaches are provided in Appendix
C. In particular, the Ehrenfest dynamics is equivalent to choosing $\gamma=0$
in the mapping theory (see Eq. 34) and an initial action-angle variables
condition (see Eq. 37 and Eq. 38) of $\varepsilon_{j}=\delta_{ij}$ (for the
initially occupied state $|i\rangle$) and $\theta_{j}=0$ (for all state
$|j\rangle$). One can thus use the same QD scheme and the mapping equation to
obtain the results of the Ehrenfest dynamics.[86]

Figure 2: The population dynamics of the adiabatic-Fock states in Shin-Metiu-
cavity model obtained from (a) Ehrenfest dynamics, (b) FSSH approach, (c)
$\gamma$-SQC method, and (d) spin-LSC dynamics. The population dynamics are
obtained with the approximate methods (open circles) and exact quantum
propagation (solid lines). Two electronic states and two Fock states are
considered in the simulation, and the light-matter coupling strength
$g_{\mathrm{c}}=0.001$ a.u.

##  IV Results

Fig. 2 presents the population dynamics of the adiabatic-Fock states simulated
using the approximate methods (open circles), including the MQC approaches
(Ehrenfest and FSSH) and the mapping dynamics methods ($\gamma$-SQC and spin-
LSC), compared to the numerically exact results (solid lines). The light-
matter coupling strength is chosen to be $g_{\mathrm{c}}=0.001$ a.u. The
system is initially prepared in the $|e0\rangle$ state and decays quickly into
the $|g1\rangle$ state during the first $\sim 12$ fs due to the large light-
matter coupling strength from the large transition dipole moment
($\bm{\mu}_{eg}$) between adiabatic electronic state $|g\rangle$ and
$|e\rangle$ (as shown in Fig. 1d). Then, the system starts to oscillate
between $|e0\rangle$ and $|g1\rangle$ until about $20$ fs. All of the MQC and
mapping dynamics methods can describe the above process reasonably well
compared to the exact results. After that, all the dynamics results (including
the exact one) show a fast population increase of the $|g0\rangle$ state,
which is due to the electronic NAC $d_{eg}$ that directly couples the
$|e0\rangle$ state to $|g0\rangle$ state (gold lines). All of the approximate
methods can qualitatively describe such a trend, but the MQC methods (panels
a-b) are less accurate compared to the mapping-based methods (panels c-d), in
terms of the rising of the $|g0\rangle$ population as well as its long time
plateau. Moreover, both Ehrenfest and FSSH dynamics predict a significant
population transfer from $|g1\rangle$ to $|e1\rangle$ state (panels a-b) as an
artifact that is not shown in the exact dynamics results. In contrast, the
$\gamma$-SQC and spin-LSC methods perform much better, where the population
transfer process from $|g1\rangle$ to $|e1\rangle$ state is largely suppressed
(panels c-d). Overall, the mapping methods outperform the MQC methods in this
small light-matter coupling case. It is worth mentioning that the population
dynamics results obtained with the FSSH method can be significantly improved
if one uses the proper estimator.[87] We have provided details of this
approach and numerical results in Appendix C. Even so, the FSSH method is
still facing many challenges from the improper treatment of the quantum
coherence and frustrated hop problems, which have been widely discussed.[27,
25, 88, 89]

Figure 3: The population dynamics of the adiabatic-Fock states in Shin-Metiu-
cavity model obtained from (a) Ehrenfest dynamics, (b) FSSH approach, (c)
$\gamma$-SQC method, and (d) spin-LSC dynamics. The population dynamics are
obtained with the approximate methods (open circles) and exact quantum
propagation (solid lines). Two electronic states and two Fock states are
considered in the simulation, and the light-matter coupling strength
$g_{\mathrm{c}}=0.005$ a.u.

Fig. 3 presents the polariton population dynamics with the coupling strength
$g_{\mathrm{c}}=0.005$ a.u. The oscillation between $|e0\rangle$ and
$|g1\rangle$ state population appears much earlier and faster compared to the
$g_{\mathrm{c}}=0.001$ results due to the larger light-matter coupling between
$|e0\rangle$ and $|g1\rangle$ states (see Eq. 5). Further, the $|g0\rangle$
and $|e1\rangle$ states are also getting populated at an earlier time, due to
the permanent dipole $\mu_{gg}$ and $\mu_{ee}$ that couples $|g1\rangle$ state
to $|g0\rangle$ state and $|e0\rangle$ state to $|e1\rangle$ state,
respectively. Similar to the $g_{\mathrm{c}}=0.001$ case, all the MQC and
mapping dynamics provide a reasonable accuracy for the population dynamics at
a short time, while the mapping methods perform much better than the MQC
methods at a longer time. In addition, the spin-LSC method outperforms
$\gamma$-SQC method in the description of $|g0\rangle$ and $|e0\rangle$ state
population after $t=20$ fs, as shown in Fig. 3c-d.

Until now, all of our simulations are restricted in the Hilbert subspace
formed by two electronic states ($|g\rangle$ and $|e\rangle$) and two photonic
Fock states ($|0\rangle$ and $|1\rangle$). The system could explore a larger
Hilbert space due to the increasing light-matter coupling strength. Thus, we
systematically check the polariton dynamics using the exact wavepacket
dynamics method with a larger number of electronic adiabatic states and Fock
states, as shown in Fig. S1 and S2 in the Supplementary Material. The results
show that, for the small light-matter coupling strength case
($g_{\mathrm{c}}=0.001$ a.u.), truncation to the Hilbert subspace formed by
two electronic states and two Fock states is enough to give an accurate
description of the population dynamics for the SM model studied in this work.
However, for the larger coupling strength case ($g_{\mathrm{c}}=0.005$ a.u.),
the polariton dynamics will converge when including four adiabatic electronic
states ($|g\rangle$, $|e\rangle$, $|f\rangle$, $|h\rangle$ with energies in
ascending order) and four Fock states ($|0\rangle$, $|1\rangle$, $|2\rangle$,
$|3\rangle$ with photon number in ascending order).

Figure 4: The population dynamics of the adiabatic-Fock states with (a,b)
Ehrenfest dynamics, (c,d) FSSH approach, (e,f) $\gamma$-SQC method, and (g,h)
spin-LSC dynamics. The population dynamics are obtained with the approximate
methods (open circles) and exact quantum propagation (solid lines). Four
adiabatic electronic states ($|g\rangle$, $|e\rangle$, $|f\rangle$,
$|h\rangle$ with energies in ascending order) and four Fock states
($|0\rangle$, $|1\rangle$, $|2\rangle$, $|3\rangle$ with photon number in
ascending order) are considered in the simulations and the light-matter
coupling strength is $g_{\mathrm{c}}=0.005$ a.u. Only the adiabatic-Fock
states with observable populations of more than $0.01$ are plotted.

To further test the performance of the mapping methods ($\gamma$-SQC and spin-
LSC) as well as the MQC methods (Ehrenfest and FSSH) in such a large Hilbert
subspace, which includes sixteen states formed by the tensor product of four
electronic states ($|g\rangle$, $|e\rangle$, $|f\rangle$, $|h\rangle$) and
four Fock states ($|0\rangle$, $|1\rangle$, $|2\rangle$, $|3\rangle$). Fig. 4
presents the results of using Ehrenfest dynamics (a-b), FSSH (c-d),
$\gamma$-SQC (e-f) and spin-LSC (g-h). Besides the adiabatic-Fock states
already appear in the four-states subspace ($|g0\rangle$, $|e0\rangle$,
$|g1\rangle$ and $|e1\rangle$), we can also see some other states
($|h0\rangle$, $|f0\rangle$, $|g2\rangle$ and $|f1\rangle$) are populated due
to the increasing light-matter coupling strength. All of the MQC (Ehrenfest,
FSSH) and mapping ($\gamma$-SQC and spin-LSC) dynamics results provide
accurate agreement with the exact one in the short time ($<20$ fs). After
that, all of the methods start to generate less accurate results (especially
for the $|g0\rangle$ population). Note that both $\gamma$-SQC and spin-LSC
perform less accurately compared to the situation in a smaller Hilbert
subspace (Fig. 3). This is because both $\gamma$-SQC and spin-LSC are
sensitive to the number of states of the system. For $\gamma$-SQC, more states
means less trajectory landed in the population action window.[90] For spin-
LSC, the ZPE correction $\Gamma$ (Eq. 42) explicitly depends on the number of
states $\mathcal{N}$. This suggests a need for the future development of more
accurate dynamics approaches. The current work, nevertheless, paves the way
for those future methods to be directly used for simulating on-the-fly
polariton quantum dynamics.

##  V Conclusions

In this work, we generalize the quasi-diabatic (QD) propagation scheme[44, 47,
45] to simulate the non-adiabatic polariton dynamics in molecule-cavity hybrid
systems. The adiabatic-Fock states, which are the tensor product states of the
adiabatic electronic states of the molecule and photon Fock states, are used
as the locally well-defined diabatic states for the dynamics propagation.[47,
45] These locally well-defined diabatic states allow using any diabatic
quantum dynamics methods for dynamics propagation, and the definition of these
states will be updated at every nuclear time step. The benefit of using such
adiabatic-Fock states is that one can conveniently obtain the electronic
adiabatic states energies, the nuclear gradient, the dipole moments and NACs
between these states, which are necessary ingredients in molecular cavity QED
simulations.

We use the recently developed non-adiabatic mapping dynamics approaches,
$\gamma$-SQC[35] and spin-LSC[37] to investigate polariton dynamics of a Shin-
Metiu model coupled to an optical cavity.[49, 23] To benchmark the results of
the obtained polariton dynamics, we performed simulations using the Ehrenfest
dynamics and the FSSH approaches, as well as the numerically exact polariton
wavepacket propagation. The results show that the mapping methods can
accurately describe the population dynamics of the molecule-cavity system both
at a short time and a longer time when compared to exact dynamics results. In
addition, the mapping methods outperform the Ehrenfest and FSSH approaches at
a long time dynamics. The numerical results also demonstrate that the
performance of the mapping methods ($\gamma$-SQC and spin-LSC) becomes less
accurate with an increased number of states in the simulation, indicating the
need for future theoretical development.

We envision that the theoretical development in this work will provide the
emerging polariton chemistry field with a general theoretical tool that
enables direct ab initio on-the-fly simulations of polariton photochemical
processes. We also anticipate that the theoretical developments in this work
will enable many recently developed diabatic quantum dynamics approaches to
directly simulate polariton quantum dynamics.

## Acknowledgments

This work was supported by the National Science Foundation CAREER Award under
Grant No. CHE-1845747 and by a Cottrell Scholar award (a program by Research
Corporation for Science Advancement). Computing resources were provided by the
Center for Integrated Research Computing (CIRC) at the University of
Rochester. D.H. appreciates valuable discussions with Wanghuai Zhou. We also
appreciate valuable suggestions by Prof. Yihan Shao.

## Conflict of Interest

The authors have no conflicts to disclose.

## Supplementary Material

See Supplementary Material for additional results of the convergence test for
the number of trajectories and the number of adiabatic electronic states and
Fock states.

## Availability of Data

The data that support the findings of this study are available from the
corresponding author upon a reasonable request.

##  Appendix A Details of the $\gamma$-SQC Approach

For practical reasons, the delta functions in Eq. 39 are broadened using well-
explored window functions, which can be used to bin the electronic action
variables in action-space. The triangle window[71, 35] is expressed as

$W_{j}({\bm{\varepsilon}})=w_{1}(\varepsilon_{j})\prod_{j^{\prime}\neq
j}^{\mathcal{N}}w_{0}(\varepsilon_{j},\varepsilon_{j^{\prime}}),$ (47)

where the window functions are defined as

$w_{1}(\varepsilon)=\begin{cases}(2-\varepsilon)^{2-\mathcal{N}},&1<\varepsilon<2\\\
0,&\mathrm{else}\end{cases}$ (48)

and

$w_{0}(\varepsilon,\varepsilon^{\prime})=\begin{cases}1,&\varepsilon^{\prime}<2-\varepsilon\\\
0,&\mathrm{else},\end{cases}$ (49)

and trajectories are assigned to state $j$ at time $t$ if $\varepsilon_{j}\geq
1$ and $\varepsilon_{j^{\prime}}<1$ for all $j^{\prime}\neq j$. The ZPE
parameter in the standard SQC method with triangle windowing is $\gamma=1/3$.

The time-dependent population of the state $|j\rangle$ is computed with Eq.
39. Using the window function estimator, the total population is no longer
properly normalized due to the fraction of trajectories that are outside of
any window region at any given time.[31] Thus, the total population must be
normalized[31] with the following procedure

${{\rho}_{jj}(t)}/{\sum_{i=1}^{N}{\rho}_{ii}(t)}\rightarrow{\rho}_{jj}(t).$
(50)

In the $\gamma$-SQC approach,[35] it was proposed that the mapping ZPE should
be chosen in such a way as to constrain the initial force to be composed
purely from the initially occupied state.[35] This new scheme has shown to
provide a significant improvement for photo-dissociation problems with coupled
Morse potentials[35] and has been combined with the kinematic momentum
approach[91] to carry out on-the-fly simulations of the methaniminium
cation.[70] The basic logic of $\gamma$-SQC is to choose an $\gamma_{j}$ for
each state $|j\rangle$ in every given individual trajectory, such that the
initial population is forced to respect the initial electronic excitation
focused onto a single excited state. If the initial electronic state is
$|i\rangle$, then

$\gamma_{j}=\varepsilon_{j}-\delta_{ji},$ (51)

or equivalently,

$\delta_{ji}=\varepsilon_{j}-\gamma_{j},$ (52)

where the $\\{\varepsilon_{j}\\}$ are uniformly sampled inside the window
function (Eq. 47), and then the $\gamma_{j}$ are chosen to satisfy Eq. 52.

These $\gamma_{j}$ will be explicitly used in the EOMs in Eqs. 35-36, and in
particular, the nuclear forces are now

${\bf F}=-{\frac{1}{2}}\sum_{ij}\nabla V_{ij}({\bf
R})\big{(}p_{i}p_{j}+q_{i}q_{j}-2\gamma_{j}\delta_{ji}\big{)},$ (53)

ensuring the initial forces (at $t=0$) are simply ${\bf F}=-\nabla V_{ii}({\bf
R})$. Previously, without any adjustments to $\gamma_{j}$, the chosen values
for $\gamma_{j}$ were only dependent on the windowing function itself, i.e.,
$\gamma_{j}=0.366$ for the square Windows and $\gamma_{j}=1/3$ for the
triangle windows. With the above $\gamma$-correction method,[35] each
individual trajectory will have its own state-specific $\gamma_{j}$ for state
$|j\rangle$ that is completely independent of the choice of window function.

##  Appendix B Details of the Adiabatic Electronic Calculation and the Exact
Quantum Dynamics Simulations

To calculate the electronic properties of the SM model, we use the Sinc DVR
basis[74] to represent the electronic adiabatic wavefunction and solve Eq. 3.
The grid of DVR is uniform with spacing $\Delta x=0.147$ in the range
$[-22,22]$ a.u. To test the convergence of grid points, we doubled the number
of grid points and the results were identical. The matrix elements of the
electronic Hamiltonian $\hat{H}_{\mathrm{el}}$ in this grid basis
$\\{|r_{i}\rangle\\}$ are given by

$\displaystyle\langle r_{i}|\hat{H}_{\mathrm{el}}|r_{j}\rangle$
$\displaystyle=\langle
r_{i}|\hat{T}_{r}+\hat{V}_{\mathrm{eN}}(\hat{r},R)+\hat{V}_{\mathrm{NN}}(R)|r_{j}\rangle$
(54) $\displaystyle=\langle
r_{i}|\hat{T}_{r}|r_{j}\rangle+[\hat{V}_{\mathrm{eN}}(r_{j},R)+\hat{V}_{\mathrm{NN}}(R)]\delta_{ij},$

where the $\langle r_{i}|\hat{T}_{r}|r_{j}\rangle$ is given analytically[74]
as follows

$\displaystyle\langle r_{i}|\hat{T}_{r}|r_{j}\rangle=$
$\displaystyle\frac{\hbar^{2}}{2m}\cdot\frac{\pi^{2}}{3(\Delta
r)^{2}}\Big{(}1+\frac{2}{N^{2}}\Big{)}\delta_{ij}$ (55)
$\displaystyle+\frac{\hbar^{2}}{2m}\cdot\frac{2(-1)^{j-i}\pi^{2}}{(\Delta
rN\sin(\frac{\pi(j-i)}{N}))^{2}}(1-\delta_{ij}).$

Directly diagonalizing the matrix of $\hat{H}_{\mathrm{el}}$ (in Eq. 54) at a
given nuclear (proton) position $R$ in this grid basis gives the accurate
adiabatic electronic states

$|\phi_{\alpha}(R)\rangle=\sum_{i}c^{\alpha}_{i}(R)|r_{i}\rangle,$ (56)

where $c^{\alpha}_{i}(R)$ is the expansion coefficient, which is purely real
for the adiabatic electronic states considered here.

A key ingredient for the QD propagation scheme is the overlap integral in Eq.
23a, Eq. 23b, and Eq. 32, which involve the overlaps between two adiabatic
states associated with two different reference geometries. These integrals are
conveniently calculated because all of the adiabatic states are represented
with the common DVR grids basis as follows

$\langle\phi_{\alpha}({R}_{0})|\phi_{\beta}({R}_{1})\rangle=\sum_{i,j}c^{\alpha}_{i}({R}_{0})\cdot
c^{\beta}_{j}({R}_{1})\langle r_{i}|r_{j}\rangle.$ (57)

Using these bases, the matrix elements for the dipole moment operator (Eq. 19)
are calculated as

$\displaystyle{\mu}_{\gamma\nu}({R})$
$\displaystyle=\sum_{ij}c^{\gamma}_{i}({R})\cdot c^{\nu}_{j}({R})\langle
r_{i}|(R-\hat{r})|r_{j}\rangle$ (58)
$\displaystyle=\sum_{ij}c^{\gamma}_{i}({R})\cdot
c^{\nu}_{j}({R})\cdot(R-r_{i})\delta_{ij}.$

Using $|\phi_{\alpha}({R})\rangle$ and $|\phi_{\beta}({R})\rangle$ in the grid
basis, we can directly evaluate the nuclear gradient
$\langle\phi_{\lambda}|\nabla\hat{H}_{\mathrm{en}}|\phi_{\nu}\rangle$ (Eq. 25)
as follows

$\langle\phi_{\lambda}({R})|\nabla\hat{H}_{\mathrm{en}}({R})|\phi_{\nu}({R})\rangle=\sum_{ij}c^{\lambda}_{i}({R})\cdot
c^{\nu}_{j}({R})\langle r_{i}|\nabla\hat{H}_{\mathrm{en}}({R})|r_{j}\rangle,$
(59)

where the $\nabla\hat{H}_{\mathrm{en}}({R})$ is evaluated analytically using
the expression of $\hat{H}_{\mathrm{en}}({R})$ in Eq. 44. This gives the
adiabatic gradient $\nabla E_{\lambda}$ (for $\lambda=\nu$) and derivative
coupling (for $\lambda\neq\nu$) as
$d_{\lambda\nu}=\langle\phi_{\lambda}|\nabla\hat{H}_{\mathrm{en}}|\phi_{\nu}\rangle/(E_{\nu}-E_{\lambda})$,
as indicated in Eq. 44. The nuclear gradient
$\langle\phi_{\lambda}|\nabla\hat{H}_{\mathrm{en}}|\phi_{\nu}\rangle$ is also
one of the key ingredient in the QD propagation.

Further, the nuclear gradient expression in Eq. II.2) for a polariton system
requires the derivative on the dipole matrix (Eq. 58). This requires the
evaluation of the derivative of the expansion coefficients $\nabla
c_{i}^{\gamma}(R)$ in Eq. 58. Instead of evaluating these derivatives, as been
commonly done in electronic structure calculations,[22] here, we evaluate this
derivative on dipole numerically as follows

$\nabla{\mu}_{\gamma\nu}(R)\approx\frac{{\mu}_{\gamma\nu}(R+\Delta
R)-{\mu}_{\gamma\nu}(R-\Delta R)}{2\Delta R}.$ (60)

To solve the exact quantum dynamics, we represent the total wavefunction of
the Hybrid system as
$|\Psi_{\xi}\rangle=\sum_{i,k}c^{\xi}_{ik}|\psi_{i}(R_{k})\rangle\otimes|R_{k}\rangle$,
where $\\{|R_{k}\rangle\\}$ is the DVR grid basis for the nucleus, and the
$|\psi_{i}(R_{k})\rangle=|\phi_{\alpha}(R_{k})\rangle\otimes|n\rangle$ is the
adiabatic-Fock basis (Eq. 8), where the electronic adiabatic basis
$|\phi_{\alpha}(R_{k})\rangle$ is obtained by solving the electronic
eigenequation (Eq. 3) using the DVR basis for the electronic DOF at the
nuclear configuration $R_{k}$. The coefficients for the total wavefunction
$c^{\xi}_{ik}$ and the eigenvalue of the total Hamiltonian $\hat{H}$ (Eq. 1)
will be obtained by solving the time-independent Schödinger equation
$\hat{H}|\Psi_{\xi}\rangle=E_{\xi}|\Psi_{\xi}\rangle$. We use the Sinc DVR
basis[74] for the nuclear DOF and solve the above eigenvalue problem to obtain
all the eigenvalues and eigenstates. We use finer grid points for nucleus
$\Delta x=0.016$ in the range $[-8,8]$. To test the convergence of grid
points, we doubled the number of grid points and the results were identical.
The time evolution dynamics is obtained by unitary evolution
$|\Phi(t)\rangle=\sum_{\xi}C_{\xi}\exp\left(-\frac{i}{\hbar}E_{\xi}t\right)|\Psi_{\xi}\rangle$,
where $C_{\xi}$ is the projection of initial total wavefunction onto the
$|\Psi_{\xi}\rangle$ as $C_{\xi}=\langle\Psi_{\xi}|\Phi(0)\rangle$, with the
initial wavefunction $|\Phi(0)\rangle$ expressed in Eq. 45. The details of the
exact polariton dynamics calculation can also be found in our recent work.[23]

##  Appendix C Details of the Ehrenfest and Surface Hopping Simulations

Besides the mapping dynamics methods ($\gamma$-SQC and spin-LSC), we also
apply the commonly used Ehrenfest and Tully’s FSSH [17, 18] algorithms to run
polariton dynamics. Details can be found in our previous work that develops
new gradient expressions for QED simulation with the MQC methods.[23] Here, we
briefly present these approaches for the completeness of this work.

In the Ehrenfest dynamics, the wavefunction of the quantum subsystem
(electronic-photonic DOFs) is written as

$|\Psi({\bf r};{\bf R}(t))\rangle=\sum_{i}c_{i}(t)|\psi_{i}({\bf
R}(t))\rangle,$ (61)

where $|\psi_{i}({\bf R}(t))\rangle=|\phi_{\alpha}({\bf
R}(t))\rangle\otimes|n\rangle$ is the adiabatic-Fock state basis (Eq. 8). The
quantum subsystem is described by the time-dependent Schödinger equation
(TDSE):

$i\hbar\frac{\partial}{\partial t}|\Psi({\bf r};{\bf
R}(t))\rangle=\hat{V}|\Psi({\bf r};{\bf R}(t))\rangle.$ (62)

The classical subsystem (nuclear DOF) is propagated using the Newton’s EOM,
where the nuclear force is evaluated from the time-dependent average potential
(mean field)

${\bf F}=-{\bf c}^{\dagger}[\nabla V]{\bf c},$ (63)

where ${\bf c}^{\dagger}$ is the transpose of the coefficient column vector
${\bf c}$ expressed as follows

${\bf c}^{\dagger}=(c_{1}(t),c_{2}(t),\ldots,c_{N}(t)).$ (64)

The nuclear gradient matrix is expressed as

$[\nabla V]\equiv\nabla[V]-[V][{\bf d}]+[{\bf d}][V],$ (65)

where $[V]$ and $[{\bf d}]$ are the matrix of $\hat{V}$ and derivative
coupling operator in the adiabatic-Fock state basis, respectively. The full
derivation of the gradient can be found in our previous work.[23]

In the FSSH dynamics, we expand the time-dependent wave function in the
polaritonic basis $|\mathcal{E}_{I}({\bf R}(t))\rangle$ (see definition in Eq.
11)

$|\Psi({\bf r};{\bf R}(t))\rangle=\sum_{I}c_{I}(t)|\mathcal{E}_{I}({\bf
R}(t))\rangle,$ (66)

where $c_{I}$ is the expansion coefficient, which will be used to compute the
fewest switching probability. Here, the nuclear force comes from only one
specific polariton state $|\mathcal{E}_{I}({\bf R}(t))\rangle$ (eigenstate of
$\hat{V}$, see Eq. 11) as follows

${\bf F}=-\nabla\mathcal{E}_{I},$ (67)

where $\mathcal{E}_{I}$ is the energy of the active adiabatic polariton state,
and $I$ is the active state index, which will be determined at every nuclear
propagation step. The nuclear gradient is

$\nabla\mathcal{E}_{I}=\sum_{jk}\langle\mathcal{E}_{I}|\psi_{j}\rangle\langle\psi_{j}|\nabla
V|\psi_{k}\rangle\langle\psi_{k}|\mathcal{E}_{I}\rangle,$ (68)

where the matrix element of $\nabla V$ is expressed in Eq. 65. The details of
the nuclear gradient in the polariton basis can be found in our previous
work.[23]

According to the “fewest switches” algorithm,[17] the probability of switching
(probability flux) from the active polariton state $|\mathcal{E}_{I}\rangle$
to any other polariton state $|\mathcal{E}_{J}\rangle$ during the time
interval between $t$ and $t+\delta t$ is

$f_{IJ}=-\frac{2\text{Re}[(\rho^{\mathrm{pl}}_{JI})^{*}\cdot{\bf{\dot{R}}}\cdot{\bf
d}_{JI}({\bf R})]\,\delta t}{\rho^{\mathrm{pl}}_{II}},$ (69)

where $\rho^{\mathrm{pl}}_{IJ}(t)$ is the reduced density matrix element in
the polariton basis expressed as follows

$\rho^{\mathrm{pl}}_{IJ}(t)=c_{I}(t)c_{J}^{*}(t).$ (70)

Since the probability should be positive definite, one sets[92] $f_{IJ}$ to 0
if $f_{IJ}<0$. The non-adiabatic transition, i.e. stochastic switch from the
currently occupied state $|\mathcal{E}_{I}\rangle$ to another state
$|\mathcal{E}_{K}\rangle$, occurs if the following condition is satisfied

$\sum_{J=1}^{K}{f_{IJ}}<\zeta<\sum_{J=1}^{K+1}{f_{IJ}},$ (71)

where $\zeta$ is a uniform randomly generated number between 0 and 1 at each
nuclear time step. If the transition is accepted, the active state is set to
the new adiabatic state $|\mathcal{E}_{K}\rangle$, while the velocities of the
nuclei are rescaled along the direction of the NAC vector ${\bf d}_{IK}({\bf
R})$ in order to conserve the total energy.[18] More details of performing
FSSH simulation of the polariton dynamics can be found in our previous
work.[23]

For the Ehrenfest dynamics and FSSH approach, we use the fourth-order Runge-
Kutta method to integrate the TDSE and the velocity Verlet algorithm to
integrate Newton’s EOM. The time step for the nuclear motion is $0.1$ a.u. and
the sub-step for solving the TDSE of the electronic-photonic subsystem is
$0.001$ a.u.. We have carefully checked that the total energy is well
conserved for all the trajectories. The initial condition is described by Eq.
45, where the nuclear DOF is sampled from the corresponding Wigner density
described in Eq. 46. For Ehrenfest dynamics, the initial coefficients
$c_{i}(0)$ for the state $|\psi_{i}\rangle=|e0\rangle$ is set to be one, and
the rest of the coefficients are set to be zero. These $\\{c_{j}(0)\\}$ can be
unitary transformed into the coefficients $\\{c_{I}(0)\\}$ for each nuclear
initial configuration described in Eq. 46. For the FSSH simulation, one needs
to choose an initial active state, and the initial state of the quantum
subsystem $|e(R),0\rangle$ is not one of the eigenstates
$|\mathcal{E}_{I}(R)\rangle$. We thus follow the previous work[93, 87] and use
the Monte-Carlo scheme to randomly choose the initial active state
$|\mathcal{E}_{I}(R)\rangle$ for each trajectory, based on the magnitude of
$|\langle\mathcal{E}_{I}(R)|e(R),0\rangle|^{2}$ for a given trajectory that
has the nuclear configuration at $R$ sampled from the Wigner density (Eq. 46).

When computing the population dynamics in a representation that is not the
adiabatic states of $\hat{V}$, there is no unique way to calculate them in the
FSSH approach.[87] In the main text, we present the populations of the
adiabatic-Fock states using the expansion coefficients $c_{I}(t)$ in Eq. 66.
There are, of course, alternative ways to compute populations of these
adiabatic-Fock states.[87] Below, we explore the alternative ways to compute
them.

For clarity, we denote the reduced density matrix in the adiabatic-Fock basis
as $\rho^{\mathrm{af}}_{ij}(t)$, and $\rho^{\mathrm{pl}}_{IJ}(t)$ is the
reduced density matrix in the polariton basis (expressed in Eq. 70). To get
the adiabatic-Fock state population of the $|\psi_{i}\rangle$ state
$\rho^{\mathrm{af}}_{ii}$ from the FSSH simulation, the most straightforward
way (as the results presented in the main text) is through following unitary
transformation

$[\rho^{\mathrm{af}}({\bf R}_{l}(t))]=\mathbf{U}[\rho^{\mathrm{pl}}({\bf
R}_{l}(t))]\mathbf{U}^{\dagger},$ (72)

where $[\rho^{\mathrm{pl}}{\bf R}_{l}(t))]$ is the reduced density matrix in
the polariton basis along a given nuclear trajectory ${\bf R}_{l}(t)$, with
$l$ as the label of the trajectory, and the elements as

$\rho^{\mathrm{pl}}_{IJ}({\bf R}_{l}(t))=c_{I}(t)c_{J}^{*}(t).$ (73)

Further, $\mathbf{U}({\bf R}_{l}(t))$ is the matrix that diagonalize the
matrix $[V({\mathbf{R}}_{l}(t))]$ as shown in Eq. 13, along the same
trajectory ${\bf R}_{l}(t)$. The adiabatic-Fock state population is then
obtained from trajectory average as follows

$P_{i}(t)=\frac{1}{N}\sum_{l}^{N}\left[\mathbf{U}[\rho^{\mathrm{pl}}({\bf
R}_{l}(t))]\mathbf{U}^{\dagger}\right]_{ii},$ (74)

where $N$ is the total number of the trajectories. This is the estimator used
in the FSSH calculation presented in the main text.

For FSSH, there are two other commonly used choices[87] to calculate the
populations that are not in an adiabatic representation. These methods vary on
how to calculate the polaritonic state density matrix
$[\rho^{\mathrm{pl}}({\bf R}_{l}(t))]$. The first choice is based on the
active state index and ignores the polaritonic state coefficients
$\\{c_{I}(t)\\}$ and the density matrix elements are written as

$\rho^{\mathrm{pl}}_{IJ}({\bf R}_{l}(t))=\begin{cases}\delta_{IK},&I=J\\\
0,&I\neq J,\end{cases}$ (75)

where $K$ is the active polaritonic state. This method explicitly assumes that
the off-diagonal elements of the polaritonic state density matrix are zero,
which is often not a good one. The adiabatic-Fock state population is then
obtained from the same transformation described in Eq. 72, and the ensemble
average over all trajectories is computed as described in Eq. 74.

The other choice[87] (motivated by the mixed quantum-classical Liouville
approach) is to calculate the diagonal elements of $\rho^{\mathrm{pl}}$ using
the active state index, and calculate the off-diagonal elements using the
polaritonic state expansion coefficients $\\{c_{I}(t)\\}$

$\rho^{\mathrm{pl}}_{IJ}({\bf R}_{l}(t))=\begin{cases}\delta_{IK},&I=J\\\
c_{I}c_{J}^{*},&I\neq J,\end{cases}$ (76)

where $K$ is the active index. The adiabatic-Fock state population is then
obtained from the same transformation described in Eq. 72, and the ensemble
average over all trajectories as described in Eq. 74.

Following the same notation as used in Ref. 87, we refer to the choice in Eq.
75 as Method 1, the choice in Eq. 73 as Method 2 (same as the FSSH results
presented in the main text), and the choice in Eq. 76 as Method 3. The FSSH
dynamics results based on these three methods are presented in Fig. 5. We can
see that Method 3 performs much better than Method 1 and Method 2, consisting
with the conclusion in Ref. 87.

Figure 5: Population dynamics of the adiabatic-Fock states obtained from the
FSSH method (open circles) and the exact quantum dynamics propagation (solid
lines), using different population estimators. (a,b) Method 1; (c,d) Method 2;
(e,f) Method 3. The results obtained using method 2 are presented in the main
text of this work. The light-matter coupling strength is set to be (a,c,e)
$g_{\mathrm{c}}=0.001$ a.u. and (b,d,f) $g_{\mathrm{c}}=0.005$ a.u.

## References

  * Ebbesen [2016] T. W. Ebbesen, Acc. Chem. Res. 49, 2403 (2016).
  * Kowalewski and Mukamel [2017] M. Kowalewski and S. Mukamel, Proc. Natl. Acad. Sci. U.S.A. 114, 3278 (2017).
  * Flick  _et al._ [2017a] J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, Proc. Natl Acad. Sci. USA 114, 3026 (2017a).
  * Ribeiro  _et al._ [2018] R. F. Ribeiro, L. A. Martínez-Martínez, M. Du, J. Camos-Gonzalez-Angulo, and J. Yuen-Zhou, Chem. Sci. 9, 6325 (2018).
  * Feist  _et al._ [2018] J. Feist, J. Galego, and F. J. Garcia-Vidal, ACS Photonics 5, 205 (2018).
  * Mandal and Huo [2019] A. Mandal and P. Huo, J. Phys. Chem. Lett. 10, 5519 (2019).
  * Hutchison  _et al._ [2012] J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux, and T. W. Ebbesen, Angew. Chem. Int. Ed. 51, 1592 (2012).
  * Thomas  _et al._ [2019] A. Thomas, L. Lethuillier-Karl, K. Nagarajan, R. M. A. Vergauwe, J. George, T. Chervy, A. Shalabney, E. Devaux, C. Genet, J. Moran, and T. W. Ebbesen, Science 363, 615 (2019).
  * Mandal  _et al._ [2020a] A. Mandal, T. D. Krauss, and P. Huo, J. Phys. Chem. B 124, 6321 (2020a).
  * Luk  _et al._ [2017] H. L. Luk, J. Feist, J. J. Toppari, and G. Groenhof, J. Chem. Theory Comput. 13, 4324 (2017).
  * Groenhof and Toppari [2018] G. Groenhof and J. J. Toppari, J. Phys. Chem. Lett. 9, 4848 (2018).
  * Groenhof  _et al._ [2019] G. Groenhof, C. Climent, J. Feist, D. Morozov, and J. J. Toppari, J. Phys. Chem. Lett. 10, 5476 (2019).
  * Tichauer  _et al._ [2021] R. H. Tichauer, J. Feist, and G. Groenhof, J. Chem. Phys. 154, 104112 (2021).
  * Tully [2012a] J. C. Tully, J. Chem. Phys. 137, 22A301 (2012a).
  * Barbatti [2011] M. Barbatti, Wiley Int. Rev. Comp. Mol. Sci. 1, 620 (2011).
  * Mai  _et al._ [2015] S. Mai, P. Marquetand, and L. Gonzalez, Int. J. Quantum Chem. 115, 1215 (2015).
  * Tully [1990] J. C. Tully, J. Chem. Phys. 93, 1061 (1990).
  * Hammes-Schiffer and Tully [1994] S. Hammes-Schiffer and J. C. Tully, J. Chem. Phys. 12, 4657 (1994).
  * Fregoni  _et al._ [2018] J. Fregoni, G. Granucci, E. Coccia, M. Persico, and S. Corni, Nat. Commun. 9, 4688 (2018).
  * Fregoni  _et al._ [2020a] J. Fregoni, S. Corni, M. Persico, and G. Granucci, J. Comput. Chem. 41, 2033 (2020a).
  * Fregoni  _et al._ [2020b] J. Fregoni, G. Granucci, M. Persico, and S. Corni, Chem 6, 250 (2020b).
  * Zhang  _et al._ [2019] Y. Zhang, T. Nelson, and S. Tretiak, J. Chem. Phys. 151, 154109 (2019).
  * Zhou  _et al._ [2022] W. Zhou, D. Hu, A. Mandal, and P. Huo, J. Chem. Phys. 157, 104118 (2022).
  * Parandekar and Tully [2006] P. V. Parandekar and J. C. Tully, J. Chem. Theory Comput. 2, 229 (2006).
  * Subotnik  _et al._ [2016] J. E. Subotnik, A. Jain, B. Landry, A. Petit, W. Ouyang, and N. Bellonzi, Annu. Rev. Phys. Chem. 67, 387 (2016).
  * Meyer and Miller [1979] H. Meyer and W. H. Miller, J. Chem. Phys. 70, 3214 (1979).
  * Stock and Thoss [1997] G. Stock and M. Thoss, Phys. Rev. Lett. 78, 578 (1997).
  * Thoss and Stock [1999] M. Thoss and G. Stock, Phys. Rev. A 59, 64 (1999).
  * Huo and Coker [2011] P. Huo and D. F. Coker, J. Chem. Phys. 135, 201101 (2011).
  * Huo and Coker [2012] P. Huo and D. F. Coker, Mol. Phys. 110, 1035 (2012).
  * Cotton and Miller [2013a] S. J. Cotton and W. H. Miller, J. Chem. Phys. 139, 234112 (2013a).
  * Cotton and Miller [2013b] S. J. Cotton and W. H. Miller, J. Phys. Chem. A 117, 7190 (2013b).
  * Hsieh and Kapral [2012] C.-Y. Hsieh and R. Kapral, J. Chem. Phys. 137, 22A507 (2012).
  * Hsieh and Kapral [2013] C.-Y. Hsieh and R. Kapral, J. Chem. Phys. 138, 134110 (2013).
  * Cotton and Miller [2019a] S. J. Cotton and W. H. Miller, J. Chem. Phys. 150, 194110 (2019a).
  * Runeson and Richardson [2019] J. E. Runeson and J. O. Richardson, J. Comp. Phys. 151, 044119 (2019).
  * Runeson and Richardson [2020] J. E. Runeson and J. O. Richardson, J. Comp. Phys. 152, 084110 (2020).
  * Bossion  _et al._ [2022] D. Bossion, W. Ying, S. Chowdhury, and P. Huo, J. Chem. Phys. 157, 084105 (2022).
  * Sun  _et al._ [1998] X. Sun, H. Wang, and W. H. Miller, J. Chem. Phys. 109, 7064 (1998).
  * Shi and Geva [2004] Q. Shi and E. Geva, J. Phys. Chem. A. 108 (2004).
  * Weight  _et al._ [2021] B. M. Weight, A. Mandal, and P. Huo, J. Chem. Phys. 155, 084106 (2021).
  * Chowdhury  _et al._ [2021] S. N. Chowdhury, A. Mandal, and P. Huo, J. Chem. Phys. 154, 044109 (2021).
  * Mandal  _et al._ [2018a] A. Mandal, S. S. Yamijala, and P. Huo, J. Chem. Theory Comput. 14, 1828 (2018a).
  * Mandal  _et al._ [2018b] A. Mandal, F. A. Shakib, and P. Huo, J. Chem. Phys. 148, 244102 (2018b).
  * Mandal  _et al._ [2019] A. Mandal, J. S. Sandoval C., F. A. Shakib, and P. Huo, J. Phys. Chem. A 123, 2470 (2019).
  * Sandoval  _et al._ [2018] J. S. Sandoval, A. Mandal, and P. Huo, J. Chem. Phys. 149, 044115 (2018).
  * Zhou  _et al._ [2019] W. Zhou, A. Mandal, and P. Huo, J. Phys. Chem. Lett. 10, 7062 (2019).
  * Shin and Metiu [1995] S. Shin and H. Metiu, J. Chem. Phys. 102, 9285 (1995).
  * Hoffmann  _et al._ [2020] N. M. Hoffmann, L. Lacombe, A. Rubio, and N. T. Maitra, J. Chem. Phys. 153, 104103 (2020).
  * Doltsinis and Marx [2002] N. L. Doltsinis and D. Marx, J. Chem. Theory Comput. 1, 319 (2002).
  * Schäfer  _et al._ [2018] C. Schäfer, M. Ruggenthaler, and A. Rubio, Phys. Rev. A 98, 043801 (2018).
  * Marx and Hutter [2009] D. Marx and J. Hutter,  _Ab initio Molecular Dynamics: Basic Theory and Advanced Methods_ (Cambridge University Press, 2009) p. 11.
  * Kowalewski  _et al._ [2016] M. Kowalewski, K. Bennett, and S. Mukamel, J. Phys. Chem. Lett. 7, 2050 (2016).
  * Taylor  _et al._ [2020] M. A. D. Taylor, A. Mandal, W. Zhou, and P. Huo, Phys. Rev. Lett. 125, 123602 (2020).
  * Rokaj  _et al._ [2018] V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and A. Rubio, J. Phys. B: At. Mol. Opt. Phys. 51, 034005 (2018).
  * Schäfer  _et al._ [2020] C. Schäfer, M. Ruggenthaler, V. Rokaj, and A. Rubio, ACS Photonics 7, 975 (2020).
  * Koessler  _et al._ [2022] E. R. Koessler, A. Mandal, and P. Huo, J. Chem. Phys. 157, 064101 (2022).
  * Flick  _et al._ [2017b] J. Flick, H. Appel, M. Ruggenthaler, and A. Rubio, J. Chem. Theory Comput. 13, 1616 (2017b).
  * Webster  _et al._ [1991] F. Webster, P.J.Rossky, and R.A.Friesner, Comput. Phys. Commun. 63, 494 (1991).
  * Meek and Levine [2014] G. A. Meek and B. G. Levine, J. Phys. Chem. Lett. 5, 2351 (2014).
  * Jain  _et al._ [2016] A. Jain, E. Alguire, and J. E. Subotnik, J. Chem. Theory Comput. 12, 5256 (2016).
  * Joubert-Doriol and Izmaylov [2018] L. Joubert-Doriol and A. F. Izmaylov, J. Chem. Phys. 148, 114102 (2018).
  * Ben-Nun and Martínez [1998] M. Ben-Nun and T. J. Martínez, J. Chem. Phys. 108, 7244 (1998).
  * Saita and Shalashilin [2012] K. Saita and D. V. Shalashilin, J. Chem. Phys. 137, 22A506 (2012).
  * Granucci  _et al._ [2001] G. Granucci, M. Persico, and A. Toniolo, J. Chem. Phys. 114, 10608 (2001).
  * Plasser  _et al._ [2012] F. Plasser, G. Granucci, J. Pittner, M. Barbatti, M. Persico, and H. Lischka, J. Chem. Phys. 137, 22A514 (2012).
  * Löwdin [1950] P. O. Löwdin, J. Chem. Phys. 18, 365 (1950).
  * Miller and Cotton [2016] W. H. Miller and S. J. Cotton, Faraday Discuss. 195, 9 (2016).
  * Müller and Stock [1999] U. Müller and G. Stock, J. Chem. Phys. 111, 77 (1999).
  * Hu  _et al._ [2021] D. Hu, Y. Xie, J. Peng, and Z. Lan, J. Chem. Theory Comput. 17, 3267 (2021).
  * Cotton and Miller [2016] S. J. Cotton and W. H. Miller, J. Chem. Phys. 145, 144108 (2016).
  * Lacombe  _et al._ [2019] L. Lacombe, N. M. Hoffmann, and N. T. Maitra, Phys. Rev. Lett. 123, 083201 (2019).
  * Martinez  _et al._ [2021] P. Martinez, B. Rosenzweig, N. M. Hoffmann, L. Lacombe, and N. T. Maitra., J. Chem. Phys. 154, 014102 (2021).
  * Colbert and Miller [1992] D. T. Colbert and W. H. Miller, J. Chem. Phys. 96, 1982 (1992).
  * Farag  _et al._ [2021] M. H. Farag, A. Mandal, and P. Huo, Phys. Chem. Chem. Phys. 23, 16868 (2021).
  * Kockum  _et al._ [2019] A. F. Kockum, A. Miranowicz, S. D. Liberato, S. Savasta, and F. Nori, Nature Rev. Phys. 1, 19 (2019).
  * Stefano  _et al._ [2019] O. D. Stefano, A. Settineri, V. Macri, L. Garziano, R. Stassi, S. Savasta, and F. Nori, Nature Phys. 15, 803 (2019).
  * Baranov  _et al._ [2020] D. G. Baranov, B. Munkhbat, E. Zhukova, A. Bisht, A. Canales, B. Rousseaux, G. Johansson, T. J. Antosiewicz, and T. Shegai, Nat. Commun. 11, 1 (2020).
  * Santhosh  _et al._ [2016] K. Santhosh, O. Bitton, L. Chuntonov, and G. Haran, Nat. Comm. 7, 11823 (2016).
  * Mandal  _et al._ [2020b] A. Mandal, S. M. Vega, and P. Huo, J. Phys. Chem. Lett. 11, 9215 (2020b).
  * Power and Zienau [1959] E. A. Power and S. Zienau, Philos. Trans. R. Soc. London, Ser. A 251, 427 (1959).
  * Cohen-Tannoudji  _et al._ [1989] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, John Wiley & Sons, Inc. (1989).
  * Forn-Díaz  _et al._ [2019] P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, Rev. Mod. Phys. 91, 025005 (2019).
  * Kelly  _et al._ [2012] A. Kelly, R. van Zon, J. Schofield, and R. Kapral, J. Chem. Phys. 136, 084101 (2012).
  * Church  _et al._ [2018] M. S. Church, T. J. H. Hele, G. S. Ezra, and N. Ananth, J. Chem. Phys. 148, 102326 (2018).
  * Bellonzi  _et al._ [2016] N. Bellonzi, A. Jain, and J. E. Subotnik, J. Chem. Phys. 144, 154110 (2016).
  * Landry  _et al._ [2013] B. R. Landry, M. J. Falk, and J. E. Subotnik, J. Chem. Phys. 139, 211101 (2013).
  * Wang  _et al._ [2016] L. Wang, A. Akimov, and O. V. Prezhdo, J. Phys. Chem. Lett. 7, 2100 (2016).
  * Granucci and Persico [2007] G. Granucci and M. Persico, J. Chem. Phys. 126, 134114 (2007).
  * Cotton and Miller [2019b] S. J. Cotton and W. H. Miller, J. Chem. Phys. 150, 104101 (2019b).
  * Cotton  _et al._ [2017] S. J. Cotton, R. Liang, and W. H. Miller, J. Chem. Phys. 147, 064112 (2017).
  * Tully [2012b] J. C. Tully, J. Chem. Phys. 137, 22A301 (2012b).
  * Hazra  _et al._ [2010] A. Hazra, A. V. Soudackov, and S. Hammes-Schiffer, J. Phys. Chem. B 114, 12319 (2010).