††thanks<EMAIL_ADDRESS>
# The origin of loose bound of the thermodynamic uncertainty relation in a
dissipative two-level quantum system

Davinder Singh  Korea Institute for Advanced Study, Seoul 02455, Korea
Changbong Hyeon  Korea Institute for Advanced Study, Seoul 02455, Korea

###### Abstract

The thermodynamic uncertainty relations (TURs), originally discovered for
classical systems, dictate the trade-off between dissipation and fluctuations
of irreversible current, specifying a minimal bound that constrains the two
quantities. In a series of effort to extend the relation to the one under more
generalized conditions, it has been noticed that the bound is less tight in
open quantum processes. To study the origin of the loose bounds, we consider
an external field-driven transition dynamics of two-level quantum system
weakly coupled to the bosonic bath as a model of open quantum system. The
model makes it explicit that the imaginary part of quantum coherence, which
contributes to dissipation to the environment, is responsible for loosening
the TUR bound by suppressing the relative fluctuations in the irreversible
current of transitions, whereas the real part of the coherence tightens it.
Our study offers a better understanding of how quantum nature affects the TUR
bound.

##  I Introduction

Thermodynamic uncertainty relations (TURs) offer quantitative ideas of how
much the irreversible current of a certain dynamical process and its
fluctuations constrain the entropy production, which makes the lower bound of
total entropy production dictated by the second law of thermodynamics more
precise [1, 2, 3, 4]. The relation was originally derived for classical
systems, particularly for continuous time Markov jump processes on networks
and for overdamped Langevin dynamics in nonequilibrium steady states (NESS)
[1, 2, 4], and it has been extended to the one at finite-time [5, 6, 7] and
discrete-time Markov processes [8]. The significance of TURs has been
illuminated in specific contexts of biological processes [9, 10, 11, 12, 13,
14, 15], heat engines [16, 17, 18, 19], and other dynamical processes [20, 21,
22, 23]. More recently, the universal bound of TUR has been used to infer the
dissipation rate from the fluctuating currents of dynamical processes [24].
The originally proposed relation has been extended to broader range of
nonequilibrium processes [25, 26, 5, 6, 27, 28, 29, 30], including those
generated in underdamped conditions [31], under periodic drives [8, 32], and
driven by velocity-dependent forces [33, 34]. However, these processes are
featured with the bounds less tight than that of the original TUR.

The TUR has been explored for quantum systems as well [35, 36, 37, 38, 39, 40,
41]. The TUR bound for an observable, such as the integrated charge current in
quantum transport, was found smaller than the value of the original TUR over a
certain parameter range [37, 38, 39, 40]. In these studies, the quantum nature
was suspected to enhance the precision of the current, and thus lowers the TUR
bound; however, expressed in the language of the transmission function of the
nonequilibrium Green’s function formalism, the physical origin of the loose
TUR bound was left unclear [37, 38, 39, 40].

Here we consider a simple example of open quantum process which allows us to
dissect the physical origin of the loose TUR bound. The model is a standard
two-level system (TLS) of ground ($\ket{g}$) and excited ($\ket{e}$) states
coupled to bosonic bath or radiation field maintained at temperature $T$. The
transitions between the two levels are stimulated by an external driving field
via dipole-electric field interaction [42, 43, 44] (Fig. 1). Without driving,
it is expected that the absorption ($\ket{g}\rightarrow\ket{e}$) and emission
($\ket{e}\rightarrow\ket{g}$) satisfy the detailed balance (DB) condition with
zero mean current ($\langle j\rangle=0$), and in this case the population
ratio of the ground and excited states must maintain the canonical
distribution, as a result of exchanging energy with the thermal bath via
system-bath coupling. On the other hand, the external field, perturbing the
system via the dipole-field coupling ($-\vec{d}\cdot\vec{E}$), generates a
transition current ($\langle j\rangle>0$) between the two levels, breaking the
DB condition, which generates heat current dissipating to the bath. The excess
number of transitions stimulated by the field is countered by the bath that
invokes time irreversible relaxation of the system to its NESS.

The expressions for the mean current of net transitions, current fluctuations,
and entropy production of the model enable us to clarify the detail of quantum
nature contributing to loosening the TUR bound.

Figure 1: A non-degenerate TLS coupled to bosonic bath, perturbed by the
incident field with the angular frequency $\omega$. The resonant frequency is
denoted by $\omega_{0}$. The difference between $\omega$ and $\omega_{0}$,
$\delta\omega=\omega-\omega_{0}$, is the detuning.

##  II Model

The total Hamiltonian in the presence of an external radiation field is given
by [43, 45]

$\displaystyle H(t)$ $\displaystyle=H_{S}+H_{\text{ext}}(t)+H_{B}+H_{SB}$ (1)

with

$\displaystyle H_{S}$ $\displaystyle=(\hbar/2)\omega_{0}\sigma_{z}$
$\displaystyle H_{\text{ext}}(t)$ $\displaystyle=-\vec{d}\cdot\vec{E}({\bf
r},t)$ $\displaystyle H_{B}$
$\displaystyle=\sum_{\textbf{k},\xi}\hbar\omega_{\textbf{k}}b_{\textbf{k},\xi}^{\dagger}b_{\textbf{k},\xi}$
$\displaystyle H_{SB}$
$\displaystyle=\sum_{\textbf{k},\xi}\hbar\left(g_{\textbf{k},\xi}^{\ast}b_{\textbf{k},\xi}^{\dagger}\sigma_{-}+g_{\textbf{k},\xi}b_{\textbf{k},\xi}\sigma_{+}\right).$

Here $\sigma$’s are the Pauli spin operators:
$\sigma_{z}=\ket{e}\bra{e}-\ket{g}\bra{g}$, $\sigma_{+}=\ket{e}\bra{g}$, and
$\sigma_{-}=\ket{g}\bra{e}$. $H_{\rm ext}(t)$ is defined by the dipole
operator $\vec{d}=\vec{d}_{eg}\sigma_{+}+\vec{d}_{ge}\sigma_{-}$ with the
transition dipole $\vec{d}_{eg}$ ($\vec{d}_{ge}$) from $\ket{g}$ to $\ket{e}$
(from $\ket{e}$ to $\ket{g}$), interacting with the electric field
$\vec{E}({\bf r},t)=\vec{\varepsilon}e^{-i({\bf k}\cdot{\bf r}-i\omega
t)}+\vec{\varepsilon}^{\ast}e^{i({\bf k}\cdot{\bf r}-i\omega t)}$, where ${\bf
k}$ is the wavevector interacting with the system at ${\bf r}$, and $\omega$
is the angular frequency of the field. In $H_{B}$ and $H_{SB}$, the summation
$\sum_{\textbf{k},\xi}$ runs over the wavevector k and polarization $\xi$. The
symbols, $b_{\textbf{k},\xi}^{\dagger}$ and $b_{\textbf{k},\xi}$ are the
creation and annihilation operators for the bath modes represented by the
harmonic oscillators with angular frequency $\omega_{k}$. The strength of
system-bath coupling is quantified by $g_{\textbf{k},\xi}\equiv-ie^{i{\bf
k}\cdot{\bf r}}\sqrt{2\pi\omega_{\bf k}/\hbar
V}\hat{\varepsilon}_{\textbf{k},\xi}\cdot\vec{d}_{eg}$, which is defined in
terms of the polarization vector of electric field
$\hat{\varepsilon}_{\textbf{k},\xi}$, quantization volume $V$.

The total density matrix $\rho_{\text{tot}}(t)$ undergoes a unitary evolution,
obeying the von Neumann equation,

$\displaystyle\frac{d\rho_{\text{tot}}(t)}{dt}=-\frac{i}{\hbar}\left[H(t),\rho_{\text{tot}}(t)\right].$
(2)

Next, tracing out the bath degrees of freedom from the total density matrix,
we obtain the evolution equation of reduced density matrix for the system
$\rho(t)=\text{Tr}_{B}(\rho_{\text{tot}}(t))$, by assuming that
$\rho_{\text{tot}}(t)$ can be described with the uncorrelated product state of
the system $\rho(t)$ and environment in equilibrium ($\rho_{B}^{eq}$), namely
$\rho_{\text{tot}}(t)\approx\rho(t)\otimes\rho_{B}^{eq}$.

We employ the Lindblad approach [43, 42] (see Appendix A for details) which
ensures that the density matrix is self-adjoint ($\rho_{eg}=\rho_{ge}^{*}$),
trace-preserving ($\rho_{ee}+\rho_{gg}=1$), and positive semi-definite, while
the system evolves and relaxes to the steady state. The evolution equation for
$\rho(t)$ reads

$\displaystyle\frac{d\rho(t)}{dt}$
$\displaystyle=-\frac{i}{\hbar}[H_{S}+H_{\text{ext}}(t),\rho(t)]+\mathcal{D}(\rho(t)).$
(3)

In evaluating the effect of $H_{\text{ext}}(t)$ on $\rho(t)$, we consider the
dipole approximation that the driving field is nearly constant over the
molecular scale ($e^{i{\bf k}\cdot{\bf r}}\simeq 1$) [45], such that
$\vec{E}(t)\simeq\vec{\varepsilon}e^{i\omega
t}+\vec{\varepsilon}^{\ast}e^{-i\omega t}$. Together with
$\vec{d}=\vec{d}_{eg}\sigma_{+}+\vec{d}_{ge}\sigma_{-}$, one obtains
$H_{\text{ext}}(t)=-\hbar\Omega(e^{i\omega t}+e^{-i\omega
t})(\sigma_{+}+\sigma_{-})$ with the _driving frequency_ ,
$\Omega=\vec{d}_{eg}\cdot\vec{\varepsilon}/\hbar=\vec{d}_{ge}\cdot\vec{\varepsilon}^{\ast}/\hbar$.
The range of the angular frequency $\omega$ relevant to our model is
$\omega=\omega_{0}+\delta\omega$ with $|\delta\omega|/\omega_{0}\ll 1$, where
$\delta\omega$ is the detuning that reflects the extent of off-resonance (see
Fig. 1). The time-irreversible dynamics of system, invoked by a weak system-
bath coupling ($\omega\gg\gamma$), is delineated by the Lindblad _dissipator_
(see Appendix A for the derivation),

$\displaystyle\mathcal{D}(\rho(t))$
$\displaystyle\equiv\gamma(\bar{n}+1)\left(\sigma_{-}\rho(t)\sigma_{+}-\frac{1}{2}\Big{\\{}\sigma_{+}\sigma_{-},\rho(t)\Big{\\}}\right)$
$\displaystyle+\gamma\bar{n}\left(\sigma_{+}\rho(t)\sigma_{-}-\frac{1}{2}\Big{\\{}\sigma_{-}\sigma_{+},\rho(t)\Big{\\}}\right),$
(4)

where $\bar{n}=(e^{\beta\hbar\omega_{0}}-1)^{-1}$ denotes the mean occupation
number of bosonic bath, and
$\\{\hat{A},\hat{B}\\}=\hat{A}\hat{B}+\hat{B}\hat{A}$ is the anti-commutator.
The parameter $\gamma$, involving the spontaneous emission frequency, sets the
timescale of relaxation of the system to the steady state, and $\gamma\bar{n}$
represents the rate of bath-induced transition. Thus, the first and second
terms of Eq. 4 arise from emissions and absorptions, respectively.

Eqs. 3 and 4 yield the evolution equations for the reduced density matrix
elements as follows

$\displaystyle\dot{\rho}_{ee}=$
$\displaystyle-\gamma(\bar{n}+1)\rho_{ee}-i\Omega\left(e^{i\omega
t}+e^{-i\omega t}\right)\rho_{eg}$ $\displaystyle+i\Omega\left(e^{i\omega
t}+e^{-i\omega t}\right)\rho_{ge}+\gamma\bar{n}\rho_{gg},$
$\displaystyle\dot{\rho}_{eg}=$ $\displaystyle-i\Omega\left(e^{i\omega
t}+e^{-i\omega
t}\right)\rho_{ee}-\left[i\omega_{0}+\frac{\gamma}{2}(2\bar{n}+1)\right]\rho_{eg}$
$\displaystyle+i\Omega\left(e^{i\omega t}+e^{-i\omega t}\right)\rho_{gg}$
$\displaystyle\dot{\rho}_{ge}=$ $\displaystyle i\Omega\left(e^{i\omega
t}+e^{-i\omega
t}\right)\rho_{ee}+\left[i\omega_{0}-\frac{\gamma}{2}(2\bar{n}+1)\right]\rho_{ge}$
$\displaystyle-i\Omega\left(e^{i\omega t}+e^{-i\omega t}\right)\rho_{gg}$
$\displaystyle\dot{\rho}_{gg}=$
$\displaystyle\gamma(\bar{n}+1)\rho_{ee}+i\Omega\left(e^{i\omega
t}+e^{-i\omega t}\right)\rho_{eg}$ $\displaystyle-i\Omega\left(e^{i\omega
t}+e^{-i\omega t}\right)\rho_{ge}-\gamma\bar{n}\rho_{gg}$ (5)

Transformation of the density matrix into the one in the rotating frame,
$\rho_{ee}\rightarrow\tilde{\rho}_{ee}$,
$\rho_{gg}\rightarrow\tilde{\rho}_{gg}$, and
$\rho_{eg}\rightarrow\tilde{\rho}_{eg}e^{-i\omega t}$ (see Appendix B)
produces the terms retaining $(1+e^{\pm 2i\omega t})$ in Eq. 5. With the
_rotating wave approximation_ (RWA), the highly oscillatory terms can be
ignored in comparison with the driving frequency $\Omega$, which approximates
$(1+e^{\pm 2i\omega t})$ to 1. To be specific, the term $e^{\pm 2i\omega t}$
vanishes when it is averaged over a time interval $\sim\Omega^{-1}$ with the
condition of $\omega/\Omega\gg 1$:

$\displaystyle\Big{|}\frac{1}{\Omega^{-1}}\int_{-\frac{1}{2\Omega}}^{\frac{1}{2\Omega}}e^{\pm
2i\omega
t}dt\Big{|}=\Big{|}\frac{\sin{(\omega/\Omega)}}{\omega/\Omega}\Big{|}\leq\Omega/\omega\ll
1.$ (6)

Thus, the RWA simplifies Eq. 5 to

$\displaystyle\frac{\partial\tilde{\rho}_{ee}}{\partial\tau}=$
$\displaystyle-(\bar{n}+1)\tilde{\rho}_{ee}-i\Omega_{\gamma}\tilde{\rho}_{eg}+i\Omega_{\gamma}\tilde{\rho}_{ge}+\bar{n}\tilde{\rho}_{gg}$
$\displaystyle\frac{\partial\tilde{\rho}_{eg}}{\partial\tau}=$
$\displaystyle-i\Omega_{\gamma}\tilde{\rho}_{ee}+\left[i\delta\omega_{\gamma}-\frac{(2\bar{n}+1)}{2}\right]\tilde{\rho}_{eg}+i\Omega_{\gamma}\tilde{\rho}_{gg}$
$\displaystyle\frac{\partial\tilde{\rho}_{ge}}{\partial\tau}=$ $\displaystyle
i\Omega_{\gamma}\tilde{\rho}_{ee}-\left[i\delta\omega_{\gamma}+\frac{(2\bar{n}+1)}{2}\right]\tilde{\rho}_{ge}-i\Omega_{\gamma}\tilde{\rho}_{gg}$
$\displaystyle\frac{\partial\tilde{\rho}_{gg}}{\partial\tau}=$
$\displaystyle(\bar{n}+1)\tilde{\rho}_{ee}+i\Omega_{\gamma}\tilde{\rho}_{eg}-i\Omega_{\gamma}\tilde{\rho}_{ge}-\bar{n}\tilde{\rho}_{gg}$
(7)

where the whole equation is rescaled with $\gamma^{-1}$, with the definitions
of dimensionless parameters $\tau\equiv\gamma t$,
$\Omega_{\gamma}\equiv\Omega/\gamma$, and
$\delta\omega_{\gamma}\equiv\delta\omega/\gamma$. It is noteworthy that once
the RWA is taken under the condition of $\omega/\Omega\gg 1$, the dynamics of
our TLS model studied in the rotating frame with the angular frequency
$\omega$ depends only on $\delta\omega_{\gamma}$, and is impervious to the
frequency ($\omega$) of the driving field ($\vec{E}(t)$) (see Eq. 7).

The steady state values of density matrix elements are obtained from Eq. 7 by
setting $d\tilde{\rho}/d\tau=0$ with the constraints of
$\tilde{\rho}_{ee}^{ss}+\tilde{\rho}_{gg}^{ss}=1$ and
$\tilde{\rho}_{eg}^{ss}=(\tilde{\rho}_{ge}^{ss})^{\ast}$

$\displaystyle\tilde{\rho}_{ee}^{ss}=\dfrac{\bar{n}[(2\bar{n}+1)^{2}+4\delta\omega_{\gamma}^{2}]+4\Omega_{\gamma}^{2}(2\bar{n}+1)}{(2\bar{n}+1)[(2\bar{n}+1)^{2}+4\delta\omega_{\gamma}^{2}+8\Omega_{\gamma}^{2}]}$
$\displaystyle\tilde{\rho}_{gg}^{ss}=\dfrac{(\bar{n}+1)[(2\bar{n}+1)^{2}+4\delta\omega_{\gamma}^{2}]+4\Omega_{\gamma}^{2}(2\bar{n}+1)}{(2\bar{n}+1)[(2\bar{n}+1)^{2}+4\delta\omega_{\gamma}^{2}+8\Omega_{\gamma}^{2}]}$
$\displaystyle\tilde{\rho}_{eg}^{ss}=\frac{-2\Omega_{\gamma}\left[2\delta\omega_{\gamma}-i(2\bar{n}+1)\right]}{(2\bar{n}+1)[(2\bar{n}+1)^{2}+4\delta\omega_{\gamma}^{2}+8\Omega_{\gamma}^{2}]}.$
(8)

For $\Omega_{\gamma}=0$, the TLS coupled only to the thermal bath should be at
equilibrium, satisfying the DB condition [46], which can be confirmed from Eq.
8. The ratio of excited and ground state populations obeys the Boltzmann
distribution,
$\tilde{\rho}^{ss}_{ee}/\tilde{\rho}^{ss}_{gg}=\bar{n}/(\bar{n}+1)=e^{-\beta\hbar\omega_{0}}$,
and the quantum coherence vanishes ($\tilde{\rho}_{eg}^{ss}=0$). In contrast,
the model with $\Omega_{\gamma}\neq 0$ results in the breakdown of the DB
condition and nonvanishing coherence ($\tilde{\rho}_{eg}^{ss}\neq 0$).

We introduce simplified notations for the real
($\rho_{R}\equiv\text{Re}(\tilde{\rho}_{eg}^{ss})$) and imaginary parts
($\rho_{I}\equiv\text{Im}(\tilde{\rho}_{eg}^{ss})$) of the coherence
($\tilde{\rho}^{ss}_{eg}=\rho_{R}+i\rho_{I}$):

$\displaystyle\rho_{R}$
$\displaystyle\equiv\frac{-4\Omega_{\gamma}\delta\omega_{\gamma}}{\coth{\frac{\mathcal{A}}{2}}[\coth^{2}{\frac{\mathcal{A}}{2}}+4\delta\omega_{\gamma}^{2}+8\Omega_{\gamma}^{2}]}$
$\displaystyle\rho_{I}$
$\displaystyle\equiv\frac{2\Omega_{\gamma}}{[\coth^{2}{\frac{\mathcal{A}}{2}}+4\delta\omega_{\gamma}^{2}+8\Omega_{\gamma}^{2}]},$
(9)

where we have used $\coth{\frac{\mathcal{A}}{2}}=(2\bar{n}+1)$ with
$k_{B}T\mathcal{A}\equiv\hbar\omega_{0}$ which corresponds to the energy gap
between the ground and excited states. The gain of energy eventually
dissipates into the bath along a single cycle of absorption and emission. A
couple of remarks on $\rho_{R}$ and $\rho_{I}$ are in place. (i) Both
$\rho_{R}$ and $\rho_{I}$ vanish with the driving frequency as $\sim
1/\Omega_{\gamma}$. (ii) $\rho_{I}$ takes the form of Lorentzian with respect
to $\delta\omega_{\gamma}$, whereas $\rho_{R}$ is an odd function of
$\delta\omega_{\gamma}$. In fact, the quantum coherence can be related to the
response function (or susceptibility) of the system between the external
electric field ($\vec{E}$) and the polarization ($\vec{P}$) in the linear
response regime ($\vec{P}=\chi\vec{E}$). More specifically, $\rho_{R}$ and
$\rho_{I}$ are related to the real and imaginary parts of the susceptibility,
which inform about the dispersion and absorption profile of the light-matter
interaction, respectively (see Appendix C).

##  III Mean current and current fluctuations

The mean current and the current fluctuations involving the net number of
transitions, $n(\tau)$, i.e., the difference between the total numbers of
emissions and absorptions for time $\tau$, at steady states, elicited by the
irradiation are calculated by employing the method of generating function [47,
48] (see Appendix D):

$\displaystyle\langle j\rangle$
$\displaystyle\equiv\lim_{\tau\rightarrow\infty}\frac{\langle
n(\tau)\rangle}{\tau}=\frac{4\Omega_{\gamma}^{2}}{\left[\coth^{2}{\frac{\mathcal{A}}{2}}+4\delta\omega_{\gamma}^{2}+8\Omega_{\gamma}^{2}\right]}$
$\displaystyle=2\Omega_{\gamma}\rho_{I},$ (10)

and

$\displaystyle\text{Var}[j]$
$\displaystyle\equiv\lim_{\tau\rightarrow\infty}\frac{\text{Var}[n(\tau)]}{\tau}$
$\displaystyle=\langle
j\rangle\coth{\left(\frac{\mathcal{A}}{2}\right)}f(\mathcal{A},\Omega_{\gamma},\delta\omega_{\gamma})$
(11)

with

$\displaystyle f(\mathcal{A},\Omega_{\gamma},\delta\omega_{\gamma})$
$\displaystyle\equiv\left(1+\dfrac{[32\delta\omega_{\gamma}^{2}-24\coth^{2}{\frac{\mathcal{A}}{2}}]\Omega_{\gamma}^{2}}{\coth^{2}{\frac{\mathcal{A}}{2}}\left[\coth^{2}{\frac{\mathcal{A}}{2}}+4\delta\omega_{\gamma}^{2}+8\Omega_{\gamma}^{2}\right]^{2}}\right)$
$\displaystyle=\left[1+2\rho_{R}^{2}-6\rho_{I}^{2}\right].$ (12)

The expression of mean current (Eq. 10) conforms to the _coherence-current
relation_ [49, 50] which states that the imaginary part of the coherence (Eq.
9) is responsible for inducing a quantum current between the states. Notice
that $\langle j\rangle\geq 0$ for the entire parameter space with $\langle
j\rangle=0$ realized when the driving field is absent ($\Omega_{\gamma}=0$) or
the energy gap (dissipation to the bath) is zero ($\mathcal{A}=0$) and that
$\langle j\rangle$ increases monotonically with $\Omega_{\gamma}$ and upper-
bounded to $\langle j\rangle\rightarrow 1/2$ for $\Omega_{\gamma}\gg 1$
($y\neq 0$). Next, Eq. 11 clarifies that $\text{Var}[j]=0$ when
$\Omega_{\gamma}=0$. For $\Omega_{\gamma}\neq 0$ and
$|\delta\omega_{\gamma}|\gg 1$, the current fluctuations are upper bounded as
$\text{Var}[j]\leq(1/2)\coth{(\frac{\mathcal{A}}{2})}$; and for
$\delta\omega_{\gamma}=0$, it is simplified to $\text{Var}[j]/\langle
j\rangle=\coth{(\frac{\mathcal{A}}{2})}(1-6\rho_{I}^{2})$, where
$\rho_{I}=\rho_{I}(\Omega_{\gamma})$ (Eq.9) is a non-monotonic function of
$\Omega_{\gamma}$.

##  IV Steady-state entropy production rate

The presence of irreversible current is conducive to the entropy production.
The total entropy production rate in the open quantum system is given by [42,
51]

$\displaystyle\Sigma$ $\displaystyle=\dfrac{dS}{d\tau}+J_{S}.$ (13)

The entropy of the system is defined with the von Neumann entropy,
$S(\tilde{\rho})=-k_{B}\text{Tr}(\tilde{\rho}\log{\tilde{\rho}})$, and $J_{S}$
corresponds to the heat current into the bath. At steady state
($\dot{\tilde{\rho}}=0$), $dS(\tilde{\rho})/d\tau=0$, so that $\Sigma$ is
contributed by $J_{S}$ part only. In the Lindblad framework, $J_{S}$ is
expressed as [42, 51]
$J_{S}=(\hbar\omega_{0}/T)\text{Tr}\left[(\bar{n}+1)\sigma_{+}\sigma_{-}\tilde{\rho}-\bar{n}\sigma_{-}\sigma_{+}\tilde{\rho}\right]$.
Therefore, the total entropy production rate at steady states is

$\displaystyle\Sigma^{ss}/k_{B}$
$\displaystyle=\beta\hbar\omega_{0}\text{Tr}\left[(\bar{n}+1)\sigma_{+}\sigma_{-}\tilde{\rho}^{ss}-\bar{n}\sigma_{-}\sigma_{+}\tilde{\rho}^{ss}\right]$
$\displaystyle=\beta\hbar\omega_{0}\left[(\bar{n}+1)\tilde{\rho}^{ss}_{ee}-\bar{n}\tilde{\rho}^{ss}_{gg}\right]$
$\displaystyle=\beta\hbar\omega_{0}\dfrac{4\Omega_{\gamma}^{2}}{\left[(2\bar{n}+1)^{2}+4\delta\omega_{\gamma}^{2}+8\Omega_{\gamma}^{2}\right]}$
$\displaystyle=\mathcal{A}\times\langle j\rangle\geq 0.$ (14)

Despite $dS(\tilde{\rho})/d\tau=0$ at NESS, the non-vanishing heat current due
to the interaction with the external field ($J_{S}^{ss}\neq 0$) gives rise to
$\Sigma^{ss}\neq 0$. We also note that the value of $\Sigma^{ss}$ is
unaffected by the transformation into the rotating frame
$\tilde{\rho}_{ij}=e^{i\omega\sigma_{z}t/2}\rho_{ij}e^{-i\omega\sigma_{z}t/2}$.

Figure 2: The uncertainty product $\mathcal{Q}$. (A) $\mathcal{Q}$ as a
function of $\Omega_{\gamma}$ and $\mathcal{A}$ for the case of perfect
resonance ($\delta\omega_{\gamma}=0$). The condition of
$(\Omega_{\gamma},\mathcal{A})$ that minimizes $\mathcal{Q}$ is marked with
the star symbol. (B) The diagram of
$\mathcal{Q}=\mathcal{Q}(\Omega_{\gamma},\delta\omega_{\gamma})$ for
$\mathcal{A}=3.6$ and $\mathcal{A}=8$. (C) Dissections of $\mathcal{Q}$.
(left) $\mathcal{Q}(\Omega_{\gamma}=0.37,\delta\omega_{\gamma})$ at
$\mathcal{A}=3.8$ and 8. (right)
$\mathcal{Q}(\Omega_{\gamma},\delta\omega_{\gamma}=0)$ for various values of
$\mathcal{A}$. The bound of the original TUR ($\mathcal{Q}=2$) is marked with
the blue dashed line.

##  V Thermodynamic uncertainty relation

The uncertainty product of TUR at NESS [1, 4] for the net transitions of the
TLS is written as

$\displaystyle\mathcal{Q}$
$\displaystyle=\lim_{\tau\rightarrow\infty}\dfrac{\Sigma^{ss}\tau}{k_{B}}\dfrac{\langle\delta
n(\tau)^{2}\rangle}{\langle
n(\tau)\rangle^{2}}=\underbrace{\left(\dfrac{\Sigma^{ss}}{k_{B}\langle
j\rangle}\right)}_{=\mathcal{A}}\underbrace{\left(\dfrac{\text{Var}[j]}{\langle
j\rangle}\right)}_{=\mathcal{F}}$
$\displaystyle=\mathcal{A}\coth{\left(\frac{\mathcal{A}}{2}\right)}f(\mathcal{A},\Omega_{\gamma},\delta\omega_{\gamma}).$
(15)

The original TUR derived for classical system states that the uncertainty
product $\mathcal{Q}$ between the cost ($\Sigma^{ss}/k_{B}$) and the precision
dictated by the square of the relative fluctuations of current
($\text{Var}[j]/\langle j\rangle^{2}$) cannot be smaller than $2$ [1, 2, 4].
However, when TUR is extended to quantum systems, the uncertainty product can
have a loose bound. The expression given in Eq. 12 points to the possibility
of $f(\mathcal{A},\Omega_{\gamma},\delta\omega_{\gamma})<1$, which in turn
leads to $\mathcal{Q}<2$, when the imaginary part of coherence is dominant. In
fact, $\rho_{R}=0$ for the case of perfect resonance
($\delta\omega_{\gamma}=0$) simplifies $\mathcal{Q}$ to
$\mathcal{Q}(\mathcal{A},\Omega_{\gamma})=\mathcal{A}\coth{(\frac{\mathcal{A}}{2})}\left[1-\frac{24\Omega_{\gamma}^{2}}{(\coth^{2}{(\frac{\mathcal{A}}{2})}+8\Omega_{\gamma}^{2})^{2}}\right]\geq\mathcal{Q}(\mathcal{A}^{\ast},\Omega_{\gamma}^{\ast})\approx
1.25$ at $(\mathcal{A}^{\ast},\Omega_{\gamma}^{\ast})\approx(3.61,0.37)$ (see
Fig. 2A). The minimal uncertainty product, $\mathcal{Q}\geq\mathcal{Q}_{\rm
min}=1.25$ coincides with the value obtained in a recent work that has studied
TUR of a spin in a rotating magnetic field [52].

For the case of finite detuning (off-resonance) $\delta\omega_{\gamma}\neq 0$,
the real part of coherence ($\rho_{R}\neq 0$) contributes to increasing the
size of $\mathcal{Q}$ (Fig. 2B, C). Thus, for a given set of parameter
$(\mathcal{A},\Omega_{\gamma},\delta\omega_{\gamma})$, the value of
$\mathcal{Q}$ is determined as a result of competition between the real and
imaginary parts of the coherence. Notably, $\rho_{R}^{2}<3\rho_{I}^{2}$, i.e.,
$|\delta\omega_{\gamma}|<(\sqrt{3}/2)\coth{\frac{\mathcal{A}}{2}}$, is a
necessary, but not a sufficient condition for the loose bound of TUR.

TUR can also be recast into the form of an inequality of the Fano factor of
transitions [30]. For our system, Eqs. 10 and 11 define the Fano factor of
transitions,

$\displaystyle\mathcal{F}$
$\displaystyle=\lim_{\tau\rightarrow\infty}\frac{\text{Var}[n(\tau)]}{\langle
n(\tau)\rangle}=\frac{\text{Var}[j]}{\langle j\rangle}$
$\displaystyle=\coth{(\frac{\mathcal{A}}{2})}f(\mathcal{A},\Omega_{\gamma},\delta\omega_{\gamma})\geq\phi(\mathcal{A}).$
(16)

For the original TUR derived for the case of constant driving [1, 2, 4],
$\phi(\mathcal{A})=\phi_{o}(\mathcal{A})=2/\mathcal{A}$. For the case of
_time-symmetric periodic driving_ , Proesmans and Van den Broeck have shown
that the steady state current fluctuations per period in the long-time limit
is bounded by $\text{Var}[j]/\langle j\rangle^{2}\geq 2\Delta
t/(e^{\Sigma^{ss}\Delta t/k_{B}}-1)$ [8, 4]. Our electromagnetic field-driven
quantum TLS is, in principle, under time-symmetric periodic driving with a
constant period of $\Delta t=2\pi/\omega$; however, such a periodicity in
driving field under weak field limit ($\Omega\ll\omega$) is reduced to a
constant driving under RWA, as detailed in Eqs. 6 and 7. Thus, it is pertinent
to ask whether or not the uncertainty product of our quantum TLS satisfies the
inequality dictated by the classical version of TUR under a constant driving.
Since $\langle j\rangle\Delta t=1$ from the definition given in Eq. 10 and
$\Sigma^{ss}\Delta t=\mathcal{A}$, the expression can be cast into
$\text{Var}[j]/\langle
j\rangle\geq\phi(\mathcal{A})=\phi_{p}(\mathcal{A})=2/(e^{\mathcal{A}}-1)$.
Mathematically, $\phi_{o}(\mathcal{A})>\phi_{p}(\mathcal{A})$ for all
$\mathcal{A}\geq 0$, but $\phi_{p}(\mathcal{A})$ converges to the form of
$\phi_{o}(\mathcal{A})$ when $\mathcal{A}\ll 1$.

$\mathcal{F}$, plotted against $\mathcal{A}$ with randomly varying set of
parameters $(\mathcal{A},\Omega_{\gamma},\delta\omega_{\gamma})$, is upper
bounded by $(5/4)\coth{\frac{\mathcal{A}}{2}}$ (the cyan line in Fig. 3),
which is obtained from
$f(\mathcal{A},\Omega_{\gamma},\delta\omega_{\gamma})\approx
1+\frac{32\delta\omega_{\gamma}^{2}\Omega_{\gamma}^{2}}{[4\delta\omega_{\gamma}^{2}+8\Omega_{\gamma}^{2}]^{2}}\leq
5/4$ for $\mathcal{A}\gg 1$, $|\delta\omega_{\gamma}|\gg 1$ and
$\Omega_{\gamma}\gg 1$. Of particular note is that there are data points that
satisfy
$\phi_{p}(\mathcal{A})<\mathcal{F}(\mathcal{A},\Omega_{\gamma},\delta\omega_{\gamma})<\phi_{o}(\mathcal{A})$
at $1\lesssim\mathcal{A}\lesssim 6$, signifying that the TUR of our quantum
TLS is characterized with the bound less tight than that of the original TUR
(Fig. 3).

Figure 3: Fano factor ($\mathcal{F}$) of a bath-coupled TLS under a weak
driving plotted against the net entropy production over a cycle of absorption
and emission, $\mathcal{A}(=\beta\hbar\omega_{0})$. For $\mathcal{F}$
calculated with a randomly generated set of parameters
$(\mathcal{A},\Omega_{\gamma},\delta\omega_{\gamma})$, there are cases whose
Fano factor is smaller than the lower bound of the original TUR (blue), but
still greater than the bound specified by the generalized uncertainty relation
(red) [8].

##  VI Discussions

The relation between the parameters $\omega$, $\omega_{0}$, $\gamma$, and
$\Omega$,

$\displaystyle\omega\sim\omega_{0}\gg\gamma\sim\Omega$ (17)

is the essential condition to study the time-irreversible evolution of our
TLS. First, for transitions induced by dipole-electric field coupling, Fermi’s
golden rule offers an expression for the emission rate,
$\gamma=4\omega_{0}^{3}|d_{eg}|^{2}/(3\hbar c^{3})$ [42]. Provided that the
dipole moment of a molecule (or an atom) is given by $|d_{eg}|=ea_{0}$ with
its electric charge $e$ and size $a_{0}(\lesssim 1$ nm), and that
$\omega_{0}/c=2\pi/\lambda$ where $c$ is the speed of light, and
$\lambda(\approx 500$ nm) is the wavelength of the field in the optical
domain, one gets $\omega_{0}/\gamma\approx 6.5\times 10^{5}$ from
$\gamma=(4\alpha/3)\left(2\pi a_{0}/\lambda\right)^{2}\omega_{0}$, with
$\alpha\equiv e^{2}/\hbar c=1/137$ and $2\pi a_{0}/\lambda\sim 10^{-2}$. The
separation of timescales ($\omega/\gamma\gg 1$) justifies decomposition of the
total density matrix into the system and equilibrated bath parts, and
rationalizes the truncation of quantum optical master equation at the second
order (_Born-Markov approximation_. See Appendix A). Second, the condition of
$\omega\gg\Omega$ was used in taking the RWA in Eq. 6 to obtain Eq. 7.

Since the Eq. 7 is independent of $\omega$, the system is effectively under a
constant driving. As a result, the system at long time limit, if classical,
should obey the original TUR [1, 31], instead of the generalized TUR under
time-symmetric periodic driving [8]. Thus, the loose TUR bound of our TLS
model ($\mathcal{Q}<2$ in Fig.2) stems from the quantum nature of the
dynamics. Our model further elaborates the contribution from the quantum
nature to TUR: The imaginary part of quantum coherence contributes to
suppressing the Fano factor of net transitions and lowering the TUR bound,
whereas the real part of the coherence always tightens the TUR bound. The
bound of TUR is determined by the competition between the dissipation and
dispersion of the input light source.

For the cases the system-bath coupling ($\gamma$) and/or the strength of
driving field ($\Omega$) are strong and become comparable to the energy scale
of the system, one should reevaluate the quantum master equation presented
here based on generalized versions of open quantum dynamics that are under
active development in recent years [53, 54, 55]. However, as long as the
energy gap of the TLS is in the optical domain ($\lambda\sim 500$ nm), one of
the conditions identified for the loose TUR bound in Fig. 2
($\Omega_{\gamma}\lesssim 0.5$) ensures the relation specified in Eq. 17,
which should reduce a generalized quantum master equation into the one studied
here. Thus, we anticipate that the condition leading to the loose TUR bound
will remain unaltered. Finally, the other condition $\mathcal{A}\sim(2-8)$
(Fig. 2) becomes accessible in the range of high temperature
$T=(2-8)^{-1}(hc/k_{B}\lambda)\sim\mathcal{O}(10^{3})$ K for $\lambda\sim 500$
nm.

Despite the apparent simplicity of the model presented here, it retains all
the essential components that enables us to explore the TUR of quantum
processes. The model also offers us a physically straightforward understanding
of the condition that leads to loosening the TUR bound.

###### Acknowledgements.

This work was supported by the KIAS Individual Grants CG077602 (D.S.) and
CG035003 (C.H.) from the Korea Institute for Advanced Study. We thank the
Center for Advanced Computation in KIAS for providing computing resources.

## Appendix

###  VI.1 Derivation of Lindblad dissipator

Introducing the transformations

$\displaystyle H^{\prime}_{I}(t)$
$\displaystyle=e^{\frac{i}{\hbar}H_{o}t}H_{I}(t)e^{-\frac{i}{\hbar}H_{o}t}$
$\displaystyle\rho^{\prime}_{\text{tot}}(t)$
$\displaystyle=e^{\frac{i}{\hbar}H_{o}t}\rho_{\text{tot}}(t)e^{-\frac{i}{\hbar}H_{o}t}$
(A1)

with the definitions of

$\displaystyle H_{o}$ $\displaystyle\equiv H_{S}+H_{B}$ $\displaystyle
H_{I}(t)$ $\displaystyle\equiv H_{\text{ext}}(t)+H_{SB},$ (A2)

one can cast the von Neumann equation into the one in the interaction picture
[43]

$\displaystyle\frac{d\rho^{\prime}_{\text{tot}}(t)}{dt}=-\frac{i}{\hbar}[H^{\prime}_{I}(t),\rho^{\prime}_{\text{tot}}(t)].$
(A3)

After integrating Eq.(A3) and expanding it up to the second order, we obtain

$\displaystyle\frac{d{\rho^{\prime}_{\text{tot}}}(t)}{dt}$
$\displaystyle=-\frac{i}{\hbar}[H^{\prime}_{SB}(t),\rho^{\prime}_{\text{tot}}(0)]-\frac{i}{\hbar}[H^{\prime}_{\text{ext}}(t),\rho^{\prime}_{\text{tot}}(0)]$
$\displaystyle-\frac{1}{\hbar^{2}}\int_{0}^{t}dt^{\prime}[H^{\prime}_{I}(t),[H^{\prime}_{I}(t^{\prime}),\rho^{\prime}_{\text{tot}}(t^{\prime})]].$
(A4)

Next, we trace out the bath modes, assuming that the density matrix
$\rho_{\text{tot}}(t)$ for the total system at time $t$ can be described with
the uncorrelated product state of system ($\rho(t)$) and environment in
equilibrium ($\rho_{B}^{eq}$), such that
$\rho_{\text{tot}}(t)\approx\rho(t)\otimes\rho_{B}^{eq}$. Then, it follows
that

$\displaystyle\frac{d{\rho^{\prime}}(t)}{dt}=-\frac{i}{\hbar}[H^{\prime}_{\text{ext}}(t),\rho^{\prime}(0)]$
$\displaystyle-\frac{1}{\hbar^{2}}\int_{0}^{t}dt^{\prime}\Tr_{B}\left([H^{\prime}_{\text{ext}}(t),[H^{\prime}_{\text{ext}}(t^{\prime}),\rho^{\prime}(t^{\prime})\otimes\rho_{B}^{eq}]]\right)$
$\displaystyle-\frac{1}{\hbar^{2}}\int_{0}^{t}dt^{\prime}\Tr_{B}\left([H^{\prime}_{SB}(t),[H^{\prime}_{SB}(t^{\prime}),\rho^{\prime}(t^{\prime})\otimes\rho_{B}^{eq}]]\right),$
(A5)

where $\rho^{\prime}=\Tr_{B}\left(\rho^{\prime}_{\text{tot}}\right)$ was used
in the left hand side of the equation and an assumption that the system is
decoupled from the environment at $t=0$ leads to
$\Tr_{B}[H^{\prime}_{SB}(t),\rho^{\prime}_{\text{tot}}(0)]=0$.

Next, to proceed further, it is convenient to define
$H_{SB}=\hbar\sum\limits_{m}s_{m}q_{m}$, where $s_{m}$ are the system
operators and $q_{m}$ are the bath operators, with $s_{1}\equiv\sigma_{-}$,
$s_{2}\equiv\sigma_{+}$,
$q_{1}\equiv\sum\limits_{\textbf{k},\xi}g_{\textbf{k},\xi}^{\ast}b_{\textbf{k},\xi}^{\dagger}$,
and
$q_{2}\equiv\sum\limits_{\textbf{k},\xi}g_{\textbf{k},\xi}b_{\textbf{k},\xi}$.
Then

$\displaystyle H^{\prime}_{SB}(t)$
$\displaystyle=\hbar\sum\limits_{m=1,2}e^{\frac{i}{\hbar}(H_{S}+H_{B})t}s_{m}q_{m}e^{-\frac{i}{\hbar}(H_{S}+H_{B})t}$
$\displaystyle=\hbar\sum\limits_{m=1,2}\left(e^{\frac{i}{\hbar}H_{S}t}s_{m}e^{-\frac{i}{\hbar}H_{S}t}\right)\left(e^{\frac{i}{\hbar}H_{B}t}q_{m}e^{-\frac{i}{\hbar}H_{B}t}\right)$
$\displaystyle=\hbar\sum\limits_{m=1,2}s_{m}(t)q_{m}(t)$ (A6)

where one can show using Baker-Campbell-Hausdorff formula that
$s_{1}(t)=\sigma_{-}e^{-i\omega_{0}t}$, $s_{2}(t)=\sigma_{+}e^{i\omega_{0}t}$,
$q_{1}(t)=\sum\limits_{{\bf k},\xi}g^{\ast}_{{\bf k},\xi}b^{\dagger}_{{\bf
k},\xi}e^{i\omega_{\bf k}t}$, and $q_{2}(t)=\sum\limits_{{\bf k},\xi}g_{{\bf
k},\xi}b_{{\bf k},\xi}e^{-i\omega_{\bf k}t}$. The master equation can be
written as

$\displaystyle\frac{d\rho^{\prime}(t)}{dt}=-\frac{i}{\hbar}[H^{\prime}_{\text{ext}}(t),\rho^{\prime}(0)]$
$\displaystyle-\frac{1}{\hbar^{2}}\int_{0}^{t}dt^{\prime}\Tr_{B}\left([H^{\prime}_{\text{ext}}(t),[H^{\prime}_{\text{ext}}(t^{\prime}),\rho^{\prime}(t^{\prime})\otimes\rho_{B}^{eq}]]\right)$
$\displaystyle-\frac{1}{\hbar^{2}}\int_{0}^{t}dt^{\prime}\left[\left(s_{1}(t)s_{2}(t^{\prime})\rho^{\prime}(t^{\prime})-s_{2}(t^{\prime})\rho^{\prime}(t^{\prime})s_{1}(t)\right)\langle
q_{1}(t)q_{2}(t^{\prime})\rangle_{B}\right.$
$\displaystyle\left.\qquad\qquad+\left(s_{2}(t)s_{1}(t^{\prime})\rho^{\prime}(t^{\prime})-s_{1}(t^{\prime})\rho^{\prime}(t^{\prime})s_{2}(t)\right)\langle
q_{2}(t)q_{1}(t^{\prime})\rangle_{B}\right.$
$\displaystyle\left.\qquad\qquad+\left(\rho^{\prime}(t^{\prime})s_{1}(t^{\prime})s_{2}(t)-s_{2}(t)\rho^{\prime}(t^{\prime})s_{1}(t^{\prime})\right)\langle
q_{1}(t^{\prime})q_{2}(t)\rangle_{B}\right.$
$\displaystyle\left.\qquad\qquad+\left(\rho^{\prime}(t^{\prime})s_{2}(t^{\prime})s_{1}(t)-s_{1}(t)\rho^{\prime}(t^{\prime})s_{2}(t^{\prime})\right)\langle
q_{2}(t^{\prime})q_{1}(t)\rangle_{B}\right],$ (A7)

where $\langle
q_{m}(t)q_{n}(t^{\prime})\rangle_{B}=\text{Tr}_{B}\left[\rho_{B}^{eq}q_{m}(t)q_{n}(t^{\prime})\right]$
represents the bath correlation function. For example,

$\displaystyle\langle q_{1}(t)q_{2}(t^{\prime})\rangle_{B}$
$\displaystyle=\sum_{{\bf k},\xi}\sum_{{\bf k}^{\prime},\xi^{\prime}}g_{{\bf
k},\xi}^{\ast}g_{{\bf k}^{\prime},\xi^{\prime}}e^{i\omega_{\bf
k}t}e^{-i\omega_{{\bf
k}^{\prime}}t^{\prime}}\text{Tr}_{B}(\rho_{B}^{eq}b_{{\bf
k},\xi}^{\dagger}b_{{\bf k}^{\prime},\xi^{\prime}})$ $\displaystyle=\sum_{{\bf
k}}|g_{{\bf k},\xi}|^{2}e^{i\omega_{\bf k}(t-t^{\prime})}\bar{n}(\omega_{\bf
k}).$ (A8)

Likewise

$\displaystyle\langle q_{2}(t)q_{1}(t^{\prime})\rangle_{B}=\sum_{{\bf
k}}|g_{{\bf k}}|^{2}e^{-i\omega_{\bf k}(t-t^{\prime})}(\bar{n}(\omega_{\bf
k})+1).$ (A9)

By substituting $s_{1,2}$, $q_{1,2}$ back, and transforming the time
$u=t-t^{\prime}$, we obtain

$\displaystyle\frac{d\rho^{\prime}(t)}{dt}=-\frac{i}{\hbar}[H^{\prime}_{\text{ext}}(t),\rho^{\prime}(0)]$
$\displaystyle-\frac{1}{\hbar^{2}}\int_{0}^{t}dt^{\prime}\Tr_{B}\left([H^{\prime}_{\text{ext}}(t),[H^{\prime}_{\text{ext}}(t^{\prime}),\rho^{\prime}(t^{\prime})\otimes\rho_{B}^{eq}]]\right)$
$\displaystyle-\sum_{{\bf k},\xi}|g_{{\bf
k},\xi}|^{2}\int_{0}^{t}du\left[\left(\sigma_{-}\sigma_{+}\rho^{\prime}(t-u)-\sigma_{+}\rho^{\prime}(t-u)\sigma_{-}\right)e^{i(\omega_{\bf
k}-\omega_{0})u}\bar{n}(\omega_{\bf k})\right.$
$\displaystyle\left.\qquad\qquad+\left(\sigma_{+}\sigma_{-}\rho^{\prime}(t-u)-\sigma_{-}\rho^{\prime}(t-u)\sigma_{+}\right)e^{-i(\omega_{\bf
k}-\omega_{0})u}(\bar{n}(\omega_{\bf k})+1)\right.$
$\displaystyle\left.\qquad\qquad+\left(\rho^{\prime}(t-u)\sigma_{-}\sigma_{+}-\sigma_{+}\rho^{\prime}(t-u)\sigma_{-}\right)e^{-i(\omega_{\bf
k}-\omega_{0})u}\bar{n}(\omega_{\bf k})\right.$
$\displaystyle\left.\qquad\qquad+\left(\rho^{\prime}(t-u)\sigma_{+}\sigma_{-}-\sigma_{-}\rho^{\prime}(t-u)\sigma_{+}\right)e^{i(\omega_{\bf
k}-\omega_{0})u}(\bar{n}(\omega_{\bf k})+1)\right].$ (A10)

To proceed further, _Markov approximation_ is often considered with
$\rho^{\prime}(t-u)\rightarrow\rho^{\prime}(t)$ and
$\int_{0}^{t}du\rightarrow\int_{0}^{\infty}du$ by assuming that the value of
density matrix is determined without memory, and
$\lim_{t\rightarrow\infty}\int_{-t}^{t}du\sum_{{\bf k},\xi}|g_{{\bf
k},\xi}|^{2}e^{\pm i(\omega_{\bf k}-\omega_{0})u}\rightarrow\gamma\equiv
2\pi\int d^{3}kJ({\bf k})|g_{{\bf k},\xi}|^{2}\delta(\omega_{\bf
k}-\omega_{0})$ where $J({\bf k})$ is the spectral density, determining the
decay rate $\gamma$ [43].

Transformed back to the Schrödinger picture, the master equation takes the
following form [43]

$\displaystyle\frac{d\rho(t)}{dt}=-\frac{i}{\hbar}[H_{S}+H_{\text{ext}}(t),\rho(t)]+\mathcal{D}(\rho(t))$
(A11)

with the Lindblad dissipator

$\displaystyle\mathcal{D}(\rho)$
$\displaystyle=\dfrac{\gamma}{2}(\bar{n}+1)\left(2\sigma_{-}\rho\sigma_{+}-\sigma_{+}\sigma_{-}\rho-\rho\sigma_{+}\sigma_{-}\right)$
$\displaystyle+\dfrac{\gamma}{2}\bar{n}\left(2\sigma_{+}\rho\sigma_{-}-\sigma_{-}\sigma_{+}\rho-\rho\sigma_{-}\sigma_{+}\right),$
(A12)

where $\bar{n}=\bar{n}(\omega_{0})=(e^{\beta\hbar\omega_{0}}-1)^{-1}$ is the
average number of thermal photons.

###  VI.2 Transformation to the rotating frame

The density matrix $\tilde{\rho}=\ket{\phi}\bra{\phi}$ in the rotating frame
is transformed into the one in the stationary frame
$\rho=\ket{\psi}\bra{\psi}$ via the following operation,

$\underbrace{\ket{\psi}\bra{\psi}}_{=\rho}=e^{-i\omega
t\sigma_{z}/2}\underbrace{\ket{\phi}\bra{\phi}}_{=\tilde{\rho}}e^{i\omega
t\sigma_{z}/2},$ (B1)

which allows us to express the excited state population in the stationary
frame into the one in the rotating frame as

$\tilde{\rho}_{ee}=e^{i\omega
t\sigma_{z}/2}\underbrace{\ket{e}\bra{e}}_{=\rho_{ee}}e^{-i\omega
t\sigma_{z}/2}.$ (B2)

Employing the Baker-Campbell-Hausdorff formula,
$e^{s\hat{A}}\hat{B}e^{-s\hat{A}}=\hat{B}+(s/1!)[\hat{A},\hat{B}]+(s^{2}/2!)[\hat{A},[\hat{A},\hat{B}]]+\cdots$,
one can show that $\tilde{\rho}_{ee}=\rho_{ee}$,
$\tilde{\rho}_{gg}=\rho_{gg}$, $\tilde{\rho}_{eg}=\rho_{eg}e^{i\omega t}$, and
$\tilde{\rho}_{ge}=\rho_{ge}e^{-i\omega t}$.

Figure 4: The susceptibility of the system as a function of scaled detuning
($\delta\omega_{\gamma}$) calculated at $\Omega_{\gamma}=0.5$ and
$\mathcal{A}=\beta\hbar\omega_{0}=2$. (A) the absorption profile (i.e.,
$\chi^{\prime\prime}$) of the TLS. (B) the dispersion profile (i.e.,
$\chi^{\prime}$) of the medium.

###  VI.3 Relation between the coherence and dielectric susceptibility

Inside a dielectric medium, the electric polarization, defined as the total
dipole moment per volume, i.e., $\vec{P}=N\langle\vec{d}\rangle$ with $N$ the
number of molecules per volume, is induced by the incident radiation field,
i.e., $\vec{P}=\chi\vec{E}$, where $\chi$ is the linear susceptibility of the
medium [45]. With
$\langle\vec{d}\rangle=\text{Tr}(\tilde{\rho}\vec{d})=\tilde{\rho}_{eg}\vec{d}_{ge}+\tilde{\rho}_{ge}\vec{d}_{eg}=e^{i\omega
t}(\rho_{eg}\vec{d}_{ge}+\rho_{ge}\vec{d}_{eg}e^{-2i\omega
t})\approx\tilde{\rho}_{eg}\vec{d}_{ge}$, the susceptibility can be expressed
as $\chi=N_{d}\tilde{\rho}_{eg}$ with $N_{d}=N|d_{ge}|/|\vec{\varepsilon}|$.
In the light-matter interactions, the properties of the medium can be
investigated by studying $\chi$, which is decomposed into the real
($\chi^{\prime}$) and imaginary part ($\chi^{\prime\prime}$), i.e.,
$\chi=\chi^{\prime}+i\chi^{\prime\prime}$. For the medium with $|\chi|\ll 1$,
since the refractive index ($\underline{n}$) is related with the dielectric
constant as $\underline{n}^{2}=\epsilon=1+4\pi\chi$ (in Gaussian units),

$\displaystyle\underline{n}\approx 1+2\pi\chi^{\prime}.$ (C1)

The incident radiation field is absorbed to the medium is quantified by the
attenuation coefficient $\alpha$, which is related to the wave vector as
$k=\beta+i\alpha/2$, which leads to $\vec{E}\sim e^{ikz}\sim e^{i\beta
z-\alpha z/2}$ and $|\vec{E}|^{2}\sim e^{-\alpha z}$. From
$k=(\omega/c)\underline{n}=(\omega/c)\sqrt{1+4\pi\chi^{\prime}+4\pi
i\chi^{\prime\prime}}\approx(\omega/c)(1+2\pi\chi^{\prime}+2\pi
i\chi^{\prime\prime})$, it follows that

$\displaystyle\alpha=4\pi\left(\frac{\omega}{c}\right)\chi^{\prime\prime}.$
(C2)

This implies that the real part of the susceptibility (i.e.,
$\chi^{\prime}=N_{d}\tilde{\rho}_{R}$) informs about the dispersion profile of
the medium, while the imaginary part of the susceptibility (i.e.,
$\chi^{\prime\prime}=N_{d}\tilde{\rho}_{I}$) describes the absorption profile
of the medium.

The degree of the absorption and dispersion varies with detuning of the field
(Fig. 4). For the perfect resonance (i.e., when $\delta\omega=0$) the
excitation is completely absorbed by the medium ($\chi^{\prime\prime}$
maximized in Fig. 4) without any dispersion ($\chi^{\prime}=0$ in Fig. 4).
With increasing $|\delta\omega|$, the absorption (dispersion) in the medium
decreases (increases).

###  VI.4 Method of generating function

The dynamical equations (Eq. 7) can be cast into the matrix form,

$\displaystyle\partial_{\tau}\tilde{\rho}(\tau)=\mathcal{L}\tilde{\rho}(\tau).$
(D1)

where
$\tilde{\rho}=(\tilde{\rho}_{ee},\tilde{\rho}_{eg},\tilde{\rho}_{ge},\tilde{\rho}_{gg})^{T}$
is a vector in Liouville space, and $\mathcal{L}$ denotes the $4\times 4$
Liouvillian super-operator. We aim to count the net number of photon transfer
to the surrounding bath for time interval $\tau$ ($n(\tau)$) and calculate its
mean current and current fluctuations at steady state. Since the light-induced
absorption and emission is a cyclic process for the system, we assume that the
system is periodic in state space and that the superoperator $\mathcal{L}$ is
periodic with its period by $L(=4)$, satisfying
$\mathcal{L}_{ij}=\mathcal{L}_{i+L,j+L}$ for all $i$, $j=0$, $1$, $2$, $3$
representing the quantum state $ee$, $eg$, $ge$, $gg$, respectively, and
introduce a generalized coordinate $\mu\in\mathbb{Z}$, which is related with
the numeric index $i=0$, $1$, $2$, $3$ as $\mu=i\text{ }(\text{mod}\text{
}L)\text{ with }L=4$. This allows us to count the cycles of absorptions
followed by emissions via
$\tilde{\rho}_{i}(\mu,\tau)\equiv\tilde{\rho}(\mu,\tau)\delta_{\mu,i}^{L}$
with the generalized Kronecker delta, $\delta_{\mu,i}^{L}=1$ if $\mu=i$ (mod
$L$) and $\delta_{\mu,i}^{L}=0$ otherwise [47, 48]. Here,
$\tilde{\rho}_{i}(\mu,\tau)$ denotes the probability of finding a system at a
given quantum state $i\in\\{ee,eg,ge,gg\\}$ at site $\mu$ at time $\tau$.

We consider a generating function for each element $\tilde{\rho}_{i}$ of the
vector $\tilde{\rho}$,

$\displaystyle\mathcal{G}_{i}(z,\tau)=\sum_{\mu=-\infty}^{\infty}e^{zn_{\mu}}\tilde{\rho}_{i}(\mu,\tau).$
(D2)

Multiplying the factor $e^{zn_{\mu}}$ to both sides of a generalized version
of Eq. D1, namely
$\partial_{\tau}\tilde{\rho}_{i}(\mu,\tau)=\sum_{j}\mathcal{L}_{ij}\tilde{\rho}_{j}(\mu,\tau)$,
and summing over $\mu$, one gets

$\displaystyle\partial_{\tau}\mathcal{G}_{i}(z,\tau)=\sum_{j}[\Gamma(z)]_{ij}\mathcal{G}_{j}(z,\tau)$
(D3)

where

$\displaystyle\Gamma(z)$
$\displaystyle\equiv\begin{bmatrix}-(\bar{n}+1)&-i\Omega_{\gamma}&i\Omega_{\gamma}&\bar{n}e^{-z}\\\
-i\Omega_{\gamma}&i\delta\omega_{\gamma}-(2\bar{n}+1)/2&0&i\Omega_{\gamma}\\\
i\Omega_{\gamma}&0&-i\delta\omega_{\gamma}-(2\bar{n}+1)/2&-i\Omega_{\gamma}\\\
(\bar{n}+1)e^{z}&i\Omega_{\gamma}&-i\Omega_{\gamma}&-\bar{n}\end{bmatrix}.$
(D4)

In obtaining the matrix $\Gamma(z)$ where the factor $e^{zn_{\mu}}$ is
multiplied to each element, we have set $n_{\mu}=-1$ for the matrix element
$[\Gamma(z)]_{ee,gg}$ corresponding to absorption, $n_{\mu}=1$ for
$[\Gamma(z)]_{gg,ee}$ corresponding to emission, and $n_{\mu}=0$ for the rest.
The formal solution of Eq. D3 can be written as

$\displaystyle\mathcal{G}_{i}(z,\tau)=\sum_{\alpha}T_{\alpha
i}(z)e^{\lambda_{\alpha}(z)\tau}$ (D5)

where $\lambda_{\alpha}(z)$ and $T_{\alpha i}(z)$ are the $\alpha$-th
eigenvalue and the corresponding eigenstate of $\Gamma(z)$, satisfying
$\sum_{j}\left[\Gamma(z)\right]_{ij}T_{\alpha j}=\lambda_{\alpha}(z)T_{\alpha
i}(z)$ with $\lambda_{0}(z)>\lambda_{1}(z)>\lambda_{2}(z)>\lambda_{3}(z)$
($\alpha=0,1,2,3$). At steady state, $\mathcal{G}_{i}(z,\tau)$ is dominated by
the term with the largest eigenvalue

$\displaystyle\lim_{\tau\rightarrow\infty}\mathcal{G}_{i}(z,\tau)\sim
T_{0i}(z,\tau)e^{\lambda_{0}(z)\tau}.$ (D6)

The eigenvalues $\lambda_{\alpha}(z)$ of $\Gamma(z)$ are obtained from the
characteristic polynomial

$\displaystyle\det|\lambda(z)\mathcal{I}-\Gamma(z)|=\sum_{n=0}^{4}a_{n}(z)\lambda^{n}(z)=0,$
(D7)

where

$\displaystyle a_{4}$ $\displaystyle=1$ $\displaystyle a_{3}$
$\displaystyle=2(2\bar{n}+1)$ $\displaystyle a_{2}$
$\displaystyle=\frac{1}{4}\left[5(2\bar{n}+1)^{2}+4\delta\omega_{\gamma}^{2}+16\Omega_{\gamma}^{2}\right]$
$\displaystyle a_{1}(z)$
$\displaystyle=\frac{1}{4}\left[(2\bar{n}+1)^{3}+4\delta\omega_{\gamma}^{2}(2\bar{n}+1)+16\Omega_{\gamma}^{2}(2\bar{n}+1)-8\Omega_{\gamma}^{2}\left(e^{z}+\bar{n}e^{-z}+\bar{n}e^{z}\right)\right]$
$\displaystyle a_{0}(z)$
$\displaystyle=-\Omega_{\gamma}^{2}\left[e^{z}+\bar{n}e^{-z}+2\bar{n}^{2}e^{-z}+3\bar{n}e^{z}+2\bar{n}^{2}e^{z}-(2\bar{n}+1)^{2}\right].$

Now, summing $\mathcal{G}_{i}(z,\tau)$ over the index $i\in\\{ee,eg,ge,gg\\}$,
we obtain the following expression that can be used as the moment generating
function:

$\displaystyle\mathcal{G}(z,\tau)\equiv\sum_{i}\mathcal{G}_{i}(z,\tau)=\sum_{i}\sum_{\mu=-\infty}^{\infty}e^{zn_{\mu}}\tilde{\rho}_{i}(\mu,\tau)$
(D8)

which allows us to calculate the first and second cumulant of $n(\tau)$ as
follows.

$\displaystyle\langle
n(\tau)\rangle=\frac{\sum_{i}\sum_{\mu}n_{\mu}\tilde{\rho}_{i}(\mu,\tau)}{\sum_{i}\sum_{\mu}\tilde{\rho}_{i}(\mu,\tau)}=\partial_{z}\log{\mathcal{G}(z,\tau)}|_{z=0}$
(D9)

and

$\displaystyle\langle n(\tau)^{2}\rangle-\langle
n(\tau)\rangle^{2}=\partial_{z}^{2}\log{\mathcal{G}(z,\tau)}|_{z=0}.$ (D10)

Finally, from Eqs. D9 and D10, and the asymptotic expression of
$\mathcal{G}(z,\tau)$,

$\displaystyle\lim_{\tau\rightarrow\infty}\mathcal{G}(z,\tau)\sim
h_{0}(z,\tau)e^{\lambda_{0}(z)\tau}.$ (D11)

with $h_{0}(z,\tau)\equiv\sum_{i}T_{0i}(z,\tau)$, we can obtain $\langle
j\rangle$ and $\text{Var}[j]$,

$\displaystyle\langle j\rangle\equiv\lim_{\tau\rightarrow\infty}\frac{\langle
n(\tau)\rangle}{\tau}=\lambda_{0}^{\prime}(0)$ (D12)

and

$\displaystyle\text{Var}[j]\equiv\lim_{\tau\rightarrow\infty}\frac{\langle
n(\tau)^{2}\rangle-\langle
n(\tau)\rangle^{2}}{\tau}=\lambda_{0}^{\prime\prime}(0).$ (D13)

These key values of $\lambda^{\prime}_{0}(0)$ and
$\lambda^{\prime\prime}_{0}(0)$ can be evaluated in terms of the coefficients
of the characteristic polynomial (Eq. D7) differentiated with respect to $z$
at $z=0$ as follows:

$\displaystyle a_{0}^{\prime}(0)+a_{1}(0)\lambda_{0}^{\prime}(0)=0$ (D14)

and

$\displaystyle
a_{0}^{\prime\prime}(0)+2a_{1}^{\prime}(0)\lambda_{0}^{\prime}(0)+a_{1}(0)\lambda_{0}^{\prime\prime}(0)+2a_{2}(0)(\lambda_{0}^{\prime}(0))^{2}=0.$
(D15)

Therefore,

$\displaystyle\langle
j\rangle=\lambda_{0}^{\prime}(0)=-\frac{a_{0}^{\prime}(0)}{a_{1}(0)}$ (D16)

and

$\displaystyle\text{Var}[j]=\lambda_{0}^{\prime\prime}(0)$
$\displaystyle=-\frac{1}{a_{1}(0)}\left(a_{0}^{\prime\prime}(0)+2a_{1}^{\prime}(0)\lambda_{0}^{\prime}(0)+2a_{2}(0)(\lambda_{0}^{\prime}(0))^{2}\right).$
(D17)

## References

  * Barato and Seifert [2015] A. C. Barato and U. Seifert, Thermodynamic Uncertainty Relation for Biomolecular Processes, Phys. Rev. Lett. 114, 158101 (2015).
  * Gingrich  _et al._ [2016] T. R. Gingrich, J. M. Horowitz, N. Perunov, and J. L. England, Dissipation Bounds All Steady-State Current Fluctuations, Phys. Rev. Lett. 116, 120601 (2016).
  * Dechant and Sasa [2018a] A. Dechant and S.-i. Sasa, Entropic bounds on currents in langevin systems, Phys. Rev. E 97, 062101 (2018a).
  * Horowitz and Gingrich [2020] J. M. Horowitz and T. R. Gingrich, Thermodynamic uncertainty relations constrain non-equilibrium fluctuations, Nature Physics 16, 15 (2020).
  * Pietzonka  _et al._ [2017] P. Pietzonka, F. Ritort, and U. Seifert, Finite-time generalization of the thermodynamic uncertainty relation, Phys. Rev. E. 96, 012101 (2017).
  * Horowitz and Gingrich [2017] J. M. Horowitz and T. R. Gingrich, Proof of the finite-time thermodynamic uncertainty relation for steady-state currents, Phys. Rev. E 96, 020103(R) (2017).
  * Dechant and Sasa [2018b] A. Dechant and S.-i. Sasa, Current fluctuations and transport efficiency for general langevin systems, J. Stat. Mech. , 063209 (2018b).
  * Proesmans and Van den Broeck [2017] K. Proesmans and C. Van den Broeck, Discrete-Time Thermodynamic-like Uncertainty Relation, EPL 119, 20001 (2017).
  * Pietzonka  _et al._ [2016] P. Pietzonka, A. C. Barato, and U. Seifert, Universal bound on the efficiency of molecular motors, J. Stat. Mech. Theory Exp. , 124004 (2016).
  * Hwang and Hyeon [2018] W. Hwang and C. Hyeon, Energetic Costs, Precision, and Transport Efficiency of Molecular Motors, J. Phys. Chem. Lett. 9, 513 (2018).
  * Marsland III  _et al._ [2019] R. Marsland III, W. Cui, and J. M. Horowitz, The thermodynamic uncertainty relation in biochemical oscillations, J. Roy. Soc. Interface 16, 20190098 (2019).
  * Song and Hyeon [2020] Y. Song and C. Hyeon, Thermodynamic cost, speed, fluctuations, and error reduction of biological copy machines, J. Phys. Chem. Lett. 11, 3136 (2020).
  * Mugnai  _et al._ [2020] M. L. Mugnai, C. Hyeon, M. Hinczewski, and D. Thirumalai, Theoretical perspectives on biological machines, Rev. Mod. Phys. 92, 025001 (2020).
  * Kim and Hyeon [2021] P. Kim and C. Hyeon, Thermodynamic optimality of glycolytic oscillations, J. Phys. Chem. B 154, 214901 (2021).
  * Song and Hyeon [2021] Y. Song and C. Hyeon, Thermodynamic uncertainty relation to assess biological processes, J. Chem. Phys. 154, 130901 (2021).
  * Dechant [2018] A. Dechant, Multidimensional thermodynamic uncertainty relations, J. Phys. A: Math. Theor. 52, 035001 (2018).
  * Pietzonka and Seifert [2018] P. Pietzonka and U. Seifert, Universal Trade-Off between Power, Efficiency, and Constancy in Steady-State Heat Engines, Phys. Rev. Lett. 120, 190602 (2018).
  * Holubec and Ryabov [2018] V. Holubec and A. Ryabov, Cycling Tames Power Fluctuations near Optimum Efficiency, Phys. Rev. Lett. 121, 120601 (2018).
  * Dechant and Sasa [2018c] A. Dechant and S.-i. Sasa, Current fluctuations and transport efficiency for general Langevin systems, J. Stat. Mech. 2018, 063209 (2018c).
  * Pigolotti  _et al._ [2017] S. Pigolotti, I. Neri, E. Roldan, and F. Julicher, Generic Properties of Stochastic Entropy Production, Phys. Rev. Lett. 119, 140604 (2017).
  * Hyeon and Hwang [2017] C. Hyeon and W. Hwang, Physical insight into the thermodynamic uncertainty relation using Brownian motion in tilted periodic potentials, Phys. Rev. E 96, 012156 (2017).
  * Brandner  _et al._ [2018] K. Brandner, T. Hanazato, and K. Saito, Thermodynamic Bounds on Precision in Ballistic Multiterminal Transport, Phys. Rev. Lett. 120, 090601 (2018).
  * Lee  _et al._ [2018] S. Lee, C. Hyeon, and J. Jo, Thermodynamic uncertainty relation of interacting oscillators in synchrony, Phys. Rev. E 98, 032119 (2018).
  * Li  _et al._ [2019] J. Li, J. M. Horowitz, T. R. Gingrich, and N. Fakhri, Quantifying dissipation using fluctuating currents, Nature Commun. 10, 1666 (2019).
  * Macieszczak  _et al._ [2018] K. Macieszczak, K. Brandner, and J. P. Garrahan, Unified Thermodynamic Uncertainty Relations in Linear Response, Phys. Rev. Lett. 121, 130601 (2018).
  * Timpanaro  _et al._ [2019] A. M. Timpanaro, G. Guarnieri, J. Goold, and G. T. Landi, Thermodynamic Uncertainty Relations from Exchange Fluctuation Theorems, Phys. Rev. Lett. 123, 090604 (2019).
  * Chiuchiu and Pigolotti [2018] D. Chiuchiu and S. Pigolotti, Mapping of uncertainty relations between continuous and discrete time, Phys. Rev. E 97, 032109 (2018).
  * Koyuk  _et al._ [2019] T. Koyuk, U. Seifert, and P. Pietzonka, A generalization of the thermodynamic uncertainty relation to periodically driven systems, J. Phys. A: Math. Theor. 52, 02LT02 (2019).
  * Barato  _et al._ [2018] A. C. Barato, R. Chetrite, A. Faggionato, and D. Gabrielli, Bounds on current fluctuations in periodically driven systems, New J. Phys. 20, 103023 (2018).
  * Hasegawa and Van Vu [2019] Y. Hasegawa and T. Van Vu, Fluctuation theorem uncertainty relation, Phys. Rev. Lett. 123, 110602 (2019).
  * Fischer  _et al._ [2020] L. P. Fischer, H.-M. Chun, and U. Seifert, Free diffusion bounds the precision of currents in underdamped dynamics, Phys. Rev. E 102, 012120 (2020).
  * Koyuk and Seifert [2019] T. Koyuk and U. Seifert, Operationally accessible bounds on fluctuations and entropy production in periodically driven systems, Phys. Rev. Lett. 122, 230601 (2019).
  * Chun  _et al._ [2019] H.-M. Chun, L. P. Fischer, and U. Seifert, Effect of a magnetic field on the thermodynamic uncertainty relation, Phys. Rev. E 99, 042128 (2019).
  * Lee  _et al._ [2019] J. S. Lee, J.-M. Park, and H. Park, Thermodynamic uncertainty relation for underdamped langevin systems driven by a velocity-dependent force, Phys. Rev. E 100, 062132 (2019).
  * Carollo  _et al._ [2019] F. Carollo, R. L. Jack, and J. P. Garrahan, Unraveling the large deviation statistics of markovian open quantum systems, Phys. Rev. Lett. 122, 130605 (2019).
  * Hasegawa [2021] Y. Hasegawa, Thermodynamic uncertainty relation for general open quantum systems, Phys. Rev. Lett. 126, 010602 (2021).
  * Ptaszynski [2018] K. Ptaszynski, Coherence-enhanced constancy of a quantum thermoelectric generator, Phys. Rev. B 98, 085425 (2018).
  * Agarwalla and Segal [2018] B. K. Agarwalla and D. Segal, Assessing the validity of the thermodynamic uncertainty relation in quantum systems, Phys. Rev. B 98, 155438 (2018).
  * Liu and Segal [2019] J. Liu and D. Segal, Thermodynamic uncertainty relation in quantum thermoelectric junctions, Phys. Rev. E 99, 062141 (2019).
  * Saryal  _et al._ [2021] S. Saryal, O. Sadekar, and B. K. Agarwalla, Thermodynamic uncertainty relation for energy transport in a transient regime: A model study, Phys. Rev. E 103, 022141 (2021).
  * Rignon-Bret  _et al._ [2021] A. Rignon-Bret, G. Guarnieri, J. Goold, and M. T. Mitchison, Thermodynamics of precision in quantum nanomachines, Phys. Rev. E 103, 012133 (2021).
  * Breuer and Petruccione [2010] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, USA, 2010).
  * Carmichael [2002] H. J. Carmichael,  _Statistical Methods in Quantum optics 1_ (Springer, 2002).
  * Leggett  _et al._ [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. Fisher, A. Garg, and W. Zwerger, Dynamics of the dissipative two-state system, Rev. Mod. Phys. 59, 1 (1987).
  * Scully and Zubairy [1997] M. O. Scully and M. S. Zubairy,  _Quantum Optics_ (Cambridge University Press, 1997).
  * Tscherbul and Brumer [2014] T. V. Tscherbul and P. Brumer, Long-Lived Quasistationary Coherences in a V-type System Driven by Incoherent Light, Phys. Rev. Lett. 113, 113601 (2014).
  * Koza [1999] Z. Koza, General technique of calculating the drift velocity and diffusion coefficient in arbitrary periodic systems, J. Phys. A. 32, 7637 (1999).
  * Bruderer  _et al._ [2014] M. Bruderer, L. D. Contreras-Pulido, M. Thaller, L. Sironi, D. Obreschkow, and M. B. Plenio, Inverse counting statistics for stochastic and open quantum systems: the characteristic polynomial approach, New J. Phys. 16, 033030 (2014).
  * Yang and Cao [2020] P.-Y. Yang and J. Cao, Steady-State Analysis of Light-Harvesting Energy Transfer Driven by Incoherent Light: From Dimers to Networks, J. Phys. Chem. Lett. 11, 7204 (2020).
  * Wu  _et al._ [2012] J. Wu, F. Liu, J. Ma, R. J. Silbey, and J. Cao, Efficient energy transfer in light-harvesting systems: Quantum-classical comparison, flux network, and robustness analysis, J. Chem. Phys. 137, 174111 (2012).
  * Breuer [2003] H.-P. Breuer, Quantum jumps and entropy production, Phys. Rev. A 68, 032105 (2003).
  * Menczel  _et al._ [2021] P. Menczel, E. Loisa, K. Brandner, and C. Flindt, Thermodynamic uncertainty relations for coherently driven open quantum systems, J. Phys. A: Math. Theor. 54, 314002 (2021).
  * Thingna  _et al._ [2012] J. Thingna, J.-S. Wang, and P. Hänggi, Generalized Gibbs state with modified Redfield solution: Exact agreement up to second order, J. Chem. Phys. 136, 194110 (2012).
  * Tanimura [2020] Y. Tanimura, Numerically “exact” approach to open quantum dynamics: The hierarchical equations of motion (HEOM), J. Chem. Phys. 153, 020901 (2020).
  * Ikeda and Sato [2020] T. N. Ikeda and M. Sato, General description for nonequilibrium steady states in periodically driven dissipative quantum systems, Sci. Adv. 6, eabb4019 (2020).