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.

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 approaches14–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 balance24 (incorrect long time population) in Ehrenfest dynamics and the creation of artificial electronic coherence25 or incorrect chemical kinetics25 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 matrix29,30 (PLDM), symmetrical quasi-classical31,32 (SQC), and the quantum-classical Liouville equation (QCLE) dynamics.33,34 In particular, the recently developed γ-SQC has been shown35 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 representation37 for the electronic DOF as well as the linearization approximation39,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 systems36 and many-state exciton Hamiltonians of light-harvesting complexes37). The γ-SQC approach has already demonstrated41 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 scheme41,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 analytically43 and numerically45,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) model48,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 γ-SQC35 and spin-LSC36,37 to obtain non-adiabatic polariton dynamics, which outperforms widely used MQC approaches.

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

(1)

where T̂n represents the nuclear kinetic energy operator, and Ĥen is the electronic Hamiltonian that describes electron–nucleus interactions. Furthermore, Ĥp, Ĥenp, and Ĥd represent the photonic Hamiltonian, electronic–nuclear–photonic interactions, and the dipole self-energy (DSE) term, respectively. A full derivation of this Hamiltonian, as well as its connection with the various atomic cavity QED models, can be found in the appendix of Ref. 42.

The electronic–nuclear potential Ĥen, which describes the common molecular Hamiltonian excluding the nuclear kinetic energy, is described as follows:

(2)

The above expression includes electronic kinetic energy, electron–electron interaction, electron–nucleus interaction, and nucleus–nucleus interaction. The expressions of these four terms can be found in previous work.50–52 Modern electronic structure theory have been developed around solving the eigenvalue problem of Ĥen, providing the following electronically adiabatic energy and its corresponding state:

(3)

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

(4)

where q̂c=/2ωc(â+â) and p̂c=iωc/2(ââ) are photon field operators, â and â 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

(5)

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 λ=1/ϵ0V0, 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=ωc/2λ. Note that the common notation used in the literature,6,53 the definition of gc also includes λμ̂. Furthermore, the total dipole operator of both electrons and nuclei is defined as

(6)

where −e is the charge of the electron and Zje is the charge of the jth nucleus.

Finally, the DSE term is expressed as

(7)

This is a necessary term in the PF Hamiltonian and ensures both gauge invariance of the Hamiltonian9,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

(8)

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

(9)

where T̂n represents the nuclear kinetic energy operator, and V̂ represents the rest of the Hamiltonian. For a bare molecular system, V̂=Ĥen expressed in Eq. (2). For a molecule–cavity hybrid system,

(10)

which is commonly referred to as the polariton Hamiltonian,3,58 also denoted as Ĥpl. In a similar way that electronic adiabatic states are defined in Eq. (3), one can further define the polaritonic state3,58 as the eigenstate of V̂=Ĥpl [see definition in Eq. (10)] through the following eigenequation:

(11)

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

(12)

where cα,nJ(R)=ϕα(R),n|EJ(R) and EJ(R) can be obtained by diagonalizing the matrix of V̂=Ĥpl [constructed from the adiabatic-Fock state basis in Eq. (8)] as

(13)

where

(14)

Note that the R-dependence of |EJ(R) 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) is the eigenstate of V̂, whereas the adiabatic state |ϕα(R)⟩ is only the eigenstate of Ĥen, and not for V̂.

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) [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 ∈ [t0, t1], where the nuclear positions evolve from R(t0) to R(t1), and the corresponding adiabatic-Fock basis [defined in Eq. (8)] are {|ψi(R(t0))⟩} and {|ψjR(t1)}. We uses the basis {|ψi(R0)⟩ ≡|ϕα(R0), n⟩} at the reference nuclear geometry R(t0) as the diabatic basis during this short-time propagation such that

(15)

With the above QD basis defined independently of R(t) within each propagation segment, the electronic derivative couplings vanish while V̂(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 ∈ [t0, t1].

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

(16)

For on-the-fly simulations, this quantity is obtained from a linear interpolation59 between Vαβ,mn(R0) and Vαβ,mn(R(t1)) as follows:

(17)

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(R0) and Vαβ,mn(R(t1)) separately for the molecule–cavity hybrid system, as discussed below. Using electronic ab initio calculation, as well as the properties of â and â for the photonic DOF, we can explicitly evaluate each term of Vαβ,mn(R0) [see Eq. (10)] as follows:

(18a)
(18b)
(18c)
(18d)

where Ĥp [see its definition in Eq. (4)] is an R-independent operator, the sum γ in the matrix element of Ĥd runs over the diabatic states, and Dαβ2 denotes the elements of DSE. Furthermore, the matrix element of the dipole operator under the diabatic representation is expressed as

(19)

Similarly, at time t1, the matrix element Vαβ,mn(R(t1))=ϕα(R0),m|V̂(R(t1))|ϕβ(R0),n can also be written explicitly, with each term expressed as follows:

(20a)
(20b)
(20c)
(20d)

where Hαβen(R(t1))ϕα(R0)|Ĥen(R(t1))|ϕβ(R0), and μαβ(R(t1))ϕα(R0)|μ̂(R(t1))|ϕβ(R0).

To conveniently calculate Hαβen(R(t1)) and μαβ(R(t1)), we use the following relations:

(21a)
(21b)

where the matrix elements at R(t1) are expressed as

(22a)
(22b)

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

(23a)
(23b)

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(t1)). First, let us focus on the gradient term from the electron–nuclear coupling term47 as follows:

(24)

Here, from the first to the second line, we have used the fact that neither ⟨ϕα(R0)| nor ⟨m| are R-dependent, which allows ∇ to bypass both and directly act on V̂(R(t1)). We have also used the fact that Ĥen does not contain any photonic operators. The gradient term ϕλ(R(t1))|Ĥen(R(t1))|ϕν(R(t1)) can be evaluated using the following well-known equality:18 

(25)

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

(26)

For the gradient on the matrix Hαβ,mnp, because there is no nuclear DOF in Ĥp, thus

(27)

For the gradient on the light–matter interaction term Hαβ,mnenp, we have

(28)

where Sαλ and Sβν are defined in Eqs. (23a) and (23b), respectively. To evaluate the term ϕλ(R)|μ̂(R)|ϕν(R) that appears in Eq. (28), we use a simple relation based on the chain rule as follows:

(29)

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̂=γ|ϕγ(R)ϕγ(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:

(30)

where the term ϕα(R(t0))|μ̂(R(t1))|ϕγ(R(t0)) can be evaluated using Eq. (21b), and the ϕα(R0)|μ̂(R(t1))|ϕβ(R0) 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 trajectories62–64 to propagate the quantum dynamics in the time step dt = t1t0 for t ∈ [t0, t1]. During the next short-time propagation segment t ∈ [t1, t2], the QD scheme adopts a new reference geometry R0R(t1) and new diabatic basis |ψj(R0)|ψj(R(t1))=|ϕβ(R(t1)),m. Between the t ∈ [t0, t1] propagation and the t ∈ [t1, t2] propagation segments, all of the necessary quantities will be transformed from {|ψi(R0)⟩} to the {|ψj(R0)} basis using the relation

(31)

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

(32)

where the electronic adiabatic state overlaps ⟨ϕα(R(t0))|ϕβ(R(t1))⟩ 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(t0))⟩} or {|ψj(R(t1))⟩}. 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(t0))|ψj(R(t1))⟩, as been done in our previous work44 for simulating photo-induced charge transfer dynamics. Here, we perform the Löwdin orthogonalization procedure68 as commonly used in the local diabatization approach66 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 ϕβ(R(t))|tϕα(R(t))=dβα(R)Ṙ, 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(t1))|ϕα(R(t2))⟩ instead of ϕβ(R(t))|tϕα(R(t)). 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(t1))|ϕα(R(t2))⟩.

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(R0)⟩} 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 R0 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.

The Meyer–Miller–Stock–Thoss (MMST) formalism26–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

(33)

where Vab(R̂)=a|V̂(r̂,R̂)|b 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 R0 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:

(34)

where 2γb is viewed as a parameter69 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)

(35a)
(35b)

with the nuclear force expressed as

(36)

Overall, the MMST mapping provides a consistent classical footing for both electronic and nuclear DOFs, and the non-adiabatic transitions between electronic states are captured through the classical motion of the fictitious harmonic oscillators. The non-adiabatic dynamics obtained from this formalism have shown good performance in the ab initio on-the-fly dynamics.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 {pb, qb} through

(37)

and the inverse relations

(38)

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 as69 

(39)

where ρ̂(0)=ρ̂R|aa| is the initial density operator, ρW(P, R) is the Wigner transform of ρ̂R operator for the nuclear DOFs, ε={ε1,ε2,,εN} is the positive-definite action variable vector for N electronic states,35,Wa(ɛ) = δ(ɛa − (1 + γa))abδ(ɛbγb) is the Wigner transformed action variables,72 and dτdP · dR · dɛ · dθ. For practical reasons, the above delta functions in Wa(ɛ) 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

(40)

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

(41a)
(41b)

The parameter Γ is related to the radius of the generalized Bloch sphere rs through Γ=2N(rs1), where s and s̄ are complementary indices in the Stratonovich–Weyl transform. Among the vast parameter space, one of the best performing choices36,37 is when rs=rs̄=N+1, which is referred to as s = W, leading to a ZPE parameter

(42)

as well as the identical expression of [|a⟩⟨a|]s and [|bb|]s̄ in Eq. (41). We further use the focused initial condition36,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Γ)=1 for initially occupied state |a⟩ and 12(qb2+pb2Γ)=0 for the initially unoccupied states |b⟩. The angle variables {θb} [Eq. (37)] are randomly sampled37 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 ψi(R0)|V̂(R)|ψj(R0) [see Eqs. (16) and (17)] and nuclear gradient ψi(R0)|V̂(R)|ψj(R0) [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(t0))⟩ ≡ |ϕα(R(t0)), n⟩} to {|ψj(R(t1))⟩ ≡ |ϕβ(R(t1)), m⟩}. This leads to the corresponding transform of mapping variables between the two consecutive QD bases as follows:43,47

(43a)
(43b)

where the overlaps between the two steps are evaluated using Eq. (32) (see discussion under that equation). More computational details for the γ-SQC and spin-LSC are provided in Sec. III B.

In this work, we use the asymmetrical Shin–Metiu model48,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 Ĥen [cf. Eq. (2)] is expressed as

(44)

where V̂p(R) represents the potential of the transferring proton, V̂e(r) represents the potential of the transferring electron, and V̂ep(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., af = 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) basis75 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.

FIG. 1.

(a) The schematic illustration of the asymmetrical SM model, where one proton and one electron can transfer between two fixed ions (donor and acceptor). The distance between the donor and acceptor is 19 a.u. Here, Vp and Ve [see Eq. (44)] are the potential of the transferring proton and electron, respectively. (b) The first two adiabatic electronic states, (c) NAC between these two electronic states, and (d) their permanent and transition dipole moments of the SM model. The PESs of the polaritonic states inside the cavity are obtained with the light–matter coupling strength (e) gc = 0.001 and (f) gc = 0.005. The color used in (e) and (f) is coded according to ââ, as shown in the upper position of panel (e).

FIG. 1.

(a) The schematic illustration of the asymmetrical SM model, where one proton and one electron can transfer between two fixed ions (donor and acceptor). The distance between the donor and acceptor is 19 a.u. Here, Vp and Ve [see Eq. (44)] are the potential of the transferring proton and electron, respectively. (b) The first two adiabatic electronic states, (c) NAC between these two electronic states, and (d) their permanent and transition dipole moments of the SM model. The PESs of the polaritonic states inside the cavity are obtained with the light–matter coupling strength (e) gc = 0.001 and (f) gc = 0.005. The color used in (e) and (f) is coded according to ââ, as shown in the upper position of panel (e).

Close modal

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 (0.1 a.u). Furthermore, we assume that the cavity field polarization direction ϵ is always aligned with the direction of the dipole operator μ̂ such that ϵμγν = μγν (for {ν, γ} = {e, g}) where μγν is the magnitude of μ̂. Explicitly considering the angle between ϵ and μ̂ 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 gc = 0.001 a.u. and gc = 0.005 a.u. in this work. The normalized coupling strength is often defined as77ηgc   |ϵμ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 gc = 0.001 a.u.) and η = 0.3 (for gc = 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) [see their definition in Eq. (11)] are presented in Fig. 1(e) with the light–matter coupling strength gc = 0.001 a.u. and in Fig. 1(f) with the light–matter coupling strength gc = 0.005. These polariton potentials are color coded [as shown in the inset of panel (e)] based on the expectation value of ââ 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 Hamiltonian56,81 because the rigorous photon number operator should be obtained by applying the Power–Zienau–Woolley (PZW) gauge transformation54,82,83 on the photon number operator ââ. 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

(45)

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 χ(R)=R|χexp[Mω0(RR0)2/2], where M is the mass of the proton (nucleus in the SM model), and R0 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 R0 under the harmonic approximation, with the harmonic oscillation frequency being ω0. We use the parameters in the original Ref. 49 for R0 = −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.

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 {pb, qb} 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 conditions37 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

(46)

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 |e0⟩. 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 variables85,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 gc 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 

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 gc = 0.001 a.u. The system is initially prepared in the |e0⟩ state and decays quickly into the |g1⟩ state during the first 12 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 |e0⟩ and |g1⟩ until about 20 fs. All of the MQC and mapping dynamics methods can describe the above process reasonably well compared to the exact results. After that, all the dynamics results (including the exact one) show a fast population increase of the |g0⟩ state, which is due to the electronic NAC deg that directly couples the |e0⟩ state to |g0⟩ 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 |g0⟩ population as well as its long time plateau. Moreover, both Ehrenfest and FSSH dynamics predict a significant population transfer from |g1⟩ to |e1⟩ 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 |g1⟩ to |e1⟩ 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

FIG. 2.

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

FIG. 2.

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

Close modal

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

FIG. 3.

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

FIG. 3.

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

Close modal

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 (gc = 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 (gc = 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 (|g0⟩, |e0⟩, |g1⟩, and |e1⟩), we can also see some other states (|h0⟩, |f0⟩, |g2⟩, and |f1⟩) 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 |g0⟩ 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.

FIG. 4.

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

FIG. 4.

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

Close modal

In this work, we generalize the quasi-diabatic (QD) propagation scheme44,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, γ-SQC35 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.

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.

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.

The authors have no conflicts to disclose.

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).

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

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 window35,72 is expressed as

(A1)

where the window functions are defined as

(A2)

and

(A3)

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 normalized31 with the following procedure:

(A4)

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 potentials35 and has been combined with the kinematic momentum approach92 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⟩,

(A5)

or equivalently,

(A6)

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

(A7)

ensuring the initial forces (at t = 0) are simply F = −∇Vaa(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.

To calculate the electronic properties of the SM model, we use the Sinc DVR basis75 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 Ĥel in this grid basis {|ri⟩} are given by

(B1)

where the ri|T̂r|rj is given analytically75 as follows:

(B2)

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

(B3)

where ciα(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:

(B4)

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

(B5)

Using |ϕα(R)⟩ and |ϕβ(R)⟩ in the grid basis, we can directly evaluate the nuclear gradient ϕλ|Ĥen|ϕν [Eq. (25)] as follows:

(B6)

where the Ĥen(R) is evaluated analytically using the expression of Ĥen(R) in Eq. (44). This gives the adiabatic gradient ∇Eλ (for λ = ν) and derivative coupling (for λν) as dλν=ϕλ|Ĥen|ϕν/(EνEλ), as indicated in Eq. (44). The nuclear gradient ϕλ|Ĥen|ϕν 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 ciγ(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:

(B7)

To solve the exact quantum dynamics, we represent the total wavefunction of the hybrid system as |Ψξ=i,kcikξ|ψi(Rk)|Rk, where {|Rk⟩} is the DVR grid basis for the nucleus, and the |ψi(Rk)⟩ = |ϕα(Rk)⟩ ⊗|n⟩ is the adiabatic-Fock basis [Eq. (8)], where the electronic adiabatic basis |ϕα(Rk)⟩ is obtained by solving the electronic eigenequation (Eq. (3)) using the DVR basis for the electronic DOF at the nuclear configuration Rk. The coefficients for the total wavefunction cikξ and the eigenvalue of the total Hamiltonian Ĥ [Eq. (1)] will be obtained by solving the time-independent Schödinger equation Ĥ|Ψξ=Eξ|Ψξ. We use the Sinc DVR basis75 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 |Φ(t)=ξCξexpiEξt|Ψξ, 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 

Besides the mapping dynamics methods (γ-SQC and spin-LSC), we also apply the commonly used Ehrenfest and Tully’s FSSH17,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

(C1)

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)

(C2)

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)

(C3)

where c is the transpose of the coefficient column vector c expressed as follows:

(C4)

The nuclear gradient matrix is expressed as

(C5)

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

In the FSSH dynamics, we expand the time-dependent wave function in the polaritonic basis |EI(R(t)) [see definition in Eq. (11)]

(C6)

where cI 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)) [eigenstate of V̂, see Eq. (11)] as follows:

(C7)

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

(C8)

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 to any other polariton state |EJ during the time interval between t and t + δt is

(C9)

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

(C10)

Since the probability should be positive definite, one sets14 fIJ to 0 if fIJ < 0. The non-adiabatic transition, i.e., stochastic switch from the currently occupied state |EI to another state |EK, occurs if the following condition is satisfied:

(C11)

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, while the velocities of the nuclei are rescaled along the direction of the NAC vector dIK(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 ci(0) for the state |ψi⟩ = |e0⟩ is set to be one, and the rest of the coefficients are set to be zero. These {cj(0)} can be unitary transformed into the coefficients {cI(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). We thus follow the previous work88,93 and use the Monte Carlo scheme to randomly choose the initial active state |EI(R) for each trajectory, based on the magnitude of |EI(R)|e(R),0|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̂, 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 cI(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 ρijaf(t), and ρ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 ρiiaf from the FSSH simulation, the most straightforward way (as the results presented in the main text) is through following unitary transformation:

(C12)

where [ρplRl(t)] is the reduced density matrix in the polariton basis along a given nuclear trajectory Rl(t), with l as the label of the trajectory, and the elements as

(C13)

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

(C14)

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 choices88 to calculate the populations that are not in an adiabatic representation. These methods vary on how to calculate the polaritonic state density matrix [ρpl(Rl(t))]. The first choice is based on the active state index and ignores the polaritonic state coefficients {cI(t)} and the density matrix elements are written as

(C15)

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 choice88 (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 {cI(t)}

(C16)

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.

FIG. 5.

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

FIG. 5.

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

Close modal
1.
2.
M.
Kowalewski
and
S.
Mukamel
,
Proc. Natl. Acad. Sci. U. S. A.
114
,
3278
(
2017
).
3.
J.
Flick
,
M.
Ruggenthaler
,
H.
Appel
, and
A.
Rubio
,
Proc. Natl. Acad. Sci. U. S. A.
114
,
3026
(
2017
).
4.
R. F.
Ribeiro
,
L. A.
Martínez-Martínez
,
M.
Du
,
J.
Campos-Gonzalez-Angulo
, and
J.
Yuen-Zhou
,
Chem. Sci.
9
,
6325
(
2018
).
5.
J.
Feist
,
J.
Galego
, and
F. J.
Garcia-Vidal
,
ACS Photonics
5
,
205
(
2018
).
6.
A.
Mandal
and
P.
Huo
,
J. Phys. Chem. Lett.
10
,
5519
(
2019
).
7.
J. A.
Hutchison
,
T.
Schwartz
,
C.
Genet
,
E.
Devaux
, and
T. W.
Ebbesen
,
Angew. Chem., Int. Ed.
51
,
1592
(
2012
).
8.
A.
Thomas
,
L.
Lethuillier-Karl
,
K.
Nagarajan
,
R. M. A.
Vergauwe
,
J.
George
,
T.
Chervy
,
A.
Shalabney
,
E.
Devaux
,
C.
Genet
,
J.
Moran
, and
T. W.
Ebbesen
,
Science
363
,
615
(
2019
).
9.
A.
Mandal
,
T. D.
Krauss
, and
P.
Huo
,
J. Phys. Chem. B
124
,
6321
(
2020
).
10.
H. L.
Luk
,
J.
Feist
,
J. J.
Toppari
, and
G.
Groenhof
,
J. Chem. Theory Comput.
13
,
4324
(
2017
).
11.
G.
Groenhof
and
J. J.
Toppari
,
J. Phys. Chem. Lett.
9
,
4848
(
2018
).
12.
G.
Groenhof
,
C.
Climent
,
J.
Feist
,
D.
Morozov
, and
J. J.
Toppari
,
J. Phys. Chem. Lett.
10
,
5476
(
2019
).
13.
R. H.
Tichauer
,
J.
Feist
, and
G.
Groenhof
,
J. Chem. Phys.
154
,
104112
(
2021
).
14.
J. C.
Tully
,
J. Chem. Phys.
137
,
22A301
(
2012
).
15.
M.
Barbatti
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
1
,
620
(
2011
).
16.
S.
Mai
,
P.
Marquetand
, and
L.
González
,
Int. J. Quantum Chem.
115
,
1215
(
2015
).
17.
J. C.
Tully
,
J. Chem. Phys.
93
,
1061
(
1990
).
18.
S.
Hammes-Schiffer
and
J. C.
Tully
,
J. Chem. Phys.
12
,
4657
(
1994
).
19.
J.
Fregoni
,
G.
Granucci
,
E.
Coccia
,
M.
Persico
, and
S.
Corni
,
Nat. Commun.
9
,
4688
(
2018
).
20.
J.
Fregoni
,
S.
Corni
,
M.
Persico
, and
G.
Granucci
,
J. Comput. Chem.
41
,
2033
(
2020
).
21.
J.
Fregoni
,
G.
Granucci
,
M.
Persico
, and
S.
Corni
,
Chem
6
,
250
(
2020
).
22.
Y.
Zhang
,
T.
Nelson
, and
S.
Tretiak
,
J. Chem. Phys.
151
,
154109
(
2019
).
23.
W.
Zhou
,
D.
Hu
,
A.
Mandal
, and
P.
Huo
,
J. Chem. Phys.
157
,
104118
(
2022
).
24.
P. V.
Parandekar
and
J. C.
Tully
,
J. Chem. Theory Comput.
2
,
229
(
2006
).
25.
J. E.
Subotnik
,
A.
Jain
,
B.
Landry
,
A.
Petit
,
W.
Ouyang
, and
N.
Bellonzi
,
Annu. Rev. Phys. Chem.
67
,
387
(
2016
).
26.
H. D.
Meyer
and
W. H.
Miller
,
J. Chem. Phys.
70
,
3214
(
1979
).
27.
G.
Stock
and
M.
Thoss
,
Phys. Rev. Lett.
78
,
578
(
1997
).
28.
M.
Thoss
and
G.
Stock
,
Phys. Rev. A
59
,
64
(
1999
).
29.
P.
Huo
and
D. F.
Coker
,
J. Chem. Phys.
135
,
201101
(
2011
).
30.
P.
Huo
and
D. F.
Coker
,
Mol. Phys.
110
,
1035
(
2012
).
31.
S. J.
Cotton
and
W. H.
Miller
,
J. Chem. Phys.
139
,
234112
(
2013
).
32.
S. J.
Cotton
and
W. H.
Miller
,
J. Phys. Chem. A
117
,
7190
(
2013
).
33.
C.-Y.
Hsieh
and
R.
Kapral
,
J. Chem. Phys.
137
,
22A507
(
2012
).
34.
C.-Y.
Hsieh
and
R.
Kapral
,
J. Chem. Phys.
138
,
134110
(
2013
).
35.
S. J.
Cotton
and
W. H.
Miller
,
J. Chem. Phys.
150
,
194110
(
2019
).
36.
J. E.
Runeson
and
J. O.
Richardson
,
J. Comput. Phys.
151
,
044119
(
2019
).
37.
J. E.
Runeson
and
J. O.
Richardson
,
J. Comput. Phys.
152
,
084110
(
2020
).
38.
D.
Bossion
,
W.
Ying
,
S. N.
Chowdhury
, and
P.
Huo
,
J. Chem. Phys.
157
,
084105
(
2022
).
39.
X.
Sun
,
H.
Wang
, and
W. H.
Miller
,
J. Chem. Phys.
109
,
7064
(
1998
).
40.
Q.
Shi
and
E.
Geva
,
J. Phys. Chem. A
108
,
6109
6116
(
2004
).
41.
B. M.
Weight
,
A.
Mandal
, and
P.
Huo
,
J. Chem. Phys.
155
,
084106
(
2021
).
42.
S. N.
Chowdhury
,
A.
Mandal
, and
P.
Huo
,
J. Chem. Phys.
154
,
044109
(
2021
).
43.
A.
Mandal
,
S. S.
Yamijala
, and
P.
Huo
,
J. Chem. Theory Comput.
14
,
1828
(
2018
).
44.
A.
Mandal
,
F. A.
Shakib
, and
P.
Huo
,
J. Chem. Phys.
148
,
244102
(
2018
).
45.
A.
Mandal
,
J. S.
Sandoval C
,
F. A.
Shakib
, and
P.
Huo
,
J. Phys. Chem. A
123
,
2470
(
2019
).
46.
J. S.
Sandoval C
,
A.
Mandal
, and
P.
Huo
,
J. Chem. Phys.
149
,
044115
(
2018
).
47.
W.
Zhou
,
A.
Mandal
, and
P.
Huo
,
J. Phys. Chem. Lett.
10
,
7062
(
2019
).
48.
S.
Shin
and
H.
Metiu
,
J. Chem. Phys.
102
,
9285
(
1995
).
49.
N. M.
Hoffmann
,
L.
Lacombe
,
A.
Rubio
, and
N. T.
Maitra
,
J. Chem. Phys.
153
,
104103
(
2020
).
50.
N. L.
Doltsinis
and
D.
Marx
,
J. Chem. Theory Comput.
01
,
319
(
2002
).
51.
C.
Schäfer
,
M.
Ruggenthaler
, and
A.
Rubio
,
Phys. Rev. A
98
,
043801
(
2018
).
52.
D.
Marx
and
J.
Hutter
,
Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods
(
Cambridge University Press
,
2009
), p.
11
.
53.
M.
Kowalewski
,
K.
Bennett
, and
S.
Mukamel
,
J. Phys. Chem. Lett.
7
,
2050
(
2016
).
54.
M. A. D.
Taylor
,
A.
Mandal
,
W.
Zhou
, and
P.
Huo
,
Phys. Rev. Lett.
125
,
123602
(
2020
).
55.
V.
Rokaj
,
D. M.
Welakuh
,
M.
Ruggenthaler
, and
A.
Rubio
,
J. Phys. B: At., Mol. Opt. Phys.
51
,
034005
(
2018
).
56.
C.
Schäfer
,
M.
Ruggenthaler
,
V.
Rokaj
, and
A.
Rubio
,
ACS Photonics
7
,
975
(
2020
).
57.
E. R.
Koessler
,
A.
Mandal
, and
P.
Huo
,
J. Chem. Phys.
157
,
064101
(
2022
).
58.
J.
Flick
,
H.
Appel
,
M.
Ruggenthaler
, and
A.
Rubio
,
J. Chem. Theory Comput.
13
,
1616
(
2017
).
59.
F.
Webster
,
P. J.
Rossky
, and
R. A.
Friesner
,
Comput. Phys. Commun.
63
,
494
(
1991
).
60.
G. A.
Meek
and
B. G.
Levine
,
J. Phys. Chem. Lett.
5
,
2351
(
2014
).
61.
A.
Jain
,
E.
Alguire
, and
J. E.
Subotnik
,
J. Chem. Theory Comput.
12
,
5256
(
2016
).
62.
L.
Joubert-Doriol
and
A. F.
Izmaylov
,
J. Chem. Phys.
148
,
114102
(
2018
).
63.
M.
Ben-Nun
and
T. J.
Martínez
,
J. Chem. Phys.
108
,
7244
(
1998
).
64.
K.
Saita
and
D. V.
Shalashilin
,
J. Chem. Phys.
137
,
22A506
(
2012
).
65.
F.
Plasser
,
M.
Ruckenbauer
,
S.
Mai
,
M.
Oppel
,
P.
Marquetand
, and
L.
González
,
J. Chem. Theory Comput.
12
,
1207
(
2016
).
66.
G.
Granucci
,
M.
Persico
, and
A.
Toniolo
,
J. Chem. Phys.
114
,
10608
(
2001
).
67.
F.
Plasser
,
G.
Granucci
,
J.
Pittner
,
M.
Barbatti
,
M.
Persico
, and
H.
Lischka
,
J. Chem. Phys.
137
,
22A514
(
2012
).
68.
P. O.
Löwdin
,
J. Chem. Phys.
18
,
365
(
1950
).
69.
W. H.
Miller
and
S. J.
Cotton
,
Faraday Discuss.
195
,
9
(
2016
).
70.
U.
Müller
and
G.
Stock
,
J. Chem. Phys.
111
,
77
(
1999
).
71.
D.
Hu
,
Y.
Xie
,
J.
Peng
, and
Z.
Lan
,
J. Chem. Theory Comput.
17
,
3267
(
2021
).
72.
S. J.
Cotton
and
W. H.
Miller
,
J. Chem. Phys.
145
,
144108
(
2016
).
73.
L.
Lacombe
,
N. M.
Hoffmann
, and
N. T.
Maitra
,
Phys. Rev. Lett.
123
,
083201
(
2019
).
74.
P.
Martinez
,
B.
Rosenzweig
,
N. M.
Hoffmann
,
L.
Lacombe
, and
N. T.
Maitra
,
J. Chem. Phys.
154
,
014102
(
2021
).
75.
D. T.
Colbert
and
W. H.
Miller
,
J. Chem. Phys.
96
,
1982
(
1992
).
76.
M. H.
Farag
,
A.
Mandal
, and
P.
Huo
,
Phys. Chem. Chem. Phys.
23
,
16868
(
2021
).
77.
A. F.
Kockum
,
A.
Miranowicz
,
S.
De Liberato
,
S.
Savasta
, and
F.
Nori
,
Nat. Rev. Phys.
1
,
19
(
2019
).
78.
O.
Di Stefano
,
A.
Settineri
,
V.
Macrì
,
L.
Garziano
,
R.
Stassi
,
S.
Savasta
, and
F.
Nori
,
Nat. Phys.
15
,
803
(
2019
).
79.
D. G.
Baranov
,
B.
Munkhbat
,
E.
Zhukova
,
A.
Bisht
,
A.
Canales
,
B.
Rousseaux
,
G.
Johansson
,
T. J.
Antosiewicz
, and
T.
Shegai
,
Nat. Commun.
11
,
1
(
2020
).
80.
K.
Santhosh
,
O.
Bitton
,
L.
Chuntonov
, and
G.
Haran
,
Nat. Commun.
7
,
11823
(
2016
).
81.
A.
Mandal
,
S.
Montillo Vega
, and
P.
Huo
,
J. Phys. Chem. Lett.
11
,
9215
(
2020
).
82.
E. A.
Power
and
S.
Zienau
,
Philos. Trans. R. Soc. London, Ser. A
251
,
427
(
1959
).
83.
C.
Cohen-Tannoudji
,
J.
Dupont-Roc
, and
G.
Grynberg
, in
Photons and Atoms: Introduction to Quantum Electrodynamics
(
John Wiley & Sons
,
1989
).
84.
P.
Forn-Díaz
,
L.
Lamata
,
E.
Rico
,
J.
Kono
, and
E.
Solano
,
Rev. Mod. Phys.
91
,
025005
(
2019
).
85.
A.
Kelly
,
R.
van Zon
,
J.
Schofield
, and
R.
Kapral
,
J. Chem. Phys.
136
,
084101
(
2012
).
86.
M. S.
Church
,
T. J. H.
Hele
,
G. S.
Ezra
, and
N.
Ananth
,
J. Chem. Phys.
148
,
102326
(
2018
).
87.
N.
Bellonzi
,
A.
Jain
, and
J. E.
Subotnik
,
J. Chem. Phys.
144
,
154110
(
2016
).
88.
B. R.
Landry
,
M. J.
Falk
, and
J. E.
Subotnik
,
J. Chem. Phys.
139
,
211101
(
2013
).
89.
L.
Wang
,
A.
Akimov
, and
O. V.
Prezhdo
,
J. Phys. Chem. Lett.
7
,
2100
(
2016
).
90.
G.
Granucci
and
M.
Persico
,
J. Chem. Phys.
126
,
134114
(
2007
).
91.
S. J.
Cotton
and
W. H.
Miller
,
J. Chem. Phys.
150
,
104101
(
2019
).
92.
S. J.
Cotton
,
R.
Liang
, and
W. H.
Miller
,
J. Chem. Phys.
147
,
064112
(
2017
).
93.
A.
Hazra
,
A. V.
Soudackov
, and
S.
Hammes-Schiffer
,
J. Phys. Chem. B
114
,
12319
(
2010
).

Supplementary Material