Quantum transport under AC drive from the leads : A Quantum Master Equation approach

TL;DR

This paper develops a comprehensive quantum master equation framework to analyze the dynamics of driven open quantum systems, enabling calculation of physical observables like current and revealing effects such as resonance.

Contribution

It introduces a general, transparent method based on Redfield equations for non-interacting particles under periodic driving, applicable to various geometries and particle types.

Findings

The theory accurately predicts current dynamics in driven quantum dot systems.

Time-dependent currents can surpass average currents, aiding experimental detection.

The approach is validated against exact numerical simulations.

Abstract

Evaluating the time-dependent dynamics of driven open quantum systems is relevant for a theoretical description of many systems, including molecular junctions, quantum dots, cavity-QED experiments, cold atoms experiments and more. Here, we formulate a rigorous microscopic theory of an out-of-equilibrium open quantum system of non-interacting particles on a lattice weakly coupled to multiple baths and driven by periodically varying thermodynamic parameters like temperature and chemical potential of the bath. The particles can be either bosonic or fermionic and the lattice can be of any dimension and geometry. Based on Redfield quantum master equation under Born-Markov approximation, we derive a complete set of linear differential equations for equal time two-point correlation functions from which various physical observables, for example, current, can be calculated. Various interesting physical effects, such as resonance, can be directly read-off from the equations. Thus, our theory is quite general gives quite transparent and easy-to-calculate results. We validate our theory by comparing with exact numerical simulations. We apply our method to a generic open quantum system, namely a double-quantum dot coupled to leads with modulating chemical potentials. Two most important experimentally relevant insights from this are : (i) time-dependent measurements of current for symmetric oscillating voltages (with zero instantaneous voltage bias) can point to the degree of asymmetry in the system, and (ii) under certain conditions, time-dependent currents can exceed time-averaged currents by several orders of magnitude, and can therefore be detected even when the average current is below the measurement threshold.

Figures

Equations

$$\mathcal{H}$$

$${\mathcal{H}}_S$$

$${\mathcal{H}}_{SB}$$

$$J_{\ell }(\omega )=2\pi \sum _r\mid {\kappa }_{\ell r}{\mid }^2\delta (\omega -{\Omega }_{\ell r})$$

$${\tau }_B\ll {\tau }_D\ll {\tau }_{expt}\ll T$$

$$c^{†}{H}^{(S)}c={\omega }^{(D)},$$

$${\hat{a}}_{\ell }=\sum _{\alpha =1}^N{c}_{\ell \alpha }{\hat{A}}_{\alpha }.$$

$$\displaystyle\frac{dC_{\alpha\nu}}{dt}=\Big{[}i\omega_{\alpha}C_{\alpha\nu}(t)-\varepsilon^{2}\sum_{\sigma=1}^{N}C_{\alpha\sigma}(t)v_{\nu\sigma}+(\alpha\leftrightarrow\nu)^{\dagger}\Big{]}$$

$$+{\epsilon }^2{Q}_{\alpha \nu }(t),$$

$$\mathrm{where},$$

$$v_{\alpha \nu }=f_{\alpha \nu }({\omega }_{\nu })-i{f}_{\alpha \nu }^{\Delta }({\omega }_{\nu })$$

$$Q_{\alpha \nu }(t)=[F_{\nu \alpha }({\omega }_{\alpha },t)-i{F}_{\nu \alpha }^{\Delta }({\omega }_{\alpha },t)+(\alpha \leftrightarrow \nu {)}^{\ast }]$$

$$f_{\alpha \nu }(\omega )=\sum _{\ell =1}^N{c}_{\ell \alpha }^{\ast }{c}_{\ell \nu }\frac{J_{\ell }(\omega )}{2},$$

$$F_{\alpha \nu }(\omega ,t)=\sum _{\ell =1}^N{c}_{\ell \alpha }^{\ast }{c}_{\ell \nu }\frac{J_{\ell }(\omega )n_{\ell }(\omega ,t)}{2}$$

$$\frac{dC}{dt}=-MC+{\epsilon }^{2}Q(t),$$

$$\mathrm{where}$$

$$M={\mathbb{I}}_N\otimes G^{\ast }+G\otimes {\mathbb{I}}_N,$$

$$G_{\alpha \nu }={\epsilon }^2{v}_{\alpha \nu }^{\ast }-i{\omega }_{\nu }{\delta }_{\alpha \nu }$$

$$C(t)=e^{-Mt}C(0)+{\epsilon }^2{\int }_0^{t}d{t}^{\prime }{e}^{-M(t-t^{\prime })}Q(t^{\prime })$$

$$Q(t+rT)=Q(t),r\in \mathbb{Z}$$

$$t=mT+\tau m\in \mathbb{Z}$$

$$\displaystyle C(mT+\tau)=\varepsilon^{2}\sum_{r=1}^{m}\Big{(}e^{-MrT}\Big{)}\int_{0}^{T}dt^{\prime}e^{-M(\tau-t^{\prime})}Q(t^{\prime})$$

$$+{\epsilon }^2{\int }_0^{\tau }d{t}^{\prime }{e}^{-M(\tau -t^{\prime })}Q(t^{\prime })+e^{-M(mT+\tau )}C(0)$$

$$C(mT+\tau )={\epsilon }^2(e^{MT}-{\mathbb{I}}_{N^2}{)}^{-1}{\int }_0^{T}d{t}^{\prime }{e}^{-M(\tau -t^{\prime })}Q(t^{\prime })$$

$$+{\epsilon }^2{\int }_0^{\tau }d{t}^{\prime }{e}^{-M(\tau -t^{\prime })}Q(t^{\prime })$$

$$C_{\alpha \nu }(t)=\sum _{p=-\infty }^{\infty }{C}_{\alpha \nu }^p{e}^{ip{\Omega }_{0}t},Q_{\alpha \nu }(t)=\sum _{p=-\infty }^{\infty }{Q}_{\alpha \nu }^p{e}^{ip{\Omega }_{0}t}$$

$$i({\omega }_{\alpha }-{\omega }_{\nu }-p{\Omega }_0)C_{\alpha \nu }^p+{\epsilon }^2{Q}_{\alpha \nu }^p$$

$$\displaystyle-\varepsilon^{2}\Big{[}\sum_{\sigma=1}^{N}C_{\alpha\sigma}^{p}v_{\nu\sigma}+C_{\sigma\nu}^{p}v_{\alpha\sigma}^{*}~{}\Big{]}=0$$

$$\displaystyle\frac{d}{dt}\Big{(}\sum_{\alpha=1}^{N}C_{\alpha\alpha}\Big{)}=\sum_{\ell}I_{B^{(\ell)}\rightarrow\ell}$$

$$\displaystyle I_{B^{(\ell)}\rightarrow\ell}=\varepsilon^{2}\Big{[}\sum_{\alpha=1}^{N}Q_{\alpha}^{(\ell)}(t)-\sum_{\sigma,\alpha=1}^{N}C_{\alpha\sigma}(t)\big{(}v_{\alpha\sigma}^{(\ell)}+v_{\sigma\alpha}^{(\ell)*}\big{)}\Big{]}$$

$$Q_{\alpha }^{(\ell )}(t)=|c_{\alpha \ell }{|}^2{\mathfrak{J}}_{\ell }({\omega }_{\alpha })n_{\ell }({\omega }_{\alpha })$$

$$\displaystyle v_{\alpha\sigma}^{(\ell)}=c_{\alpha\ell}^{*}c_{\sigma\ell}\Big{(}\frac{\mathfrak{J}_{\ell}(\omega_{\sigma})}{2}-i\mathcal{P}\int\frac{d\omega}{2\pi}\frac{\mathfrak{J}_{\ell}(\omega)}{\omega-\omega_{\sigma}}\Big{)}$$

$$\displaystyle I_{\ell\rightarrow\ell+1}=2\hskip 2.0ptIm\Big{(}H^{(s)}_{\ell\hskip 2.0pt\ell+1}\hat{a}_{\ell}^{\dagger}\hat{a}_{\ell+1}\Big{)}$$

$$\frac{{\omega }_{\alpha }-{\omega }_{\nu }}{p}-\frac{{\omega }_{\alpha }-{\omega }_{\nu }}{p+1}=\frac{{\omega }_{\alpha }-{\omega }_{\nu }}{p(p+1)}\sim \frac{1}{p^2}$$

$${\hat{\mathcal{H}}}_S$$

$${\hat{\mathcal{H}}}_B^{(\ell )}$$

$${\hat{\mathcal{H}}}_{SB}$$

$${\mathfrak{J}}_{\ell }(\omega )={\Gamma }_{\ell }\sqrt{1-{\left(\frac{\omega }{2t_B}\right)}^2},{\Gamma }_{\ell }=\frac{2{\gamma }_{\ell }^2}{t_B}.$$

$$n_{\ell }(\omega ,t)=[e^{\beta (\omega -{\mu }_{\ell }(t))}+1{]}^{-1},{\mu }_{\ell }(t+rT)={\mu }_{\ell }(t)\forall r\in \mathbb{Z}$$

Full text

.1 Introduction

The dynamics of a quantum system coupled to an external environment (so-called “open” quantum system) are a result of the intricate interplay between the coherent evolution of the quantum degrees of freedom, the structure of the environment and the form and strength of the system-environment coupling. From quantum cavities Aspelmeyer et al. (2014); Ludwig and Marquardt (2013); Brunelli et al. (2016) and superconducting qubits Fink et al. (2008); Hofheinz et al. (2008); Bishop et al. (2009); Schmidt and Koch (2013); Hoffman et al. (2011); Houck et al. (2012) , through quantum dots van der Wiel et al. (2002); Keller et al. (2016); Hartmann et al. (2015); Venturelli et al. (2013); Gruss et al. (2016); Ajisaka and Barra (2013), molecular junctions Peskin (2016, 2010); Thingna et al. (2016); Segal and Agarwalla (2016); Ajisaka et al. (2015); Wang et al. (2015); Zelovich et al. (2015); Elenewski et al. (2017); Esposito and Galperin (2010); Hod et al. (2016); Rudge and Kosov (2016); Dubi (2013) and cold atoms Chien et al. (2015); Baumann et al. (2010); Krinner et al. (2015); Kulkarni et al. (2013) , to excitons traveling in photosynthetic complexes Mohseni et al. (2014); Lambert et al. (2013); Fleming et al. (2011); Scholes and Fleming (2005), open quantum systems show dynamics which can be far richer and more surprising than their coherent (environment-free) counterparts.

Recent ideas of designing a quantum state by engineering a specific environment Kraus et al. (2008); Barreiro et al. (2010); Kienzler et al. (2015); Verstraete et al. (2009); Diehl et al. (2008); Bardyn et al. (2013); Vorberg et al. (2013) or by a periodic modulation of the open system’s Hamiltonian Volokitin et al. (2016); Hartmann et al. (2016); Peskin (2016) open a path to new forms of control over quantum systems. Combining these two concepts of time-periodic modulations and environment (bath) engineering, here we study the dynamics of open quantum systems where the environment thermodynamic parameters are periodically modulated. Even though there exists formally exact methods of treating such set-ups Jauho et al. (1994); Zhu et al. (2005); Maciejko et al. (2006); Tu and Zhang (2008); Jinshuang et al. (2010); Ridley et al. (2015, 2016a, 2016b); Yang and Zhang (2016), they are generally quite difficult to use, bein computationally demanding. To simplify calculations, often the adiabatic approximation is used Perfetto et al. (2015); Yadalam and Harbola (2016); Perfetto et al. (2016). The adiabatic approximation assumes the environment parameters are modulated very slowly such that the system is always at a non-equilibrium steady state with the instantaneous environment. This approximation fails to capture dynamical effects which arise due to the time-periodic modulations. In this work, we therefore seek to go beyond the adiabatic approximation.

For weak system-bath coupling, a combination of Quantum Master Equation methods (typically of the Lindblad form) with Floquet theory is often used to go beyond the adiabatic approximation Kohler et al. (2005); Ponomarev et al. (2009); Wu and Timm (2010); Hausinger and Grifoni (2010); Stace et al. (2013); Szczygielski et al. (2013). However, Floquet theory has the drawback of converting a finite-dimensional problem to an infinite dimensional one. Consequently, most interesting results from Floquet theory are often obtained as a perturbation expansion in terms of inverse drive frequency, giving results primarily for high frequency driving.

Moreover, from calculations in the absence of time-periodic modulations, the commonly used phenomenological Lindblad Quantum Master Equations are known to have several drawbacks, especially in an out-of-equilibrium situation Purkayastha et al. (2016); Wichterich et al. (2007); Bauer (2011); Levy and Kosloff (2014). It has been recently shown that microscopically derived Redfield Quantum Master Equation (RQME) under Born-Markov approximation can be used to overcome the drawbacks of phenomenological Lindblad equations and to obtain correct results up-to leading order in system-bath coupling Purkayastha et al. (2016). To our knowledge, there has been no work in extending the RQME to the case when thermodynamic parameters of the baths are periodically modulated.

In this work, we derive a Born-Markov approximated Redfield theory for a quantum system consisting of non-interacting particles (bosons or fermions) on a lattice of arbitrary dimension and geometry, weakly coupled to multiple periodically-modulated baths. The baths are modeled by non-interacting particles (bosons or fermions) with infinite degrees of freedom. For a lattice of $N$ sites, this approach leads to a closed system of $N^{2}$ linear differential equations which can be easily solved numerically in time domain, thereby avoiding the drawback of infinite dimensional theory like Floquet theory. We apply our formalism to a generic example of an open quantum system, namely a double quantum-dot coupled to two fermionic reservoirs, characterized by periodically modulated chemical potentials. The current through the system depends on both the chemical potential modulation properties and the chemical potential difference (i.e. voltage bias) between the electrodes. We show that even when there is no bias between the electrodes the dynamics of the system are non-trivial and reveal the internal structure of the junction (i.e. the internal spectrum of the dot and the asymmetry in the junction). When there is a voltage bias between the electrodes, the modulation induces unusual hysteresis behavior. These predictions can in principle be tested experimentally, demonstrating the applicability of our formalism to realistic open quantum systems.

.2 Set-up and protocol

Consider a general non-interacting tight-binding Hamiltonians for the whole system+bath set-up :

[TABLE]

$H^{(s)}$ is a Hermitian matrix and $\hat{a}_{\ell}$ correspond to fermionic (bosonic) annihilation operators defined respectively on $\ell$th lattice point of system and $\hat{B}_{\ell r}$ to those of baths attached to the $\ell$th point. The baths have infinite degrees of freedom. $\varepsilon$ is a dimensionless parameter that controls system bath coupling, so that $\{\kappa_{\ell r}\}$ have dimensions of energy. Here each lattice point is connected to a bath. By setting system-bath coupling to zero at an arbitrary number of lattice points, any geometry of system-bath connections can be achieved. We also introduce the bath spectral functions:

[TABLE]

Note that the full system+bath Hamiltonian (1) is time-independent. In the conventional dc bias case, for weak system-bath coupling, the baths are always assumed to be in thermal equilibrium with their own corresponding temperature and chemical potential. In other words, $Tr_{B}(\hat{B}_{\ell r}^{\dagger}\hat{B}_{\ell r})$, where $Tr_{B}(..)$ implies trace taken over bath degrees of freedom, is given by the corresponding fermi (bose) distribution \Big{[}e^{-\beta_{\ell}(\Omega_{\ell r}-\mu_{\ell})}\pm 1\Big{]}^{-1}. Here, we are interested in the case when the inverse temperature $\beta_{\ell}$ and chemical potential $\mu_{\ell}$ are periodic functions of time. At this stage, one may want to model such a set-up by deriving a quantum master equation for the density matrix of the system for the time independent case and making the fermi (bose) distributions time periodic ‘by hand’. However, the regime of validity of such an approach will not be clear. So, instead, in the following, we assume a protocol, and make physical assumptions so that such an equation may be reached. This shows how and under what conditions such a situation can arise physically.

We now describe the protocol. In the following, $\chi$ is the full density matrix of system+bath, $\rho_{B}$ is the density matrix of the bath, and $\rho$ is density matrix of the system obtained by tracing $\chi$ over bath variables. Let us, for the time being, assume one bath. The protocol can readily be generalized to multiple baths.

The protocol for a single bath is as follows :

a) At time $t=0$, $\chi=\rho\otimes\rho_{B}(0)$, with $\rho=\rho_{0}(0)$, which is some arbitrary initial system state, and $\rho_{B}(0)={exp[-\beta(0)(\mathcal{H}_{B}-\mu(0)\mathcal{N})]}/{Z(\tau_{S})}$. That is, the inital state is a product state of an arbitrary system state and a thermal bath state.

b) We evolve the system for a time $\tau_{D}$. After time $\tau_{D}$, the system reaches the state $\rho_{0}(\tau_{D})$. Note that during this time, the temperatures and chemical potentials of the bath has not changed.

c) At $t=\tau_{D}$, we restart the entire system+bath setup with the initial state $\chi=\rho\otimes\rho_{B}(\tau_{D})$, with $\rho=\rho_{1}(0)=\rho_{0}(\tau_{D})$ and $\rho_{B}(\tau_{D})={exp[-\beta(\tau_{D})(\mathcal{H}_{B}-\mu(\tau_{D})\mathcal{N})]}/{Z(\tau_{D})}$. That is, at time $\tau_{D}$, the bath is changed into the thermal state with new inverse temperature and chemical potential $\beta(\tau_{D}),\mu(\tau_{D})$, and the full system+bath state is again taken as the product state. This step implicitly assumes weak-system bath coupling. This is because, under weak system-bath coupling, we can assume that the bath is hardly affected by the system, and to leading order in system-bath coupling, the full density matrix is in product form.

d) We again let the system density matrix evolve under this new bath for time $\tau_{D}$ starting from the new initial state. Again, after time $\tau_{D}$, we restart the entire set-up with initial state $\chi=\rho(0)\otimes\rho_{B}(2\tau_{D})$, with $\rho=\rho_{2}(0)=\rho_{1}(\tau_{D})$ and $\rho_{B}(2\tau_{D})={exp[-\beta(2\tau_{D})(\mathcal{H}_{B}-\mu(2\tau_{D})\mathcal{N})]}/{Z(2\tau_{D})}$. This protocol continues for a very long time and we will be mostly interested in the long time dynamics of this process.

In the limit of very small $\tau_{D}$, the above protocol gives a nearly continuous evolution of temperatures and chemical potentials of the baths. In particular, if the temperature and chemical potential vary periodically with a period $T$, and $\tau_{D}\ll T$, then the above protocol describes dynamics of a system evolving under a continuous periodic drive from the temperature and chemical potential of the bath.

To make analytical progress with this problem, we need to make one more assumption, the Markov assumption. Let $\tau_{B}$ be the characteristic time-scale of relaxation of the effects of system-bath coupling on the bath. Then, we assume that $\tau_{B}\ll\tau_{D}$. It is only under this condition that microscopically derived Markovian Quantum Master Equations can be applied. So, if $\tau_{expt}$ be experimentally the smallest time scale, most of our following discussion holds if the following condition on time-scales is maintained :

[TABLE]

It can be shown that $\tau_{B}$ depends on the temperature of the bath, $\tau_{B}\sim\beta$ (see Appendix B). Hence, our analytical discussion below will not be valid at extremely low temperatures. It is only valid when temperature is large enough so that above condition on time scales can be satisfied.

For multiple baths, the above protocol is followed for each bath. $\tau_{B}$ is then taken as the largest of the relaxation times of the baths, so that $\tau_{D}$ is much larger than relaxation time scales of all baths. When the system is driven by multiple baths according to our protocol, we call it an ac drive process. Correspondingly, the usual case of having a time independent temperature or chemical potential difference between baths is called a dc bias process.

The above protocol, along with the condition (Eq. 4) on time scales, breaks down the ac drive process into steps of time independent processes where Born-Markov approximation can be applied. For such time independent processes, the Born-Makov approximated Redfield Qunatum Master Equation (RQME), as well as the evolution equation for two point correlation functions from the RQME have been derived for Hamiltonian (1) recently Purkayastha et al. (2016) (see Appendix A). Going through the above protocol with the RQME only has the effect of making the fermi (bose) distribution functions time-dependent with the instantaneous temperatures and chemical potentials (see Appendix C). The final result, as desired, will be same as that obtained by ‘putting by hand’ time dependent fermi of bose distribution functions in the RQME derived without such time dependence. However, it is important to note that, the description via such an equation would not be possible if each step of our protocol would not satisfy the Markov condition $\tau_{B}\ll\tau_{D}$. Hence, only when system-bath coupling is weak and the condition on time scales in Eq. 4 is satisfied, can our set-up and protocol be described by such an equation. This crucial fact would not be clear in a ‘putting by hand’ approach. Also, for small systems, the protocol can easily be exactly (i.e, without Markov approx) simulated with finite but large baths, and hence allows for numerical validation of results.

.3 Redfield equation for the correlation matrix (single particle density matrix)

Since the system is quadratic, various physical observables like current and occupation can be directly calculated from the two-point correlation functions. Using the results in Ref. Peskin (2016) (which are also re-derived in Appendix), we can readily write down the evolution equation for two-point correlation functions for our set-up. For this, it is convenient to go to the eigenbasis of the system Hamiltonian. Let $c$ be the unitary matrix which diagonalizes $H^{(S)}$, i.e.,

[TABLE]

where $c^{\dagger}c=I$ and $\omega^{(D)}$ is a diagonal matrix with elements $\omega_{\nu}$. Then we also define new operators $\{A_{\alpha}\}$ through the transformation

[TABLE]

Thus $\hat{A}_{\alpha}$ is the annihilation operator for the $\alpha$th eigen-mode with energy $\omega_{\alpha}$.

The evolution equation for the equal time two point correlation functions C_{\alpha\nu}(t)=\langle\hat{A}_{\alpha}^{\dagger}(t)\hat{A}_{\nu}(t)\rangle=Tr\big{(}\rho\hat{A}_{\alpha}^{\dagger}(t)\hat{A}_{\nu}(t)\big{)} can be derived by following the derivation of the Redfield equation for the density matrix $\rho$ Breuer and Petruccione (2007); Carmichael (2002); Nitzan (2006) and substituting it back to the definition of the correlation functions Purkayastha et al. (2016)(see Appendix A). The protocol, along with Born-Markov approximation, only makes the fermi (bose) distributions time-dependent (see Appendix B). The resulting equation is

[TABLE]

with $f_{\alpha\nu}^{\Delta}(\omega)=\mathcal{P}\int\frac{}{}\frac{d\omega^{\prime}f_{\alpha\nu}(\omega^{\prime})}{\pi(\omega^{\prime}-\omega)},\hskip 2.0ptF_{\alpha\nu}^{\Delta}(\omega,t)=~{}\mathcal{P}\int\frac{}{}\frac{d\omega^{\prime}F_{\alpha\nu}(\omega^{\prime},t)}{\pi(\omega^{\prime}-\omega)}$, where $\mathcal{P}$ denotes principal value, n_{\ell}(\omega,t)~{}=~{}\Big{[}e^{-\beta_{\ell}(t)(\omega-\mu_{\ell}(t))}\pm 1\Big{]}^{-1} is the instantaneous fermi or bose distribution function and $J_{\ell}(\omega)$ is the spectral function (Eq. 3) of bath (lead) attached to $\ell$th site of the system. Weak system-bath coupling requires $\varepsilon^{2}J_{\ell}(\omega)\ll\{\omega_{\alpha}\}$.

Eq. .3 gives a closed set of $N^{2}$ linear differential equations. The matrix with elements $C_{\alpha\nu}$ is sometimes also called the single particle density matrix of the system. Writing $C$ and $Q$ as vectors, this equation can be cast in the form

[TABLE]

$\mathbb{I}_{N}$ is the $N$ dimensional identity matrix and $\otimes$ denotes Kronecker product. Eq. 8 has the formal solution

[TABLE]

Note that the matrix $M$ has no time dependence and hence is the same matrix as would appear in the dc problem. If under dc bias the system reaches a steady state after a long time, the real part of eigenvalues of the matrix $M$ has to be positive. We will assume this henceforth.

We wish to look at the long time properties of Eq. 9. Let our ac drive be periodic with a time period $T$. Then, for integer $r$,

[TABLE]

We break up the observation time $t$ into steps of $T$, so that, for integer $m$,

[TABLE]

Then Eq. 9 can be written as

[TABLE]

Now, since we have assumed that the real part of all the eigenvalues of $M$ are positive, we can perform the sum in above equation. Also, we are interested in the long time limit, $m\gg 1$. Hence we have

[TABLE]

Note that, in the long time limit, the RHS is independent of $m$, and is also independent of the initial condition. This means that in the long time limit, the correlation functions settle down in a periodic state, with period same as the ac drive. This is consistent with Floquet theory. However, note that, Floquet theory was not used explicitly to derive this result. Eq. .3 is the central result of the manuscript.

.4 Resonances and Currents

Now, let us go back to Eq. .3. Since in the long time limit the periods of $C(t)$ and $Q(t)$ are same, we can perform a Fourier series expansion :

[TABLE]

where $\Omega_{0}=\frac{2\pi}{T}$. Substituting in Eq. .3, we obtain the following equation for each Fourier mode,

[TABLE]

This immediately gives us two important results. First, we note that for the time-average correlation functions, $\overline{C}_{\alpha\nu}(t)=\frac{1}{T}\int_{t}^{t+T}dt^{\prime}C_{\alpha\nu}(t^{\prime})=C_{\alpha\nu}^{0}$, and similarly for $Q_{\alpha\nu}$. With this we find that the equation averaged over one time period gives exactly the steady state equation for a dc bias given by the time-period averaged bose or fermi distribution functions $\overline{n_{\ell}(\omega,t)}$.

Second, we see that when $\omega_{\alpha}-\omega_{\nu}=p\Omega_{0}$, Eq. .4 becomes independent of system-bath coupling $\varepsilon^{2}$. This is the phenomenon of resonance. At resonance, the leading term of system correlation functions do not depend on system bath coupling. So all system properties, like current between two sites inside the system, which were otherwise proportional to $\varepsilon^{2}$ in the leading term, become larger by orders of magnitude.

Even though system properties become large at resonance, the current from any of the baths still remain small. To see this, we need the expression for current from the baths in terms of the correlation functions. The expression for particle current is obtained from the continuity equation for conservation of particles

[TABLE]

$I_{B^{(\ell)}\rightarrow\ell}$ is the particle current from bath attached to the $\ell$th site. Under dc bias, in steady state, the LHS of above equation is zero, and hence the currents from the baths are equal. However, under ac drive, even at long time, the LHS is not zero, and thus, instantaneous currents from the baths can be different. The summation on the right runs over all lattice sites connected to bath. The expressions for the currents from the baths obtained from above continuity equations and Eq. .3 are

[TABLE]

where

[TABLE]

These expressions show that the currents from the baths are explicitly proportional to $\varepsilon^{2}$. So even when the system correlation functions are independent of $\varepsilon^{2}$ in the leading order, the currents from the baths are still $O(\varepsilon^{2})$. However, this is not true for particle current between two sites which reside within the system. The current between $\ell$th and $\ell+1$th lattice sites of the system is given by

[TABLE]

This current is not explicitly proportional to $\varepsilon^{2}$. At resonance, it is independent of $\varepsilon^{2}$ in the leading order. Thus, at resonance, the current between two adjacent lattice points in the system can be much larger than the current from the baths. On the other hand, the time-period averaged current, which corresponds to steady state of a dc bias, is same inside the system as from the baths.

The difference in frequency between adjacent resonances, which is given by

[TABLE]

decreases as $1/p^{2}$. Hence small driving frequencies are always close to resonance with one of the higher modes (large $p$) of the steady state oscillations. So, for small driving frequencies, there is not much difference between the resonance and off-resonance condition. The values of the correlation functions, however, may not be so large as the first resonance ($p=1$). This is because, the weight of the driving signal at higher modes may decrease.

It is also interesting to note that none of the expressions for currents can be reduced to the form of difference between fermi (bose) distributions of the various baths. Therefore, even when all the baths are driven by the exactly same time dependent temperature or chemical potential (symmetric ac drive), so that there is no instantaneous temperature or chemical potential difference, there can be an instantaneous current, both between the system and the baths and inside the system. At resonance, the internal currents may have a large amplitude (compared to the system-bath currents). The time-period averaged current is, of course, zero in this case.

It is important to state that we have not made the ubiquitous adiabatic approximation $T\gg t_{steady}$, where $t_{steady}$ is the time to reach steady state, which corresponds to the smallest real part of eigenvalues of the matrix $M$ in Eq. 8. If adiabatic approximation were made, then the expressions for the correlation functions at any time would be given by the dc-bias steady state results with the fermi (bose) distributions given by the instantaneous temperatures and chemical potentials. In that case, the expressions for currents would have reduced to the form of difference between fermi (bose) distributions of the various baths, and no instantaneous current would have been seen in the symmetric ac drive case.

Under dc bias, the steady state can be quite easily obtained for non-interacting open systems described by Hamiltonian of the form Eq. 1 by exact methods (i.e, not under Born-Markov approximation). Transient dynamics of approach to steady state, on the other hand, is generally quite difficult to calculate from such methods. From our protocol, it can be seen that when, $T\lesssim t_{steady}$ (i.e, when adiabatic approximation is not valid), it is the transient behaviour of approach to steady state for dc bias that becomes important to describe the ac driven case. Therefore, such exact methods are difficult to use in our ac drive set-up. However, as has been recently established, RQME gives the correct time dynamics of approach to steady state under dc bias of two-point correlation functions for our set-up as long as Born-Markov approximation is satisfied Purkayastha et al. (2016). This fact has allowed us to use RQME for ac drive case to go beyond the adiabatic approximation easily.

Recently, another method based on a quantum master equation called the Discrete Lioville von-Neumann (DLvN) method has been formulated Zelovich et al. (2014, 2017); Chen et al. (2014), which has been shown to give accurate results for time-dependent behaviour Chen et al. (2014). However, in this method one models the leads via finite but large size Hamiltonians, and needs to solve explicitly for the leads along with the system. As a result, if there are $N_{L}$ sites in left lead, $N_{R}$ sites in the right lead and $N$ sites in the system, then one needs to solve for a closed set of $(N_{L}+N_{R}+N)^{2}$ linear differential equations. In comparison, using our RQME method, the leads are considered implicitly in the spectral functions (Eq. 3) and one needs to solve only for a closed set of $N^{2}$ linear differential equations (see Eqs. .3 and 8).

RQME has the drawback of not being completely positive because of not being of the Linblad form. As a result, under certain initial conditions, it may lead to a non-physical state at short times. In fact, the rigorous microscopic derivation of the dc bias Redfield equation for quadratic baths requires the observation time $t$ to satisfy $t\gg\tau_{B}$ (see Appendix A). Thus, it is only in this regime of observation times that the results should be valid. As a consequence, one expects that, problems regarding non-positivity show up only at short times ($t\leq\tau_{B}$), where the Redfield equation itself is not justified Eastham et al. (2016). Based on this, it is expected that, in our ac drive set-up, when the condition on time scales in Eq. 4 is satisfied, problems regarding positivity will not arise. This, however, is not a completely rigorous statement, and requires further investigation. In Ref Purkayastha et al. (2016), the long time behaviour of the correlation functions for the dc bias case of our set-up have been shown to be given correctly by the Redfield equation. Since our protocol breaks down the ac drive process into a series of such dc drive steps, this gives further support regarding validity of our results. But, as mentioned before, since $\tau_{B}$ is inversely proportional to temperature (see Appendix B), our theory may not be valid at extremely low temperatures.

.5 A detailed example: the $N=2$ case

All the above results are quite general and hold for non-interacting system of bosons or fermions in lattice of any dimension and geometry. In the following, to validate the theory, as well as to better understand the physics of such ac drive, we apply this theory to a simple model of two fermionic lattice sites, each connected to their own baths. This generic model can be used to describe physical systems such as driven double quantum dotsStehlik et al. (2016) or single-molecule junctions with two molecular moieties (e.g. biphenyl-dithiol molecular junctions Venkataraman et al. (2006)).

We consider the following specific two-site system coupled to baths which are one-dimensional chains:

[TABLE]

where the operators are all fermionic and $\hat{b}_{s}^{\ell}$ is the annihilation operator of the $s$th bath site of the $\ell$th bath. The eigenmodes of the system are given by $\hat{A}_{1}={(\hat{a}_{1}-\hat{a}_{2})}/{\sqrt{2}}$, $\hat{A}_{2}={(\hat{a}_{1}+\hat{a}_{2})}/{\sqrt{2}}$ with eigenvalues $\omega_{1}=\omega_{0}-g$, $\omega_{2}=\omega_{0}+g$. The bath spectral functions, defined in Eq. 3, can be obtained explicitly by going to eigenmodes of the baths and are given by Purkayastha et al. (2016); Peskin (2016) (see Appendix D)

[TABLE]

$\omega_{0}\gg(\varepsilon^{2}/t_{B})$ so that QME can be applied while the parameter $g$ can be varied freely. The two baths are taken at the same temperature. Inverse temperature $\beta$ is taken to be constant and the ac drive is given by periodic chemical potential :

[TABLE]

We choose the drive frequency as

[TABLE]

If $p$ is integer, the system is at one of the resonances. If $p$ is not an integer, it is away from resonance.

We look at the particle current in the ac driven steady state. The expressions for the currents are given in Eq. 16, 19. For the $N=2$ case, the expression for $I_{1\rightarrow 2}$ simplifies to $I_{1\rightarrow 2}=2g\hskip 2.0ptIm(C_{12})$. Two different ac drives are considered:

(i) symmetric ac drive or zero voltage drive :

[TABLE]

(ii) asymmetric ac drive :

[TABLE]

.5.1 Comparison between analytic formula and exact numerics

Before proceeding to elucidate the physics dominating the driven system, we wish to validate the analytic derivation. For this we compare results obtained from Eq. 9 with that obtained from full numerical simulation of our protocol. For small systems like the two-site case considered here, the protocol for our ac drive set-up can be simulated exactly with finite but large baths. For each time-independent step of our protocol, we choose a bath of finite size with bath correlations satisfying fermi distributions and evolve the full system+bath Hamiltonian $\hat{\mathcal{H}}$ using unitary Hamiltonian dynamics. Let us collectively denote by “$d$” a column vector with all annihilation operators of both system and baths. The full Hamiltonian can be written as $\hat{\mathcal{H}}=\sum_{i,j}H_{ij}d^{\dagger}_{i}d_{j}$ where $i$ now refers to either system or bath sites. If $D=\langle dd^{\dagger}\rangle$ denotes the full correlation matrix of system and baths, its time evolution is given by $D(t)=e^{iHt}De^{-iHt}$. The bath correlations are then changed according to the protocol, and the process is repeated. Various observables like current inside the system calculated using this exact numerical simulation can be compared with that obtained from Eq. 9, thus providing a way to validate our theory. Note that the numerical simulation does not take into consideration the ‘Markov’ condition on time scales given in Eq. 4. It holds even when Eq. 4 is not respected. Thus it allows us to check the validity of the crucial assumption on time scales required for our analytical treatment.

Fig. 1 shows the results for current inside the system $I_{1\rightarrow 2}^{s}$ and $I_{1\rightarrow 2}^{a}$ for cases (i) Eq. 25 and (ii) Eq. 26 respectively (superscripts in $I$ stand for symmetric and asymmetric currents), as obtained from exact numerics as well as from our theory. Eq. .3 has been used to obtain the long time result, while Eq. 9 has been used to get the result at all times, showing approach to the long time dynamics. Numerical simulation has been done with baths of size $256$ sites, which are large enough to have negligible finite-size effects. The near perfect agreement with exact numerical simulations validate our theory. For this plot, the drive frequency is chosen to be $\Omega_{0}=\omega_{1}-\omega_{2}=2g$, so that the set-up is at the first resonance. The near perfect match occurs for other frequencies also as long as Eq. 4 is satisfied. For our choice of parameters $\beta\sim 0.1$ satisfies Eq. 4. The initial condition for plots shown in the figure corresponds to no particle in the system, but the agreement with numerical simulations has been checked for other initial conditions (like randomly chosen initial values of correlation functions etc.) also.

Note that numerical validation was only possible owing to the small size of the system, which allowed for using reasonably sized finite baths. For larger system sizes, much larger baths will be required and the set-up will not be amenable to numerical simulation. However, the theory can be easily used for much larger systems connected to infinite baths.

Having validated the theory, we now look at the physics of the long time dynamics for both the cases.

.5.2 Symmetric - or zero voltage - ac drive

First, we look at long time dynamics of the symmetric drive (Eq. 25). In this case, both chemical potentials are the exact same sinusoidal so that there is no instantaneous voltage difference, i.e, $V(t)=\mu_{1}(t)-\mu_{2}(t)=0$. Even though there is no voltage difference, as discussed before, because of being a driven system, there can still be an instantaneous current. Moreover, because of presence of displacement current, the current from the left bath, current in the system, current from the right bath can be different.

Panel (a) of Fig. 2 shows currents from the left (right) bath to site 1(2), $I^{s}_{B^{(1)}\rightarrow 2}~{}(I^{s}_{B^{(2)}\rightarrow 2})$, and current inside the system over one time period for the case of first resonance ($p=1$ in Eq. 24). As expected from our discussion of resonance, the current inside the system is much larger than the currents from the baths. A physical explanation for the non-zero instantaneous current for zero voltage drive can be given as follows.

Let $V_{0}>\omega_{2}$ (see Eq. 25). As the chemical potential varies, the Fermi energy of the particles in the bath varies. When it exceeds $\omega_{2}$, particles flow into the system from the baths. When the Fermi energy of the particles in the bath is smaller than $\omega_{1}$, particles flow out of the system into the baths. If the time to reach steady state is larger than the time period of the drive, this transient behavior is observed, which leads to the instantaneous current.

It is clear from this argument, that the instantaneous current increases with increase in $V_{0}$. Particularly, if $V_{0}<\omega_{1}$, there will be a small instantaneous current. This behaviour is shown in panel (b) of Fig. 2 which shows the maximum instantaneous currents from the baths and inside the system as a function of $V_{0}$. The current increases continuously with $V_{0}$ and finally saturates. The saturation occurs because there are only two eigen-energy levels of the system. Actually, because of the fermionic nature of the set-up, one would expect steps or kinks at the positions $V_{0}=\omega_{1}$ and $V_{0}=\omega_{2}$. However, such behaviour is not observed because temperature is not low enough. On the other hand, it is important to note that the time-period averaged current is zero always because it is proportional to the difference between time-period averaged fermi distributions of the two baths.

We note that the internal current, while strictly speaking cannot be measured directly, nevertheless may have physical consequences. Specifically, such internal current may lead to local heating of the junction (due to Joule heating). This can lead for instance, to a breakdown of the system (in molecular junctions) or to heating-induced observable changes in current (for a double quantum dot).

Panel (c) of Fig. 2 shows a more interesting effect. The plot shows the maximum instantaneous current from the baths and inside the system as a function of asymmetry $\Gamma_{1}/\Gamma_{2}$ of the system-bath coupling (note that the chemical potentials are still symmetric). The set-up is still at first resonance. The maximum current inside the system decreases with decrease in asymmetry, and vanishes for the fully symmetric junction, $\Gamma_{1}/\Gamma_{2}=1$. At the same time, in the symmetric point the maximum currents from left and right baths are equal. This can be physically explained as follows. If the baths are identically coupled as well as identically driven, the particles come into the system at exactly same rate leading to same current from both baths. Since the particles are fermions, if one particle occupy each site, there can be no current in between the two sites due to Pauli exclusion principle. This leads to the fact that if rate of inflow of particles from both baths is same, there is no current between the two sites. Thus, current between the two sites of the system comes from a mismatch between rate of inflow of particles from the baths. Therefore, current inside the system increases with increase in asymmetry of system-bath coupling.

Panel (d) of Fig. 2 shows plots of the maximum current inside the system for asymmetric system-bath coupling with $p$ (as defined in Eq. 24) for two different temperatures. When $p$ is equal to an integer, the set-up is at resonance. Thus, from our previous discussion, the maximum current inside the system should show peaks at integer values of $p$. The first resonance peak is strong. However, lower frequency resonance peaks are much weaker. This is because the driving signal itself has a small contribution from such modes. Moreover, we see that higher temperature washes out the lower frequency resonances.

We point out that asymmetry in molecular junctions may be a key ingredient in deciphering its electronic transport properties and in designing single-molecule devices (see, e.g., Warner et al. (2015); Zhang et al. (2015); Wang et al. (2014)). Our results demonstrate that measuring the time-dependent current (for zero voltage bias but time-dependent voltages) is a direct way to measure asymmetry in molecular junctions, which can serve (in parallel to usual transport measurements) for junction characterization.

.5.3 Asymmetric ac drive

Having discussed the physics of long time dynamics of the symmetric ac drive, we move to the case of the asymmetric ac drive. Consider the situation where the left chemical potential and right chemical potentials vary sinusoidally out of phase with same frequency but have different amplitude in general.

Unlike the symmetric ac drive case, here the voltage across the system $V=\mu_{1}-\mu_{2}$ is non-zero. For $V_{1}\neq V_{2}$. The difference between time-period averaged Fermi distributions of the two baths is not zero in this case. Hence, there is a net time-period averaged current through the system. Fig. 3 shows plots of currents from the baths and current in the system at first resonance. Panel (a) shows current vs voltage ($V=\mu_{1}-\mu_{2}$) curve over one time period. We observe an interesting hysteresis behavior. Also, because of resonance, the current inside the system is much greater than current from the baths. This effect would not have been seen in the adiabatic limit. Panel (b) shows the time-period averaged current as a function of $V_{1}-V_{2}$. The current increases, and finally saturates (similar to $I_{max}$ vs $V_{0}$ curve for the symmetric drive case). Note that the time-averaged current is much smaller than the maximal instantaneous currents. This implies that perhaps time-dependent signals can be measured even when the average currents are small and below the noise level.

The hysteresis in I-V curve for asymmetric system bath coupling can potentially have device applications. This behavior depends intricately on the phase difference between the two drives, the amplitude difference of the two drives, as well as on the asymmetry of system-bath coupling. These dependences are quite complicated and a detailed investigation of the hysteresis behavior will be taken up in a future work.

The two site set-up described above is experimentally realizable using quantum dots. However, the bottleneck experimental parameter for observation of most of the above effects is frequency of the drive $\Omega_{0}$. Our most interesting results are close to resonance $\Omega_{0}=2g$. The maximum frequency of ac drive currently experimentally realizable is $\sim 10GHz$ Trasobares et al. (2016). This means to observe the above effects $g\sim 0.01meV$. The baths must have much wider bandwidths than system energy scales, which can be easily arranged experimentally. As evidenced by Eq. 22, having wide bandwidths will automatically realize the weak system-bath coupling. Temperature $\sim 1mK$ (which is equivalent to $\beta\sim 0.1$) will be consistent with the Markov approximation. These experimental parameters may be challenging but not impossible. Our study thus points to new rich physics in such experimental domain.

.6 Conclusion

We have formulated a theory of an out-of-equilibrium open quantum system set-up where the system is connected to multiple baths and is driven by periodically varying thermodynamic parameters (temperature or chemical potential) of the baths. The key assumptions in our theory are weak system-bath coupling (Born approx) and that the relaxation time scale of the bath is much smaller than the time period of the drive (Markov approx). We do not make any assumption on the relaxation time scales of the system. This has allowed us to go beyond the adiabatic approximation where the system is always assumed to be in the steady state with the instantaneous bath. Our theory works for systems of non-interacting particles (i.e, system Hamiltonian is quadratic) bilinearly coupled to multiple baths (also quadratic Hamiltonian) in a lattice of any dimension or geometry (the particles may be all bosonic or all fermionic, boson-fermion couplings are not bilinear). Even though there exists formally exact methods of treating such periodically driven microscopic models, they are generally quite difficult to calculate beyond the adiabatic regime. On the other hand, phenomenological models based on Lindblad equations have been shown to be either quite restrictive or to have severe shortcomings. Unlike these, under Born-Markov approximation, our theory gives a much simpler, more transparent and much less restrictive results, providing a complete set of linear differential equations for equal time two-point correlation functions from which various physical observables (for example currents and populations) can be directly obtained. Also various physical effects can be directly read-off from such evolution equation for correlation functions. The long time dynamics of the system are readily shown to be periodic with the same period as the drive. The occurrence of resonance when the frequency of the periodic drive equals the difference of two eigen-energies of the system is also easily seen from the equations.

We have tested our theory by applying it to the case of two fermionic sites with hopping between them, weakly coupled to two different baths at same constant temperature but with sinusoidally driven chemical potentials. In this we have considered two cases. We have validated our theory by comparing with exact numerics done with finite but large baths for both the cases. The first case is when the sinusoidal drive is symmetric, i.e, both chemical potentials vary in an exactly same manner so that there is no instantaneous voltage difference. In this case, even though there is no net time-period averaged current, there is an instantaneous current, which can become quite large at resonance. Furthermore, inside the system, the instantaneous current depends on asymmetry of the system-bath coupling and becomes zero when system-bath coupling is symmetric. This gives an experimental way of detecting the asymmetry of system-bath coupling. The maximum instantaneous current increases with amplitude of drive and shows resonance peaks. Increase in temperature washes out the higher resonance peaks.

In the other case, the chemical potentials of the baths have an amplitude difference, as well as a phase difference. Due to amplitude difference, the net time-period averaged current is non-zero and increases with increase in amplitude difference. We find interesting hysteresis in the I-V curves, a detailed investigation of which will be taken up in a future work.

Two most important experimentally relevant insights from our results, especially for studying electronic transport in molecular junctions or through quantum dots, are : (i) time-dependent measurements of current for symmetric oscillating voltages (with zero instantaneous voltage bias) can point to the degree of asymmetry in the system, and (ii) under certain conditions, time-dependent currents can exceed time-averaged currents by several orders of magnitude, and can therefore be detected even when the average current is below the measurement threshold.

Since the formalism is for non-interacting fermions (bosons) on a lattice, bilinearly coupled to fermionic (bosonic) baths, interactions can be built in at a mean field level. Thus our theory can be used to describe a wide range of experimental set-ups (such as molecular junctions, quantum dots, cavity-QED expts, cold atoms, cavity optomechanics etc.). However, the theory cannot be applied to extremely low temperatures, because of the Markov approximation. The relaxation time scale of a bath can be shown to be inversely proportional to the temperature of the bath. Hence, at extremely low temperatures, the Markov approximation would be violated. Also, our formalism cannot treat cases where system-bath coupling is not bilinear, such as fermions connected to bosonic baths and vice-versa. Further work is required for go beyond these limitations. Other future work includes treating interacting systems beyond mean-field, as well as incorporating the effects of dephasing or local vibrations.

Acknowledgements:

Y.D. acknowledges support from the Israel Science Fund (Grant No. 1256/14). A.P would like to thank Prof Abhishek Dhar for very insightful discussions, and also for funding support from his grant under the Indo-Israel joint research project No.6-8/2014(IC).

Appendix A Derivation of RQME for time independent case

The Redfield equation for our set-up in the dc bias case has been rigorously deriven in Ref Purkayastha et al. (2016). We include the same derivation here for completion.

The standard Redfield equation is derived assuming weak-system bath coupling (Born approximation) Carmichael (2002); Breuer and Petruccione (2007). For bilinear system-bath coupling $\hat{\mathcal{H}}_{SB}=\sum_{\ell}\hat{S}_{\ell}\hat{B}_{\ell}$ where $\hat{S}_{\ell}$ operates on system and $\hat{B}_{\ell}$ operates on bath, system-bath coupling being turned on at time $t=0$, we have the master equation Carmichael (2002)

[TABLE]

where we use the interaction representation $\hat{O}^{I}(t)=e^{i(\hat{\mathcal{H}}_{S}+\hat{\mathcal{H}}_{B})t}\hat{O}e^{-i(\hat{\mathcal{H}}_{S}+\hat{\mathcal{H}}_{B})t}$ and $\langle...\rangle_{B}$ refers to the average taken only with respect to bath. For our set-up in Eq. 1, we have the corresponding system and bath operators in interaction picture $\hat{a}_{\ell}^{I}(t)=\sum_{\alpha=1}^{N}c_{\ell\alpha}\hat{A}_{\alpha}e^{-\omega_{\alpha}t}$, and $\hat{B}_{\ell r}^{I}(t)=\hat{B}_{\ell r}e^{-\Omega_{\ell r}t}$, where we have used definitions in Eq. 6. Putting these in Eq. A and going back to Schroedinger picture, one obtains

[TABLE]

where $\Omega_{min}^{\ell},\Omega_{max}^{\ell}$ are the minimum and the maximum energy levels of the bath coupled to the $\ell$th site. This equation still has an integration over time from [math] to observation time $t$ and hence has time-dependent rates. To obtain a QME with time independent rates, one assumes $t\gg\tau_{B}$ (Markov approximation). This allows one to take the $t\rightarrow\infty$ in the integration limit. Upon doing this, we obtain the microscopically derived Born-Markov RQME with time independent rates. Multiplying this equation by $\hat{A}_{\alpha}^{\dagger}\hat{A}_{\nu}$ and taking trace, one can obtain the closed set of linear differential equations for the equal time two point correlation functions. This equation will be exactly same as Eq. .3 except that $Q$ will be time independent.

Appendix B Finding $\tau_{B}$

The condition on time scales Eq. 4 is one of the major requirements for validity of our theory. This entails that we must find a good estimate of bath relaxation time $\tau_{B}$. We have seen above that to derive the RQME with time independent rates, we need observation time $t\gg\tau_{B}$, so that we can take upper limits of certain integrals to infinity. We can thus look more carefully at the integrals to check under what conditions of bath parameters such an assumption is valid. Looking at Eq. 28, and using the spectral function in Eq. 22 and fermi distribution, we find that the relevant integrals are $\int_{0}^{t}\mathcal{I}_{1}(t^{\prime})dt^{\prime}$ and $\int_{0}^{t}\mathcal{I}_{2}(t^{\prime})dt^{\prime}$, where

[TABLE]

The bath relaxation time scale $\tau_{B}$ is essentially the time in which $\mathcal{I}_{1}(t)$ and $\mathcal{I}_{2}(t)$ decay. If the observation time is much larger than $\tau_{B}$, while doing the time integrals, the upper limit can be extended to infinity, because, the contribution to the time integral from time greater than $\tau_{B}$ will be negligible. Using contour integrals, $\mathcal{I}_{1}(t)$ and $\mathcal{I}_{2}(t)$ can be analytically evaluated. If $t\gg(2t_{B})^{-1},(2t_{B}-\mu)^{-1}$, the leading order term for $\mathcal{I}_{1}(t)$ is given by

[TABLE]

where $r$ are integers and $r^{*}=\frac{\beta}{2\pi}(2t_{B}-\mu)-\frac{1}{2}$. The validity of this result is shown in Fig. 4, where $\mathcal{I}_{1}(t)$ calculated via exact numerical integration and via above equation are plotted for two different temperatures. Similar result with opposite sign holds for $\mathcal{I}_{2}(t)$. From this result, it is clear that, if $2t_{B}$ is quite large, then $\tau_{B}\sim\beta$, and is independent of $\mu$. However, in our AC voltage drive case, we want to vary $\mu$ as a function of time, and make the Markov approximation always. So, in our case, $\tau_{B}\sim\beta$. This means, as mentioned in the main text, the condition on time scales (Eq. 4) render our theory inapplicable at very low temperatures.

We would like to mention that while this result was calculated explicitly for our chosen spectral function and fermi distribution, similar results yielding $\tau_{B}\sim\beta$ can also be obtained for other popular choices of spectral functions (for example, $J(\omega)~{}=~{}const$, $J(\omega)~{}\propto~{}\frac{\omega}{\omega^{2}+\Lambda^{2}}$), as well as for bose distribution.

Appendix C Deriving Eq. .3 rigorously from the protocol

Now, we would like to demonstrate that, following through our protocol in Section .2 rigorously using the time-independent RQME leads to Eq. .3. Between the $(r-1)$th and the $r$th steps of the protocol, the evolution equation for correlation matrix is as given by time-independent RQME

[TABLE]

The correlation matrix at time $r\tau_{D}$ is then given by the formal solution of this equation which can be written as

[TABLE]

If $\tau_{D}$ is small, $t^{\prime}$ is also small. But $r\tau_{D}$ can still be large for $r\gg 1$. In such case, we can expand each term to the linear order in $\tau_{D}$ to obtain

[TABLE]

This equation is the discrete time version of Eq. .3 and for $\tau_{D}\rightarrow 0$ and $r\tau_{D}\rightarrow t$, exactly gives Eq. .3.

Eq. .3 may also be obtained by ‘putting by hand’ time dependence in the fermi or bose distributions that appear in the dc bias RQME. However, note that this form would not have been possible without the Markov approx $\tau_{D}\gg\tau_{B}$. Without Markov assumption, the matrix $M$ would be time dependent which will not lead to a formal solution of the form Eq. 33, and as a consequence would not lead to Eq. .3.

Appendix D Derivation of spectral function of bath Eq. 22

We now derive the bath spectral function Eq. 22 for the explicit bath Hamiltonian in Eq. 21. In Eq. 1, the system couples to eigenmodes of the bath. Hence we can obtain the spectral function by using the eigen-energies and the eigenvectors of the bath. Since the bath has infinite number of modes, the spectrum can be assumed continuous. With this, the bath eigen-energies are $\Omega(q_{\ell})=-2t_{B}\cos q_{\ell}$ and system-bath coupling to bath modes are given by $\kappa(q_{\ell})=\gamma_{\ell}\sqrt{\frac{2}{\pi}}\sin q_{\ell}$, with $0\leq q_{\ell}\leq\pi$. Thus,

[TABLE]

We also need to the following result $\mathcal{P}\int\frac{}{}\frac{d\omega^{\prime}J_{\ell}(\omega^{\prime})}{2\pi(\omega^{\prime}-\omega)}=-\frac{\gamma_{\ell}^{2}\omega}{2t_{B}^{2}}$ to calculate $f_{\alpha\nu}^{\Delta}(\omega)$ [in Eq.  .3]. $F_{\alpha\nu}^{\Delta}(\omega)$ cannot be written in a simple closed form and is calculated numerically.