On the Proper Derivation of the Floquet-based Quantum Classical Liouville Equation and Surface Hopping Describing a Molecule or Material Subject to an External Field
aa r X i v : . [ phy s i c s . c h e m - ph ] M a y On the Proper Derivation of the Floquet-based Quantum ClassicalLiouville Equation Describing a Molecule or Material Subject to anExternal Field
Hsing-Ta Chen, ∗ Zeyu Zhou, and Joseph E. Subotnik Department of Chemistry, University of Pennsylvania,Philadelphia, Pennsylvania 19104, U.S.A.
Abstract
We investigate two approaches to derive the proper Floquet-based quantum-classical Liouvilleequation (F-QCLE) for laser-driven electron-nuclear dynamics. The first approach projects the op-erator form of the standard QCLE onto the diabatic Floquet basis, then transforms to the adiabaticrepresentation. The second approach directly projects the QCLE onto the Floquet adiabatic basis.Both approaches yield a form which is similar to the usual QCLE with two modifications: 1. Theelectronic degrees of freedom are expanded to infinite dimension. 2. The nuclear motion followsFloquet quasi-energy surfaces. However, the second approach includes an additional cross deriva-tive force due to the dual dependence on time and nuclear motion of the Floquet adiabatic states.Our analysis and numerical tests indicate that this cross derivative force is a fictitious artifact,suggesting that one cannot safely exchange the order of Floquet state projection with adiabatictransformation. Our results are in accord with similar findings by Izmaylov et al ., who found thattransforming to the adiabatic representation must always be the last operation applied, thoughnow we have extended this result to a time-dependent Hamiltonian. This paper and the properderivation of the F-QCLE should lay the basis for further improvements of Floquet surface hopping. ∗ [email protected] . INTRODUCTION A computational understanding light-matter interactions for a molecular system in alaser field is useful key for interpreting spectroscopy and photochemistry, where the dy-namical interplay between electronic non-adiabatic transitions and photon excitation playsan important role for many exciting phenomena, such as molecular photodissociation[1–4]and coherent X-ray diffraction[5–8]. These phenomena usually involve dynamical processesin which electrons in a molecular system can make a transition through either (a) non-adiabatic coupling associated with the reorganization of nuclear configurations or (b) radia-tive coupling in conjunction with absorption or emission of photons. Thus, simulating theseprocesses concurrently requires accurate theoretical treatments of both non-adiabatic molec-ular dynamics and light-matter interactions[9–11]. Over the past decades, many successfulsimulation schemes have been developed based on mixed quantum−classical frameworks inwhich the electronic wavefunction evolves according to quantum mechanics while the nu-clear degrees of freedom and the laser excitation are treated as classical parameters in atime-dependent electronic Hamiltonian[12–18].Among the myraid of semiclassical dynamics, Floquet-based fewest switch surface hop-ping (F-FSSH) has emerged as one of the most powerful methods especially for simulat-ing photodissociation and ionization in a monochromatic laser field[19–21]. In a nutshell,F-FSSH integrates Floquet theory with Tully’s FSSH algorithm[22]. The general idea isto expand the electronic wavefunction in a Floquet state basis (with the electronic statesdressed by e imωt for an integer m and the laser frequency ω ), so that one can recast anexplicitly time-dependent Hamiltonian into a time-independent Floquet Hamiltonian, albeitof infinite dimension. With the Floquet Hamiltonian, one can simply employ Tully’s FSSHmethod in the Floquet state representation with a minimal modification[19, 20]. In addi-tion to the standard advantages of the usual surface hopping algorithm—stability, efficiency,and ease of incorporation with electronic structure calculations—F-FSSH also yield a bet-ter estimate for both electronic and nuclear observables than other FSSH-based methodsrelying on instantaneous adiabatic surfaces[12]. Furthermore, given the time-independentnature of the Floquet Hamiltonian, many techniques designed to improve standard FSSHmethod, such as velocity reversal and decoherence[23–30], should be applicable within theFloquet formalism. That being said, due to the fact that one cannot directly derive Tully’s2SSH—but rather only indirectly connect the equations of motion for FSSH dynamics withthe quantum-classical Liouville equation (QCLE)[31, 32]—a proper understanding of thecorrect F-QCLE is crucial if we aim to make future progress in semiclassical modeling oflight-matter interactions.Unfortunately, even without a light field, the proper derivation of the correct QCLE istricky—one can find two different versions of QCLE if one invokes slightly different formalderivations. Following Kapral’s approach[33, 34], the proper derivation of the standardQCLE includes two operations: (i) Wigner transformation and (ii) projection onto theadiabatic electronic state basis. First, one performs a partial Wigner transformation withrespect to the nuclear degrees of freedom to obtain the operator form of the QCLE. Wignertransformation provides an exact framework to interpret the full quantum density matrixin terms of the joint electronic-nuclear probability density in the phase space of the nuclearconfiguration while retaining the quantum operator character of the electronic subsystem.Second, one projects the operator form of the QCLE onto the adiabatic electronic states basisobtained by diagonalizing the electronic Hamiltonian. This adiabatic representation allowsthe connection to electronic structure calculations in a mixed quantum-classical sense[35–39]. This approach is called the Wigner-then-Adiabatic (WA) approach. As Izmaylov andco-workers have shown[40], however, exchanging these two operations (the Adiabatic-then-Wigner (AW) approach) leads to a different QCLE that cannot capture geometric phaseeffects arising from a conical intersection.With this background in mind, the proper derivation of the F-QCLE is now even morechallenging. In addition to the two operations above, there is a third step: one needs to dressthe electronic states and expand the density matrix in the Floquet state basis. In the litera-ture to date, the F-QCLE has been derived via the AW approach (projecting in the Floquetadiabatic representation, and then performing partial Wigner transformation)[21]. Never-theless, as shown by Izmaylov and co-workers[34, 40], even the limit of a time-independentHamiltonian, such an (incorrect) AW approach will lead to a QCLE that neglects geometricphase related features in the nuclear dynamics or introduces artificial nuclear effects. Despiterecent progress, the proper derivation of the F-QCLE is still an open question.In this paper, our goal is to explore different approaches to derive the F-QCLE as we shuf-fle the three key operations and quantify their differences in the context of driven electron-nuclear dynamics. By isolating the correct F-QCLE, our work will not only validate F-FSSH3ethods, it should also provide means to improve surface hopping methods. This paper isorganized as follows. In Sec. II, we formulate the three operations that are required to derivethe F-QCLE. In Sec. III, we derive F-QCLEs via approaches with different ordering of oper-ations. In Sec. IV, we implement F-FSSH calculations corresponding to these F-QCLEs andanalyze their results for a modified avoided crossing model. We conclude with an outlookfor the future in Sec. V.For notation, we denote a quantum operator by ˆ H and use bold font for matrix H . Weuse f H to denote the corresponding matrix in expanded Floquet basis (infinite dimensional).The nuclear position and momentum are ~R = { R α } , ~P = { P α } where α is the nuclearcoordinate index. We use a shorthand notation for dot product: X α · Y α = P α X α Y α . II. THREE OPERATIONS
In the context of driven electron-nuclear dynamics, let us formulate the three necessaryoperations for deriving F-QCLE. Consider a coupled electron-nuclei system driven by anexternal field of frequency ω , the total Hamiltonian takes the form of ˆ H = ˆ T ( ˆ P α ) + ˆ V ( ˆ R α , t ) where ˆ T ( ˆ P α ) = P α ( ˆ P α ) M α is the nuclear kinetic energy and ˆ V ( ˆ R α , t ) is the electronic Hamil-tonian with explicit time periodicity ˆ V ( t ) = ˆ V ( t + τ ) with τ = 2 π/ω . Formally, the dynamicsof the total system can be described by the time-dependent Schrödinger equation (TDSE) i ~ ∂∂t | Ψ i = ˆ H | Ψ i of the total wavefunction | Ψ i or the quantum Liouville equation (QLE) ∂∂t ˆ ρ = − i ~ [ ˆ H, ˆ ρ ] of the total density matrix ˆ ρ = | Ψ ih Ψ | . To derive F-QCLE, we need toapply the following three operations to the QLE. A. Partial Wigner transformation of the nuclear degrees of freedom
To describe the dynamics in a mixed quantum-classical sense, we will follow Kapral’sapproach and perform a partial Wigner transformation with respect to the nuclear degreesof freedom ˆ ρ W (cid:16) ~R, ~P , t (cid:17) = 1(2 π ~ ) N Z d ~S * ~R + ~S (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˆ ρ ( t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ~R − ~S + e − i ~P · ~S/ ~ (1)where N is the dimension of the nuclear coordinate. A nuclear position eigenstate is definedas ˆ R α | R α i = R α | R α i . In Eq. (1), the density matrix operator has been transformed intoa Wigner wavepacket in phase space with coordinates ( ~R, ~P ) . In what follows, we will4enote the partial Wigner transformed operator by the subscript W [for example ˆ V ( ˆ R α , t ) → ˆ V W ( R α , t ) ]. Note that, after the partial Wigner transformation, ˆ ρ W and ˆ V W remain electronicoperators while R α and P α are parameters.The equation of motion of the Wigner wavepacket can be obtained by transformingthe QLE by ∂∂t ˆ ρ W = − i ~ [( ˆ H ˆ ρ ) W − ( ˆ ρ ˆ H ) W ] . The Wigner transform of operator productscan be expanded further by the Wigner–Moyal operator ( ˆ H ˆ ρ ) W = ˆ H W e − i ~ ←→ Λ / ˆ ρ W with ←→ Λ = P α ←−− ∂∂P α −−→ ∂∂R α − ←−− ∂∂R α −−→ ∂∂P α [34]. Then, if we expand the the Wigner–Moyal operator inTaylor series and truncate to the first order of ~ , we obtain the operator form of the QCLE ∂∂t ˆ ρ W = − i ~ h ˆ V W , ˆ ρ W i − P α M α ∂ ˆ ρ W ∂R α + 12 ( ∂ ˆ V W ∂R α , ∂ ˆ ρ W ∂P α ) . (2)Here, the commutator is h ˆ A, ˆ B i = ˆ A ˆ B − ˆ B ˆ A and the anti-commutator is n ˆ A, ˆ B o = ˆ A ˆ B + ˆ B ˆ A .Note that Eq. (2) is exact if the partial Wigner transformed Hamiltonian is quadratic in R α ,for example harmonic oscillators.To propagate the Wigner wavepacket in Eq. (2), one must project the operator form ofthe QCLE in an electronic basis. One straightforward choice is to use a complete set ofdiabatic states for the electronic subsystem {| µ i} ; such a set does not depend on any nu-clear configuration. With this electronic diabatic basis, one can derive equations of motionfor the density matrix ( A dia νm = h ν | ˆ ρ W | µ i ) using matrix elements of the electronic Hamilto-nian, V νµ ( ~R, t ) = h µ | ˆ V W ( ~R, t ) | µ i . However, for many realistic electron-nuclei systems (andcertainly any ab initio calculations), this diabatic QCLE cannot be solved since finding acomplete set of exactly diabatic electronic states over a large set of nuclear geometries isrigorously impossible and quite demanding in practice even for approximate diabats. B. Dress the electronic basis in the Floquet formalism
Let us now focus on the Floquet formalism, according to which one solves the TDSE bytransforming the time-dependent Hamiltonian into a time-independent Floquet Hamiltonianin an extended Hilbert space of infinite dimension. For the moment let us ignore all nuclearmotion and focus on the electronic exclusively. According to Floquet theory, we utilizethe time periodicity of the electronic Hamiltonian and dress the electronic diabatic states {| µ i} by a time-periodic function e imωt where m is an integer formally from −∞ to ∞ .We denote the dressed state as the the Floquet diabatic state | µm i ≡ e imωt | µ i . In terms5f the Floquet diabatic basis, a time periodic electronic wavefunction can be expressedas | Ψ i = P µm ˜ c µm | µm i where ˜ c µm is an infinite dimensional state vector. The electronicwavefunction coefficient must satisfy the electronic TDSE i ~ X µm ∂ ˜ c µm ∂t | µm i = X µm ˆ V F ( t ) | µm i ˜ c µm (3)where the Floquet Hamiltonian operator is defined as ˆ V F ( t ) ≡ ˆ V ( t ) − i ~ ∂∂t . (4)Next, we close Eq. 3 by multiplying both sides by h ν | and write h ν | ˆ V F ( t ) | µm i = P n e V F νn,µm e inωt as a Fourier series: i ~ X m ∂ ˜ c νm ∂t e imωt = X µm X n e V F νn,µm e inωt ˜ c µm . (5)Thus, the TDSE in Eq. (3) can be solved by grouping together all terms with the same timedependence, leading to the equation of motion for ˜ ci ~ ∂∂t ˜ c νn = X µm e V F νn,µm ˜ c µm . (6)The matrix elements of the Floquet Hamiltonian can be obtained by performing a Fouriertransformation on the matrix elements e V F νn,µm = hh νn | ˆ V F | µm ii = 1 τ Z τ dt h ν | ˆ V F ( t ) | µ i e i ( m − n ) ωt . (7)Here, we define the double bracket projection by hh νn | · · · | µm ii = τ R τ dt h ν | · · · | µ i e i ( m − n ) ωt . Given that the electronic Hamiltonian operator ˆ V ( t ) is periodic in time, the double-bracketprojection eliminates all time dependence and the Floquet Hamiltonian matrix reads e V F νn,µm = hh νn | ˆ V | µm ii + δ µν δ mn m ~ ω. (8)In the end, with this time-independent Hamiltonian, Eq. (6) can be formally solved by theexponential operator exp( − i e V F t/ ~ ) with an arbitrary initial state.At this point, we will allow nuclei to move and turn out attention to the equation ofmotion for the density matrix b ρ W ( ~R, ~P , t ) within the Floquet diabatic basis. The Wigner-transformed density matrix in the Floquet diabatic representation is e A dia νn,µm ( ~R, ~P , t ) = h νn | b ρ W ( ~R, ~P , t ) | µm i . (9)6or a proper F-QCLE, we will need to calculate the time derivative of e A dia ∂∂t e A dia νn,µm = h νn | ∂ b ρ W ∂t | µm i − i ( n − m ) ~ ω e A dia νn,µm (10)where the Floquet diabatic states depends on time explicitly. We begin by using Eq. (2) toproject ∂∂t ˆ ρ W into a Floquet diabatic basis. For the commutator term in Eq. (2), we candivide the operator product into matrix multiplication: h νn | h ˆ V W , ˆ ρ W i | µm i = X λl hh νn | ˆ V W | λl ii e A dia λl,µm − e A dia νn,λl hh λl | ˆ V W | µm ii (11) = [ e V W , e A dia ] νn,µm by inserting the identity of the diabatic electronic basis: ˆ1 = P λ | λ ih λ | and expanding thetime-dependent coefficients in terms of a Fourier series; see Appendix A for more details.Furthermore, if we combine Eq. (11) with the second term on the RHS of Eq. (10), we canwrite the sum of both terms as [ e V F , e A dia ] , i.e. we can replace e V W with e V F .For the anti-commutator term, we can use the same procedure to divide the operatorproduct h νn | ( ∂ ˆ V W ∂R α , ∂ ˆ ρ W ∂P α ) | µm i = X λl hh νn | ∂ ˆ V W ∂R α | λl ii ∂ e A dia λl,µm ∂P α + ∂ e A dia νn,λl ∂P α hh λl | ∂ ˆ V W ∂R α | µm ii (12)where h νn | ∂ b ρ W ∂P α | µm i = ∂∂P α e A dia νn,µm . Note that, since the Floquet diabatic states do not dependon the nuclear coordinate, we can rewrite the derivative of the electronic Hamiltonian interms of the Floquet Hamiltonian hh νn | ∂ ˆ V W ∂R α | µm i = ∂∂R α e V F µm,νn . (13)In the end, we may combine the above expressions to write down a complete diabatic F-QCLE ∂∂t e A dia = − i ~ h e V F , e A dia i − P α M ∂ e A dia ∂R α + 12 ( ∂ e V F ∂R α , ∂ e A dia ∂P α ) (14)As a final remark, we emphasize that the Floquet Hamiltonian e V F = e V F ( ~R ) is a time-independent matrix, so Eq. (14) is simply the diabatic QCLE corresponding to an infinitelylarge electronic Hamiltonian e V F . 7 . Transformation to the adiabatic representation To recast the diabatic F-QCLE in an adiabatic representation, we diagonalize the FloquetHamiltonian matrix by solving the eigenvalue problem: X νn e V F µm,νn ( ~R ) G Jνn ( ~R ) = E F J ( ~R ) G Jµm ( ~R ) . (15)The eigenvalues E F J = E F J ( ~R ) are the so-called Floquet quasi-energies with correspondingeigenvectors G Jµm ( ~R ) . Since e V F is Hermitian, we can choose the eigenvectors G Jµm to beothornormal so that we have the identities G † G = GG † = I , i.e. P λℓ G J ∗ λℓ G Kλℓ = δ JK and P L G Lµm G L ∗ νn = δ µm,νn . The Floquet adiabatic state corresponding to the quasi-energy E F J are | φ J ( ~R, t ) i = X µm G Jµm ( ~R ) | µm i . (16)As a practical matter, although e V F is infinite dimensional, we can truncate highly oscillatingFloquet states by replacing P ∞ m = −∞ with P Mm = − M .With this Floquet adiabatic state basis, the probability density can be obtained by adiabatic-to-adiabatic transformation e A adi JK ( ~R, ~P , t ) = h Φ J | b ρ W (cid:16) ~R, ~P , t (cid:17) | Φ K i = (cid:16) G † e A dia G (cid:17) JK (17)in the Floquet adiabatic representation. Note that, since the eigenvectors G Jµm ( ~R ) do notdepend on time explicitly, the time-derivative of the adiabatic probability density can becalculated simply to be: ∂ e A adi ∂t = G † ∂ e A dia ∂t G . (18)We are now ready to derive the adiabatic F-QCLE in the following section. III. DIFFERENT APPROACHES TO DERIVE F-QCLE
In this section, we present two approaches with different orders for the three operationsabove; as will be shown, different orders will result in different adiabatic F-QCLEs. Wesummarize these approaches and the corresponding F-QCLEs in Table I. Note that, inpractice, even more approaches are possible, but we will ignore all AW approaches giventhat Izmaylov and co-workers have shown that such an ordering is inappropriate[40].8 . Wigner-Floquet-Adiabatic (WFA) approach
Our first approach follows the order presented above: we first perform partial Wignertransformation, we second project to a Floquet diabatic basis, we third transform to anadiabatic representation. For the last step, following Eq. (18), we transform the diabaticF-QCLE by sandwiching the diabatic F-QCLE (Eq. (14)) with G † and G . The first term ofEq. (14) (the commutator term) becomes X νn X µm G J ∗ νn h e V F , e A dia i νn,µm G Kµm = (cid:0) E F J − E F K (cid:1) e A adi JK (19)For the second term, since the Floquet adiabatic states depend on the nuclear coordinate ~R , the R α derivative of the density must yield X νn X µm G J ∗ νn ∂ e A dia νn,µm ∂R α G Kµm = ∂ e A adi JK ∂R α + X L (cid:16) D αJL e A adi LK − e A adi JL D αLK (cid:17) (20)where the derivative coupling is D αJK = hh Φ J | ∂∂R α | Φ K ii = P µm G J ∗ µm ∂G Kµm ∂R α corresponding tothe change of the Floquet adiabatic states with respect to the nuclear coordinate R α . Notethat, if the Floquet Hamiltonian is real, the diagonal element of the derivative coupling iszero ( D αJJ = 0 ). For the third term (the anti-commutator term), the R α derivative of theFloquet Hamiltonian can be written in terms of the force matrix X νn X µm G J ∗ νn ∂ e V F νn,µm ∂R α G Kµm = − F αJK (21)explicitly, F αJK = − ∂ E F J ∂R α δ JK + ( E F J − E F K ) D αJK . The force matrix accounts for direct changes in the nuclear momentum associated with theelectronic coupling. One can understand the diagonal element F αJJ = − ∂ E F J ∂R α as the classicalforce for nuclear dynamics moving along the J -th Floquet quasi-energy surface in the phasespace.Finally, the F-QCLE via the WFA approach reads ∂∂t e A adi JK = − i ~ (cid:0) E F J − E F K (cid:1) e A adi JK − P α M α ∂ e A adi JK ∂R α − P α M α X L (cid:16) D αJL e A adi LK − e A adi JL D αLK (cid:17) − X L F αJL ∂ e A adi LK ∂P α + ∂ e A adi JL ∂P α F αLK ! (22)9e find that Eq. 22 takes exactly the same form as the standard QCLE in the adiabaticrepresentation (for electron-nuclear dynamics without a driving laser). B. Wigner-Adiabatic-Floquet (WAF) approach
For the second approach, we exchange the “adiabatic” and “Floquet” operations after thepartial Wigner transformation. In this case, we directly project Eq. (2) onto the Floquetadiabatic basis | φ J ( ~R, t ) i . Namely, we make the diabatic-to-adiabatic transform of theFloquet electronic basis prior to the projection onto the dressed electronic states. Thus, weconsider this path as the Wigner-Adiabatic-Floquet (WAF) approach.Overall, we apply a procedure similar to what was used in Eq. 22. For the commutatorterm, we divide operator products using the same technique as in Appendix A. The R α derivative term yields a derivative coupling term as the Floquet adiabatic basis dependson R α explicitly. In the end, the WAF approach includes the first three terms exactly asEq. (22). However, from the anti-commutator term of Eq. (2), the WAF approach leads toan additional cross derivative force. To see this, we focus on the derivative of the electronicHamiltonian in the adiabatic representation hh Φ J | ∂ ˆ V∂R α | Φ K ii = hh Φ J | ∂ ˆ V F ∂R α | Φ K ii + i ~ hh Φ J | ∂∂R α ∂∂t | Φ K ii (23)where we have used the definition of the Floquet Hamiltonian ∂ ˆ V F ∂R α = ∂ ˆ V∂R α − i ~ ∂∂R α ∂∂t . Thederivative of the electronic Hamiltonian yields two terms: first, the same force matrix weobtained in Eq. (21): hh Φ J | ∂ ˆ V F ∂R α | Φ K ii = − F αJK (24)second, a cross derivative force comes from the explicit dependence of the Floquet adiabaticstates on both the nuclear coordinate and time i ~ hh Φ J | ∂∂R α ∂∂t | Φ K ii = − χ αJK , and explicitly, χ αJK = X µm m ~ ωG J ∗ µm ∂G Kµm ∂R α . (25)Note that, unlike the derivative coupling D αJK , the diagonal element of χ αJK is non-zero andreal. 10 -QCLE effective force matrixWFA Eq. (22) − ∂ E F J ∂R α δ JK + ( E F J − E F K ) D αJK WAF Eq. (26) − ∂ E F J ∂R α δ JK + ( E F J − E F K ) D αJK − P µm m ~ ωG J ∗ µm ∂G Kµm ∂R α Table I. The F-QCLEs obtained via different approaches differ in the effective force matrices fornuclear motions.
Finally, the F-QCLE via the WAF approach reads: ∂∂t e A adi = − i ~ (cid:0) E F J − E F K (cid:1) e A adi JK − P α M α ∂ e A adi JK ∂R α − P α M α X L (cid:16) D αJL e A adi LK − e A adi JL D αLK (cid:17) − X L ( F αJL + χ α ∗ LJ ) ∂ e A adi LK ∂P α + ∂ e A adi JL ∂P α ( F αLK + χ αLK ) (26)We observe that, while Eq. (22) and Eq. 26 take the same form, the “effective” force matrix(defined as the coefficients of ∂ e A adi ∂P α ) includes an additional cross derivative force that indicatesthe difference between these two equations. The similarities suggest that the time evolutionof the electronic degrees of freedom ( ∂ e A adi ∂R α = ∂ e A adi ∂P α = 0 ) should follow the same equation forboth the WFA and WAF. However, the difference in the “effective” force matrix will affectthe nuclear dynamics in phase space. Specifically, from Eq. 26, the classical force on the J -th quasi-energy surface is given by F αJJ + χ αJJ , implying that the nuclear dynamics shouldexperience the additional cross derivative force on top of the Floquet quasi-energy surface. IV. RESULTSA. Shifted avoided crossing model
To analyze the difference between these approaches, we consider a shifted avoided crossingmodel composed of a two-level electronic system coupled to a 1D nuclear motion. Thediabatic electronic states are denoted as | g i and | e i and the electronic Hamiltonian is givenby ˆ V ( R, t ) = V gg ( R ) V ge ( R, t ) V eg ( R, t ) V ee ( R ) (27)11here the diabatic energy is given by V gg ( R ) = A (1 − e − BR ) R > − A (1 − e BR ) R < (28) V ee ( R ) = − V gg ( R ) + ~ ω (29)and the diabatic coupling is periodic in time V ge ( R, t ) = V eg ( R, t ) = Ce − DR cos ωt. (30)The parameters are A = 0 . , B = 1 . , C = 0 . , D = 1 . , and the nuclear mass is M =2000 . Note that this model is Tully’s simple avoided crossing model with two modifications:the diabatic coupling becomes time-periodic and the excited potential energy surface isshifted by ~ ω . For simplicity, we choose the laser frequency ~ ω = 0 . large enough( ~ ω > A ) so that the diabatic Floquet states | gm i have an avoided crossing only with | e ( m − i (at R = 0 ) and do not have any trivial crossings.We assume the initial wavepacket is a Gaussian centered at the initial position R andmomentum P on diabat | g i : | Ψ i = 1 N exp (cid:18) − ( R − R ) σ + i ~ P ( R − R ) (cid:19) | g i (31)where the normalization factor is N = πσ . The width of the Gaussian is chosen to be σ = 20 ~ /P . The wavepacket can be propagated exactly in the diabatic representation. B. F-FSSH based on WFA and WAF
To show the difference between the dynamics as obtained by the different F-QCLEs, wewill simulate F-FSSH results for both the WFA and WFA approaches. Within F-FSSH, wedescribe the propagation of the Floquet wavepacket by a swarm of trajectories, each withits own electronic amplitudes ˜ c J satisfying ∂ ˜ c J ∂t = − i ~ E F J ˜ c J − PM X L D JL ˜ c L (32)where E F J = E F J ( R ) , D JL = D JL ( R ) and R = R ( t ) , P = P ( t ) representing nuclear trajectory.All nuclear trajectories move classically along an active Floquet state ( J ) obeying dRdt = PM (33)12 Pdt = − ∂ E F J ∂R for WFA − ∂ E F J ∂R − χ JJ for WAF (34)Here, based on the connection between the QCLE and FSSH, the nuclear force in Eq. 34 isdetermined according to the diagonal element of the effective force matrix in the F-QCLEs.Consistent with the standard FSSH technique, the hopping probability from active Flo-quet state J to state K is given byProb ( J → K ) = − Re (cid:18) P α M α · D αKJ ˜ c J ˜ c ∗ K | ˜ c J | (cid:19) dt where dt is the classical time step. After each successful hop, the velocity is adjusted toconserve the total Floquet quasi-energy. If a frustrated hop occurs, we implement veloc-ity reversal[25]. Note that we neglect the decoherence correction since the over-coherenceproblem should not be severe for a simple avoided crossing model[41].In the end, the probability to measure diabatic state µ can be evaluated by the densitymatrix interpretation[42] P µ = X m N ( µm ) N traj + X n = m P N ( µm ) l P N ( µn ) k ˜ c ( l ) µm ˜ c ( k ) ∗ µn e i ( m − n ) ωt N ( µm ) N ( µn ) (35)where N ( µm ) = P N traj l δ J ( l ) µm is the number of the trajectories that have the active surface J ( l ) end up on the Floquet state | µm i . Here l and k are the trajectory indices. We propagate N traj trajectories for an amount of time long enough for each trajectory to pass through thecoupling region ( | R | < for this parameter set). C. Effective Floquet quasi-energy surfaces for nuclear dynamics
First, we analyze the effective potential energy surfaces for nuclear dynamics by integrat-ing Eq. 34 over R for the WFA and WAF approaches respectively. For the WFA approach,the effective PES simply recovers the Floquet quasi-energy surfaces E F J . For the WAF ap-proach, the effective PES is V eff ( R ) = E F J ( R ) + R R −∞ χ JJ ( R ′ ) dR ′ where the quasi-energysurface is modified by the integration of the cross derivative force. We find that includingthe cross derivative force in the WAF approach increases the crossing barrier for the nucleardynamics on the lower adiabat (see Fig. 1). Note that, in terms of an F-FSSH calcula-tion, these changes will have a direct effect on the nuclear dynamics, but not the electronicamplitudes. 13 igure 1. The effective potential energy surfaces for the WFA approach (solid lines) and the WAFapproach (dash lines). The lower (upper) quasi-energy surface is in red (blue). The diabatic Floquetstate energies e V F µm,µm are plotted for | g i (green) and | e i (orange) in dotted lines. Note that thebarrier height and the equilibrium energy ratio of the WAF surface is significantly modified (relativeto the WFA approach) by the presence of the additional cross derivative force. D. Transmission and reflection
Next, we turn our attention to the transmission and reflection probabilities producedby the F-FSSH calculations. Overall, the WFA results are more accurate than the WAFresults. In Fig. 2, we find that, according to the WFA approach, there should be a risein transmission on the lower adiabat around P ≈ . , which is the momentum for whichtransmission should be allowed classically; see the barrier height ( ≈ . ) in Fig. 1. Indeed,such a threshold at P ≈ . is found according to exact wavefunction simulation as well.However, for the WAF approach with the cross derivative force, one find a higher crossingbarrier energy ( ≈ . ), and the transmission on the lower adiabat occurs (incorrectly) at P ≈ . . This result suggests that the cross derivative force is a fictitious term: the WAFsemiclassical derivation is spurious.Let us now focus on the WFA results in more detail. Several points are worth mentioning.First, the transmission to the upper adiabat occurs after P ≈ . , which agrees with theclassical energy difference A = 0 . (see Eq. (28) for the definition of A ). Second, forhigh initial momentum ( P > . ), the F-FSSH-WFA can almost recover the correct nucleardynamics. Third, in the intermediate momentum region P ∈ (6 , , unfortunately, the F-14 Initial Momentum P Transmission
Initial Momentum P P r obab ili t y P r obab ili t y Reflection
ExactF-FSSH (WFA)F-FSSH (WAF)
UpperLower LowerUpper
Figure 2. The probability of transmission (right) and reflection (left) on the upper and loweradiabats as a function of the initial momentum. The F-FSSH dynamics is implemented usingthe effective nuclear forces of the WFA (red) and WAF (blue). Overall, the WFA result is moreaccurate than the WAF result. The WFA result almost recover the correct nuclear dynamics. Dueto the cross derivative force, the nuclear trajectory of the F-FSSH(WAF) experiences a much highercrossing barrier, requiring larger P for transmission. FSSH wavepacket exhibits less transmission than the exact calculation. This discrepancymay be attributed to FSSH’s incapability to capture nuclear tunneling effects.
V. CONCLUSION AND OUTLOOK
We have analyzed two approaches for deriving the QCLE within a Floquet formalism,and found two different F-QCLEs. While these F-QCLEs take similar forms, the differencein the effective force matrix can lead to large discrepancies in nuclear dynamics. As such,in the context of driven electron-nuclear dynamics, our results reiterate the fact that onecannot change the order of the operations in the derivation of the correct QCLE. Specifically,as opposed to the WFA approach, the WAF approach [exchanging Floquet electronic basisdressing (F) and adiabatic transformation (A)] is spurious. Overall, our results are veryconsistent with the results of Izmaylov and Kapral[40] who find that one must be carefulwhen deriving the QCLE even without time dependence; in the end, with or without a15ime-dependent Hamiltonian, it appears that one will always derive the correct semiclassicalequation of motion provided that one moves to the adiabatic representation as the very laststep .Looking forward, the derivation of the F-QCLE presented here validates the F-FSSHmethod and paves the way to further improvements in the future. With regard to coherenceand decoherence, given the time-independent nature of the Floquet Hamiltonian, we can im-mediately apply many decoherence schemes, including augmented moment decoherence[28,41, 43–45], to the F-FSSH algorithm. As far as geometric phase effect is concerned, it isknown that Berry phases are already included within the QCLE[46] for a time-independentelectronic Hamiltonian, and so we would expect that similar effects should already be in-cluded within this proper F-QCLE for periodic (time-dependent) electronic Hamiltonians.Nevertheless, however, there is one nuance which we have conveniently neglected in thepresent paper. Note that, according to Eq. (7), we have every reason to believe that theF-QCLE formalism (especially for a non-monocrhomatic driving field with more than oneFourier mode in the time-dependent Hamiltonian) will necessarily introduce a complex (i.e.not real) Floquet Hamiltonian. In such a case, we should find not just Berry phases, but alsoBerry force[47]. Future research into the nature of this intrinsic magnetic Berry force—howor if it appears in the context of F-QCLE and surface hopping dynamics—is currently un-derway and represents an exciting new direction for non-adiabatic theory.Lastly, it has been recently reported that novel control schemes, such as Floquetengineering[48, 49], can enhance the excitation energy transfer rate even in the presenceof strong fluctuations and dissipation. Given so many potential applications for F-QCLEsimulations, we believe the present manuscript should find immediate use in the physicalchemistry and chemical physics community.
ACKNOWLEDGMENT
This material is based upon work supported by the U.S. Department of Energy, Office ofScience, Office of Basic Energy Sciences under Award Number de-sc0019397. This researchalso used resources of the National Energy Research Scientific Computing Center (NERSC),a U.S. Department of Energy Office of Science User Facility operated under Contract No.DE-AC02-05CH11231. We thank Abraham Nitzan for very helpful discussions.16
ATA AVAILABLITY
The data that support the findings of this study are available from the correspondingauthor upon reasonable request.
Appendix A: Rewriting the Product of Wignerized Operators As Matrix Multipli-cation in the Floquet Representation
Throughout this paper, we have constantly used one trick. Namely, we have consistentlyrewritten the product of two Wigner transformed operators (one of which must be time-periodic) into a non-standard matrix product in the Floquet representation. To see how thistrick works in practice, we consider (for example) the operator product ˆ ρ W ˆ H W .We insert the electronic identity operator ˆ1 = P λ | λ ih λ | in between ˆ ρ W and ˆ H W : h νn | ˆ ρ W ˆ H W | µm i = X λ h νn | ˆ ρ W | λ ih λ | ˆ H W | µm i (A1)Next, as in Sec. II, we express h λ | ˆ H W | µm i in the form of Fourier series: h λ | ˆ H W | µm i = X l hh λl | ˆ H W | µm ii e ilωt (A2)where the double bracket projection is defined by Fourier transform. With the Fourier series,we write | λ i e ilωt = | λl i and obtain h νn | ˆ ρ W ˆ H W | µm i = X λl h νn | ˆ ρ W | λl ihh λl | ˆ H W | µm ii (A3) [1] P. M. Regan, S. R. Langford, A. J. Orr-Ewing, and M. N. R. Ashfold,The Journal of Chemical Physics , 281 (1999).[2] K. J. Franks, H. Li, and W. Kong, The Journal of Chemical Physics , 11779 (1999).[3] K. I. Hilsabeck, J. L. Meiser, M. Sneha, N. Balakrishnan, and R. N. Zare,Physical Chemistry Chemical Physics (2019), 10.1039/C8CP06107F.[4] K. I. Hilsabeck, J. L. Meiser, M. Sneha, N. Balakrishnan, and R. N. Zare,Phys. Chem. Chem. Phys. , 14195 (2019), publisher: The Royal Society of Chemistry.
5] J. M. Glownia, A. Natan, J. P. Cryan, R. Hartsock, M. Kozina, M. P. Minitti, S. Nelson,J. Robinson, T. Sato, T. van Driel, G. Welch, C. Weninger, D. Zhu, and P. H. Bucksbaum,Physical Review Letters (2016), 10.1103/PhysRevLett.117.153003.[6] J. M. Glownia, A. Natan, J. P. Cryan, R. Hartsock, M. Kozina, M. P. Minitti, S. Nelson,J. Robinson, T. Sato, T. van Driel, G. Welch, C. Weninger, D. Zhu, and P. H. Bucksbaum,Phys. Rev. Lett. , 069302 (2017).[7] H. T. Lemke, K. S. Kjær, R. Hartsock, T. B. van Driel, M. Chollet, J. M. Glownia, S. Song,D. Zhu, E. Pace, S. F. Matar, M. M. Nielsen, M. Benfatto, K. J. Gaffney, E. Collet, andM. Cammarata, Nature Communications , 15342 (2017).[8] F. D. Fuller, S. Gul, R. Chatterjee, E. S. Burgie, I. D. Young, H. Lebrette, V. Srinivas, A. S.Brewster, T. Michels-Clark, J. A. Clinger, B. Andi, M. Ibrahim, E. Pastor, C. d. Lichtenberg,R. Hussein, C. J. Pollock, M. Zhang, C. A. Stan, T. Kroll, T. Fransson, C. Weninger, M. Kubin,P. Aller, L. Lassalle, P. Bräuer, M. D. Miller, M. Amin, S. Koroidov, C. G. Roessler, M. Allaire,R. G. Sierra, P. T. Docker, J. M. Glownia, S. Nelson, J. E. Koglin, D. Zhu, M. Chollet, S. Song,H. Lemke, M. Liang, D. Sokaras, R. Alonso-Mori, A. Zouni, J. Messinger, U. Bergmann, A. K.Boal, J. M. Bollinger, C. Krebs, M. Högbom, G. N. Phillips, R. D. Vierstra, N. K. Sauter,A. M. Orville, J. Kern, V. K. Yachandra, and J. Yano, Nat Methods , 443 (2017), number:4 Publisher: Nature Publishing Group.[9] J. J. Bajo, G. Granucci, and M. Persico, The Journal of Chemical Physics , 044113 (2014).[10] N. M. Hoffmann, H. Appel, A. Rubio, and N. T. Maitra, Eur. Phys. J. B , 180 (2018).[11] A. Abedi, N. T. Maitra, and E. K. U. Gross, Physical Review Letters (2010), 10.1103/PhysRevLett.105.123002.[12] Z. Zhou, H.-T. Chen, A. Nitzan, and J. E. Subotnik,J. Chem. Theory Comput. , 821 (2020).[13] M. Richter, P. Marquetand, J. González-Vázquez, I. Sola, and L. González,Journal of Chemical Theory and Computation , 1253 (2011).[14] R. Mitrić, J. Petersen, and V. Bonačić-Koutecký, Physical Review A (2009), 10.1103/PhysRevA.79.053416.[15] S. Mai, P. Marquetand, and L. González, Wiley Interdisciplinary Reviews: Computational Molecular Science , e1370 (2018).[16] M. Thachuk, M. Y. Ivanov, and D. M. Wardlaw, The Journal of Chemical Physics , 4094 (1996).[17] J. J. Bajo, J. González-Vázquez, I. R. Sola, J. Santamaria, M. Richter, P. Marquetand, andL. González, The Journal of Physical Chemistry A , 2800 (2012).[18] P. G. Lisinetskaya and R. Mitrić, Physical Review A (2011), 10.1103/PhysRevA.83.033408.
19] T. Fiedlschuster, J. Handt, and R. Schmidt, Physical Review A , 053409 (2016).[20] T. Fiedlschuster, J. Handt, E. K. U. Gross, and R. Schmidt,Physical Review A (2017), 10.1103/PhysRevA.95.063424.[21] I. Horenko, B. Schmidt, and C. Schütte, The Journal of Chemical Physics , 5733 (2001).[22] J. C. Tully, J. Chem. Phys. , 22A301 (2012).[23] E. R. Bittner and P. J. Rossky, J. Chem. Phys. , 8130 (1995).[24] A. W. Jasper and D. G. Truhlar, J. Chem. Phys. , 64103 (2005).[25] A. W. Jasper and D. G. Truhlar, Chem. Phys. Lett. , 60 (2003).[26] T. Nelson, S. Fernandez-Alberti, A. E. Roitberg, and S. Tretiak,J. Chem. Phys. , 224111 (2013).[27] J. E. Subotnik, A. Jain, B. Landry, A. Petit, W. Ouyang, and N. Bellonzi,Annual Review of Physical Chemistry , 387 (2016).[28] A. Jain, E. Alguire, and J. E. Subotnik, Journal of Chemical Theory and Computation , 5256 (2016).[29] J. Y. Fang and S. Hammes-Schiffer, J. Phys. Chem. A , 9399 (1999).[30] M. J. Bedard-Hearn, R. E. Larsen, and B. J. Schwartz, J. Chem. Phys. , 234106 (2005).[31] J. E. Subotnik, W. Ouyang, and B. R. Landry, J. Chem. Phys. , 214107 (2013).[32] R. Kapral, Chemical Physics , 77 (2016).[33] R. Kapral, Annu. Rev. Phys. Chem. (2006), 10.1146/annurev.physchem.57.032905.104702.[34] R. Kapral and G. Ciccotti, J. Chem. Phys. , 8919 (1999), publisher: American Institute ofPhysics.[35] K. F. Wong and P. J. Rossky, J. Chem. Phys. , 8429 (2002).[36] K. F. Wong and P. J. Rossky, J. Chem. Phys. , 8418 (2002).[37] F. Webster, P. Rossky, and R. Friesner, Comput. Phys. Commun. , 494 (1991).[38] A. Kelly and T. E. Markland, J. Chem. Phys. , 014104 (2013).[39] H. W. Kim, A. Kelly, J. W. Park, and Y. M. Rhee, J. Am. Chem. Soc. , 11640 (2012),publisher: American Chemical Society.[40] I. G. Ryabinkin, C.-Y. Hsieh, R. Kapral, and A. F. Izmaylov,The Journal of Chemical Physics , 084104 (2014).[41] B. R. Landry and J. E. Subotnik, J. Chem. Phys. , 0 (2012).[42] B. R. Landry, M. J. Falk, and J. E. Subotnik, The Journal of Chemical Physics , 211101 (2013).[43] A. S. Petit and J. E. Subotnik, The Journal of Chemical Physics , 014107 (2014).
44] B. J. Schwartz, E. R. Bittner, O. V. Prezhdo, and P. J. Rossky,J. Chem. Phys. , 5942 (1996).[45] O. Prezhdo and P. Rossky, Phys. Rev. Lett. , 5294 (1998).[46] J. Subotnik, G. Miao, N. Bellonzi, H.-H. Teh, and W. Dou,J. Chem. Phys. , 074113 (2019).[47] G. Miao, N. Bellonzi, and J. Subotnik, J. Chem. Phys. , 124101 (2019).[48] N. Thanh Phuc and A. Ishizaki, The Journal of Physical Chemistry Letters , 1243 (2018).[49] K. Schwennicke and J. Yuen-Zhou, The Journal of Physical Chemistry C (2020), 10.1021/acs.jpcc.9b10030., 1243 (2018).[49] K. Schwennicke and J. Yuen-Zhou, The Journal of Physical Chemistry C (2020), 10.1021/acs.jpcc.9b10030.