# Coarse-graining master equation for periodically driven systems

Ronja Hotz1 Gernot Schaller1,2<EMAIL_ADDRESS>1 Institut für Theoretische
Physik, Technische Universität Berlin, Hardenbergstr. 36, 10623 Berlin,
Germany  2 Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstraße 400,
01328 Dresden, Germany

###### Abstract

We analyze Lindblad-Gorini-Kossakowski-Sudarshan-type generators for selected
periodically driven open quantum systems. All these generators can be obtained
by temporal coarse-graining procedures, and we compare different coarse-
graining schemes. Similar as for undriven systems, we find that a dynamically
adapted coarse-graining time, effectively yielding non-Markovian dynamics by
interpolating through a series of different but invididually Markovian
solutions, gives the best results among the different coarse-graining schemes,
albeit at highest computational cost.

##  I Introduction

While the propagation of undriven closed quantum systems via Schrödingers
equation may be challenging for many degrees of freedom, it can be calculated
with standard methods. This changes when the Hamiltonian is subject to an
external driving, when time-ordering becomes relevant. Even for the simple
case of periodic driving, where Floquet theory Floquet (1883); Shirley (1965)
applies, the exact calculation of the unitary propagator may be notoriously
difficult Grifoni and Hänggi (1998); Eckardt and Anisimovas (2015); Bukov  _et
al._ (2015). In addition, it is then often hidden that the time-dependent
variation of external parameters requires that the system is coupled in some
way to an outside world, which to be complete would actually require to model
the system as open. While open quantum systems in absence of driving are well-
studied and interesting-subject on their own Weiss (1993); Breuer and
Petruccione (2002); Schlosshauer (2007), the joint discussion of the effects
of system-reservoir coupling and periodic driving is challenging.

However, with in particular periodic types of driving being experimentally
quite feasible with current standards Chu (1991); Rice and Zhao (2000);
Aidelsburger  _et al._ (2013); Atala  _et al._ (2014); Jotzu  _et al._ (2014),
the derivation of proper dissipators for periodically driven open systems is
of great concern. In principle, it is formally possible to perform the same
approximations that lead to Lindblad-Gorini-Kossakowski-Sudarshan (LGKS)
master equations Lindblad (1976); Gorini  _et al._ (1976) for an undriven
system also for its driven version Kohler  _et al._ (1997); Grifoni and Hänggi
(1998); Gelbwaser-Klimovsky  _et al._ (2013); Szczygielski  _et al._ (2013);
Szczygielski (2014). In the interaction picture, the application of Born-,
Markov-, and secular approximations generically leads to LGKS master equations
that come in handy due to their stability and – despite the lack of
thermalization Shirai  _et al._ (2015, 2016) – their transparent thermodynamic
interpretation Cuetara  _et al._ (2015). However, it has also frequently been
noted that in particular the involved secular approximation may lead to strong
artifacts Restrepo  _et al._ (2019); Luo  _et al._ (2020). An intuitive
explanation for this is that the level spacing in the extended (Floquet) space
may become very small, which conflicts with the secular approximation Hone
_et al._ (2009). Therefore, traditional Floquet-LGKS approaches have been
questioned in many works Restrepo  _et al._ (2019); Schnell  _et al._ (2020).

Already for undriven systems, it can be argued that the (non-LGKS) Redfield
equation applied in the proper regime only leads to small violations of
density matrix properties in trade for a closer description of physical
reality Hartmann and Strunz (2020). As a practical benefit, improper density
matrices can then be used as an indicator for leaving the region of validity
of such perturbative schemes, whereas such a witness is missing for LGKS-based
approaches.

Figure 1: (Color Online) Sketch illustrating the different coarse-graining
timescales to derive a master equation in the weak coupling limit. For period
coarse graining (PCG) (orange) the coarse-graining extends over one period of
the driving field $\tau=T$. For dynamical coarse-graining (DCG) (purple) the
coarse-graining time is chosen dynamically as the physical time. For
$\tau\to\infty$ all methods recover the Born-Markov-secular (BMS) limit (red).

Another viewpoint is to preserve the formal LGKS form of the generator while
improving in some ways the involved approximations Whitney (2008). There are
various routes to achieve this, e.g. formally restoring the LGKS form by
discarding its subspace with negative eigenvalues Kiršanskas  _et al._ (2018);
Ptaszyński and Esposito (2019), by performing only a partial secular
approximation Farina and Giovannetti (2019); Cattaneo  _et al._ (2020), or by
taking the effects of the reservoir on the system pointer basis into account
Schultz and von Oppen (2009); Schultz (2010); Teretenkov (2021); Trushechkin
(2021). While these approaches apply to undriven systems, we would like to
focus on an approach that can be easily combined with periodic driving.

Such an approach that can also be well-motivated using a microscopic
derivation is temporal coarse-graining Vacchini (2000); Lidar  _et al._
(2001); Schaller and Brandes (2008); Schaller  _et al._ (2009); Benatti  _et
al._ (2009, 2010); Royer (1996); Schaller  _et al._ (2009); Majenz  _et al._
(2013); Rivas (2017, 2019); Hartmann and Strunz (2020); Mozgunov and Lidar
(2020). Intuitively, it can be understood as an approximate Markovian
representation of the reduced exact solution over a chosen coarse-graining
time interval. The freedom in choosing the coarse-graining interval is
exemplified in Fig. 1 and shall be the subject of this paper with a focus on
periodically driven systems. We begin with a brief introduction into the
method in Sec. II and then compare various dissipators for three examples of
driven open systems in Secs. III, IV, and V.

##  II Coarse-graining methods for periodically driven open systems

The starting point of our considerations is a representation of the
Hamiltonian of the full universe in the form

$\displaystyle H(t)=H_{S}(t)+\sum_{\alpha}A_{\alpha}\otimes B_{\alpha}+H_{B}$
(1)

with periodically driven system Hamiltonian $H_{S}(t+T)=H_{S}(t)$, system and
bath coupling operators $A_{\alpha}$ and $B_{\alpha}$, respectively, and
reservoir Hamiltonian $H_{B}$. We do not necessarily require the coupling
operators to be individually hermitian, but the total interaction Hamiltonian
is hermitian $H_{I}=\sum_{\alpha}A_{\alpha}\otimes
B_{\alpha}=H_{I}^{\dagger}$. Furthermore, we assume the system coupling
operators $A_{\alpha}$ to be normalized with respect to some norm, such that
unit and strength of the interaction is carried by the bath coupling operators
$B_{\alpha}$. Since system and bath operators act on different Hilbert spaces,
we assume from the beginning that $[H_{S},H_{B}]=0$ and also
$[A_{\alpha},B_{\beta}]=0$, $[A_{\alpha},H_{B}]=0$, $[B_{\alpha},H_{S}]=0$. In
an interaction picture with respect to $H_{0}=H_{S}(t)+H_{B}$ (marked by bold
symbols), we can approximate the time evolution operator
$\mbox{\boldmath$U$}(t_{0}+\tau,t_{0})$. The coarse-graining approach
originates from the attempt to match the reduced evolution of the open quantum
system with a time-local master equation

$\displaystyle e^{\mbox{\boldmath${\cal
L}_{\tau}$}\tau}\rho_{S}^{0}\stackrel{{\scriptstyle!}}{{=}}{\rm
Tr_{B}}\left\\{\mbox{\boldmath$U$}(t_{0}+\tau,t_{0})\rho_{S}^{0}\otimes\rho_{B}^{0}\mbox{\boldmath$U^{\dagger}$}(t_{0}+\tau,t_{0})\right\\}\,.$
(2)

Here, we have already put an index ${\cal L}_{\tau}$ to the dissipation
superoperator to highlight the fact that with a constant superoperator ${\cal
L}$ it will in general not be possible to capture the exact reduced dynamics
for all times, which would require a Kraus representation Kraus (1971). For
weak system-reservoir interaction (or small propagation times $\tau$) we can
use that $\mbox{\boldmath$U$}(t_{0}+\tau,t_{0})$ is close to the identity and
$\mbox{\boldmath${\cal L}$}_{\tau}$ is small. One can then expand the
expressions on both sides and assuming furthermore that ${\rm
Tr}\left\\{B_{\alpha}\rho_{B}^{0}\right\\}=0$, one arrives at a defining
equation for the coarse-graining generator, see App. A.1 for details. Choosing
for simplicity $t_{0}=0$ it reads

${\cal L}_{\tau}$$\rho_{S}$
$\displaystyle=-\mathrm{i}\Big{[}\frac{1}{2\mathrm{i}\tau}\int\limits^{\tau}_{0}dt_{1}dt_{2}\sum_{\alpha\bar{\alpha}}C_{\alpha\bar{\alpha}}(t_{1},t_{2})$
(3) $\displaystyle\qquad\qquad\times{\rm
sgn}(t_{1}-t_{2})\boldsymbol{A}_{\alpha}(t_{1})\boldsymbol{A}_{\bar{\alpha}}(t_{2}),\boldsymbol{\rho_{S}}\Big{]}$
$\displaystyle\qquad+\frac{1}{\tau}\int\limits^{\tau}_{0}dt_{1}dt_{2}\sum_{\alpha\bar{\alpha}}C_{\alpha\bar{\alpha}}(t_{1},t_{2})$
$\displaystyle\qquad\times\Big{[}\boldsymbol{A}_{\bar{\alpha}}(t_{2})\boldsymbol{\rho_{S}}\boldsymbol{A}_{\alpha}(t_{1})-\frac{1}{2}\\{\boldsymbol{A}_{\alpha}(t_{1})\boldsymbol{A}_{\bar{\alpha}}(t_{2}),\boldsymbol{\rho_{S}}\\}\Big{]}\,,$

with bath correlation functions $C_{\alpha\bar{\alpha}}(t_{1},t_{2})={\rm
Tr_{B}}\\{\boldsymbol{B_{\alpha}}(t_{1})\boldsymbol{B_{\bar{\alpha}}}(t_{2})\rho_{B}^{0}\\}={\rm
Tr}\left\\{\mbox{\boldmath$B_{\alpha}(t_{1}-t_{2})$}B_{\beta}\rho_{B}^{0}\right\\}$,
where the last equality holds when $[H_{B},\rho_{B}^{0}]=0$. Here, the
interaction picture representation of the system coupling operators is
formally given by
$\mbox{\boldmath$A_{\alpha}$}(t)=U_{S}^{\dagger}(t)A_{\alpha}U_{S}(t)$ with
$U_{S}(t)=\mathcal{T}\left\\{e^{-\mathrm{i}\int_{0}^{t}H_{S}(t^{\prime})dt^{\prime}}\right\\}$
and $\mathcal{T}$ denoting the time ordering operator. One can show that for
any kind of driving and also for any fixed choice of $t_{0}$ and $\tau$, the
above coarse-graining dissipator ${\cal L}_{\tau}$ is of LGKS form Mozgunov
and Lidar (2020); Schaller and Ablaßmayer (2020). This contrasts the above
generator from other coarse-graining approaches that arise from a temporal
averaging of Redfield equations, which require subsequent secular-type
approximations to reach an LGKS form Elouard  _et al._ (2020).

Thus, we can generally split the generator into a unitary part and a
dissipative one

$\displaystyle\boldsymbol{\dot{\rho}_{S}}$
$\displaystyle=-\mathrm{i}\left[\mbox{\boldmath$H_{\rm
LS}^{\tau}$},\boldsymbol{\rho_{S}}\right]+\mbox{\boldmath$\mathcal{D}^{\tau}$}\mbox{\boldmath$\rho_{S}$}\,,$
(4)

which specifically for periodic driving have the single-integral expressions

$H_{\rm LS}^{\tau}$ $\displaystyle=\int
d\omega\sum_{\alpha\bar{\alpha}}\sigma_{\alpha\bar{\alpha}}(\omega)\frac{1}{2\mathrm{i}}\sum_{nn^{\prime}}\sum_{abcd}$
$\displaystyle\qquad\times
f_{0}^{\tau}(\bar{E}_{a}-\bar{E_{b}}+n\Omega,\bar{E}_{c}-\bar{E}_{d}+n^{\prime}\Omega,\omega)$
$\displaystyle\qquad\times
A_{\alpha,ab}^{n}A_{\bar{\alpha},dc}^{-n^{\prime}}L_{ab}L_{dc}\,,$
$\displaystyle\mbox{\boldmath$\mathcal{D}^{\tau}$}\boldsymbol{\rho_{S}}$
$\displaystyle=\int
d\omega\sum_{\alpha\bar{\alpha}}\gamma_{\alpha\bar{\alpha}}(\omega)\sum_{nn^{\prime}}\sum_{abcd}$
$\displaystyle\qquad\times
f_{0}^{\tau}(\bar{E}_{a}-\bar{E}_{b}+n\Omega,\bar{E}_{c}-\bar{E}_{d}+n^{\prime}\Omega,\omega)$
$\displaystyle\qquad\times
A_{\alpha,ab}^{+n}A_{\bar{\alpha},dc}^{-n^{\prime}}$
$\displaystyle\qquad\times{\left[L_{dc}\boldsymbol{\rho_{S}}L_{ab}-\frac{1}{2}\left\\{L_{ab}L_{dc},\boldsymbol{\rho_{S}}\right\\}\right]}\,.$
(5)

Here, the first term converges for large $\tau$ to the LGKS Lambshift term,
and the second term denotes the dissipative influence of the reservoir. The
functions $\sigma_{\alpha\bar{\alpha}}(\omega)$ as well as
$\gamma_{\alpha\bar{\alpha}}(\omega)$ denote odd and even Fourier transforms
of the reservoir correlation functions (A.1),
$f_{0}^{\tau}(E_{a},E_{b},\omega)$ is a nascent $\delta$-function, the
$A_{\alpha,ab}^{n}$ are Fourier components of the system coupling operator
matrix elements, the quasienergies $\bar{E}_{a}\in[-\Omega/2,+\Omega/2)$ are
chosen in the first Brillouin zone defined by the period of the driving
$\Omega=2\pi/T$, and the Lindblad operators $L_{ab}$ generate transitions
between Floquet Hamiltonian eigenstates, see App A.1 for the precise
definitions.

The different coarse-graining schemes sketched in Fig. 1 may lead to different
LGKS-type dissipators. For later reference, we coin the coarse-graining
schemes dynamical coarse-graining Schaller and Brandes (2008) (DCG) with
$t_{0}=0$ and $\tau=t$ and period-coarse-graining (PCG) with $t_{0}=0$ and
$\tau=T$. Details for these schemes are found in Appendices A.3 and A.4,
respectively. Furthermore, we remark that the limit of very large coarse-
graining times $\tau\to\infty$ allows for significant simplifications, see
App. A.5, which however do not always coincide with the frequently-used Born-
Markov-secular (BMS) Floquet master equation exposed in App. A.6 and a further
simplified variant thereof – the Born-Markov-ultrasecular (BMU) master
equation – in App. A.7.

Our problem is defined by three physical timescales. We will be concerned with
driven two-level systems, where the natural internal time scale of the system
dynamics is defined by their level splitting $\Delta E$ via the relation
$\tau_{\rm int}=2\pi/\Delta E$. Then, we have the obvious timescale given by
the period of the external driving $\tau_{\rm drv}=2\pi/\Omega=T$. Finally,
the reservoir timescales are reflected in the decay timescale of the reservoir
correlation function $\tau_{\rm dec}=2\pi/\omega_{\rm c}$, where $\omega_{\rm
c}$ is a cutoff in the reservoir spectral density. We are going to investigate
how the technical timescale – the coarse-graining time $\tau$ – can be
optimally chosen in dependence of the physical ones.

In what follows, we will compare the solutions to various periodically driven
problems that are based on these different coarse-graining approaches, i.e.,
for the DCG approach we evaluate Eq. (II) with $\tau=t$, for the PCG approach
we use Eq. (II) with $\tau=T$, whereas for the BMS and BMU results we use
Eqns. (A.6) and (A.7), respectively. Furthermore, in our calculations we will
for simplicity neglect the Hermitian correction term to the Hamiltonian
$H_{\rm LS}^{\ldots}$ throughout. Whereas for the particular pure-dephasing
model considered below this term has no effect anyways, in general it can only
be neglected in the weak-coupling regime where ${\left|H_{\rm
LS}^{\ldots}\right|}\ll{\left|H_{S}\right|}$.

##  III Pure-dephasing models

By the term pure-dephasing models we summarize models where the interaction
commutes with the system Hamiltonian at all times, i.e., system and reservoir
cannot exchange energy. The Hamiltonian of a general driven-system pure-
dephasing model with a reservoir of bosonic oscillators is then given by:

$\displaystyle H(t)$
$\displaystyle=H_{S}(t)+A\otimes\sum_{k}\left(h_{k}b_{k}+h_{k}^{*}b_{k}^{\dagger}\right)+\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}\,,$
$\displaystyle\left[H_{S}(t),A\right]=0\,.$ (6)

The assumption of commuting coupling operator $A$ and system Hamiltonian
$H_{S}(t)$ allows to solve the dynamics exactly, see App. B.

###  III.1 Exact solution

Specifically for a driven two level system with
$H_{S}=\sigma^{z}\left[\frac{\Delta}{2}+\lambda\cos(\Omega t)\right]$ with the
single system coupling operator $A=\sigma^{z}$ we obtain that the populations
in the $\sigma^{z}$ eigenbasis remain constant

$\displaystyle\rho_{S,00}(t)=\rho_{S,00}^{0}\,,\qquad\rho_{S,11}(t)=\rho_{S,11}^{0}\,,$
(7)

whereas the coherences decay oscillatorily according to

$\displaystyle\rho_{S,10}(t)$ $\displaystyle=e^{+\mathrm{i}\left(\Delta\cdot
t+2\frac{\lambda}{\Omega}\sin(\Omega t)\right)}$ (8)
$\displaystyle\qquad\times\exp\left\\{-\frac{4}{\pi}\int\limits_{-\infty}^{\infty}d\omega\gamma(\omega)\frac{\sin^{2}\left(\frac{\omega
t}{2}\right)}{\omega^{2}}\right\\}\rho_{S,10}^{0}\,.$

Here, the Fourier transform of the single reservoir correlation function is
$\gamma(\omega)=\Gamma(\omega)[1+n_{B}(\omega)]$ with spectral density
$\Gamma(\omega)\equiv
2\pi\sum_{k}{\left|h_{k}\right|}^{2}\delta(\omega-\omega_{k})$ (analytically
continued via $\Gamma(-\omega)=-\Gamma(\omega)$ to the complete real axis) and
the Bose distribution is $n_{B}(\omega)=[e^{\beta\omega}-1]^{-1}$. In absence
of driving ($\lambda\to 0$), this falls back to the known pure-dephasing
solution of the spin-boson model, see e.g. Ref. Lidar  _et al._ (2001).
Additionally, we remark that one may also obtain the above result from the
undriven model by a gauge transform.

###  III.2 Master equation solutions

We can now compare this exact solution with the various approximate approaches
discussed before by dropping the indices as above, i.e., using
$A_{1}=A=\sigma^{z}$ and
$\gamma_{11}(\omega)=\gamma(\omega)=\Gamma(\omega)[1+n_{B}(\omega)]$ for
coupling operator and Fourier transform of the reservoir correlation function.
Notably, the matrix elements of the coupling operators become rather trivial
$A_{1,ab}^{n}=\delta_{n,0}\delta_{ab}\bra{\bar{a}}\sigma^{z}\ket{\bar{a}}$,
where $\ket{\bar{a}}$ are Floquet Hamiltonian eigenstates. To compute them, it
suffices to realize that for this simple problem the system Floquet
Hamiltonian is just $\bar{H}=\frac{\Delta}{2}\sigma^{z}$, such that its
eigenstates are trivial. All master equations capture the constant populations
(7).

The coarse-graining solution for fixed coarse-graining time $\tau$ predicts
for the coherences a decay according to:

$\displaystyle\rho_{cg,10}(t)$ $\displaystyle=e^{+\mathrm{i}\left(\Delta\cdot
t+2\frac{\lambda}{\Omega}\sin(\Omega t)\right)}$ (9)
$\displaystyle\qquad\times\exp\left\\{-\frac{4}{\pi}{\frac{t}{\tau}}\int\limits_{-\infty}^{\infty}\gamma(\omega)\frac{\sin^{2}(\frac{\omega{\tau}}{2})}{\omega^{2}}d\omega\right\\}\rho^{0}_{10}\,,$

where the PCG approach is recovered for $\tau=2\pi/\Omega$.

Clearly, for DCG ($\tau=t$), the above equation exactly reproduces the exact
solution (8).

Furthermore, for this model the BMS and BMU solutions coincide:

$\displaystyle\rho_{\rm BM,10}(t)=e^{+\mathrm{i}\left(\Delta\cdot
t+2\frac{\lambda}{\Omega}\sin(\Omega
t)\right)}e^{-2\gamma(0)t}\rho^{0}_{10}\,,$ (10)

which can also be obtained by virtue of (40) from Eq. (9) in the limit
$\tau\to\infty$.

###  III.3 Comparison

We analytically find that the DCG approach (adaptively choosing $\tau=t$ and
$t_{0}=0$ in Eq. (II)) always reproduces the exact solution. We also find
analytically from the trivial structure of the model that BMS (A.6) and BMU
(A.7) solutions must concide. In Fig. 2 one can see that the additional
driving modifies the decay of coherences by super-imposing small additional
oscillations to the undriven solution (green dashed).

Figure 2: (Color Online) Plot of $\left<\sigma^{x}\right>/2$ comparing coarse-
graining PCG for $\tau=T=\frac{2\pi}{\Omega}$ (orange) and the BMS/BMU
solution (red) to the exact analytical solution (blue) and the exact solution
in absence of driving (dashed green). BMS and BMU solutions are identical, and
the DCG approach ($\tau=t$) reproduces the exact solution. Inset: By
construction, PCG (orange) and DCG solutions (blue) intersect at $t=T$
(vertical dashed line). Parameters: $\Omega=10\;\Delta$,
$\lambda=\frac{1}{2}\;\Delta$, $\Gamma(\omega)=\Gamma_{0}\omega
e^{-{\left|\omega\right|}/\omega_{\rm c}}$ with $\Gamma_{0}=0.05$,
$\omega_{c}=20\;\Delta$, $\beta\,\Delta=1$, $\rho_{\rm
S}^{0}=1/2(\mbox{\boldmath$1$}+\sigma^{x})$.

Apart from that, we see that the exact reduced dynamics cannot be reproduced
by a single Markovian generator (like PCG, BMS or BMU). Fig. 2 has been
computed for parameters with $\omega_{c}>\Omega\gg\Delta$ ($\tau_{\rm
dec}<\tau_{\rm drv}\ll\tau_{\rm int})$.

When we consider slower drivings $\omega_{c}\gg\Omega=\Delta$ ($\tau_{\rm
dec}\ll\tau_{\rm drv}=\tau_{\rm int}$), we observe that also the PCG method
converges towards the exact solution, see Fig. 3, whereas the BMS approach
fails to capture the exact dynamics.

Figure 3: (Color Online) Analogous to Fig. 2, but for slower driving
$\Omega=\Delta$. PCG and exact solution agree much better.

The failure of the BMS approach in this limit is not too surprising, given
that it is typically derived for fast driving. The fact that for slow driving
it suffices to coarse-grain over the first period (PCG) to get a sufficiently
accurate dissipator is a consequence of the relatively high temperatures
chosen. This can be formally seen by performing a variable transform
$\omega=\omega^{\prime}t/\tau$ in (9) and using that in the considered range
of parameters $\gamma(\Omega
t/(2\pi)\omega)\approx\Gamma_{0}/\beta\approx\gamma(\omega)$. For lower
temperatures (not shown) and other models (see below) we do not see such a
clear tendency.

Since pure-dephasing-type models are a very specific model class, we consider
models that allow the exchange of energy between system and reservoir below.

##  IV Circular driving

As a second model, we consider a periodically driven open two level system of
the form

$\displaystyle H(t)$
$\displaystyle=\frac{\Delta}{2}\sigma^{z}+P\sigma^{+}e^{-\mathrm{i}\Omega
t}+P^{*}\sigma^{-}e^{+\mathrm{i}\Omega t}$
$\displaystyle\quad+\sum_{k}(\sigma^{+}h_{k}b_{k}+\sigma^{-}h_{k}^{*}b_{k}^{\dagger})+\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}\,,$
(11)

where $P$ denotes a driving amplitude. For $P=0$ and a vacuum reservoir, the
model is exactly solvable Garraway (1997); Breuer and Petruccione (2002), but
we are here interested in the case of finite driving amplitudes but weak
system-reservoir coupling. We therefore interpret the DCG solution as
benchmark for this model, see App. A.2.

Before treating the model with our standard approaches, we demonstrate that
the model can be mapped to a different frame where the Hamiltonian becomes
time-independent. In this time-independent frame, a standard master equation
can be derived that is equivalent to the Floquet (BMS) master equation in the
original frame.

###  IV.1 Mapping to a time-independent frame

By employing the unitary gauge transform

$\displaystyle V(t)=\exp\left[-\mathrm{i}\left(\frac{\Omega
t}{2}\sigma^{z}+\Omega t\sum_{k}b_{k}^{\dagger}b_{k}\right)\right]$ (12)

one can see from the observations

$\displaystyle V^{\dagger}(t)\sigma^{z}V(t)$
$\displaystyle=\sigma^{z}\,,\qquad
V^{\dagger}(t)\sigma^{\pm}V(t)=\sigma^{\pm}e^{\pm\mathrm{i}\Omega t}\,,$
$\displaystyle V^{\dagger}(t)b_{k}V(t)$ $\displaystyle=e^{-\mathrm{i}\Omega
t}b_{k}\,,\qquad V^{\dagger}(t)b_{k}^{\dagger}V(t)=e^{+\mathrm{i}\Omega
t}b_{k}^{\dagger}\,,$ (13)

that the total Hamiltonian becomes time-independent under the transformation

$\displaystyle V^{\dagger}(t)H(t)V(t)$
$\displaystyle=\frac{\Delta}{2}\sigma^{z}+P\sigma^{+}+P^{*}\sigma^{-}$ (14)
$\displaystyle\quad+\sum_{k}(\sigma^{+}h_{k}b_{k}+\sigma^{-}h_{k}^{*}b_{k}^{\dagger})+\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}\,.$

Accordingly, by transforming into the corresponding picture, the expectation
value of an observable $O$ can be written as

$\displaystyle\left<O\right>={\rm
Tr}\left\\{OU(t)\rho_{0}U^{\dagger}(t)\right\\}\equiv{\rm
Tr}\left\\{\tilde{O}(t)\tilde{U}(t)\rho_{0}\tilde{U}^{\dagger}(t)\right\\}$
(15)

with the effective time evolution operator $\tilde{U}(t)=V^{\dagger}(t)U(t)$
and transformed observable $\tilde{O}(t)=V^{\dagger}(t)OV(t)$. The effective
time evolution operator evolves according to a time-independent picture

$\displaystyle\frac{d}{dt}\tilde{U}(t)$
$\displaystyle=-\mathrm{i}\tilde{H}\tilde{U}(t)\,,$ (16)
$\displaystyle\tilde{H}$
$\displaystyle=V^{\dagger}(t)H(t)V(t)+\mathrm{i}\dot{V}V$
$\displaystyle=\frac{\Delta-\Omega}{2}\sigma^{z}+P\sigma^{+}+P^{*}\sigma^{-}$
$\displaystyle\quad+\sum_{k}(\sigma^{+}h_{k}b_{k}+\sigma^{-}h_{k}^{*}b_{k}^{\dagger})+\sum_{k}(\omega_{k}-\Omega)b_{k}^{\dagger}b_{k}\,,$

such that we can write $\tilde{U}(t)=e^{-\mathrm{i}\tilde{H}t}$.

Therefore, under the gauge transform (14), the time evolution is just that of
a time-independent open two-level system. In this picture, we can derive the
BMS master equation for an undriven system. Formally, this master equation
looks like Eq. (A.6) with keeping only the $n=0$ term and instead of the
Floquet Hamiltonian eigenstates and eigenvalues using the eigenstates and
eigenvalues of
$\tilde{H}_{S}=\frac{\Delta-\Omega}{2}\sigma^{z}+P\sigma^{+}+P^{*}\sigma^{-}$.

In the last line of the above equation, we can see that the reservoir single-
particle energies are shifted, which eventually also affects the Kubo-Martin-
Schwinger Kubo (1957); Martin and Schwinger (1959); Kossakowski  _et al._
(1977) relation obeyed by the reservoir correlation function. Therefore, the
resulting master equation does not thermalize even in the time-independent
picture but reaches some nonequilibrium steady state, which – after
transforming back to the original picture – exhibits a time-dependence with
the period of the driving. This becomes more evident in a simpler long-term
analysis: For non-degenerate eigenvalues of the effective system Hamiltonian
$\tilde{H}_{S}$, the populations decouple from the coherences and the
evolution is given by a simple rate equation:
$\dot{\tilde{\rho}}_{S,aa}=\sum_{b}\gamma_{ab,ab}\tilde{\rho}_{S,bb}-\sum_{b}\gamma_{ba,ba}\tilde{\rho}_{S,aa}$
with transition rates from energy level $b$ to $a$ given by $\gamma_{ab,ab}$.
The steady state solution is then given by

$\displaystyle\bar{\tilde{\rho}}_{S}$
$\displaystyle=P_{-}\ket{-}\bra{-}+(1-P_{-})\ket{+}\bra{+}\,,$ $\displaystyle
P_{-}$ $\displaystyle=\frac{\gamma_{-+,-+}}{\gamma_{-+,-+}+\gamma_{+-,+-}}\,,$
(17)

where $\ket{-}$ and $\ket{+}$ denote the ground and excited state of
$\tilde{H}_{S}$, respectively. For our example the two transition rates read

$\displaystyle\gamma_{-+,-+}$
$\displaystyle=\tilde{\gamma}_{12}(E_{+}-E_{-})|\bra{-}\sigma^{-}\ket{+}|^{2}$
$\displaystyle\qquad+\tilde{\gamma}_{21}(E_{+}-E_{-})|\bra{-}\sigma^{+}\ket{+}|^{2}\,,$
$\displaystyle\gamma_{+-,+-}$
$\displaystyle=\tilde{\gamma}_{12}(E_{-}-E_{+})|\bra{+}\sigma^{-}\ket{-}|^{2}$
$\displaystyle\qquad+\tilde{\gamma}_{21}(E_{-}-E_{+})|\bra{+}\sigma^{+}\ket{-}|^{2}\,,$
(18)

with
$\tilde{\gamma}_{12}(\omega)=\Theta(\omega+\Omega)\Gamma(\omega+\Omega)[1+n_{B}(\omega+\Omega)]$
and
$\tilde{\gamma}_{21}(\omega)=\Theta(-\omega+\Omega)\Gamma(-\omega+\Omega)n_{B}(-\omega+\Omega)$
exemplifying the broken Kubo-Martin-Schwinger relation
$\frac{\tilde{\gamma}_{21}(-\omega)}{\tilde{\gamma}_{12}(+\omega)}=e^{-\beta(\omega+\Omega)}$.
In this picture, system observables can be computed via

$\displaystyle\left<O\right>\to{\rm Tr_{S}}\left\\{e^{+\mathrm{i}\Omega
t/2\sigma^{z}}Oe^{-\mathrm{i}\Omega
t/2\sigma^{z}}\tilde{\rho}_{S}(t)\right\\}\,,$ (19)

where for large times we can insert $\tilde{\rho}(t)\to\bar{\tilde{\rho}}_{S}$
as given by Eq. (IV.1).

We find numerically that the solutions derived in this way fully agree with
the BMS solutions in the original frame (solid red curves in Figs. 4 and 5).
This is not surprising as the gauge transformation that we used should not
change the general outcome. In addition, we checked that the simplified long-
term dynamics agrees with the full time-dependent solution for large times.

###  IV.2 Master equation solutions

We now go back to the original Schrödinger picture representation of Eq. (IV),
where we can identify system coupling operators $A_{1}=\sigma^{+}$ and
$A_{2}=\sigma^{-}$. Additionally, the non-vanishing correlation functions
become $\gamma_{12}(\omega)=\Theta(\omega)\Gamma(\omega)[1+n_{B}(\omega)]$ and
$\gamma_{21}(\omega)=\Theta(-\omega)\Gamma(-\omega)n_{B}(-\omega)$. With a
gauge transform analogous to (12) but only involving the system parts one can
find the exact time evolution operator of the driven system alone, and by
demanding periodicity of the kick operator one finds that kick-operator and
Floquet Hamiltonian are given by

$\displaystyle\bar{H}$
$\displaystyle=\frac{\Delta-\Omega}{2}\sigma^{z}+P\sigma^{+}+P^{*}\sigma^{-}-\frac{\Omega}{2}\mbox{\boldmath$1$}\,,$
$\displaystyle U_{\rm kick}(t)$
$\displaystyle=\exp\left(-\mathrm{i}\frac{\Omega
t}{2}(\mbox{\boldmath$1$}+\sigma^{z})\right)\,.$ (20)

This allows us to compute the coarse-graining dissipators DCG (App. A.3) and
PCG (App. A.4) as well as BMS (App. A.6) and BMU (App. A.7) dissipators.

###  IV.3 Comparison

We compare the resulting dynamics with the previously discussed asymptotic
solution in Fig. 4.

Figure 4: (Color Online) Plot of the time-dependent expectation value
$\left<\sigma^{x}\right>$ according to the undriven DCG solution (dashed
green), driven DCG (purple crosses), PCG (orange), and BMS (solid red) and BMU
(dashed red) solutions. The gridlines mark multiples of the driving period
$T$. At long times, DCG, BMS, and BMU solutions all agree (right inset). For
short times, the BMS and BMU solutions differ from each other and also from
the DCG (benchmark) solution. Parameters: $\Omega=2\;\Delta$,
$P=\frac{1}{2}\;\Delta$ (with Floquet energies
($\bar{E}_{\pm}=\pm\frac{2-\sqrt{2}}{2}\Delta$) in the first Brillouin zone),
$\beta\,\Delta=0.1$, $\Gamma_{0}=0.05$, $\omega_{c}=15\;\Delta$,
$\rho_{S}^{0}=1/2[\mbox{\boldmath$1$}+0.6\sigma^{x}+0.4\sigma^{z}]$.

At steady state, we see a convincing agreement of the BMS and BMU (solid and
dashed red) approaches with the DCG approach. Additionally, we checked that
the BMS approach fully agrees with the master equation solution from the time-
independent picture (not shown). In the short-term limit however, some
differences become visible (left inset). In contrast, the PCG solution
(orange) captures the short-term evolution but fails to match the DCG solution
for large times.

Analogous results hold for other expectation values, see Fig. 5.

Figure 5: (Color Online) Analogous to Fig. 4 but for
$\left<\sigma^{z}\right>$. By construction, for short times the DCG solution
is superior to the other perturbative methods, but again DCG, BMS, and BMU
solutions converge to for large times.

Here one finds strong differences only transiently, where for short times the
benchmark DCG approach by construction must approximate the unknown exact
solution for the chosen parameters.

These figures have been computed for $\omega_{c}\gg\Omega>\Delta$, i.e., for
$\tau_{\rm dec}\ll\tau_{\rm drv}<\tau_{\rm int}$, which explains the good
long-term agreement of the BMS and BMU solutions.

We also explored the regime of faster driving ($\omega_{c}>\Omega\gg\Delta$),
where the BMS solution performs also transiently better and the regime of
slower driving ($\omega_{c}\gg\Omega\approx\Delta$), where the BMS solution
performed (transiently) worse. However, unlike Sec. III, the PCG approach did
not approximate the DCG results very well for all times and observables at
slow drivings.

##  V Fast-Driving solution

As a last example, we consider systems where the driving breaks the pure-
dephasing character of the system

$\displaystyle H(t)$ $\displaystyle=H_{S}^{0}+\lambda\cos(\Omega
t)C+A\otimes\sum_{k}(h_{k}b_{k}+h_{k}^{*}b_{k}^{\dagger})$
$\displaystyle\qquad+\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}\,,$ (21)

with system Hamiltonian $H_{S}^{0}$ and system coupling operator
$A=A^{\dagger}$ fulfilling only the constraint $[H_{S}^{0},A]=0$ and system
driving operator $C=C^{\dagger}$ with in general $[H_{S}^{0},C]\neq 0$. Under
a naive rotating wave approximation (RWA), the contribution of the driving
vanishes completely, and in case of a qubit the exact solution (8) from the
previous section in absence of driving ($\lambda=0$) applies for all values of
the system-reservoir coupling strength. However, to go beyond the naive RWA,
we transform into an interaction picture with respect to the driving

$\displaystyle U_{1}(t)=e^{-\mathrm{i}\lambda/\Omega\sin(\Omega t)C}\,,$ (22)

and in this picture (marked by a tilde e.g. via
$\tilde{A}(t)=U_{1}^{\dagger}(t)AU_{1}(t)$) the transformed Hamiltonian
assumes a time-dependent pure-dephasing form

$\displaystyle\tilde{H}(t)$
$\displaystyle=\tilde{H}_{S}^{0}(t)+\tilde{A}(t)\otimes\sum_{k}(h_{k}b_{k}+h_{k}^{*}b_{k}^{\dagger})+\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}\,,$
(23)

where we still have the equal-time commutation
$\left[\tilde{H}_{S}^{0}(t),\tilde{A}(t)\right]=U_{1}^{\dagger}(t)\left[H_{S}^{0},A\right]U_{1}(t)=0$.
Nevertheless, we cannot apply the naive polaron treatment of App. B, since now
also the coupling operator has picked up a periodic time-dependence. Then, the
polaron transform would become time-dependent as well, leading to an
additional term in the effective Hamiltonian, such that the system-reservoir
interaction would just be moved elsewhere.

However, depending on the problem one may have the situation that also the
time averages of the operators

$\displaystyle\bar{H}_{S}^{0}\equiv\frac{1}{T}\int_{0}^{T}\tilde{H}_{S}^{0}(t)dt\,,\qquad\bar{A}\equiv\frac{1}{T}\int_{0}^{T}\tilde{A}(t)dt$
(24)

commute with each other $[\bar{H}_{S}^{0},\bar{A}]=0$, and in this case the
Hamiltonian after an RWA in the interaction picture

$\displaystyle\bar{H}_{\rm
tot}=\bar{H}_{S}^{0}+\bar{A}\otimes\sum_{k}(h_{k}b_{k}+h_{k}^{*}b_{k}^{\dagger})+\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}$
(25)

is of pure-dephasing type and can be solved with standard methods, see App. C
for details.

The transformation (22) is of course not equivalent to the interaction picture
representation. However, under the RWA in this transformed frame, the time
evolution operator of the system becomes $U_{S}(t)\approx
U_{1}(t)e^{-\mathrm{i}\bar{H}_{S}^{0}t}$, such that we can identify
approximations to kick operator $U_{\rm kick}(t)\approx U_{1}(t)$ and system
Floquet Hamiltonian $\bar{H}\approx\bar{H}_{S}^{0}$, respectively.

###  V.1 Fast driving benchmark

Specifically, for a two level system with

$\displaystyle H_{S}^{0}=\frac{\Delta}{2}\sigma^{z}\,,\quad
A=\sigma^{z}\,,\quad C=\sigma^{x},$ (26)

we obtain with the Fourier decomposition

$\displaystyle\tilde{A}(t)$
$\displaystyle=U_{1}^{\dagger}(t)\sigma^{z}U_{1}(t)$
$\displaystyle=\sigma^{z}\sum_{n=-\infty}^{\infty}{\cal
J}_{2n}\left(\frac{2\lambda}{\Omega}\right)e^{\mathrm{i}2n\Omega t}$
$\displaystyle\qquad-\mathrm{i}\sigma^{y}\sum_{n=-\infty}^{\infty}{\cal
J}_{2n+1}\left(\frac{2\lambda}{\Omega}\right)e^{\mathrm{i}(2n+1)\Omega t}$
(27)

that this method – under the RWA in the transformed frame – yields the
following dynamics

$\displaystyle\left<\sigma^{x}\right>$
$\displaystyle=\big{[}\cos\left(\mu_{1}\Delta\cdot
t\right)\left<\sigma^{x}\right>_{0}-\sin\left(\mu_{1}\Delta\cdot
t\right)\left<\sigma^{y}\right>_{0}\big{]}\Sigma(t)\,,$
$\displaystyle\left<\sigma^{y}\right>$
$\displaystyle=-\sin\left(\tfrac{2\lambda}{\Omega}\sin(\Omega
t)\right)\left<\sigma^{z}\right>_{0}+\cos\left(\tfrac{2\lambda}{\Omega}\sin(\Omega
t)\right)$ $\displaystyle\qquad\times\big{[}\sin\left(\mu_{1}\Delta\cdot
t\right)\left<\sigma^{x}\right>_{0}+\cos\left(\mu_{1}\Delta\cdot
t\right)\langle\sigma^{y}\rangle_{0}\big{]}\Sigma(t)\,,$
$\displaystyle\left<\sigma^{z}\right>$
$\displaystyle=+\cos\left(\tfrac{2\lambda}{\Omega}\sin(\Omega
t)\right)\left<\sigma^{z}\right>_{0}+\sin\left(\tfrac{2\lambda}{\Omega}\sin(\Omega
t)\right)$ $\displaystyle\qquad\times\big{[}\sin\left(\mu_{1}\Delta\cdot
t\right)\left<\sigma^{x}\right>_{0}+\cos\left(\mu_{1}\Delta\cdot
t\right)\left<\sigma^{y}\right>_{0}\big{]}\Sigma(t)\,,$ (28)

with
$\Sigma(t)\equiv\exp\left\\{-\frac{4}{\pi}\left(\mu_{1}\right)^{2}\int\limits_{-\infty}^{\infty}\gamma(\omega)\frac{\sin^{2}\left(\frac{\omega
t}{2}\right)}{\omega^{2}}d\omega\right\\}$ and
$\mu_{1}\equiv\mathcal{J}_{0}\left(\tfrac{2\lambda}{\Omega}\right)$ being
given by the Bessel function of the first kind Abramowitz and Stegun (1970).
Details on this can be found in App. C by inserting the expressions for a two
level system in Eq. (67). In absence of of driving $\lambda\to 0$ one can
check with $\mu_{1}\to 1$ (as
$\mathcal{J}_{0}\left(\frac{2\lambda}{\Omega}\right)\to 1$ for
$\lambda\rightarrow 0$) that this falls back to the exact pure-dephasing
solution Lidar  _et al._ (2001); Leggett  _et al._ (1987) with constant
populations (7) and decaying coherences (8).

Furthermore, at vanishing system-reservoir coupling we have $\Sigma(t)\to 1$,
where one can immediately verify that the purity of an initial state is
preserved. Finally, we can see that for reasonable spectral densities we find
that $\lim_{t\to\infty}\Sigma(t)=0$, such that the long-term asymptotics of
(V.1) is given by

$\displaystyle\left<\sigma^{y}\right>$
$\displaystyle\to-\sin\left(\frac{2\lambda}{\Omega}\sin(\Omega
t)\right)\left<\sigma^{z}\right>_{0}\,,$
$\displaystyle\left<\sigma^{z}\right>$
$\displaystyle\to+\cos\left(\frac{2\lambda}{\Omega}\sin(\Omega
t)\right)\left<\sigma^{z}\right>_{0}\,.$ (29)

Importantly, we stress that the fast-driving solution is valid also for
stronger system-reservoir couplings Restrepo  _et al._ (2016). One can extend
this perturbative scheme by transforming Eq. (23) into yet another picture and
perform the RWA then in this frame. Doing so has for the parameters considered
not led to significant changes, which makes us believe that the above fast-
driving solution can for large $\Omega$ indeed be considered a benchmark
solution.

###  V.2 Master equation solutions

Since we are interested in the fast driving regime here, we consider kick
operator $U_{\rm kick}(t)\approx\exp\\{-\mathrm{i}\lambda/\Omega\sin(\Omega
t)\sigma^{x}\\}$ and system Floquet Hamiltonian
$\bar{H}\approx\frac{\Delta}{2}{\cal
J}_{0}\left(\frac{2\lambda}{\Omega}\right)\sigma^{z}$ also by applying the RWA
in the transformed frame. Considering the scaling of the Bessel functions at
small arguments, we therefore consistently only keep the lowest order terms
with $n\in\\{-1,0,+1\\}$ in Eq. (V.1). Then, the coupling operator in the
interaction picture can be approximately written as

$\displaystyle\mbox{\boldmath$A$}(t)$ $\displaystyle\approx{\cal
J}_{0}\left(\frac{2\lambda}{\Omega}\right)\sigma^{z}+2\sin(\Omega t){\cal
J}_{1}\left(\frac{2\lambda}{\Omega}\right)e^{+\mathrm{i}\bar{H}t}\sigma^{y}e^{-\mathrm{i}\bar{H}t}$
$\displaystyle\qquad+{\cal
O}\left\\{\frac{\lambda^{2}}{\Omega^{2}}\right\\}\,,$ (30)

which suffices to set up all master equations. We have numerically checked
convergence in the shown parameter regime by considering also the case
$n\in\\{-2,\ldots,+2\\}$ (not shown).

###  V.3 Comparison

In the limit of both fast driving and weak coupling we can compare the
benchmark solution with the perturbative PCG/DCG/BMS/BMU approaches. The
expectation value $\left<\sigma^{x}\right>$ is shown in Fig. 6.

Figure 6: (Color Online) Plot of the time-dependent expectation value of Pauli
matrix $\sigma^{x}$ comparing the no-driving solution (dashed green), the
analytical fast-driving solution (blue), PCG (orange), BMS (solid red), BMU
(dashed red) and DCG (purple cross). The gridlines mark multiples of the
driving period $T$. A higher order approximation for the driving yields
similar results (not shown). Parameters: $\Omega=25\;\Delta$,
$\lambda=\frac{1}{2}\;\Delta$, $\beta\,\Delta=1$, $\Gamma_{0}=0.05$,
$\omega_{c}=15\;\Delta$,
$\rho_{0}=1/2(\mbox{\boldmath$1$}+0.2\sigma^{x}+0.4\sigma^{z})$.

One can see that only the BMU (dashed red) and PCG (orange) solutions fail to
capture the correct dynamics, whereas BMS (solid red) and DCG (symbols)
solutions closely approximate it. It should be kept in mind however that for
such fast drivings, already the naive RWA in the original picture – yielding
the solution (8) with $\lambda=0$ (dashed green) – would yield acceptable
results for $\left<\sigma^{x}\right>$. This is different for the expectation
value $\left<\sigma^{z}\right>$, which we show in Fig. 7.

Figure 7: (Color Online) Analogous to Fig. 6, but for
$\left<\sigma^{z}\right>$. Whereas for the transient dynamics DCG, PCG, and
BMS solutions are close to the benchmark fast driving solution (inset),
significant differences arise for longer times (main plot).

For this expectation value, a naive RWA treatment would just predict a
constant evolution – which is indeed best captured by the long-term fast-
driving evolution and not by the perturbative approaches.

In both figures we observe that the DCG approach is closest to the FD
solution. It has however the disadvantage that the generator has to be re-
computed for any desired time. Still, all perturbative approaches fail to
correctly capture the long-term dynamics of the fast-driving solution (V.1),
which is not surprising as the latter is non-perturbative in the system-
reservoir coupling strength.

So far, we have considered the regime $\Omega\gg\omega_{c}\gg\Delta$, which
corresponds to $\tau_{\rm drv}\ll\tau_{\rm dec}\ll\tau_{\rm int}$. We see that
BMS and BMU solution fail to describe the long-term dynamics correctly, but
this is also the case for the DCG approach, which captures the short-term
dynamics by construction. In the regime $\Omega\gg\omega_{c}\approx\Delta$,
the BMS, BMU, and PCG solutions perform significantly worse, since then the
Markovian approximation is not valid anymore. In contrast, the DCG solution
still matches the fast-driving benchmark very well, as it is not subject to a
Markovian approximation. We show this in Fig. 8 for the expectation value
$\left<\sigma^{z}\right>$ (the expectation value of $\left<\sigma^{x}\right>$
looks similar, but there the driving has little effect anyways).

Figure 8: (Color Online) Analogous to Fig. 7, but for smaller cutoffs implying
a slower decay of the reservoir correlation function $\omega_{c}=\Delta$.
While DCG and fast-driving benchmark still agree perfectly, the Markovian
approaches (BMS, BMU, and PCG) all deviate strongly.

##  VI Conclusions

Our intention in this study was to find an improved time-independent Lindblad
(Markovian Breuer (2012)) master equation applicable to periodically driven
systems. For this, we analyzed driven qubit systems coupled to thermal
reservoirs with coarse-graining approaches. All methods used had the formal
advantage of being in LGKS form, which unconditionally preserves all density
matrix properties and leads to a stable numerical propagation. The negative
result is that the analyzed coarse-graining approach over one period of the
driving (PCG) did not match the expectations, i.e., it was inferior to the DCG
approach, even in the special case of pure dephasing.

On the positive side, we found that the BMS approach can yield quite reliable
results in the long term limit provided the secular approximation is performed
in a proper way. The DCG approach performed at least as good as the BMS
variant in the long-term limit and approaches the exact solution by
construction in the short-term limit. As the corresponding systematic
expansion converges in the weak-coupling regime (see App. A.2), we consider it
as the most accurate second-order perturbative solution among the schemes
tested. Its drawback is the computational cost, since for each time a suitable
dissipator needs to be calculated numerically. To avoid this, we note that by
using spectral densities with a simpler polynomial structure, the involved
integrals in (II) can be performed analytically. The fact that it is not
possible to capture the full dynamics faithfully with a single coarse-graining
time could be taken as a hint that the restriction to LGKS dynamics is too
severe. Using a Redfield approach however may be highly unstable for driven
systems even in the weak-coupling limit, such that one should rather aim at
deriving Kraus representations from first principles.

###### Acknowledgements.

R.H. acknowledges financial support by the Japanese Student Services
Organization JASSO. The authors gratefully acknowledge discussions with J.
Ablaßmayer, T. Becker, A. Eckardt, and A. Schnell.

##  Appendix A Coarse-graining generators

###  A.1 General coarse-graining derivation

To evaluate the r.h.s. of the defining equation (2), we first calculate the
time evolution operator in the interaction picture, where it follows the
equation

$\displaystyle\frac{d}{dt}\mbox{\boldmath$U$}(t,t_{0})=-\mathrm{i}\mbox{\boldmath$H$}_{I}(t,t_{0})\mbox{\boldmath$U$}(t,t_{0})\,,$
(31)

together with the initial condition
$\mbox{\boldmath$U$}(t_{0},t_{0})=\mbox{\boldmath$1$}$. Formally integrating
this equation and repeatedly inserting the solution into the r.h.s. yields the
representation

$\displaystyle\mbox{\boldmath$U$}(t,t_{0})$
$\displaystyle=\sum_{n=0}^{\infty}(-\mathrm{i})^{n}\int\limits_{t_{0}}^{t}dt_{1}\int\limits_{t_{0}}^{t_{1}}dt_{2}\ldots\int\limits_{t_{0}}^{t_{n-1}}dt_{n}$
$\displaystyle\qquad\times\mbox{\boldmath$H$}_{I}(t_{1})\mbox{\boldmath$H$}_{I}(t_{2})\ldots\mbox{\boldmath$H$}_{I}(t_{n})\,,$
(32)

where the nested integration bounds make the time-ordering explicit and the
$n=0$-term is understood as identity. Now, inserting a generic tensor product
decomposition of the interaction Hamiltonian

$\displaystyle\mbox{\boldmath$H$}_{I}(t,t_{0})$
$\displaystyle=\sum_{\alpha}\mbox{\boldmath$A$}_{\alpha}(t,t_{0})\otimes\mbox{\boldmath$B$}_{\alpha}(t,t_{0})$
$\displaystyle=\sum_{\alpha}U_{S}^{\dagger}(t,t_{0})A_{\alpha}U_{S}(t,t_{0})$
$\displaystyle\qquad\otimes
e^{+\mathrm{i}H_{B}(t-t_{0})}B_{\alpha}e^{-\mathrm{i}H_{B}(t-t_{0})}$ (33)

for system and bath coupling operators $A_{\alpha}$ and $B_{\alpha}$ 111Our
derivation also allows for driven system and bath coupling operators that are
time-dependent already in the Schrödinger picture, but we will drop this for
brevity., respectively, we see that the isolated time evolution operator of
the (undriven) reservoir can be given explicitly, whereas the explicit
calculation is required for the isolated time evolution operator of the
(driven) system $\mbox{\boldmath$U$}_{S}(t,t_{0})$.

Applied to the periodically-driven Schrödinger equation, Floquet theory
predicts that the time evolution operator can be expressed as a product of a
unitary and periodic kick operator (that has the period of the driving) and
the exponential of a Floquet Hamiltonian Eckardt and Anisimovas (2015).
Specifically, for $t_{0}=0$ one may write

$\displaystyle U_{S}(t)$ $\displaystyle=U_{\rm
kick}(t)e^{-\mathrm{i}\bar{H}t}\,,$ (34) $\displaystyle U_{\rm kick}(t)$
$\displaystyle=U_{\rm kick}(t+T)\,,\qquad U_{\rm kick}(nT)=\mathds{1}\,,$

with $n\in\mathds{Z}$, such that the stroboscopic evolution between full
periods is given by the Floquet Hamiltonian
$U_{S}(nT)=e^{-\mathrm{i}\bar{H}nT}$. It should be stressed here that even in
this simple case the above decomposition is not unique. For example,
introducing the eigenbasis of the Floquet Hamiltonian via
$\bar{H}\ket{\bar{a}}=\bar{E}_{a}\ket{\bar{a}}$, we can see that the time
evolution operator remains invariant under the simultaneous transformations
$\bar{H}\to\bar{H}+m\Omega\ket{\bar{a}}\bra{\bar{a}}$ and $U_{\rm kick}(t)\to
U_{\rm kick}(t)e^{+\mathrm{i}m\Omega t\ket{\bar{a}}\bra{\bar{a}}}$, which
maintains the periodicity of the kick operator. This gauge degree of freedom
allows one to deliberately choose the convention that all eigenenergies of the
Floquet Hamiltonian should lie in the first Brillouin zone
$\bar{E}_{a}\in[-\Omega/2,+\Omega/2)$. Hence, the system coupling operators in
the interaction picture
$\mbox{\boldmath$A$}_{\alpha}(t)=\mbox{\boldmath$A$}_{\alpha}(t,0)$ can be
expressed as:

$\displaystyle\boldsymbol{A}_{\alpha}(t)$
$\displaystyle=U_{S}^{\dagger}(t)A_{\alpha}U_{S}(t)=e^{\mathrm{i}\bar{H}t}U_{\rm
kick}^{\dagger}(t)A_{\alpha}U_{\rm kick}(t)e^{-\mathrm{i}\bar{H}t}$
$\displaystyle=\sum_{n=-\infty}^{+\infty}e^{\mathrm{i}\bar{H}t}\hat{A}_{\alpha}^{n}e^{\mathrm{i}n\Omega
t}e^{-\mathrm{i}\bar{H}t}\,,$ (35)

where due to the periodicity of the kick operator we have expanded $U_{\rm
kick}^{\dagger}(t)A_{\alpha}U_{\rm kick}(t)$ in a Fourier series

$\displaystyle\hat{A}_{\alpha}^{n}=\frac{\Omega}{2\pi}\int\limits^{T/2}_{-T/2}U_{\rm
kick}^{\dagger}(t)A_{\alpha}U_{\rm kick}(t)e^{-\mathrm{i}n\Omega t}dt,$
$\displaystyle T=\frac{2\pi}{\Omega}\,.$ (36)

Now, in the eigenbasis of the Floquet Hamiltonian we can also write this as

$\displaystyle\boldsymbol{A}_{\alpha}(t_{1})$
$\displaystyle=\sum_{ab}\sum_{n}\bra{\bar{a}}\hat{A}_{\alpha}^{n}\ket{\bar{b}}e^{\mathrm{i}(\bar{E}_{a}-\bar{E}_{b}+n\Omega)t_{1}}\ket{\bar{a}}\bra{\bar{b}}$
$\displaystyle\equiv\sum_{ab}\sum_{n}A_{\alpha,ab}^{+n}e^{+\mathrm{i}(\bar{E}_{a}-\bar{E}_{b}+n\Omega)t_{1}}L_{ab}\,,$
$\displaystyle\boldsymbol{A}_{\bar{\alpha}}(t_{2})$
$\displaystyle=\sum_{cd}\sum_{n^{\prime}}\bra{\bar{c}}\hat{A}_{\bar{\alpha}}^{n^{\prime}}\ket{\bar{d}}e^{\mathrm{i}(\bar{E}_{c}-\bar{E}_{d}+n^{\prime}\Omega)t_{2}}\ket{\bar{c}}\bra{\bar{d}}$
$\displaystyle\equiv\sum_{cd}\sum_{n^{\prime}}A_{\bar{\alpha},dc}^{-n^{\prime}}e^{-\mathrm{i}(\bar{E}_{c}-\bar{E}_{d}+n^{\prime}\Omega)t_{2}}L_{dc}\,,$
(37)

where $L_{ab}\equiv\ket{\bar{a}}\bra{\bar{b}}$ and we exchanged
$c\leftrightarrow d$ and $n^{\prime}\to-n^{\prime}$ in the last line.

The reservoir contributions to the full time evolution operator in the
interaction picture can be re-organized into reservoir correlation functions,
which for a reservoir at equilibrium depend only on the difference of time
arguments. Therefore, we can introduce their even and odd Fourier transforms
222 We note that higher-order correlation functions can be related to two-
point correlation functions for linear operators and thermal harmonic
reservoirs.

$\displaystyle C_{\alpha\bar{\alpha}}(t_{1}-t_{2})$
$\displaystyle=\frac{1}{2\pi}\int
d\omega\gamma_{\alpha\bar{\alpha}}(\omega)e^{-\mathrm{i}\omega(t_{1}-t_{2})}\,,$
$\displaystyle C_{\alpha\bar{\alpha}}(t_{1}-t_{2}){\rm sgn}(t_{1}-t_{2})$
$\displaystyle=\frac{1}{2\pi}\int
d\omega\sigma_{\alpha\bar{\alpha}}(\omega)e^{-\mathrm{i}\omega(t_{1}-t_{2})}\,.$
(38)

Inserting both the system and reservoir contributions above into (2),
expanding also its l.h.s. to first order and solving for ${\cal L}_{\tau}$,
the temporal integrals can be analytically performed with the help of

$\displaystyle f_{t_{0}}^{\tau}(\alpha,\beta,\omega)$
$\displaystyle\equiv\frac{1}{2\pi\tau}\int\limits_{t_{0}}^{t_{0}+\tau}\int\limits_{t_{0}}^{t_{0}+\tau}e^{-\mathrm{i}\omega(t_{1}-t_{2})}e^{\mathrm{i}\alpha
t_{1}}e^{-\mathrm{i}\beta t_{2}}dt_{1}dt_{2}$
$\displaystyle=\frac{\tau}{2\pi}e^{\mathrm{i}(\alpha-\beta)t_{0}}e^{\mathrm{i}(\alpha-\beta)\tau/2}{\rm
sinc}\left[(\alpha-\omega)\frac{\tau}{2}\right]$
$\displaystyle\qquad\times{\rm
sinc}\left[(\beta-\omega)\frac{\tau}{2}\right]\,.$ (39)

These definitions are used in the generators in Eqn. (II) in the main text and
in Eqns. (A.5), (A.6), and (A.7) below.

The function defined above has additional useful properties. For example, for
discrete $\alpha$ and $\beta$ one finds

$\displaystyle\lim_{\tau\to\infty}f_{t_{0}}^{\tau}(\alpha,\beta,\omega)=\delta_{\alpha\beta}\delta(\alpha-\omega)\,,$
(40)

which is used to derive the long-term limit of the coarse-graining generator
in App. A.5. Additionally, also for finite $\tau$, integrals involving this
function remain bounded as we show below.

###  A.2 Convergence of the expansion

Inserting the expansion for the time-evolution operator (A.1) into the
defining coarse-graining equation (2), we additionally obtain many temporal
integrals over higher-order correlation functions – some of them equipped with
Heaviside-$\Theta$ functions depending on time differences that enforce the
time-ordering. For linear couplings and harmonic reservoirs, only correlation
functions with an even number of operators remain, and furthermore, these can
then be expressed as products of two-point correlation functions. Thus,
defining Fourier transforms as before, we also obtain for higher-order terms
the nascent $\delta$-function (39) together with a function $g(\omega)$ that
contains Fourier transforms in analogy to
$\gamma_{\alpha\bar{\alpha}}(\omega)$ and
$\sigma_{\alpha\bar{\alpha}}(\omega)$ from Eq. (A.1). We can bound integrals
over products of this nascent $\delta$-function with other bounded functions
$g(\omega)$ by the essential supremum of the latter

$\displaystyle D$ $\displaystyle\equiv{\left|\int d\omega
g(\omega)f_{t_{0}}^{\tau}(\omega_{\alpha},\omega_{\beta},\omega)\right|}$
$\displaystyle\leq\int
d\omega{\left|g(\omega)f_{t_{0}}^{\tau}(\omega_{\alpha},\omega_{\beta},\omega)\right|}$
$\displaystyle\leq\left[{\rm
ess\;sup}_{\omega}{\left|g(\omega)\right|}\right]\int{\left|f_{t_{0}}^{\tau}(\omega_{\alpha},\omega_{\beta},\omega)\right|}d\omega$
$\displaystyle\leq\left[{\rm
ess\;sup}_{\omega}{\left|g(\omega)\right|}\right]\sqrt{\int\frac{\tau}{2\pi}{\rm
sinc}^{2}\left[(\omega_{\alpha}-\omega)\frac{\tau}{2}\right]d\omega}\times$
$\displaystyle\qquad\times\sqrt{\int\frac{\tau}{2\pi}{\rm
sinc}^{2}\left[(\omega_{\beta}-\omega)\frac{\tau}{2}\right]d\omega}$
$\displaystyle={\rm ess\;sup}_{\omega}{\left|g(\omega)\right|}\,.$ (41)

Here, in the second inequality we have used Hölders inequality for the
$1$-norm and the $\infty$-norm, and in the third inequality we have first
inserted (39) and then used Hölders inequality for the $2$-norm. Thus, if the
maximum of even and odd Fourier transforms of the correlation functions (that
enter as $g(\omega)$) is small, i.e., in the weak-coupling limit, we may
expect convergence of the series for all $\tau$.

For the bosonic reservoirs that we consider, a finite Fourier transform of the
reservoir correlation function requires at least ohmic spectral densities. If,
for example, the Fourier transform of the reservoir correlation function is
given by $\gamma(\omega)=\Gamma_{0}\omega
e^{-{\left|\omega\right|}/\omega_{c}}[1+1/(e^{\beta\omega}-1)]$, we find at
zero temperatures $\gamma(\omega)\leq\Gamma_{0}\omega_{\rm c}e^{-1}$ and at
high temperatures $\gamma(\omega)\leq\Gamma_{0}/\beta$.

###  A.3 Dynamical coarse-graining (DCG)

For dynamical coarse-graining, the time-dependent system density matrix is
expressed as

$\displaystyle\mbox{\boldmath$\rho_{S}$}(t)=\exp\left\\{\mbox{\boldmath${\cal
L}_{t}$}\cdot t\right\\}\mbox{\boldmath$\rho_{S}$}(t_{0})\,.$ (42)

This density matrix then formally evolves according to
$\frac{d}{dt}\mbox{\boldmath$\rho_{S}$}=\left[\left(\frac{d}{dt}e^{\mbox{\boldmath${\cal
L}_{t}$}\cdot t}\right)e^{-\mbox{\boldmath${\cal L}_{t}$}\cdot
t}\right]\mbox{\boldmath$\rho_{S}$}(t)$, where the term in square brackets can
be seen as a non-Markovian generator of the dynamical coarse-graining
evolution Rivas (2017). The full term in square brackets is not of LGKS form.
Nevertheless, the scheme preserves all density matrix properties, which
follows as ${\cal L}_{t}$ is of LGKS form for fixed $t$, such that the scheme
(42) is an interpolation through physically valid density matrices.

###  A.4 Initial-period coarse-graining (PCG)

Coarse-graining only over the initial period ($t_{0}=0$ and $\tau=T$), the
PCG-coarse-graining generator is obtained from (II) by inserting

$\displaystyle f_{0}^{T}(\omega_{1},\omega_{1},\omega)$
$\displaystyle=\frac{T}{2\pi}e^{\mathrm{i}(\omega_{1}-\omega_{2})T/2}$ (43)
$\displaystyle\qquad\times{\rm
sinc}\left[(\omega_{1}-\omega)\frac{T}{2}\right]{\rm
sinc}\left[(\omega_{2}-\omega)\frac{T}{2}\right]\,.$

The generator then still contains an integration over $\omega$, which –
depending on the form of the spectral density – may need to be solved
numerically.

###  A.5 Long-term coarse-graining

When $\tau\to\infty$, we get from Eq. (40) that the remaining integral in (II)
can be fully resolved, leading to

$\displaystyle{\bf H_{\rm LS}^{\infty}}$
$\displaystyle=\frac{1}{2\mathrm{i}}\sum_{\alpha\bar{\alpha}}\sum_{nn^{\prime}}\sum_{abcd}\sigma_{\alpha\bar{\alpha}}(\bar{E}_{a}-\bar{E}_{b}+n\Omega)$
$\displaystyle\qquad\times\delta_{\bar{E}_{a}-\bar{E}_{b}+n\Omega,\bar{E}_{c}-\bar{E}_{d}+n^{\prime}\Omega}A_{\alpha,ab}^{n}A_{\bar{\alpha},dc}^{-n^{\prime}}{L_{ab}L_{dc}}\,,$
$\displaystyle\mathcal{D}^{\infty}\boldsymbol{\rho_{S}}$
$\displaystyle=\sum_{\alpha\bar{\alpha}}\sum_{nn^{\prime}}\sum_{abcd}\gamma_{\alpha\bar{\alpha}}(\bar{E}_{a}-\bar{E}_{b}+n\Omega)$
$\displaystyle\qquad\times\delta_{\bar{E}_{a}-\bar{E}_{b}+n\Omega,\bar{E}_{c}-\bar{E}_{d}+n^{\prime}\Omega}A_{\alpha,ab}^{+n}A_{\bar{\alpha},dc}^{-n^{\prime}}$
$\displaystyle\qquad\times{\left[L_{dc}\boldsymbol{\rho_{S}}L_{ab}-\frac{1}{2}\left\\{L_{ab}L_{dc},\boldsymbol{\rho_{S}}\right\\}\right]}\,.$
(44)

Thus, the corresponding generator can be readily evaluated, and in particular
the long-term results from App. A.3 will coincide with the long-term results
from the equation above.

###  A.6 Floquet-BMS master equation

In Eq. (A.5), the evaluation of the resonance condition resulting from the
Kronecker-$\delta$

$\displaystyle\bar{E}_{a}-\bar{E}_{b}+n\Omega=\bar{E}_{c}-\bar{E}_{d}+n^{\prime}\Omega$
(45)

may require some care. The typical argument is that for fast driving, the
above resonance can only be met if separately $n^{\prime}=n$ and
$\bar{E}_{a}-\bar{E}_{b}=\bar{E}_{c}-\bar{E}_{d}$. However, the applicability
of this argument critically depends on the Floquet spectrum. If for a two-
level system we by chance have Floquet energies well within the first
Brillouin zone $\bar{E}_{a}\in\\{-\Omega/4,+\Omega/4\\}$, this generates
energy differences $\bar{E}_{a}-\bar{E}_{b}\in\\{-\Omega/2,0,+\Omega/2\\}$. In
this case, the above resonance condition can also be met with
$\bar{E}_{a}-\bar{E}_{b}=+\Omega/2$, $\bar{E}_{c}-\bar{E}_{d}=-\Omega/2$ and
$n^{\prime}=n+1$, demonstrating that in general, the long-term limit of
periodically driven coarse-graining master equation need not coincide with the
significantly simpler Floquet-BMS master equation that results from demanding
the resonance separately (i.e., by setting $n^{\prime}=n$ and
$\bar{E}_{a}-\bar{E}_{b}=\bar{E}_{c}-\bar{E}_{d}$)

$H_{\rm LS}^{\rm BMS}$
$\displaystyle=\frac{1}{2\mathrm{i}}\sum_{\alpha\bar{\alpha}}\sum_{n}\sum_{abcd}\sigma_{\alpha\bar{\alpha}}(\bar{E}_{a}-\bar{E}_{b}+n\Omega)$
$\displaystyle\qquad\times\delta_{\bar{E}_{a}-\bar{E}_{b},\bar{E}_{c}-\bar{E}_{d}}A_{\alpha,ab}^{n}A_{\bar{\alpha},dc}^{-n}{L_{ab}L_{dc}}\,,$
$\displaystyle\mathcal{D}^{\rm BMS}\boldsymbol{\rho_{S}}$
$\displaystyle=\sum_{\alpha\bar{\alpha}}\sum_{n}\sum_{abcd}\gamma_{\alpha\bar{\alpha}}(\bar{E}_{a}-\bar{E}_{b}+n\Omega)$
$\displaystyle\qquad\times\delta_{\bar{E}_{a}-\bar{E}_{b},\bar{E}_{c}-\bar{E}_{d}}A_{\alpha,ab}^{+n}A_{\bar{\alpha},dc}^{-n}$
$\displaystyle\qquad\times{\left[L_{dc}\boldsymbol{\rho_{S}}L_{ab}-\frac{1}{2}\left\\{L_{ab}L_{dc},\boldsymbol{\rho_{S}}\right\\}\right]}\,.$
(46)

###  A.7 The ultra-secular limit (BMU)

Even when (A.6) is applicable, we note that the remaining Kronecker-$\delta$
leaves many terms that are often neglected. For example, the above resonance
can always be trivially fulfilled with $\bar{E}_{a}=\bar{E}_{b}$ and
$\bar{E}_{c}=\bar{E}_{d}$, even for the undriven case. For many models, one
has that $\gamma_{\alpha\bar{\alpha}}(0)\to 0$, such that such terms would not
contribute anyways for undriven systems, where $\Omega\to 0$. But in general
they will have to be kept. For comparison we therefore also state the ultra-
secular approximation (BMU), where only terms with $a=c$ and $b=d$ are kept

$\displaystyle{\bf H_{\rm LS}^{\rm BMU}}$
$\displaystyle=\frac{1}{2\mathrm{i}}\sum_{\alpha\bar{\alpha}}\sum_{n}\sum_{ab}\sigma_{\alpha\bar{\alpha}}(\bar{E}_{a}-\bar{E}_{b}+n\Omega)$
$\displaystyle\quad\times
A_{\alpha,ab}^{n}A_{\bar{\alpha},ba}^{-n}{L_{ab}L_{ba}}\,,$
$\displaystyle\mathcal{D}^{\rm BMU}\boldsymbol{\rho_{S}}$
$\displaystyle=\sum_{\alpha\bar{\alpha}}\sum_{n}\sum_{ab}\gamma_{\alpha\bar{\alpha}}(\bar{E}_{a}-\bar{E}_{b}+n\Omega)A_{\alpha,ab}^{+n}A_{\bar{\alpha},ba}^{-n}$
$\displaystyle\quad\times\left[L_{ba}\boldsymbol{\rho_{S}}L_{ab}-\frac{1}{2}\left\\{L_{ab}L_{ba},\boldsymbol{\rho_{S}}\right\\}\right]\,.$
(47)

##  Appendix B Exact pure-dephasing solution

###  B.1 General solution

For a driven Hamiltonian of pure dephasing form (III), the exact dynamics can
be calculated by using a polaron (or Lang-Firsov Lang and Firsov (1963))
transform:

$\displaystyle U_{\rm
p}=e^{A\otimes\left(\frac{h_{k}}{\omega_{k}}b_{k}-\frac{h_{k}^{*}}{\omega_{k}}b_{k}^{\dagger}\right)}.$
(48)

Applying this transformation to system and bath operators yields

$\displaystyle U_{\rm p}^{\dagger}H_{S}(t)U_{\rm p}$
$\displaystyle=H_{S}(t)\,,\qquad U_{\rm p}^{\dagger}AU_{\rm p}=A\,,$ (49)
$\displaystyle U_{\rm p}^{\dagger}b_{k}U_{\rm p}$
$\displaystyle=b_{k}-\frac{h^{*}_{k}}{\omega_{k}}A\,,\qquad U_{\rm
p}^{\dagger}b_{k}^{\dagger}U_{\rm
p}=b_{k}^{\dagger}-\frac{h_{k}}{\omega_{k}}A\,,$

which leads to decoupled system and bath Hamiltonians

$\displaystyle U_{\rm p}^{\dagger}H(t)U_{\rm p}$
$\displaystyle=H_{S}(t)-\left(\sum_{k}\frac{|h_{k}|^{2}}{\omega_{k}}\right)A^{2}+\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}$
$\displaystyle\equiv\bar{H}_{S}(t)+\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k},$
(50)

with $\bar{H}_{S}(t)$ denoting an effective system Hamiltonian. Hence, the
time evolution operator in the polaron picture $\tilde{U}(t)$ is given by the
system and bath evolution separately

$\displaystyle\tilde{U}(t)=\mathcal{T}\left\\{e^{-\mathrm{i}\int_{0}^{t}\bar{H}_{S}(t^{\prime})dt^{\prime}}\right\\}e^{-\mathrm{i}\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}t}\,,$
(51)

where $\mathcal{T}$ denotes the time-ordering operator of the system. For a
general unitary transformation $V(t)$ into another frame via
$\ket{\Psi(t)}=V(t)\ket{\tilde{\Psi}(t)}$, the time evolution operator in the
original frame $U(t)$ can be expressed in terms of the time evolution operator
in the new frame $\tilde{U}(t)$ as $U(t)=V(t)\tilde{U}(t)V^{\dagger}(0)$.
Thus, for the time-independent polaron transform this relation yields

$\displaystyle U(t)=U_{\rm p}\tilde{U}(t)U_{\rm p}^{\dagger}.$ (52)

From this, the expectation value of any system observable $O_{S}$,

$\displaystyle\left<O_{S}\right>$ $\displaystyle={\rm
Tr}\left\\{(O_{S}\otimes\mbox{\boldmath$1$}_{B})\rho(t)\right\\}$
$\displaystyle={\rm
Tr}\left\\{U^{\dagger}(t)(O_{S}\otimes\mbox{\boldmath$1$}_{B})U(t)\rho^{0}\right\\},$
(53)

can be calculated as:

$\displaystyle\left<O_{S}\right>={\rm Tr}\left\\{U_{\rm
p}\tilde{U}^{\dagger}(t)U_{\rm
p}^{\dagger}(O_{S}\otimes\mbox{\boldmath$1$}_{B})U_{\rm p}\tilde{U}(t)U_{\rm
p}^{\dagger}\rho^{0}\right\\}\,.$ (54)

###  B.2 Two-level system

For the example of a two level system,

$\displaystyle H_{S}(t)=\frac{\Delta}{2}\sigma^{z}{+\lambda\cos\Omega
t}\,,\qquad A=\sigma^{z},$ (55)

the populations
$\rho_{S,00}(t)=\frac{1}{2}\langle\mathds{1}+\sigma^{z}\rangle$ and
$\rho_{S,11}(t)=\frac{1}{2}\langle\mathds{1}-\sigma^{z}\rangle$ are constant
as $\mathds{1}$ and $\sigma^{z}$ commute with $U_{\rm p}$ and $\tilde{U}(t)$,
and we obtain (7) in the main text.

For the coherences $\rho_{01}=\left<\sigma^{-}\right>$ and
$\rho_{10}=\left<\sigma^{+}\right>$, the calculations are more involved. Using
that

$\displaystyle U_{\rm p}^{\dagger}\sigma^{\pm}U_{\rm p}=e^{\pm
2\sum_{k}\left(\frac{h_{k}^{*}}{\omega_{k}}b_{k}^{\dagger}-\frac{h_{k}}{\omega_{k}}b_{k}\right)}\sigma^{\pm},$
(56)

and organising bath and system parts together, one gets

$\displaystyle\left<\sigma^{\pm}\right>$ $\displaystyle={\rm
Tr}\Big{\\{}U_{\rm p}U_{B}^{\dagger}(t)e^{\pm
2\sum_{k}\left(\frac{h_{k}^{*}}{\omega_{k}}b_{k}^{\dagger}-\frac{h_{k}}{\omega_{k}}b_{k}\right)}U_{B}(t)U_{S}^{\dagger}(t)\sigma^{\pm}$
$\displaystyle\qquad\times U_{S}(t)U_{\rm p}^{\dagger}\rho^{0}\Big{\\}},$ (57)

with $U_{B}(t)=e^{-\mathrm{i}\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}t}$ and
$U_{S}(t)=e^{-\mathrm{i}\left(\frac{\Delta}{2}t+\frac{\lambda}{\Omega}\sin(\Omega
t)\right)\sigma^{z}}$. Now, one can use that

$\displaystyle e^{\mathrm{i}\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}t}e^{\pm
2\sum_{k}\left(\frac{h_{k}^{*}}{\omega_{k}}b_{k}^{\dagger}-\frac{h_{k}}{\omega_{k}}b_{k}\right)}e^{-\mathrm{i}\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}t}$
$\displaystyle\qquad=e^{\pm
2\sum_{k}\left(\frac{h_{k}^{*}}{\omega_{k}}b_{k}^{\dagger}e^{\mathrm{i}\omega_{k}t}-\frac{h_{k}}{\omega_{k}}b_{k}e^{-\mathrm{i}\omega_{k}t}\right)}\,,$
$\displaystyle
e^{\mathrm{i}x\sigma^{z}}\sigma^{\pm}e^{-\mathrm{i}x\sigma^{z}}=e^{\pm
2\mathrm{i}x}\sigma^{\pm}\,,$ (58)

which leads by inserting the identity $U_{\rm p}^{\dagger}U_{\rm
p}=\mathds{1}$ to:

$\displaystyle\left<\sigma^{\pm}\right>$ $\displaystyle=e^{\pm
2\mathrm{i}\left(\frac{\Delta}{2}t+\frac{\lambda}{\Omega}\sin(\Omega
t)\right)}$ (59) $\displaystyle\quad\times{\rm Tr}\left\\{U_{\rm p}e^{\pm
2\sum_{k}\left(\frac{h_{k}^{*}}{\omega_{k}}b_{k}^{\dagger}e^{\mathrm{i}\omega_{k}t}-{\rm
h.c.}\right)}U_{\rm p}^{\dagger}U_{\rm p}\sigma^{\pm}U_{\rm
p}^{\dagger}\rho^{0}\right\\}\,.$

Next, the polaron transformation is applied to the system and bath parts
separately:

$\displaystyle U_{\rm p}e^{\pm
2\sum\limits_{k}\frac{h_{k}^{*}}{\omega_{k}}b_{k}^{\dagger}e^{\mathrm{i}\omega_{k}t}-{\rm
h.c.}}U_{\rm p}^{\dagger}$ $\displaystyle=e^{\pm
2\sum\limits_{k}\frac{h_{k}^{*}}{\omega_{k}}\left(b_{k}^{\dagger}+\frac{h_{k}}{\omega_{k}}\sigma^{z}\right)e^{\mathrm{i}\omega_{k}t}-{\rm
h.c.}}\,,$ $\displaystyle U_{\rm p}\sigma^{\pm}U_{\rm p}^{\dagger}$
$\displaystyle=e^{\mp
2\sum_{k}\left(\frac{h_{k}^{*}}{\omega_{k}}b_{k}^{\dagger}-\frac{h_{k}}{\omega_{k}}b_{k}\right)}\sigma^{\pm}.$
(60)

Separating system and bath parts and inserting the initial condition
$\rho^{0}=\rho^{0}_{S}\otimes\bar{\rho}_{B}$ yields

$\displaystyle\left<\sigma^{\pm}\right>$ $\displaystyle=e^{\pm
2\mathrm{i}\left(\frac{\Delta}{2}t+\frac{\lambda}{\Omega}\sin(\Omega
t)\right)}$ $\displaystyle\qquad\times{\rm Tr}\left\\{e^{\pm
4\mathrm{i}\sum_{k}\frac{|h_{k}|^{2}}{\omega_{k}^{2}}\sin(\omega_{k}t)\sigma^{z}}\sigma^{\pm}\rho_{S}^{0}\right\\}B_{\pm}(t)\,,$
$\displaystyle B_{\pm}(t)$ $\displaystyle={\rm Tr}\Big{\\{}e^{\pm
2\sum_{k}\left(\frac{h_{k}^{*}}{\omega_{k}}b_{k}^{\dagger}e^{\mathrm{i}\omega_{k}t}-\frac{h_{k}}{\omega_{k}}b_{k}e^{-\mathrm{i}\omega_{k}t}\right)}$
$\displaystyle\qquad\times e^{\mp
2\sum_{k}\left(\frac{h_{k}^{*}}{\omega_{k}}b_{k}^{\dagger}-\frac{h_{k}}{\omega_{k}}b_{k}\right)}\bar{\rho}_{B}\Big{\\}}\,.$
(61)

The trace over the bath parts $B_{\pm}(t)$ gives:

$\displaystyle B_{\pm}(t)$ $\displaystyle=\exp\left\\{\mp
4\mathrm{i}\sum_{k}\frac{|h_{k}|^{2}}{\omega_{k}^{2}}\sin(\omega_{k}t)\right\\}$
(62)
$\displaystyle\times\exp\left\\{-4\sum_{k}\frac{|h_{k}|^{2}}{\omega_{k}^{2}}[1-\cos(\omega_{k}t)][1+2n_{B}(\omega_{k})]\right\\}\,.$

Here, $n_{B}(\omega_{k})=\frac{1}{e^{\beta\omega_{k}}-1}$ denotes the Bose
distribution. Plugging the expression for $B_{\pm}(t)$ into Eq. (61) and using

$\displaystyle{\rm Tr}\left\\{e^{\pm
4\mathrm{i}\sum_{k}\frac{|h_{k}|^{2}}{\omega_{k}^{2}}\sin(\omega_{k}t)\sigma^{z}}\sigma^{\pm}\rho_{S}^{0}\right\\}$
$\displaystyle=e^{\pm
4i\sum_{k}\frac{|h_{k}|^{2}}{\omega_{k}^{2}}\sin(\omega_{k}t)}$
$\displaystyle\qquad\times{\rm Tr}\left\\{\sigma^{\pm}\rho_{S}^{0}\right\\},$
(63)

yields – after inserting the spectral coupling density
$\Gamma(\omega)=2\pi\sum_{k}|h_{k}|^{2}\delta(\omega-\omega_{k})$ – for the
expectation value of $\left<\sigma^{+}\right>=\rho_{10}$ Eq. (8) in the main
text.

##  Appendix C Analytic fast-driving approximation

We consider the simplified form of Eq. (V) with $H_{S}^{0}=\frac{\Delta}{2}A$.
First, we move into an interaction picture with respect to the driving by
performing the transformation
$U_{1}(t)=e^{-\mathrm{i}\frac{\lambda}{\Omega}\sin(\Omega t)C}$, which yields
after applying the rotating wave approximation

$\displaystyle\tilde{H}_{\rm RWA}$
$\displaystyle=\tilde{A}\left(\frac{\Delta}{2}+\sum_{k}(h_{k}b_{k}+h_{k}^{*}b_{k}^{\dagger})\right)+\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}\,,$
$\displaystyle\tilde{A}$
$\displaystyle=\frac{\Omega}{2\pi}\int_{0}^{\frac{2\pi}{\Omega}}U_{1}^{\dagger}(t)AU_{1}(t)dt\,.$
(64)

This is the Hamiltonian of a simple pure dephasing system, see Eq. (III), with
time independent system Hamiltonian.

Applying the polaron transform (48), system and bath part of the Hamiltonian
can be decoupled and (50) becomes

$\displaystyle H_{\rm p}$ $\displaystyle=U_{\rm p}^{\dagger}\tilde{H}_{\rm
RWA}U_{\rm p}$
$\displaystyle=\frac{\Delta}{2}\tilde{A}-\sum_{k}\tfrac{|h_{k}|^{2}}{\omega_{k}}\tilde{A}^{2}+\sum_{k}h_{k}b_{k}^{\dagger}b_{k}\,,$
(65)

which yields for the time evolution operator in this picture:

$\displaystyle\tilde{U}(t)\approx
e^{-\mathrm{i}\left(\frac{\Delta}{2}\tilde{A}-\sum_{k}\tfrac{|h_{k}|^{2}}{\omega_{k}}\tilde{A}^{2}\right)t}e^{-\mathrm{i}\sum_{k}h_{k}b_{k}^{\dagger}b_{k}t}.$
(66)

In an anologous way as for the pure dephasing system, we can write the
expectation value of an observable as:

$\displaystyle\langle O_{S}\rangle={\rm Tr}\left\\{U_{\rm
p}\tilde{U}^{\dagger}(t)U_{\rm
p}^{\dagger}U_{1}^{\dagger}(O_{S}\otimes\mbox{\boldmath$1$}_{B})U_{1}U_{\rm
p}\tilde{U}(t)U_{\rm p}^{\dagger}\rho^{0}\right\\}.$ (67)

For a two level system Eq. (26), we get

$\displaystyle\tilde{A}=\mathcal{J}_{0}\left(\tfrac{2\lambda}{\Omega}\right)\sigma^{z},$
(68)

with $\mathcal{J}_{n}(x)$ denoting the Bessel function of the first kind,
which is defined as the solution of the differential equation
$x^{2}\mathcal{J}_{n}^{\prime\prime}(x)+x\mathcal{J}_{n}^{\prime}(x)-(z^{2}+n^{2})\mathcal{J}_{n}(x)=0$,
eventually leading to (V.1) in the main text.

## References

  * Floquet (1883) G. Floquet, Annales scientifiques de l’École Normale Supérieure 12, 47 (1883).
  * Shirley (1965) J. H. Shirley, Phys. Rev. 138, B979 (1965).
  * Grifoni and Hänggi (1998) M. Grifoni and P. Hänggi, Physics Reports 304, 229 (1998).
  * Eckardt and Anisimovas (2015) A. Eckardt and E. Anisimovas, New Journal of Physics 17, 093039 (2015).
  * Bukov  _et al._ (2015) M. Bukov, L. D’Alessio, and A. Polkovnikov, Advances in Physics 64, 139 (2015).
  * Weiss (1993) U. Weiss,  _Quantum Dissipative Systems_ , Series of Modern Condensed Matter Physics, Vol. 2 (World Scientific, Singapore, 1993).
  * Breuer and Petruccione (2002) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2002).
  * Schlosshauer (2007) M. Schlosshauer,  _Decoherence and the Quantum-To-Classical Transition_ (Springer-Verlag Berlin Heidelberg, 2007).
  * Chu (1991) S. Chu, Science 253, 861 (1991).
  * Rice and Zhao (2000) S. Rice and M. Zhao,  _Optical Control of Molecular Dynamics_ (Wiley, 2000).
  * Aidelsburger  _et al._ (2013) M. Aidelsburger, M. Atala, S. Nascimbène, S. Trotzky, Y.-A. Chen, and I. Bloch, Applied Physics B 113, 1–11 (2013).
  * Atala  _et al._ (2014) M. Atala, M. Aidelsburger, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Nature Physics 10, 588–593 (2014).
  * Jotzu  _et al._ (2014) G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, and T. Esslinger, Nature 515, 237–240 (2014).
  * Lindblad (1976) G. Lindblad, Communications in Mathematical Physics 48, 119 (1976).
  * Gorini  _et al._ (1976) V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Journal of Mathematical Physics 17, 821 (1976).
  * Kohler  _et al._ (1997) S. Kohler, T. Dittrich, and P. Hänggi, Physical Review E 55, 300–313 (1997).
  * Gelbwaser-Klimovsky  _et al._ (2013) D. Gelbwaser-Klimovsky, R. Alicki, and G. Kurizki, Physical Review E 87, 012140 (2013).
  * Szczygielski  _et al._ (2013) K. Szczygielski, D. Gelbwaser-Klimovsky, and R. Alicki, Physical Review E 87, 012120 (2013).
  * Szczygielski (2014) K. Szczygielski, Journal of Mathematical Physics 55, 083506 (2014).
  * Shirai  _et al._ (2015) T. Shirai, T. Mori, and S. Miyashita, Physical Review E 91 (2015), 10.1103/physreve.91.030101.
  * Shirai  _et al._ (2016) T. Shirai, J. Thingna, T. Mori, S. Denisov, P. Hänggi, and S. Miyashita, New Journal of Physics 18, 053008 (2016).
  * Cuetara  _et al._ (2015) G. B. Cuetara, A. Engel, and M. Esposito, New Journal of Physics 17, 055002 (2015).
  * Restrepo  _et al._ (2019) S. Restrepo, S. Böhling, J. Cerrillo, and G. Schaller, Phys. Rev. B 100, 035109 (2019).
  * Luo  _et al._ (2020) J. Luo, Y. Yan, H. Wang, J. Hu, X.-L. He, and G. Schaller, Phys. Rev. B 101, 125410 (2020).
  * Hone  _et al._ (2009) D. W. Hone, R. Ketzmerick, and W. Kohn, Phys. Rev. E 79, 051129 (2009).
  * Schnell  _et al._ (2020) A. Schnell, A. Eckardt, and S. Denisov, Phys. Rev. B 101, 100301 (2020).
  * Hartmann and Strunz (2020) R. Hartmann and W. T. Strunz, Phys. Rev. A 101, 012103 (2020).
  * Whitney (2008) R. S. Whitney, Journal of Physics A: Mathematical and General 41, 175304 (2008).
  * Kiršanskas  _et al._ (2018) G. Kiršanskas, M. Franckié, and A. Wacker, Phys. Rev. B 97, 035432 (2018).
  * Ptaszyński and Esposito (2019) K. Ptaszyński and M. Esposito, Phys. Rev. Lett. 122, 150603 (2019).
  * Farina and Giovannetti (2019) D. Farina and V. Giovannetti, Phys. Rev. A 100, 012107 (2019).
  * Cattaneo  _et al._ (2020) M. Cattaneo, G. L. Giorgi, S. Maniscalco, and R. Zambrini, Phys. Rev. A 101, 042108 (2020).
  * Schultz and von Oppen (2009) M. G. Schultz and F. von Oppen, Physical Review B 80, 033302 (2009).
  * Schultz (2010) M. G. Schultz, Physical Review B 82, 155408 (2010).
  * Teretenkov (2021) A. E. Teretenkov, Journal of Physics A: Mathematical and Theoretical 54, 265302 (2021).
  * Trushechkin (2021) A. Trushechkin, Phys. Rev. A 103, 062226 (2021).
  * Vacchini (2000) B. Vacchini, Physical Review Letters 84, 1374–1377 (2000).
  * Lidar  _et al._ (2001) D. A. Lidar, Z. Bihary, and K. Whaley, Chemical Physics 268, 35–53 (2001).
  * Schaller and Brandes (2008) G. Schaller and T. Brandes, Physical Review A 78 (2008), 10.1103/physreva.78.022106.
  * Schaller  _et al._ (2009) G. Schaller, P. Zedler, and T. Brandes, Physical Review A 79, 032110 (2009).
  * Benatti  _et al._ (2009) F. Benatti, R. Floreanini, and U. Marzolino, epl 88 (2009), 10.1209/0295-5075/88/20011.
  * Benatti  _et al._ (2010) F. Benatti, R. Floreanini, and U. Marzolino, Phys. Rev. A 81, 012105 (2010).
  * Royer (1996) A. Royer, Physical Review Letters 77 (1996), 10.1103/PhysRevLett.77.3272.
  * Majenz  _et al._ (2013) C. Majenz, T. Albash, H.-P. Breuer, and D. A. Lidar, Physical Review A 88, 012103 (2013).
  * Rivas (2017) A. Rivas, Physical Review A 95, 042104 (2017).
  * Rivas (2019) A. Rivas, Entropy (2019), 10.3390/e21080725.
  * Mozgunov and Lidar (2020) E. Mozgunov and D. Lidar, Quantum 4, 227 (2020).
  * Kraus (1971) K. Kraus, Annals of Physics 64, 311 (1971).
  * Schaller and Ablaßmayer (2020) G. Schaller and J. Ablaßmayer, Entropy 22, 525 (2020).
  * Elouard  _et al._ (2020) C. Elouard, D. Herrera-Martí, M. Esposito, and A. Auffèves, New Journal of Physics 22, 103039 (2020).
  * Garraway (1997) B. M. Garraway, Phys. Rev. A 55, 2290 (1997).
  * Kubo (1957) R. Kubo, Journal of the Physical Society of Japan 12, 570 (1957).
  * Martin and Schwinger (1959) P. C. Martin and J. Schwinger, Physical Review 115, 1342 (1959).
  * Kossakowski  _et al._ (1977) A. Kossakowski, A. Frigerio, V. Gorini, and M. Verri, Communications in Mathematical Physics 57, 97 (1977).
  * Abramowitz and Stegun (1970) M. Abramowitz and I. A. Stegun, eds.,  _Handbook of Mathematical Functions_ (National Bureau of Standards, 1970).
  * Leggett  _et al._ (1987) A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Reviews of Modern Physics 59, 1 (1987).
  * Restrepo  _et al._ (2016) S. Restrepo, J. Cerrillo, V. M. Bastidas, D. G. Angelakis, and T. Brandes, Phys. Rev. Lett. 117, 250401 (2016).
  * Breuer (2012) H.-P. Breuer, Journal of Physics B: Atomic, Molecular and Optical Physics 45, 154001 (2012).
  * Note (1) Our derivation also allows for driven system and bath coupling operators that are time-dependent already in the Schrödinger picture, but we will drop this for brevity.
  * Note (2) We note that higher-order correlation functions can be related to two-point correlation functions for linear operators and thermal harmonic reservoirs.
  * Lang and Firsov (1963) I. G. Lang and Y. A. Firsov, Soviet Physics JETP 16, 1301 (1963).