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–6} which have been shown to facilitate new chemical reactivities.^{1,6–9} Theoretical investigations play a crucial role in understanding the fundamental limit and basic principles in this emerging field,^{5,6,10–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–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–13,19–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), are valid for 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 that are based on the Meyer–Miller–Stock–Thoss (MMST) mapping formalism.^{26–28} These methods include partial linearized density matrix^{29,30} (PLDM), symmetrical quasi-classical^{31,32} (SQC), and the quantum-classical Liouville equation (QCLE) dynamics.^{33,34} In particular, the recently developed *γ*-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 (spin-LSC) approach,^{36–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 *γ*-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,9,42}

Recently, we have developed the quasi-diabatic (QD) propagation scheme^{41,43–47} 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 *γ*-SQC^{35} and spin-LSC^{36,37} to obtain non-adiabatic polariton dynamics, which outperforms widely used MQC approaches.

## II. THEORY AND METHODS

### A. 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

where $T\u0302n$ represents the nuclear kinetic energy operator, and $H\u0302en$ is the electronic Hamiltonian that describes electron–nucleus interactions. Furthermore, $H\u0302p$, $H\u0302enp$, and $H\u0302d$ 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 $H\u0302en$, which describes the common molecular Hamiltonian excluding the nuclear kinetic energy, is described as follows:

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–52} Modern electronic structure theory have been developed around solving the eigenvalue problem of $H\u0302en$, providing the following electronically adiabatic energy and its corresponding state:

Here, |*ϕ*_{α}(**R**)⟩ represents the *α*_{th} many-electron adiabatic state for a given molecular system, with the adiabatic energy *E*_{α}(**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

where $q\u0302c=\u210f/2\omega c(a\u0302\u2020+a\u0302)$ and $p\u0302c=i\u210f\omega c/2(a\u0302\u2020\u2212a\u0302)$ are photon field operators, $a\u0302\u2020$ and $a\u0302$ are the photonic creation and annihilation operators, respectively, and *ω*_{c} is the photon frequency.

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

where ** λ** =

*λ*·

**characterizes the cavity photon field strength, and**

*ϵ***is the direction of the field polarization. The photon field strength is determined by the volume of the cavity as $\lambda =1/\u03f50V0$, where**

*ϵ**ϵ*

_{0}is the permittivity inside the cavity and $V0$ is the effective quantization volume inside the cavity. Another way to characterize the light–matter coupling strength is using $gc=\u210f\omega c/2\lambda $. Note that the common notation used in the literature,

^{6,53}the definition of

*g*

_{c}also includes $\lambda \u22c5\mu \u0302$. Furthermore, the total dipole operator of both electrons and nuclei is defined as

where −*e* is the charge of the electron and *Z*_{j}*e* is the charge of the *j*_{th} nucleus.

Finally, the DSE term is expressed as

This is a necessary term in the PF Hamiltonian and ensures both gauge invariance of the Hamiltonian^{9,54} and a bounded ground state.^{9,55,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

where quantum number *i* ≡ {*α*, *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 the 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

where $T\u0302n$ represents the nuclear kinetic energy operator, and $V\u0302$ represents the rest of the Hamiltonian. For a bare molecular system, $V\u0302=H\u0302en$ expressed in Eq. (2). For a molecule–cavity hybrid system,

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

where $|EJ(R)\u3009$ is the polariton state with polariton energy $EJ(R)$. The polariton eigenstate can be expressed as

where $c\alpha ,nJ(R)=\u27e8\varphi \alpha (R),n|EJ(R)\u27e9$ and $EJ(R)$ can be obtained by diagonalizing the matrix of $V\u0302=H\u0302pl$ [constructed from the adiabatic-Fock state basis in Eq. (8)] as

where

Note that the **R**-dependence of $|EJ(R)\u3009$ is entirely coming from the **R**-dependence of the adiabatic states |*ϕ*_{α}(**R**)⟩, and the Fock state |*n*⟩ is completely **R**-independent. Meanwhile, the $|EJ(R)\u3009$ is the eigenstate of $V\u0302$, whereas the adiabatic state |*ϕ*_{α}(**R**)⟩ is only the eigenstate of $H\u0302en$, and not for $V\u0302$.

### B. 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.^{41,43–47}

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 *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 |*ψ*_{i}(**R**)⟩ = |*ϕ*_{α}(**R**), *n*⟩ [Eq. (8)], which is not the same as the polariton states $|EJ(R)\u3009$ [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* ∈ [*t*_{0}, *t*_{1}], where the nuclear positions evolve from **R**(*t*_{0}) to **R**(*t*_{1}), and the corresponding adiabatic-Fock basis [defined in Eq. (8)] are {|*ψ*_{i}(**R**(*t*_{0}))⟩} and ${|\psi jR(t1)\u3009}$. We uses the basis {|*ψ*_{i}(**R**_{0})⟩ ≡|*ϕ*_{α}(**R**_{0}), *n*⟩} at the reference nuclear geometry **R**(*t*_{0}) as the *diabatic* basis during this short-time propagation such that

With the above QD basis defined independently of **R**(*t*) within each propagation segment, the electronic derivative couplings vanish while $V\u0302(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* ∈ [*t*_{0}, *t*_{1}].

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

For on-the-fly simulations, this quantity is obtained from a linear interpolation^{59} between *V*_{αβ,mn}(**R**_{0}) and *V*_{αβ,mn}(**R**(*t*_{1})) as follows:

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*_{αβ,mn}(**R**_{0}) and *V*_{αβ,mn}(**R**(*t*_{1})) separately for the molecule–cavity hybrid system, as discussed below. Using electronic *ab initio* calculation, as well as the properties of $a\u0302\u2020$ and $a\u0302$ for the photonic DOF, we can explicitly evaluate each term of *V*_{αβ,mn}(**R**_{0}) [see Eq. (10)] as follows:

where $H\u0302p$ [see its definition in Eq. (4)] is an **R**-independent operator, the sum *∑*_{γ} in the matrix element of $H\u0302d$ runs over the diabatic states, and $D\alpha \beta 2$ denotes the elements of DSE. Furthermore, the matrix element of the dipole operator under the diabatic representation is expressed as

Similarly, at time *t*_{1}, the matrix element $V\alpha \beta ,mn(R(t1))=\u27e8\varphi \alpha (R0),m|V\u0302(R(t1))|\varphi \beta (R0),n\u27e9$ can also be written explicitly, with each term expressed as follows:

where $H\alpha \beta en(R(t1))\u2261\u27e8\varphi \alpha (R0)|H\u0302en(R(t1))|\varphi \beta (R0)\u27e9$, and $\mu \alpha \beta (R(t1))\u2261\u27e8\varphi \alpha (R0)|\mu \u0302(R(t1))|\varphi \beta (R0)\u27e9$.

To conveniently calculate $H\alpha \beta en(R(t1))$ and *μ*_{αβ}(**R**(*t*_{1})), we use the following relations:

where the matrix elements at **R**(*t*_{1}) are expressed as

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

Using the above information, as well as Eq. (17), we can obtain each term of *V*_{αβ,mn}(**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 ∇*V*_{αβ,mn}(**R**(*t*_{1})). First, let us focus on the gradient term from the electron–nuclear coupling term^{47} as follows:

Here, from the first to the second line, we have used the fact that neither ⟨*ϕ*_{α}(**R**_{0})| nor ⟨*m*| are **R**-dependent, which allows ∇ to bypass both and directly act on $V\u0302(R(t1))$. We have also used the fact that $H\u0302en$ does not contain any photonic operators. The gradient term $\u27e8\varphi \lambda (R(t1))|\u2207H\u0302en(R(t1))|\varphi \nu (R(t1))\u27e9$ can be evaluated using the following well-known equality:^{18}

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

For the gradient on the matrix $H\alpha \beta ,mnp$, because there is no nuclear DOF in $H\u0302p$, thus

For the gradient on the light–matter interaction term $H\alpha \beta ,mnenp$, we have

where *S*_{αλ} and $S\beta \nu \u2020$ are defined in Eqs. (23a) and (23b), respectively. To evaluate the term $\u27e8\varphi \lambda (R)|\u2207\mu \u0302(R)|\varphi \nu (R)\u27e9$ that appears in Eq. (28), we use a simple relation based on the chain rule as follows:

where **d**_{λγ}(**R**) = ⟨*ϕ*_{λ}(**R**)|∇*ϕ*_{γ}(**R**)⟩ = −⟨∇*ϕ*_{λ}(**R**)|*ϕ*_{γ}(**R**)⟩ and **d**_{γν}(**R**) = ⟨*ϕ*_{γ}(**R**)|∇*ϕ*_{ν}(**R**)⟩ = −⟨∇*ϕ*_{γ}(**R**)|*ϕ*_{ν}(**R**)⟩ are the electronic derivative couplings [defined in Eq. (26)]. Note that from the first line to the second line, we have inserted $P\u0302=\u2211\gamma |\varphi \gamma (R)\u3009\u3008\varphi \gamma (R)|$, which is the resolution of identity in the electronic subspace. As one can clearly see, this term requires the evaluation of derivative coupling **d**_{λγ}(**R**) and **d**_{γν}(**R**), and the derivative on dipole matrix element ∇*μ*_{λν}(**R**). For most of the electronic structure methods, the derivatives on dipole matrix elements ∇*μ*_{λν}(**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:

where the term $\u27e8\varphi \alpha (R(t0))|\mu \u0302(R(t1))|\varphi \gamma (R(t0))\u27e9$ can be evaluated using Eq. (21b), and the $\u27e8\varphi \alpha (R0)|\u2207\mu \u0302(R(t1))|\varphi \beta (R0)\u27e9$ type of derivative can be computed in the same fashion as elaborated in Eq. (28).

Using the matrix elements *V*_{αβ,mn}(**R**(*t*)) [Eq. (17)] and the nuclear gradient ∇*V*_{αβ,mn}(**R**) [as outlined in Eqs. (24)–(30)], one can, in principle, use any trajectory-based approaches or wavepacket approaches with guiding trajectories^{62–64} to propagate the quantum dynamics in the time step *dt* = *t*_{1} − *t*_{0} for *t* ∈ [*t*_{0}, *t*_{1}]. During the next short-time propagation segment *t* ∈ [*t*_{1}, *t*_{2}], the QD scheme adopts a new reference geometry $R0\u2032\u2261R(t1)$ and new *diabatic* basis $|\psi j(R0\u2032)\u3009\u2261|\psi j(R(t1))\u3009=|\varphi \beta (R(t1)),m\u3009$. Between the *t* ∈ [*t*_{0}, *t*_{1}] propagation and the *t* ∈ [*t*_{1}, *t*_{2}] propagation segments, all of the necessary quantities will be transformed from {|*ψ*_{i}(**R**_{0})⟩} to the ${|\psi j(R0\u2032)\u3009}$ basis using the relation

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

where the electronic adiabatic state overlaps ⟨*ϕ*_{α}(**R**(*t*_{0}))|*ϕ*_{β}(**R**(*t*_{1}))⟩ are directly calculated using the discrete variable representation (DVR) basis, and the Fock states are orthonormal to each other ⟨*n*|*m*⟩ = *δ*_{n,m}. For *ab initio* on-the-fly simulations, these overlap matrices can be computed based on the approach outlined in Ref. 65, which we have extensively used in our previous work on QD based on-the-fly simulation with the CASSCF method.^{41,47}

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.^{66,67} Nevertheless, the orthogonality remains to be well satisfied among {|*ψ*_{i}(**R**(*t*_{0}))⟩} or {|*ψ*_{j}(**R**(*t*_{1}))⟩}. 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.^{66,67} This potential issue can be easily resolved by using orthonormalization procedure among the vectors of the overlap matrix composed by ⟨*ψ*_{i}(**R**(*t*_{0}))|*ψ*_{j}(**R**(*t*_{1}))⟩, as been done in our previous work^{44} for simulating photo-induced charge transfer dynamics. Here, we perform the Löwdin orthogonalization procedure^{68} as commonly used in the local diabatization approach^{66} to ensure unitary propagation.

Furthermore, the QD scheme ensures a stable propagation of the quantum dynamics compared to directly solving it in the adiabatic representation. This is due to the fact that directly solving electronic dynamics in the adiabatic state requires the non-adiabatic coupling $\u27e8\varphi \beta (R(t))|\u2202\u2202t\varphi \alpha (R(t))\u27e9=d\beta \alpha (R)\u22c5R\u0307$, which might exhibit highly peaked values and cause large numerical errors. The QD scheme explicitly alleviates this difficulty by using the well behaved transformation matrix elements ⟨*ϕ*_{β}(**R**(*t*_{1}))|*ϕ*_{α}(**R**(*t*_{2}))⟩ instead of $\u27e8\varphi \beta (R(t))|\u2202\u2202t\varphi \alpha (R(t))\u27e9$. A detailed numerical example can be found in Fig. 4 of Ref. 46. In this work, the time step *dt* is a fixed constant. We have carefully checked the convergence of the dynamics at the single trajectory level with *dt*. As explained in our previous work,^{46} the QD propagation scheme makes the dynamics less sensitive due to the well defined quantities ⟨*ϕ*_{β}(**R**(*t*_{1}))|*ϕ*_{α}(**R**(*t*_{2}))⟩.

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 {|*ψ*_{i}(**R**_{0})⟩} 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 **R**_{0} is close to the nuclear geometry **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.

### C. Non-adiabatic mapping dynamics methods

The Meyer–Miller–Stock–Thoss (MMST) formalism^{26–28} maps the discrete electronic DOFs onto continuous phase space variables. In the strict diabatic basis {|*a*⟩} (in the sense that ⟨*a*|∇|*b*⟩ = 0 for all |*a*⟩ and |*b*⟩), the total Hamiltonian in Eq. (9) is expressed as

where $Vab(R\u0302)=\u27e8a|V\u0302(r\u0302,R\u0302)|b\u27e9$ are the matrix elements of the electronic Hamiltonian. Note that here |*a*⟩ is used to represent the strict diabatic basis, and not to be confused with the adiabatic-Fock state |*ψ*_{i}(**R**)⟩ = |*ϕ*_{α}(**R**), *n*⟩ introduced in Eq. (8). Nevertheless, based on the QD scheme, these adiabatic-Fock states with a reference geometry **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:

where 2*γ*_{b} is viewed as a parameter^{69} that specifies the ZPE of the mapping oscillators.^{36,37,69,70} In principle, 2*γ*_{b} 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)

with the nuclear force expressed as

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.^{41,47,71}

To sample the initial electronic condition and estimate the population, it is also convenient to use the action-angle variables, {*ɛ*_{b}, *θ*_{b}}, which are related to the canonical mapping variables {*p*_{b}, *q*_{b}} through

and the inverse relations

where *ɛ*_{b} 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 |*b*⟩, which is to be evaluated as^{69}

where $\rho \u0302(0)=\rho \u0302R\u2297|a\u3009\u3008a|$ is the initial density operator, *ρ*_{W}(**P**, **R**) is the Wigner transform of $\rho \u0302R$ operator for the nuclear DOFs, $\epsilon ={\epsilon 1,\epsilon 2,\u2026,\epsilon N}$ is the positive-definite action variable vector for $N$ electronic states,^{35}^{,}*W*_{a}(** ɛ**) =

*δ*(

*ɛ*

_{a}− (1 +

*γ*

_{a}))

*∏*

_{a≠b}

*δ*(

*ɛ*

_{b}−

*γ*

_{b}) is the Wigner transformed action variables,

^{72}and

*d*

**≡**

*τ**d*

**P**·

*d*

**R**·

*d*

**·**

*ɛ**d*

**. For practical reasons, the above delta functions in**

*θ**W*

_{a}(

**) are broadened using a distribution function (so-called window function) that can be used to bin the resulting electronic action variables in action-space.**

*ɛ*^{69}Furthermore, we use the

*γ*-SQC approach,

^{35}which uses a

*state-specific*and

*trajectory-specific*

*γ*

_{b}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.

^{41,71}The details of

*γ*-SQC are provided in Appendix A.

For the spin-LSC approach,^{36,37} one chooses a universal ZPE parameter 2*γ*_{b} = Γ for all states and trajectories. The spin-LSC population dynamics is calculated as

where the population estimators are obtained from the Stratonovich–Weyl transformed electronic projection operators, with the expressions as follows:^{37}

The parameter Γ is related to the radius of the generalized Bloch sphere *r*_{s} through $\Gamma =2N(rs\u22121)$, where s and $s\u0304$ are complementary indices in the Stratonovich–Weyl transform. Among the vast parameter space, one of the best performing choices^{36,37} is when $rs=rs\u0304=N+1$, which is referred to as s = W, leading to a ZPE parameter

as well as the identical expression of [|*a*⟩⟨*a*|]_{s} and $[|b\u3009\u3008b|]s\u0304$ in Eq. (41). We further use the focused initial condition^{36,37} that replaces the sampling of the mapping variables in the *d*** τ** integral of Eq. (40) with specific values of the mapping variables, such that $12(qa2+pa2\u2212\Gamma )=1$ for initially occupied state |

*a*⟩ and $12(qb2+pb2\u2212\Gamma )=0$ for the initially unoccupied states |

*b*⟩. The angle variables {

*θ*

_{b}} [Eq. (37)] are randomly sampled

^{37}in the range of [0, 2

*π*).

Using the QD propagation scheme, one can directly perform non-adiabatic using both *γ*-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 B, one can obtain the polariton coupling $\u27e8\psi i(R0)|V\u0302(R)|\psi j(R0)\u27e9$ [see Eqs. (16) and (17)] and nuclear gradient $\u2207\u27e8\psi i(R0)|V\u0302(R)|\psi j(R0)\u27e9$ [see Eqs. (24)–(30)], which are the necessary ingredients to solve the MMST mapping EOMs in Eqs. (35) and (36). Between two propagation steps, the QD basis is transformed from {|*ψ*_{i}(**R**(*t*_{0}))⟩ ≡ |*ϕ*_{α}(**R**(*t*_{0})), *n*⟩} to {|*ψ*_{j}(**R**(*t*_{1}))⟩ ≡ |*ϕ*_{β}(**R**(*t*_{1})), *m*⟩}. This leads to the corresponding transform of mapping variables between the two consecutive QD bases as follows:^{43,47}

## III. COMPUTATIONAL DETAILS

### A. 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. 1(a). 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.^{49,73,74} The electron–nuclear interaction potential operator $H\u0302en$ [cf. Eq. (2)] is expressed as

where $V\u0302p(R)$ represents the potential of the transferring proton, $V\u0302e(r)$ represents the potential of the transferring electron, and $V\u0302ep(r,R)$ 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^{75} to represent the electronic adiabatic states. These adiabatic states |*ϕ*_{α}(*R*)⟩ 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) 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)]. Figure 1(b) 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. Figure 1(c) 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. 1(d).

When coupling the SM molecular model with the cavity, the photon frequency of the cavity mode is chosen as *ℏω*_{c} = 2.721 eV ($\u22480.1$ a.u). Furthermore, we assume that the cavity field polarization direction ** ϵ** is always aligned with the direction of the dipole operator $\mu \u0302$ such that

**⋅**

*ϵ*

*μ*_{γν}=

*μ*

_{γν}(for {

*ν*,

*γ*} = {

*e*,

*g*}) where

*μ*

_{γν}is the magnitude of $\mu \u0302$. Explicitly considering the angle between

**and $\mu \u0302$ will generate a polariton induced conical intersection (even for a diatomic molecule), which will induce geometric phase effects.**

*ϵ*^{76}We consider two different light–matter coupling strengths

*g*

_{c}= 0.001 a.u. and

*g*

_{c}= 0.005 a.u. in this work. The normalized coupling strength is often defined as

^{77}

*η*≡

*g*

_{c}|

*ϵμ*_{eg}|/

*ω*

_{c}, where |

**⋅**

*ϵ*

*μ*_{eg}| is the typical magnitude of the transition dipole projected along the field polarization direction. For the coupling strength considered above (and taking

*ℏω*

_{c}≈ 2.721 eV for the model calculation), the normalized coupling strength is

*η*= 0.06 (for

*g*

_{c}= 0.001 a.u.) and

*η*= 0.3 (for

*g*

_{c}= 0.005 a.u). When 0.1 <

*η*< 1, the light and matter interaction achieves the ultra-strong coupling regime,

^{77,78}which is difficult to achieve but still within the reach of the current experimental setup.

^{79,80}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 $EJ(R)$ associated with polariton states $|EJ(R)\u3009$ [see their definition in Eq. (11)] are presented in Fig. 1(e) with the light–matter coupling strength *g*_{c} = 0.001 a.u. and in Fig. 1(f) with the light–matter coupling strength *g*_{c} = 0.005. These polariton potentials are color coded [as shown in the inset of panel (e)] based on the expectation value of $\u27e8a\u0302\u2020a\u0302\u27e9$ 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^{56,81} because the rigorous photon number operator should be obtained by applying the Power–Zienau–Woolley (PZW) gauge transformation^{54,82,83} on the photon number operator $a\u0302\u2020a\u0302$. 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.^{84} Here, we use it to characterize the photonic character for a given polariton state.

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

which corresponds to a Franck–Condon excitation of the hybrid system to the |*e*, 0⟩ state, with |*χ*⟩ as the initial nuclear wavefunction. For the SM model in this work, we use $\chi (R)=\u27e8R|\chi \u27e9\u223cexp[\u2212M\omega 0(R\u2212R0)2/2\u210f]$, where *M* is the mass of the proton (nucleus in the SM model), and *R*_{0} is the position with a minimum potential energy of the ground electronic state. Here, *χ*(*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 *ω*_{0}. We use the parameters in the original Ref. 49 for *R*_{0} = −4 and *ω*_{0} = 0.000 382 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.

### B. Details of *γ*-SQC and spin-LSC dynamics

To perform the *γ*-SQC dynamics, we need to sample the initial condition for the quantum subsystem. In this work, we first sample the action-angle variables {*ɛ*_{b}, *θ*_{b}} then transform them to the mapping variables {*p*_{b}, *q*_{b}} using Eq. (38). Among them, the action variables {*ɛ*_{b}} are sampled according to the window function in Eq. (A1), and the angle variables {*θ*_{b}} are randomly sampled from [0, 2*π*). 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 Sec. II C, where the action variable *ɛ*_{a} is set to be 1 + Γ/2 for the initially occupied state and Γ/2 for the initially unoccupied state, with Γ expressed in Eq. (42). The angle variables {*θ*_{b}} are randomly generated between [0, 2*π*) as in the *γ*-SQC method. The canonical mapping variables are obtained from Eq. (38).

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

which is the Wigner transformation of the nuclear wavefunction *χ*(*R*) = ⟨*R*|*χ*⟩ 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 |*e*0⟩. The nuclear time step used in the QD-*γ*-SQC and QD-spin-LSC is *dt* = 0.1 a.u., with 100 equally spaced electronic time steps for the mapping variables’ integration during each nuclear time step. The equation of motion in Eqs. (35) and (36) are integrated using a second-order symplectic integrator for the MMST variables^{85,86} 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*_{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 *γ* = 0 in the mapping theory [see Eq. (34)] and an initial action-angle variables condition [see Eqs. (37) and (38)] of *ɛ*_{b} = *δ*_{ab} (for the initially occupied state |*a*⟩) and *θ*_{b} = 0 (for all state |*b*⟩). One can thus use the same QD scheme and the mapping equation to obtain the results of the Ehrenfest dynamics.^{87}

## IV. RESULTS

Figure 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 (*γ*-SQC and spin-LSC), compared to the numerically exact results (solid lines). The light–matter coupling strength is chosen to be *g*_{c} = 0.001 a.u. The system is initially prepared in the |*e*0⟩ state and decays quickly into the |*g*1⟩ state during the first $\u223c12$ fs due to the large light–matter coupling strength from the large transition dipole moment (*μ*_{eg}) between adiabatic electronic state |*g*⟩ and |*e*⟩ [as shown in Fig. 1(d)]. Then, the system starts to oscillate between |*e*0⟩ and |*g*1⟩ 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 |*g*0⟩ state, which is due to the electronic NAC *d*_{eg} that directly couples the |*e*0⟩ state to |*g*0⟩ state (gold lines). All of the approximate methods can qualitatively describe such a trend, but the MQC methods [(a) and (b)] are less accurate compared to the mapping-based methods [(c) and (d)], in terms of the rising of the |*g*0⟩ population as well as its long time plateau. Moreover, both Ehrenfest and FSSH dynamics predict a significant population transfer from |*g*1⟩ to |*e*1⟩ state [(a) and (b)] as an *artifact* that is not shown in the exact dynamics results. In contrast, the *γ*-SQC and spin-LSC methods perform much better, where the population transfer process from |*g*1⟩ to |*e*1⟩ state is largely suppressed [(c) and (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 estimator specifically for computing diabatic populations.^{88} 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.^{25,27,89,90}

Figure 3 presents the polariton population dynamics with the coupling strength *g*_{c} = 0.005 a.u. The oscillation between |*e*0⟩ and |*g*1⟩ state population appears much earlier and faster compared to the *g*_{c} = 0.001 results due to the larger light–matter coupling between |*e*0⟩ and |*g*1⟩ states [see Eq. (5)]. Furthermore, the |*g*0⟩ and |*e*1⟩ states are also getting populated at an earlier time, due to the permanent dipole *μ*_{gg} and *μ*_{ee} that couples |*g*1⟩ state to |*g*0⟩ state and |*e*0⟩ state to |*e*1⟩ state, respectively. Similar to the *g*_{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 *γ*-SQC method in the description of |*g*0⟩ and |*e*0⟩ state population after *t* = 20 fs, as shown in Figs. 3(c) and 3(d).

Until now, all of our simulations are restricted in the Hilbert subspace formed by two electronic states (|*g*⟩ and |*e*⟩) and two photonic Fock states (|0⟩ and |1⟩). 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 Figs. S1 and S2 in the supplementary material. The results show that, for the small light–matter coupling strength case (*g*_{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*_{c} = 0.005 a.u), the polariton dynamics will converge when including four adiabatic electronic states (|*g*⟩, |*e*⟩, |*f*⟩, and |*h*⟩ with energies in ascending order) and four Fock states (|0⟩, |1⟩, |2⟩, and |3⟩ with photon number in ascending order).

To further test the performance of the mapping methods (*γ*-SQC and spin-LSC) as well as the MQC methods (Ehrenfest and FSSH) in such a large Hilbert subspace, which includes 16 states formed by the tensor product of four electronic states (|*g*⟩, |*e*⟩, |*f*⟩, and |*h*⟩) and four Fock states (|0⟩, |1⟩, |2⟩, and |3⟩). Figure 4 presents the results of using Ehrenfest dynamics [(a) and (b)], FSSH [(c) and (d)], *γ*-SQC [(e) and (f)] and spin-LSC [(g) and (h)]. Besides the adiabatic-Fock states already appear in the four-state subspaces (|*g*0⟩, |*e*0⟩, |*g*1⟩, and |*e*1⟩), we can also see some other states (|*h*0⟩, |*f*0⟩, |*g*2⟩, and |*f*1⟩) are populated due to the increasing light–matter coupling strength. All of the MQC (Ehrenfest, FSSH) and mapping (*γ*-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 |*g*0⟩ population). Note that both *γ*-SQC and spin-LSC perform less accurately compared to the situation in a smaller Hilbert subspace (Fig. 3). This is because both *γ*-SQC and spin-LSC are sensitive to the number of states of the system. For *γ*-SQC, more states means less trajectory landed in the population action window.^{91} For spin-LSC, the ZPE correction Γ [Eq. (42)] explicitly depends on the number of states $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,45,47} 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.^{45,47} 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, *γ*-SQC^{35} and spin-LSC,^{37} to investigate polariton dynamics of a Shin–Metiu model coupled to an optical cavity.^{23,49} 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 system at both short- and long-time dynamics compared to the exact results. In addition, the mapping methods outperform the Ehrenfest and FSSH approaches for long-time dynamics. The numerical results also demonstrate that the performance of the mapping methods (*γ*-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.

## SUPPLEMENTARY MATERIAL

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

## 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 Professor Yihan Shao.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Deping Hu**: Conceptualization (equal); Data curation (lead); Formal analysis (equal); Investigation (lead); Methodology (lead); Project administration (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). **Arkajit Mandal**: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – review & editing (supporting). **Braden M. Weight**: Conceptualization (equal); Investigation (supporting); Methodology (equal); Software (equal); Writing – review & editing (equal). **Pengfei Huo**: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (lead); Supervision (lead); Validation (equal); Writing – original draft (equal); Writing – review & editing (lead).

## DATA AVAILABILITY

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

### APPENDIX A: DETAILS OF THE *γ*-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^{35,72} is expressed as

where the window functions are defined as

and

and trajectories are assigned to state *b* at time *t* if *ɛ*_{b} ≥ 1 and *ɛ*_{b′} < 1 for all *b*′ ≠ *b*. The ZPE parameter in the standard SQC method with triangle windowing is *γ* = 1/3.

The time-dependent population of the state |*b*⟩ 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:

In the *γ*-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^{92} to carry out on-the-fly simulations of the methaniminium cation.^{71} The basic logic of *γ*-SQC is to choose an *γ*_{b} for each state |*b*⟩ 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 |*a*⟩,

or equivalently,

where the {*ɛ*_{b}} are uniformly sampled inside the window function [Eq. (A1)], and then the *γ*_{b} are chosen to satisfy Eq. (A6).

These *γ*_{b} will be explicitly used in the EOMs in Eqs. (35) and (36), and, in particular, the nuclear forces are now

ensuring the initial forces (at *t* = 0) are simply **F** = −∇*V*_{aa}(**R**). Previously, without any adjustments to *γ*_{b}, the chosen values for *γ*_{b} were only dependent on the windowing function itself, i.e., *γ*_{b} = 0.366 for the square Windows and *γ*_{b} = 1/3 for the triangle windows. With the above *γ*-correction method,^{35} each individual trajectory will have its own state-specific *γ*_{b} for state |*b*⟩ 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^{75} to represent the electronic adiabatic wavefunction and solve Eq. (3). The grid of DVR is uniform with spacing Δ*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 $H\u0302el$ in this grid basis {|*r*_{i}⟩} are given by

where the $\u27e8ri|T\u0302r|rj\u27e9$ is given analytically^{75} as follows:

Directly diagonalizing the matrix of $H\u0302el$ [in Eq. (B1)] at a given nuclear (proton) position *R* in this grid basis gives the accurate adiabatic electronic states

where $ci\alpha (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 Eqs. (23a), (23b), and (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:

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

Using |*ϕ*_{α}(*R*)⟩ and |*ϕ*_{β}(*R*)⟩ in the grid basis, we can directly evaluate the nuclear gradient $\u27e8\varphi \lambda |\u2207H\u0302en|\varphi \nu \u27e9$ [Eq. (25)] as follows:

where the $\u2207H\u0302en(R)$ is evaluated analytically using the expression of $H\u0302en(R)$ in Eq. (44). This gives the adiabatic gradient ∇*E*_{λ} (for *λ* = *ν*) and derivative coupling (for *λ* ≠ *ν*) as $d\lambda \nu =\u27e8\varphi \lambda |\u2207H\u0302en|\varphi \nu \u27e9/(E\nu \u2212E\lambda )$, as indicated in Eq. (44). The nuclear gradient $\u27e8\varphi \lambda |\u2207H\u0302en|\varphi \nu \u27e9$ is also one of the key ingredient in the QD propagation.

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

To solve the exact quantum dynamics, we represent the total wavefunction of the hybrid system as $|\Psi \xi \u3009=\u2211i,kcik\xi |\psi i(Rk)\u3009\u2297|Rk\u3009$, where {|*R*_{k}⟩} is the DVR grid basis for the nucleus, and the |*ψ*_{i}(*R*_{k})⟩ = |*ϕ*_{α}(*R*_{k})⟩ ⊗|*n*⟩ is the adiabatic-Fock basis [Eq. (8)], where the electronic adiabatic basis |*ϕ*_{α}(*R*_{k})⟩ 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 $cik\xi $ and the eigenvalue of the total Hamiltonian $H\u0302$ [Eq. (1)] will be obtained by solving the time-independent Schödinger equation $H\u0302|\Psi \xi \u3009=E\xi |\Psi \xi \u3009$. We use the Sinc DVR basis^{75} for the nuclear DOF and solve the above eigenvalue problem to obtain all the eigenvalues and eigenstates. We use finer grid points for nucleus Δ*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)\u3009=\u2211\xi C\xi \u2061exp\u2212i\u210fE\xi t|\Psi \xi \u3009$, where *C*_{ξ} is the projection of initial total wavefunction onto the |Ψ_{ξ}⟩ as *C*_{ξ} = ⟨Ψ_{ξ}|Φ(0)⟩, with the initial wavefunction |Φ(0)⟩ 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 (*γ*-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

where |*ψ*_{i}(**R**(*t*))⟩ = |*ϕ*_{α}(**R**(*t*))⟩ ⊗|*n*⟩ is the adiabatic-Fock state basis [Eq. (8)]. The quantum subsystem is described by the time-dependent Schödinger equation (TDSE)

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)

where **c**^{†} is the transpose of the coefficient column vector **c** expressed as follows:

The nuclear gradient matrix is expressed as

where [*V*] and [**d**] are the matrix of $V\u0302$ 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 $|EI(R(t))\u3009$ [see definition in Eq. (11)]

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 $|EI(R(t))\u3009$ [eigenstate of $V\u0302$, see Eq. (11)] as follows:

where $EI$ 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

where the matrix element of ∇*V* is expressed in Eq. (C5). 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 $|EI\u3009$ to *any other* polariton state $|EJ\u3009$ during the time interval between *t* and *t* + *δt* is

where $\rho IJpl(t)$ is the reduced density matrix element in the polariton basis expressed as follows:

Since the probability should be positive definite, one sets^{14} *f*_{IJ} to 0 if *f*_{IJ} < 0. The non-adiabatic transition, i.e., stochastic switch from the currently occupied state $|EI\u3009$ to another state $|EK\u3009$, occurs if the following condition is satisfied:

where *ζ* 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 $|EK\u3009$, while the velocities of the nuclei are rescaled along the direction of the NAC vector **d**_{IK}(**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 |*ψ*_{i}⟩ = |*e*0⟩ 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⟩ is not one of the eigenstates $|EI(R)\u3009$. We thus follow the previous work^{88,93} and use the Monte Carlo scheme to randomly choose the initial active state $|EI(R)\u3009$ for each trajectory, based on the *magnitude* of $|\u27e8EI(R)|e(R),0\u27e9|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 $V\u0302$, there is no unique way to calculate them in the FSSH approach.^{88} In the main text, we present the populations of the adiabatic-Fock states using the expansion coefficients *c*_{I}(*t*) in Eq. (C6). There are, of course, alternative ways to compute populations of these adiabatic-Fock states.^{88} Below, we explore the alternative ways to compute them.

For clarity, we denote the reduced density matrix in the adiabatic-Fock basis as $\rho ijaf(t)$, and $\rho IJpl(t)$ is the reduced density matrix in the polariton basis [expressed in Eq. (C10)]. To get the adiabatic-Fock state population of the |*ψ*_{i}⟩ state $\rho iiaf$ from the FSSH simulation, the most straightforward way (as the results presented in the main text) is through following unitary transformation:

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

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

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^{88} to calculate the populations that are not in an adiabatic representation. These methods vary on how to calculate the polaritonic state density matrix [*ρ*^{pl}(**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

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. (C12), and the ensemble average over all trajectories is computed as described in Eq. (C14).

The other choice^{88} (motivated by the mixed quantum-classical Liouville approach) is to calculate the diagonal elements of *ρ*^{pl} using the active state index, and calculate the off-diagonal elements using the polaritonic state expansion coefficients {*c*_{I}(*t*)}

where *K* is the active index. The adiabatic-Fock state population is then obtained from the same transformation described in Eq. (C12), and the ensemble average over all trajectories as described in Eq. (C14).

Following the same notation as used in Ref. 88, we refer to the choice in Eq. (C15) as method 1, the choice in Eq. (C13) as method 2 (same as the FSSH results presented in the main text), and the choice in Eq. (C16) 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 methods 1 and 2, consisting with the conclusion in Ref. 88.