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

To describe the dynamics in a mixed quantum-classical sense, we will follow Kapral’s approach and perform a partial Wigner transformation with respect to the nuclear degrees of freedom

$$\hat{\rho}_{W}\left(\vec{R},\vec{P},t\right)=\frac{1}{(2\pi\hbar)^{N}}\int d\vec{S}\left\langle\vec{R}+\frac{\vec{S}}{2}\right|\hat{\rho}(t)\left|\vec{R}-\frac{\vec{S}}{2}\right\rangle e^{-i\vec{P}\cdot\vec{S}/\hbar}$$

(1)

where $N$ is the dimension of the nuclear coordinate. A nuclear position eigenstate is defined as $\hat{R}^{\alpha}|R^{\alpha}\rangle=R^{\alpha}|R^{\alpha}\rangle$. In Eq. (1), the density matrix operator has been transformed into a Wigner wavepacket in phase space with coordinates $(\vec{R},\vec{P})$. In what follows, we will denote the partial Wigner transformed operator by the subscript $W$ [for example $\hat{V}(\hat{R}^{\alpha},t)\rightarrow\hat{V}_{W}(R^{\alpha},t)$]. Note that, after the partial Wigner transformation, $\hat{\rho}_{W}$ and $\hat{V}_{W}$ remain electronic operators while $R^{\alpha}$ and $P^{\alpha}$ are parameters.

The equation of motion of the Wigner wavepacket can be obtained by transforming the QLE by $\frac{\partial}{\partial t}\hat{\rho}_{W}=-\frac{i}{\hbar}[(\hat{H}\hat{\rho})_{W}-(\hat{\rho}\hat{H})_{W}]$. The Wigner transform of operator products can be expanded further by the Wigner–Moyal operator $(\hat{H}\hat{\rho})_{W}=\hat{H}_{W}e^{-i\hbar\overleftrightarrow{\Lambda}/2}\hat{\rho}_{W}$ with $\overleftrightarrow{\Lambda}=\sum_{\alpha}\overleftarrow{\frac{\partial}{\partial P^{\alpha}}}\overrightarrow{\frac{\partial}{\partial R^{\alpha}}}-\overleftarrow{\frac{\partial}{\partial R^{\alpha}}}\overrightarrow{\frac{\partial}{\partial P^{\alpha}}}$(Kapral and Ciccotti, 1999). Then, if we expand the the Wigner–Moyal operator in Taylor series and truncate to the first order of $\hbar$, we obtain the operator form of the QCLE

$$\frac{\partial}{\partial t}\hat{\rho}_{W}=-\frac{i}{\hbar}\left[\hat{V}_{W},\hat{\rho}_{W}\right]-\frac{P^{\alpha}}{M^{\alpha}}\frac{\partial\hat{\rho}_{W}}{\partial R^{\alpha}}+\frac{1}{2}\left\{\frac{\partial\hat{V}_{W}}{\partial R^{\alpha}},\frac{\partial\hat{\rho}_{W}}{\partial P^{\alpha}}\right\}.$$

(2)

Here, the commutator is $\left[\hat{A},\hat{B}\right]=\hat{A}\hat{B}-\hat{B}\hat{A}$ and the anti-commutator is $\left\{\hat{A},\hat{B}\right\}=\hat{A}\hat{B}+\hat{B}\hat{A}$. Note that Eq. (2) is exact if the partial Wigner transformed Hamiltonian is quadratic in $R^{\alpha}$, for example harmonic oscillators.

To propagate the Wigner wavepacket in Eq. (2), one must project the operator form of the QCLE in an electronic basis. One straightforward choice is to use a complete set of diabatic states for the electronic subsystem $\{|\mu\rangle\}$; such a set does not depend on any nuclear configuration. With this electronic diabatic basis, one can derive equations of motion for the density matrix ($A_{\nu\mu}^{\text{dia}}=\langle\nu|\hat{\rho}_{W}|\mu\rangle$) using matrix elements of the electronic Hamiltonian, $V_{\nu\mu}(\vec{R},t)=\langle\nu|\hat{V}_{W}(\vec{R},t)|\mu\rangle$. However, for many realistic electron-nuclei systems (and certainly any ab initio calculations), this diabatic QCLE cannot be solved since finding a complete set of exactly diabatic electronic states over a large set of nuclear geometries is rigorously impossible and quite demanding in practice even for approximate diabats.

Let us now focus on the Floquet formalism, according to which one solves the TDSE by transforming the time-dependent Hamiltonian into a time-independent Floquet Hamiltonian in an extended Hilbert space of infinite dimension. For the moment let us ignore all nuclear motion and focus on the electronic subsystem exclusively. According to Floquet theory, we utilize the time periodicity of the electronic Hamiltonian and dress the electronic diabatic states $\{|\mu\rangle\}$ by a time-periodic function $e^{im\omega t}$ where $m$ is an integer formally from $-\infty$ to $\infty$. We denote the dressed state as the the Floquet diabatic state $|\mu m\rangle\equiv e^{im\omega t}|\mu\rangle$. In terms of the Floquet diabatic basis, a time periodic electronic wavefunction can be expressed as $|\Psi\rangle=\sum_{\mu m}\tilde{c}_{\mu m}|\mu m\rangle$ where $\tilde{c}_{\mu m}$ is an infinite dimensional state vector. The electronic wavefunction coefficient must satisfy the electronic TDSE

$$i\hbar\sum_{\mu m}\frac{\partial\tilde{c}_{\mu m}}{\partial t}|\mu m\rangle=\sum_{\mu m}\hat{V}^{\text{F}}(t)|\mu m\rangle\tilde{c}_{\mu m}$$

(3)

where the Floquet Hamiltonian operator is defined as

$$\hat{V}^{\text{F}}(t)\equiv\hat{V}(t)-i\hbar\frac{\partial}{\partial t}.$$

(4)

Next, we close Eq. 3 by multiplying both sides by $\langle\nu|$ and write $\langle\nu|\hat{V}^{\text{F}}(t)|\mu m\rangle=\sum_{n}\widetilde{V}_{\nu n,\mu m}^{\text{F}}e^{in\omega t}$ as a Fourier series:

$$i\hbar\sum_{m}\frac{\partial\tilde{c}_{\nu m}}{\partial t}e^{im\omega t}=\sum_{\mu m}\sum_{n}\widetilde{V}_{\nu n,\mu m}^{\text{F}}e^{in\omega t}\tilde{c}_{\mu m}.$$

(5)

Thus, the TDSE in Eq. (3) can be solved by grouping together all terms with the same time dependence, leading to the equation of motion for $\tilde{c}$

$$i\hbar\frac{\partial}{\partial t}\tilde{c}_{\nu n}=\sum_{\mu m}\widetilde{V}_{\nu n,\mu m}^{\text{F}}\tilde{c}_{\mu m}.$$

(6)

The matrix elements of the Floquet Hamiltonian can be obtained by performing a Fourier transformation on the matrix elements

$$\widetilde{V}_{\nu n,\mu m}^{\text{F}}=\langle\langle\nu n|\hat{V}^{\text{F}}|\mu m\rangle\rangle=\frac{1}{\tau}\int_{0}^{\tau}dt\left\langle\nu\right|\hat{V}^{\text{F}}(t)\left|\mu\right\rangle e^{i(m-n)\omega t}.$$

(7)

Here, we define the double bracket projection of an electronic operator by $\langle\langle\nu n|\cdots|\mu m\rangle\rangle=\frac{1}{\tau}\int_{0}^{\tau}dt\left\langle\nu\right|\cdots\left|\mu\right\rangle e^{i(m-n)\omega t}$. Given that the electronic Hamiltonian operator $\hat{V}(t)$ is periodic in time, the double-bracket projection eliminates all time dependence and the Floquet Hamiltonian matrix reads

$$\widetilde{V}_{\nu n,\mu m}^{\text{F}}=\langle\langle\nu n|\hat{V}|\mu m\rangle\rangle+\delta_{\mu\nu}\delta_{mn}m\hbar\omega.$$

(8)

In the end, with this time-independent Hamiltonian, Eq. (6) can be formally solved by the exponential operator $\exp(-i\widetilde{\mathbf{V}}^{\text{F}}t/\hbar)$ with an arbitrary initial state.

At this point, we will allow nuclei to move and turn out attention to the equation of motion for the density matrix $\widehat{\rho}_{W}(\vec{R},\vec{P},t)$ within the Floquet diabatic basis. The Wigner-transformed density matrix in the Floquet diabatic representation is

$$\widetilde{A}_{\nu n,\mu m}^{\text{dia}}(\vec{R},\vec{P},t)=\langle\nu n|\widehat{\rho}_{W}(\vec{R},\vec{P},t)|\mu m\rangle.$$

(9)

For a proper F-QCLE, we will need to calculate the time derivative of $\widetilde{\mathbf{A}}^{\text{dia}}$

$$\frac{\partial}{\partial t}\widetilde{A}_{\nu n,\mu m}^{\text{dia}}=\left\langle\nu n\right|\frac{\partial\widehat{\rho}_{W}}{\partial t}\left|\mu m\right\rangle-i\left(n-m\right)\hbar\omega\widetilde{A}_{\nu n,\mu m}^{\text{dia}}$$

(10)

where the Floquet diabatic states depends on time explicitly. We begin by using Eq. (2) to project $\frac{\partial}{\partial t}\hat{\rho}_{W}$ into a Floquet diabatic basis. For the commutator term in Eq. (2), we can divide the operator product into matrix multiplication:

	$\displaystyle\langle\nu n|\left[\hat{V}_{W},\hat{\rho}_{W}\right]|\mu m\rangle$	$\displaystyle=$	$\displaystyle\sum_{\lambda l}\langle\langle\nu n|\hat{V}_{W}|\lambda l\rangle\rangle\widetilde{A}_{\lambda l,\mu m}^{\text{dia}}-\widetilde{A}_{\nu n,\lambda l}^{\text{dia}}\langle\langle\lambda l|\hat{V}_{W}|\mu m\rangle\rangle$		(11)
		$\displaystyle=$	$\displaystyle[\widetilde{\mathbf{V}}_{W},\widetilde{\mathbf{A}}^{\text{dia}}]_{\nu n,\mu m}.$

Here, we have inserted the identity of the diabatic electronic basis ($\hat{1}=\sum_{\lambda}|\lambda\rangle\langle\lambda|$) and expanded the time-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 can write the sum of both terms as $[\widetilde{\mathbf{V}}^{\text{F}},\widetilde{\mathbf{A}}^{\text{dia}}]$, i.e. we can replace $\widetilde{\mathbf{V}}_{W}$ with $\widetilde{\mathbf{V}}^{\text{F}}$.

For the anti-commutator term, we can use the same procedure to divide the operator product

	$\displaystyle\langle\nu n|\left\{\frac{\partial\hat{V}_{W}}{\partial R^{\alpha}},\frac{\partial\hat{\rho}_{W}}{\partial P^{\alpha}}\right\}|\mu m\rangle$	$\displaystyle=$	$\displaystyle\sum_{\lambda l}\langle\langle\nu n|\frac{\partial\hat{V}_{W}}{\partial R^{\alpha}}|\lambda l\rangle\rangle\frac{\partial\widetilde{A}_{\lambda l,\mu m}^{\text{dia}}}{\partial P^{\alpha}}$		(12)
			$\displaystyle+\frac{\partial\widetilde{A}_{\nu n,\lambda l}^{\text{dia}}}{\partial P^{\alpha}}\langle\langle\lambda l|\frac{\partial\hat{V}_{W}}{\partial R^{\alpha}}|\mu m\rangle\rangle$

where $\langle\nu n|\frac{\partial\widehat{\rho}_{W}}{\partial P^{\alpha}}|\mu m\rangle=\frac{\partial}{\partial P^{\alpha}}\widetilde{A}_{\nu n,\mu m}^{\text{dia}}$. Note that, since the Floquet diabatic states do not depend on the nuclear coordinate, we can rewrite the derivative of the electronic Hamiltonian in terms of the Floquet Hamiltonian

$$\langle\langle\nu n|\frac{\partial\hat{V}_{W}}{\partial R^{\alpha}}|\mu m\rangle\rangle=\frac{\partial}{\partial R^{\alpha}}\widetilde{V}_{\nu n,\mu m}^{\text{F}}.$$

(13)

In the end, we may combine the above expressions to write down a complete diabatic F-QCLE

	$\displaystyle\frac{\partial}{\partial t}\widetilde{\mathbf{A}}^{\text{dia}}$	$\displaystyle=$	$\displaystyle-\frac{i}{\hbar}\left[\widetilde{\mathbf{V}}^{\text{F}},\widetilde{\mathbf{A}}^{\text{dia}}\right]-\frac{P^{\alpha}}{M}\frac{\partial\widetilde{\mathbf{A}}^{\text{dia}}}{\partial R^{\alpha}}$		(14)
			$\displaystyle+\frac{1}{2}\left\{\frac{\partial\widetilde{\mathbf{V}}^{\text{F}}}{\partial R^{\alpha}},\frac{\partial\widetilde{\mathbf{A}}^{\text{dia}}}{\partial P^{\alpha}}\right\}$

As a final remark, we emphasize that the Floquet Hamiltonian $\widetilde{\mathbf{V}}^{\text{F}}=\widetilde{\mathbf{V}}^{\text{F}}(\vec{R})$ is a time-independent matrix, so Eq. (14) is simply the diabatic QCLE corresponding to an infinitely large electronic Hamiltonian $\widetilde{\mathbf{V}}^{\text{F}}$.

To recast the diabatic F-QCLE in an adiabatic representation, we diagonalize the Floquet Hamiltonian matrix by solving the eigenvalue problem:

$$\sum_{\nu n}\widetilde{V}_{\mu m,\nu n}^{\text{F}}(\vec{R})G_{\nu n}^{J}(\vec{R})={\cal E}_{J}^{\text{F}}(\vec{R})G_{\mu m}^{J}(\vec{R}).$$

(15)

The eigenvalues ${\cal E}_{J}^{\text{F}}={\cal E}_{J}^{\text{F}}(\vec{R})$ are the so-called Floquet quasi-energies with corresponding eigenvectors $G_{\mu m}^{J}(\vec{R}).$ Since $\widetilde{\mathbf{V}}^{\text{F}}$ is Hermitian, we can choose the eigenvectors $G_{\mu m}^{J}$ to be othornormal so that we have the identities $\mathbf{G}^{\dagger}\mathbf{G}=\mathbf{G}\mathbf{G}^{\dagger}=\mathbf{I}$, i.e. $\sum_{\lambda\ell}G_{\lambda\ell}^{J*}G_{\lambda\ell}^{K}=\delta_{JK}$ and $\sum_{L}G_{\mu m}^{L}G_{\nu n}^{L*}=\delta_{\mu m,\nu n}$. The Floquet adiabatic state corresponding to the quasi-energy ${\cal E}_{J}^{\text{F}}$ are

$$|\Phi^{J}(\vec{R},t)\rangle=\sum_{\mu m}G_{\mu m}^{J}(\vec{R})\left|\mu m\right\rangle.$$

(16)

As a practical matter, although $\widetilde{\mathbf{V}}^{\text{F}}$ is infinite dimensional, we can truncate highly oscillating Floquet states by replacing $\sum_{m=-\infty}^{\infty}$ with $\sum_{m=-m_{\text{max}}}^{m_{\text{max}}}$.

With this Floquet adiabatic state basis, the probability density can be obtained by a diabatic-to-adiabatic transformation

$$\widetilde{A}_{JK}^{\text{adi}}(\vec{R},\vec{P},t)=\langle\Phi^{J}|\widehat{\rho}_{W}\left(\vec{R},\vec{P},t\right)|\Phi^{K}\rangle=\left(\mathbf{G}^{\dagger}\widetilde{\mathbf{A}}^{\text{dia}}\mathbf{G}\right)_{JK}$$

(17)

in the Floquet adiabatic representation. Note that, since the eigenvectors $G_{\mu m}^{J}(\vec{R})$ do not depend on time explicitly, the time-derivative of the adiabatic probability density can be calculated simply to be:

$$\frac{\partial\widetilde{\mathbf{A}}^{\text{adi}}}{\partial t}=\mathbf{G}^{\dagger}\frac{\partial\widetilde{\mathbf{A}}^{\text{dia}}}{\partial t}\mathbf{G}.$$

(18)

We are now ready to derive the adiabatic F-QCLE in the following section.

Our first approach follows the order presented above: we first perform partial Wigner transformation, we second project to a Floquet diabatic basis, we third transform to an adiabatic representation. For the last step, following Eq. (18), we transform the diabatic F-QCLE by sandwiching the diabatic F-QCLE (Eq. (14)) with $\mathbf{G}^{\dagger}$ and $\mathbf{G}$. The first term on the right hand side of Eq. (14) (the commutator term) becomes

$\displaystyle\sum_{\nu n}\sum_{\mu m}G_{\nu n}^{J*}\left[\widetilde{V}^{\text{F}},\widetilde{A}^{\text{dia}}\right]_{\nu n,\mu m}G_{\mu m}^{K}$

$\displaystyle=$

$\displaystyle\left({\cal E}_{J}^{\text{F}}-{\cal E}_{K}^{\text{F}}\right)\widetilde{A}_{JK}^{\text{adi}}$

(19)

For the second term on the right hand side of Eq. (14), since the Floquet adiabatic states depend on the nuclear coordinate $\vec{R}$, the $R^{\alpha}$ derivative of the density must yield

$$\sum_{\nu n}\sum_{\mu m}G_{\nu n}^{J*}\frac{\partial\widetilde{A}_{\nu n,\mu m}^{\text{dia}}}{\partial R^{\alpha}}G_{\mu m}^{K}=\frac{\partial\widetilde{A}_{JK}^{\text{adi}}}{\partial R^{\alpha}}+\sum_{L}\left(D_{JL}^{\alpha}\widetilde{A}_{LK}^{\text{adi}}-\widetilde{A}_{JL}^{\text{adi}}D_{LK}^{\alpha}\right)$$

(20)

where the derivative coupling is $D_{JK}^{\alpha}=\langle\langle\Phi^{J}|\frac{\partial}{\partial R^{\alpha}}|\Phi^{K}\rangle\rangle=\sum_{\mu m}G_{\mu m}^{J*}\frac{\partial G_{\mu m}^{K}}{\partial R^{\alpha}}$ corresponding to the change of the Floquet adiabatic states with respect to the nuclear coordinate $R^{\alpha}$. Note that, if the Floquet Hamiltonian is real-valued, the diagonal element of the derivative coupling is zero ($D_{JJ}^{\alpha}=0$). For the third term on the right hand side of Eq. (14) (the anti-commutator term), the $R^{\alpha}$ derivative of the Floquet Hamiltonian can be written in terms of the force matrix

$$\sum_{\nu n}\sum_{\mu m}G_{\nu n}^{J*}\frac{\partial\widetilde{V}_{\nu n,\mu m}^{\text{F}}}{\partial R^{\alpha}}G_{\mu m}^{K}=-F_{JK}^{\alpha}$$

(21)

explicitly,

$$F_{JK}^{\alpha}=-\frac{\partial{\cal E}_{J}^{\text{F}}}{\partial R^{\alpha}}\delta_{JK}+({\cal E}_{J}^{\text{F}}-{\cal E}_{K}^{\text{F}})D_{JK}^{\alpha}.$$

The force matrix accounts for direct changes in the nuclear momentum associated with the electronic coupling. One can understand the diagonal element $F_{JJ}^{\alpha}=-\frac{\partial{\cal E}_{J}^{\text{F}}}{\partial R^{\alpha}}$ as the classical force for nuclear dynamics moving along the $J$-th Floquet quasi-energy surface in the phase space. The off-diagonal force term $F_{JK}^{\alpha}=({\cal E}_{J}^{\text{F}}-{\cal E}_{K}^{\text{F}})D_{JK}^{\alpha}$ for $J\neq K$ corresponds to nuclear velocity rescaling associated with non-adiabatic transitions between Floquet quasi-energy surfaces(Kapral and Ciccotti, 1999; Kapral, 2006).

Finally, the F-QCLE via the WFA approach reads

	$\displaystyle\frac{\partial}{\partial t}\widetilde{A}_{JK}^{\text{adi}}$	$\displaystyle=$	$\displaystyle-\frac{i}{\hbar}\left({\cal E}_{J}^{\text{F}}-{\cal E}_{K}^{\text{F}}\right)\widetilde{A}_{JK}^{\text{adi}}$		(22)
			$\displaystyle-\frac{P^{\alpha}}{M^{\alpha}}\frac{\partial\widetilde{A}_{JK}^{\text{adi}}}{\partial R^{\alpha}}-\frac{P^{\alpha}}{M^{\alpha}}\sum_{L}\left(D_{JL}^{\alpha}\widetilde{A}_{LK}^{\text{adi}}-\widetilde{A}_{JL}^{\text{adi}}D_{LK}^{\alpha}\right)$
			$\displaystyle-\frac{1}{2}\sum_{L}\left(F_{JL}^{\alpha}\frac{\partial\widetilde{A}_{LK}^{\text{adi}}}{\partial P^{\alpha}}+\frac{\partial\widetilde{A}_{JL}^{\text{adi}}}{\partial P^{\alpha}}F_{LK}^{\alpha}\right)$

We find that Eq. 22 takes exactly the same form as the standard QCLE in the adiabatic representation (for electron-nuclear dynamics without a driving laser).

For the second approach, we exchange the “adiabatic” and “Floquet” operations after the partial Wigner transformation. In this case, we directly project Eq. (2) onto the Floquet adiabatic basis $|\phi^{J}(\vec{R},t)\rangle$. Namely, we make the diabatic-to-adiabatic transform of the Floquet electronic basis prior to the projection onto the dressed electronic states. Thus, we consider this path as the Wigner-Adiabatic-Floquet (WAF) approach.

Overall, we apply a procedure similar to what was used in Eq. 22. For the commutator term, we divide operator products using the same technique as in Appendix A. The $R^{\alpha}$ derivative term yields a derivative coupling term as the Floquet adiabatic basis depends on $R^{\alpha}$ explicitly. In the end, the WAF approach includes the first three terms exactly as Eq. (22). However, from the anti-commutator term of Eq. (2), the WAF approach leads to an additional cross derivative force. To see this, we focus on the derivative of the electronic Hamiltonian in the adiabatic representation

$$\langle\langle\Phi^{J}|\frac{\partial\hat{V}}{\partial R^{\alpha}}|\Phi^{K}\rangle\rangle=\langle\langle\Phi^{J}|\frac{\partial\hat{V}^{\text{F}}}{\partial R^{\alpha}}|\Phi^{K}\rangle\rangle+i\hbar\langle\langle\Phi^{J}|\frac{\partial}{\partial R^{\alpha}}\frac{\partial}{\partial t}|\Phi^{K}\rangle\rangle$$

(23)

where we have used the definition of the Floquet Hamiltonian $\frac{\partial\hat{V}^{\text{F}}}{\partial R^{\alpha}}=\frac{\partial\hat{V}}{\partial R^{\alpha}}-i\hbar\frac{\partial}{\partial R^{\alpha}}\frac{\partial}{\partial t}$. The derivative of the electronic Hamiltonian yields two terms: first, the same force matrix we obtained in Eq. (21):

$$\langle\langle\Phi^{J}|\frac{\partial\hat{V}^{\text{F}}}{\partial R^{\alpha}}|\Phi^{K}\rangle\rangle=-F_{JK}^{\alpha}$$

(24)

second, a cross derivative force comes from the explicit dependence of the Floquet adiabatic states on both the nuclear coordinates and time

$$i\hbar\langle\langle\Phi^{J}|\frac{\partial}{\partial R^{\alpha}}\frac{\partial}{\partial t}|\Phi^{K}\rangle\rangle=-\chi_{JK}^{\alpha},$$

and explicitly,

$$\chi_{JK}^{\alpha}=\sum_{\mu m}m\hbar\omega G_{\mu m}^{J*}\frac{\partial G_{\mu m}^{K}}{\partial R^{\alpha}}.$$

(25)

Note that, unlike the derivative coupling $D_{JK}^{\alpha}$, the diagonal element of $\chi_{JK}^{\alpha}$ is non-zero and real-valued.

Finally, in the same Floquet adiabatic basis as given by Eq. (16), the F-QCLE reads:

	$\displaystyle\frac{\partial}{\partial t}\widetilde{A}_{JK}^{\text{adi}}$	$\displaystyle=$	$\displaystyle-\frac{i}{\hbar}\left({\cal E}_{J}^{\text{F}}-{\cal E}_{K}^{\text{F}}\right)\widetilde{A}_{JK}^{\text{adi}}$		(26)
			$\displaystyle-\frac{P^{\alpha}}{M^{\alpha}}\frac{\partial\widetilde{A}_{JK}^{\text{adi}}}{\partial R^{\alpha}}-\frac{P^{\alpha}}{M^{\alpha}}\sum_{L}\left(D_{JL}^{\alpha}\widetilde{A}_{LK}^{\text{adi}}-\widetilde{A}_{JL}^{\text{adi}}D_{LK}^{\alpha}\right)$
			$\displaystyle-\frac{1}{2}\sum_{L}\left(F_{JL}^{\alpha}+\chi_{LJ}^{\alpha*}\right)\frac{\partial\widetilde{A}_{LK}^{\text{adi}}}{\partial P^{\alpha}}+\frac{\partial\widetilde{A}_{JL}^{\text{adi}}}{\partial P^{\alpha}}\left(F_{LK}^{\alpha}+\chi_{LK}^{\alpha}\right)$

We observe that, while Eq. (22) and Eq. (26) take the same form, the “effective” force matrix (defined as the coefficients of $\frac{\partial\widetilde{\mathbf{A}}^{\text{adi}}}{\partial P^{\alpha}}$) includes an additional cross derivative force that indicates the difference between these two equations. In other words, even though these two QCLEs will propagate electrons equivalently, the difference in the “effective” force matrix will lead to different nuclear dynamics in phase space. Specifically, from Eq. 26, the classical force on the $J$-th quasi-energy surface is given by $F_{JJ}^{\alpha}+\chi_{JJ}^{\alpha}$, implying that the nuclear dynamics should experience the additional cross derivative force on top of the Floquet quasi-energy surface.

Finally, for completeness, we should emphasize that the first derivation of an F-QCLE was published by Horenko, Schmidt, and Schütte who used an approach that was different from either of the approaches above. In Ref. (Horenko, Schmidt, and Schütte, 2001), these authors derived their F-QCLE using an Adiabatic-then-Wigner (AW) approach. Their final equation of motion was:

	$\displaystyle\frac{\partial}{\partial t}\widetilde{A}_{JK}^{\text{adi}}$	$\displaystyle=$	$\displaystyle-\frac{i}{\hbar}\left({\cal E}_{J}^{\text{F}}-{\cal E}_{K}^{\text{F}}\right)\widetilde{A}_{JK}^{\text{adi}}$		(27)
			$\displaystyle-\frac{P^{\alpha}}{M^{\alpha}}\frac{\partial\widetilde{A}_{JK}^{\text{adi}}}{\partial R^{\alpha}}-\frac{P^{\alpha}}{M^{\alpha}}\sum_{L}\left(D_{JL}^{\alpha}\widetilde{A}_{LK}^{\text{adi}}-\widetilde{A}_{JL}^{\text{adi}}D_{LK}^{\alpha}\right)$
			$\displaystyle+\frac{1}{2}\sum_{L}\left(\frac{\partial{\cal E}_{J}^{\text{F}}}{\partial R^{\alpha}}+\frac{\partial{\cal E}_{K}^{\text{F}}}{\partial R^{\alpha}}\right)\frac{\partial\widetilde{A}_{JK}^{\text{adi}}}{\partial P^{\alpha}}$

As mentioned above, in the limit of a time-independent Hamiltonian, Izmaylov and co-workers have shown that the AW approach results in the loss of geometric phase related features for a 2D problem(Ryabinkin et al., 2014); the WA approach is far more accurate (see Figs in Ref. (Ryabinkin et al., 2014)). Moreover, in contrast to Eq. (22) and Eq. (26), the last term in Eq. 27 does not include any off-diagonal force term which implies no momentum rescaling. Within a surface hopping context, a lack of momentum rescaling will inevitably lead to a violation of detailed balance(Parandekar and Tully, 2006, 2005; Schmidt, Parandekar, and Tully, 2008). Note that early numerical assessments of the AW approach did not investigate either geometric phase effects or detailed balance(Ando, 2002; Ando and Santer, 2003). For these reasons, we expect the AW approach to the F-QCLE to be less accurate than the WFA and WAF approaches. Nevertheless, to satisfy the readers who is curious about the original AW approach, we will present the relevant transmission and reflection probabilities and final momentum distribution data in Appendix B which should remind one why momentum jumps are so important for scattering calculations¹¹1Interestingly, the AW approach to the time-independent QCLE was found to perform well for short times for small closed systems in 2003(Ando, 2002; Ando and Santer, 2003), but we are unaware of more systematic tests of the method beyond Ref. (Ryabinkin et al., 2014), especially within the context of long time dynamics..

To analyze the difference between these approaches, we consider a shifted avoided crossing model composed of a two-level electronic system coupled to a 1D nuclear motion. The diabatic electronic states are denoted as $|g\rangle$ and $|e\rangle$ and the electronic Hamiltonian is given by

$$\hat{V}(R,t)=\left(\begin{array}[]{cc}V_{gg}(R)&V_{ge}(R,t)\\ V_{eg}(R,t)&V_{ee}(R)\end{array}\right)$$

(28)

where the diabatic energy is given by

$$V_{gg}(R)=\begin{cases}A(1-e^{-BR})&R>0\\ -A(1-e^{BR})&R<0\end{cases}$$

(29)

$$V_{ee}(R)=-V_{gg}(R)+\hbar\omega$$

(30)

and the diabatic coupling is periodic in time

$$V_{ge}(R,t)=V_{eg}(R,t)=Ce^{-DR^{2}}\cos\omega t.$$

(31)

The parameters are $A=0.01$, $B=1.6$, $C=0.01$, $D=1.0$, 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 is shifted by $\hbar\omega$. For simplicity, we choose the laser frequency $\hbar\omega=0.024$ large enough ($\hbar\omega>2A$) so that the diabatic Floquet states $|gm\rangle$ have an avoided crossing only with $|e(m-1)\rangle$ (at $R=0$) and do not have any trivial crossings. The size of the Floquet basis is truncated by $m_{\text{max}}=4$ as the laser field is weak ($C/\hbar\omega<1$).

We assume the initial wavepacket is a Gaussian centered at the initial position $R_{0}$ and momentum $P_{0}$ on diabat $|g\rangle$:

$$|\Psi_{0}\rangle=\frac{1}{{\cal N}}\exp\left(-\frac{(R-R_{0})^{2}}{2\sigma^{2}}+\frac{i}{\hbar}P_{0}(R-R_{0})\right)|g\rangle$$

(32)

where the normalization factor is ${\cal N}^{4}=\pi\sigma^{2}$. The width of the Gaussian is chosen to be $\sigma=20\hbar/P_{0}$. The wavepacket can be propagated exactly in the diabatic representation.

To show the difference between the dynamics as obtained by the different F-QCLEs, we will simulate F-FSSH results for both the WFA and WAF approaches. Within F-FSSH, we describe the propagation of the Floquet wavepacket by a swarm of trajectories, each with its own electronic amplitudes $\tilde{c}_{J}$ satisfying

$$\frac{\partial\tilde{c}_{J}}{\partial t}=-\frac{i}{\hbar}{\cal E}_{J}^{\text{F}}\tilde{c}_{J}-\frac{P}{M}\sum_{L}D_{JL}\tilde{c}_{L}$$

(33)

where ${\cal E}_{J}^{\text{F}}={\cal E}_{J}^{\text{F}}(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

$$\frac{dR}{dt}=\frac{P}{M}$$

(34)

$$\frac{dP}{dt}=\begin{cases}-\frac{\partial{\cal E}_{J}^{\text{F}}}{\partial R}&\text{for WFA}\\ -\frac{\partial{\cal E}_{J}^{\text{F}}}{\partial R}-\chi_{JJ}&\text{for WAF}\end{cases}$$

(35)

Here, based on the connection between the QCLE and FSSH, the nuclear force in Eq. 35 is determined 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 Floquet state $J$ to state $K$ is given by

$$\text{Prob}(J\rightarrow K)=-2\text{Re}\left(\frac{P^{\alpha}}{M^{\alpha}}\cdot D_{KJ}^{\alpha}\frac{\tilde{c}_{J}\tilde{c}_{K}^{*}}{|\tilde{c}_{J}|^{2}}\right)dt$$

where $dt$ is the classical time step. After each successful hop, the velocity is adjusted to conserve the total Floquet quasi-energy. If a frustrated hop occurs, we implement velocity reversal(Jasper and Truhlar, 2003). Note that we neglect the decoherence correction since the over-coherence problem should not be severe for a simple avoided crossing model(Landry and Subotnik, 2012).

In the end, the probability to measure diabatic state $\mu$ can be evaluated by the density matrix interpretation(Landry, Falk, and Subotnik, 2013)

$$P_{\mu}=\sum_{m}\frac{N(\mu m)}{N_{\text{traj}}}+\sum_{n\neq m}\frac{\sum_{l}^{N(\mu m)}\sum_{k}^{N(\mu n)}\tilde{c}_{\mu m}^{(l)}\tilde{c}_{\mu n}^{(k)*}e^{i(m-n)\omega t}}{N(\mu m)N(\mu n)}$$

(36)

where $N(\mu m)=\sum_{l}^{N_{\text{traj}}}\delta_{J^{(l)}\mu m}$ is the number of the trajectories that have the active surface $J^{(l)}$ end up on the Floquet state $|\mu m\rangle$. Here $l$ and $k$ are the trajectory indices. We propagate $N_{\text{traj}}$ trajectories for an amount of time long enough for each trajectory to pass through the coupling region ($|R|<3$ for this parameter set).

First, we analyze the effective potential energy surfaces for nuclear dynamics by integrating Eq. 35 over $R$ for the WFA and WAF approaches respectively. For the WFA approach, the effective PES simply recovers the Floquet quasi-energy surfaces ${\cal E}_{J}^{\text{F}}$. For the WAF approach, the effective PES is obtained by $V_{eff}(R)={\cal E}_{J}^{\text{F}}(R)+\int_{-\infty}^{R}\chi_{JJ}(R^{\prime})dR^{\prime}$ where the quasi-energy surface is modified by the integration of the cross derivative force. We find that including the cross derivative force in the WAF approach increases the crossing barrier for the nuclear dynamics on the lower adiabat (see Fig. 1). Note that, in terms of an F-FSSH calculation, these changes will have a direct effect on the nuclear dynamics, but not the electronic amplitudes.

Next, we turn our attention to the transmission and reflection probabilities produced by the F-FSSH calculations. Overall, the WFA results are more accurate than the WAF results. In Fig. 2, we find that, according to the WFA approach, there should be a rise in transmission on the lower adiabat around $P_{0}\approx 5.3$, which is the momentum for which transmission should be allowed classically; see the barrier height ($\approx 0.007$) in Fig. 1. Indeed, such a threshold at $P_{0}\approx 5.3$ is found according to exact wavefunction simulation as well. However, for the WAF approach with the cross derivative force, one find a higher crossing barrier energy ($\approx 0.015$), and the transmission on the lower adiabat occurs (incorrectly) at $P_{0}\approx 7.8$. This result suggests that the cross derivative force is a fictitious term: the WAF semiclassical 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_{0}\approx 8.9$, which agrees with the classical energy difference $2A=0.02$ (see Eq. (29) for the definition of $A$). Second, for high initial momentum ($P_{0}>8.0$), the F-FSSH-WFA can almost recover the correct nuclear dynamics. Third, in the intermediate momentum region $P_{0}\in(6,8)$, unfortunately, the F-FSSH wavepacket exhibits less transmission than the exact calculation. This discrepancy may be attributed to FSSH’s inability to capture nuclear tunneling effects.