We derive a rigorous nuclear gradient for a molecule-cavity hybrid system using the quantum electrodynamics Hamiltonian. We treat the electronic–photonic degrees of freedom (DOFs) as the quantum subsystem and the nuclei as the classical subsystem. Using the adiabatic basis for the electronic DOF and the Fock basis for the photonic DOF and requiring the total energy conservation of this mixed quantum–classical (MQC) system, we derived the rigorous nuclear gradient for the molecule–cavity hybrid system, which is naturally connected to the approximate gradient under the Jaynes–Cummings approximation. The nuclear gradient expression can be readily used in any MQC simulations and will allow one to perform the non-adiabatic on-the-fly simulation of polariton quantum dynamics. The theoretical developments in this work could significantly benefit the polariton quantum dynamics community with a rigorous nuclear gradient of the molecule–cavity hybrid system and have a broad impact on the future non-adiabatic simulations of polariton quantum dynamics.

## 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} These polariton states have hybridized potential energy surfaces (PESs) (and their curvatures) from both the ground and the excited electronic states,^{5,6} which have been shown to facilitate new chemical reactivities.^{1,6–9} Thus, polariton chemistry provides a new paradigm for chemical transformations. Theoretical investigations play a crucial role in understanding the 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 that does not include quantized photons or quantum optics that does not have a well-defined theory to include the influence of nuclear vibrations.^{2}

Historically, the mixed quantum–classical (MQC) approaches^{14–16} have played an important role in simulating the non-adiabatic dynamics of the coupled electronic–nuclear DOFs. Two of the most commonly used MQC methods are Ehrenfest dynamics and fewest switches surface hopping (FSSH) approaches.^{17,18} Both methods describe the electronic subsystem quantum mechanically and treat the nuclear DOF classically. It is thus a natural idea for the theoretical chemistry community to extend these MQC approaches to investigate polariton chemistry by treating the electronic–photonic DOF (or the so-called polariton subsystem) quantum mechanically and describing the nuclear DOF classically. Indeed, incorporating the description of the photon field into the MQC methods has become a basic problem of crucial importance to study polariton chemistry.^{10–13,19–22}

There has been enormous progress during the past few years in developing MQC approaches to simulate polariton dynamics. Kowalewski *et al.*^{23} derived the expression of the derivative couplings using the Jaynes–Cummings (JC) model Hamiltonian^{24} where the rotating wave approximation (RWA) is assumed for the molecule–cavity coupling term. Groenhof and co-workers developed a multi-scale simulation approach combining Tavis–Cummings (TC) model^{25} using FSSH approach^{10,11} or Ehrenfest dynamics^{12,13} to simulate an ensemble of molecules coupled to a cavity. Zhang *et al.* extended the JC and TC models to include multiple molecular excited states,^{22} derived the corresponding nuclear gradient, and performed MQC simulations for polariton dynamics. Fregoni *et al.* developed the MQC simulations with the JC-type model (that excludes the permanent dipole moment) to perform FSSH simulations of azobenzene photoisomerization in cavity.^{19–21} Foley and co-workers have incorporated cavity loss through non-Hermitian approach^{26} and performed FSSH simulations of a model photoisomerization reaction coupled to an optical cavity.^{27} Furthermore, Li and Hammes-Schiffer^{28} have recently developed a semi-classical approach for molecular polaritons by self-consistently propagating the real-time dynamics of classical cavity modes and a quantum molecular subsystem described by the nuclear–electronic orbital (NEO) method.

The key ingredient in the MQC simulations of polariton dynamics is the expression of the nuclear gradient, which is a necessary component for propagating the motion of nuclei. However, previous MQC simulations often use the rotating wave approximation, exclude the dipole self-energy (DSE) terms,^{29} or neglect the permanent dipole moment,^{20} all of which could change the fundamental polariton quantum dynamics in the strong and ultra-strong coupling regimes.^{30–34} Nuclear gradient expressions beyond the JC models have been derived under the full Configuration-Interaction (CI) expansion framework.^{20} Rigorous nuclear gradient expressions under the quantum electrodynamics Linear-Response time-dependent Density Functional Theory (QED-LR-TDDFT)^{35} framework using Pauli–Fierz (PF) type QED Hamiltonian.^{3,33} However, these gradient expressions lack a clear physical picture of light–matter interactions, as well as a clear connection to the more intuitive (but less accurate) gradient of the JC-type Hamiltonian.

In this paper, we derive a rigorous expression of the nuclear gradient using the PF QED Hamiltonian without making unnecessary approximations. We treat the molecule–cavity hybrid system as a MQC system, where the electronic–photonic DOFs are treated quantum mechanically, and the nuclear DOFs are treated classically. In this work, the electronic DOFs are described using the electronic adiabatic states, and the photonic DOFs are described with the Fock states. Requiring the total energy conservation of this MQC system, we derive the rigorous nuclear gradient expression [Eq. (29)], which is intuitively connected to those approximate gradient expressions under the JC approximation.^{23} Our gradient expressions, and the corresponding MQC simulation approaches, can in principle include any number of electronic states and photon states, which is a desired property for investigating the dynamical process of polariton chemistry. The nuclear gradient expression derived here can be readily used in any MQC simulations, such as Ehrenfest and FSSH approaches, as demonstrated in the Result section.

## II. THEORY AND METHODS

### A. The Pauli–Fierz QED Hamiltonian

The Pauli–Fierz (PF) QED Hamiltonian in the dipole gauge describing one molecule coupled to quantized radiation field inside a 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 DSE, 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. 36. A detailed discussion of Hamiltonian under different gauges can be found in Ref. 33.

The electronic–nuclear potential, $H\u0302en$, is expressed as follows:

which describes the molecular Hamiltonian (without the nuclear kinetic energy). The above expression in Eq. (2) includes electronic kinetic energy $T\u0302e$, electron–electron interaction $V\u0302ee$, electronic–nuclear interaction $V\u0302en$, and nuclei–nuclei interaction $V\u0302nn$. The expressions of these four terms can be found in standard textbooks.^{37–39} Modern electronic structure theory is focused on solving the eigenvalue problem of $H\u0302en$, providing the following *adiabatic* energy and corresponding adiabatic states as follows:

Here, |*ϕ*_{ν}(**R**)⟩ represents the *ν*_{th} many-electron adiabatic state for a given molecular system, with the adiabatic energy *E*_{ν}(**R**).

The photonic Hamiltonian is written as

where $a\u0302\u2020$ and $a\u0302$ are the photonic creation and annihilation operator, respectively, $q\u0302c=\u210f/2\omega c(a\u0302\u2020+a\u0302)$ and $p\u0302c=\u2212i\u210f\omega c/2(a\u0302\u2212a\u0302\u2020)$, and *ω*_{c} is the photon frequency inside the cavity. For clarity, we restrict our discussions to the cavity with only one photonic mode, and all of the theoretical expressions presented here can be easily generalized into a cavity with many modes.^{13,40,41}

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

where *λ*_{c} = *λ*_{c} ⋅ ** ϵ** characterizes the photon–matter coupling strength, and

**is the direction of the field polarization. Furthermore, the total dipole operator of both electrons and nuclei is defined as**

*ϵ*where −*e* is the charge of the electron and *Z*_{α}*e* is the charge of the *α*_{th} nucleus. The coupling light–matter strength is determined by the volume of the cavity as $\lambda c=1/\u03f50V0$, where *ϵ*_{0} is the permittivity inside the cavity and *V*_{0} is the effective quantization volume inside the cavity. Another way to define the light–matter coupling strength is using $gc=\u210f\omega c/2\lambda c$. Note that the *g*_{c} used in this work differs from the common notation used in the literature based on the quantum optics models,^{6,42} which often include the magnitude of the dipole inside the definition of *g*_{c}. This is because, in the quantum optics models, the magnitude of the dipole is treated as a constant, whereas in molecular cavity QED, the dipole matrix elements explicitly depend on **R**.

Finally, the dipole self-energy is expressed as

This is a necessary term of the PF Hamiltonian, in order to make sure gauge invariance of the PF Hamiltonian^{9,33} and a bounded ground polariton state.^{9,43,44}

Recent investigations of polariton photochemistry have been mainly focused on using the JC model^{22,29} or TC model^{12} to describe the quantum light–matter interactions. These models usually only consider two electronic states {|*g*⟩, |*e*⟩} and the transition dipole $\mu ge(R)=\u27e8g|\mu \u0302|e\u27e9$ among them, where the permanent dipole is often ignored. In this context, one can further define the creation and annihilation operators for molecular excitation as $\sigma \u0302\u2020\u2261|e\u3009\u3008g|$ and $\sigma \u0302\u2261|g\u3009\u3008e|$, and thus $\mu \u0302=\mu eg(R)\u22c5(\sigma \u0302\u2020+\sigma \u0302)$. The molecule–cavity interaction term in Eq. (5) can now be expressed as

Assuming the RWA by ignoring the counter-rotating terms (CRTs) $a\u0302\u2020\sigma \u0302\u2020$ and $a\u0302\sigma \u0302$, and explicitly dropping the DSE term $H\u0302d$ in Eq. (7), one arrives at the following JC model:

The Rabi model, on the other hand, ignores the DSE and the permanent dipole.^{45–47}

In this work, we do not consider the cavity loss due to the interactions of the cavity mode and non-cavity modes. Modeling such an effect can be accomplished by incorporating Lindblad type decay dynamics with MQC simulations.^{48}

### B. Nuclear gradient in the mixed-quantum classical simulations of polariton dynamics

For the molecule–cavity hybrid system, a convenient basis for quantum dynamics simulations could be the dressed states

where the quantum number *i* ≡{*ν*, *n*} indicates both the adiabatic electronic state of the molecule and the Fock states. Note that we have introduced a shorthand notation in Eq. (10), which will be used throughout the rest of this paper. This is one of the most straightforward choices of the basis for the hybrid system because of the readily available adiabatic electronic states from the electronic structure calculations, as well as the dipole moments we need to construct the elements of Hamiltonian.

Here, we provide a rigorous derivation of general nuclear gradient expression used in MQC simulations in a real orthogonal basis, following the similar procedure of Tully.^{49} In the MQC simulation, such as the Ehrenfest dynamics or the FSSH approach, 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$ is expressed in Eq. (2). For a molecule–cavity hybrid system,

which is commonly referred to as the polariton Hamiltonian,^{3,50} also denoted as $H\u0302pl$. For the molecule–cavity hybrid system, we define the polaritonic state^{3,50} as the eigenstate of $V\u0302=H\u0302pl$ [see definition in Eq. (12)] through the following eigenequation:

where $|EI(R)\u3009$ is the polariton state with polariton energy $EI(R)$.

In MQC dynamics simulations, one treats the nuclear DOF classically, such that the Hamiltonian in Eq. (11) becomes

The electronic–photonic wave function is expanded in the basis |*ψ*_{i}(**r**; **R**(*t*))⟩,

where *N*_{b} is the number of basis we use and *c*_{i} is the expansion coefficients. The wave function satisfies the time-dependent Schödinger equation (TDSE),

where $Vij=\u27e8\psi i(R)|V\u0302|\psi j(R)\u27e9$. The time derivative coupling (or non-adiabatic coupling, NAC) between two basis states is

where $dij\alpha $ is the NAC associated with the *α*_{th} atom, defined as follows:

Using the above notations, Eq. (17) can be rewritten as

The total energy for the MQC system expressed in Eq. (14) can be expressed as

where *M*_{α} is the nuclear mass of *α*_{th} atom and **P**_{α} is the corresponding momentum. In order to get the equation of motion for the classical nuclei, we resort to the conservation of the above total energy.^{49} Setting the time derivative of the above total energy in Eq. (21) to zero, i.e. *dE*/*dt* = 0, we obtain

where we have used the chain rule with respect to Hamiltonian matrix elements. As shown in Appendix A, using Eq. (20), we can prove the following identity:

These results are reduced to the familiar expression for an isolated molecule, as shown in Appendix B.

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

Equation (24) can then be written as a more compact form

where we have used [ ⋅⋅⋅ ] to denote a matrix, and the gradient of potential matrix is defined as

where [*V*] and [**d**^{α}] are the matrix of $V\u0302$ and derivative coupling operator, respectively, and we have defined the matrix

We can write the matrix elements of the nuclear gradient as follows:

This is the most general expression of the nuclear gradient with a *real orthogonal basis* that explicitly depends upon **R**. Note that although we introduced the short-hand notation [∇_{α}*V*] in Eq. (27), it can be justified as the matrix of the operator $\u2207\alpha V\u0302$ using the chain rule as follows:

which is indeed equivalent to Eq. (29).

### C. Nuclear gradient of molecule-cavity hybrid systems

In this section, we provide the detailed expression of the nuclear gradient for molecule–cavity hybrid systems. With the adiabatic-Fock basis |*ϕ*_{ν}*n*⟩ and |*ϕ*_{γ}*m*⟩ introduced in Eq. (10), the matrix elements of every term in $V\u0302$ [Eq. (12)] can be explicitly expressed as follows (using the properties of creation and annihilation operators of photonic DOF):

where *E*_{ν}(**R**) and *μ*_{γν}(**R**) explicitly depend on nuclear position **R**. In Eq. (31d), $D\gamma \nu 2$ denotes the matrix elements of DSE, and the sum *∑*_{ξ} in the matrix element of $H\u0302d$ runs over the adiabatic states |*ϕ*_{ξ}(**R**)⟩ considered in the calculation (as opposed to including all possible adiabatic states of the molecule). Note that this effectively projects $\mu \u0302$ inside the DSE term within the matter subspace, with the projection operator $P\u0302=\u2211\xi |\varphi \xi \u3009\u3008\varphi \xi |$. This matter state truncation scheme^{33,51} makes sure that all operators are properly confined in the same truncated electronic subspace^{33} in order to generate consistent and meaningful results. Indeed, if the electronic Hilbert space identity is $I\u0302=P\u0302+Q\u0302$, where $Q\u0302$ contains the adiabatic electronic states outside the subspace defined by $P\u0302$, $P\u0302\mu \u0302P\u0302$ is properly confined in the subspace $P\u0302$, whereas $P\u0302\mu \u03022P\u0302=P\u0302\mu \u0302(P\u0302+Q\u0302)\mu \u0302P\u0302$ contains the terms outside the subspace $P\u0302$. More details of this discussion can be found in Refs. 52 and 33, and a numerical example that demonstrates the validity of the truncation scheme can be found in Fig. S2 of Ref. 33. Note that if solving the entire polariton eigenvalue problem [Eq. (13)] using QED-electronic structure theory, the DSE can be expressed in a compact form that explicitly contains quadruple contributions.^{26} Furthermore, in Eq. (31), we have used the matrix element of the dipole operator $\mu \u0302$ [Eq. (6)] under the adiabatic representation

which parametrically depends on **R**.

To express the nuclear gradients in Eq. (29), we need the gradients ∇_{α}*V*_{ij} and derivative coupling matrix [**d**^{α}] that often can be directly obtained from the *ab initio* electronic structure calculations or numerical differential techniques.^{22} To calculate ∇_{α}*V*_{ij}, we need to take the derivative on each term of *V*_{ij}(**R**) expressed in Eq. (31), including ∇_{α}*E*_{ν}(**R**) and ∇_{α}(** ϵ** ⋅

*μ*_{γν}(

**R**)). To construct the derivative coupling matrix [

**d**

^{α}], one can simply use

because the Fock states do not explicitly depend upon **R** and are orthonormal to each other. In the above equation, $d\gamma \nu \alpha $ is the regular derivative couplings among adiabatic states of the molecule. Using the above information, we can explicitly write down all matrix elements of the PF Hamiltonian and the nuclear gradients [using Eq. (29)].

To summarize, the most general results of the nuclear gradient $[\u2207\alpha V]ij$ for a molecule–cavity hybrid system is presented in Eq. (29), with the matrix elements of *V*_{ij} expressed in Eq. (31) and the matrix elements of **d**^{α} in Eq. (33). This is the *central* results of this paper.

Below, we give more detailed expressions of these gradient expressions in two specific subspaces. For clarity and concreteness, we will only consider two adiabatic electronic states |*g*⟩ and |*e*⟩, which are the electronic ground and first excited states of the molecule, respectively. The eigenvalues corresponding to these two states |*g*⟩ and |*e*⟩ are *E*_{g} and *E*_{e}, respectively. The photon Fock basis considered here is from vacuum state |0⟩ to some finite number that assures the convergence of the properties we calculate. Unless explicitly stated, we always assume that the photon number in the basis is in the ascending order. For the state with the same photon number, the electronic states are also in the ascending order. Under the above settings, the two simplest basis one can choose are {|*e*0⟩, |*g*1⟩} and {|*g*0⟩, |*e*0⟩, |*g*1⟩, |*e*1⟩}, which are referenced hereafter as two-state and four-state basis, respectively.

Under the two-state subspace {|*e*0⟩, |*g*1⟩}, the potential matrix for the molecule–cavity hybrid system under this two-state basis is

The derivative coupling between |*g*1⟩ and |*e*0⟩ states is

due to the fact that ⟨1|0⟩ = 0, even though that ⟨*g*(**R**)|∇_{α}|*e*(**R**)⟩ ≠ 0. Thus, in the {|*e*, 0⟩, |*g*, 1⟩} subspace, the derivative coupling among these two states is 0, making them effectively “diabatic” even though |*e*⟩ and |*g*⟩ themselves are adiabatic states. Furthermore, using the general results in Eqs. (31)–(33), we have [**d**^{α}] = 0 because of ⟨*g*1|∇_{α}|*g*1⟩ = 0, ⟨*g*1|∇_{α}|*e*0⟩ = 0, and ⟨*e*0|∇_{α}|*e*0⟩ = 0.

The corresponding nuclear gradient thus only comes from ∇_{α}[*V*], which is explicitly expressed as

The above formula (without the DSE related terms) can be directly derived^{23} using the JC model [Eq. (9)] and has been used to compute the nuclear force of two-level molecules coupled to a cavity.^{10}

In the four-state basis {|*g*0⟩, |*e*0⟩, |*g*1⟩, |*e*1⟩}, the potential matrix is

The above matrix is Hermitian because *μ*_{eg} = *μ*_{ge} and $Deg2=Dge2$. The derivative coupling matrix [**d**^{α}] is

where $dge\alpha =\u27e8g|\u2207\alpha |e\u27e9=\u2212deg\alpha $. The elements are non-zero only when the corresponding basis (ket and bra) have the same photon number, according to Eq. (33). The nuclear gradients can be decomposed into two matrices as indicated by Eq. (29), with the ∇_{α}[*V*] expressed as

as well as the [**X**^{α}] matrix as follows:

All of the quantities appear in [*V*] and the gradient matrix [Eqs. (37)–(40)] can be obtained when solving the electronic structure problem [Eq. (3)] by calculating the adiabatic states of the molecular Hamiltonian $H\u0302en$ [Eq. (2)].

Furthermore, one can clearly see that the DSE term is a necessary component in $H\u0302PF$ [Eq. (7)], which is originated from the Power–Zienau–Woolley (PZW) gauge transformation on the minimum coupling Hamiltonian. These DSE terms also explicitly show up in the matrix elements of $V\u0302$ [Eq. (31d)] as well as the nuclear gradient terms [Eqs. (39) and (40)]. Without DSE, the gauge invariance will explicitly break down.^{30,43,52,53} This is a well-known result in QED as well as revisited in the current literature.^{9,30,51–53} The importance of including the DSE terms in the PF Hamiltonian as well as in the nuclear gradient expression is discussed in Appendix C, with numerical examples.

### D. Gradient expression in the polaritonic basis

Until now, the nuclear gradient expressions were formulated in the adiabatic-Fock basis (i.e., photon-dressed electronic states), which is not the eigenbasis of $V\u0302$ [defined in Eq. (12)]. Some MQC methods, such as FSSH,^{17} are formulated specifically in the eigenstates of $V\u0302$ and require the calculation of the nuclear gradients on the eigenstates of $V\u0302$. This means that one needs to obtain the polariton eigenstates ${|EI(R)\u3009}$ by solving Eq. (13). The transformation between {|*ψ*_{i}(**R**)⟩} = {|*ϕ*_{ν}(**R**)⟩ ⊗|*n*⟩} and ${|EI(R)\u3009}$ basis can be accomplished by applying the transformation matrix **U** that diagonalizes the matrix [*V*] as

where the diagonal elements $EI(R)$ are the eigenvalues of the polariton states $|EI(R)\u3009$. The adiabatic gradient (diagonal terms) and NAC between polaritonic states can thus be obtained by transforming the gradient matrix [∇_{α}*V*] using **U**.

To connect with the nuclear gradient expressions in the previous literature, we derive the gradient explicitly in polaritonic basis. The nuclear gradient associated with the $|EI(R)\u3009$ polaritonic state is

where we have used the Hellman–Feynman theorem. Assuming the completeness relation $\u2211i|\psi i\u3009\u3008\psi i|=I\u0302$ (where |*ψ*_{i}⟩ = |*ϕ*_{ν}(**R**)⟩ ⊗|*n*⟩) and inserting it into Eq. (42), we have

where the transformation matrix element is $UkI=\u27e8\psi k|EI\u27e9$. The gradient term ⟨*ψ*_{j}|∇_{α}*V*|*ψ*_{k}⟩ can be evaluated using Eq. (29).

The non-adiabatic coupling (derivative coupling) between two *polaritonic states* (when *I* ≢ *J*) can be expressed as^{18}

Note that this should not be confused with the molecular derivative coupling $dij\alpha $ defined in Eq. (33). One can further express Eq. (44) by inserting the completeness relation as

where $UlJ=\u27e8\psi l|EJ\u27e9$ and ⟨*ψ*_{k}|∇_{α}*V*|*ψ*_{l}⟩ can be evaluated using Eq. (30).

To give a more concrete example, for the two-state basis {|*e*0⟩, |*g*1⟩}, we have

where $\u2207\alpha E1$ and $\u2207\alpha E2$ are gradients of two adiabatic (polaritonic) states and $d12\alpha $ is the NAC between them. Since one can obtain the analytical formula of the transformation matrix **U** to diagonalize a 2 × 2 matrix, we can explicitly write down the analytical formula for the NAC for this special case. To diagonalize the [*V*] matrix in Eq. (34), we have the following **U** matrix:

where the mixing angle *θ* satisfies the condition^{54}

with $\mu \u0303ge=gc\u03f5\u22c5\mu ge$ and $\Delta V=Ve+Dee2\u2212Vg\u2212Dgg2\u2212\u210f\omega c$, which is the difference between diagonal elements in Eq. (34).

The two polaritonic states can be expressed as

which are commonly referred to as the upper polariton state (for | + ⟩) and lower polariton state (for | − ⟩).

The NAC between these two states is $d12\alpha =\u27e8E1|\u2207\alpha |E2\u27e9$. Using Eq. (49), one can find that

Using the expression of the mixing angle *θ* in Eq. (48), we have

The above $d12\alpha $ expression is identical to the derivative couplings derived by Kowalewski and Mukamel^{23} for a JC light–matter interaction model [Eq. (9)], subject to a difference with the presence of DSE terms in the expression of $\Delta V$. This is not surprising because, within the {|*e*0⟩, |*g*1⟩} subspace, the light–matter interaction Hamiltonian reduces to a JC type Hamiltonian [see Eq. (34)].

## III. DETAILS OF MODEL CALCULATIONS

### A. Shin–Metiu model for the molecule

In this work, we use the Shin–Metiu (SM) model^{55} as the “*ab initio*” model molecular system to investigate strong and ultra-strong light–matter interactions between a molecule and an optical cavity. This simple and idealized model allows us to solve the polariton dynamics exactly, thus providing an ideal test case for assessing the validity of the nuclear gradient expressions. This model, on the other hand, also contains realistic electron–nuclear interactions and goes beyond the simple diabatic system-bath type models. We choose two different parameter sets of the SM model, and hereafter we refer to them as SM1 and SM2, respectively. The SM1 is the model-I documented in Ref. 55. The electronic–nuclear potential reads as

where the *r* and *R* are the position of the electron and the nucleus, respectively. The distance between two fixed ions is *L* = 18.897 a.u., and the cut-off for the modified Coulomb interaction *R*_{c} = 2.8345 a.u.

To calculate the electronic properties of the SM model, we use the Sinc discrete variable representation (DVR) basis^{56} to represent the electronic adiabatic wavefunction and solve Eq. (3). The grid of DVR is uniform with spacing Δ*r* = 0.147 in the range [ −22, 22] a.u. for the electronic DOF. 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^{56} as follows:

Directly diagonalizing the matrix of $H\u0302el$ [in Eq. (53)] 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.

Furthermore, the matrix elements for dipole moment operator [Eq. (32)] are calculated as

Figure 1(a) presents the two lowest adiabatic electronic states [defined in Eq. (3)] of SM1 model (red and blue curves), and Fig. 1(c) presents the NAC between them (green curve). The energy gap at *R* = 0 a.u. between these two states is about 1.281 eV. The NAC is in the order of 0.2 a.u. and has a peak at *R* = 0 a.u. as well. The matrix elements of the dipole moment under the adiabatic representation [Eq. (32)] of SM1 are presented in Fig. 1(e). Note that, in some nuclear configurations, the magnitude of the permanent dipoles (*μ*_{gg} and *μ*_{ee}) could be even larger than the transition dipole *μ*_{eg}, making the atomic JC model invalid. This can significantly change the polaritonic potential energy surfaces (PESs) when the light–matter coupling is strong enough and cause some strange dynamical behaviors as we will discuss later.

The SM2 model, which is an asymmetrical proton-coupled electron transfer model, is directly adapted from Ref. 40. The electron–nucleus interaction potential operator is

We choose the same parameters used in Ref. 40, 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. For the polariton dynamics simulation when coupling this molecule with cavity, the photon energy of the cavity mode is chosen as *ℏω*_{c} = 2.721 eV ($\u22480.1$ a.u.). Figure 1(b) presents the two lowest adiabatic electronic states of SM2 (red and blue curves). Figure 1(d) presents the NAC between them (green curve). The matrix elements of the dipole moment under the adiabatic representation [Eq. (32)] of SM2 are presented in Fig. 1(f).

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 cause geometric phase effects.**

*ϵ*^{57}We consider three different light–matter coupling strengths

*g*

_{c}= 0.001 a.u.,

*g*

_{c}= 0.005 a.u., and

*g*

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

^{30}

*η*≡

*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.7 eV for the model calculation), the normalized coupling strength is

*η*= 0.06 (for

*g*

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

*η*= 0.3 (for

*g*

_{c}= 0.005 a.u.), and

*η*= 0.6 (for

*g*

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

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

^{30,52}which is difficult to achieve but within the reach of the current experimental setup.

^{58,59}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.

For both models, the initial states (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 both two models in this work, we use $\chi (R)=\u27e8R|\chi \u27e9\u223cexp[\u2212M\omega (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 is *ω*. For SM1, *R*_{0} = −4.156 and *ω* = 0.00270 a.u.; for SM2, we use^{40} *R*_{0} = −4 and *ω* = 0.000382 a.u.

### B. Details of MQC simulations

For the MQC simulations, all of the electronic properties of SM model are calculated on-the-fly using the Sinc DVR method,^{56} including the adiabatic electronic eigenstate, dipole moments, and derivative coupling. In particular, using |*ϕ*_{α}(*R*)⟩ in the grid basis, we can directly evaluate the adiabatic nuclear gradient $\u27e8\varphi \lambda |\u2207H\u0302en|\varphi \nu \u27e9$ as follows:

where the $\u2207H\u0302en(R)$ is evaluated analytically using the expression of $V\u0302en(R)$ in Eqs. (52) and (57). This gives the adiabatic gradient ∇*E*_{γ} (when *γ* = *ν*) and derivative coupling (for *γ* ≠ *ν*) as

Furthermore, the nuclear gradient expression [Eq. (39)] for a polariton system requires the derivative on the dipole matrix [Eq. (56)]. This requires the evaluation of the derivative of the expansion coefficients $\u2207ci\gamma (R)$ in Eq. (56). 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 perform the Ehrenfest dynamics simulation, we use the TDSE in Eq. (20) for quantum subsystem (electronic–photonic DOFs), and equation of motion in Eq. (24) for the classical subsystem (nuclear DOF). We use the fourth-order Runge–Kutta method to integrate the TDSE and the velocity Verlet algorithm to integrate Newton’s equation of motion. 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. The total energy is well conserved for all the trajectories.

We use Tully’s FSSH^{17,18} algorithm to perform surface hopping simulations for polariton dynamics. Note that a similar simulation has been performed recently,^{19,21,22} and our focus here is to test the importance of the new gradient expressions. The equation of motion for electronic wavefunction in FSSH is described in Eq. (20), which is the same as Ehrenfest dynamics. The main difference is how the nuclear force is treated. In FSSH, the nuclear force is contributed from *only one* specific eigenstate of $V\u0302$, the so-called active state as follows:

where $|EI\u3009$ is the active adiabatic polariton state (which is the eigenstate of $V\u0302$). The active state index will be determined at every single nuclear propagation step. According to the “fewest switches” algorithm,^{17} the probability of switching (probability flux) from the active polariton state |*I*⟩ to any other polariton state |*J*⟩ during the time interval between *t* and *t* + *δt* is

where $\rho IJ=cI*cJ$ is the electronic density matrix element. Since the probability should not be negative, we set *f*_{IJ} to 0 if *f*_{IJ} < 0. Then, the non-adiabatic transition, i.e., stochastic switch from the current occupied state |*I*⟩ to another state |*K*⟩, occurs if the following expression is met:

where *ξ* is a 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 |*K*⟩, while the velocities of the nuclei are rescaled along the direction of the non-adiabatic coupling vector **d**_{IK}(**R**) in order to conserve the total energy

where the universal scaling constant *τ*_{IK} is calculated with the smallest magnitude obtained from the following expression:^{18}^{,}$\tau IK=1AIKBIK\xb1BIK2+2AIK\u22c5CIK$, where $AIK=dIK2(R)/M$, $bIK=R\u0307\u22c5dIK(R)$, and $CIK=EI(R)\u2212EK(R)$. When the nuclear kinetic energy is not large enough to compensate the polariton potential difference, a frustrated hop occurs,^{60} meaning that the hop is rejected and the active state is set to the original adiabatic state |*I*⟩. In addition, we do not modify the nuclear velocity in this case. In the Result section, we compare population dynamics in the basis of |*ψ*_{i}(**R**)⟩ = |*ϕ*_{α}(*R*)⟩ ⊗|*n*⟩, instead of in the adiabatic polariton basis $|EI(R)\u3009$. When computing the population in a representation that is rather than the adiabatic states of $V\u0302$, there is no unique way to calculate them in FSSH. Here, we choose the simplest possible choice of using the electronic expansion coefficients [in Eqs. (15) and (17)] to compute the population dynamics. Other choices as suggested in Ref. 61 might provide more accurate results for FSSH. These approaches will be explored in future studies.

The initial nuclear distribution of MQC simulations are generated by sampling the Wigner density $[\u27e8R|\chi \u27e9]w=1\u210f\pi e\u2212M(P2+\omega 2(R\u2212R0)2)/\omega \u210f$, which is the Wigner transformation of the nuclear wavefunction *χ*(*R*) = ⟨*R*|*χ*⟩ in the initial state [see Eq. (58)]. Here, *R* and *P* are the nuclear coordinate and momentum, respectively. The initial state for the electronic–photonic subsystem is set to be |*e*0⟩. The population dynamics were averaged over 2000 trajectories, although 500 trajectories was good enough to produce the basic trend of the polariton dynamics. The light–matter coupling strength *g*_{c} was chosen to be 0.001, 0.005, and 0.01. For FSSH simulation, we need to choose an initial active state, and the initial state for the polariton subsystem |*e*(*R*), 0⟩ is not one of the eigenstate $|E(R)\u3009$. We follow the previous work^{62} and use the canonical 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$.

### C. Exact polariton quantum dynamics

To benchmark the performance of the MQC polariton dynamics simulations, we also solve the exact quantum dynamics for the molecule-cavity hybrid system. We express the total wavefunction of the hybrid system through the Born–Huang expansion as follows:

where *χ*_{νn}(**R**) is the nuclear wavefunction. In our cases, electronic |*ϕ*_{ν}(**R**)⟩ and photonic |*n*⟩ basis states are truncated to some finite number of states. For example, in the four-level basis, *ν* = *g*, *e*, and *n* = 0, 1, and the total wavefunction can be written out explicitly as

To simplify our notation, we omit **R**-dependence of the adiabatic electronic states. Plugging the above expansion to time-independent Schödinger equation (TISE) $H\u0302|\Psi \u3009=\epsilon |\Psi \u3009$, we can obtain the eigenvalue equation for the nuclear wavefunction [see Eq. (70)]. The main obstacle is how to write the nuclear kinetic energy operator correctly in the dressed states. We discard the nuclear index *α* because we only have one nucleus in our model system (i.e., ∇_{α} = ∇). The nuclear gradient operator is applied to the total wavefunction [Eq. (67)], and we arrive at the following two equations:

Multiply ⟨*g*0|, ⟨*e*0|, ⟨*g*1|, and ⟨*e*1| to the above equation, respectively, we can rewrite the above equation in matrix form

The matrix in right hand side can be seen as the nuclear kinetic operator in the adiabatic-Fock states basis except a constant −*ℏ*^{2}/2*M*. Since we already derived the formula of $V\u0302$ in the same basis in Eq. (37), we can directly plug these expressions into the TISE and obtain the following eigenequation:

We use the Sinc DVR basis^{56} to represent the nuclear wavefunction {*χ*} 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

where *ɛ*_{j} is the *j*_{th} eigenvalue and *c*_{j} is the projection of initial total wavefunction onto the *j*_{th} eigenstate |*ɛ*_{j}⟩,

where |Ψ(*t* = 0)⟩ = |Φ⟩ defined in Eq. (58).

The potential technical challenge is to obtain the second order derivative coupling, e.g., ⟨*g*|∇^{2}*g*⟩ and ⟨*g*|∇^{2}*e*⟩, terms in Eq. (69), which are not always available in electronic structure methods. However, for a two-electronic-state problem, they can be calculated explicitly using **d**_{ge} as

## IV. RESULTS AND DISCUSSIONS

Figure 2(a) presents the first four polaritonic surfaces of the SM1 model coupled to a resonant optical cavity, where the potential is color coded 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^{34,53} because the rigorous photon number operator should be obtained by applying the Power–Zienau–Woolley (PZW) Gauge transformation^{33,65,66} 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.^{67} We use a light–matter coupling strength of *g*_{c} = 0.005 a.u., and the photon energy of the cavity mode is chosen as *ℏω*_{c} = 1.281 eV ($\u22480.047$ a.u.).

Figure 2(b) presents the exact quantum dynamics within the JC subspace {|*e*0⟩, |*g*1⟩} (dashed lines) or in the larger {|*g*0⟩, |*e*0⟩, |*g*1⟩, |*e*1⟩} (solid lines) sub-space, under which the polaritonic eigenstates and dynamics are converged. The initial condition is described in Eq. (58). Clearly, additional states beyond the JC subspace will be explored by the quantum dynamics of the hybrid system due to the NACs among these states. The population of the |*e*1⟩ state is mainly contributed from the population transfer from the |*e*0⟩ state due to the light–matter coupling originated from the permanent dipole *μ*_{ee}. In addition, the population transition between |*e*0⟩ and |*g*1⟩ is mediated by the cavity-induced coupling near *R* = 0, where the PES exhibits a strong mixing between these two states as shown in Fig. 2(a) at R = 0.

Figure 2(c) presents several components of the **X**_{ij} matrix in the general gradient expression [Eq. (40)] in the subspace of {|*g*0⟩, |*e*0⟩, |*g*1⟩, |*e*1⟩}. Comparing the gradient expression in the two-state [Eq. (36)] and four-state subspaces [Eqs. (39) and (40)], the [**X**] terms [Eq. (40)] do not show up in the two-level JC case [Eq. (36)]. The JC gradient component $\u2207Ve0,g1=gc\u03f5\u22c5\u2207\mu ge$ and the regular gradient from the electronic derivative couplings $Xg0,e0=dge(Ve\u2212Vg+Dee2\u2212Dgg2)$ are compared with the rest of the non-JC type of gradients, such as $Xg0,g0=2dgeDeg2$ (where $Deg2$ is the matrix of the DSE operator), $Xg0,e1=dgegc\u03f5\u22c5\mu ee\u2212\mu gg$, and **X**_{g0,g1} = 2**d**_{ge}*g*_{c}** ϵ** ⋅

*μ*_{ge}. As we show in Fig. 2(c), all of the gradients we discussed above have a similar magnitude in this Shin–Metiu-cavity model with a coupling strength

*g*

_{c}= 0.005 a.u. Some elements can be even larger than ∇[

*V*]

_{e0,g1}, such as

**X**

_{g0,e0}and

**X**

_{g0,g1}. All of these non-JC type of gradient terms, unfortunately, are not included in the previous MQC studies of polariton quantum dynamics based upon the JC models.

^{22,23,42}When the system starts to explore all of the states and generate sizable population and coherences [$\rho ij=ci*cj$, where

*c*

_{i}is the expansion coefficients in Eq. (17)], the nuclear forces [for example, in Eq. (24)] from these new gradient terms will have a non-negligible contribution that is required to be explicitly and correctly included.

Figure 2(d) presents the polariton dynamics obtained from the Ehrenfest MQC simulations (dotted lines) in the four-level subspace, compared to the numerically exact results (solid lines). It shows that the population dynamics of MQC agrees reasonably well with the exact results. For the two-level case, the Ehrenfest MQC simulation results are also consistent with the exact quantum simulation shown in Fig. 2(b) (dashed line), which are not presented here.

Figure 3 presents the MQC population dynamics with light–matter coupling strength of *g*_{c} = 0.001 [panels (a), (d), and (g)], *g*_{c} = 0.005 [panels (b), (e), and (h)], and *g*_{c} = 0.01 [panels (c), (f), and (i)]. The polariton population dynamics are obtained from both Ehrenfest dynamics [open circles in panels (d), (e), and (f)] as well as the FSSH approach [dots in panels (g), (h), and (i)], and compared to the numerically exact results (solid lines).

When the coupling is weak (*g*_{c} = 0.001), the polaritonic PESs [panel (a)] have nearly identical curvatures of the photon dressed states |*g*1⟩ and |*e*0⟩. The light-induced mixing only occurs near *R* = 0 between the |*e*0⟩ and |*g*1⟩ states. The population dynamics presented in Figs. 3(d) and 3(g) suggest that the population transfer mainly occurs between the |*g*1⟩ and |*e*0⟩ states. Both Ehrenfest dynamics and the FSSH method provide an accurate description of the population dynamics. When the coupling strength increases, the polaritonic PES starts to change its curvatures and characters, which can be clearly seen from the color coding of the polariton PES in Figs. 3(b) and 3(c). As can be seen from Fig. 3(b) (*g*_{c} = 0.005) and Fig. 3(c) (*g*_{c} = 0.01), the gap between the 2nd and the 3rd polariton PES becomes larger and the state |*e*1⟩ starts to mix with the states within the JC subspace. The |*g*0⟩ state, on the other hand, only slightly mixes with the other photon dressed states, even with the largest light–matter couplings we considered.

The population dynamics for *g*_{c} = 0.005 are presented in Figs. 3(e) and 3(h). The main transitions occur between the |*e*0⟩ and |*g*1⟩ states, and the |*e*1⟩ state also shows a finite oscillating population, which is mainly caused by the light–matter interaction carried by the permanent dipole *μ*_{ee}. When increasing the coupling strength to *g*_{c} = 0.01, the main population transition now happens between |*e*0⟩ and |*e*1⟩ states, due to the permanent dipole *μ*_{ee} that dominates the light–matter coupling. For all cases, both Ehrenfest and FSSH approaches generate reasonably accurate polariton population dynamics, with FSSH slightly outperforming the Ehrenfest dynamics.

Figure 4 presents the polariton dynamics of the SM2 model coupled to a cavity.^{40} This model has been used to investigate how cavity can influence proton-coupled electron transfer reaction with the exact factorization approach.^{40,68,69} The PES of SM2 is asymmetric and there is an avoided crossing near *R* = 2.0 a.u. The photon energy is *ℏω*_{c} = 0.1 a.u., causing a light-induced avoid crossing near *R* = −4.0 a.u. [see Fig. 4(a)]. Figure 4 presents the MQC population dynamics with the light–matter coupling strength of *g*_{c} = 0.001 [panels (a), (d), and (g)], *g*_{c} = 0.005 [panels (b), (e), and (h)], and *g*_{c} = 0.01 [panels (c), (f), and (i)]. The polariton population dynamics are obtained from both Ehrenfest dynamics [open circles in panels (d), (e), and (f)] as well as the FSSH approach [dots in panels (g), (h), and (i)] and compared to the numerically exact results (solid lines).

When the coupling is *g*_{c} = 0.001, the polaritonic PESs [Fig. 4(a)] have nearly identical curvatures of the photon dressed states |*g*1⟩ and |*e*0⟩ (which is indicated by the color coding of the surfaces), except at *R* ≈ − 4 a.u., where the |*e*, 0⟩ and |*g*, 1⟩ states mix. In Figs. 4(d) and 4(g), even for a relatively small light–matter coupling strength *g*_{c} = 0.001, the system will have a finite population for |*g*0⟩ state. This is different compared to the case of SM1 [Figs. 3(d) and 3(g)], where the dynamics is largely confined with in the JC subspace {|*g*1⟩, |*e*0⟩} under the same coupling strength. Although SM2 has a large permanent dipole, the population of |*g*0⟩ only starts to grow later in time (20 fs) after the initial excitation. Thus, the population transfer to |*g*0⟩ and |*e*1⟩ states are mainly due to the electronic NAC *d*_{eg} that directly couples the |*g*1⟩ state to |*e*1⟩ state, as well as the |*e*0⟩ state to |*g*0⟩ state. Note that the |*e*1⟩ state does not have a significant population (see solid lines for the exact results); however, both Ehrenfest dynamics [open circles in Fig. 4(d)] and FSSH approach [dots in Fig. 4(g)] incorrectly predict a larger |*e*1⟩ population due to their less accurate MQC approximation.

When the light–matter coupling strength increases, the polaritonic PES starts to change its curvatures and characters, which can be clearly seen from the color coding of the polariton PES. The population of |*e*1⟩ starts to grow at an earlier time, due to the permanent dipole *μ*_{ee} that couples |*e*0⟩ state to the |*e*1⟩ state, whereas the electronic NAC *d*_{eg} still contribute to the later time population transfer to the |*g*0⟩ state. Through out all coupling strengths investigated here, both Ehrenfest and FSSH approaches provide a reasonable accuracy for the population dynamics for the |*g*1⟩ and |*e*0⟩ states at a short time, with a less satisfied accuracy at a longer time compared to the case of SM1 coupled to the cavity in Fig. 3.

The canonical picture of the JC model for cavity QED has provided invaluable insights into the field of polariton physics and chemistry. In this model, the cavity photons dress the molecular energy level by the energy of one photon, providing the photon-dressed ground state |*g*1⟩ state, which can be resonant with the |*e*0⟩ state and then hybridize with it. The transition dipole *μ*_{eg} causes the superposition between |*e*0⟩ and |*g*1⟩ states [through the light–matter coupling term in Eq. (9)]. The NAC between polaritonic PES *d*_{12} [Eq. (51)] will be large due to the near degeneracy of the polariton states. When the NAC from other states are smaller than *d*_{12} [Eq. (51)], the polariton dynamics is largely confined within the Hilbert subspace {|*e*0⟩, |*g*1⟩} states. However, this approximation will breakdown under several scenarios, and we just enumerate two commonly encountered ones that we have already explored in our numerical results. First, the electronic NAC *d*_{eg} [Eq. (19)] between |*e*0⟩ and |*g*0⟩ could play an important role in the dynamics (as shown in Fig. 4) and could have a large contribution in the [**X**] expression [Eq. (40)]. Note that the effect of electronic NAC *d*_{eg} is not included in the JC subspace [see Eq. (35)]. Correctly including the intrinsic NAC is essential to investigate molecular cavity QED. The second scenario is when molecules have large permanent dipole moments. In our case, both SM1 and SM2 have a large permanent dipole *μ*_{ee} that causes the coupling between |*e*0⟩ and |*e*1⟩ through the light–matter coupling in Eq. (5), and this coupling is proportional to the light–matter coupling strength *g*_{c}. Thus, with an increasing *g*_{c}, the |*e*1⟩ state could get a larger population (as shown in Fig. 3). Unlike the atomic cavity QED, the dynamical interplay among the electronic, photonic, and nuclear DOFs are often more complicated; thus, the truncation to any Hilbert subspace (including the JC subspace) must be done with caution.^{33,34,41}

## V. CONCLUSIONS

In this paper, we derive a rigorous expression of the nuclear gradient using the QED Hamiltonian without making unnecessary approximations. Under the adiabatic-Fock representation, and requiring the total energy conservation of this MQC system, we derived the rigorous nuclear gradient, and it is intuitively connected to those approximate gradients under the JC approximation.^{23} This rigorous expression has additional terms that go beyond the expression of the JC Hamiltonian.^{23} We have numerically demonstrated the importance of these terms in the MQC simulations of polariton dynamics.

The nuclear gradient expression can be readily used in any MQC simulations, such as Ehrenfest and FSSH approaches demonstrated in this work, as well as nuclear wavepacket approaches that require a nuclear gradient.^{70,71} The theoretical developments in this work could significantly benefit the polariton quantum dynamics community with a rigorous nuclear gradient of the molecule-cavity hybrid system and could have a broad impact on the future non-adiabatic simulations of 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). W.Z. appreciates the support from the China Scholarship Council (CSC) during his visit at the University of Rochester. Computing resources were provided by the Center for Integrated Research Computing (CIRC) at the University of Rochester. The authors appreciate valuable discussions and comments from Braden M. Weight.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Wanghuai Zhou**: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (equal); Software (lead); Validation (equal); Visualization (equal); Writing – original draft (lead). **Deping Hu**: Data curation (supporting); Formal analysis (supporting); Validation (equal); Writing – original draft (supporting); Writing – review & editing (equal). **Arkajit Mandal**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Pengfei Huo**: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (lead); Resources (lead); Software (supporting); Supervision (lead); Validation (equal); Visualization (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 reasonable request.

### APPENDIX A: PROOF OF EQ. (23)

Here, we provide a simple proof of Eq. (23) that was used to derive the nuclear gradient expression in Eq. (24). Starting from the equation of coefficients Eq. (20), we have

Combine the above two equations together, we have

Using the above result, the left-hand side of Eq. (23) becomes

The order of summation over the index *i*, *j*, *k* does not affect the final result. By relabeling the index *i* → *l*, *k* → *j*, and *j* → *i*, we can see that

Thus, the second term of the right-hand side of Eq. (A3) is equal to zero. Similarly, by relabeling the index *l* → *i*, *j* → *k*, and *i* → *j*, one can show that

Plugging the above equation into Eq. (A3), we have

which is Eq. (23) of the main text. In the last step, we have used the property $dij\alpha =\u2212dji\alpha $.

### APPENDIX B: NUCLEAR GRADIENTS OF MOLECULES

For the case of an isolated molecular system, it is easy to check that the nuclear forces in Eq. (24) reduce to the familiar expressions. For a diabatic basis, the derivative coupling between states is strictly zero, i.e., **d**_{ij} = 0, and Eq. (24) becomes the Ehrenfest equation of motion

For the adiabatic basis, the electronic potential matrix [*V*] now is diagonal, i.e., *V*_{ij} = *E*_{i}(**R**)*δ*_{ij}, where *E*_{i} is *i*_{th} eigenvalue satisfying $V\u0302|\varphi i\u3009=Ei|\varphi i\u3009$. Under this case, Eq. (24) becomes

where we have used the chain rule $\u27e8\varphi i|\u2207\alpha V\u0302|\varphi j\u27e9=\u2207\alpha \u27e8\varphi i|V\u0302|\varphi j\u27e9\u2212\u27e8\u2207\alpha \varphi i|V\u0302|\varphi j\u27e9\u2212\u27e8\varphi i|V\u0302|\u2207\alpha \varphi j\u27e9$ as well as the following well-known equality:^{18}

### APPENDIX C: ADDITIONAL NUMERICAL RESULTS

#### 1. Convergence of the Fock states

All simulations presented in the main text are restricted in the Fock state basis |0⟩ and |1⟩. Under the strong light–matter interactions, more Fock states are required to reach convergence.^{9,34} In Fig. 5, we provide the convergence test of the Fock states with the intermediate coupling case (*g*_{c} = 0.005) for SM1.

Figure 5(a) presents the population dynamics of the |e0⟩ state with the exact quantum dynamics. The results obtained using *n* = 6 and *n* = 8 Fock states are nearly identical and only display a small difference after 30 fs, and the *n* = 10 Fock state results are identical to the *n* = 8 result (which is on top of the *n* = 10 result). Figure 5(b) presents the same calculations using the Ehrenfest dynamics, with the general nuclear gradient expression in Eqs. (29) and (31). As we already presented in the main text, the Ehrenfest simulation results (or FSSH results) are very close to the exact simulation results for the SM1-cavity coupling model investigated here, the Ehrenfest dynamics thus presents almost identical trend of the convergence for the population dynamics.

#### 2. Importance of the dipole self-energy

We demonstrate the importance of the DSE terms in polariton quantum dynamics simulations, especially in the strong and ultra-strong light–matter coupling regimes.

Figure 6 presents the value of the DSE term $Dge2$ along the nuclear coordinate (proton coordinate) for the SM1 model coupled to a cavity with different light–matter coupling strengths. For comparison, we also present other off-diagonal light–matter coupling terms, such as *g*_{c}** ϵ** ⋅

*μ*_{gg}(red),

*g*

_{c}

**⋅**

*ϵ*

*μ*_{ge}(green), and

*g*

_{c}

**⋅**

*ϵ*

*μ*_{ee}(blue), which are the matrix elements of the light–matter coupling operator [see Eq. (37)]. Note that, according to the definition of $Dge2$ in Eq. 31(d), $Dge2=gc2\u210f\omega c(\u03f5\u22c5\mu gg)(\u03f5\u22c5\mu ge)+(\u03f5\u22c5\mu ge)(\u03f5\u22c5\mu ee)$, which depends on both the magnitude of the dipoles as well as the light–matter coupling strength

*g*

_{c}.

Figure 6 clearly shows that by increasing light–matter coupling strength, the magnitude of DSE term $Dge2$ (purple lines) becomes comparable with light–matter coupling terms and thus is non-negligible in the QED simulations. Other DSE terms, such as $Dgg2$ and $Dee2$, have a similar behavior.

Figure 7 presents the population dynamics of the SM1 model coupled to a cavity, obtained from the numerically exact approach [panels (a), (c), and (e)] at different coupling strengths. For comparison, the results obtained based on the PF Hamiltonian with the presence of the DSE terms (solid lines) and the dynamics without all DSE terms (open circles) are presented. One can clearly see that when the light–matter coupling strength is relatively weak (*g*_{c} = 0.001), the population dynamics is not significantly influenced by omitting the DSE term [Fig. 7(a)]. At an intermediate coupling strength (*g*_{c} = 0.005), or a large coupling strengths (*g*_{c} = 0.01), the discrepancy of the population dynamics becomes larger [Fig. 7(c)], and the population dynamics stars to show different features [Fig. 7(e)] if DSE terms are not included.

We have further performed the Ehrenfest dynamics [panels (b), (d), and (f)] based on the PF Hamiltonian with the DSE terms (solid lines) or exclude the presence of the DSE (open circles) for both couplings [*V*] as well as in the nuclear gradient terms [∇_{α}*V*]. The results obtained from the Ehrenfest dynamics simulations are in a good agreement with the numerically exact simulations (panels a, c, e). This demonstrates the accuracy of the Ehrenfest simulations for this particular model system, as well as the importance of the DSE terms for an accurate description of the light–matter interactions.