# Enhancing energy transfer in quantum systems via periodic driving: Floquet
master equations

Z. Nafari Qaleh  Department of Physics, Sharif University of Technology,
Tehran 14588, Iran  A. T. Rezakhani  Department of Physics, Sharif University
of Technology, Tehran 14588, Iran

###### Abstract

We provide a comprehensive study of the energy transfer phenomenon—populating
a given energy level—in $3$\- and $4$-level quantum systems coupled to two
thermal baths. In particular, we examine the effects of an external periodic
driving and the coherence induced by the baths on the efficiency of the energy
transfer. We consider the Floquet-Lindblad and the Floquet-Redfield scenarios,
which both are in the Born-Markov, weak-coupling regime but differ in the
treatment of the secular approximation, and for the latter we develop an
appropriate Floquet-type master equation by employing a partial secular
approximation. Throughout the whole analysis we keep Lamb-shift corrections in
the master equations. We observe that, especially in the Floquet-Redfield
scenario, the driving field can enhance the energy transfer efficiency
compared to the nondriven scenario. In addition, unlike degenerate systems
where Lamb-shift corrections do not contribute significantly on the energy
transfer, in the Redfield and the Floquet-Redfield scenarios these corrections
have nonnegligible effects.

##  I Introduction

The process of energy transfer is a fundamental subject in physics, chemistry,
and biology, which involves energy exchange among parts of physical, chemical,
and biological systems. In particular, it is interesting to study such
processes in quantum regimes [1, 2]. There exist a vast literature of
pertinent studies in various systems such as crystals with impurities [3],
quantum-dot nanostructures [4], polymer chains [5, 6], and light-harvesting
complexes in photosynthesis [7, 8, 9, 10]. Increasing efficiency (or power) of
energy transfer processes has been one of the main objectives in these
investigations [11, 12]. Energy transfer usually occurs in two forms: (i) in
real space, where energy in the form of excitation (or exciton) is to move
from one particular “site” of a lattice to another target site [13, 14, 15,
16]; and (ii) in energy space, where energy is to move from an eigenstate of
the system Hamiltonian to populate another target eigenstate or energy level
[17, 18]. The latter event in an open quantum system in the presence of baths
is of our interest in this paper.

In real systems, interaction between the system and its environment (or bath)
is inevitable. This results in dissipative dynamics of the system in the sense
that the system evolution is given by nonunitary master equations, usually
after applying relevant approximations such as the Born-Markov, weak-coupling,
and secular approximations [19, 20, 21, 22]. Under the first two
approximations, the Redfield master equation is obtained; the three
approximations yield the Lindblad master equation [21]. The Redfield equation
may not necessarily guarantee the complete-positivity of the system density
matrix (nonpositivity issue), and evolution of populations (diagonal terms of
the density matrix) through this equation may involve dependence on coherence
(offdiagonal) terms. Applying the secular approximation (if holds) retrieves
the complete-positivity in the resulting Lindblad equation. If the Hamiltonian
of the system is nondegenerate, the Lindblad equation yields an evolution for
populations fully independently from coherence terms [23, 24, 25]. Although
there are some arguments against the use of the Redfield equation [26],
several recent studies have shown the utility of this equation in some
applications [27, 28]. In fact, the main nonpositivity issue with the Redfield
equation may be partially alleviated without need to apply the full secular
approximation, e.g., by applying a suitable “partial” secular approximation
[29, 28, 27]. This way the Redfield equation does not generate unphysical
states yet it maintains its appeal to couple populations and coherences and
gives physically relevant results.

In this paper, we consider a model which includes a multilevel system
$\mathsf{S}$ in contact with two thermal baths at two different temperatures.
The hot bath (e.g., sunlight) is set to excite the system to its highest
energy level(s) and then the exciton moves through the energy space in the
presence of another cold bath. We work in the local master equation
approximation where the couplings are so weak that the baths can be assumed to
act only locally on the part of the system they directly interact with [30,
31], and thus effectively do not see each other. We consider four scenarios–as
in Fig. 1—and their impacts on the energy transfer efficiency: (i) Lindblad
scenario: assuming the Born-Markov and secular approximations, with no applied
driving field on the system, (ii) Floquet-Lindblad scenario: the Born-Markov
and secular approximations and driving the system by a time-periodic applied
field, (iii) Redfield scenario: the Born-Markov approximation, no secular
approximation, and no applied field on the system, and (iv) Floquet-Redfield
scenario: the Born-Markov approximation, a partial secular approximation, and
driving the system by a time-periodic applied field.

Figure 1: Four scenarios of the energy transfer. (a) Lindblad scenario: time-
independent system by applying the Born-Markov and secular approximations.
Here the standard Lindblad equation [Eq. (6)] applies. (b) Redfield scenario:
time-independent system and performing only the Born-Markov approximation.
Here the Redfield equation [Eq. (27)] applies. (c) Floquet-Lindblad scenario:
time-dependent system with the Born-Markov and secular approximations. Here
the Floquet-Lindblad equation [Eq. (11)] applies. (d) Floquet-Redfield
scenario: time-dependent system with the Born-Markov approximation and a
partial secular approximation. Here the Floquet-Redfield equation [Eq. (34)]
applies. In time-dependent cases the system is driven by a time-periodic
applied field.

Each of these scenarios require its own dynamical equation for the system. In
particular, in the time-dependent driven cases since the applied field is
periodic, one needs to use Floquet theory. Scenarios (i), (ii), and (iii) can
be described by, respectively, Lindblad equation, Floquet-Lindblad equation,
and Redfield equation [32, 33]. However, for scenario (iv) we develop an
appropriate Floquet-type Redfield equation where we employ a particular
partial secular approximation (which is different from that of Ref. [29]).

Our objective is to populate a particular energy level of the system as the
target state (denoted by $|b\rangle$). Performance of this task can be
quantified by the time averaging of the population of the target state [34,
35, 36, 37],

$\displaystyle\eta(t_{\mathrm{f}})=(1/t_{\mathrm{f}})\textstyle{\int_{0}^{t_{\mathrm{f}}}}ds\,\langle
b|\varrho_{\mathsf{S}}(s)|b\rangle,$ (1)

where $t_{\mathrm{f}}$ is the final time and $\varrho_{\mathsf{S}}(s)$ denotes
the state of the system at time $0\leqslant s\leqslant t_{\mathrm{f}}$. Our
results show that driving the system by a periodic external field leads to
increasing the population of the the target level relative to the case where
no external field is applied. In the Redfield and Floquet-Redfield master
equations we include the Lamb-shift terms, which are computed by appropriate
QED considerations. According to our analysis, it is seen that at least for
degenerate systems the Lamb-shift terms do not have a significant impact on
the dynamics. In addition, we observe that in the presence of a periodic
external field, generating coherence in the system through the dynamics can
enhance the energy transfer efficiency (in the Floquet-Lindblad and Floquet-
Redfield scenarios).

The structure of this paper is as follows. In Sec. II we investigate the
Lindblad and Floquet-Lindblad scenarios and their associated dynamical
equations. In addition, we consider the case of a $3$-level quantum system and
calculate the energy transfer efficiency for each of the scenarios. In Sec.
III we move to the Redfield and Floquet-Redfield scenarios. In particular, we
delineate how one can derive an appropriate dynamical equation which holds for
a driven quantum system under a time-periodic external field under the Born-
Markov condition. To do so, we introduce a specific partial secular
approximation. Next, we calculate the energy transfer efficiency for these
scenarios for a $4$-level system, which also enables us to see how coherence
may also affect the energy transfer efficiency. We conclude the paper with a
summary in Sec. IV. Several appendixes are also included to elaborate further
details and steps of derivations.

##  II Lindblad and Floquet-Lindblad scenarios

Here we sketch derivation of the quantum master equation for the state of an
open quantum quantum system which is driven by a periodic external field.
Before presenting the Floquet-Lindblad equation we set up the problem and
briefly review the standard Lindblad master equation.

###  II.1 Lindblad scenario

Figure 2: A $3$-level system ($\\{|0\rangle,|1\rangle,|b\rangle\\}$) coupled
to two baths and driven by an external periodic field $V(t)$. In (a) the
external field couples $|0\rangle$ and $|b\rangle$, while in (b) it couples
$|1\rangle$ and $|b\rangle$.

Consider a nondriven, closed bipartite quantum system with the total
Hamiltonian

$\displaystyle
H_{\mathsf{tot}}=H_{\mathsf{S}}+H_{\mathsf{B}}+H_{\mathsf{int}},$ (2)

where
$H_{\mathsf{S}}=\sum_{\epsilon}\epsilon|\epsilon\rangle\langle\epsilon|$,
$H_{\mathsf{B}}$, and $H_{\mathsf{int}}=\sum_{\alpha}S_{\alpha}\otimes
B_{\alpha}$ are, respectively, the system, bath, and interaction Hamiltonians.
It is more convenient to derive the master equation in the interaction picture
which is defined by

$\displaystyle\bm{\varrho}_{\mathsf{tot}}(t)=U^{\dagger}_{\mathsf{S}}(t,0)\otimes
U^{\dagger}_{\mathsf{B}}(t,0)\varrho_{\mathsf{tot}}(t)U_{\mathsf{S}}(t,0)\otimes
U_{\mathsf{B}}(t,0),$ (3)
$\displaystyle\bm{H}_{\mathsf{int}}(t)=U^{\dagger}_{\mathsf{S}}(t,0)\otimes
U^{\dagger}_{\mathsf{B}}(t,0)H_{\mathsf{int}}U_{\mathsf{S}}(t,0)\otimes
U_{\mathsf{B}}(t,0),$ (4)

where $U_{\mathsf{S}}(t,0)=e^{-itH_{\mathsf{S}}/\hbar}$ and
$U_{\mathsf{B}}(t,0)=e^{-itH_{\mathsf{B}}/\hbar}$. We use boldface letters for
the interaction picture throughout this paper. The dynamics of the total
system in this picture is given by

$\displaystyle\dfrac{d\bm{\varrho}_{\mathsf{tot}}(t)}{dt}=-\dfrac{i}{\hbar}\left[\bm{H}_{\mathsf{int}}(t),\bm{\varrho}_{\mathsf{tot}}(t)\right].$
(5)

It has been argued that under the Born-Markov and secular approximations the
dynamics of the open quantum system is given by the Lindblad quantum master
equation [21],

$\displaystyle\dfrac{d\bm{\varrho}_{\mathsf{S}}(t)}{dt}=-\dfrac{i}{\hbar}\left[H_{\mathrm{lamb}},\bm{\varrho}_{\mathsf{S}}(t)\right]+\mathpzc{D}\left[\bm{\varrho}_{\mathsf{S}}(t)\right],$
(6)

where

$\displaystyle
H_{\mathrm{lamb}}=(1/\hbar)\textstyle{\sum_{\alpha,\alpha^{\prime}}}\sum_{{\omega}}\xi_{\alpha\alpha^{\prime}}(\omega)S_{\alpha}^{\dagger}(\omega)S_{\alpha^{\prime}}(\omega),$
(7)
$\displaystyle\mathpzc{D}[\bm{\varrho}_{\mathsf{S}}(t)]=(1/\hbar^{2})\textstyle{\sum_{\alpha,\alpha^{\prime}}}\sum_{{\omega}}\gamma_{\alpha\alpha^{\prime}}(\omega)$
$\displaystyle\,\times\big{[}S_{\alpha^{\prime}}(\omega)\bm{\varrho}_{\mathsf{S}}(t)S^{\dagger}_{\alpha}(\omega)-\dfrac{1}{2}\left\\{S_{\alpha}^{\dagger}(\omega)S_{\alpha^{\prime}}(\omega),\bm{\varrho}_{\mathsf{S}}(t)\right\\}\big{]},$
(8)

are the Lamb-shift Hamiltonian and the dissipator, respectively. Here the
$\gamma_{\alpha\alpha^{\prime}}(\omega)$ and
$\xi_{\alpha\alpha^{\prime}}(\omega)$ coefficients are related to the bath
correlation functions (Appendix A), and the Lindblad operators
$S_{\alpha}(\omega)$ are defined based on the system Hamiltonian as

$\displaystyle
S_{\alpha}(\omega)=\textstyle{\sum_{\epsilon-\epsilon^{\prime}=\hbar\omega}}|\epsilon\rangle\langle\epsilon|S_{\alpha}|\epsilon^{{}^{\prime}}\rangle\langle\epsilon^{{}^{\prime}}|.$
(9)

###  II.2 Floquet-Lindblad scenario

Now we study a driven open quantum system whose Hamiltonian is periodic in
time, $H_{\mathsf{S}}(t+\tau)=H_{\mathsf{S}}(t)$ with period $\tau$. It is
known from the Floquet theorem [38, 39] that the unitary operator
$U_{\mathsf{S}}(t,0)$ satisfying the equation
$i\hbar\frac{d}{dt}U_{\mathsf{S}}(t,0)=H_{\mathsf{S}}(t)U_{\mathsf{S}}(t,0)$
can also be represented in the form

$\displaystyle U_{\mathsf{S}}(t,0)=P(t,0)e^{-i\bar{H}t/\hbar},$ (10)

where $\bar{H}$ is a Hermitian time-independent operator, referred to as the
Floquet Hamiltonian, and $P(t,0)$ is a time-periodic unitary operator,
referred to shortly as the periodic operator, such that
$P(t+n\tau,n\tau)=P(t+n\tau,0)=P(t,0)$ and $P(n\tau,0)=\mathbbmss{I}$,
$\forall n\in\mathbbmss{N}$.

Starting from the dynamical equation (5) and employing the Floquet theorem, it
has been argued that under appropriate Born-Markov and secular approximations
the dynamical equation of an open system $\mathsf{S}$ is given by [21, 33]

$\displaystyle\dfrac{d\bm{\varrho}_{\mathsf{S}}(t)}{dt}=-\dfrac{i}{\hbar}\big{[}H_{\mathrm{lamb}}^{(\mathrm{F})},\bm{\varrho}_{\mathsf{S}}(t)\big{]}+\mathpzc{D}^{(\mathrm{F})}\left[\bm{\varrho}_{\mathsf{S}}(t)\right],$
(11)

where

$\displaystyle
H_{\mathrm{lamb}}^{(\mathrm{F})}=\textstyle{\sum_{\alpha,\alpha^{\prime}}}\sum_{{\omega}}\sum_{q\in\mathbbmss{Z}}\dfrac{1}{\hbar}\xi_{\alpha\alpha^{\prime}}(\omega+q\Omega)S_{\alpha}^{\dagger}(q,\omega)S_{\alpha^{\prime}}(q,\omega),$
(12)
$\displaystyle\mathpzc{D}^{(\mathrm{F})}[\bm{\varrho}_{\mathsf{S}}(t)]=\dfrac{1}{\hbar^{2}}\textstyle{\sum_{\alpha,\alpha^{\prime}}}\sum_{{\omega}}\sum_{q\in\mathbbmss{Z}}\gamma_{\alpha\alpha^{\prime}}(\omega+q\Omega)$
$\displaystyle\times\big{[}S_{\alpha^{\prime}}(q,\omega)\bm{\varrho}_{\mathsf{S}}(t)S^{\dagger}_{\alpha}(q,\omega)-\dfrac{1}{2}\left\\{S_{\alpha}^{\dagger}(q,\omega)S_{\alpha^{\prime}}(q,\omega),\bm{\varrho}_{\mathsf{S}}(t)\right\\}\big{]},$
(13)

are the Lamb-shift Hamiltonian and the dissipator, respectively. Here the
Floquet-Lindblad operators $S_{\alpha}(q,\omega)$ are defined as follows.
First, note that we can Fourier transform the periodic factor
$P^{{\dagger}}(t,0)S_{\alpha}P(t,0)$ as
$P^{{\dagger}}(t,0)S_{\alpha}P(t,0)=\textstyle{\sum_{q\in\mathbbmss{Z}}}S_{\alpha}(q)e^{iq\Omega
t}$, with $\Omega=2\pi/\tau$ being the frequency of the applied field. Next,
we define

$S_{\alpha}(q,\omega)=\textstyle{\sum_{\bar{\epsilon}-\bar{\epsilon}^{\prime}=\hbar\omega}}|\bar{\epsilon}\rangle\langle\bar{\epsilon}|S_{\alpha}(q)|\bar{\epsilon}^{\prime}\rangle\langle\bar{\epsilon}^{\prime}|,$
(14)

where $\bar{H}|\bar{\epsilon}\rangle=\bar{\epsilon}|\bar{\epsilon}\rangle$. In
addition, the coefficients $\gamma_{\alpha\alpha^{\prime}}(\omega+q\Omega)$
and $\xi_{\alpha\alpha^{\prime}}(\omega+q\Omega)$ are the real and imaginary
parts of the one-sided Fourier transformation of the bath correlation
functions. For details, see Appendix A. An immediate observation is that,
unlike the Lindblad scenario, in this case the Lamb-shift Hamiltonian (12)
does not necessarily commute with the system Hamiltonian.

###  II.3 Case study: $3$-level system

Figure 3: Time evolution of the reduced density matrix elements of a $3$-level
system (see Sec. II). The system is coupled to two baths with
$\beta_{\mathsf{c}}/\beta_{\mathsf{h}}=30/4$. (a) Applying no external driving
field on the system, (b) Applying an external driving field $V_{0}(t)$. Note
that here $\Omega/\omega_{b0}=0.9$. (c) Applying an external driving field
$V_{1}(t)$. Note also that $\Omega/\omega_{1b}=4.5$. The properties of the
system, baths, and the external fields are listed in Table 1 and the paragraph
after Eq. (23).

Consider a driven $3$-level system which is weakly coupled to two bosonic
baths $\mathsf{B}_{\mathsf{h}}$ and $\mathsf{B}_{\mathsf{c}}$ with
temperatures $T_{\mathsf{h}}$ and $T_{\mathsf{c}}$ ($T_{\mathsf{h}}\geqslant
T_{\mathsf{c}}$); see Fig. 2. The baths are assumed sufficiently large and
initially in their associated equilibrium states
$\varrho_{\mathsf{B}}=e^{-\beta_{\mathsf{B}}H_{\mathsf{B}}}/Z_{\mathsf{B}}$,
where $\beta_{\mathsf{B}}=1/(k_{B}T_{\mathsf{B}})$ and
$Z_{\mathsf{B}}=\mathrm{Tr}[e^{-\beta_{\mathsf{B}}H_{\mathsf{B}}}]$, with
$\mathsf{B}\in\\{\mathsf{B}_{\mathsf{c}},\mathsf{B}_{\mathsf{h}}\\}$. One can
simulate the baths with electromagnetic radiation fields which interact with a
system through dipole interaction. A careful analysis of this interaction is
reviewed in Appendix B. As a result, the total Hamiltonian of a driven
$3$-level system interacting with a radiation field can be written as

$\displaystyle
H_{\mathsf{tot}}(t)=H_{\mathsf{S}}(t)+H_{\mathsf{B}_{\mathsf{h}}}+H_{\mathsf{B}_{\mathsf{c}}}+H_{\mathsf{S}\mathsf{B}_{\mathsf{h}}}+H_{\mathsf{S}\mathsf{B}_{\mathsf{c}}},$
(15)

where

$\displaystyle H_{\mathsf{S}}(t)=H_{0}+V_{i}(t),$ (16) $\displaystyle
H_{0}=\epsilon_{0}|\epsilon_{0}\rangle\langle\epsilon_{0}|+\epsilon_{b}|\epsilon_{b}\rangle\langle\epsilon_{b}|+\epsilon_{1}|\epsilon_{1}\rangle\langle\epsilon_{1}|,$
(17) $\displaystyle V_{i}(t)=\mu\cos(\Omega
t)\big{(}|\epsilon_{b}\rangle\langle\epsilon_{i}|+|\epsilon_{i}\rangle\langle\epsilon_{b}|\big{)},\qquad
i\in\\{0,1\\}$ (18) $\displaystyle
H_{\mathsf{B}_{\mathsf{h}}}=\textstyle{\sum_{k}}\hbar\zeta_{k}\big{(}\hat{a}_{k}^{\dagger}\hat{a}_{k}+1/2\big{)},$
(19) $\displaystyle
H_{\mathsf{B}_{\mathsf{c}}}=\textstyle{\sum_{q}}\hbar\nu_{q}\big{(}\hat{b}_{q}^{\dagger}\hat{b}_{q}+1/2\big{)},$
(20) $\displaystyle
H_{\mathsf{S}\mathsf{B}_{\mathsf{h}}}=\textstyle{\sum_{k}}f_{k}\big{(}|\epsilon_{0}\rangle\langle\epsilon_{1}|\otimes\hat{a}^{\dagger}_{k}+|\epsilon_{1}\rangle\langle\epsilon_{0}|\otimes\hat{a}_{k}\big{)},$
(21) $\displaystyle
H_{\mathsf{S}\mathsf{B}_{\mathsf{c}}}=\textstyle{\sum_{q}}g_{q}\big{(}|\epsilon_{b}\rangle\langle\epsilon_{1}|\otimes\hat{b}^{\dagger}_{q}+|\epsilon_{1}\rangle\langle\epsilon_{b}|\otimes\hat{b}_{q}\big{)}.$
(22)

Here $H_{0}$ denotes the free Hamiltonian of the nondriven system, with energy
levels $\epsilon_{0}<\epsilon_{b}<\epsilon_{1}$. For brevity, hereafter and
throughout the paper we use the shorthand $|i\rangle$ rather than
$|\epsilon_{i}\rangle$. (For an $n$-level system, the free Hamiltonian is
extended to

$H_{0}=\textstyle{\sum_{i}}\epsilon_{i}|i\rangle\langle
i|,\,\,i\in\\{0,b,1,2,\ldots,n-2\\},$ (23)

with $\epsilon_{0}<\epsilon_{b}<\epsilon_{1}<\ldots<\epsilon_{n-2}$.) The
external field $V_{0}(t)$ ($V_{1}(t)$) couples the ground level $|0\rangle$
and the target level $|b\rangle$ (level $|1\rangle$ to the target level
$|b\rangle$). The operators $\hat{a}^{\dagger}_{k}$ and
$\hat{b}^{\dagger}_{k}$ are the creation operators of the $k$th mode of the
hot and the cold baths, respectively, which satisfy the bosonic commutation
relations
$[\hat{a}_{l},\hat{a}_{k}^{\dagger}]=[\hat{b}_{l},\hat{b}_{k}^{\dagger}]=\delta_{lk}$.
Moreover, the interaction Hamiltonian between the system and the hot and the
cold baths are given, respectively, by Eqs. (21) and (22), where we assume
$f_{k}$s and $g_{k}$s to be real-valued parameters.

Table 1: Temperatures and constants of the spectral density functions of the cold and hot baths [denoted by $(\mathsf{c})$ and $(\mathsf{h})$], in the form $J(\omega)=J_{0}\omega e^{-\omega^{2}/\omega^{2}_{\mathrm{cutoff}}}$, in natural units. $1/\beta_{\mathsf{h}}$ | $1/\beta_{\mathsf{c}}$ | $J_{0}^{(\mathsf{h})}$ | $J_{0}^{(\mathsf{c})}$ | $\omega_{\mathrm{cutoff}}^{(\mathsf{h})}$ | $\omega_{\mathrm{cutoff}}^{(\mathsf{c})}$
---|---|---|---|---|---
$30$ | $4$ | $4\times 10^{-4}$ | $4\times 10^{-3}$ | $\sqrt{2}$ | $\sqrt{0.2}$

In the Born-Markov and secular approximations, the evolution of system
$\mathsf{S}$ is governed by Eq. (11). Specifically, here we consider $H_{0}$
with $(\epsilon_{0},\epsilon_{b},\epsilon_{1})=(0,2.5,3)$ and
$\\{|0\rangle=(1,0,0)^{T},|1\rangle=(0,1,0)^{T},|b\rangle=(0,0,1)^{T}\\}$. The
system is assumed to be initially in the state
$\varrho_{\mathsf{S}}(0)=|0\rangle_{\mathsf{S}}\langle 0|$, and is coupled to
two thermal baths with $\beta_{\mathsf{c}}/\beta_{\mathsf{h}}=30/4$. We also
assume $(\mu,\Omega)=(0.1,2.25$). In all simulations throughout the paper, we
consider natural units where we set $\hbar\equiv c\equiv k_{\mathrm{B}}\equiv
1$, with $c$ being the speed of light and $k_{\mathrm{B}}$ the Boltzmann
constant. The properties of the baths are listed in Table 1. We consider three
cases: (i) there is no driving field; (ii) there is an external field
$V_{0}(t)$ which couples the $|0\rangle$ and $|b\rangle$ levels; and (iii) the
external field $V_{1}(t)$ couples the $|1\rangle$ and $|b\rangle$ levels.
(Table 2 in Appendix C shows the parameters of the Floquet-Lindblad master
equation for this example.) We also study the effect of the Lamb-shift in the
Lindblad and Floquet-Lindblad scenarios for these cases.

Figure 4: The comparison of energy transfer for the three cases of Fig. 3,
when $\beta_{\mathsf{c}}/\beta_{\mathsf{h}}=30/4$; see Sec. II. Figure 5:
Lamb-shift vs. no-Lamb-shift in the dynamics (see Sec. II). (a) As Fig. 3 a,
with no applied field. (b) As Fig. 3 b, with a near-resonant applied field
$V_{0}(t)$. (c) As Fig. 3 c, with an offresonant applied field $V_{1}(t)$. All
dashed lines indicate the cases with no-Lamb-shift Hamiltonian.

Figure 3 shows the evolution of the state of the system. When the system is
not driven, it is observed that it reaches its stationary state monotonically;
Fig. 3 a. Since there is no coherence in the initial state and also noting
that the master equation of the system is the Lindblad equation, no coherence
is generated during this evolution either. However, when the system is driven
by a time-dependent periodic field (Figs. 3 b and 3 c), it is observed that
the system eventually reaches a stationary state, experiencing some
fluctuations in the meanwhile. The amplitudes of these transient fluctuations
depend on the ratio $\Omega/\omega_{b0}$ (which for Fig. 3 b is in the near-
resonance regime) and $\Omega/\omega_{1b}$ (which for Fig. 3 c is in the
offresonance regime). In the former case, at earlier times the amplitude of
the fluctuations are considerable, whereas eventually the baths would suppress
them. In the latter case, in contrast, the system effectively experiences an
effective static potential, which implies negligible fluctuations [38]. We
note that at long times $|\varrho_{0b}|=0$ [in case (ii)] and
$|\varrho_{b1}|=0$ [in case (iii)], which implies a bath-induced decoherence.
In addition, by comparing Figs. 3 a – 3 c, it is observed that in case (iii)
the driven system reaches its stationary state relatively sooner than in the
other cases. See Appendix C for details of calculations.

Figure 4 shows the time evolution of the population of the target level
($\varrho_{bb}$) and the efficiency of energy transfer [as defined in Eq. (1)]
for the three cases. It is observed that the external field improves this
efficiency. However, at earlier times for case (ii) the efficiency is higher
compared to the two other cases. At long times, at $t\approx 6000$, when
$V_{0}(t)$ is applied, the efficiency grows by $\approx 7\%$ compared to the
nondriven case. Similarly, when $V_{1}(t)$ is applied, the efficiency
increases by $\approx 8\%$ compared to the nondriven case. As a remark, note
that for two baths with equal temperature ratios, their interactions with the
corresponding energy levels can be guided through appropriate frequency
filters [40, 41, 42, 43].

In practical studies of open quantum systems through master equations, it is
often assumed that the Lamb-shift Hamiltonian $H_{\mathrm{lamb}}$ is
negligible [44, 45, 46]. However, a remark is in order. We note that
$H_{\mathrm{lamb}}$ in the Lindblad and Floquet-Lindblad scenarios is of the
second-order with respect to $H_{\mathsf{int}}$—so is the dissipator [21, 33].
One can better analyze the effect of $H_{\mathrm{lamb}}$ on the system
evolution by removing this term from the master equation and comparing the
results with the complete equation. In the Lindblad scenario, since the master
equation is a rate-type equation for the population, then the Lamb term does
not affect the population evolutions for the nondegenerate case (see Fig. 5
a). However, the evolution of the offdiagonal terms of the state of the system
$\varrho_{\mathsf{S}}$ (in the eigenbasis of $H_{\mathsf{S}}$) depends on the
Lamb-shift Hamiltonian; that is, $H_{\mathrm{lamb}}$ contributes to the
coherence effects of $H_{\mathsf{S}}$ in the dynamics. Hence, due to its
smallness, neglecting $H_{\mathrm{lamb}}$ compared to $H_{\mathsf{S}}$ seems
justifiable.

Nevertheless, note that unlike the Lindblad case, in the Floquet-Lindblad
scenario keeping the Lamb-shift Hamiltonian affects the system evolution
because $H_{\mathrm{lamb}}^{(\mathrm{F})}$ does not commute with the bare
system Hamiltonian $H_{\mathsf{S}}$. Thus, care must be taken when in
particular one deals with such cases. From Figs. 5 b and 5 c we observe that
neglecting the Lamb-shift Hamiltonian in the Lindblad scenario as well as the
offresonance Floquet-Lindblad scenario are plausible assumptions, because in
these cases coherent oscillations are already negligible. Nevertheless, in the
near-resonance regime of the Floquet-Lindblad scenario this assumption does
not hold unless at the long-time limit of the dynamics.

##  III Redfield and Floquet-Redfield scenarios

Here we study the effect of the coherence on the energy transfer efficiency.
To do so, we consider two levels $|1\rangle$ and $|2\rangle$ for the system
Hamiltonian and investigate the effect of their coherence on the evolution of
the population of the target level $|b\rangle$. We note that the Lindblad and
the Floquet-Lindblad master equations [Eqs. (6) and (11)] are not appropriate
for this particular purpose because neither of them can couple the evolution
of the populations to the coherences of the system state. However, if we relax
the full secular approximation and retain only the Born-Markov approximation,
the resulting Redfield equations for the nondriven and driven cases can offer
a physically relevant framework to study the coherence effect. As we explain
later, in the driven case (the Floquet-Redfield scenario) we shall still need
a suitable partial secular approximation to guarantee physicality of the
system states.

###  III.1 Redfield scenario

Consider an open quantum system with a time-independent Hamiltonian
$H_{\mathsf{S}}$. When the system and the bath interact weakly so that we can
apply the Born-Markov approximation, then the dynamical master equation of the
system is given by the Redfield equation, whose general form is as follows
[21]:

$\displaystyle\dfrac{d\bm{\varrho}_{\mathsf{S}}(t)}{dt}=$
$\displaystyle\dfrac{1}{\hbar^{2}}\textstyle{\int^{\infty}_{0}}ds\,\mathrm{Tr}_{\mathsf{B}}\big{[}\bm{H}_{\mathsf{int}}(t-s)\bm{\varrho}_{\mathsf{S}}(t)\otimes\bm{\varrho}_{\mathsf{B}}\bm{H}_{\mathsf{int}}(t)$
$\displaystyle-\bm{H}_{\mathsf{int}}(t)\bm{H}_{\mathsf{int}}(t-s)\bm{\varrho}_{\mathsf{S}}(t)\otimes\bm{\varrho}_{\mathsf{B}}+\mathrm{h.c.}\big{]}.$
(24)

This is to some extent similar to the Lindblad equation but without invoking
the secular approximation. To remind the derivation, see Appendix A (where we
need to consider $H_{\mathsf{S}}$ time-independent and all discussions up to
Eq. (50) go through). This equation has been vastly used in various
applications to study problems such as light-harvesting [27].

Figure 6: A $4$-level system which is in contact with two thermal baths. Two
pair of the states are coupled with hot bath, and two other pair with the cold
bath. (a) The system is not driven. (b) An external field couples the
$|0\rangle$ and $|b\rangle$ levels to each other ($V_{0}$ field).

To be specific, let us consider a multilevel system interacting with a bath of
radiation field, where the total Hamiltonian is given by Eq. (2) with
$H_{\mathsf{S}}=H_{0}$ as in Eq. (23) and the system and the field are coupled
through the electric dipole interaction in the “$\vec{E}\cdot\vec{r}$”
representation [47, 48, 49] (see also Appendix B),

$\displaystyle
H_{\mathsf{B}}=\textstyle{\sum_{k,\lambda}}\hbar\nu_{k}\big{(}\hat{a}_{k,\lambda}^{\dagger}\hat{a}_{k,\lambda}+1/2\big{)},$
(25) $\displaystyle
H_{\mathsf{int}}=\textstyle{\sum_{i,j}}\textstyle{\sum_{k,\lambda}}g_{k\lambda}^{(ij)}|i\rangle\langle
j|\otimes\hat{a}_{k,\lambda}+\mathrm{h.c.}$ (26)

Here “h.c.” denotes Hermitian conjugate, $\hat{a}^{{\dagger}}_{k\lambda}$ is
the creation operator of a bath bosonic mode (harmonic oscillator) with
frequency $\nu_{k}$, wave vector $\\{k\\}$, and polarization $\lambda$, and
$g_{k\lambda}^{(ij)}$ is the coupling constant of the dipole interaction of
the $|i\rangle\leftrightarrow|j\rangle$ transition of the system with the
$(k,\lambda)$th mode of the bath; in particular,
$g_{k\lambda}^{(ij)}=\vec{\mu}_{ij}\cdot\hat{e}_{k\lambda}$, with
$\hat{e}_{k\lambda}$ being the unit polarization vector,
$\vec{\mu}_{ij}=\langle i|\vec{\mu}|j\rangle$ being the electric-dipole
transition matrix elements in the system energy basis [50].

After some algebra one obtains (see Appendix D) the master equation in the
Schrödinger picture as

$\displaystyle\dfrac{d\varrho_{\mathsf{S}}(t)}{dt}=-\dfrac{i}{\hbar}\left[H_{\mathsf{S}},\varrho_{\mathsf{S}}(t)\right]-\dfrac{1}{6\pi
c^{3}\hbar\epsilon_{0}}\textstyle{\sum_{i}}\sum_{i^{\prime}j^{\prime}}$
$\displaystyle\big{[}(\vec{\mu}_{i^{\prime}i}\cdot\vec{\mu}_{i^{\prime}j^{\prime}})N_{1}(\omega_{i^{\prime}j^{\prime}},\beta)\left(|i\rangle\langle
j^{\prime}|\varrho_{\mathsf{S}}(t)+\varrho_{\mathsf{S}}(t)|j^{\prime}\rangle\langle
i|\right)+$
$\displaystyle(\vec{\mu}_{ij^{\prime}}\cdot\vec{\mu}_{i^{\prime}j^{\prime}})N_{2}(\omega_{i^{\prime}j^{\prime}},\beta)\left(|i\rangle\langle
i^{\prime}|\varrho_{\mathsf{S}}(t)+\varrho_{\mathsf{S}}(t)|i^{\prime}\rangle\langle
i|\right)\big{]}$ $\displaystyle+\dfrac{1}{6\pi
c^{3}\hbar\epsilon_{0}}\textstyle{\sum_{ij}}\sum_{i^{\prime}j^{\prime}}\vec{\mu}_{ij}\cdot\vec{\mu}_{i^{\prime}j^{\prime}}\times$
$\displaystyle\big{[}N_{1}(\omega_{i^{\prime}j^{\prime}},\beta)\left[|i\rangle\langle
j|\varrho_{\mathsf{S}}(t)|j^{\prime}\rangle\langle
i^{\prime}|+|i^{\prime}\rangle\langle
j^{\prime}|\varrho_{\mathsf{S}}(t)|j\rangle\langle i|\right]+$ $\displaystyle
N_{2}(\omega_{i^{\prime}j^{\prime}},\beta)\left[|j\rangle\langle
i|\varrho_{\mathsf{S}}(t)|i^{\prime}\rangle\langle
j^{\prime}|+|j^{\prime}\rangle\langle
i^{\prime}|\varrho_{\mathsf{S}}(t)|i\rangle\langle j|\right]\big{]}$
$\displaystyle-\dfrac{i}{\hbar}\sum_{ii^{\prime}j}\dfrac{\vec{\mu}_{ij}\cdot\vec{\mu}_{i^{\prime}j}}{6\epsilon_{0}\pi^{2}c^{3}}C_{2}(\omega_{i^{\prime}j},\beta)\big{(}|i\rangle\langle
i^{\prime}|\varrho_{\mathsf{S}}(t)-\varrho_{\mathsf{S}}(t)|i^{\prime}\rangle\langle
i|\big{)}$
$\displaystyle-\dfrac{i}{\hbar}\sum_{ijj^{\prime}}\dfrac{\vec{\mu}_{ij}\cdot\vec{\mu}_{ij^{\prime}}}{6\epsilon_{0}\pi^{2}c^{3}}C_{1}(\omega_{ij^{\prime}},\beta)\big{(}\varrho_{\mathsf{S}}(t)|j^{\prime}\rangle\langle
j|-|j\rangle\langle j^{\prime}|\varrho_{\mathsf{S}}(t)\big{)}$
$\displaystyle+\dfrac{i}{\hbar}\sum_{ii^{\prime}jj^{\prime}}\dfrac{\vec{\mu}_{ij}\cdot\vec{\mu}_{i^{\prime}j^{\prime}}}{6\epsilon_{0}\pi^{2}c^{3}}\Big{[}C_{1}(\omega_{i^{\prime}j^{\prime}},\beta)\big{(}|i\rangle\langle
j|\varrho_{\mathsf{S}}(t)|j^{\prime}\rangle\langle i^{\prime}|-$
$\displaystyle|i^{\prime}\rangle\langle
j^{\prime}|\varrho_{\mathsf{S}}(t)|j\rangle\langle
i|\big{)}+C_{2}(\omega_{i^{\prime}j^{\prime}},\beta)\big{(}|j^{\prime}\rangle\langle
i^{\prime}|\varrho_{\mathsf{S}}(t)|i\rangle\langle j|-$
$\displaystyle|j\rangle\langle
i|\varrho_{\mathsf{S}}(t)|i^{\prime}\rangle\langle j^{\prime}|\big{)}\Big{]},$
(27)

where

$\displaystyle
C_{1}(x,\beta)=\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}\overline{n}(\nu,\beta)}{x-\nu},$
$\displaystyle
C_{2}(x,\beta)=\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}[\overline{n}(\nu,\beta)+1]}{x-\nu},$
$\displaystyle N_{1}(x,\beta)=x^{3}\bar{n}(x,\beta),$ $\displaystyle
N_{2}(x,\beta)=x^{3}[\bar{n}(x,\beta)+1].$ (28)

Here $\bar{n}(\omega,\beta)=(e^{\beta\hbar\omega}-1)^{-1}$. We see that in the
time-independent case, there is not any oscillatory terms on the right-hand
side (RHS) of the master equation in the Schrödinger picture. Hence it seems
that this equation [Eq. (27)] may generate nonnegative populations.

Before studying examples, there are two points to take into account. (i)
Unlike the Lindblad master equation, in general in the Redfield master
equation the spectral density function does not appear. This complicates the
comparison of the results of these two approaches. To alleviate this issue and
make a fairer comparison, we can study a specific example where the two
approaches can yield the similar master equations and then by comparing the
equations one can read corresponding parameters. For example, for a qubit
interacting with a thermal bath, with an odd spectral density function
[$J(-\omega)=-J(\omega)$], comparing the Lindblad and the Redfield equations
yield

$\displaystyle(\mu_{10}/\sqrt{6\pi c^{3}\hbar\epsilon_{0}})^{2}=\pi
J(\omega_{10})/\omega_{10}^{3},$ (29)

where $\omega_{10}$ is the qubit Hamiltonian gap. We remind that in all
examples in this paper, we choose our parameters according to the explanation
after Eq. (23), Table 1, and Eq. (29).

(ii) Note the integrals in Eq. (28) which appear in the Redfield and Floquet-
Redfield master equations. As explained in Appendix D, these integrals are
associated to the Lamb-shift terms. It has been known that the Lamb-shift
terms give rise to modifications of the energy levels. Since the corrections
of the energy levels due to the Lamb shift are typically small, in some
applications they can be safely ignored [44, 45, 46, 29]. In contrast to the
Lindblad scenarios, in the Redfield ones without introducing the spectral
density functions, some Lamb-shift terms mathematically diverge. However, by
introducing a suitable cutoff frequency (which is compatible with the dipole
approximation) and employing appropriate QED considerations [51, 52, 48, 53,
54], these terms take finite values and the divergence issue can be removed.
We also follow a similar approach to investigate the Lamb-shift effect in Sec.
III.3. In some situations such as nondegenerate or time-dependent system
Hamiltonians, we observe that the Lamb-shift terms lead to considerable
effects on the dynamics.

Now we study a specific example. To see the effect of neglecting the secular
approximation, we examine a $4$-level system which is coupled to two thermal
baths as in Fig. 6 a. The Hamiltonian of the system is

$H_{\mathsf{S}}=H_{0}=\textstyle{\sum_{i=0}^{2}}\epsilon_{i}|i\rangle\langle
i|+\epsilon_{b}|b\rangle\langle b|.$ (30)

We assume that the hot bath leads to the two transitions
$|0\rangle\leftrightarrow|1\rangle$ and $|0\rangle\leftrightarrow|2\rangle$,
whereas the cold bath causes the $|b\rangle\leftrightarrow|1\rangle$ and
$|b\rangle\leftrightarrow|2\rangle$ transitions. The system-bath interaction
Hamiltonian reads as

$\displaystyle
H_{\mathsf{int}}=H_{\mathsf{S}\mathsf{B}_{\mathrm{h}}}+H_{\mathsf{S}\mathsf{B}_{\mathrm{c}}},$
(31)

where

$\displaystyle
H_{\mathsf{S}\mathsf{B}_{\mathsf{h}}}=\textstyle{\sum_{i=1,2}}\textstyle{\sum_{k,\lambda}}g_{k\lambda}^{(i0)}|i\rangle\langle
0|\otimes\hat{a}_{k,\lambda}+\mathrm{h.c.},$ $\displaystyle
H_{\mathsf{S}\mathsf{B}_{\mathrm{c}}}=\textstyle{\sum_{i=1,2}}\textstyle{\sum_{q,\nu}}f_{q\nu}^{(ib)}|i\rangle\langle
b|\otimes\hat{b}_{q,\nu}+\mathrm{h.c.}$ (32)

Figure 7: Time evolution of a $4$-level nondriven system in contact with two
thermal baths with $\beta_{\mathsf{c}}/\beta_{\mathsf{h}}=30/4$. The Lamb-
shift term is not included here. (a) The levels $|1\rangle$ and $|2\rangle$
are degenerate. (b) The gap between $|1\rangle$ and $|2\rangle$ is assumed to
be $0.05$. All quantities are in natural units.

We study two cases according to the value of the gap between $|1\rangle$ and
$|2\rangle$. First, we consider the nonzero gap case
$\omega_{21}=(\epsilon_{2}-\epsilon_{1})/\hbar=0.05$, and next we focus on the
gapless case where $|1\rangle$ and $|2\rangle$ are degenerate states. At
first, let us disregard the integrals in Eq. (27) (i.e., $C_{1}(x,\beta)$ and
$C_{2}(x,\beta)$ coefficients) and analyze the solutions of this dynamical
equation. In Fig. 7 the evolution of the system in both degenerate and
nondegenerate cases are shown (for the corresponding analysis of the dynamics
with Lamb-shift terms, see Sec. III.3). In Figs. 7 a and 7 b the temperature
ratio of the two baths is $30/4$. When there is a gap between $|1\rangle$ and
$|2\rangle$, the final configuration of the system has populations in all
levels, which yields a relative reduction of the population of the target
level. As expected, according to Eq. (27) even in the nondegenerate case a
nonvanishing coherence between the levels $|1\rangle$ and $|2\rangle$ is
generated because the two $|0\rangle\leftrightarrow|1\rangle$ and
$|0\rangle\leftrightarrow|2\rangle$ ($|b\rangle\leftrightarrow|1\rangle$ and
$|b\rangle\leftrightarrow|2\rangle$) transitions are due to the same bath,
i.e., there are more than one way to populate the levels $|1\rangle$ and
$|2\rangle$. The role of coherence in the energy transfer efficiency is
studied in Sec. III.3, where the Lamb-shift terms are included in the master
equation.

###  III.2 Floquet-Redfield scenario

Now consider a quantum system which interacts with a bath and is driven by a
periodic external field. The total system Hamiltonian is

$\displaystyle
H_{\mathsf{tot}}(t)=H_{\mathsf{S}}(t)+H_{\mathsf{B}}+H_{\mathsf{int}},$ (33)

where the terms on the RHS are given by Eqs. (16), (23), (25), and (26).

Recalling Sec. II.2, one can obtain a related master equation for this
scenario. However, completely positive of this dynamics is not necessarily
guaranteed. To alleviate this issue, we apply a partial secular approximation
dictated by the assumption $\Omega^{-1}\ll\mathpzc{T}_{R}$, where
$\mathpzc{T}_{R}$ is the relaxation time of the driven system. This
approximation neglects all terms in the master equation where $q^{\prime}\neq
q$. This removes some fast oscillating terms. Note that this approximation is
different from that of Ref. [29]. After some algebra (see Appendix D for
details), we obtain the following Floquet-Redfield master equation in the
Schrödinger picture:

$\displaystyle\dfrac{d\varrho_{\mathsf{S}}(t)}{dt}=$
$\displaystyle-\dfrac{i}{\hbar}\left[H_{\mathsf{S}}(t),\varrho_{\mathsf{S}}(t)\right]-\dfrac{1}{6\pi
c^{3}\hbar\epsilon_{0}}\textstyle{\sum_{ij}}\sum_{i^{\prime}j^{\prime}}\sum_{q}\sum_{\omega\omega^{\prime}}(\vec{\mu}_{ij}\cdot\vec{\mu}_{i^{\prime}j^{\prime}})\times$
$\displaystyle\Big{[}\big{[}N_{2}(\omega^{\prime}+q\Omega\,,\beta)+C_{2}(\omega^{\prime}+q\Omega\,,\beta)\big{]}\times\big{[}L_{ij}(q,\omega;t)L^{\dagger}_{i^{\prime}j^{\prime}}(q,\omega^{\prime};t)\varrho_{\mathsf{S}}(t)-L^{\dagger}_{i^{\prime}j^{\prime}}(q,\omega^{\prime};t)\varrho_{\mathsf{S}}(t)L_{ij}(q,\omega;t)\big{]}$
$\displaystyle+\big{[}N_{2}(\omega^{\prime}+q\Omega\,,\beta)-C_{2}(\omega^{\prime}+q\Omega\,,\beta)\big{]}\times\big{[}\varrho_{\mathsf{S}}(t)L_{i^{\prime}j^{\prime}}(q,\omega^{\prime};t)L^{\dagger}_{ij}(q,\omega;t)-L^{\dagger}_{ij}(q,\omega;t)\varrho_{\mathsf{S}}(t)L_{i^{\prime}j^{\prime}}(q,\omega^{\prime};t)\big{]}$
$\displaystyle+\big{[}N_{1}(\omega^{\prime}+q\Omega\,,\beta)+C_{1}(\omega^{\prime}+q\Omega\,,\beta)\big{]}\times\big{[}\varrho_{\mathsf{S}}(t)L^{\dagger}_{i^{\prime}j^{\prime}}(q,\omega^{\prime};t)L_{ij}(q,\omega;t)-L_{ij}(q,\omega;t)\varrho_{\mathsf{S}}(t)L^{\dagger}_{i^{\prime}j^{\prime}}(q,\omega^{\prime};t)\big{]}$
$\displaystyle+\big{[}N_{1}(\omega^{\prime}+q\Omega\,,\beta)-C_{1}(\omega^{\prime}+q\Omega\,,\beta)\big{]}\times\big{[}L^{\dagger}_{ij}(q,\omega;t)L_{i^{\prime}j^{\prime}}(q,\omega^{\prime};t)\varrho_{\mathsf{S}}(t)-L_{i^{\prime}j^{\prime}}(q,\omega^{\prime};t)\varrho_{\mathsf{S}}(t)L^{\dagger}_{ij}(q,\omega;t)\big{]}\Big{]},$
(34)

where the coefficients are defined in Eqs. (28) and

$\displaystyle
L_{ij}(q,\omega;t)=P(t,0)\sigma_{ij}(q,\omega)P^{\dagger}(t,0).$ (35)

with $\sigma_{ij}(q,\omega)$ defined through

$\displaystyle\sigma_{ij}\equiv|i\rangle\langle j|,$ $\displaystyle
P^{\dagger}(t,0)\sigma_{ij}P(t,0)=\textstyle{\sum_{q\in\mathbbmss{Z}}}\sigma_{ij}(q)e^{iq\Omega
t},$
$\displaystyle\sigma_{ij}(q,\omega)=\textstyle{\sum_{\bar{\epsilon}-\bar{\epsilon}^{\prime}=\hbar\omega}}|\bar{\epsilon}\rangle\langle\bar{\epsilon}|\sigma_{ij}(q)|\bar{\epsilon}^{\prime}\rangle\langle\bar{\epsilon}^{\prime}|,$

Note that Eq. (34) includes the Lamb-shift terms.

Figure 8: Time evolution of the density matrix for the $4$-level driven system
in two different temperature ratios. The external field is given by Eq. (18).

As an example we consider a $4$-level system coupled to two thermal baths and
assume that an external field $V_{0}(t)$ couples the $|0\rangle$ and
$|b\rangle$ energy levels—see Fig. 6 b. The total system Hamiltonian is
described by Eq. (33), where

$\displaystyle H_{\mathsf{S}}(t)=H_{0}+V_{0}(t)$ (36)
$\displaystyle=\textstyle{\sum_{i=0}^{2}}\epsilon_{i}|i\rangle\langle
i|+\epsilon_{b}|b\rangle\langle b|+\mu\cos(\Omega t)\left(|0\rangle\langle
b|+|b\rangle\langle 0|\right).$

The Hamiltonians of the baths and their interaction are also given by Eqs.
(25) and (31), respectively.

At first, we ignore the contribution of the Lamb-shift terms in the dynamics,
i.e., we remove the $C_{1}(x,\beta)$ and $C_{2}(x,\beta)$ coefficients from
Eq. (34). In the next section, by considering Lamb-shift terms we will observe
that in the Redfield scenario of Sec. III.1, the energy transfer efficiency in
the degenerate case is larger than the nondegenerate case (Fig. 10). Hence in
the case of the driven system, we limit ourselves only to the degenerate case,
where the two levels $|1\rangle$ and $|2\rangle$ are considered to be
degenerate. For the temperature ratio
$\beta_{\mathsf{c}}/\beta_{\mathsf{h}}=30/4$ the time evolution of the state
of the system is represented in Fig. 8. It is observed that the coherence
$|\varrho_{12}|$, which is also equal to the populations of the degenerate
levels, is not destroyed by the bath. However, the coherence $|\varrho_{0b}|$,
which is due to the external field, is affected significantly by the baths. By
comparing Figs. 8 and 7, it is seen that when the dynamics is driven and given
by the Floquet-Redfield equation (34), the system approaches an almost
stationary configuration sooner than in the nondriven Redfield dynamics (27).

In the next section, we consider the full master equations (27) and (34) with
the Lamb-shift terms included and investigate the effect of the noise-induced
coherence as well as the Lamb-shift terms on the energy transfer efficiency.

###  III.3 Effect of the Lamb-shift terms in the Redfield master equations

Figure 9: Dynamics of a nondriven nondegenerate $4$-level system evolving
through Eq. (27). The solid curves correspond to the Redfield scenario with
the Lamb-shift correction, while the dotted curves correspond to the dynamics
without the Lamb-shift terms. Note that the evolution of the nondriven
degenerate system in the presence of the Lamb-shift terms is similar to Fig. 7
a.

In the derivation of the Redfield-type master equations we substitute
$\int_{0}^{\infty}ds\,e^{-i(x+\nu)s}=\pi\delta(x+\nu)-i\mathbbmss{P}1/(x+\nu)$—see
Appendix D.2. The delta function leads to terms which appear in the dissipator
and the Cauchy principal value results in the Lamb-shift term, which
contributes to the coherent evolution. Although in some studies the principal
value part is disregarded (which yield to ignoring the Lamb-shift Hamiltonians
in these master equations), a careful analysis of the Lamb-shift term is still
in order. To do so, we start from the Hamiltonian of a system interacting with
a radiation field. This Hamiltonian is written in the $\vec{E}\cdot\vec{r}$
representation and in the long-wavelength approximation. That is, only the
modes with frequency less than a cutoff frequency value $W<\infty$ are taken
into account. Hence the upper limit of all integrals in Eq. (28) should be
replaced with the parameter $W$. In addition, we assume a low-intensity and
low-energy radiation, which respectively yield to disregarding the second-
order terms $O(\vec{A}^{2})$ term and the interaction of the magnetic field of
the radiation and the spin of the system. Now consider the expression for
$C_{2}(\omega,\beta)$, which includes integrals of the following types:

$\displaystyle\mathbbmss{P}\int_{0}^{W}d\nu\,\dfrac{\nu^{3}}{\omega-\nu}+\mathbbmss{P}\int_{0}^{W}d\nu\,\dfrac{\nu^{3}\overline{n}(\nu,\beta)}{\omega-\nu},$
(37)

where

$\displaystyle\mathbbmss{P}\int_{0}^{W}d\nu\,\dfrac{\nu^{3}}{\omega-\nu}=-\int_{0}^{W}d\nu\,\nu^{2}-\omega\int_{0}^{W}d\nu\,\nu$
$\displaystyle-\omega^{2}\int_{0}^{W}d\nu+\mathbbmss{P}\int_{0}^{W}d\nu\dfrac{\omega^{3}}{\omega-\nu}.$
(38)

By employing appropriate QED arguments (see Appendix E), one can find that the
first term on the RHS of the above relation is completely compensated with
another term with opposite sign, which is attributed to a self-energy term. In
addition, the second term of the above relation can be concluded to be
irrelevant within the low-intensity approximation. Thus, keeping only the
relevant terms of Eq. (37), in the master equation we can effectively replace

$\displaystyle\mathbbmss{P}\int_{0}^{W}d\nu\,\dfrac{\nu^{3}\big{(}\overline{n}(\nu,\beta)+1\big{)}}{\omega-\nu}\to$
$\displaystyle\mathbbmss{P}\int_{0}^{W}d\nu\,\dfrac{\nu^{3}\overline{n}(\nu,\beta)}{\omega-\nu}-\omega^{2}W$
$\displaystyle+\omega^{3}\ln(W/\omega).$ (39)

Note that to establish the validity of the long-wavelength approximation, the
cutoff frequency $W$ is chosen large compared to the characteristic frequency
of the system and much smaller than the relativistic modes, $\omega_{0}\ll
W\ll mc^{2}/\hbar$, where $m$ is the mass of the charged particle [47].

Figure 10: Energy transfer efficiency vs. time for the nondriven cases of the
Lindblad and the Redfield scenarios. All quantities are in natural units.

At first we discuss the time-independent Redfield scenario in the presence as
well as absence of the Lamb-shift terms. Our simulations show that for the
time-independent degenerate $4$-level system (with $|1\rangle$ and $|2\rangle$
degenerate levels) the Lamb-shift terms do not affect the dynamics of the
system and the results are the same as Fig. 7 a. Thus, for the nondriven
scenario (Fig. 6 a), we only discuss the nondegenerate system. In all plots of
Fig. 9, we consider a nondriven nondegenerate $4$-level system with a $0.05$
gap between the $|1\rangle$ and $|2\rangle$ levels. This figure shows the
evolution of the state of the system in the presence as well as absence of the
Lamb-shift terms. It is evident that the Lamb-shift terms are nonnegligible.
Here we have taken $W\approx 4\times 10^{4}$ (in natural units). Note that the
correct values of the cutoff frequency $W$ in a model depends on its
specifics. Here, however, we have estimated this frequency simply by trial and
error such that the dynamics generates positive density matrices during the
evolution—for a similar analysis see Ref. [55].

To see the effect of coherence on the efficiency of energy transfer see Fig.
10, where we compare the Redfield scenario [Eq. (27)] with the Lindblad
scenario [Eq. (6)] for two degenerate and nondegenerate cases with the
temperature ratio $\beta_{\mathsf{c}}/\beta_{\mathsf{h}}=30/4$. At long times
the efficiencies of the Redfield scenario for both cases are greater than the
similar efficiencies of the Lindblad scenario. For example, at $t\approx
2800$, the Redfield dynamics yields $\approx 8\%$ and $\approx 24\%$ increase
in the efficiency over the Lindblad dynamics for the degenerate and
nondegenerate ($\omega_{12}=0.05$) cases, respectively. These increases of the
efficiency can be attributed to the coupling between the populations and
coherence terms in the Redfield scenarios.

Now let us consider the Floquet-Redfield scenario for a degenerate $4$-level
system whose dynamics is given by Eq. (34). Figure 11 compares the effects of
the full and the partial secular approximations, that is, it shows how the
Floquet-Lindblad and Floquet-Redfield master equations compare, note that the
Lamb-shift terms are maintained. From this figure, despite remarkable
difference between the two sets of plots at short times, it is seems that in
the long-time limit these equations yield similar solutions. In Fig. 12 the
energy transfer efficiency of these two cases are plotted. Since in the
Floquet-Redfield scenario the system reaches its stationary state relatively
sooner than the Floquet-Lindblad case, from short times the efficiency takes
its maximum value and remains so. Figure 13 compares the Floquet-Redfield
scenarios with and without the Lamb-shift terms. We observe that in the
Floquet-Redfield scenario for the degenerate system, the Lamb-shift terms have
an insignificant impact on the dynamics. It appears that for the driven
nondegenerate system considering the Lamb-shifts terms has nonnegligible
impact on the time evolution of the system. This is similar to the time-
independent case (Redfield scenario) where including the Lamb-shift terms
yielded nonnegligible effect on the time evolution for the nondegenerate
system (see Fig. 9),

Figure 11: State of the $4$-level system. Solid curves are for the Floquet-
Redfield dynamics (34), whereas the dotted curves correspond to the Floquet-
Lindblad dynamics (11). Both equations include the Lamb-shift terms. Figure
12: Effect of the coherence on the energy transfer efficiency in the presence
of the external field. The orange curve corresponds to the Floquet-Redfield
master equation (with the lamb-shift terms), and the blue one is related to
the Floquet-Lindblad master equation (with the Lamb-shift term).

##  IV Summary and conclusions

We have studied the problem of energy transfer in a multilevel system in
contact with two thermal baths at two different temperatures. The hot bath
excites the system to its highest energy level(s) and then the excitation
moves through the energy space in the presence of another cold bath. We have
considered four different dynamical scenarios within the Born-Markov regime:
(i) Lindblad scenario: assuming the secular approximation, with no applied
driving field on the system, (ii) Floquet-Lindblad scenario: assuming the
secular approximation and driving the system by a time-periodic driving field,
(iii) Redfield scenario: assuming no secular approximation and no driving
field on the system, and (iv) Floquet-Redfield scenario: assuming a particular
partial secular approximation and driving the system by a time-periodic
applied field. We have developed a particular dynamical equation for the last
scenario. Our overall objective has been to study the efficiency of energy
transfer to a particular energy level from other energy levels of the system
Hamiltonian in each of the dynamical scenarios.

Figure 13: Effect of the Lamb-shift terms in the Floquet-Redfield scenario.
Solid curves are related to the dynamics with the Lamb-shift terms, and the
dotted curves are associated with the no-Lamb-shift dynamics.

We have shown that in the two Lindblad scenarios, Lindblad and Floquet-
Lindblad, applying an external periodic field on a $3$-level system which
couples the target energy level to another energy level has yielded higher
efficiency for the energy transfer phenomenon. We have investigated the effect
of the Lamb-shift correction to the system Hamiltonian. In particular, as is
known we see that for the no-field Lindblad scenario the Lamb-shit term does
not have any effect on the dynamics of the system—and hence on the energy
transfer. A similar behavior also has manifested for the case in which an
applied field has coupled the target energy level to another level in the
offresonance regime. However, when the applied field is in the near-resonance
regime, it affects the dynamics at short times, while this effect diminishes
as time increases.

For the two Redfield scenarios, Redfield and Floquet-Redfield, we have
considered a $4$-level system, in which a coherence between the two upper
levels can be generated through the dynamics. In the Redfield scenario we have
considered both nondegenerate and degenerate cases. Since in the Redfield
scenario the energy transfer efficiency in the degenerate system grows more
than the nondegenerate case, we have only studied the degenerate system in the
Floquet-Redfield scenario. Moreover, we have investigate the effect of the
Lamb-shift terms in the Redfield scenarios and observed that at least in the
degenerate case the lamb-shift terms do not have any significant impact on the
dynamics.

Our careful and comprehensive study can be useful in the analysis of energy
transfer in synthetic systems, and it offers new outlooks about how one can
enhance such processes through suitable driving fields. As a next step, for
example, one can study how correlation effects may play a role in enhancing or
reducing energy transfer efficiency [56].

###### Acknowledgements.

The authors acknowledge helpful discussions with V. Rezvani. This work was
partially supported by Sharif University of Technology’s Office of Vice
President for Research and Technology through Contract No. QA960512.

##  Appendix A Derivation of the Floquet-Lindblad master equation

Here based on Ref. [33], we reproduce the derivation of the Floquet-Lindblad
master equation. Consider a quantum system which is in contact with a bath and
is driven by an external periodic field. We want to obtain a master equation
which describes the dynamic of the system. The total Hamiltonian is given by

$\displaystyle
H_{\mathsf{tot}}(t)=H_{\mathsf{S}}(t)+H_{\mathsf{B}}+H_{\mathsf{int}},$ (40)

where $H_{\mathsf{S}}(t)$ is the periodic driving Hamiltonian of the system
($H_{\mathsf{S}}(t+\tau)=H_{\mathsf{S}}(t)$ with period $\tau$),
$H_{\mathsf{B}}$ is the bath Hamiltonian, and $H_{\mathsf{int}}$ describes the
interaction between the system and the bath. Evolution of the total system in
the interaction picture is governed by

$\displaystyle\dfrac{d\bm{\varrho}_{\mathsf{tot}}(t)}{dt}=-\dfrac{i}{\hbar}\left[\bm{H}_{\mathsf{int}}(t),\bm{\varrho}_{\mathsf{tot}}(t)\right],$
(41)

where the interaction Hamiltonian in the interaction picture is defined as

$\displaystyle\bm{H}_{\mathsf{int}}(t)=U^{\dagger}_{\mathsf{S}}(t,0)\otimes
U^{\dagger}_{\mathsf{B}}(t,0)H_{\mathsf{int}}U_{\mathsf{S}}(t,0)\otimes
U_{\mathsf{B}}(t,0).$ (42)

From Eq. (41) we obtain

$\displaystyle\dfrac{d\bm{\varrho}_{\mathsf{tot}}(t)}{dt}=-\dfrac{i}{\hbar}\left\\{\left[\bm{H}_{\mathsf{int}}(t),\bm{\varrho}_{\mathsf{tot}}(0)\right]-\dfrac{i}{\hbar}\big{[}\bm{H}_{\mathsf{int}}(t),\textstyle{\int^{t}_{0}}ds\,\left[\bm{H}_{\mathsf{int}}(s),\bm{\varrho}_{\mathsf{tot}}(s)\right]\big{]}\right\\}.$
(43)

By considering the Markov approximation, the dynamical equation reduces to

$\displaystyle\dfrac{d\bm{\varrho}_{\mathsf{tot}}(t)}{dt}=-\dfrac{i}{\hbar}\left[\bm{H}_{\mathsf{int}}(t),\bm{\varrho}_{\mathsf{tot}}(0)\right]-\dfrac{1}{\hbar^{2}}\textstyle{\int^{t}_{0}}ds\,\big{[}\bm{H}_{\mathsf{int}}(t),\left[\bm{H}_{\mathsf{int}}(s),\bm{\varrho}_{\mathsf{tot}}(t)\right]\big{]},$
(44)

from whence

$\displaystyle\dfrac{d\bm{\varrho}_{\mathsf{S}}(t)}{dt}=-\dfrac{i}{\hbar}\mathrm{Tr}_{\mathsf{B}}\left[\bm{H}_{\mathsf{int}}(t),\bm{\varrho}_{\mathsf{tot}}(0)\right]-\dfrac{1}{\hbar^{2}}\textstyle{\int^{t}_{0}}ds\,\mathrm{Tr}_{\mathsf{B}}\big{[}\bm{H}_{\mathsf{int}}(t),\left[\bm{H}_{\mathsf{int}}(s),\bm{\varrho}_{\mathsf{tot}}(t)\right]\big{]}.$
(45)

Consider the general form of the interaction Hamiltonian in the Schrödinger
picture as

$\displaystyle H_{\mathsf{int}}=\textstyle{\sum_{\alpha}}S_{\alpha}\otimes
B_{\alpha},$ (46)

where $S_{\alpha}^{\dagger}=S_{\alpha}$ and $B_{\alpha}^{\dagger}=B_{\alpha}$.
Substituting Eq. (46) into Eq. (42) yields

$\displaystyle\bm{H}_{\mathsf{int}}(t)=\textstyle{\sum_{\alpha}}U^{\dagger}_{\mathsf{S}}(t,0)S_{\alpha}U_{\mathsf{S}}(t,0)\otimes
U^{\dagger}_{\mathsf{B}}(t,0)B_{\alpha}U_{\mathsf{B}}(t,0)=\sum_{\alpha}\bm{S}_{\alpha}(t)\otimes\bm{B}_{\alpha}(t),$
(47)

where
$U_{\mathsf{S}}(t,0)=\mathbbmss{T}e^{(-i/\hbar)\int_{0}^{t}ds\,H_{\mathsf{S}}(s)}$,
$U_{\mathsf{B}}(t,0)=e^{(-it/\hbar)H_{\mathsf{B}}}$,
$\bm{S}_{\alpha}(t)=U^{\dagger}_{\mathsf{S}}(t,0)S_{\alpha}U_{\mathsf{S}}(t,0)$,
and
$\bm{B}_{\alpha}(t)=U^{\dagger}_{\mathsf{B}}(t,0)B_{\alpha}U_{\mathsf{B}}(t,0)$.

If we assume that the initial state of the total system is an uncorrelated
state
$\varrho_{\mathsf{tot}}(0)=\varrho_{\mathsf{S}}(0)\otimes\varrho_{\mathsf{B}}(0)$
and the bath is prepared in a thermal state, then the first term in Eq. (45)
vanishes,

$\displaystyle\mathrm{Tr}_{\mathsf{B}}\left[\bm{H}_{\mathsf{int}}(t),\bm{\varrho}_{\mathsf{tot}}(0)\right]$
$\displaystyle=\mathrm{Tr}_{\mathsf{B}}\left[\textstyle{\sum_{\alpha}}\bm{S}_{\alpha}(t)\otimes\bm{B}_{\alpha}(t),\varrho_{\mathsf{S}}(0)\otimes\varrho_{\mathsf{B}}(0)\right]$
$\displaystyle=\textstyle{\sum_{\alpha}}\left[\bm{S}_{\alpha}(t),\varrho_{\mathsf{S}}(0)\right]\mathrm{Tr}_{\mathsf{B}}\left[\bm{B}_{\alpha}(t)\varrho_{\mathsf{B}}(0)\right]=0,$
(48)

because
$\mathrm{Tr}_{\mathsf{B}}\left[\bm{B}_{\alpha}(t)\varrho_{\mathsf{B}}(0)\right]=0$.
In addition, we assume the Born approximation, where the coupling between the
system and the bath is considered to be weak such that one can neglect the
effects of the system on the bath. This approximation is often incorporated by
writing the instantaneous state of the total system approximately as a tensor
product of the system state $\varrho_{\mathsf{S}}(t)$ and the bath state is
assumed almost unchanged in the course of time,
$\bm{\varrho}_{\mathsf{B}}(t)\approx\bm{\varrho}_{\mathsf{B}}(0)\equiv\bm{\varrho}_{\mathsf{B}}$;
that is,
$\bm{\varrho}_{\mathsf{tot}}(t)\approx\bm{\varrho}_{\mathsf{S}}(t)\otimes\bm{\varrho}_{\mathsf{B}}$.
Thus, by applying this approximation, Eq. (45) becomes

$\displaystyle\dfrac{d\bm{\varrho}_{\mathsf{S}}(t)}{dt}=-\dfrac{1}{\hbar^{2}}\textstyle{\int^{t}_{0}}ds\,\mathrm{Tr}_{\mathsf{B}}\big{[}\bm{H}_{\mathsf{int}}(t),\left[\bm{H}_{\mathsf{int}}(s),\bm{\varrho}_{\mathsf{S}}(t)\otimes\bm{\varrho}_{\mathsf{B}}\right]\big{]}.$
(49)

If we substitute $s\rightarrow t-s$ in Eq. (49) and put the upper limit of the
integral to be infinity (presuming that the integrand, which is proportional
to the bath correlation functions, does not have a substantial contribution
for $s\gg\mathpzc{T}_{\mathsf{B}}$) [21], this reduces Eq. (49) to the
following Markovian master equation:

$\displaystyle\dfrac{d\bm{\varrho}_{\mathsf{S}}(t)}{dt}=\dfrac{1}{\hbar^{2}}\textstyle{\int^{\infty}_{0}}ds\,\mathrm{Tr}_{\mathsf{B}}\big{[}\bm{H}_{\mathsf{int}}(t-s)\bm{\varrho}_{\mathsf{S}}(t)\otimes\bm{\varrho}_{\mathsf{B}}\bm{H}_{\mathsf{int}}(t)-\bm{H}_{\mathsf{int}}(t)\bm{H}_{\mathsf{int}}(t-s)\bm{\varrho}_{\mathsf{S}}(t)\otimes\bm{\varrho}_{\mathsf{B}}+\mathrm{h.c.}\big{]}.$
(50)

We now need to incorporate $\bm{H}_{\mathsf{int}}(t)$ in the above equation
from Eq. (47). To do so, we need to calculate $U_{\mathsf{S}}(t,0)$, which in
general requires a time-ordering as
$U_{\mathsf{S}}(t,0)=\mathbbmss{T}e^{-(i/\hbar)\int_{0}^{t}H_{\mathsf{S}}(s)\,ds}$.
However, noting that the system Hamiltonian $H_{\mathsf{S}}(t)$ is assumed to
be periodic, one can employ the Floquet theorem to partially simplify the form
of $U_{\mathsf{S}}(t,0)$.

The Floquet theorem states that [38] the time evolution operator associated to
a periodic Hamiltonian $H_{\mathsf{S}}(t)$,
$i\hbar\frac{d}{dt}U_{\mathsf{S}}(t,t_{0})=H_{\mathsf{S}}(t)U_{\mathsf{S}}(t,t_{0})$,
can be factorized as

$\displaystyle
U_{\mathsf{S}}(t,t_{0})=P(t,t_{0})e^{-i\bar{H}(t-t_{0})/\hbar},$ (51)

where $\bar{H}$ is a Hermitian time-independent operator, referred to as the
Floquet Hamiltonian, and $P(t,t_{0})$ is a unitary and periodic operator such
that $P(t+n\tau,n\tau)=P(t+n\tau,0)=P(t,0)$ and $P(n\tau,0)=\mathbbmss{I}$,
with $\tau$ being the period of the Hamiltonian and $n\in\mathbbmss{N}$.

The Floquet theorem helps us rewrite $\bm{S}_{\alpha}(t)$ operators in Eq.
(47) as

$\displaystyle\bm{S}_{\alpha}(t)=e^{i\bar{H}t/\hbar}P^{\dagger}(t,0)S_{\alpha}P(t,0)e^{-i\bar{H}t/\hbar}.$
(52)

Since the operator $P(t,0)$ is periodic, then
$P^{\dagger}(t,0)S_{\alpha}P(t,0)$ is also periodic and hence has a Fourier
expansion as

$\displaystyle
P^{{\dagger}}(t,0)S_{\alpha}P(t,0)=\textstyle{\sum_{q\in\mathbbmss{Z}}}S_{\alpha}(q)e^{iq\Omega
t},$ (53)

where $\Omega=2\pi/\tau$ is the frequency of the periodic Hamiltonian
$H_{\mathsf{S}}(t)$. The Fourier coefficients read

$\displaystyle
S_{\alpha}(q)=\dfrac{1}{\tau}\textstyle{\int^{\tau}_{0}}dt\,e^{-iq\Omega
t}P^{{\dagger}}(t,0)S_{\alpha}P(t,0),$ (54)

whose calculation needs the knowledge of the periodic operator $P(t,0)$ only
in the first period $0\leqslant t\leqslant\tau$.

One can construct another set of related operators as

$\displaystyle
S_{\alpha}(q,\omega)=\textstyle{\sum_{\bar{\epsilon}-\bar{\epsilon}^{\prime}=\hbar\omega}}|\bar{\epsilon}\rangle\langle\bar{\epsilon}|S_{\alpha}(q)|\bar{\epsilon}^{\prime}\rangle\langle\bar{\epsilon}^{\prime}|,$
(55)

where $\bar{\epsilon}$ and $|\bar{\epsilon}\rangle$ are, respectively, the
eigenvalues and eigenvectors of the Floquet Hamiltonian,
$\bar{H}|\bar{\epsilon}\rangle=\bar{\epsilon}|\bar{\epsilon}\rangle$. The
summation in Eq. (55) is over all eigenvalues with the same (quasi) energy
difference $\hbar\omega$, and whenever we encounter a $\sum_{\omega}$ it
denotes the summation over all possible (quasi) energy gaps.

It is useful to identify several symmetries here. By comparing Eqs. (53) and
(55) with their Hermitian conjugates, we find
$S_{\alpha}^{\dagger}(q)=S_{\alpha}(-q)$ and
$S_{\alpha}^{\dagger}(q,\omega)=S_{\alpha}(-q,-\omega)$. By summing
$S_{\alpha}(q,\omega)$ over all $\omega$s we have

$\displaystyle\textstyle{\sum_{\omega}}S_{\alpha}(q,\omega)=S_{\alpha}(q).$
(56)

Substituting Eq. (53) into Eq. (52) and then using Eq. (55), gives

$\displaystyle\bm{S}_{\alpha}(t)=\textstyle{\sum_{q\in\mathbbmss{Z}}}\sum_{\omega}e^{it(\omega+q\Omega)}S_{\alpha}(q,\omega),$
(57)

which in turn yields

$\displaystyle\bm{H}_{\mathsf{int}}(t)=\textstyle{\sum_{\alpha}}\textstyle{\sum_{q\in\mathbbmss{Z}}}\textstyle{\sum_{\omega}}e^{it(\omega+q\Omega)}S_{\alpha}(q,\omega)\otimes\bm{B}_{\alpha}(t).$
(58)

Inserting this form in Eq. (50) gives the following dynamical equation:

$\displaystyle\dfrac{d\bm{\varrho}_{\mathsf{S}}(t)}{dt}=$
$\displaystyle\dfrac{1}{\hbar^{2}}\textstyle{\sum_{\alpha,\alpha^{\prime}}}\sum_{{\omega},{\omega^{\prime}}}\sum_{q,q^{\prime}}e^{it\left[(\omega-\omega^{\prime})+(q-q^{\prime})\Omega\right]}\Big{[}\left(S_{\alpha^{\prime}}(q,\omega)\bm{\varrho}_{\mathsf{S}}(t)S_{\alpha}^{\dagger}(q^{\prime},\omega^{\prime})-S_{\alpha}^{\dagger}(q^{\prime},\omega^{\prime})S_{\alpha^{\prime}}(q,\omega)\bm{\varrho}_{\mathsf{S}}(t)\right)$
$\displaystyle\times\textstyle{\int^{\infty}_{0}}ds\,e^{-is(\omega+q\Omega)}\mathrm{Tr}_{\mathsf{B}}\big{[}\bm{B}_{\alpha}^{\dagger}(t)\bm{B}_{\alpha^{\prime}}(t-s)\varrho_{\mathsf{B}}\big{]}\Big{]}+\mathrm{h.c.}$
(59)

This from can still be simplified further by introducing the one-sided Fourier
transformation of the bath correlation functions as

$\displaystyle\Gamma_{\alpha\alpha^{\prime}}(\omega+q\Omega)\equiv\textstyle{\int^{\infty}_{0}}ds\,e^{-is(\omega+q\Omega)}\langle\bm{B}_{\alpha}^{\dagger}(t)\bm{B}_{\alpha^{\prime}}(t-s)\rangle,$
(60)

where
$\langle\bm{B}_{\alpha}^{\dagger}(t)\bm{B}_{\alpha^{\prime}}(t-s)\rangle=\mathrm{Tr}_{\mathrm{B}}\big{[}\bm{B}_{\alpha}^{\dagger}(t)\bm{B}_{\alpha^{\prime}}(t-s)\varrho_{\mathsf{B}}\big{]}$.
If we assume the thermal equilibrium state as $\varrho_{\mathsf{B}}=e^{-\beta
H_{\mathsf{B}}}/Z_{\mathsf{B}}$, with $\beta=1/k_{B}T$ and
$Z_{\mathsf{B}}=\mathrm{Tr}[e^{-\beta H_{\mathsf{B}}}]$, then one can see that
the bath correlation functions become homogeneous in time,

$\displaystyle\left\langle\bm{B}_{\alpha}^{\dagger}(t)\bm{B}_{\alpha^{\prime}}(t-s)\right\rangle=\left\langle\bm{B}_{\alpha}^{\dagger}(s)\bm{B}_{\alpha^{\prime}}(0)\right\rangle.$
(61)

Thus, we obtain

$\displaystyle\dfrac{d\bm{\varrho}_{\mathsf{S}}(t)}{dt}=$
$\displaystyle\dfrac{1}{\hbar^{2}}\textstyle{\sum_{\alpha,\alpha^{\prime}}}\sum_{{\omega},{\omega^{\prime}}}\sum_{q,q^{\prime}}\Gamma_{\alpha\alpha^{\prime}}(\omega+q\Omega)e^{it\left[(\omega-\omega^{\prime})+(q-q^{\prime})\Omega\right]}$
$\displaystyle\times\left[S_{\alpha^{\prime}}(q,\omega)\bm{\varrho}_{\mathsf{S}}(t)S_{\alpha}^{\dagger}(q^{\prime},\omega^{\prime})-S_{\alpha}^{\dagger}(q^{\prime},\omega^{\prime})S_{\alpha^{\prime}}(q,\omega)\bm{\varrho}_{\mathsf{S}}(t)\right]+\mathrm{h.c.}$
(62)

If the time scale
$\mathpzc{T}_{\mathsf{S}}=\max_{\begin{subarray}{c}\omega\neq\omega^{\prime},m\end{subarray}}|\omega-\omega^{\prime}+m\Omega|^{-1}$
is considerably small compared to the relaxation time scale of the driven
system ($\mathpzc{T}_{R}$), i.e.,
$\max_{\begin{subarray}{c}\omega\neq\omega^{\prime},m\end{subarray}}|\omega-\omega^{\prime}+m\Omega|^{-1}\ll\mathpzc{T}_{R}$,
then in the summation $\sum_{\alpha,\alpha^{\prime}}$ and
$\sum_{\omega,\omega^{\prime}}$ only the terms for which
$\omega=\omega^{\prime}$ and $q=q^{\prime}$ will have a nonvanishing
contribution. This condition is called the secular approximation, and gives

$\displaystyle\dfrac{d\bm{\varrho}_{\mathsf{S}}(t)}{dt}=$
$\displaystyle\dfrac{1}{\hbar^{2}}\textstyle{\sum_{\alpha,\alpha^{\prime}}}\sum_{{\omega}}\sum_{q}\Gamma_{\alpha\alpha^{\prime}}(\omega+q\Omega)\big{[}S_{\alpha^{\prime}}(q,\omega)\bm{\varrho}_{\mathsf{S}}(t)S_{\alpha}^{\dagger}(q,\omega)-S_{\alpha}^{\dagger}(q,\omega)S_{\alpha^{\prime}}(q,\omega)\bm{\varrho}_{\mathsf{S}}(t)\big{]}+\mathrm{h.c.}$
(63)

We can decompose $\Gamma_{\alpha\alpha^{\prime}}(x)$ to a real part and an
imaginary part as follows:

$\displaystyle\Gamma_{\alpha\alpha^{\prime}}(x)=\dfrac{1}{2}\gamma_{\alpha\alpha^{\prime}}(x)+i\xi_{\alpha\alpha^{\prime}}(x).$
(64)

By substituting the above expression in Eq. (63), we obtain the master
equation of the system in the interaction picture as follows:

$\displaystyle\dfrac{d\bm{\varrho}_{\mathsf{S}}(t)}{dt}=-\dfrac{i}{\hbar}\big{[}H_{\mathrm{lamb}}^{(\mathrm{F})},\bm{\varrho}_{\mathsf{S}}(t)\big{]}+\mathpzc{D}^{(\mathrm{F})}\left[\bm{\varrho}_{\mathsf{S}}(t)\right],$
(65)

where

$\displaystyle
H_{\mathrm{lamb}}^{(\mathrm{F})}=\dfrac{1}{\hbar}\textstyle{\sum_{\alpha,\alpha^{\prime}}}\sum_{{\omega}}\sum_{q}\xi_{\alpha\alpha^{\prime}}(\omega+q\Omega)S_{\alpha}^{\dagger}(q,\omega)S_{\alpha^{\prime}}(q,\omega),$
(66)
$\displaystyle\mathpzc{D}^{(\mathrm{F})}[\bm{\varrho}_{\mathsf{S}}(t)]=\dfrac{1}{\hbar^{2}}\textstyle{\sum_{\alpha,\alpha^{\prime}}}\sum_{{\omega}}\sum_{q}\gamma_{\alpha\alpha^{\prime}}(\omega+q\Omega)\big{[}S_{\alpha^{\prime}}(q,\omega)\bm{\varrho}_{\mathsf{S}}(t)S^{\dagger}_{\alpha}(q,\omega)-\dfrac{1}{2}\big{\\{}S_{\alpha}^{\dagger}(q,\omega)S_{\alpha^{\prime}}(q,\omega),\bm{\varrho}_{\mathsf{S}}(t)\big{\\}}\big{]}.$
(67)

Here $H_{\mathrm{lamb}}^{(\mathrm{F})}$ is referred to as the Lamb-shift
Hamiltonian, which is due to the interaction between the system and the bath.
The second term on the RHS of Eq. (11)
$\left(\mathpzc{D}^{(\mathrm{F})}[\bm{\varrho}_{\mathsf{S}}(t)]\right)$ is a
decoherence term induced by the bath.

For the sake of comparison, we close this appendix with recalling the standard
Lindblad master equation for a time-independent system in the interaction
picture,

$\displaystyle\dfrac{d\bm{\varrho}_{\mathsf{S}}(t)}{dt}=-\dfrac{i}{\hbar}\left[H_{\mathrm{lamb}},\bm{\varrho}_{\mathsf{S}}(t)\right]+\mathpzc{D}\left[\bm{\varrho}_{\mathsf{S}}(t)\right],$
(68)

where

$\displaystyle
H_{\mathrm{lamb}}=\textstyle{\sum_{\alpha,\alpha^{\prime}}}\sum_{{\omega}}\dfrac{1}{\hbar}\xi_{\alpha\alpha^{\prime}}(\omega)S_{\alpha}^{\dagger}(\omega)S_{\alpha^{\prime}}(\omega),$
(69)
$\displaystyle\mathpzc{D}[\bm{\varrho}_{\mathsf{S}}(t)]=\dfrac{1}{\hbar^{2}}\textstyle{\sum_{\alpha,\alpha^{\prime}}}\sum_{{\omega}}\gamma_{\alpha\alpha^{\prime}}(\omega)\big{[}S_{\alpha^{\prime}}(\omega)\bm{\varrho}_{\mathsf{S}}(t)S^{\dagger}_{\alpha}(\omega)-\dfrac{1}{2}\big{\\{}S_{\alpha}^{\dagger}(\omega)S_{\alpha^{\prime}}(\omega),\bm{\varrho}_{\mathsf{S}}(t)\big{\\}}\big{]}.$
(70)

Unlike the time-dependent case, here the Lindblad operators
$S_{\alpha}(\omega)$ are obtained according to the system Hamiltonian, i.e.,
$H_{\mathsf{S}}=H_{0}$. In addition, here the $\omega$’s are the gaps of the
system Hamiltonian. Note that the $\gamma_{\alpha\alpha^{\prime}}(\omega)$ and
$\xi_{\alpha\alpha^{\prime}}(\omega)$ coefficients are obtained by Eq. (64),
too.

##  Appendix B Hamiltonian of the charge-field interaction

Consider a particle of charge $q$ and mass $m$ (an atom or a molecule)
localized in the origin and coupled to a (quantized) electromagnetic radiation
field. Let us also assume the long-wavelength approximation, which neglects
the spatial variation of the electromagnetic field on the length scale of the
charged particle, and the Coulomb gauge for the electromagnetic field. Based
on elementary physical principles the total charge-field Hamiltonian is given
by [47]

$\displaystyle
H_{\mathrm{tot}}=\dfrac{1}{2m}[\vec{p}-q\vec{A}(0)]^{2}+V_{\mathrm{Coul}}(\vec{r})+\textstyle{\sum_{j}}\hbar\omega_{j}(\hat{a}^{\dagger}_{j}\hat{a}_{j}+1/2),$
(71)

where $\vec{p}$ is the momentum operator of the charged particle, $\vec{A}$ is
the vector potential operator associated with the electromagnetic field,
$V_{\mathrm{Coul}}$ is the Coulomb potential, $\\{\omega_{j}\\}$ are the
quantized frequencies or modes of the field, and $\hat{a}_{j}$ is the
annihilation operator corresponding to mode $\omega_{j}$ of the field. This
form of the Hamiltonian is usually referred to as the “$\vec{A}\cdot\vec{p}$”
representation. One can rewrite this Hamiltonian such that wherein three
separate terms can be discerned,

$\displaystyle H_{\mathrm{tot}}=H_{\mathsf{P}}+H_{\mathsf{R}}+H_{I},$ (72)

where

$\displaystyle H_{\mathsf{P}}=\dfrac{\vec{p}^{2}}{2m}+V_{\mathrm{Coul}},$ (73)
$\displaystyle
H_{\mathsf{R}}=\textstyle{\sum_{j}}\hbar\omega_{j}(\hat{a}^{\dagger}_{j}\hat{a}_{j}+1/2),$
(74) $\displaystyle H_{I}=H_{I1}+H_{I2},$ (75)

are the Hamiltonians associated to the particle, the radiation field, and the
interaction, respectively. Here $H_{I1}$ is linear in terms of the field
$\vec{A}$ and $H_{I2}$ is quadratic,

$\displaystyle H_{I1}=-\dfrac{q}{m}\vec{p}\cdot\vec{A}(0),$ (76)
$\displaystyle H_{I2}=\dfrac{q^{2}}{2m}[\vec{A}(0)]^{2}.$ (77)

For practical purposes, it might be more useful to transform the
$\vec{A}\cdot\vec{p}$ representation to another equivalent
“$\vec{E}\cdot\vec{r}$” representation which explicitly depends on the
associated electric field $\vec{E}$ rather than the vector potential
$\vec{A}$. The transformation between these two representations is carried out
by the unitary operator $T=e^{-(i/\hbar)q\vec{r}\cdot\vec{A}(0)}$, which
yields

$\displaystyle H^{\prime}_{\mathrm{tot}}$
$\displaystyle=TH_{\mathrm{tot}}T^{\dagger}$
$\displaystyle=\dfrac{\vec{p}^{2}}{2m}+V_{\mathrm{Coul}}+\textstyle{\sum_{j}}\hbar\omega_{j}(\hat{a}^{\dagger}_{j}\hat{a}_{j}+1/2)-\vec{d}\cdot\vec{E}(0)+\varepsilon_{\mathrm{dipole}},$
(78)

where $\varepsilon_{\mathrm{dipole}}$ is the dipole self-energy related to the
charged particle Hamiltonian (which is usually ignored in the literature), and
the fourth term is the expected dipole interaction, where $\vec{d}=q\vec{r}$
is the dipole operator of the charged particle. It is seen that another appeal
of the $\vec{E}\cdot\vec{r}$ representation is that in this representation the
Hamiltonian only depends on the electric field linearly. One can note that the
charge-field Hamiltonian (78) can be recast in the well-known form

$H^{\prime}_{\mathrm{tot}}=H_{\mathsf{P}}+H_{\mathsf{R}}-\vec{d}\cdot\vec{E}(0),$
(79)

where the first, second, and last terms denote the particle, the reservoir,
and the particle-reservoir dipole interaction, respectively. In this paper we
work in the $\vec{E}\cdot\vec{r}$ representation and simply denote the
associated Hamiltonian by $H_{\mathrm{tot}}$ (without prime).

######  Remark 1.

In some physical applications, one may assume the low-intensity approximation
which corresponds to ignoring the $O(\vec{A}^{2})$ term from the Hamiltonian
(71). In the $\vec{E}\cdot\vec{r}$ representation this approximation leads to
the subtraction of $H_{I2}$ from the Hamiltonian (78).

##  Appendix C Floquet Hamiltonian and periodic operator for a driven system

###  C.1 General analysis

As it was shown in the Born-Markov and secular approximations, the evolution
of a driven system with a periodic external field is governed by Eq. (11), and
for using this equation we should calculate the Floquet Hamiltonian. Thus, we
present a way to obtain the Floquet Hamiltonian approximately.

On the one hand, we have

$U_{\mathsf{S}}(t,0)=\mathbbmss{T}e^{-(i/\hbar)\int_{0}^{t}ds\,H_{\mathsf{S}}(s)},$
(80)

where $H_{\mathsf{S}}(t)=H_{0}+V(t)$ and $\mathbbmss{T}$ denotes time-
ordering. On the other hand, from the Floquet theory, we have

$\displaystyle U_{\mathsf{S}}(\tau,0)=e^{-i\bar{H}\tau/\hbar}.$ (81)

According to Eq. (81) to obtain the Floquet Hamiltonian $\bar{H}$ it is better
that we have the evolution operator in the exponential form, to do this notice
the following equation which relates the evolution operator in the both
Schrödinger and interaction picture together

$\displaystyle U_{\mathsf{S}}(t,0)=U_{0}(t,0)\bm{U}_{I}(t,0),$ (82)

where $U_{0}(t,0)=e^{-iH_{0}t/\hbar}$ and $\bm{U}_{I}(t,0)$ is the evolution
operator in the interaction picture which evolves as

$\displaystyle\dfrac{d\bm{U}_{I}(t,0)}{dt}=-\dfrac{i}{\hbar}\bm{V}_{I}(t)\bm{U}_{I}(t,0),$
(83)

where $\bm{V}_{I}(t)=e^{iH_{0}t/\hbar}V(t)e^{-iH_{0}t/\hbar}$. By
straightforward calculation the operator $\bm{V}_{I}(t)$ takes the following
form for the operator defined as Eq. (18) [see Fig. 2]

$\displaystyle\bm{V}_{I}(t)=\mu\cos(\Omega
t)\left[e^{-it\omega_{b0}}|0\rangle\langle
b|+e^{it\omega_{b0}}|b\rangle\langle 0|\right],$ (84)

where $\omega_{b0}\equiv(\epsilon_{b}-\epsilon_{0})/\hbar$ is the gap of the
two states which are couple together through the external field with strength
$\mu$. Thus the RHS of the differential Eq. (83) is proportional to $\mu$ and
we can apply the Magnus expansion to obtain the $\bm{U}_{I}(t,0)$ as follows:

$\displaystyle\bm{U}_{I}(t,0)=e^{\bm{\Lambda}(t)},$ (85)

where the operator $\bm{\Lambda}(t)$ obtains by the following series
expansion:

$\displaystyle\bm{\Lambda}(t)=\textstyle{\sum_{n=1}^{\infty}}(-i\mu/\hbar)^{n}\bm{\Lambda}_{n}(t).$
(86)

The three first terms of the series Eq. (86) is as follows [57]:

$\displaystyle\bm{\Lambda}_{1}(t)=\textstyle{\int_{0}^{t}}dt_{1}\,\bm{V}(t_{1}),$
$\displaystyle\bm{\Lambda}_{2}(t)=(1/2!)\textstyle{\int_{0}^{t}}dt_{1}\,\int_{0}^{t_{1}}dt_{2}\left[\bm{V}(t_{1}),\bm{V}(t_{2})\right],$
$\displaystyle\bm{\Lambda}_{3}(t)=(1/3!)\textstyle{\int_{0}^{t}}dt_{1}\,\int_{0}^{t_{1}}dt_{2}\int_{0}^{t_{2}}dt_{3}\big{\\{}\left[\bm{V}(t_{1}),\left[\bm{V}(t_{2}),\bm{V}(t_{3})\right]\right]+\left[\left[\bm{V}(t_{1}),\bm{V}(t_{2})\right],\bm{V}(t_{3})\right]\big{\\}},$
(87)

where $\bm{V}(t)=\cos(\Omega t)\big{(}e^{-it\omega_{b0}}|0\rangle\langle
b|+e^{it\omega_{b0}}|b\rangle\langle 0|\big{)}$. Here we keep $\bm{U}_{I}(t)$
up to $O(\mu^{3})$. By direct calculation the evolution operator in the
interaction picture is obtained as

$\displaystyle\bm{U}_{I}(t,0)=e^{x(t)|0\rangle\langle b|+y(t)|b\rangle\langle
0|+h(t)|0\rangle\langle 0|-h(t)|b\rangle\langle b|},$ (88)

where

$\displaystyle x(t)\equiv g(t)+a(t)+b(t)+c(t)+d(t)+e(t)+f(t),$ $\displaystyle
y(t)\equiv-[g^{*}(t)+a^{*}(t)+b^{*}(t)-c(t)+d^{*}(t)+e^{*}(t)+f^{*}(t)],$
$\displaystyle
h(t)\equiv\dfrac{\mu^{2}}{\hbar^{2}}\dfrac{i\omega_{b0}}{\omega^{2}_{b0}-\Omega^{2}}\Big{[}\dfrac{\Omega\cos(\omega_{b0}t)\sin(\Omega
t)-\omega_{b0}\cos(\Omega
t)\sin(\omega_{b0}t)}{\omega^{2}_{b0}-\Omega^{2}}+\dfrac{2t\Omega+\sin(2t\Omega)}{4\Omega}\Big{]},$

and the functions $a(t)$, $b(t)$, $c(t)$, $d(t)$, $e(t)$, $f(t)$, and $g(t)$
are defined as follows:

$\displaystyle a(t)\equiv$
$\displaystyle\dfrac{\mu^{3}\omega_{b0}}{12\hbar^{3}\Omega(\Omega^{2}-\omega^{2}_{b0})}\Big{[}\dfrac{4\Omega^{2}\sin(t\Omega)e^{-it\omega_{b0}}\left(t\Omega^{2}+2i\omega_{b0}-t\omega^{2}_{b0}\right)-4\Omega\left(\Omega^{2}+\omega^{2}_{b0}\right)}{(\Omega^{2}-\omega^{2}_{b0})^{2}}$
$\displaystyle+\dfrac{e^{-it\omega_{b0}}\left(\Omega\cos(\Omega
t)+i\omega_{b0}\sin(\Omega
t)\right)-\Omega}{\omega_{b0}^{2}-\Omega^{2}}+\dfrac{e^{-it\omega_{b0}}\big{(}3\Omega\cos(3\Omega
t)+i\omega_{b0}\sin(3\Omega t)\big{)}-3\Omega}{\omega_{b0}^{2}-9\Omega^{2}}$
$\displaystyle+\dfrac{4\Omega
e^{-it\omega_{b0}}\cos(t\Omega)\left(\Omega^{2}-it\Omega^{2}\omega_{b0}+\omega^{2}_{b0}+it\omega_{b0}^{3}\right)}{(\Omega^{2}-\omega^{2}_{b0})^{2}}\Big{]},$
(90) $\displaystyle b(t)\equiv$
$\displaystyle-\dfrac{2\mu^{3}\omega_{b0}}{3\hbar^{3}(\Omega^{2}-\omega^{2}_{b0})^{2}}\times\dfrac{e^{-2it\omega_{b0}}}{8\Omega(\omega^{2}_{b0}-\Omega^{2})}\Big{[}\Omega\cos(2t\Omega)\left(\Omega^{2}(1+e^{2it\omega_{b0}})-\omega^{2}_{b0}(e^{2it\omega_{b0}}-1)\right)$
$\displaystyle-\Omega^{3}\left(e^{2it\omega_{b0}}(1+2it\omega_{b0})+1\right)+\Omega\omega^{2}_{b0}\left(1+e^{2it\omega_{b0}}(2it\omega_{b0}-1)\right)$
$\displaystyle+i\omega_{b0}\sin(2t\Omega)\left(\Omega^{2}(2-e^{2it\omega_{b0}})+\omega_{b0}^{2}e^{2it\omega_{b0}}\right)\Big{]},$
$\displaystyle c(t)\equiv$
$\displaystyle\dfrac{2i\mu^{3}\omega_{b0}^{2}}{3\hbar^{3}(\Omega^{2}-\omega^{2}_{b0})^{2}}\Big{[}\dfrac{\Omega\cos(t\omega_{b0})\sin(t\Omega)-\omega_{b0}\cos(t\Omega)\sin(t\omega_{b0})}{\Omega^{2}-\omega^{2}_{b0}}-\dfrac{2t\Omega+\sin(2t\Omega)}{4\Omega}\Big{]},$
$\displaystyle d(t)\equiv$
$\displaystyle\dfrac{i\mu^{3}}{12\hbar^{3}}\dfrac{\Omega^{2}-2\omega^{2}_{b0}}{(\Omega^{2}-\omega^{2}_{b0})^{3}}\left[-i\omega_{b0}+e^{it\omega_{b0}}\big{(}i\omega_{b0}\cos(t\Omega)+\Omega\sin(t\Omega)\big{)}\right],$
(91) $\displaystyle e(t)\equiv$
$\displaystyle\dfrac{i\mu^{3}}{12\hbar^{3}\Omega(\Omega^{2}-\omega^{2}_{b0})^{2}(9\Omega^{4}-10\Omega^{2}\omega^{2}_{b0}+\omega^{4}_{b0})}\Big{[}i\Omega\omega_{b0}(9\Omega^{2}-\omega^{2}_{b0})(\omega_{2}+4\omega^{2}_{b0})$
$\displaystyle+\sin(t\Omega)\Bigl{(}-\omega_{b0}^{6}-2\Omega^{2}\omega_{b0}^{4}(1+it\omega_{b0})+4\Omega^{4}\omega_{b0}^{2}(13+5it\omega_{b0})-9\Omega^{6}(1+2it\omega_{b0})$
$\displaystyle-\omega_{b0}^{2}(\omega_{b0}^{4}-3\Omega^{4}+2\Omega^{2}\omega_{b0}^{2})\Bigr{)}+2i\Omega\omega^{3}_{b0}(\Omega^{2}-\omega^{2}_{b0})-2i\Omega\omega^{3}_{b0}\cos(3t\Omega)(\Omega^{2}-\omega^{2}_{b0})$
$\displaystyle+e^{-it\omega_{b0}}\bigl{(}-\Omega\omega_{b0}\cos(t\Omega)(9\Omega^{2}-\omega^{2}_{b0})(-2\omega^{2}_{b0}(-2i+t\omega_{b0})+\Omega^{2}(i+2t\omega_{b0})\bigr{)}\Big{]},$
(92) $\displaystyle f(t)\equiv$
$\displaystyle\dfrac{-\mu^{3}e^{-it\omega_{b0}}}{24\hbar^{3}(\Omega^{2}-\omega^{2}_{b0})^{2}(9\Omega^{4}-10\Omega^{2}\omega_{b0}^{2}+\omega_{b0}^{4})}\Big{[}\cos(3t\Omega)\left(\omega_{b0}^{5}-3\Omega^{4}\omega_{b0}+2\Omega^{2}\omega_{b0}^{3}\right)$
$\displaystyle+2\Omega\Bigl{(}-8\Omega\omega_{b0}e^{it\omega_{b0}}(\omega_{b0}^{2}-3\Omega^{2})-i\Omega^{2}\sin(t\Omega)(e^{2it\omega_{b0}}-1)(\omega_{b0}^{2}-9\Omega^{2})$
$\displaystyle+2i\omega_{b0}^{2}\sin(3t\Omega)(\omega_{b0}^{2}-\Omega^{2})\Bigr{)}-\omega_{b0}\cos(t\Omega)\,(\omega_{b0}^{2}-9\omega_{2})\bigl{(}\omega_{b0}^{2}-\Omega^{2}(3+2e^{2it\omega_{b0}})\bigr{)}\Big{]},$
$\displaystyle g(t)\equiv$
$\displaystyle\dfrac{-i\mu}{\hbar}\left[-i\omega_{b0}+ie^{-it\omega_{b0}}\big{(}\omega_{b0}\cos(t\Omega)+i\Omega\sin(t\Omega)\big{)}\right].$
(93)

Now we want to calculate the evolution operator in the Schrödinger picture by
applying the Baker-Campbell-Hausdorff (BCH) formula to the RHS of Eq. (82). By
keeping only first twelve terms of this formula, we obtain the operator
$\mathpzc{E}$ as [58]

$\displaystyle\mathpzc{E}(t)=$
$\displaystyle\log[e^{\Theta(t)}e^{\bm{\Lambda}(t)}]\approx\Theta(t)+\bm{\Lambda}(t)+\dfrac{1}{2}\left[\Theta(t),\bm{\Lambda}(t)\right]$
$\displaystyle+\dfrac{1}{12}\left(\left[\Theta(t),\left[\Theta(t),\bm{\Lambda}(t)\right]\right]+\left[\bm{\Lambda}(t),\left[\bm{\Lambda}(t),\Theta(t)\right]\right]\right)-\dfrac{1}{24}\left[\bm{\Lambda}(t),\left[\Theta(t),\left[\Theta(t),\bm{\Lambda}(t)\right]\right]\right]$
$\displaystyle-\dfrac{1}{720}\left(\left[\bm{\Lambda}(t),\left[\bm{\Lambda}(t),\left[\bm{\Lambda}(t),\left[\bm{\Lambda}(t),\Theta(t)\right]\right]\right]\right]+\left[\Theta(t),\left[\Theta(t),\left[\Theta(t),\left[\Theta(t),\bm{\Lambda}(t)\right]\right]\right]\right]\right)$
$\displaystyle+\dfrac{1}{360}\left(\left[\Theta(t),\left[\bm{\Lambda}(t),\left[\bm{\Lambda}(t),\left[\bm{\Lambda}(t),\Theta(t)\right]\right]\right]\right]+\left[\bm{\Lambda}(t),\left[\Theta(t),\left[\Theta(t),\left[\Theta(t),\bm{\Lambda}(t)\right]\right]\right]\right]\right)$
$\displaystyle+\dfrac{1}{120}\left(\left[\bm{\Lambda}(t),\left[\Theta(t),\left[\bm{\Lambda}(t),\left[\Theta(t),\bm{\Lambda}(t)\right]\right]\right]\right]+\left[\Theta(t),\left[\bm{\Lambda}(t),\left[\Theta(t),\left[\bm{\Lambda}(t),\Theta(t)\right]\right]\right]\right]\right)+\ldots,$
(94)

where

$\displaystyle\Theta(t)\equiv(-it/\hbar)H_{0}=-it\big{(}\omega_{0}|0\rangle\langle
0|+\omega_{b}|b\rangle\langle
b|+\textstyle{\sum_{j=1}^{n-2}}\omega_{j}|j\rangle\langle j|\big{)},$
$\displaystyle\bm{\Lambda}(t)\equiv x(t)|0\rangle\langle
b|+y(t)|b\rangle\langle 0|+h(t)|0\rangle\langle 0|-h(t)|b\rangle\langle b|,$
(95) $\displaystyle\omega_{j}\equiv\epsilon_{j}/\hbar,$ $\displaystyle
U_{\mathsf{S}}(t,0)=e^{\mathpzc{E}(t)}.$ (96)

Equation (95) is related for $n$-level atoms. In the following a general
expression for the evolution operator in the Schrödinger picture for an
$n$-level system which its two levels are coupled together via an external
field is obtained,

$\displaystyle\mathpzc{E}(t)=u_{00}(t)|0\rangle\langle
0|+u_{0b}(t)|0\rangle\langle b|+u_{b0}(t)|b\rangle\langle
0|+u_{bb}(t)|b\rangle\langle
b|-it\textstyle{\sum_{j=1}^{n-2}}\omega_{j}|j\rangle\langle j|.$ (97)

where the functions $u_{00}(t)$, $u_{0b}(t)$, $u_{b0}(t)$ and $u_{bb}(t)$ have
been obtained by BCH expansion as follows:

$\displaystyle u_{00}(t)\equiv$ $\displaystyle-
it\omega_{0}+h(t)+\dfrac{i}{6}tx(t)y(t)\omega_{b0}-\dfrac{i}{90}tx(t)y(t)\omega_{b0}\bigl{(}h(t)^{2}+x(t)y(t)\bigr{)}$
$\displaystyle+\dfrac{i}{180}t^{3}\omega_{b0}^{3}x(t)y(t)+\dfrac{1}{30}t^{2}\omega_{b0}^{2}h(t)x(t)y(t)$
$\displaystyle u_{0b}(t)\equiv$ $\displaystyle
x(t)+\dfrac{i}{2}tx(t)\omega_{b0}-\dfrac{1}{12}\omega_{b0}^{2}t^{2}x(t)-\dfrac{i}{6}tx(t)h(t)\omega_{b0}+\dfrac{1}{12}h(t)\omega_{b0}^{2}t^{2}x(t)$
$\displaystyle+\dfrac{i}{90}tx(t)h(t)\omega_{b0}\bigl{(}h^{2}(t)+x(t)y(t)\bigr{)}-\dfrac{1}{720}t^{4}\omega_{b0}^{4}x(t)+\dfrac{i}{60}t^{3}\omega_{b0}^{3}x(t)h(t)$
$\displaystyle-\dfrac{i}{180}t^{3}\omega_{b0}^{3}x(t)h(t)-\dfrac{1}{30}t^{2}\omega_{b0}^{2}x(t)h^{2}(t)+\dfrac{1}{90}t^{2}\omega_{b0}^{2}x(t)\bigl{(}h^{2}(t)+x(t)y(t)\bigr{)}$
$\displaystyle u_{b0}(t)\equiv$ $\displaystyle
y(t)-\dfrac{i}{2}ty(t)\omega_{b0}-\dfrac{1}{12}\omega_{b0}^{2}t^{2}y(t)-\dfrac{i}{6}ty(t)h(t)\omega_{b0}-\dfrac{1}{12}h(t)\omega_{b0}^{2}t^{2}y(t)$
$\displaystyle+\dfrac{i}{90}ty(t)h(t)\omega_{b0}\bigl{(}h^{2}(t)+x(t)y(t)\bigr{)}-\dfrac{1}{720}t^{4}\omega_{b0}^{4}y(t)+\dfrac{i}{60}t^{3}\omega_{b0}^{3}y(t)h(t)$
$\displaystyle-\dfrac{i}{180}t^{3}\omega_{b0}^{3}y(t)h(t)-\dfrac{1}{30}t^{2}\omega_{b0}^{2}y(t)h^{2}(t)+\dfrac{1}{90}t^{2}\omega_{b0}^{2}y(t)\bigl{(}h^{2}(t)+x(t)y(t)\bigr{)}$
$\displaystyle u_{bb}(t)\equiv$ $\displaystyle-
it\omega_{\mathsf{B}}-h(t)-\dfrac{i}{6}tx(t)y(t)\omega_{b0}+\dfrac{i}{90}tx(t)y(t)\omega_{b0}\bigl{(}h(t)^{2}+x(t)y(t)\bigr{)}$
$\displaystyle-\dfrac{i}{180}t^{3}\omega_{b0}^{3}x(t)y(t)-\dfrac{1}{30}t^{2}\omega_{b0}^{2}h(t)x(t)y(t).$
(98)

Thus, from Eq. (81) the Floquet Hamiltonian is

$\displaystyle\bar{H}=(i\hbar/\tau)\mathpzc{E}(\tau),$ (99)

from which we can also obtain the eigensystem
$\\{\bar{\epsilon},|\bar{\epsilon}\rangle\\}$. By having the Floquet
Hamiltonian and using Eq. (51) with $t_{0}=0$ the periodic operator
$P(t,t_{0})$ is as follows:

$\displaystyle P(t,0)=U_{\mathsf{S}}(t,0)e^{it\bar{H}/\hbar},$ (100)

where $U_{\mathsf{S}}(t,0)$ is obtained by solving Eq. (80) numerically. To
construct the master equation, the next step is to calculate the correlation
functions of the baths. In doing so we should rewrite the interaction
Hamiltonian between the system and the hot and the cold baths [Eqs. (21) and
(22)] as in the following:

$\displaystyle H_{\mathsf{S}\mathsf{B}_{\mathsf{h}}}$
$\displaystyle=\sigma_{x}^{(\mathsf{h})}\otimes\textstyle{\sum_{k}}f_{k}(\hat{a}_{k}+\hat{a}_{k}^{\dagger})+\sigma_{y}^{(\mathsf{h})}\otimes\textstyle{\sum_{k}}if_{k}(\hat{a}_{k}-\hat{a}_{k}^{\dagger}),$
(101) $\displaystyle H_{\mathsf{S}\mathsf{B}_{\mathsf{c}}}$
$\displaystyle=\sigma_{x}^{(\mathsf{c})}\otimes\textstyle{\sum_{q}}g_{q}(\hat{b}_{q}+\hat{b}_{q}^{\dagger})+\sigma_{y}^{(\mathsf{c})}\otimes\textstyle{\sum_{q}}ig_{q}(\hat{b}_{q}-\hat{b}_{q}^{\dagger}),$
(102)

where

$\displaystyle\sigma_{x}^{(\mathsf{h})}\equiv\dfrac{1}{2}\left(|0\rangle\langle
1|+|1\rangle\langle 0|\right)\qquad
B_{1}^{(\mathsf{h})}=\textstyle{\sum_{k}}f_{k}\big{(}\hat{a}_{k}+\hat{a}_{k}^{\dagger}\big{)},$
$\displaystyle\sigma_{y}^{(\mathsf{h})}\equiv\dfrac{i}{2}\left(|0\rangle\langle
1|-|1\rangle\langle 0|\right)\qquad
B_{2}^{(\mathsf{h})}=\textstyle{\sum_{k}}if_{k}\big{(}\hat{a}_{k}-\hat{a}_{k}^{\dagger}\big{)},$
(103)
$\displaystyle\sigma_{x}^{(\mathsf{c})}\equiv\dfrac{1}{2}\left(|b\rangle\langle
1|+|1\rangle\langle b|\right)\qquad
B_{1}^{(\mathsf{c})}=\textstyle{\sum_{q}}g_{q}\big{(}\hat{b}_{q}+\hat{b}_{q}^{\dagger}\big{)},$
$\displaystyle\sigma_{y}^{(\mathsf{c})}\equiv\dfrac{i}{2}\left(|b\rangle\langle
1|-|1\rangle\langle b|\right)\qquad
B_{2}^{(\mathsf{c})}=\textstyle{\sum_{q}}ig_{q}\big{(}\hat{b}_{q}-\hat{b}_{q}^{\dagger}\big{)}.$
(104)

According to Eq. (60), we also need to calculate the bath operators in the
interaction picture,
$\bm{B}(t)=e^{iH_{\mathsf{B}}t/\hbar}Be^{-iH_{\mathsf{B}}t/\hbar}$, which are
obtained as

$\displaystyle\bm{B}_{1}^{(\mathsf{h})}(t)=\textstyle{\sum_{k}}f_{k}\big{(}e^{-it\zeta_{k}}\hat{a}_{k}+e^{it\zeta_{k}}\hat{a}_{k}^{\dagger}\big{)},$
$\displaystyle\bm{B}_{2}^{(\mathsf{h})}(t)=\textstyle{\sum_{k}}if_{k}\big{(}e^{-it\zeta_{k}}\hat{a}_{k}-e^{it\zeta_{k}}\hat{a}_{k}^{\dagger}\big{)},$
(105)
$\displaystyle\bm{B}_{1}^{(\mathsf{c})}(t)=\textstyle{\sum_{q}}g_{q}\big{(}e^{-it\nu_{q}}\hat{b}_{q}+e^{it\nu_{q}}\hat{b}_{q}^{\dagger}\big{)},$
$\displaystyle\bm{B}_{2}^{(\mathsf{c})}(t)=\textstyle{\sum_{q}}ig_{q}\big{(}e^{-it\nu_{q}}\hat{b}_{q}-e^{it\nu_{q}}\hat{b}_{q}^{\dagger}\big{)}.$
(106)

All $\Gamma_{ij}(x)$ coefficients can be obtained similarly. For example,

$\displaystyle\Gamma^{(\mathsf{h})}_{11}(x)$
$\displaystyle=\textstyle{\int_{0}^{\infty}}ds\,e^{-isx}\big{\langle}\bm{B}^{(\mathsf{h})\dagger}_{1}(s)\bm{B}^{(\mathsf{h})}_{1}(0)\big{\rangle}=\int_{0}^{\infty}ds\,e^{-isx}\textstyle{\sum_{k}}f_{k}^{2}\left[\bar{n}(\zeta_{k},\beta)e^{is\zeta_{k}}+[\bar{n}(\zeta_{k},\beta)+1]e^{-is\zeta_{k}}\right]$
$\displaystyle=\textstyle{\sum_{k}}f_{k}^{2}\left\\{\bar{n}(\zeta_{k},\beta)\int_{0}^{\infty}ds\,e^{-i(x-\zeta_{k})s}+[\bar{n}(\zeta_{k},\beta)+1]\int_{0}^{\infty}ds\,e^{-i(x+\zeta_{k})s}\right\\}$
$\displaystyle=\textstyle{\int_{\infty}^{\infty}}d\nu\,\textstyle{\sum_{k}}\delta(\nu-\zeta_{k})f^{2}_{k}\left\\{\bar{n}(\nu,\beta)\int_{0}^{\infty}ds\,e^{-i(x-\nu)s}+[\bar{n}(\nu,\beta)+1]\int_{0}^{\infty}ds\,e^{-i(x+\nu)s}\right\\}$
$\displaystyle=\textstyle{\int_{-\infty}^{\infty}}d\nu\,J^{(\mathsf{h})}(\nu)\left\\{\bar{n}(\nu,\beta)\left[\pi\delta(x-\nu)-i\mathbbmss{P}\dfrac{1}{x-\nu}\right]+[\bar{n}(\nu,\beta)+1]\left[\pi\delta(x+\nu)-i\mathbbmss{P}\dfrac{1}{x+\nu}\right]\right\\}$
$\displaystyle=\pi\bar{n}(x,\beta)\big{[}J^{(\mathsf{h})}(x)-J^{(\mathsf{h})}(-x)\big{]}-i\mathbbmss{P}\textstyle{\int_{-\infty}^{\infty}}d\nu\,J^{(\mathsf{h})}(\nu)\left[\dfrac{\bar{n}(\nu,\beta)}{x-\nu}+\dfrac{\bar{n}(\nu,\beta)+1}{x+\nu}\right],$
(107)

where we have introduced the spectral density function of the hot bath as

$\displaystyle
J^{(\mathsf{h})}(x)=\textstyle{\sum_{k}}f^{2}_{k}\delta(x-\zeta_{k}),$ (108)

and also we have used the following formula:

$\displaystyle\textstyle{\int_{0}^{\infty}}ds\,e^{-i(x+\nu)s}=\pi\delta(x+\nu)-i\mathbbmss{P}\dfrac{1}{x+\nu},$
(109)

with $\mathbbmss{P}$ being the Cauchy principal value. Moreover, the following
relations for the bosonic creation and annihilation operators have been used:

$\displaystyle\langle\hat{a}^{\dagger}_{k}\hat{a}_{l}\rangle=\delta_{kl}\,\bar{n}(\zeta_{k},\beta),$
(110)
$\displaystyle\langle\hat{a}_{l}\hat{a}^{\dagger}_{k}\rangle=\delta_{kl}\,[\bar{n}(\zeta_{k},\beta)+1],$
(111) $\displaystyle\langle\hat{a}_{k}\hat{a}_{l}\rangle=0,$ (112)
$\displaystyle\langle\hat{a}_{k}^{\dagger}\hat{a}_{l}^{\dagger}\rangle=0,$
(113)

where $\beta$ is the inverse temperature of the bath. Finally according to Eq.
(64) by separating the real part from the imaginary part of the bath
correlation function, we find the following expressions:

$\displaystyle\gamma^{(\mathsf{h})}_{11}(x)=2\pi\bar{n}(x,\beta)\big{[}J^{(\mathsf{h})}(x)-J^{(\mathsf{h})}(-x)\big{]},\qquad\xi^{(\mathsf{h})}_{11}(x)=-\mathbbmss{P}\int_{-\infty}^{\infty}d\nu\,J^{(\mathsf{h})}({\nu})\left[\dfrac{\bar{n}(\nu,\beta)}{x-\nu}+\dfrac{\bar{n}(\nu,\beta)+1}{x+\nu}\right],$
$\displaystyle\gamma^{(\mathsf{h})}_{12}(x)=2\mathbbmss{P}\int_{-\infty}^{\infty}d\nu
J^{(\mathsf{h})}({\nu})\left[\dfrac{\bar{n}(\nu,\beta)}{x-\nu}-\dfrac{\bar{n}(\nu,\beta)+1}{x+\nu}\right],\qquad\xi^{(\mathsf{h})}_{12}(x)=\pi\bar{n}(x,\beta)\left[J^{(\mathsf{h})}(x)+J^{(\mathsf{h})}(-x)\right],$
$\displaystyle\gamma^{(\mathsf{h})}_{21}(x)=2\mathbbmss{P}\int_{-\infty}^{\infty}d\nu\,J^{(\mathsf{h})}({\nu})\left[\dfrac{\bar{n}(\nu,\beta)+1}{x+\nu}-\dfrac{\bar{n}(\nu,\beta)}{x-\nu}\right],\qquad\xi^{(\mathsf{h})}_{21}(x)=-\pi\bar{n}(x,\beta)\left[J^{(\mathsf{h})}(x)+J^{(\mathsf{h})}(-x)\right],$
$\displaystyle\gamma^{(\mathsf{h})}_{22}(x)=2\pi\bar{n}(x,\beta)\big{[}J^{(\mathsf{h})}(x)-J^{(\mathsf{h})}(-x)\big{]},\qquad\xi^{(\mathsf{h})}_{22}(x)=-\mathbbmss{P}\int_{-\infty}^{\infty}d\nu\,J^{(\mathsf{h})}({\nu})\left[\dfrac{\bar{n}(\nu,\beta)}{x-\nu}+\dfrac{\bar{n}(\nu,\beta)+1}{x+\nu}\right].$
(114)

For specificity, we choose the Ohmic spectral density function for the baths,

$\displaystyle J(x)=J_{0}xe^{-x^{2}/\omega^{2}_{\mathrm{cutoff}}},$ (115)

where $J_{0}$ is a model-dependent and also frequency-independent constant and
$\omega_{\mathrm{cutoff}}$ is the bath cutoff frequency. Hence Eqs. (114)
yield

$\displaystyle\gamma_{ij}(x)=0,\qquad i\neq j$
$\displaystyle\gamma_{ii}(x)=\begin{cases}4\pi\bar{n}(x,\beta)J(x)&\text{,
$x\neq 0$}\\\ 4\pi J_{0}/(\hbar\beta)&\text{, $x=0$},\end{cases}$ (116)
$\displaystyle\xi_{ij}(x)=0\qquad i\neq j$
$\displaystyle\xi_{ii}(x)=-2\mathbbmss{P}\textstyle{\int_{0}^{\infty}}d\omega\,G(x,\omega),$
(117)

where

$\displaystyle G(x,\omega)\equiv\begin{cases}J_{0}\omega
e^{-\omega^{2}/\omega_{\mathrm{cutoff}}^{2}}\frac{(x-\omega)e^{\beta\hbar\omega}+x+\omega}{\left(e^{\beta\hbar\omega}-1\right)(x^{2}-\omega^{2})}&\text{,
$\omega\neq 0$}\\\ 2J_{0}/(\hbar\beta x)&\text{, $\omega=0$}.\end{cases}$
(118)

Figure 14: (a) Fidelity between the exact and approximate solutions (Magnus +
BCH expansions) of the evolution operator in the Schrödinger picture and in
the first period for the model discussed in Sec. II.3. (b) Fidelity of
$P(t,0)$ in two periods which demonstrates the periodic property of $P(t,0)$.
All quantities are in natural units.

######  Remark 2.

If rather than taking the external field as in Eq. (18), we had chosen the
field as

$\displaystyle V_{1}(t)=\mu\cos(\Omega t)\big{(}|1\rangle\langle
b|+|b\rangle\langle 1|\big{)},$ (119)

then Eqs. (84), (95), and (97) would be, respectively, replaced with

$\displaystyle{\bm{V}}_{I}(t)=\mu\cos(\Omega
t)\left(e^{it\omega_{1b}}|1\rangle\langle
b|+e^{-it\omega_{1b}}|b\rangle\langle 1|\right),$ (120)
$\displaystyle\bm{\Lambda}(t)=x(t)|b\rangle\langle 1|+h(t)|b\rangle\langle
b|-h(t)|1\rangle\langle 1|+y(t)|1\rangle\langle b|,$ (121)
$\displaystyle\mathpzc{E}(t)=v_{bb}(t)|b\rangle\langle
b|+v_{b1}(t)|b\rangle\langle 1|+v_{1b}(t)|1\rangle\langle
b|+v_{11}(t)|1\rangle\langle 1|-it\omega_{0}|0\rangle\langle
0|-it\textstyle{\sum_{j=2}^{n-1}}\omega_{j}|j\rangle\langle j|,$ (122)

where the time-dependent functions $v_{bb}(t)$, $v_{b1}(t)$, $v_{1b}(t)$, and
$v_{11}(t)$ have the same definitions as $u_{00}(t)$, $u_{0b}(t)$,
$u_{b0}(t)$, and $u_{bb}(t)$, respectively, after the replacements
$0\rightarrow b$ and $b\rightarrow 1$. Thus, the Floquet Hamiltonian
$\bar{H}$, the periodic operator $P(t,0)$, and the Lindblad operators would
similarly be determined through Eqs. (99), (100), and (55).

###  C.2 Benchmarking the results

As in Sec. II.3, here we assume that $H_{0}$ is given by
$(\epsilon_{0},\epsilon_{b},\epsilon_{1})=(0,2.5,3)$ and
$\\{|0\rangle=(1,0,0)^{T},|1\rangle=(0,1,0)^{T},|b\rangle=(0,0,1)^{T}\\}$. The
system is assumed to be initially in the state
$\varrho_{\mathsf{S}}(0)=|0\rangle_{\mathsf{S}}\langle 0|$, and is coupled to
two thermal baths with $\beta_{\mathsf{c}}/\beta_{\mathsf{h}}=30/4$. We also
assume $(\mu,\Omega=2\pi/\tau)=(0.1,0.25$).

In order to obtain the Floquet Hamiltonian we first need to employ the Magnus
+ BCH expansions to obtain the operator $U_{\mathsf{S}}(t,0)$ as an
exponential in the form of Eq. (96). Since these expansions are infinite
series, in practice we need to truncate them, which in turn renders the
obtained $U_{\mathsf{S}}(t,0)$, $\bar{H}$, and $P(t,0)$ to be approximate. As
a result, one should investigate the validity of these approximations. To
benchmark the approximate operators, here we have also calculated
$U_{\mathsf{S}}(t,0)$ numerically through the $4$th order Runge-Kutta method
[59]. Our figures-of-merit are the following fidelities [60]:

Table 2: Parameters of the master equation (11). Here we consider two cases for the external field: (i) the field couples $|0\rangle$ and $|b\rangle$ levels [Eq. (18)], and (ii) the field couples $|1\rangle$ and $|b\rangle$ states [Eq. (119)]. All matrix representations are in the $\\{|0\rangle,|1\rangle,|b\rangle\\}$ eigenbasis of $H_{0}$ (as introduced in Sec. II.3), and all quantities are in natural units. |  $V_{0}$: field couples $|0\rangle$ and $|b\rangle$ |  $V_{1}$: field couples $|1\rangle$ and $|b\rangle$
---|---|---
$\bar{H}$ | $\left(\begin{matrix}0.0317&0&-0.2939+0.3347i\\\ 0&3&0\\\ -0.2939-0.3347i&0&2.4683\end{matrix}\right)$ | $\left(\begin{matrix}0&0&0\\\ 0&2.9991&0.0250\\\ 0&0.0250&2.5009\end{matrix}\right)$
$\bar{\epsilon}$ | $\left\\{3,2.5472,-0.0472\right\\}$ | $\left\\{3.0003,2,4997,0\right\\}$
$\omega$ | $\left\\{\pm 0.4528,\pm 3.0472,\pm 2.5943,0\right\\}$ | $\left\\{\pm 0.5006,\pm 3.0003,\pm 2.4997,0\right\\}$
$q$ | $\left[-24,24\right]$ | $\left[-3,3\right]$
$H_{\mathrm{lamb}}^{(\mathrm{F}),(\mathrm{h})}$ | $\left(\begin{matrix}-0.0145&0&-0.0016+0.0018i\\\ 0&0.0166&0\\\ -0.0016-0.0018i&0&-0.0015\end{matrix}\right)$ | $\left(\begin{matrix}-0.0190&0&0\\\ 0&0.0162&0.0006\\\ 0&0.0006&0.0035\end{matrix}\right)$
$H_{\mathrm{lamb}}^{(\mathrm{F}),(\mathrm{c})}$ | $\left(\begin{matrix}-0.0073&0&0.0243-0.0277i\\\ 0&0.2392&0\\\ 0.0243+0.0277i&0&-0.2088\end{matrix}\right)$ | $\left(\begin{matrix}0&0&0\\\ 0&0.1912&0.0182\\\ 0&0.0182&-0.1711\end{matrix}\right)$

$\displaystyle\mathpzc{F}\left[U^{\mathrm{app}}_{\mathsf{S}}(t,0),U^{\mathrm{ex}}_{\mathsf{S}}(t,0)\right]=\dfrac{1}{3}\big{|}\mathrm{Tr}\big{[}U^{\mathrm{app}}_{\mathsf{S}}(t,0)U_{\mathsf{S}}^{\mathrm{ex}\,\dagger}(t,0)\big{]}\big{|},$
(123)
$\displaystyle\mathpzc{F}\left[P^{\mathrm{app}}(t,0),P^{\mathrm{app}\,{\dagger}}(t+\tau,\tau)\right]=\dfrac{1}{3}\big{|}\mathrm{Tr}\big{[}P^{\mathrm{app}}(t,0),P^{\mathrm{app}\,{\dagger}}(t+\tau,\tau)\big{]}\big{|},$
(124)

where $0\leqslant t\leqslant\tau$. The first quantity measures the fidelity
between the approximate and almost exact solutions for $U_{\mathsf{S}}(t,0)$.
The second quantity measures periodicity of the approximate $P(t,0)$ [obtained
according to Eq. (100) where we have used only the Magnus expansion in order
to calculate $U_{\mathsf{S}}(t,0)$], as expected in
$P(t+n\tau,n\tau)=P(t+n\tau,0)=P(t,0)$ and the discussion after Eq. (51), for
two periods. Figure 14 a shows that the fidelity of the approximate and
almost-exact unitaries $U_{\mathsf{S}}(t,0)$ is more than $0.97$ in
$0\leqslant t\leqslant\tau$. Figure 14 b demonstrates that the fidelity
between the periodic unitary operators $P(t,0)$ in two periods is more than
$0.96$. We remark that in numerical calculations of $S_{\alpha}(q,\omega)$, we
have offset all values less than $10^{-3}$ to $0$. With these favorable
results, we obtain the other relevant quantities of the model (the Floquet
Hamiltonian $\bar{H}$, the Lamb-shift Hamiltonian
$H^{(\mathrm{F})}_{\mathrm{lamb}}$, and the range of $q$s and $\omega$s) as in
Table 2. These results have then been used to obtain the results discussed in
Sec. II.3.

##  Appendix D Redfield master equation

In this section, we explain how the Redfield master equation for a multilevel
quantum system can be obtained, and in particular we focus on the Jaynes-
Cummings model discuss in this paper.

###  D.1 Time-independent system Hamiltonian: Redfield master equation

Consider a quantum system which is coupled to a radiation bath. The total
system Hamiltonian is described by Eq. (33), where we consider the time-
independent case with $H_{\mathsf{S}}=H_{0}$, and the bath is comprised of
bosons with wave vectors $\vec{k}$, frequencies $\nu_{k}$, and polarization
$\lambda$ [see Eq. (25)]. The interaction Hamiltonian is described as Eq.
(26), where the couplings between the system and the bath are given by

$\displaystyle
g_{k\lambda}^{(ij)}=-E_{k}\vec{\mu}_{ij}\cdot\hat{e}_{k\lambda},$ (125)

where $\hat{e}_{k\lambda}$ denotes the unit vector of $\lambda$ polarization
for mode $k$, and $\vec{\mu}_{ij}=e\langle i|\vec{r}|j\rangle$ is the
electric-dipole transition matrix elements, with $e$ being the electric charge
[49]. Here the quantity $E_{k}$ has the dimension of an electric field and is
defined as $E_{k}=\sqrt{\hbar\nu_{k}/(2\epsilon_{0}V)}$, where $V$ is the
quantization volume and $\epsilon_{0}$ is the vacuum permittivity. We stress
that in the derivation of the master equation we use Eq. (109) for the
integral $\textstyle{\int_{0}^{\infty}}ds\,e^{-i(x+\nu)s}$, including both
delta function and the Cauchy principal value terms, contributing to the
decoherence operators and the Lamb-shift Hamiltonian, respectively. We remind
the following relation for transforming the interaction picture back to the
Schrödinger picture:

$\displaystyle\bm{\varrho}_{\mathsf{S}}(t)=U_{\mathsf{S}}^{\dagger}(t,0)\varrho_{\mathsf{S}}(t)U_{\mathsf{S}}(t,0),$
(126)
$\displaystyle\dfrac{d\varrho_{\mathsf{S}}(t)}{dt}=-\dfrac{i}{\hbar}\left[H_{\mathsf{S}}(t),\varrho_{\mathsf{S}}(t)\right]+U_{\mathsf{S}}(t,0)\big{[}\dfrac{d\bm{\varrho}_{\mathsf{S}}(t)}{dt}\big{]}U_{\mathsf{S}}^{\dagger}(t,0),$
(127)

where $U_{\mathsf{S}}(t,0)$ is the unitary operator generated by
$H_{\mathsf{S}}(t)=H_{0}$. After doing some algebra, we obtain the following
Redfield master equation in the Schrödinger picture:

$\displaystyle\dfrac{d\varrho_{\mathsf{S}}(t)}{dt}=-\dfrac{i}{\hbar}\left[H_{\mathsf{S}},\varrho_{\mathsf{S}}(t)\right]-\dfrac{1}{6\pi
c^{3}\hbar\epsilon_{0}}\textstyle{\sum_{i}}\sum_{i^{\prime}j^{\prime}}\omega_{i^{\prime}j^{\prime}}^{3}\times$
$\displaystyle\times\big{[}(\vec{\mu}_{i^{\prime}i}\cdot\vec{\mu}_{i^{\prime}j^{\prime}})\bar{n}(\omega_{i^{\prime}j^{\prime}},\beta)\left(|i\rangle\langle
j^{\prime}|\varrho_{\mathsf{S}}(t)+\varrho_{\mathsf{S}}(t)|j^{\prime}\rangle\langle
i|\right)+(\vec{\mu}_{ij^{\prime}}\cdot\vec{\mu}_{i^{\prime}j^{\prime}})[\bar{n}(\omega_{i^{\prime}j^{\prime}},\beta)+1]\left(|i\rangle\langle
i^{\prime}|\varrho_{\mathsf{S}}(t)+\varrho_{\mathsf{S}}(t)|i^{\prime}\rangle\langle
i|\right)\big{]}$ $\displaystyle+\dfrac{1}{6\pi
c^{3}\hbar\epsilon_{0}}\textstyle{\sum_{ij}}\sum_{i^{\prime}j^{\prime}}\vec{\mu}_{ij}\cdot\vec{\mu}_{i^{\prime}j^{\prime}}\omega_{i^{\prime}j^{\prime}}^{3}\times$
$\displaystyle\times\big{[}\bar{n}(\omega_{i^{\prime}j^{\prime}},\beta)\left[|i\rangle\langle
j|\varrho_{\mathsf{S}}(t)|j^{\prime}\rangle\langle
i^{\prime}|+|i^{\prime}\rangle\langle
j^{\prime}|\varrho_{\mathsf{S}}(t)|j\rangle\langle
i|\right]+[\bar{n}(\omega_{i^{\prime}j^{\prime}},\beta)+1]\left[|j\rangle\langle
i|\varrho_{\mathsf{S}}(t)|i^{\prime}\rangle\langle
j^{\prime}|+|j^{\prime}\rangle\langle
i^{\prime}|\varrho_{\mathsf{S}}(t)|i\rangle\langle j|\right]\big{]}$
$\displaystyle-\dfrac{i}{\hbar}\sum_{ij}\sum_{i^{\prime}}\dfrac{\vec{\mu}_{ij}\cdot\vec{\mu}_{i^{\prime}j}}{6\epsilon_{0}\pi^{2}c^{3}}\big{[}\big{(}|i\rangle\langle
i^{\prime}|\varrho_{\mathsf{S}}(t)-\varrho_{\mathsf{S}}(t)|i^{\prime}\rangle\langle
i|\big{)}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}(\bar{n}(\nu,\beta)+1)}{\omega_{i^{\prime}j}-\nu}\big{]}$
$\displaystyle-\dfrac{i}{\hbar}\sum_{ij}\sum_{j^{\prime}}\dfrac{\vec{\mu}_{ij}\cdot\vec{\mu}_{ij^{\prime}}}{6\epsilon_{0}\pi^{2}c^{3}}\big{[}\big{(}\varrho_{\mathsf{S}}(t)|j^{\prime}\rangle\langle
j|-|j\rangle\langle
j^{\prime}|\varrho_{\mathsf{S}}(t)\big{)}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}\bar{n}(\nu,\beta)}{\omega_{ij^{\prime}}-\nu}\big{]}$
$\displaystyle+\dfrac{i}{\hbar}\sum_{ij}\sum_{i^{\prime}j^{\prime}}\dfrac{\vec{\mu}_{ij}\cdot\vec{\mu}_{i^{\prime}j^{\prime}}}{6\epsilon_{0}\pi^{2}c^{3}}\big{[}\big{(}|i\rangle\langle
j|\varrho_{\mathsf{S}}(t)|j^{\prime}\rangle\langle
i^{\prime}|-|i^{\prime}\rangle\langle
j^{\prime}|\varrho_{\mathsf{S}}(t)|j\rangle\langle
i|\big{)}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}\bar{n}(\nu,\beta)}{\omega_{i^{\prime}j^{\prime}}-\nu}\big{]}$
$\displaystyle+\dfrac{i}{\hbar}\sum_{ij}\sum_{i^{\prime}j^{\prime}}\dfrac{\vec{\mu}_{ij}\cdot\vec{\mu}_{i^{\prime}j^{\prime}}}{6\epsilon_{0}\pi^{2}c^{3}}\big{[}\big{(}|j^{\prime}\rangle\langle
i^{\prime}|\varrho_{\mathsf{S}}(t)|i\rangle\langle j|-|j\rangle\langle
i|\varrho_{\mathsf{S}}(t)|i^{\prime}\rangle\langle
j^{\prime}|\big{)}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}(\bar{n}(\nu,\beta)+1)}{\omega_{i^{\prime}j^{\prime}}-\nu}\big{]}.$
(128)

We have used this equation in Sec. III.1.

###  D.2 Time-periodic system Hamiltonian: Floquet-Redfield master equation

Now we consider a quantum system in contact with a bath and also driven by an
external periodic field. The model we use is similar to the one discussed in
the precious section. To derive the master equation for such a system system,
we need to combine both Redfield and Floquet approaches. According to Appendix
A, to derive the master equation we need $\bm{H}_{\mathsf{int}}(t)$. By using
Eqs. (47) and (51) we obtain

$\displaystyle\bm{H}_{\mathsf{int}}(t)$
$\displaystyle=-\textstyle{\sum_{i,j}}\sum_{k,\lambda}g_{k\lambda}^{(ij)}e^{i\bar{H}t/\hbar}P^{\dagger}(t,0)\sigma_{ij}P(t,0)e^{-i\bar{H}t/\hbar}\otimes
U_{\mathsf{B}}^{\dagger}(t,0)\hat{a}_{k,\lambda}U_{\mathsf{B}}(t,0)+\mathrm{h.c.},$
$\displaystyle=-\textstyle{\sum_{i,j}}\sum_{k,\lambda}\sum_{q}\sum_{\left\\{\omega\right\\}}g_{k\lambda}^{(ij)}e^{it(\omega+q\Omega-\nu_{k})}\sigma_{ij}(q,\omega)\otimes\hat{a}_{k,\lambda}+\mathrm{h.c.},$
(129)

where

$\displaystyle\sigma_{ij}\equiv|i\rangle\langle j|,$ (130) $\displaystyle
U_{\mathsf{B}}^{\dagger}(t,0)\hat{a}_{k,\lambda}U_{\mathsf{B}}(t,0)=e^{-i\nu_{k}t}\hat{a}_{k,\lambda},$
(131) $\displaystyle
P^{\dagger}(t,0)\sigma_{ij}P(t,0)=\textstyle{\sum_{q\in\mathbbmss{Z}}}\sigma_{ij}(q)e^{iq\Omega
t},\qquad\sigma_{ij}^{\dagger}(q)=\sigma_{ij}(-q)$ (132)
$\displaystyle\sigma_{ij}(q,\omega)=\textstyle{\sum_{\bar{\epsilon}-\bar{\epsilon}^{\prime}=\hbar\omega}}|\bar{\epsilon}\rangle\langle\bar{\epsilon}|\sigma_{ij}(q)|\bar{\epsilon}^{\prime}\rangle\langle\bar{\epsilon}^{\prime}|,\qquad\sigma_{ij}^{\dagger}(q,\omega)=\sigma_{ij}(-q,-\omega)$
(133)
$\displaystyle\sigma_{ij}(t)=\textstyle{\sum_{q}}\sum_{\omega}e^{it(\omega+q\Omega)}\sigma_{ij}(q,\omega).$
(134)

By inserting $\bm{H}_{\mathsf{int}}(t)$ [Eq. (129)] in Eq. (49) and after some
algebra, we obtain the master equation of the system in the interaction
picture,

$\displaystyle\dfrac{d\bm{\varrho}_{\mathsf{S}}(t)}{dt}=$
$\displaystyle-\dfrac{1}{6\pi
c^{3}\hbar\epsilon_{0}}\textstyle{\sum_{i,j}}\sum_{i^{\prime},j^{\prime}}\sum_{q,\omega}\sum_{q^{\prime},\omega^{\prime}}(\vec{\mu}_{ij}\cdot\vec{\mu}_{i^{\prime}j^{\prime}})\times$
$\displaystyle\Big{(}e^{it(\omega-\omega^{\prime}+(q-q^{\prime})\Omega)}\big{[}(\omega^{\prime}+q^{\prime}\Omega)^{3}[\bar{n}(\omega^{\prime}+q^{\prime}\Omega,\beta)+1]+\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}(\bar{n}(\nu,\beta)+1)}{\omega^{\prime}+q^{\prime}\Omega-\nu}\big{]}\sigma_{ij}(q,\omega)\sigma_{i^{\prime}j^{\prime}}^{\dagger}(q^{\prime},\omega^{\prime})\bm{\varrho}_{\mathsf{S}}(t)$
$\displaystyle+e^{-it(\omega-\omega^{\prime}+(q-q^{\prime})\Omega)}\big{[}(\omega^{\prime}+q^{\prime}\Omega)^{3}\bar{n}(\omega^{\prime}+q^{\prime}\Omega,\beta)-\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}\bar{n}(\nu,\beta)}{\omega^{\prime}+q^{\prime}\Omega-\nu}\big{]}\sigma_{ij}^{\dagger}(q,\omega)\sigma_{i^{\prime}j^{\prime}}(q^{\prime},\omega^{\prime})\bm{\varrho}_{\mathsf{S}}(t)$
$\displaystyle-e^{it(\omega-\omega^{\prime}+(q-q^{\prime})\Omega)}\big{[}(\omega^{\prime}+q^{\prime}\Omega)^{3}\bar{n}(\omega^{\prime}+q^{\prime}\Omega,\beta)+\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}\bar{n}(\nu,\beta)}{\omega^{\prime}+q^{\prime}\Omega-\nu}\big{]}\sigma_{ij}(q,\omega)\bm{\varrho}_{\mathsf{S}}(t)\sigma_{i^{\prime}j^{\prime}}^{\dagger}(q^{\prime},\omega^{\prime})$
$\displaystyle-e^{-it(\omega-\omega^{\prime}+(q-q^{\prime})\Omega)}\big{[}(\omega^{\prime}+q^{\prime}\Omega)^{3}[\bar{n}(\omega^{\prime}+q^{\prime}\Omega,\beta)+1]-\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}(\bar{n}(\nu,\beta)+1)}{\omega^{\prime}+q^{\prime}\Omega-\nu}\big{]}\sigma_{ij}^{\dagger}(q,\omega)\bm{\varrho}_{\mathsf{S}}(t)\sigma_{i^{\prime}j^{\prime}}(q^{\prime},\omega^{\prime})$
$\displaystyle-e^{-it(\omega-\omega^{\prime}+(q-q^{\prime})\Omega)}\big{[}(\omega^{\prime}+q^{\prime}\Omega)^{3}\bar{n}(\omega^{\prime}+q^{\prime}\Omega,\beta)-\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}\bar{n}(\nu,\beta)}{\omega^{\prime}+q^{\prime}\Omega-\nu}\big{]}\sigma_{i^{\prime}j^{\prime}}(q^{\prime},\omega^{\prime})\bm{\varrho}_{\mathsf{S}}(t)\sigma_{ij}^{\dagger}(q,\omega)$
$\displaystyle-e^{it(\omega-\omega^{\prime}+(q-q^{\prime})\Omega)}\big{[}(\omega^{\prime}+q^{\prime}\Omega)^{3}[\bar{n}(\omega^{\prime}+q^{\prime}\Omega,\beta)+1]+\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}(\bar{n}(\nu,\beta)+1)}{\omega^{\prime}+q^{\prime}\Omega-\nu}\big{]}\sigma_{i^{\prime}j^{\prime}}^{\dagger}(q^{\prime},\omega^{\prime})\bm{\varrho}_{\mathsf{S}}(t)\sigma_{ij}(q,\omega)$
$\displaystyle+e^{-it(\omega-\omega^{\prime}+(q-q^{\prime})\Omega)}\big{[}(\omega^{\prime}+q^{\prime}\Omega)^{3}[\bar{n}(\omega^{\prime}+q^{\prime}\Omega,\beta)+1]-\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}(\bar{n}(\nu,\beta)+1)}{\omega^{\prime}+q^{\prime}\Omega-\nu}\big{]}\bm{\varrho}_{\mathsf{S}}(t)\sigma_{i^{\prime}j^{\prime}}(q^{\prime},\omega^{\prime})\sigma_{ij}^{\dagger}(q,\omega)$
$\displaystyle+e^{it(\omega-\omega^{\prime}+(q-q^{\prime})\Omega)}\big{[}(\omega^{\prime}+q^{\prime}\Omega)^{3}\bar{n}(\omega^{\prime}+q^{\prime}\Omega,\beta)+\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}\bar{n}(\nu,\beta)}{\omega^{\prime}+q^{\prime}\Omega-\nu}\big{]}\bm{\varrho}_{\mathsf{S}}(t)\sigma_{i^{\prime}j^{\prime}}^{\dagger}(q^{\prime},\omega^{\prime})\sigma_{ij}(q,\omega)\Big{)}.$
(135)

To transform back into the Schrödinger picture, we apply Eq. (126), where now
$U_{\mathsf{S}}(t,0)$ is the unitary operator given by Eq. (51). Substituting
Eq. (135) on the RHS of Eq. (127) then yields

$\displaystyle\dfrac{d\varrho_{\mathsf{S}}(t)}{dt}=-\dfrac{i}{\hbar}\left[H_{\mathsf{S}}(t),\varrho_{\mathsf{S}}(t)\right]-\dfrac{1}{6\pi
c^{3}\hbar\epsilon_{0}}\textstyle{\sum_{ij}}\sum_{i^{\prime}j^{\prime}}\sum_{q\omega}\sum_{q^{\prime}\omega^{\prime}}(\vec{\mu}_{ij}\cdot\vec{\mu}_{i^{\prime}j^{\prime}})$
$\displaystyle\times\Big{[}e^{it(q-q^{\prime})\Omega}\big{[}(\omega^{\prime}+q^{\prime}\Omega)^{3}[\bar{n}(\omega^{\prime}+q^{\prime}\Omega,\beta)+1]+\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}(\bar{n}(\nu,\beta)+1)}{\omega^{\prime}+q^{\prime}\Omega-\nu}\big{]}P(t,0)\sigma_{ij}(q,\omega)\sigma_{i^{\prime}j^{\prime}}^{\dagger}(q^{\prime},\omega^{\prime})P^{\dagger}(t,0)\varrho_{\mathsf{S}}(t)$
$\displaystyle+e^{-it(q-q^{\prime})\Omega}\big{[}(\omega^{\prime}+q^{\prime}\Omega)^{3}\bar{n}(\omega^{\prime}+q^{\prime}\Omega,\beta)-\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}\bar{n}(\nu,\beta)}{\omega^{\prime}+q^{\prime}\Omega-\nu}\big{]}P(t,0)\sigma_{ij}^{\dagger}(q,\omega)\sigma_{i^{\prime}j^{\prime}}(q^{\prime},\omega^{\prime})P^{\dagger}(t,0)\varrho_{\mathsf{S}}(t)$
$\displaystyle-e^{it(q-q^{\prime})\Omega}\big{[}(\omega^{\prime}+q^{\prime}\Omega)^{3}\bar{n}(\omega^{\prime}+q^{\prime}\Omega,\beta)+\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}\bar{n}(\nu,\beta)}{\omega^{\prime}+q^{\prime}\Omega-\nu}\big{]}P(t,0)\sigma_{ij}(q,\omega)P^{\dagger}(t,0)\varrho_{\mathsf{S}}(t)P(t,0)\sigma_{i^{\prime}j^{\prime}}^{\dagger}(q^{\prime},\omega^{\prime})P^{\dagger}(t,0)$
$\displaystyle-e^{-it(q-q^{\prime})\Omega}\big{[}(\omega^{\prime}+q^{\prime}\Omega)^{3}[\bar{n}(\omega^{\prime}+q^{\prime}\Omega,\beta)+1]-\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}(\bar{n}(\nu,\beta)+1)}{\omega^{\prime}+q^{\prime}\Omega-\nu}\big{]}P(t,0)\sigma_{ij}^{\dagger}(q,\omega)P^{\dagger}(t,0)\varrho_{\mathsf{S}}(t)P(t,0)\sigma_{i^{\prime}j^{\prime}}(q^{\prime},\omega^{\prime})P^{\dagger}(t,0)$
$\displaystyle-e^{-it(q-q^{\prime})\Omega}\big{[}(\omega^{\prime}+q^{\prime}\Omega)^{3}\bar{n}(\omega^{\prime}+q^{\prime}\Omega,\beta)-\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}\bar{n}(\nu,\beta)}{\omega^{\prime}+q^{\prime}\Omega-\nu}\big{]}P(t,0)\sigma_{i^{\prime}j^{\prime}}(q^{\prime},\omega^{\prime})P^{\dagger}(t,0)\varrho_{\mathsf{S}}(t)P(t,0)\sigma_{ij}^{\dagger}(q,\omega)P^{\dagger}(t,0)$
$\displaystyle-e^{it(q-q^{\prime})\Omega}\big{[}(\omega^{\prime}+q^{\prime}\Omega)^{3}[\bar{n}(\omega^{\prime}+q^{\prime}\Omega,\beta)+1]+\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}(\bar{n}(\nu,\beta)+1)}{\omega^{\prime}+q^{\prime}\Omega-\nu}\big{]}P(t,0)\sigma_{i^{\prime}j^{\prime}}^{\dagger}(q^{\prime},\omega^{\prime})P^{\dagger}(t,0)\varrho_{\mathsf{S}}(t)P(t,0)\sigma_{ij}(q,\omega)P^{\dagger}(t,0)$
$\displaystyle+e^{-it(q-q^{\prime})\Omega}\big{[}(\omega^{\prime}+q^{\prime}\Omega)^{3}[\bar{n}(\omega^{\prime}+q^{\prime}\Omega,\beta)+1]-\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}(\bar{n}(\nu,\beta)+1)}{\omega^{\prime}+q^{\prime}\Omega-\nu}\big{]}\varrho_{\mathsf{S}}(t)P(t,0)\sigma_{i^{\prime}j^{\prime}}(q^{\prime},\omega^{\prime})\sigma_{ij}^{\dagger}(q,\omega)P^{\dagger}(t,0)$
$\displaystyle+e^{it(q-q^{\prime})\Omega}\big{[}(\omega^{\prime}+q^{\prime}\Omega)^{3}\bar{n}(\omega^{\prime}+q^{\prime}\Omega,\beta)+\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}\bar{n}(\nu,\beta)}{\omega^{\prime}+q^{\prime}\Omega-\nu}\big{]}\varrho_{\mathsf{S}}(t)P(t,0)\sigma_{i^{\prime}j^{\prime}}^{\dagger}(q^{\prime},\omega^{\prime})\sigma_{ij}(q,\omega)P^{\dagger}(t,0)\Big{]}.$
(136)

The above master equation contains some fast oscillating terms which may lead
to nonpositivity of the density matrices. To alleviate this problem, we remove
some fast oscillating terms by applying a partial secular approximation which
is given by the assumption $\Omega^{-1}\ll\mathpzc{T}_{R}$, where
$\mathpzc{T}_{R}$ is the relaxation time of the driven system. This
approximation neglects all terms where $q^{\prime}\neq q$. Note that this
approximation is different from that of Ref. [29]. As a result, we obtain the
following modified Floquet-Redfield master equation:

$\displaystyle\dfrac{d\varrho_{\mathsf{S}}(t)}{dt}=$
$\displaystyle-\dfrac{i}{\hbar}\left[H_{\mathsf{S}}(t),\varrho_{\mathsf{S}}(t)\right]-\dfrac{1}{6\pi
c^{3}\hbar\epsilon_{0}}\sum_{ij}\sum_{i^{\prime}j^{\prime}}\sum_{q}\sum_{\omega\omega^{\prime}}(\vec{\mu}_{ij}\cdot\vec{\mu}_{i^{\prime}j^{\prime}})\times$
$\displaystyle\Big{[}L_{ij}(q,\omega;t)L^{\dagger}_{i^{\prime}j^{\prime}}(q,\omega^{\prime};t)\varrho_{\mathsf{S}}(t)\big{[}(\omega^{\prime}+q\Omega)^{3}[\bar{n}(\omega^{\prime}+q\Omega,\beta)+1]+\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}(\bar{n}(\nu,\beta)+1)}{\omega^{\prime}+q\Omega-\nu}\big{]}$
$\displaystyle+L^{\dagger}_{ij}(q,\omega;t)L_{i^{\prime}j^{\prime}}(q,\omega^{\prime};t)\varrho_{\mathsf{S}}(t)\big{[}(\omega^{\prime}+q\Omega)^{3}\bar{n}(\omega^{\prime}+q\Omega,\beta)-\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}\bar{n}(\nu,\beta)}{\omega^{\prime}+q\Omega-\nu}\big{]}$
$\displaystyle-
L_{ij}(q,\omega;t)\varrho_{\mathsf{S}}(t)L^{\dagger}_{i^{\prime}j^{\prime}}(q,\omega^{\prime};t)\big{[}(\omega^{\prime}+q\Omega)^{3}\bar{n}(\omega^{\prime}+q\Omega,\beta)+\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}\bar{n}(\nu,\beta)}{\omega^{\prime}+q\Omega-\nu}\big{]}$
$\displaystyle-L^{\dagger}_{ij}(q,\omega;t)\varrho_{\mathsf{S}}(t)L_{i^{\prime}j^{\prime}}(q,\omega^{\prime};t)\big{[}(\omega^{\prime}+q\Omega)^{3}[\bar{n}(\omega^{\prime}+q\Omega,\beta)+1]-\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}(\bar{n}(\nu,\beta)+1)}{\omega^{\prime}+q\Omega-\nu}\big{]}$
$\displaystyle-
L_{i^{\prime}j^{\prime}}(q,\omega^{\prime};t)\varrho_{\mathsf{S}}(t)L^{\dagger}_{ij}(q,\omega;t)\big{[}(\omega^{\prime}+q\Omega)^{3}\bar{n}(\omega^{\prime}+q\Omega,\beta)-\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}\bar{n}(\nu,\beta)}{\omega^{\prime}+q\Omega-\nu}\big{]}$
$\displaystyle-L^{\dagger}_{i^{\prime}j^{\prime}}(q,\omega^{\prime};t)\varrho_{\mathsf{S}}(t)L_{ij}(q,\omega;t)\big{[}(\omega^{\prime}+q\Omega)^{3}[\bar{n}(\omega^{\prime}+q\Omega,\beta)+1]+\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}[\bar{n}(\nu,\beta)+1]}{\omega^{\prime}+q\Omega-\nu}\big{]}$
$\displaystyle+\varrho_{\mathsf{S}}(t)L_{i^{\prime}j^{\prime}}(q,\omega^{\prime};t)L^{\dagger}_{ij}(q,\omega;t)\big{[}(\omega^{\prime}+q\Omega)^{3}[\bar{n}(\omega^{\prime}+q\Omega,\beta)+1]-\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}[\bar{n}(\nu,\beta)+1]}{\omega^{\prime}+q\Omega-\nu}\big{]}$
$\displaystyle+\varrho_{\mathsf{S}}(t)L^{\dagger}_{i^{\prime}j^{\prime}}(q,\omega^{\prime};t)L_{ij}(q,\omega;t)\big{[}(\omega^{\prime}+q\Omega)^{3}\bar{n}(\omega^{\prime}+q\Omega,\beta)+\dfrac{i}{\pi}\mathbbmss{P}\int_{0}^{\infty}d\nu\,\dfrac{\nu^{3}\bar{n}(\nu,\beta)}{\omega^{\prime}+q\Omega-\nu}\big{]}\Big{]},$
(137)

where

$\displaystyle
L_{ij}(q,\omega;t)=P(t,0)\sigma_{ij}(q,\omega)P^{\dagger}(t,0).$ (138)

We have used this equation in Sec. III.2.

##  Appendix E Computing the Lamb-shift correction in the Redfield equation

In both Eqs. (128) and (137) we have included the Lamb-shift terms. These
terms may give rise to divergence. It has been known that by by introducing a
finite cutoff frequency $W<\infty$ one can make these corrections finite. In
fact, the physical reason for the introduction of the cutoff frequency is
related to the long-wavelength approximation ($\omega|\vec{r}|/c\ll 1$, where
$|\vec{r}|$ is of the order of the atom size) [47]. After introducing $W$ we
need the following integral for the computation of the Lamb-shift corrections:

$\displaystyle\mathbbmss{P}\int_{0}^{W}d\nu\,\dfrac{\nu^{3}[\bar{n}(\nu,\beta)+1]}{\omega-\nu}=\mathbbmss{P}\int_{0}^{W}d\nu\,\dfrac{\nu^{3}\bar{n}(\nu,\beta)}{\omega-\nu}+\mathbbmss{P}\int_{0}^{W}d\nu\,\dfrac{\nu^{3}}{\omega-\nu}.$
(139)

Note that the second term of the RHS can also be recast as

$\displaystyle\mathbbmss{P}\int_{0}^{W}d\nu\,\dfrac{\nu^{3}}{\omega-\nu}=-\int_{0}^{W}d\nu\,\nu^{2}-\omega\int_{0}^{W}d\nu\,\nu-\omega^{2}\int_{0}^{W}d\nu+\mathbbmss{P}\int_{0}^{W}d\nu\dfrac{\omega^{3}}{\omega-\nu}.$
(140)

A careful quantum-electrodynamical (QED) investigation reveals that the first
term of Eq. (140) is fully compensated and cancelled out by other terms
emerged from the QED treatment [48]. As in Appendix B, consider a charged
particle (an atom or a molecule) coupled to a radiation field. Based on Eq.
(78) in the “$\vec{E}\cdot\vec{r}$” representation, we have an extra term
$\varepsilon_{\mathrm{dipole}}$ corresponding to a dipole self-energy of the
particle. In particular, in the Jaynes–Cummings model by considering a cutoff
frequency this dipole self-energy term is simply a constant energy value and
thus it is usually discarded in the literature [61]. However, interestingly,
it has been known that this term contributes an extra term in the Lamb-shift
correction which is the same as the first term of the RHS of Eq. (140), now
with the opposite sign. That is, the dipole self-energy contribution in the
energy shift completely cancels out the effect of the first term of the RHS of
Eq. (140). Furthermore, as argued in Remark 1 of Appendix B, in the low-
intensity regime, the appearance of $-H_{I2}$ term in the Hamiltonian
contributes an opposite term in the Lamb-shift correction, which cancels out
the effect of the the second term of the RHS of Eq. (140).

######  Remark 3.

Note that for low-energy photons one can ignore the interaction between the
spin magnetic moment of the charged particle and the magnetic field of the
radiation field [47].

## References

  * Mohseni  _et al._ [2014] M. Mohseni, Y. Omar, G. S. Engel, and M. B. Plenio, eds.,  _Quantum Effects in Biology_ (Cambridge University Press, Cambridge, 2014).
  * Demtröder [2005] W. Demtröder,  _Molecular Physics: Theoretical Principles and Experimental Methods_ (Wiley-VCH, Weinheim, Germany, 2005).
  * Dexter [1953] D. L. Dexter, A theory of sensitized luminescence in solids, J. Chem. Phys. 21, 836 (1953).
  * Gerardot  _et al._ [2005] B. D. Gerardot, S. Strauf, M. J. A. de Dood, A. M. Bychkov, A. Badolato, K. Hennessy, E. L. Hu, D. Bouwmeester, and P. M. Petroff, Photon Statistics from Coupled Quantum Dots, Phys. Rev. Lett. 95, 137403 (2005).
  * Collini and Scholes [2009] E. Collini and G. D. Scholes, Coherent intrachain energy migration in a conjugated polymer at room temperature, Science 323, 369 (2009).
  * Herz  _et al._ [2004] L. M. Herz, C. Silva, A. C. Grimsdale, K. Müllen, and R. T. Phillips, Time-dependent energy transfer rates in a conjugated polymer guest-host system, Phys. Rev. B 70, 165207 (2004).
  * Cheng and Fleming [2009] Y.-C. Cheng and G. R. Fleming, Dynamics of light harvesting in photosynthesis, Ann. Rev. Phys. Chem. 60, 241 (2009).
  * Plenio and Huelga [2008] M. B. Plenio and S. F. Huelga, Dephasing-assisted transport: quantum networks and biomolecules, New J. Phys. 10, 113019 (2008).
  * Caruso  _et al._ [2009] F. Caruso, A. W. Chin, A. Datta, S. F. Huelga, and M. B. Plenio, Highly efficient energy excitation transfer in light-harvesting complexes: The fundamental role of noise-assisted transport, J. Chem. Phys. 131, 105106 (2009).
  * Wu  _et al._ [2013] J. Wu, R. J. Silbey, and J. Cao, Generic Mechanism of Optimal Energy Transfer Efficiency: A Scaling Theory of the Mean First-Passage Time in Exciton Systems, Phys. Rev. Lett. 110, 200402 (2013).
  * Svidzinsky  _et al._ [2011] A. A. Svidzinsky, K. E. Dorfman, and M. O. Scully, Enhancing photovoltaic power by Fano-induced coherence, Phys. Rev. A 84, 053818 (2011).
  * Dorfman  _et al._ [2011] K. E. Dorfman, M. B. Kim, and A. A. Svidzinsky, Increasing photocell power by quantum coherence induced by external source, Phys. Rev. A 84, 053829 (2011).
  * Caruso  _et al._ [2010] F. Caruso, A. W. Chin, A. Datta, S. F. Huelga, and M. B. Plenio, Entanglement and entangling power of the dynamics in light-harvesting complexes, Phys. Rev. A 81, 062346 (2010).
  * Campos Venuti and Zanardi [2011] L. Campos Venuti and P. Zanardi, Excitation transfer through open quantum networks: Three basic mechanisms, Phys. Rev. B 84, 134206 (2011).
  * Sarovar  _et al._ [2011] M. Sarovar, Y.-C. Cheng, and K. B. Whaley, Environmental correlation effects on excitation energy transfer in photosynthetic light harvesting, Phys. Rev. E 83, 011906 (2011).
  * Chen  _et al._ [2013] G.-Y. Chen, N. Lambert, C.-M. Li, Y.-N. Chen, and F. Nori, Rerouting excitation transfers in the Fenna-Matthews-Olson complex, Phys. Rev. E 88, 032120 (2013).
  * Irish  _et al._ [2014] E. K. Irish, R. Gómez-Bombarelli, and B. W. Lovett, Vibration-assisted resonance in photosynthetic excitation-energy transfer, Phys. Rev. A 90, 012510 (2014).
  * Roden and Whaley [2016] J. J. J. Roden and K. B. Whaley, Probability-current analysis of energy transport in open quantum systems, Phys. Rev. E 93, 012128 (2016).
  * Zwanzig [1960] R. Zwanzig, Ensemble method in the theory of irreversibility, J. Chem. Phys. 33, 1338 (1960).
  * Nakajima [1958] S. Nakajima, On quantum theory of transport phenomena: Steady diffusion, Prog. Theor. Phys. 20, 948 (1958).
  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2002).
  * Alipour  _et al._ [2020] S. Alipour, A. T. Rezakhani, A. P. Babu, K. Mølmer, M. Möttönen, and T. Ala-Nissila, Correlation-Picture Approach to Open-Quantum-System Dynamics, Phys. Rev. X 10, 041024 (2020).
  * Lindblad [1975] G. Lindblad, Completely positive maps and entropy inequalities, Commun. Math. Phys. 40, 147 (1975).
  * Lindblad [1976] G. Lindblad, On the generators of quantum dynamical semigroups, Commun. Math. Phys. 48, 119 (1976).
  * Gorini  _et al._ [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Completely positive dynamical semigroups of $N$-level systems, J. Math. Phys. 17, 821 (1976).
  * Palmieri  _et al._ [2009] B. Palmieri, D. Abramavicius, and S. Mukamel, Lindblad equations for strongly coupled populations and coherences in photosynthetic complexes, J. Chem. Phys. 130, 204512 (2009).
  * Jeske  _et al._ [2015] J. Jeske, D. J. Ing, M. B. Plenio, S. F. Huelga, and J. H. Cole, Bloch-Redfield equations for modeling light-harvesting complexes, J. Chem. Phys. 142, 064104 (2015).
  * Farina and Giovannetti [2019] D. Farina and V. Giovannetti, Open-quantum-system dynamics: Recovering positivity of the Redfield equation via the partial secular approximation, Phys. Rev. A 100, 012107 (2019).
  * Tscherbul and Brumer [2015] T. V. Tscherbul and P. Brumer, Partial secular Bloch-Redfield master equation for incoherent excitation of multilevel quantum systems, J. Chem. Phys. 142, 104107 (2015).
  * Hofer  _et al._ [2017] P. P. Hofer, M. Perarnau-Llobet, L. D. M. Miranda, G. Haack, R. Silva, J. B. Brask, and N. Brunner, Markovian master equations for quantum thermal machines: local versus global approach, New J. Phys. 19, 123037 (2017).
  * De Chiara  _et al._ [2018] G. De Chiara, G. Landi, A. Hewgill, B. Reid, A. Ferraro, A. J. Roncaglia, and M. Antezza, Reconciliation of quantum local master equations with thermodynamics, New J. Phys. 20, 113024 (2018).
  * Alicki  _et al._ [2006] R. Alicki, D. A. Lidar, and P. Zanardi, Internal consistency of fault-tolerant quantum error correction in light of rigorous derivations of the quantum Markovian limit, Phys. Rev. A 73, 052311 (2006).
  * Szczygielski [2014] K. Szczygielski, On the application of Floquet theorem in development of time-dependent Lindbladians, J. Math. Phys. 55, 083506 (2014).
  * Rebentrost  _et al._ [2009a] P. Rebentrost, R. Chakraborty, and A. Aspuru-Guzik, Non-Markovian quantum jumps in excitonic energy transfer, J. Chem. Phys. 131, 184102 (2009a).
  * Qin  _et al._ [2014] M. Qin, H. Z. Shen, X. L. Zhao, and X. X. Yi, Dynamics and quantumness of excitation energy transfer through a complex quantum network, Phys. Rev. E 90, 042140 (2014).
  * Rebentrost  _et al._ [2009b] P. Rebentrost, M. Mohseni, I. Kassal, S. Lloyd, and A. Aspuru-Guzik, Environment-assisted quantum transport, New J. Phys. 11, 033003 (2009b).
  * Rebentrost  _et al._ [2009c] P. Rebentrost, M. Mohseni, and A. Aspuru-Guzik, Role of quantum coherence and environmental fluctuations in chromophoric energy transport, J. Phys. Chem. B 113, 9942 (2009c).
  * Bukov  _et al._ [2015] M. Bukov, L. D’Alessio, and A. Polkovnikov, Universal high-frequency behavior of periodically driven systems: from dynamical stabilization to Floquet engineering, Adv. Phys. 64, 139 (2015).
  * Chicone [2006] C. Chicone,  _Ordinary Differential Equations with Applications_ (Springer, New York, 2006).
  * Scovil and Schulz-DuBois [1959] H. E. D. Scovil and E. O. Schulz-DuBois, Three-Level Masers as Heat Engines, Phys. Rev. Lett. 2, 262 (1959).
  * Geusic  _et al._ [1967] J. E. Geusic, E. O. Schulz-DuBios, and H. E. D. Scovil, Quantum equivalent of the Carnot cycle, Phys. Rev. 156, 343 (1967).
  * Henrich  _et al._ [2007] M. J. Henrich, G. Mahler, and M. Michel, Driven spin systems as quantum thermodynamic machines: Fundamental limits, Phys. Rev. E 75, 051118 (2007).
  * Gelbwaser-Klimovsky  _et al._ [2013] D. Gelbwaser-Klimovsky, R. Alicki, and G. Kurizki, Minimal universal quantum heat machine, Phys. Rev. E 87, 012140 (2013).
  * Eastham  _et al._ [2013] P. R. Eastham, A. O. Spracklen, and J. Keeling, Lindblad theory of dynamical decoherence of quantum-dot excitons, Phys. Rev. B 87, 195306 (2013).
  * Vogt  _et al._ [2013] N. Vogt, J. Jeske, and J. H. Cole, Stochastic Bloch-Redfield theory: Quantum jumps in a solid-state environment, Phys. Rev. B 88, 174514 (2013).
  * Bamba and Imoto [2016] M. Bamba and N. Imoto, Maxwell boundary conditions imply non-Lindblad master equation, Phys. Rev. A 94, 033802 (2016).
  * Cohen-Tannoudji  _et al._ [1989] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg,  _Photons and Atoms – Introduction to Quantum Electrodynamics_ (Wiley-VCH, Weinheim, Germany, 1989).
  * Cohen-Tannoudji  _et al._ [2004] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg,  _Atom-Photon Interactions – Basic Processes and Applications_ (Wiley-VCH, Weinheim, Germany, 2004).
  * Scully and Zubairy [1997] M. O. Scully and M. S. Zubairy,  _Quantum Optics_ (Cambridge University Press, Cambridge, 1997).
  * Carmichael [1993] H. Carmichael,  _An Open Systems Approach to Quantum Optics_ (Springer, Berlin, 1993).
  * Power  _et al._ [1959] E. A. Power, S. Zienau, and H. S. W. Massey, Coulomb gauge in non-relativistic quantum electro-dynamics and the shape of spectral lines, Phil. Trans. Roy. Soc. Lond. Ser. A Math. Phys. Sci. 251, 427 (1959).
  * Agarwal [1973] G. S. Agarwal, Rotating-wave approximation and spontaneous emission, Phys. Rev. A 7, 1195 (1973).
  * Zhang and Chen [2013] S.-Q. Zhang and Z.-D. Chen, Lamb shift of a lossy two-level emitter inside a leaky cavity: A universal description based on an equivalent dissipative model, Phys. Rev. A 87, 013804 (2013).