Coupling molecules to the quantized radiation field inside an optical cavity creates a set of new photon-matter hybrid states, which are commonly referred to as polaritons. [1][2][3][4][5][6] 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][7][8][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][11][12][13] as these polariton chemical reactions often involve a rich dynamical interplay among the electronic, nuclear, and photonic degrees of freedom (DOFs). Accurately simulating polaritonic quantum dynamics remains a challenging task and is beyond the scope of photochemistry 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][15][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][11][12][13][19][20][21][22] The Journal of Chemical Physics
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][20][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 ultrastrong coupling regimes. [30][31][32][33][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.
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 Tn 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 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, Ĥen, is expressed as follows:
which describes the molecular Hamiltonian (without the nuclear kinetic energy). The above expression in Eq. ( 2) includes electronic kinetic energy Te, electron-electron interaction Vee, electronic-nuclear interaction Ven, and nuclei-nuclei interaction Vnn. The expressions of these four terms can be found in standard textbooks. [37][38][39] Modern electronic structure theory is focused on solving the eigenvalue problem of Ĥen, 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 â † and â are the photonic creation and annihilation operator, respectively, qc = h2ωc(â † + â) and pc = -i hωc2(ââ † ), 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
The Journal of Chemical Physics
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 λc = 1 0 V 0 , 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 = hωc2λ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 ge (R) = ge 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 σ † ≡ eg and σ ≡ ge, and thus = eg (R) ⋅ (σ † + σ). The molecule-cavity interaction term in Eq. ( 5) can now be expressed as
Assuming the RWA by ignoring the counter-rotating terms (CRTs) â † σ † and âσ, and explicitly dropping the DSE term Ĥd 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][46][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
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 Tn represents the nuclear kinetic energy operator, and V represents the rest of the Hamiltonian. For a bare molecular system, V = Ĥen 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 Ĥpl . For the molecule-cavity hybrid system, we define the polaritonic state 3,50 as the eigenstate of V = Ĥpl [see definition in Eq. ( 12)] through the following eigenequation:
where EI(R) 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 ci is the expansion coefficients. The wave function satisfies the time-dependent Schödinger equation (TDSE),
Plugging Eq. ( 15) into (16), we obtain the equations of motion for ci as follows:
where Vij = ψi(R) Vψj(R). The time derivative coupling (or nonadiabatic coupling, NAC) between two basis states is
where d α ij is the NAC associated with the α th atom, defined as follows:
The Journal of Chemical Physics
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. dEdt = 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:
where the derivative coupling is defined as Eq. ( 19). Inserting Eq. ( 23) into Eq. ( 22), we have
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 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 ∇α V using the chain rule as follows:
which is indeed equivalent to Eq. ( 29).
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 [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 2 γν denotes the matrix elements of DSE, and the sum ∑ ξ in the matrix element of Ĥd 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 inside the DSE term within the matter subspace, with the projection operator P = ∑ ξ ϕ ξ ϕ ξ . 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 Î = P + Q, where Q contains the adiabatic electronic states outside the subspace defined by P, P P is properly confined in the subspace P, whereas P 2 P = P( P + Q) P contains the terms outside the subspace P. 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 [Eq. ( 6)] under the adiabatic representation
which parametrically depends on R.
To express the nuclear gradients in Eq. ( 29), we need the gradients ∇αVij 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 ∇αVij, we need to take the derivative on each term of Vij(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 α γν 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 [∇αV]ij for a molecule-cavity hybrid system is presented in Eq. ( 29), with the matrix elements of Vij 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 Eg and Ee, 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 {e0, g1} and {g0, e0, g1, e1}, which are referenced hereafter as two-state and four-state basis, respectively.
Under the two-state subspace {e0, g1}, the potential matrix for the molecule-cavity hybrid system under this two-state basis is
The derivative coupling between g1 and e0 states is
due to the fact that 10 = 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 g1∇αg1 = 0, g1∇αe0 = 0, and e0∇αe0 = 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 {g0, e0, g1, e1}, the potential matrix is
The above matrix is Hermitian because eg = ge and D 2 eg = D 2 ge . The derivative coupling matrix
where d α ge = g∇αe = -d α eg . 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 Ĥen [Eq. ( 2)]. Furthermore, one can clearly see that the DSE term is a necessary component in ĤPF [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 [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][52][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.
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 [defined in Eq. ( 12)]. Some MQC methods, such as FSSH, 17 are formulated specifically in the eigenstates of V and require the calculation of the nuclear gradients on the eigenstates of V. This means that one needs to obtain the polariton eigenstates {EI(R)} by solving Eq. ( 13). The transformation between {ψ i (R)} = {ϕ ν (R) ⊗ n} and {EI(R)} 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). 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) polaritonic state is
where we have used the Hellman-Feynman theorem. Assuming the completeness relation ∑ i ψiψi = Î (where
and inserting it into Eq. ( 42), we have
where the transformation matrix element is U kI = ψ k EI. 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 d α ij defined in Eq. (33). One can further express Eq. ( 44) by inserting the completeness relation as
where U lJ = ψ l EJ and ψ k ∇αVψ l can be evaluated using Eq. (30).
To give a more concrete example, for the two-state basis {e0, g1}, we have
where ∇αE 1 and ∇αE 2 are gradients of two adiabatic (polaritonic) states and d α 12 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:
The Journal of Chemical Physics
where the mixing angle θ satisfies the condition 54
with ge = gc ⋅ ge and V = Ve + D 2 ee -Vg -D 2 gg -hω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 d α 12 = E 1 ∇αE 2 . Using Eq. ( 49), one can find that
Using the expression of the mixing angle θ in Eq. ( 48), we have
The above d α 12 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 V. This is not surprising because, within the {e0, g1} subspace, the light-matter interaction Hamiltonian reduces to a JC type Hamiltonian [see Eq. ( 34)].
In this work, we use the Shin-Metiu (SM) model 55 as the "ab initio" model molecular system to investigate strong and ultrastrong 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 Rc = 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 Ĥel in this grid basis {ri} are given by ri Ĥel rj = ri Tr + VeN (r, R) + VNN (R)rj
where the ri Trrj is given analytically 56 as follows:
Directly diagonalizing the matrix of Ĥel [in Eq. ( 53)] at a given nuclear (proton) position R in this grid basis gives the accurate adiabatic electronic states
where c α i (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. mode is chosen as hωc = 2.721 eV (≈ 0.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 , 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 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 hω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 moleculecavity 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
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.002 70 a.u.; for SM2, we use 40 R 0 = -4 and ω = 0.000 382 a.u.
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 ϕ λ ∇ Ĥenϕν as follows:
where the ∇ Ĥen(R) is evaluated analytically using the expression of Ven(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 ∇c γ i (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, the so-called active state as follows:
where EI is the active adiabatic polariton state (which is the eigenstate of V). 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 ρIJ = c * I cJ is the electronic density matrix element. Since the probability should not be negative, we set fIJ to 0 if fIJ < 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 dIK(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
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). When computing the population in a representation that is rather than the adiabatic states of V, there is no unique way 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 [Rχ]w = 1 hπ e -M(P 2 +ω 2 (R-R ) 2 )ω h , 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 e0. 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). We follow the previous work 62 and use the canonical 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 .
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) ĤΨ = εΨ, 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 g0, e0, g1, and e1 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 -h 2 2M.
Since we already derived the formula of V 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 cj 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 dge as g∇ 2 g = e∇ 2 e = -dge 2 , (73a) g∇ 2 e = -e∇ 2 g = ∇⋅dge.
(
The above strategy to solve the exact polariton dynamics can also be found in many other previous studies, such as Refs. 63 and 64.
For more than two electronic states, these quantities are evaluated numerically.
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 â † â 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 â † â. Nevertheless, it can be viewed as an approximate estimation of the photon number when the light-matter couplings are not in the ultrastrong 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 hωc = 1.281 eV (≈0.047 a.u.).
Figure 2(b) presents the exact quantum dynamics within the JC subspace {e0, g1} (dashed lines) or in the larger {g0, e0, g1, e1} (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 e1 state is mainly contributed from the population transfer from the e0 state due to the light-matter coupling originated from the permanent dipole ee . In addition, the population transition between e0 and g1 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 Xij matrix in the general gradient expression [Eq. (40)] in the subspace of {g0, e0, g1, e1}. 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 ∇[V] e0,g1 = gc ⋅∇ ge and the regular gradient from the electronic derivative couplings X g0,e0 = dge(Ve -Vg + D 2 ee -D 2 gg ) are compared with the rest of the non-JC type of gradients, such as X g0,g0 = 2dgeD 2 eg (where The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp the matrix of the DSE operator), X g0,e1 = dgegc ⋅ eegg , and X g0,g1 = 2dge g c ⋅ ge . As we show in Fig. 2(c), all of the gradients we discussed above have a similar magnitude in this Shin-Metiucavity 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 [ρij = c * i cj, where ci 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 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 e1 starts to mix with the states within the JC subspace. The g0 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 e0 and g1 states, and the e1 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 e0 and e1 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 hωc = 0.1 a.u., causing a light-induced avoid crossing near R = -4.0 a.u. [see Fig. 4 When the coupling is g c = 0.001, the polaritonic PESs [Fig. 4(a)] have nearly identical curvatures of the photon dressed states g1 and e0 (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 g0 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 {g1, e0} under the same coupling strength. Although SM2 has a large permanent dipole, the population of g0 only starts to grow later in time (20 fs) after the initial excitation. Thus, the population transfer to g0 and e1 states are mainly due to the electronic NAC deg that directly couples the g1 state to e1 state, as well as the e0 state to g0 state. Note that the e1 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 e1 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 e1 starts to grow at an earlier time, due to the permanent dipole ee that couples e0 state to the e1 state, whereas the electronic NAC deg still contribute to the later time population transfer to the g0 state. Through out all coupling strengths investigated here, both Ehrenfest and FSSH approaches provide a reasonable accuracy for the population dynamics for the g1 and e0 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 photondressed ground state g1 state, which can be resonant with the e0 state and then hybridize with it. The transition dipole eg causes the superposition between e0 and g1 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 {e0, g1} 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 deg [Eq. (19)] between e0 and g0 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 deg 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 e0 and e1 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 e1 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.
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 d α ij = -d α ji .
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., dij = 0, and Eq. ( 24) becomes the Ehrenfest equation of motion
For the adiabatic basis, the electronic potential matrix [V] now is diagonal, i.e., Vij = Ei(R)δij, where Ei is i th eigenvalue satisfying Vϕi = Eiϕi. Under this case, Eq. ( 24) becomes
where we have used the chain rule ϕi∇α Vϕj = ∇αϕi Vϕj -∇αϕi Vϕjϕi V∇αϕj as well as the following well-known equality: 18 ϕi∇α Vϕj = ∇αEi, (i = j),
The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp APPENDIX C: ADDITIONAL NUMERICAL RESULTS
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.
We demonstrate the importance of the DSE terms in polariton quantum dynamics simulations, especially in the strong and ultrastrong light-matter coupling regimes. ge 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 37) are also presented for comparison, such as g c ⋅ gg (red), g c ⋅ ge (green), and g c ⋅ ee (blue). 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.
Published under an exclusive license by AIP Publishing