Nonadiabatic Dynamics Near Metal Surface With Periodic Drivings: A Floquet Surface Hopping Algorithm

TL;DR

This paper introduces a Floquet surface hopping method to simulate nonadiabatic molecular dynamics near metal surfaces under periodic light driving, effectively capturing fast oscillations and steady states.

Contribution

The paper develops a novel Floquet surface hopping approach based on a Floquet quantum master equation, enabling efficient simulation of driven nonadiabatic dynamics near metal surfaces.

Findings

FaSH-density algorithm accurately reproduces FQME results

Method captures fast oscillations and steady states

Applicable to strong light-matter interaction studies

Abstract

We develop a Floquet surface hopping (FSH) approach to deal with nonadiabatic dynamics of molecules near metal surfaces subjected to time-periodic drivings from strong light-matter interactions. The method is based on a Floquet classical master equation (FCME) derived from a Floquet quantum master equation (FQME), followed by a Wigner transformation to treat nuclear motion classically. We then propose different trajectory surface hopping algorithms to solve the FCME. We find that a Floquet averaged surface hopping with electron density (FaSH-density) algorithm works the best as benchmarked with the FQME, capturing both the fast oscillations due to the driving as well as the correct steady state observables. This method will be very useful to study strong light-matter interactions with a manifold of electronic state.

Equations

$$\hat{H}(t)={\hat{H}}_S(t)+{\hat{H}}_B+{\hat{H}}_T$$

$$\begin{array}{rl}{\hat{H}}_S(t)=(E_d+A\sin (\Omega t))d^{+}d+& g(a^{+}+a)d^{+}d+\\ & \hslash \omega (a^{+}a+\frac{1}{2})\end{array}$$

$${\hat{H}}_B=\sum _k{ϵ}_k{c}_k^{+}{c}_k$$

$${\hat{H}}_T=\sum _k{V}_k(d^{+}{c}_k+c_k^{+}d)$$

$${\hat{H}}_S=(E_d+A\sin (\Omega t))d^{+}d+\sqrt{\frac{2m\omega }{\hslash }}gx{d}^{+}d+\frac{p^2}{2m}+\frac{1}{2}m{\omega }^2{x}^2$$

$${\hat{H}}_S(t)={\hat{H}}_{mol}+A\sin (\Omega t)d^{+}d$$

$${\hat{H}}_{mol}=(E_d+\sqrt{\frac{2m\omega }{\hslash }}gx)d^{+}d+\frac{p^2}{2m}+\frac{1}{2}m{\omega }^2{x}^2$$

$${\hat{H}}_{mol}=H_0|0⟩⟨0|+H_1|1⟩⟨1|$$

$${\hat{H}}_0=\frac{p^2}{2m}+\frac{1}{2}m{\omega }^2{x}^2$$

$${\hat{H}}_1=\frac{p^2}{2m}+\frac{1}{2}m{\omega }^2{x}^2+\sqrt{\frac{2m\omega }{\hslash }}gx+E_d$$

$$\frac{d\rho (t)}{dt}=-\frac{i}{\hslash }[H_T(t),\rho (0)]-\frac{1}{{\hslash }^2}{\int }_0^{t}d{t}^{\prime }[H_T(t),[H_T(t^{\prime }),\rho (t^{\prime })]]$$

$$H_T(t)=U^{+}(t)H_{T}U(t)$$

$$\begin{array}{rl}U(t)& =\hat{T}\exp \left[-\frac{i}{\hslash }{\int }_0^{t}d{t}^{\prime }(H_B+H_S(t^{\prime }))\right]\\ & =\hat{T}\exp \left[-\frac{i}{\hslash }{\int }_0^{t}d{t}^{\prime }(H_B+H_{mol}+A\sin (\Omega t^{\prime })d^{+}d)\right]\\ & =\exp (-i(H_B+H_{mol})t/\hslash -ig(t)d^{+}d/\hslash )\end{array}$$

$$\begin{array}{rl}\frac{d\rho (t)}{dt}=& -\frac{i}{\hslash }[H_T(t),\rho (0)]-\\ & \frac{1}{{\hslash }^2}{\int }_0^{\infty }d\tau [H_T(t),[H_T(t-\tau ),{\rho }_B^{eq}⨂{\rho }_S(t)]]\end{array}$$

$$\begin{array}{r}\frac{d{\rho }_S(t)}{dt}=-\frac{1}{{\hslash }^2}{\int }_0^{\infty }d\tau T{r}_B[H_T(t),[H_T(t-\tau ),{\rho }_B^{eq}⨂{\rho }_S(t)]]\end{array}$$

$$\begin{array}{rl}\frac{d{\rho }_0}{dt}=& -i[H_0,{\rho }_0]-\sum _k\frac{|V_k{|}^2}{{\hslash }^2}\times \\ & {\int }_0^{\infty }d\tau [e^{i{ϵ}_k\tau /\hslash -i(g(t)-g(t-\tau ))/\hslash }f({ϵ}_k)e^{-i{H}_1\tau /\hslash }{e}^{i{H}_0\tau /\hslash }{\rho }_0\\ & -e^{i{ϵ}_k\tau /\hslash -i(g(t)-g(t-\tau ))/\hslash }(1-f({ϵ}_k)){\rho }_1{e}^{-i{H}_1\tau /\hslash }{e}^{i{H}_0\tau /\hslash }\\ & +e^{-i{ϵ}_k\tau /\hslash +i(g(t)-g(t-\tau ))/\hslash }f({ϵ}_k){\rho }_0{e}^{-i{H}_0\tau /\hslash }{e}^{i{H}_1\tau /\hslash }\\ & -e^{-i{ϵ}_k\tau /\hslash +i(g(t)-g(t-\tau ))/\hslash }(1-f({ϵ}_k))e^{-i{H}_0\tau /\hslash }{e}^{i{H}_1\tau /\hslash }{\rho }_1]\end{array}$$

$$\begin{array}{rl}\frac{d{\rho }_1}{dt}=& -i[H_1,{\rho }_1]-\sum _k\frac{|V_k{|}^2}{{\hslash }^2}\times \\ & {\int }_0^{\infty }d\tau [e^{-i{ϵ}_k\tau /\hslash +i(g(t)-g(t-\tau ))/\hslash }(1-f({ϵ}_k))e^{-i{H}_0\tau /\hslash }{e}^{i{H}_1\tau /\hslash }{\rho }_1\\ & -e^{-i{ϵ}_k\tau /\hslash +i(g(t)-g(t-\tau ))/\hslash }f({ϵ}_k){\rho }_0{e}^{-i{H}_0\tau /\hslash }{e}^{i{H}_1\tau /\hslash }\\ & +e^{i{ϵ}_k\tau /\hslash -i(g(t)-g(t-\tau ))/\hslash }(1-f({ϵ}_k)){\rho }_1{e}^{-i{H}_1\tau /\hslash }{e}^{i{H}_0\tau /\hslash }\\ & -e^{i{ϵ}_k\tau /\hslash -i(g(t)-g(t-\tau ))/\hslash }f({ϵ}_k)e^{-i{H}_1\tau /\hslash }{e}^{i{H}_0\tau /\hslash }{\rho }_0]\end{array}$$

$$\begin{array}{r}{e}^{ig(t)}=e^{\frac{i}{\hslash }\frac{A}{\Omega }}{e}^{\frac{i}{\hslash }\frac{A}{\Omega }\cos (\Omega t)}=e^{\frac{i}{\hslash }\frac{A}{\Omega }}\sum _{n=-\infty }^{+\infty }(i{)}^n{J}_n(z)e^{in\Omega t}\end{array}$$

$$\begin{array}{r}{e}^{-i(g(t)-g(t-\tau ))/\hslash }=\sum _{n,m}(i{)}^n(-i{)}^m{J}_n(z)J_m(z)e^{i(n-m)\Omega t}{e}^{im\Omega \tau }\end{array}$$

$$\begin{array}{rl}\frac{d{\rho }_0(i,j)}{dt}=& -\frac{i}{\hslash }({ϵ}_0(i)-{ϵ}_0(j)){\rho }_0(i,j)\\ & -\frac{\Gamma }{2\hslash }\sum _{i^{\prime },k}\tilde{f}({ϵ}_1(i^{\prime })-{ϵ}_0(k))F_{i\to i^{\prime }}{F}_{k\to i^{\prime }}{\rho }_0(k,j)\\ & +\frac{\Gamma }{2\hslash }\sum _{i^{\prime },j^{\prime }}(1-\tilde{f}({ϵ}_1(j^{\prime })-{ϵ}_0(j)))F_{i\to i^{\prime }}{F}_{j\to j^{\prime }}{\rho }_1(i^{\prime },j^{\prime })\\ & -\frac{\Gamma }{2\hslash }\sum _{i^{\prime },k}{\rho }_0(i,k)\tilde{f}({ϵ}_1(i^{\prime })-{ϵ}_0(k))F_{j\to i^{\prime }}{F}_{k\to i^{\prime }}\\ & +\frac{\Gamma }{2\hslash }\sum _{i^{\prime },j^{\prime }}{\rho }_1(i^{\prime },j^{\prime })(1-\tilde{f}({ϵ}_1(i^{\prime })-{ϵ}_0(i)))F_{i\to i^{\prime }}{F}_{j\to j^{\prime }}\end{array}$$

$$\begin{array}{rl}\frac{d{\rho }_1(i^{\prime },j^{\prime })}{dt}=& -\frac{i}{\hslash }({ϵ}_1(i^{\prime })-{ϵ}_1(j^{\prime })){\rho }_1(i^{\prime },j^{\prime })\\ & -\frac{\Gamma }{2\hslash }\sum _{i,k^{\prime }}(1-\tilde{f}({ϵ}_1(k^{\prime })-{ϵ}_0(i)))F_{i\to i^{\prime }}{F}_{i\to k^{\prime }}{\rho }_1(k^{\prime },j^{\prime })\\ & +\frac{\Gamma }{2\hslash }\sum _{i,j}\tilde{f}({ϵ}_1(j^{\prime })-{ϵ}_0(j))F_{i\to i^{\prime }}{F}_{j\to j^{\prime }}{\rho }_0(i,j)\\ & -\frac{\Gamma }{2\hslash }\sum _{i,k^{\prime }}{\rho }_1(i^{\prime },k^{\prime })(1-\tilde{f}({ϵ}_1(k^{\prime })-{ϵ}_0(i)))F_{i\to j^{\prime }}{F}_{i\to k^{\prime }}\\ & +\frac{\Gamma }{2\hslash }\sum _{i,j}{\rho }_0(i,j)\tilde{f}({ϵ}_1(i^{\prime })-{ϵ}_0(i))F_{i\to i^{\prime }}{F}_{j\to j^{\prime }}\end{array}$$

$$\begin{array}{r}{\overline{E}}_d\equiv E_d-E_r\end{array}$$

$$\begin{array}{r}{F}_{i\to i^{\prime }}=⟨i^{\prime }⟩i=\int dx{\varphi }_i(x+\sqrt{2}\lambda ){\varphi }_i(x),\lambda \equiv g/\hslash \omega \end{array}$$

$$\begin{array}{r}{F}_{i\to i^{\prime }}=(p!/Q!{)}^{1/2}{\lambda }^{Q-p}{e}^{-{\lambda }^2/2}{L}_p^{Q-p}({\lambda }^2)[sgn(i^{\prime }-i){]}^{i-i^{\prime }}\end{array}$$

$$\begin{array}{r}\tilde{f}(x)=\sum _{n,m}(i{)}^n(-i{)}^m{J}_n(z)J_m(z)e^{i(n-m)\Omega t}f(x-m\Omega )\end{array}$$

$$\begin{array}{r}\overline{\tilde{f}}(x)=\sum _n{J}_n(z{)}^{2}f(x-n\Omega )\end{array}$$

$$\begin{array}{r}\Gamma (ϵ)=2\pi \sum _k|V_k{|}^2\delta ({ϵ}_k-ϵ)\end{array}$$

$$\begin{array}{rl}\frac{\partial P_0(x,p)}{\partial t}=& \frac{\partial H_0(x,p)}{\partial x}\frac{\partial P_0(x,p)}{\partial p}-\frac{\partial H_0(x,p)}{\partial p}\frac{\partial P_0(x,p)}{\partial x}\\ & -{\gamma }_{0\to 1}{P}_0(x,p)+{\gamma }_{1\to 0}{P}_1(x,p)\end{array}$$

$$\begin{array}{rl}\frac{\partial P_1(x,p)}{\partial t}=& \frac{\partial H_1(x,p)}{\partial x}\frac{\partial P_1(x,p)}{\partial p}-\frac{\partial H_1(x,p)}{\partial p}\frac{\partial P_1(x,p)}{\partial x}\\ & +{\gamma }_{0\to 1}{P}_0(x,p)-{\gamma }_{1\to 0}{P}_1(x,p)\end{array}$$

$$\begin{array}{r}{\gamma }_{0\to 1}=\frac{\Gamma }{\hslash }\tilde{f}(\Delta V)\end{array}$$

Full text

††preprint: AIP/123-QED

I Introduction

Floquet engineering is a powerful tool to control the dynamics of quantum systems using time-periodic external fields (e.g. electromagnetic fields), which is account for various phenomenaOka and Kitamura (2019). In practice, surface plasmon or a Fabry-Pérot cavity can create a resonant electromagnetic (EM) field and realize Floquet engineering through strong light-matter interactions. For instance, strong light-matter coupling is achieved by placing materials in a resonance EM field and show a great promise to control properties of matter without structural modificationsHertzog et al. (2019); Garcia-Vidal, Ciuti, and Ebbesen (2021). In addition, n-type or p-type organic semiconductors deposited on top of a metal surface, which acts as an open electromagnetic cavity, have been observed a noteworthy conductivity enhancementOrgiu et al. (2015); Nagarajan et al. (2020). Moreover, strong light-matter coupling can largely extend spatial range of energy transport in organic materialsGonzalez-Ballestero et al. (2015); Hou et al. (2020); Zhong et al. (2017). Furthermore, numerous photophysical and photochemical properties, such as electronic relaxation pathwaysLidzey et al. (2002); Coles et al. (2011, 2013), quantum yieldsWang et al. (2014); Ballarini et al. (2014); Grant et al. (2016); Stranius, Hertzog, and Börjesson (2018),, and chemical reactivitiesHutchison et al. (2012); Lin et al. (2016); Feist, Galego, and Garcia-Vidal (2018); Li, Mandal, and Huo (2021), can also be manipulated in the strong coupling regime. The microscopic mechanisms leading to these modifications through strong light-matter interactions, however, remain largely unexplained.

On the one hand, without Floquet driving, when nuclei interact non-adiabatically with a manifold electronic states (e.g. an adsorbate at a metal surface), there is a drastic breakdown of the Born-Oppenheimer approximation. There are several numerically exact solutions for such non-adiabatic processes when only a few degrees of freedom (DoFs) are considered, including multiconfiguration time dependent Hartree (MCTDH)Thoss, Kondov, and Wang (2007), quantum Monte Carlo (QMC)Mühlbacher and Rabani (2008),numerical renormalization group (NRG)Hewson and Meyer (2001); Cabra, Franco, and Galperin (2020); Wilson (1975), and hierarchical quantum master equation (HQME)Schinabeck et al. (2016a, b). That being said, these methods are difficult to apply to large, realistic systems. Many attempts have been given to develop approximate methods to deal with realistic systems involving nonadiabatic dynamics in open quantum systemsGalperin, Ratner, and Nitzan (2007); Reimers et al. (2007); Bode et al. (2012); Akimov and Prezhdo (2013). On the other hand, when considering nonadiabatic dynamics under Flquet driving, there exists Floquet surface hoppingFiedlschuster, Handt, and Schmidt (2016); Zhou et al. (2020); Chen, Zhou, and Subotnik (2020), as well as coupled-trajectory mixed quantum–classicalSchirò, Eich, and Agostini (2021) methods to study dynamics for a closed system. For open quantum systems, there are far less approaches available with only a few exceptions, such as Floquet scattering matrixWu and Cao (2006); Moskalets (2014); Pantazopoulos and Stefanou (2019), Floquet Green’s functionsCapolino et al. (2000); Torres (2005); Wu and Cao (2010); Honeychurch and Kosov (2023), and Floquet master equationsKohler, Dittrich, and Hänggi (1997); Wu and Timm (2010); Kuperman, Nagar, and Peskin (2020). Again, these methods are only applicable to small systems. Accurate methods are needed to study large and realistic systems.

In this article, we attempt to offer an accurate and efficient tools to study nonadiabatic dynamics in open quantum systems with Floquet engineering. We first derive a Floquet quantum master equation (FQME) to accurately uncover the non-adiabatic processes of a periodically driven impurity level near a metal surface. Then a corresponding Floquet classical master equation (FCME) is achieved by taking the Wigner transformation of FQME, similar to the steps in Ref. Dou, Nitzan, and Subotnik (2015) for non-Floquet scenario. We then propose three Floquet surface hopping algorithms to solve the FCME. We find that an algorithm named time-averged Floquet surface hopping with density (FaSH-density) is capable to reproduce fully dynamics as well as steady state results as benchmarked against FQME.

An outline of this paper is as follows. In Sec. II, we outline how to derive the FQME and FCME, and present the FSH algorithm. In Sec. III, we present results of both electronic population and phonon relaxation dynamics under different Floquet drivings employing five different methods (FQME, FSH, FaQME, FaSH, and FaSH-density). We conclude in Sec. IV.

II Theory

II.1 Model Hamiltonian

We consider the Anderson-Holstein (AH) Hamiltonian with one molecule coupled to a Fermionic bath as well as subjected to periodic driving. The molecule consists of a single electronic level coupled to a single phonon and a manifold of electrons:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

$\hat{H}_{S}(t)$ is the system Hamiltonian with $E_{d}$ being the energy level and $\omega$ being the frequency of the harmonic oscillator. $g$ represents the electron-phonon (el-ph) coupling strength. $A$ is the driving amplitude and $\Omega$ is driving frequency. In addition, $\hat{H}_{B}$ represents the bath Hamiltonian with $\epsilon_{k}$ being the energy level of the Fermion $c_{k}$ in the bath. $\hat{H}_{T}$ is the interaction Hamiltonian with $V_{k}$ being the coupling between the impurity level $d$ and the Fermion $c_{k}$ in the bath.

It will be convenient to replace $a^{+}$ and $a$ with the position $x$ and momentum $p$ (i.e. $x=\sqrt{\frac{\hbar}{2m\omega}}(a^{+}+a)$ and $p=i\sqrt{\frac{m\hbar\omega}{2}}(a^{+}-a)$), where $m$ is the nuclear mass, such that the system Hamiltonian can be written as

[TABLE]

To derive a Floquet quantum master equation for treating the dynamics of the periodically driven open quantum system, we first separate the time independent part from the time dependent part in $H_{S}(t)$, such that

[TABLE]

[TABLE]

We can further write $H_{mol}$ as

[TABLE]

[TABLE]

[TABLE]

where $\ket{0}$ ($\ket{1}$) denotes unoccupied (occupied) state of the impurity.

II.2 Floquet Quantum Master Equation

We now sketch out the derivation of a Floquet Quantum master equation in treating the dynamics of the system. The key quantity of interest is the reduced density matrix of the system. Starting with the quantum Liouville equation in the interaction picture, we find the total density matrix can be expressed as

[TABLE]

where

[TABLE]

Here, the time evolution operator $U(t)$ is

[TABLE]

where $\hat{T}$ is the time ordering operator and $g(t)=\frac{A}{\hbar\Omega}(1-\cos(\Omega t))$

In the Born-Markovian approximation, we replace $\rho(t^{\prime})$ in Eq. (11) by $\rho_{B}^{eq}\bigotimes\rho_{S}(t)$ in the integrand that relies on the assumptions that the bath remains in equilibrium throughout the process and that bath correlation functions decay fast on the system time scale. This leads to (setting $\tau=t-t^{\prime}$):

[TABLE]

Next, we assume the initial total density matrix is a direct product of the system density matrix and the equilibrium bath density matrix, i.e., $\rho(0)=\rho_{B}^{eq}\bigotimes\rho_{S}(0)$ and take the trace of Eq. (11) over bath degrees of freedom. We also use $Tr_{B}(H_{T}(t)\rho_{B}^{eq})=0$, which yields

[TABLE]

Back to Schr$\rm\ddot{o}$dinger picture, $\rho_{S}(t)=e^{iH_{mol}t}\rho_{s}e^{-iH_{mol}t}$, and assuming $\rho_{s}=\rho_{0}\ket{0}\bra{0}+\rho_{1}\ket{1}\bra{1}$, which ensures that there will be no coherence between occupied and unoccupied states at later time. The reduced density matrix for state 0 (unocupied) and state 1 (occupied) evolves as:

[TABLE]

[TABLE]

Here, $f$ is the Fermi function ($f(x)=1/(1+e^{\beta x})$). By employing the Jacobi-Anger expansion, we can express $e^{ig(t)}$ as:

[TABLE]

where $n$ is the integer, $J_{n}(z)$ is the $n$-th Bessel function of the first kind, and $z=\frac{A}{\hbar\Omega}$. Thus, we can expand the term $e^{-i(g(t)-g(t-\tau))/\hbar}$ appears in Eqs. 16 and 17 as:

[TABLE]

We now expand the reduced density matrix in a basis of harmonic oscillator eigenstates ($H_{0}|i\rangle=\epsilon_{0}(i)|i\rangle$, $H_{1}|i^{\prime}\rangle=\epsilon_{1}(i^{\prime})|i^{\prime}\rangle$,

[TABLE]

[TABLE]

The above equations are referred to as our Floquet Quantum Master equation. Here $\epsilon_{0}(i)=\hbar\omega(i+\frac{1}{2})$, and $\epsilon_{1}(i^{\prime})=\hbar\omega(i^{\prime}+\frac{1}{2})+\overline{E}_{d}$. $\overline{E}_{d}$ is the renormalized impurity energy level.

[TABLE]

where $E_{r}\equiv g^{2}/\hbar\omega$ is the reorganization energy. $F$ is the Frank-Condon factor.

[TABLE]

where $\phi_{i}(x)$ is the $i$th eigenfunction of the harmonic oscillator. The Frank-Condon factor can be expressed as

[TABLE]

p(Q) is the minimum (maximum) of $i$ and $i^{\prime}$, and $L_{n}^{m}$ is generalized Laguerre polynomial. $\tilde{f}(x)$ in the above equations is the modified Fermi function with Floquet replicas, which is given by

[TABLE]

Here we have defined $z=\frac{A}{\hbar\Omega}$. Notice that $\tilde{f}(x)$ is time dependent. When taking time average over a cycle, we arrive at time-independent modified Fermi function

[TABLE]

Finally, $\Gamma$ in Eqs. 20 and 21 is the hybridization function,

[TABLE]

In the wide band limit, we assume that $\Gamma$ is a constant (i.e., does not change with $\epsilon$ or $x$)

II.3 Floquet Classical Master Equation

In the limit of $kT>\hbar\omega$, we can treat the nuclear motion classically, such that we will arrive at the Floquet classical master equation (FCME). The FCME is obtained by taking Wigner transform of FQME (Eqs. 20 and 21), similar to the steps taken in Ref. Dou, Nitzan, and Subotnik (2015):

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

$\tilde{f}(\Delta V)$ is given in Eq. 25. The Floquet CME is similar to the non-Floquet CME Dou, Nitzan, and Subotnik (2015), except that the Fermi function is modified with Floquet replicas. In the non-Floquet CME, $\gamma_{0\rightarrow 1}$ ( $\gamma_{1\rightarrow 0}$) is a real number, and can be interpret as the hopping rates. Such that the detailed balance is achieved. In the Floquet CME, however, the modified Fermi function $\tilde{f}(\Delta V)$ is time-dependent and complex valued, such that $\gamma_{0\rightarrow 1}$ ( $\gamma_{1\rightarrow 0}$) cannot be simply interpreted as hopping rates in trajectory based algorithm (see below). Still, in the limit of fast driving (large $\Omega$), we can take the time average of the modified Fermi function, such that $\gamma_{0\rightarrow 1}$ ( $\gamma_{1\rightarrow 0}$) will become real valued and can be interpreted as hopping rate. Below, we will use trajectories based surface hopping algorithm to solve the Floquet CME.

II.4 Floquet Surface Hopping

As stated above, the Floquet CME is similar to non-Flouqet CME in form, except the hopping rates are complexed valued. We now introduce 3 surface hopping algorithm to solve Floquet CME:

FSH: in FSH, we interpret the real positive part of $\gamma_{0\rightarrow 1}$ as our hopping rates and determine a hopping event from state $\ket{0}$ to state $\ket{1}$ based on $\xi<\Re(\gamma_{0\rightarrow 1})dt$, where $\xi$ is a random number between [math] and $1$. Similarly, a hopping event from state $\ket{1}$ to state $\ket{0}$ is determined by $\xi<\Re(\gamma_{1\rightarrow 0})dt$ 2. 2.

FaSH: in FaSH, we take the time average of $\tilde{f}(\Delta V)$ over a cycle that becomes $\overline{\tilde{f}}(\Delta V)$, such that the hopping rates are real valued and determine the hopping probability. 3. 3.

FaSH-density: in FaSH-density, we use the time-averaged $\overline{\tilde{f}}(\Delta V)$ as our hopping rates to propagate nuclear dynamics as same as FasH. In addition, we propagate electron density using $\dot{P_{0}}=-\gamma_{0\rightarrow 1}P_{0}+\gamma_{1\rightarrow 0}P_{1}$ and $\dot{P_{1}}=\gamma_{0\rightarrow 1}P_{0}(x,p)-\gamma_{1\rightarrow 0}P_{1}(x,p)$. Here $\gamma_{0\rightarrow 1}$ and $\gamma_{1\rightarrow 0}$ are non time-averaged, such that we can capture the oscillation in a cycle for electronic dynamics.

We initialize our nuclei in one well (the unoccupied electronic state) with a Boltzmann distribution. Notice that the initial nuclei temperature can be different from electronic temperature, such that we can simulate energy relaxation of the nuclear kinetic energy.

III Results

We now benchmark our trajectory based algorithms against FQME for electronic population as well as nuclear kinetic energy.

We first look at the the electronic population with periodic driving as a function of time. In Fig. 1a, we show the electronic population dynamics from FQME and FSH at different temperatures ($kT=0.25,0.5,1$). Here, we set nuclear vibration frequency as $\hbar\omega=0.3$. As expect, at high temperature ($kT>\hbar\omega$), FSH agrees with FQME very well; whereas at low temperature, FSH shows discrepancy as compared to FQME. Similar trend is observed for FaSH and FaQME as shown in Fig. 1b: in the high temperature limt, FaSH agrees with FaQME very well. Notice that, the time-averaged floquet results (FaSH or FaQME) do not show oscillations and reach to a steady state in the long time. The non time-averaged methods (FSH or FQME) do show the oscillations due to the non vanishing driving. In the long time, FSH or FQME reaches to a cycle limit instead of a stead state. The frequency of the oscillation in electronic population is equal to the frequency of the driving (we have set $\Omega=0.2$).

To further verify our methods, we show the electronic population and phonon relaxation dynamics as a function of time without any Floquet driving in Fig. 2. In such a case, we expect that the five algorithms (FQME, FSH, FaQME, FaSH, and FaSH-density) should all agree with each other at high temperature. Indeed, for the case of $kT>\hbar\omega$ ($kT=1$, $\hbar\omega=0.3$), the algorithms all agree with each other for both electronic population (Fig. 2a) as well as nuclear kinetic energy (Fig. 2b). Obviously, without any Floquet driving, both electronic population and nuclear kinetic energy reach to a steady state or equilibrium state. Indeed, the steady state kinetic energy is $E_{k}=1/2kT$ and electronic population is $N=f(\overline{E}_{d})$. Here $kT$ is the temperature of the electronic bath. Below, we will mainly look at the high temperature limit where the classical nuclei is valid.

In Fig. 3, we benchmark our algorithms for relatively small driving Amplitude with different driving frequencies. Here, the driving amplitude is $A=0.2$, which is comparable to the nuclear frequency ($\hbar\omega=0.3$). In such a case, the five methods all agree with each other regardless of the driving frequencies. That being said, the non time-averaged methods (FSH, FQME, and FaSH-density) do show small oscillations in electronic dynamics. The steady state in this case is very close to the equilibrium state with $E_{k}=1/2kT$ and $N=f(\overline{E}_{d})$.

We now turn to the intermediate driving amplitude, where we have set $A=1$. In such a case, the five methods start to show deviations as shown in Fig. 4. For smaller driving frequency (e.g. $\Omega=0.5$), the oscillation feature from FQME and FSH becomes stronger both in electronic populaiton and nuclear dynamics, whereas FaSH and FaQME miss the feature completely. FaSH-density shows oscillation in electronic dynamics and fail to reproduce the feature in nuclear dynamics. For slightly larger driving frequency ($\Omega=1$ or $10$), the oscillation feature becomes weaker. However, FSH starts to show deviations in the steady state at very large driving frequency ($\Omega=10$), where the kinetic energy is higher than the results from FQME (and others). Consequently, the electronic population from FSH is smaller than the results from other methods. This deviation under quick drivings is mainly derived from disregarding the negative part of the hopping rate in the FSH algorithm. This violates the detailed balance. As a result, the effective temperature of the system from FSH is higher than the true effective temperature from other methods. Note also that, the true effective temperature is slightly different from equilibrium temperature. This is due to the effects of the driving.

Finally, we show the case of very strong driving amplitude ($A=4$) in Fig. 5. In such a case, the oscillation feature becomes very strong, particularly for small driving amplitude. Again, FSH works very well in small driving frequency, but fails to produce the steady state electronic population and nuclear kinetic energy at large driving frequency. FaSH misses all oscillation feature at small driving frequency, but does reach to correct steady state. FaSH-density performs the best, working both in small driving frequency and large driving frequency.

IV Conclusions

In summary, we derived the Floquet quantum master equation to study the dynamics of a periodically driven system near a metal surface. In the high temperature limit, the FQME reduces to a FCME. We proposed three surface hopping algorithms (FSH, FaSH, and FaSH-density) to solve the FCME. In the limit of small driving frequency, FSH works very well, which captures the full dynamics for both electronic population and nuclear kinetic energy as benchmark against FQME. In the large driving frequency, FSH fails to produce the correct steady states. This is due to the fact that the hopping rates are complex valued and we throw out the negative hopping rates in FSH. The FaSH works well in the large driving frequency, but fails to reproduce the oscillation feature at small driving frequency. The FaSH-density performs the best, which reproduces the full electronic dynamics at small and large driving frequencies and reaches to the correct steady states for nuclear motion. We expect that our algorithms will be very useful to study non-adiabatic dynamics with periodic driving in open quantum systems.

Acknowledgements.

We thank Amikam Levy, Jacob Bätge, and Michael Thoss for useful and inspiring conversations. This material is based upon work supported by National Science Foundation of China (NSFC No. 22273075)