# Demon driven by geometric phase

Ryosuke Yoshii<EMAIL_ADDRESS>Center for Liberal Arts and Sciences,
Sanyo-Onoda City University, Yamaguchi 756-0884, Japan  Hisao Hayakawa
<EMAIL_ADDRESS>Center for Gravitational Physics and Quantum
Information, Yukawa Institute for Theoretical Physics, Kyoto University,
Kitashirakawa-oiwake cho, Sakyo-ku, Kyoto 606-8502, Japan

(August 27, 2024)

###### Abstract

We theoretically study the entropy production and the work extracted from a
system connected to two reservoirs by periodic modulations of the
electrochemical potentials of the reservoirs and the parameter of a system
Hamiltonian under isothermal conditions. We find that the modulation of the
parameters can drive a geometric state, which is away from a nonequilibrium
steady state. Using this property, we construct a demon in which the entropy
production during the first one-cycle is negative such that we can extract the
work if we start from the nonequilibrium steady state without parameter
modulations. We use the Anderson model to implement the demon in a realistic
situation.

##  I Introduction

The second law of thermodynamics is one of the most fundamental laws of
physics, specifying the upper limit of the available work that can be
extracted from reservoirs. Maxwell proposed an ideal setup to violate the
second law, in which a demon quickly opens and closes the gate, allowing only
fast-moving molecules to pass through in one direction [1]. This leads to a
decrease in entropy without any work being done and violates the second law of
thermodynamics. A simplified version of Maxwell’s demon was proposed by
Szilard [2]. It was shown that Maxwell’s demon can be used to implement the
Szilard engine, while the original Maxwell’s demon makes the temperature
gradient act like a refrigerator. Since then, the demon problem has attracted
much attention from physicists [3, 4, 5, 6, 7, 8]. Modern experiments use
photonic devices to realise Maxwell’s demon [9, 10, 11]. A comprehensive
historical review of Maxwell’s demon is given in Ref. [12].

Since Maxwell’s original idea was based on the measurement of molecules, it is
natural to combine the physical science with information thermodynamics to
realise Maxwell’s demon [3, 4, 5, 6, 7, 8, 9, 13, 14, 15, 16, 17, 18].
However, the cost of implementing the informational Maxwell’s demon is high,
although the theoretical formulation ignores this cost. Instead, we propose a
geometric demon using Berry’s phase [19] in a geometric engine as an extension
of the Thouless pumping [20, 21, 22]. We consider a small system sandwiched
between two thermal reservoirs. If two parameters in the reservoirs and one
parameter in the system Hamiltonian are controlled by an external agent, we
can extract the work from the system. This is a natural application of the
Thouless pumping [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
35, 36] and geometric thermodynamics [37, 38, 39, 40, 41, 43, 44, 42].

It is well known that the Kullback-Leibler (KL) divergence is non-negative,
being zero only when the system is in a nonequilibrium steady state (NESS)
[45, 46, 47, 48, 49, 50, 43, 44, 51]. It is also known that the KL divergence
cannot increase for all completely positive and trace-preserving (CPTP)
processes [45, 46, 47, 48, 49, 50]. As a consequence, the KL divergence
decreases monotonically towards zero in the relaxation process. This
monotonicity of the KL divergence corresponds to the second law of
thermodynamics or the H-theorem. Nevertheless, we must be careful in choosing
the correct KL divergence in the presence of modulations. Indeed, we show that
the correct choice of KL divergence under a periodic modulation is not the
relative entropy between the time-evolving density matrix and that in a steady
state but between the matrix and that in a geometric state because of the
generation of a geometric term (see Fig. 1).

We can therefore extract the work done by this geometric engine by means of an
increase in the relative entropy corresponding to a negative entropy
production. Let us call such an engine a geometric demon. The advantage of the
geometric demon is that we do not need any costs for the feedback control, but
we do need a non-adiabatic control of the system.

Figure 1: A schematic of the process used for the implementation of the
geometric demon. At $\theta=0$, we start modulating the parameters. Although
the modulation is slow enough to prevent the system from being excited, the
geometric term is suddenly added to the state vector. The state relaxes to the
geometric state as time goes. The discrepancy between the steady state in the
absence of the parameter modulation $|\hat{\rho}^{\rm SS}\rangle$ and the
geometric state $|\hat{\rho}^{\rm GS}\rangle$ can be used to determine the
work.

The organization of this paper is as follows. In Sec. II, we explain the setup
and the geometric state under the parameter modulation. We also present the
corresponding KL divergences associated including the geometric contribution.
In Sec. III, we discuss the work, the heat, and the efficiency in the present
setup. In Sec. IV, we apply our formulation to the Anderson model to
demonstrate the protocol for realizing the geometric demon. In Sec. V we
summarize our results and give their physical interpretations. We also give
future perspectives. In Appendix A, we describe the derivation of the time
evolution of the state including the geometric contribution. We also give a
detailed expression for the BSN curvature. In Appendix B, we explain the
differences between the present study and the previous studies. In Appendix C,
we present the detailed properties of the Anderson model. In Appendix D, we
present some detailed properties of the density matrix, the von Neumann
entropy, and the heat flow from the comparison of the results of the Born-like
approximation used in this study with those obtained by the numerical
solution.

##  II Geometric phase and entropy production

In the present study, we focus on a system connected to two reservoirs
depicted in Fig. 2. The left and right reservoirs are characterized by
electrochemical potentials ($\mu_{\rm L}$ and $\mu_{\rm R}$) and temperature
$T$, respectively.

Figure 2: A schematic of the model considered in the present paper. The system
is connected to two reservoirs characterized by chemical potentials and
temperatures.

The electrochemical potentials are modulated as

$\mu_{\mathrm{L}}=\overline{\mu}(1+r_{\mathrm{L}}\sin\theta),\quad\mu_{\mathrm{R}}=\overline{\mu}[1+r_{\mathrm{R}}\sin(\theta+\delta)],$
(1)

where
$\overline{\mu}:=\frac{1}{2\pi}\int_{0}^{2\pi}d\theta\mu_{\alpha}(\theta)$ is
the one-cycle average of the electrochemical potential $\mu_{\alpha}$ in a
reservoir $\alpha(=$${\mathrm{L}}$ or ${\mathrm{R}}$). We assume that
$\mu_{\alpha}$ depends only on the modulation phase $\theta:=\omega t$, where
$\omega$ and $t$ are the modulation angular frequency and time, respectively.
We also assume that the system Hamiltonian $\hat{H}(\lambda(\theta))$ is
perfectly periodic in this paper, i. e.,

$\hat{H}(\lambda(\theta))=\hat{H}(\lambda(\theta+2\pi))$ (2)

through a control parameter $\lambda(\theta)$, where
$\lambda(\theta):=1+r_{H}\cos\theta$. In this paper, we consider only the case
of $r:=r_{\mathrm{L}}=r_{\mathrm{R}}=r_{H}$ in order to reduce the number of
parameters for simplicity. To maintain the positivity of the parameters, we
assume $|r|<1$. Thus, our system is characterized by a set of fixed
parameters, such as $T$ and $\overline{\mu}$, and two control parameters, $r$
and $\delta$. To express a set of control parameters, we introduce

$\bm{\Lambda}(\theta,\delta):=\left(\lambda(\theta),\frac{\mu_{\rm
L}(\theta)}{\overline{\mu}},\frac{\mu_{\rm
R}(\theta,\delta)}{\overline{\mu}}\right)$ (3)

using $\Lambda_{\mu}$ as one of its components. Figure 3 shows a schematic of
the control parameters for $r=1$ and $\delta=\pi/4$.

Figure 3: A schematic of control parameters for $r=1$ and $\delta=\pi/4$.

We consider the master equation for the density matrix
$\hat{\rho}(\theta,\delta)$:

$\displaystyle\frac{d}{d\theta}|\hat{\rho}(\theta,\delta)\rangle=\epsilon^{-1}\hat{K}(\theta,\delta)|\hat{\rho}(\theta,\delta)\rangle,$
(4)

where $\hat{K}(\theta,\delta)$ is the evolution operator. We use the vector
notation $|\hat{\rho}(\theta,\delta)\rangle$ in Eq. (4), where the components
of the density matrix $\hat{\rho}(\theta,\delta)$ are aligned. We introduce
the dimensionless parameter $\epsilon$ in Eq. (4) as
$\epsilon:=2\pi\omega/\Gamma$ where $\Gamma$ characterizes the hopping rate of
the electrons from the reservoirs into the system. Introducing the eigenvalue
$\varepsilon_{i}$ with the subscript $i$ for $\hat{K}$, the corresponding left
and right eigenstates $\langle\ell_{i}|$ and $|r_{i}\rangle$ satisfy the
orthonormal relation $\langle\ell_{i}|r_{j}\rangle=\delta_{ij}$, if the
eigenvalues are non-degenerate. We assume that there exists a non-degenerate
largest eigenvalue $\varepsilon_{0}:=0$ corresponding to a stationary state.
Probability conservation leads to the left-zero eigenvector
$\langle\ell_{0}|$, which is defined as $\langle\ell_{0}|\hat{K}=0$, whose
diagonal components are $1$ in the matrix form and $0$ otherwise. The right
zero eigenstate $|r_{0}\rangle$ satisfying $\hat{K}|r_{0}\rangle=0$ is also
expressed as $|\hat{\rho}^{\rm SS}\rangle$ to specify the NESS.

In the absence of the modulations, the physical state is relaxed to
$|\hat{\rho}^{\rm SS}\rangle$, since all eigenvalues except for
$\varepsilon_{0}=0$ are negative. This suggests that the natural choice of the
initial state would be $|\hat{\rho}^{\rm SS}\rangle$. Therefore, we analyze
the case where the initial state is the steady state
$|\hat{\rho}(0,\delta)\rangle=|r_{0}\rangle=|\hat{\rho}^{\mathrm{SS}}\rangle$.

As soon as we start the modulation, the state relaxes to the other state,
which we call the geometric state. In the presence of modulation, the state
vector $|\hat{\rho}\rangle$ contains the contribution of a geometric phase. As
shown in Appendix A.1, the approximate expression for $|\hat{\rho}\rangle$ is
given by

$\displaystyle|\hat{\rho}(\theta,\delta)\rangle\simeq|\hat{\rho}^{\rm
SS}(\theta,\delta)\rangle+\sum_{i\neq
0}^{n}C_{i}(\theta,\delta)|r_{i}(\theta,\delta)\rangle,$ (5)

where

$\displaystyle C_{i}(\theta,\delta)$ $\displaystyle=-\int_{0}^{\theta}d\phi
e^{\epsilon^{-1}\int^{\theta}_{\phi}dz\varepsilon_{i}(z,\delta)}\langle\ell_{i}(\phi,\delta)|\frac{d}{d\phi}|r_{0}(\phi,\delta)\rangle.$
(6)

The second term on the RHS of Eq. (5) is the Berry-Sinitsyn-Nemenman (BSN)
connection [21, 22]. A similar expression is derived in Ref. [36] ( see
Appendix B for details of the relationship between our analysis and that in
Ref. [36]). Note that Eqs. (5) and (6) describe the evolution process from the
initial steady state $|\hat{\rho}^{\mathrm{SS}}\rangle$, but the expression
can be used if we are interested in the long time ($\theta\gg\epsilon$)
behavior, where the choice of the initial condition is not important. We also
present the expression for $|\hat{\rho}(\theta,\delta)\rangle$ starting from a
general initial condition in Appendix A.1 as well as the relaxation to Eqs.
(5) and (6) for $\theta\gg\epsilon$.

The state described by Eqs. (5) and (6) relaxes to the geometric state:

$\displaystyle|\hat{\rho}^{\rm{GS}}(\theta,\delta)\rangle:=\lim_{n\to\infty}\mathcal{T}\exp\left(\int^{\theta+2n\pi}_{0}d\phi\hat{K}(\phi,\delta)\right)|r_{0}(0,\delta)\rangle,$
(7)

where we have introduced the time-ordering operator $\cal{T}$ defined as
$\mathcal{T}[\hat{A}(t_{2})\hat{B}(t_{1})\cdots\hat{Z}(t_{n})]:=\hat{Z}(t_{n})\cdots\hat{A}(t_{2})\hat{B}(t_{1})$
with $t_{1}<t_{2}<\cdots<t_{n}$. Since this relaxation process is fast, we can
use the approximate expression for the geometric state as

$|\hat{\rho}^{\rm{GS}}(\theta,\delta)\rangle\simeq\mathcal{T}\exp\left(\int^{\theta+2\pi}_{0}d\phi\hat{K}(\phi,\delta)\right)|r_{0}(0,\delta)\rangle,$
(8)

Note that the trace preservation is always satisfied for an arbitrary $\theta$
from $\langle\ell_{0}|r_{i}\rangle=\delta_{0i}$ and
$\langle\ell_{0}|\hat{K}=0$ in the dynamics described by Eq. (4). In this
paper, we restrict our interest to quasi-adiabatic processes. There is an
analogy between Eq. (5) and the path-integral with the Born approximation in
the quantum theory. The phase factor
$e^{\epsilon^{-1}\int^{\theta}_{\phi}dz\varepsilon_{i}(z,\delta)}$ in Eq. (6)
can be interpreted as free relaxation instead of free propagation. If the
modulation is not slow, i.e. $\epsilon$ is larger than $|\varepsilon_{i\neq
0}|$, it would not be adequate to ignore the higher order processes. The
expressions in Eqs. (5) and (6) are compatible with the slow-modulation
approximation employed in Ref. [44]. The leading contribution of the
modulation to the entropy production is indicated by the second term on the
right-hand side (RHS) of Eq. (5).

As shown in Appendix A.2, the deviation from the initial state after a
modulation cycle satisfying
$|r_{i}(2\pi,\delta)\rangle=|r_{i}(0,\delta)\rangle$ is written as

$\displaystyle\Delta|\hat{\rho}\rangle:=|\hat{\rho}(2\pi,\delta)\rangle-|\hat{\rho}(0,\delta)\rangle=\sum_{i\neq
0}\mathcal{C}_{i}(\delta)|r_{i}(0,\delta)\rangle$ (9)

with

$\displaystyle\mathcal{C}_{i}(\delta):=\int_{0}^{2\pi}d\phi
e^{\epsilon^{-1}\int^{2\pi}_{\phi}dz\varepsilon_{i}(z,\delta)}\mathcal{A}^{\mu}_{i}\frac{\partial\Lambda_{\mu}}{\partial\phi},$
(10)

where we have introduced the BSN connection $\mathcal{A}^{\mu}_{i}$ as

$\displaystyle\mathcal{A}^{\mu}_{i}(\phi,\delta):=-\langle\ell_{i}(\Lambda(\phi,\delta))|\frac{\partial}{\partial\Lambda_{\mu}}|r_{0}(\Lambda(\phi,\delta))\rangle.$
(11)

According to $\mathcal{A}^{\mu}_{i}$, we define the BSN curvature as

$\displaystyle
F^{\mu\nu}_{i}(\theta,\delta):=\left(\frac{\partial\mathcal{A}^{\nu}_{i}}{\partial\Lambda_{\mu}}\right)_{\theta}-\left(\frac{\partial\mathcal{A}^{\mu}_{i}}{\partial\Lambda_{\nu}}\right)_{\theta}.$
(12)

Due to the damping factor in Eq. (10), the contribution of the BSN curvature
is localized in $\theta$.

Now, let us consider the KL divergence. The proper KL divergence is not
$D^{\mathrm{KL}}(\hat{\rho}(\theta,\delta)||\hat{\rho}^{\mathrm{SS}}(\theta,\delta))$
but
$D^{\mathrm{KL}}(\hat{\rho}(\theta,\delta)||\hat{\rho}^{\mathrm{GS}}(\theta,\delta))$,
where we have defined the KL divergence [45, 46, 47, 48, 49, 50, 43, 44] as

$\displaystyle D^{\mathrm{KL}}(\hat{\rho}||\hat{\sigma}):=$
$\displaystyle\mathrm{Tr}\left[\hat{\rho}(\ln\hat{\rho}-\ln\hat{\sigma})\right].$
(13)

Since we expect $|\hat{\rho}(\theta,\delta)\rangle$ to quickly relax to
$|\hat{\rho}^{\rm GS}(\theta,\delta)\rangle$, the KL divergence approximately
satisfies

$D^{\mathrm{KL}}(\hat{\rho}(\theta,\delta)||\hat{\rho}^{\mathrm{GS}}(\theta,\delta))\simeq
0$ (14)

for $\theta\geq 2\pi$. The initially positive
$D^{\mathrm{KL}}(\hat{\rho}(\theta,\delta)||\hat{\rho}^{\mathrm{GS}}(\theta,\delta))$
decreases monotonically as $\theta$ approaches zero. This is the H-theorem.

The production of the total entropy associated with the KL divergence should
be

$\displaystyle\Delta\Sigma(\theta,\delta)$
$\displaystyle\quad:=D^{\mathrm{KL}}(\hat{\rho}(0,\delta)||\hat{\rho}^{\mathrm{GS}}(0,\delta))-D^{\mathrm{KL}}(\hat{\rho}(\theta,\delta)||\hat{\rho}^{\mathrm{GS}}(\theta,\delta)).$
(15)

This $\Delta\Sigma(\theta,\delta)$ which is always non-negative could play an
important role.

On the other hand, we can consider the entropy production during the one-cycle
modulation starting from the NESS at $\theta=0$ as

$\displaystyle\Delta S(\delta):$
$\displaystyle=D^{\mathrm{KL}}(\hat{\rho}(0,\delta)||\hat{\rho}^{\mathrm{SS}}(0,\delta))-D^{\mathrm{KL}}(\hat{\rho}(2\pi,\delta)||\hat{\rho}^{\mathrm{SS}}(0,\delta))$
$\displaystyle=-D^{\mathrm{KL}}(\hat{\rho}(2\pi,\delta)||\hat{\rho}^{\mathrm{SS}}(0,\delta)).$
(16)

This $\Delta S(\delta)$ is always non-positive due to the non-negativity of
the KL divergence. We could expect $\Delta S(\delta)\geq 0$, since the KL
divergence is expected to decrease with $\theta$ [45, 46, 47, 48, 49, 50] as

$\Delta S(\delta)=-\int_{0}^{2\pi}d\theta\frac{\partial}{\partial\theta}D^{\rm
KL}(\hat{\rho}(\theta,\delta)||\hat{\rho}^{\mathrm{SS}}(\theta,\delta))\geq
0.$ (17)

However, the compatible value of $\Delta S(\delta)$ in Eq. (17) with the non-
negative KL divergence is zero under the periodic modulation. Nevertheless,
$\Delta S(\delta)$ can be negative for most $\delta$ as will be shown, due to
$\lim_{\theta\to+0}\hat{\rho}^{\mathrm{GS}}(\theta,\delta)\neq\hat{\rho}^{\mathrm{SS}}(0,\delta)$.
This implies that $D^{\rm
KL}(\hat{\rho}(\theta,\delta)||\hat{\rho}^{\mathrm{SS}}(\theta,\delta))$ is
not proper KL divergence. In fact, as will be shown, $D^{\rm
KL}(\hat{\rho}(\theta,\delta)||\hat{\rho}^{\mathrm{SS}}(\theta,\delta))$ is
not monotonic with $\theta$ but oscillates with $\theta$. Thus, we can extract
the work from the negative entropy production $\Delta S(\delta)<0$.

##  III Work and efficiency

###  III.1 First law of thermocynamics

Now, we discuss the thermodynamic relations to construct the geometric demon.
First, we consider the first law of thermodynamics [51, 52, 53]. The internal
energy is given by

$E(\theta,\delta)=\mathrm{Tr}[\hat{H}(\lambda(\theta))\hat{\rho}(\theta,\delta)].$
(18)

The rate of change of internal energy is given by

$\displaystyle\frac{\partial E}{\partial\theta}(\theta,\delta)=$
$\displaystyle\mathrm{Tr}\left[\frac{\partial\hat{H}(\lambda(\theta))}{\partial\theta}\hat{\rho}(\theta,\delta)\right]+\mathrm{Tr}\left[\hat{H}(\lambda(\theta))\frac{\partial\hat{\rho}(\theta,\delta)}{\partial\theta}\right].$
(19)

Here, the first term on the RHS of Eq. (19) corresponds to the power, which
can be rewritten as

$\displaystyle\mathscr{P}(\theta,\delta):={\rm
Tr}\left[\hat{\rho}(\theta,\delta)\frac{\partial\hat{H}(\lambda(\theta))}{\partial\lambda(\theta)}\right]\dot{\lambda}(\theta),$
(20)

where $\dot{\lambda}(\theta):=\frac{d}{d\theta}\lambda(\theta)$. Thus the work
generated from the initial state to $\theta$ becomes

$\displaystyle\mathcal{W}(\theta,\delta):=\int_{0}^{\theta}d\phi\mathscr{P}(\phi,\delta).$
(21)

By definition, the relation
$\mathcal{\dot{W}}(\theta,\delta)=\mathscr{P}(\theta,\delta)$ holds. The work
$W(\delta)$ during the one cycle is defined as

$\displaystyle W(\delta):=\mathcal{W}(2\pi,\delta).$ (22)

The work $W(\delta)$ and power $\mathscr{P}(\theta,\delta)$ can be positive or
negative depending on the situation. A positive $\mathscr{P}(\theta,\delta)$
is interpreted as a power input from the external agent, while a negative
$\mathscr{P}(\theta,\delta)$ can be interpreted as a power loss. The second
term on the RHS of Eq. (19) corresponds to the heat flow

$\displaystyle\mathcal{\dot{Q}}(\theta,\delta):={\rm
Tr}\left[\frac{\partial\hat{\rho}(\theta,\delta)}{\partial\theta}\hat{H}(\lambda(\theta))\right].$
(23)

Thus the heat generated from the initial state to $\theta$ becomes

$\displaystyle\mathcal{Q}(\theta,\delta):=\int^{\theta}_{0}d\phi\mathcal{\dot{Q}}(\phi,\delta).$
(24)

The heat generated during one cycle is defined as

$\displaystyle Q(\delta):=\mathcal{Q}(2\pi,\delta).$ (25)

Thus, we can rewrite Eq. (19) in the well-known form of the first law of
thermodynamics as

$\displaystyle\frac{\partial}{\partial\theta}E(\theta,\delta)=\frac{\partial}{\partial\theta}\mathcal{W}(\theta,\delta)+\frac{\partial}{\partial\theta}\mathcal{Q}(\theta,\delta).$
(26)

The power $\mathscr{P}(\theta,\delta)$, the work $W(\delta)$, the heat flow
$\mathcal{Q}(\theta,\delta)$, and the heat $Q(\delta)$ can be calculated once
we know the density matrix $\hat{\rho}(\theta,\delta)$ as in Eqs. (5) and (6).

###  III.2 Second law of thermocynamics

Because the KL divergence satisfies the second law of thermodynamics or the H
theorem expressed as
$D^{\mathrm{KL}}(\hat{\rho}(\theta,\delta)||\hat{\rho}^{\mathrm{GS}}(\theta,\delta))\geq
0$ and its monotonicity, we can construct the thermodynamics for the driven
system by considering the system entropy and the heat between the system and
the reservoirs. As will be shown, the system entropy differs from the von
Neumann entropy due to the contribution of the driven terms.

Let us consider $\Delta\Sigma(\theta,\delta)$ which satisfies
$\Delta\Sigma(\theta,\delta)=\Delta S_{\rm sys}+\Delta S_{\rm r}$ where
$\Delta S_{\rm sys}$ and $\Delta S_{\rm r}$ express the changes of the system
entropy and the reservoir, respectively. Because $\Delta S_{\rm r}$ is
expressed as $\Delta S_{\rm r}=-\beta\mathcal{Q}$, we obtain [51]

$\displaystyle\Delta\Sigma(\theta,\delta)=\Delta S_{\rm
sys}(\theta,\delta)-\beta\mathcal{Q}(\theta,\delta),$ (27)

where $\beta:=1/T$, $\Delta S_{\rm sys}(\theta,\delta):=S_{\rm
sys}(\theta,\delta)+\mathrm{Tr}[\hat{\rho}(0,\delta)\ln\hat{\rho}(0,\delta)]$
with

$\displaystyle S_{\rm sys}(\theta,\delta):=$ $\displaystyle
S^{\mathrm{vN}}(\theta,\delta)+\beta\mathcal{Q}(\theta,\delta)+\beta\mathcal{\tilde{Q}}(\theta,\delta),$
(28)

where we have introduced von Neumann entropy

$\displaystyle S^{\mathrm{vN}}(\theta,\delta):$
$\displaystyle=-\mathrm{Tr}[\hat{\rho}(\theta,\delta)\ln\hat{\rho}(\theta,\delta)],$
(29)

and

$\displaystyle\beta\mathcal{\tilde{Q}}(\theta,\delta):$
$\displaystyle=\mathrm{Tr}[\hat{\rho}(\theta,\delta)\ln\hat{\rho}^{\rm
GS}(\theta,\delta)]-\mathrm{Tr}[\hat{\rho}(0,\delta)\ln\hat{\rho}^{\rm
GS}(0,\delta)].$ (30)

This $S_{\rm sys}(\theta,\delta)$ is the system entropy which is deviated from
von Neumann entropy $S^{\mathrm{vN}}(\theta,\delta)$. Note that the $S_{\rm
sys}(\theta,\delta)$ is reduced to $S^{\mathrm{vN}}(\theta,\delta)$ when the
parameter in the Hamiltonian is not modulated and the system is relaxed to an
equilibrium state. Since $\hat{\rho}^{\rm GS}(\theta,\delta)$ is reduced to be
$\hat{\rho}^{\rm eq}(\theta,\delta)\propto e^{-\beta\hat{H}}$ in such a
situation, the RHS of Eq. (30) becomes
$\int_{0}^{\theta}d\phi\frac{d}{d\phi}[-\beta\hat{H}\rho]=-\beta\int_{0}^{\theta}d\phi[\hat{H}\dot{\rho}]=-\beta\mathcal{Q}(\theta,\delta)$.
Since we have already obtained the expression of $\hat{\rho}(\theta,\delta)$
as in Eqs. (5) and (6), one can obtain the heat $\mathcal{Q}(\theta,\delta)$.

###  III.3 Absorbing heat, work, and efficiency

Because $\dot{\mathcal{Q}}(\theta,\delta)$ can be negative, the heat
$\mathcal{Q}(\theta,\delta)$ can be negative. We can decompose
$\mathcal{Q}(\theta,\delta)$ into the absorbing heat and the releasing heat as

$\displaystyle Q_{A}(\delta):$
$\displaystyle=\int_{0}^{2\pi}d\theta\dot{\mathcal{Q}}(\theta,\delta)\Theta(\dot{\mathcal{Q}}(\theta,\delta)),$
(31) $\displaystyle Q_{R}(\delta):$ $\displaystyle=Q(\delta)-Q_{A}(\delta),$
(32)

where $\Theta(x)$ is the Heaviside step function that satisfies $\Theta(x)=1$
if $x\geq 0$ and $\Theta(x)=0$ otherwise. If we consider a slow modulation
process, we expect the deviation of the state from the steady state
$\hat{\rho}^{\rm SS}$ to be small and the time evolution will be nearly
periodic. Thus, $\Delta E(\delta):=E(2\pi,\delta)-E(0,\delta)$ would be much
smaller than the other variables such as $W(\delta)$ and $Q(\delta)$. (Namely,
we expect $|W(\delta)|,|Q(\delta)|\gg|\Delta E(\delta)|$). This implies that
$Q_{A}(\delta)\geq|Q(\delta)|\simeq|W(\delta)|\geq W(\delta)$ holds due to the
first law of thermodynamics. We can prove this relation by considering a slow
modulation process $\epsilon\ll 1$ under the rotational wave approximation
[54]. 111 We estimate $\Delta E(\delta)$ as follows. $\Delta E(\delta)={\rm
Tr}[\hat{H}\hat{\rho}(2\pi,\delta)]-{\rm
Tr}[\hat{H}\hat{\rho}(0,\delta)]=\sum_{i,j}[\hat{H}_{ij}\hat{\rho}^{\rm
GS}_{ji}-\hat{H}_{ij}\hat{\rho}^{\rm SS}_{ji}]$. Since both $\hat{\rho}^{\rm
GS}_{ij}$ and $\hat{\rho}^{\rm SS}_{ij}$ are diagonal under the rotational
wave approximation [54], $\Delta E(\delta)$ is estimated as $\Delta
E(\delta)=\sum_{i}[\hat{H}_{ii}\hat{\rho}^{\rm
GS}_{ii}-\hat{H}_{ii}\hat{\rho}^{\rm
SS}_{ii}]=\sum_{i}[\hat{H}_{ii}(\hat{\rho}^{\rm
SS}_{ii}+C_{1}r_{1,i}+\cdots)-\hat{H}_{ij}\hat{\rho}^{\rm SS}_{ii}]\simeq
C_{1}\sum_{i}[\hat{H}_{ii}r_{1,i}]$, where $|r_{1}\rangle$ is the right
eigenvector of $\hat{K}$ associated with the slowest eigenmode.
$\dot{\mathcal{Q}}={\rm Tr}[\hat{H}\dot{\hat{\rho}}]={\rm
Tr}[\hat{H}\epsilon^{-1}\mathcal{K}\hat{\rho}]\simeq C_{1}{\rm
Tr}[\hat{H}\epsilon^{-1}\varepsilon_{1}\hat{r}_{1}]=C_{1}\epsilon^{-1}\varepsilon_{1}\sum_{i}[\hat{H}_{ii}r_{1,i}]$,
and $\mathscr{P}=\int d\theta{\rm
Tr}[\hat{H}\dot{\hat{\rho}}]=\epsilon^{-1}\int d\theta
C_{1}\varepsilon_{1}\sum_{i}[\hat{H}_{ii}r_{1,i}]$, where $\hat{r}_{1}$ is the
matrix expression of $|r_{1}\rangle$. Thus, $\dot{\mathcal{Q}}$ and
$\mathscr{P}$ are much larger than $\Delta E$ because of the factor
$\epsilon^{-1}$ in the quasi-adiabatic situation.  If the work $W(\delta)$ is
negative, the system can be regarded as an engine in which the work done by
the system is larger than the work done by the reservoirs.

Thus, to construct a geometric demon, we require that $W(\delta)<0$ to use the
negative entropy production $\Delta S(\delta)<0$. To the best of our
knowledge, we cannot determine the sign of $W(\delta)$ in general. Therefore,
we use the Anderson model to show that $W(\delta)$ can be negative.

Using the absorbed energy, the work, and the efficiency can be introduced as

$\eta(\delta):=-\frac{W(\delta)}{Q_{\mathrm{A}}(\delta)}.$ (33)

When the work becomes positive, the efficiency becomes ill-defined.

##  IV Application to the Anderson model

Let us apply the general framework to the Anderson model for a quantum dot
(QD), where a single dot is coupled to two electron reservoirs. 222Although
higher energy levels exist in realistic quantum dot systems, for simplicity
here we consider only the single energy level for simplicity. The total
Hamiltonian $\hat{H}^{\rm tot}$ can be written as

$\displaystyle\hat{H}^{\rm tot}:$ $\displaystyle=\hat{H}+\hat{H}^{\rm
r}+\hat{H}^{\rm int},$ (34)

where $\hat{H}$, the reservoir Hamiltonian $\hat{H}^{\rm r}$ and the
interaction Hamiltonian $\hat{H}^{\rm int}$ are, respectively, given by

$\displaystyle\hat{H}$
$\displaystyle=\sum_{\sigma}\epsilon_{0}\hat{d}_{\sigma}^{\dagger}\hat{d}_{\sigma}+U(\theta)\hat{n}_{\uparrow}\hat{n}_{\downarrow},$
(35) $\displaystyle\hat{H}^{\rm r}$
$\displaystyle=\sum_{\alpha,k,\sigma}\epsilon_{k}\hat{a}_{\alpha,k,\sigma}^{\dagger}\hat{a}_{\alpha,k,\sigma},$
(36) $\displaystyle\hat{H}^{\rm int}$
$\displaystyle=\sum_{\alpha,k,\sigma}V_{\alpha}\hat{d}_{\sigma}^{\dagger}\hat{a}_{\alpha,k,\sigma}+{\rm
h.c.},$ (37)

where $\hat{a}_{\alpha,k,\sigma}^{\dagger}$ and $\hat{a}_{\alpha,k,\sigma}$
are the creation and annihilation operators, respectively, for the electrons
in the reservoirs $\alpha(=$L or R) with the wavenumber $k$, energy
$\epsilon_{k}$, and spin $\sigma(=\uparrow$ or $\downarrow$). Furthermore,
$\hat{d}^{\dagger}_{\sigma}$ and $\hat{d}_{\sigma}$ are those in the QD, and
$\hat{n}_{\sigma}=\hat{d}^{\dagger}_{\sigma}\hat{d}_{\sigma}$.
$U(\theta):=U_{0}\lambda(\theta)$ and $V_{\alpha}$ are, respectively, the
time-dependent electron-electron interaction in the QD and the transfer energy
between the QD and the reservoir $\alpha$. We adopt a model in the wide-band
limit for the reservoirs. We denote the line width in this paper as
$\Gamma=\pi\varrho(V_{\mathrm{L}}^{2}+V_{\mathrm{R}}^{2})$, where $\varrho$ is
the density of states in the reservoirs. For simplicity, we set $V_{L}=V_{R}$.

The Anderson model for the QD has four states: doubly occupied, singly
occupied with an up-spin, singly occupied with a down-spin, and empty.
Therefore, the density matrix is expressed as a $4\times 4$ matrix. However,
as shown in Ref. [44], $\hat{\rho}(\theta,\delta)$ of the Anderson model under
the wide-band approximation is reduced to a diagonal matrix, where the
diagonal elements correspond to the probability of finding the states in the
empty state $\rho_{e}$, the down-spin state $\rho_{\downarrow}$, the up-spin
state $\rho_{\uparrow}$, and the doubly occupied state $\rho_{d}$,
respectively. This implies that the model is quasi-classical.

The trace-preserving condition
$\mathrm{Tr}\hat{\rho}=\rho_{e}+\rho_{\uparrow}+\rho_{\downarrow}+\rho_{d}=1$
reduces to the probability preserving condition. The explicit forms of the
evolution matrix $\hat{K}$ and the corresponding eigenstates
$\langle\ell_{i}|$ and $|r_{i}\rangle$ are summarized in Ref. [44]. The
explicit forms of the BSN connection in Eq. (10) are presented in Appendix C.

Figure 4 shows the time evolution of the elements $\rho_{e}$,
$\rho_{\downarrow}$, $\rho_{\uparrow}$ and $\rho_{d}$ of $\hat{\rho}(\theta)$
(see also Appendix C). This figure clearly supports the positivity of all
elements, therefore the evolution dynamics satisfies CPTP properties. As shown
in Appendix D the approximate expression of the density matrix in Eqs. (5) and
(6) is very precise, where the deviation between the numerical solution of Eq.
(4) and the approximate expression in Eqs. (5) and (6) is almost invisible,

Figure 4: Scaled time ($\theta$) evolutions of the elements of the density
matrix. The parameters are set to be $\beta\bar{\mu}=0.1$, $\beta U_{0}=0.1$,
and $\beta\epsilon_{0}=0.1$.

In the following, we set $\epsilon=0.1$ to get some explicit results for the
Anderson model. Converting $|\hat{\rho}(\theta)\rangle$ to the matrix form
$\hat{\rho}(\theta,\delta)$ and inserting it into Eq. (II), we get
$D^{\mathrm{KL}}(\hat{\rho}(\theta,\delta)||\hat{\rho}^{\mathrm{GS}}(\theta,\delta))$.
The time evolution of
$D^{\mathrm{KL}}(\hat{\rho}(\theta,0)||\hat{\rho}^{\mathrm{GS}}(\theta,0))$ is
shown in Fig. 5 for $r=0.9$, and $\beta U_{0}=0.1$. This figure clearly
supports the H-theorem, where the KL divergence decreases monotonically as an
exponential function of $\theta$ to approach zero. Thus, we can evaluate
$\mathcal{Q}_{\mathrm{KL}}(\delta)$ introduced in Eq. (II). We have also
confirmed $\hat{\rho}(\theta=2\pi,\delta)=\hat{\rho}(4\pi,\delta)$ within the
numerical accuracy. This justifies the practical treatment of Eq. (8). As a
result, the trivial relation
$D^{\mathrm{KL}}(\hat{\rho}(\theta,\delta)||\hat{\rho}^{\mathrm{GS}}(\theta,\delta))=0$
holds for $\theta\geq 2\pi$.

Figure 5: Plot of
$D^{\mathrm{KL}}(\hat{\rho}(\theta,0)||\hat{\rho}^{\mathrm{GS}}(\theta,0))$
against $\theta$ for $\beta U_{0}=0.1$ and $r=0.9$. The inset is the semi-log
plot of
$D^{\mathrm{KL}}(\hat{\rho}(\theta,\delta)||\hat{\rho}^{\mathrm{GS}}(\theta,\delta))$
against $\theta$.  Figure 6: (a) Schematics of a contour of the integral of
$\mathcal{C}_{1}$, where the black solid line is the trajectory of the
parameters. The color scale at a value of $\theta$ expresses
$F^{{\mu_{\mathrm{L}}}{\mu_{\mathrm{R}}}}_{1}$. (b) and (c) The BSN curvatures
$F^{{\mu_{\mathrm{L}}}{\mu_{\mathrm{R}}}}_{1}$ and
$F^{{\mu_{\mathrm{L}}}{\mu_{\mathrm{R}}}}_{3}$ at $\theta=2\pi$ are plotted.
The parameters are set to be $\beta\bar{\mu}=0.1$, $\beta U_{0}=0.1$, and
$\beta\epsilon_{0}=0.1$ for all figures.

Figure 6(a) shows the contour of the integral given in Eq. (10) in the
parameter space $(\beta\mu_{\rm L}(\theta),\beta\mu_{\rm
R}(\theta,\delta),\theta)$. For the explicit calculation, we use Eqs. (5) and
(6) with $\theta=2\pi$ for Figs. 6 (b) and (c). As can be seen, the BSN
curvature is always present for $\theta>0$, although its magnitude decreases
with $\theta$. The BSN curvatures at specific $\theta$s are plotted in Figs.
6(b) and (c), where the peak and the dip of
$F^{\mu_{\mathrm{L}}\mu_{\mathrm{R}}}_{1}$ are approximately located at
$\beta\mu_{\mathrm{R}}=-\beta\mu_{\mathrm{L}}\approx 10$,
$\beta\mu_{\mathrm{R}}=-\beta\mu_{\mathrm{L}}\approx-10$, respectively, while
the peak and the dip of $F^{\mu_{\mathrm{L}}\mu_{\mathrm{R}}}_{3}$ are
approximately located at $\beta\mu_{\mathrm{L}}\approx\pm
10,\beta\mu_{\mathrm{R}}=0$ and $\beta\mu_{\mathrm{R}}\approx\pm
10,\beta\mu_{\mathrm{L}}=0$, respectively. Integrating by parts, we can verify
that $C_{2}(\theta,\delta)=0$ and
$F^{\mu_{\mathrm{L}}\mu_{\mathrm{R}}}_{2}=0$, because $\langle\ell_{2}|$ is
independent of $\theta$.

Figure 7 plots $D^{\mathrm{KL}}(\hat{\rho}(\theta,0)||\hat{\rho}^{\rm
SS}(\rho(\theta,0))$ versus $\theta$. This illustrates the rapid relaxation
near $\theta=0$ and the periodic oscillation of $\theta$. This non-monotonic
behavior means that $D^{\mathrm{KL}}(\hat{\rho}(\theta,0)||\hat{\rho}^{\rm
SS}(\rho(\theta,0))$ is not the proper KL divergence.

Figure 7: Plots of $D^{\mathrm{KL}}(\rho(\theta,0)|\hat{\rho}^{\rm
SS}(\hat{\rho}(\theta,0))$ versus $\theta$ for $r=0.9$ and $\beta U_{0}=0.3$.

Figure 8 plots $\Delta S(\delta)$ versus $\delta$ for various $\beta U_{0}$
for $r=0.9$. As can be seen in Eq. (II) and Fig. 7 the main contribution of
$\Delta S(\delta)$ comes from the drastic change near $\theta=0$ where the
mismatch between the geometric state and the initial steady state exists. We
also note that $\Delta S(\delta)$ approaches zero in the limit $\beta U_{0}\to
0$. This implies that the modulation of the Hamiltonian generates $\Delta
S(\delta)$. Although the periodic control of the Hamiltonian generates an
electric current and the corresponding Joule heat even though the two
reservoirs have no chemical potential difference, we have verified that their
contribution is small. Note that the expression of $\Delta S(\delta)$ does not
have room to include such contributions.

Figure 8: Plots of $\Delta S(\delta)$ versus $\delta$ for $\beta U_{0}=0.3$
(solid line), $0.5$ (dotted line), $0.7$ (dashed line) with fixing $r=0.9$.

We also plot the system entropy $S_{\rm sys}(\theta,0)$ against $\theta$ in
Fig. 9. As can be seen, $S_{\rm sys}(\theta,0)$ takes the maximum value at
$\theta\simeq\pi$. We also note that the behavior of
$S^{\mathrm{vN}}(\theta,0)$ is presented in Appendix D, which is qualitatively
different from
$D^{\mathrm{KL}}(\hat{\rho}(\theta,0)||\hat{\rho}^{\mathrm{GS}}(\theta,0))$,
$S_{\rm sys}(\theta,0)$ and ${\mathcal{Q}}(\theta,0)$.

Figure 9: Plot of $S_{\rm sys}(\theta,0)$ as a function of $\theta$ for $\beta
U_{0}=0.1$ and $r=0.9$.

We have confirmed that the work $W(\delta)$ introduced in Eq. (22) becomes
negative except for the small parameter region around $\delta=\pi$ (i.e.
$2.824\leq\delta\leq 3.569$), as shown in Fig. 10(a). This indicates that we
can extract the work by cyclic modulations of the parameters in the Anderson
model without fine-tuning in a wide range of parameters. The change of the
energy in the first cycle $\Delta E(\delta)$ is shown in Fig. 10(b). We have
verified the relationship $|\Delta E(\delta)|\ll|W(\delta)|$ and thus we have
$|W(\delta)|\simeq|Q(\delta)|$ in this set of parameters. The heat absorbed
during a modulation cycle is shown in Fig. 10(c). Figure 10(d) shows the
efficiency $\eta(\delta)$, defined in Eq. (33), for one modulation cycle. The
efficiency is not defined for $2.824\leq\delta\leq 3.569$ because the work
becomes positive in this region. These negative work and efficiency are the
result of $\Delta S(\delta)<0$. Thus, our engine is suitable for calling the
geometric demon.

Figure 10:  Plots of the work, the change of the internal energy, the
absorbing heat, and the efficiency against $\delta$ for $\beta U_{0}=0.1$ and
$r=0.9$. (a) Plot of $\beta W$ against $\delta$ (solid line). The region where
the geometric demon is not achieved is shown as the dashed line. (b) Plot of
$\beta\Delta E(\delta)$ against $\delta$ (solid line). (c) Plot of
$\beta\mathcal{Q}_{A}$ against $\delta$ (solid line). (d) Plot of $\eta$
against $\delta$ (solid line). The region $2.824\leq\delta\leq 3.569$ is kept
blank since the work becomes positive in the region.

##  V Concluding Remarks

We have implemented a geometric demon by modulating the electrochemical
potentials in the two reservoirs and the repulsive $U(\theta)$ in the system
Hamiltonian under isothermal conditions. We found that the state tends to
relax towards the geometric state instead of the steady state when the
parameter is kept modulated. The discrepancy between the steady state and the
geometric state under a non-adiabatic modulation leads to the negative
production of relative entropy. We can automatically extract the work from
this engine with an increase in relative entropy if we start from the
nonequilibrium steady state. We have verified that the work becomes negative
in the wide range of the parameters for the Anderson model and thus the
geometric demon can be implemented. Our geometric demon does not require any
observation of states to reduce entropy. In this sense, our geometric demon
can be easily implemented in realistic situations; therefore, we expect wide
applications of this demon. Note that efficiency is unchanged even if we
consider the second cycle. Indeed, we can stop the modulation after the first
cycle and wait for the realization of the steady state, and we restart the
modulation.

Our future tasks are as follows. (i) Since the present method of the argument
is restricted to the case of $\epsilon\ll 1$, we will need to extend the
analysis to the regime of larger $\epsilon$ as done in Ref. [36]. (ii)
Although we have analyzed a pure quantum system to clarify entanglement or
coherence effects [41, 55, 56], our treatment remains quasi-classical.

Acknowledgements.- The authors thank Ville M. M. Paasonen, Kazutaka Takahashi,
Kiyoshi Kanazawa, Sosuke Ito, Hikaru Watanabe, Asahi Yamaguchi and Frédéric
van Wijland for fruitful discussions and useful comments. This study was
partially supported by a Grant-in-Aid of MEXT for Scientific Research KAKENHI
(Grant Nos. 21H01006, 19K14616, and 20H01838).

##  Appendix A Time evolution of the density matrix

This section consists of two subsections. In the first part, A.1, we derive
the time evolution of $\hat{\rho}(\theta,\delta)$ from the general initial
condition to demonstrate the universality of Eqs. (5) and (6) for
$\theta\gg\epsilon$. In the second part, A.2, we present the detailed
derivations of the BSN connection and BSN curvature.

###  A.1 Derivation of the time-dependent expression of the density matrix
starting from the general initial condition

####  A.1.1 Time evolution and geometrical state

The purpose of this section is to derive Eqs. (5) and (6). Although we assume
in the main text that the initial state is given by $\hat{\rho}^{\rm
SS}(\theta=0,\delta)$, we can still derive these equations even when if we
start with the generalized initial condition:

$|\hat{\rho}_{\mathrm{ini}}\rangle=\sum_{i}a_{i}|r_{i}(0,\delta)\rangle,$ (38)

where $a_{i}$ is given by

$a_{i}:=\langle\ell_{i}|\hat{\rho}_{\mathrm{ini}}\rangle.$ (39)

The normalization of the density matrix leads to the condition:

$\mathrm{Tr}\hat{\rho}_{\mathrm{ini}}=\langle\ell_{0}|\hat{\rho}_{\mathrm{ini}}\rangle=\sum_{i}a_{i}\langle\ell_{0}|r_{i}(0,\delta)\rangle=a_{0}=1.$
(40)

The time evolution of $|\rho(\theta,\delta)\rangle$ under Eq. (4) is given
formally as

$\displaystyle|\hat{\rho}(\theta,\delta)\rangle=\mathcal{T}\exp\left(\int^{\theta}_{0}d\phi\hat{K}(\phi,\delta)\right)|\hat{\rho}_{\mathrm{ini}}\rangle,$
(41)

As shown in the main text, the formal expression Eq. (41) quickly approaches
the geometric state $|\hat{\rho}^{\rm{GS}}(\theta,\delta)\rangle$. Although
the geometric state is defined in the limit $\theta\to\infty$ in the strict
sense, the difference between Eq. (41) and the geometric state
$|\hat{\rho}^{\rm{GS}}(\theta,\delta)\rangle$ decreases exponentially as
$\theta$ increases. We note that Eq. (41) satisfies the normalization
condition

$\displaystyle\mathrm{Tr}\hat{\rho}(\theta,\delta)=\langle\ell_{0}|\hat{\rho}(\theta,\delta)\rangle=\langle\ell_{0}|\hat{\rho}_{\mathrm{ini}}\rangle=1,$
(42)

because $\langle\ell_{0}|$ is the left zero-eigenstate of
$\hat{K}(\phi,\delta)$. Before entering the detailed calculation, we briefly
explain the late time behavior of Eq. (41). As we will see below, except for
$i=0$ the term in Eq. (38) is exponentially damped, and since $a_{0}=1$ the
late time behavior hardly depends on the choice of the initial state.

####  A.1.2 Path-integral representation of the state

Let us use the pass-integral-like expression for
$|\hat{\rho}(\theta,\delta)\rangle$ as

$\displaystyle|\hat{\rho}(\theta+\Delta\theta,\delta)\rangle$
$\displaystyle=e^{\epsilon^{-1}\hat{K}(\theta)\Delta\theta}|\hat{\rho}(\theta,\delta)\rangle$
$\displaystyle=e^{\epsilon^{-1}\hat{K}(\theta)\Delta\theta}\left(\sum_{m}|r_{m}(\theta,\delta)\rangle\langle\ell_{m}(\theta,\delta)|\right)|\hat{\rho}(\theta,\delta)\rangle$
$\displaystyle=\sum_{m}|r_{m}(\theta,\delta)\rangle
e^{\varepsilon_{m}(\theta)\frac{\Delta\theta}{\epsilon}}\langle\ell_{m}(\theta,\delta)|\hat{\rho}(\theta,\delta)\rangle.$
(43)

Repeating Eq. (43) we can rewrite
$|\hat{\rho}(\theta+\Delta\theta,\delta)\rangle$ as

$\displaystyle|\hat{\rho}(\theta+\Delta\theta,\delta)\rangle$
$\displaystyle=\sum_{m,n}|r_{m}(\theta,\delta)\rangle
e^{\varepsilon_{m}(\theta)\frac{\Delta\theta}{\epsilon}}\langle\ell_{m}(\theta,\delta)|r_{m}(\theta-\Delta\theta,\delta)\rangle$
$\displaystyle\hskip
30.00005pte^{\varepsilon_{n}(\theta-\Delta\theta)\frac{\Delta\theta}{\epsilon}}\langle\ell_{m}(\theta-\Delta\theta,\delta)|\hat{\rho}(\theta-\Delta\theta,\delta)\rangle.$
(44)

Thus, the recursive calculation yields

$\displaystyle|\hat{\rho}(\theta+\Delta\theta,\delta)\rangle$
$\displaystyle=\sum_{i,j,k,\cdots,l}|r_{i}(\theta,\delta)\rangle
e^{\varepsilon_{i}(\theta,\delta)\frac{\Delta\theta}{\epsilon}}\langle\ell_{i}(\theta,\delta)|r_{j}(\theta-\Delta\theta,\delta)\rangle$
$\displaystyle\hskip 30.00005pt\times
e^{\varepsilon_{j}(\theta-\Delta\theta,\delta)\frac{\Delta\theta}{\epsilon}}\langle\ell_{j}(\theta-\Delta\theta,\delta)|r_{k}(\theta-2\Delta\theta,\delta)\rangle$
$\displaystyle\hskip 30.00005pt\cdots\times
e^{\varepsilon_{l}(\Delta\theta)\frac{\Delta\theta}{\epsilon}}\langle\ell_{l}(\Delta\theta,\delta)|\hat{\rho}(0,\delta)\rangle$
$\displaystyle=\sum_{i,j,k,\cdots,l,m}a_{n}|r_{i}(\theta,\delta)\rangle
e^{\varepsilon_{i}(\theta,\delta)\frac{\Delta\theta}{\epsilon}}\langle\ell_{i}(\theta,\delta)|r_{j}(\theta-\Delta\theta,\delta)\rangle$
$\displaystyle\hskip 30.00005pt\times
e^{\varepsilon_{j}(\theta-\Delta\theta,\delta)\frac{\Delta\theta}{\epsilon}}\langle\ell_{j}(\theta-\Delta\theta,\delta)|r_{k}(\theta-2\Delta\theta,\delta)\rangle$
$\displaystyle\hskip 30.00005pt\cdots\times
e^{\varepsilon_{l}(\Delta\theta)\frac{\Delta\theta}{\epsilon}}\langle\ell_{l}(\Delta\theta,\delta)|r_{n}(0,\delta)\rangle,$
(45)

where we have adopted the general initial condition
$|\hat{\rho}(0,\delta)\rangle=|\hat{\rho}_{\mathrm{ini}}\rangle=\sum_{n}a_{n}|r_{n}(0,\delta)\rangle$.

Next we evaluate the connection (transfer matrix)
$\langle\ell_{i}(\theta,\delta)|r_{j}(\theta-\Delta\theta,\delta)\rangle$ that
connects the states in the different $\theta$. The connection can be
calculated as follows up to the first order of $\Delta\theta$.

$\displaystyle\langle\ell_{i}(\theta,\delta)|r_{j}(\theta-\Delta\theta,\delta)\rangle$
$\displaystyle=\langle\ell_{i}(\theta,\delta)|\left(|r_{j}(\theta,\delta)\rangle-\Delta\theta\frac{d}{d\theta}|r_{j}(\theta,\delta)\rangle\right)$
$\displaystyle=\delta_{ij}-\Delta\theta\langle\ell_{i}(\theta,\delta)|\frac{d}{d\theta}|r_{j}(\theta,\delta)\rangle.$
(46)

The two terms in the last expression of Eq. (46) have a natural
interpretation; the first describing the second free propagation, and the
second describing the jump process from the state $j$ to the state $i$. If we
substitute Eq. (46) into Eq. (45), we obtain the expression:

$\displaystyle|\hat{\rho}(\theta,\delta)\rangle=|\hat{\rho}(\theta,\delta)\rangle_{0}+|\hat{\rho}(\theta,\delta)\rangle_{1}+|\hat{\rho}(\theta,\delta)\rangle_{2}+\cdots.$
(47)

The first term on the RHS of Eq. (47) is the free relaxation:

$\displaystyle|\hat{\rho}(\theta,\delta)\rangle_{0}$
$\displaystyle=\sum_{n}a_{n}|r_{n}(\theta,\delta)\rangle K_{n}(\theta,0),$
(48)

where we have introduced $K_{n}(b,a)$ as

$\displaystyle
K_{n}(b,a):=e^{\int_{a}^{b}\varepsilon_{n}(\phi,\delta)\frac{d\phi}{\epsilon}}.$
(49)

The second term on the RHS of Eq. (47) contains a jump process:

$\displaystyle|\hat{\rho}(\theta,\delta)\rangle_{1}=\sum_{m,n}C_{mn}(\theta,\delta)a_{n}|r_{m}(\theta,\delta)\rangle,$
(50)

where the coefficient $C_{mn}$ is given by

$\displaystyle C_{mn}(\theta,\delta):=\int_{0}^{\theta}d\phi
K_{m}(\theta,\phi)J_{mn}(\phi)K_{n}(\phi,0)$ (51)

with the jump operator

$\displaystyle
J_{mn}(\phi):=-\langle\ell_{m}(\phi,\delta)|\frac{d}{d\phi}|r_{n}(\phi,\delta)\rangle.$
(52)

The first and the second terms on the RHS of Eq. (47) can be understood by the
diagrams shown in Figs. 11(a) (I) and (II), where the straight line represents
the free relaxation and the cross mark represents the jump process. The first
term on the RHS of Eq. (47) represents the free relaxation without any jump
process. The second term on the RHS of Eq. (47), the initial state
$\sum_{n}a_{n}|r_{n}(0)\rangle$ evolves with $\theta$ by $K_{n}(\phi,0)$ and
then the state $|r_{n}\rangle$ jumps to the state $|r_{m}\rangle$ at $\phi$
and then it evolves with $\theta$ by $K_{m}(\theta,\phi)$. Using Fig.
11(a)(III), we can write the third term on the RHS of Eq. (47)

$\displaystyle|\hat{\rho}(\theta,\delta)\rangle_{2}=\sum_{m,n,k}C_{kmn}(\theta)a_{n}|r_{k}(\theta,\delta)\rangle,$
(53)

where

$\displaystyle C_{kmn}(\theta,\delta)$
$\displaystyle=\int_{0}^{\theta}d\phi\int_{\phi}^{\theta}d\phi^{\prime}K_{k}(\theta,\phi^{\prime})J_{km}(\phi^{\prime})$
$\displaystyle\hskip 70.0001pt\times
K_{m}(\phi^{\prime},\phi)J_{mn}(\phi)K_{n}(\phi,0).$ (54)

Figure 11: Diagrammatic expressions of the Born type expansion. Figure 11 (a)
shows the 0th order (I), 1st order (II), and 2nd order (III) processes,
respectively. The state evolves with $\theta$ by the free relaxation process
(straight line) until the scattering process (represented by the cross mark)
takes place (Fig. 11 (b)).

####  A.1.3 Born-like approximation and further simplification

In the following, we assume that the modulation speed is slow enough to ignore
higher-order jump processes. Thus, we consider only the first and the second
terms on the RHS of Eq. (47). This approximation is analogous to the Born
approximation in the scattering theory. As emphasized in the main text, the
state vector quickly approaches the geometric state with the Born-like
approximation which is independent of the initial condition but depends on the
BSN connection term.

In the Born-like approximation, the resulting state is

$\displaystyle|\hat{\rho}(\theta,\delta)\rangle_{0}+|\hat{\rho}(\theta,\delta)\rangle_{1}$
$\displaystyle=\sum_{m}a_{m}e^{\int_{0}^{\theta}\varepsilon_{m}(\phi,\delta)\frac{d\phi}{\epsilon}}|r_{m}(\theta,\delta)\rangle+\sum_{m,n}C_{mn}a_{n}|r_{m}(\theta,\delta)\rangle.$
(55)

The phase factor in the first term on RHS of Eq. (55) is the dynamical phase,
while the second term on the RHS is the term generated by the geometric phase.

Since $a_{0}=1$ and $\varepsilon_{0}=0$ we can write Eq. (55) as

$\displaystyle|\hat{\rho}(\theta,\delta)\rangle_{0}+|\hat{\rho}(\theta,\delta)\rangle_{1}$
$\displaystyle=|r_{0}(\theta,\delta)\rangle+\sum_{m\neq
0}a_{m}e^{\int_{0}^{\theta}\varepsilon_{m}(\phi,\delta)\frac{d\phi}{\epsilon}}|r_{m}(\theta,\delta)\rangle$
$\displaystyle\quad+\sum_{m}C_{m0}|r_{m}(\theta,\delta)\rangle+\sum_{m\neq
0,n\neq 0}C_{mn}a_{n}|r_{m}(\theta,\delta)\rangle$
$\displaystyle\quad+\sum_{n}C_{0n}a_{n}|r_{m}(\theta,\delta)\rangle.$ (56)

The coefficient in the last term on the RHS of Eq. 56 is zero as

$\displaystyle C_{0n}=$ $\displaystyle-\int^{\theta}_{0}d\phi
e^{\int_{0}^{\phi}\varepsilon_{n}(\zeta,\delta)\frac{d\zeta}{\epsilon}}\langle\ell_{0}(\phi,\delta)|\frac{d}{d\phi}|r_{n}(\phi,\delta)\rangle$
$\displaystyle=$
$\displaystyle-e^{\int_{0}^{\phi}\varepsilon_{n}(\zeta,\delta)\frac{d\zeta}{\epsilon}}\left\\{\langle\ell_{0}|r_{n}(\theta,\delta)\rangle-\langle\ell_{0}|r_{n}(0,\delta)\rangle\right\\}$
$\displaystyle+\int^{\theta}_{0}d\phi\frac{d}{d\phi}\left(e^{\int_{0}^{\phi}\varepsilon_{n}(\zeta,\delta)\frac{d\zeta}{\epsilon}}\langle\ell_{0}|\right)|r_{n}(\phi,\delta)\rangle$
$\displaystyle=$
$\displaystyle\delta_{0n}\int^{\theta}_{0}d\phi\frac{d}{d\phi}\left(e^{\int_{0}^{\phi}\varepsilon_{n}(\zeta,\delta)\frac{d\zeta}{\epsilon}}\right)=0,$
(57)

where we have used that $\langle\ell_{0}(\phi,\delta)|$ is independent of
$\phi$ and $\delta$, and thus
$\frac{d}{d\phi}\langle\ell_{0}(\phi,\delta)|=0$. We have also used the
orthonormal relation $\langle\ell_{m}|r_{n}\rangle=\delta_{mn}$. The last
equality in Eq. (A.1.4) holds since $\delta_{0n}=0$ for $n\neq 0$ and
$\frac{d}{d\phi}\left(e^{\int_{0}^{\phi}\varepsilon_{n}(\zeta,\delta)\frac{d\zeta}{\epsilon}}\right)=0$
for $n=0$. Thus, Eq. (56) becomes

$\displaystyle|\hat{\rho}(\theta,\delta)\rangle_{0}+|\hat{\rho}(\theta,\delta)\rangle_{1}$
$\displaystyle=|r_{0}(\theta,\delta)\rangle+\sum_{m\neq
0}\left[a_{m}e^{\int_{0}^{\theta}\varepsilon_{m}(\phi,\delta)\frac{d\phi}{\epsilon}}+C_{m0}\right]|r_{m}(\theta,\delta)\rangle$
$\displaystyle\quad+\sum_{m\neq 0,n\neq
0}C_{mn}a_{n}|r_{m}(\theta,\delta)\rangle$ (58)

Note that the coefficient $C_{m0}$ is equivalent to $C_{m}$ in the main text.
We also note that we can directly obtain the relation
$\langle\ell_{0}|\hat{\rho}(\theta,\delta)\rangle_{n\geq 1}=0$ from trace
preservation

####  A.1.4 Late time behavior

Now, let us consider the behavior for $\theta/\epsilon\gg 1$. It is obvious
that the second term on the RHS of Eq. (58) is negligible, because of the
exponential damping factor that appears for $m\neq 0$. We also note that the
sum of the first and third terms on the RHS is equivalent to
$|\hat{\rho}(\theta,\delta)\rangle$ used in Eq. (5). Thus, Eq. (58) is reduced
to

$\displaystyle|\hat{\rho}(\theta,\delta)\rangle_{0}+|\hat{\rho}(\theta,\delta)\rangle_{1}$
$\displaystyle=|\hat{\rho}^{\rm
GS}(\theta,\delta)\rangle+\sum_{m,n}C_{mn}(\theta,\delta)a_{n}|r_{m}(\theta,\delta)\rangle$
(59)

for $\theta/\epsilon\gg 1$. It is straightforward to evaluate
$C_{mn}(\theta,\delta)$ for $n\neq 0$ as

$\displaystyle|C_{mn}(\theta)|=$ $\displaystyle\left|\int_{0}^{\theta}d\phi
e^{\int_{\phi}^{\theta}\varepsilon_{m}(\xi)\frac{d\xi}{\epsilon}}e^{\int_{0}^{\phi}\varepsilon_{n}(\zeta)\frac{d\zeta}{\epsilon}}\langle\ell_{m}(\phi)|\frac{d}{d\phi}|r_{n}(\phi)\rangle\right|$
$\displaystyle\leq$ $\displaystyle\int_{0}^{\theta}d\phi
e^{\int_{\phi}^{\theta}\varepsilon_{m}(\xi)\frac{d\xi}{\epsilon}}e^{\int_{0}^{\phi}\varepsilon_{n}(\zeta)\frac{d\zeta}{\epsilon}}\left|\langle\ell_{m}(\phi)|\frac{d}{d\phi}|r_{n}(\phi)\rangle\right|$
$\displaystyle\leq$ $\displaystyle\int_{0}^{\theta}d\phi
e^{\int_{0}^{\theta}\mathrm{max}(\varepsilon_{m},\varepsilon_{n})(\xi)\frac{d\xi}{\epsilon}}\left|\langle\ell_{m}(\phi)|\frac{d}{d\phi}|r_{n}(\phi)\rangle\right|,$
(60)

where we have omitted the $\delta$ dependence in Eq. (60). Thus, $C_{mn}$ with
$n\neq 0$ is much smaller than $C_{m0}$ because of the existence of the
exponential factor in Eq. (60). It can be shown that this late time behavior
does not depend on the initial state since the coefficient
$C_{m0}(\theta,\delta)$ depends only on the steady state
$|r_{0}(\phi,\delta)\rangle$ and the states $|r_{m}(\phi,\delta)\rangle$ for
$|\phi-\theta|\ll\theta_{0}$, where $\theta_{0}$ is the largest relaxation
time of the states $|r_{m}(\phi,\delta)\rangle$. It is easy to show that the
higher Born corrections are also independent of the initial state. Thus the
geometric state can be approximated as

$\displaystyle|\hat{\rho}^{\rm GS}(\theta,\delta)\rangle$
$\displaystyle\simeq|\hat{\rho}^{\rm
SS}(\theta,\delta)\rangle+\sum_{m}C_{m0}(\theta,\delta)|r_{m}(\theta,\delta)\rangle\
(\theta\gg\theta_{0}),$ (61)

where the equality holds up to the negligibly small deviation as in the case
of Eq. (8). Here we note that (61) hardly depends on the choice of the initial
state since

$\displaystyle C_{m0}=$ $\displaystyle-\int^{\theta}_{0}d\phi
e^{\int_{\phi}^{\theta}\varepsilon_{m}(\zeta,\delta)\frac{d\zeta}{\epsilon}}\langle\ell_{m}(\phi,\delta)|\frac{d}{d\phi}|r_{0}(\phi,\delta)\rangle$
$\displaystyle\simeq$ $\displaystyle-\int^{\theta}_{\theta_{0}}d\phi
e^{\int_{\phi}^{\theta}\varepsilon_{m}(\zeta,\delta)\frac{d\zeta}{\epsilon}}\langle\ell_{m}(\phi,\delta)|\frac{d}{d\phi}|r_{0}(\phi,\delta)\rangle,$
(62)

where the difference between the RHS of the first line and that of the second
line is negligibly small due to the exponential damping factor.

###  A.2 Derivation of BSN connection and BSN curvature

For the cyclic modulation
$|r_{i}(2\pi,\delta)\rangle=|r_{i}(0,\delta)\rangle$, the time evolution from
the general initial state becomes

$\displaystyle\Delta|\hat{\rho}\rangle=$ $\displaystyle\sum_{i\neq
0}a_{i}\left(e^{\int_{0}^{2\pi}\varepsilon_{i}(\phi,\delta)\frac{d\phi}{\epsilon}}-1\right)|r_{i}(0,\delta)\rangle$
$\displaystyle+\sum_{i,j}\mathcal{C}_{ij}a_{j}e^{\int_{0}^{2\pi}\varepsilon_{i}(\phi,\delta)\frac{d\phi}{\epsilon}}|r_{i}(0,\delta)\rangle,$
(63)

where

$\displaystyle\mathcal{C}_{ij}:=$
$\displaystyle\delta_{ij}\oint_{\partial\Omega}d\phi\frac{\varepsilon_{i}(\phi,\delta)}{\epsilon}-\oint_{\partial\Omega}d\Lambda_{\mu}\mathcal{A}^{\mu}_{ij}$
(64)

with the introduction of the BSN connection $\mathcal{A}^{\mu}_{ij}$:

$\displaystyle\mathcal{A}^{\mu}_{ij}:=\langle\tilde{\ell}_{i}(\phi,\delta)|\frac{\partial}{\partial\Lambda_{\mu}}|\tilde{r}_{j}(\phi,\delta)\rangle,$
(65)

where we have introduced the gauged vectors

$\displaystyle\langle\tilde{\ell}_{i}(\phi,\delta)|:=\langle\ell_{i}(\phi,\delta)|e^{-\epsilon^{-1}\int_{0}^{\phi}\varepsilon_{i}(\zeta,\delta)\frac{d\zeta}{\epsilon}},$
(66)
$\displaystyle|\tilde{r}_{j}(\phi,\delta)\rangle:=e^{\int_{0}^{\phi}\varepsilon_{m}(\zeta,\delta)\frac{d\zeta}{\epsilon}}|\tilde{r}_{j}(\phi,\delta)\rangle.$
(67)

Here, the sum of the first term on the RHS of Eq. (63) is taken except for
$i=0$ because $\varepsilon_{0}(\phi,\delta)=0$ always gives
$e^{\int_{0}^{\theta}\varepsilon_{0}(\phi,\delta)\frac{d\phi}{\epsilon}}=1$.
Using the Stokes theorem, we can rewrite Eq. (64) as

$\displaystyle\mathcal{C}_{ij}=\delta_{ij}\oint_{\partial\Omega}d\phi\frac{\varepsilon_{i}(\phi,\delta)}{\epsilon}-\int_{\Omega}dS_{\mu\nu}F^{\mu\nu}_{ij},$
(68)

where $\Omega$ is the area enclosed by the closed trajectory $\partial\Omega$,
$S_{\mu\nu}=\frac{1}{2}d\Lambda_{\mu}\wedge d\Lambda_{\nu}$, and
$F^{\mu\nu}_{ij}$ represents the BSN curvatures defined as

$\displaystyle
F^{\mu\nu}_{ij}:=\frac{\partial\langle\tilde{\ell}_{i}|}{\partial\Lambda_{\nu}}\frac{\partial|\tilde{r}_{j}\rangle}{\partial\Lambda_{\mu}}-\frac{\partial\langle\tilde{\ell}_{i}|}{\partial\Lambda_{\mu}}\frac{\partial|\tilde{r}_{j}\rangle}{\partial\Lambda_{\nu}}.$
(69)

It is also possible to rewrite Eq. (68) as

$\displaystyle\mathcal{C}_{ij}$
$\displaystyle=\delta_{ij}\oint_{\partial\Omega}d\phi\frac{\varepsilon_{i}(\phi,\delta)}{\epsilon}+\frac{1}{2}\int_{\Omega}d\langle\tilde{\ell}_{i}|\wedge
d|\tilde{r}_{j}\rangle.$ (70)

##  Appendix B Differences between the present study and the previous studies

This section is devoted to explaining the difference between the present study
and the previous studies. In the first part, we show the consistency between
the results in Ref. [36] and those in the present study. In the second part,
we present the difference between the setup in the present study and the
previous study [44].

###  B.1 The relationship of our work with Ref. [36]

Here, we briefly comment on the superficial differences between the present
study and Ref. [36]. In the present analysis, we have fixed the gauge degrees
of freedom to fix the overall scale of states. It is known that
$\\{|r_{m}(\theta)\rangle,\langle\ell_{m}(\theta)|\\}$ and
$\\{e^{-f_{m}(\theta)}|r_{m}(\theta)\rangle,e^{f_{m}(\theta)}\langle\ell_{m}(\theta)|\\}$
yield the same results. The gauge freedom can be fixed if we choose the
scaling, e.g. $\langle\ell_{0}|=(1,1,1,1)$ for the Anderson model. We also
note that the 0th-order Born approximation in this paper is called the
adiabatic approximation in Ref. [36].

In Ref. [36], the gauge degrees of freedom for the left eigenstates
($\langle\phi_{n}|$) and the right eigenstates ($|\phi_{m}\rangle$) are fixed
by 333We note that the lowest eigenstate is labeled by $m=1$ in Ref. [36]
whereas it is labeled by $m=0$ in this paper. To avoid confusion, we relabel
the subscripts $m$ for the eigenstates in Ref. [36] from $\\{1,\cdots\\}$ to
$\\{0,\cdots\\}$

$\displaystyle|\tilde{\phi}_{m}(\theta)\rangle=e^{-\int^{\theta}_{0}d\eta\langle\phi_{m}(\eta)|\frac{d}{d\eta}|\phi_{m}(\eta)\rangle}|\phi_{m}(\theta)\rangle,$
(71)
$\displaystyle\langle\tilde{\phi}_{m}(\theta)|=e^{\int^{\theta}_{0}d\eta\langle\phi_{m}(\eta)|\frac{d}{d\eta}|\phi_{m}(\eta)\rangle}\langle\phi_{m}(\theta)|.$
(72)

Since both states $|\phi_{m}\rangle,\langle\phi_{m}|$ and
$|r_{m}\rangle,\langle\ell_{m}|$ are left and right eigenstates of $\hat{K}$
with eigenvalue $\varepsilon_{m}$, those are identical up to the scaling
factor as

$\displaystyle|\phi_{m}(\theta)\rangle=e^{g(\theta)}|r_{m}(\theta)\rangle,$
(73)
$\displaystyle\langle\phi_{m}(\theta)|=e^{-g(\theta)}\langle\ell_{m}(\theta)|.$
(74)

Substituting Eqs. (73) and (74) into Eqs. (71) and (72), we obtain the
relation between
$\\{|\tilde{\phi}_{m}(\theta)\rangle,\langle\tilde{\phi}_{m}(\theta)\\}$ in
Ref. [36] and $\\{|r_{m}(\theta)\rangle,\langle\ell_{m}(\theta)\\}$ in the
present study as

$\displaystyle|\tilde{\phi}_{m}(\theta)\rangle=e^{g(0)+\int_{0}^{\theta}d\eta\langle\ell_{m}(\eta)|\frac{d}{d\eta}|r_{m}(\eta)\rangle}|r_{m}(\theta)\rangle,$
(75)
$\displaystyle\langle\tilde{\phi}_{m}(\theta)|=e^{-g(0)-\int_{0}^{\theta}\langle\ell_{m}(\eta)|\frac{d}{d\eta}|r_{m}(\eta)\rangle}\langle\ell_{m}(\theta)|.$
(76)

By using
$\int_{0}^{\theta}d\eta\langle\ell_{0}(\eta)|\frac{d}{d\eta}|r_{0}(\eta)\rangle=\left[\langle\ell_{0}(\eta)|r_{0}(\eta)\rangle\right]^{\theta}_{0}-\int_{0}^{\theta}\frac{d\langle\ell_{0}(\eta)|}{d\eta}|r_{0}(\eta)\rangle=0$,
one can show that the lowest eigenstate in Ref. [36] is identical to
$|r_{0}\rangle,\langle\ell_{0}|$.

We note that if we ignore the last term on the RHS of Eq. (58), our expression
is reduced to Eq. (5) in Ref. [36], where they assume that the time evolution
does not make transitions to different eigenstates, which corresponds to the
Born-type approximation in our approach. The other assumption is made in Ref.
[36] which corresponds to setting $C_{mn}=0$ with $n\neq 0$ in Eq. (58). This
approximation is also justified in our approach and will be explained in the
next subsection. With these approximations, we can directly show the
consistency between our results and those in Ref. [36] by substituting Eqs.
(75) and (76) into Eq. (58). In the main part of our paper, since we start
from the steady state ($a_{0}=1,a_{m\neq 0}$), $C_{mn}$ with $n\neq 0$ does
not appear, and thus the latter assumption is not necessary,

###  B.2 Difference between Ref. [44] and the present study

A major difference between this paper and Ref. [44] is the setup. We mainly
discuss the relaxation process from NESS in this paper, but Ref. [44]
discussed the behavior after the system reaches a quasi-periodic state, in
which the geometric term can be ignored. Thus the previous study did not
consider the initial relaxation and treated the steady state $\rho^{\rm SS}$.

It is obvious that $\hat{\rho}^{\mathrm{GS}}(\theta,\delta)$ is different from
$\hat{\rho}^{\mathrm{SS}}(\theta,\delta)$ due to the existence of the BSN
connection term as shown in Eqs. (5) and (6). Moreover, as shown in the main
text, the mismatch between $\hat{\rho}^{\mathrm{SS}}(0,\delta)$ and
$\hat{\rho}^{\mathrm{GS}}(\theta,\delta)$ for small $\theta$ plays a
significant role. Thus, if we are interested in the initial relaxation process
from $\hat{\rho}^{\mathrm{SS}}(\delta)$ to
$\hat{\rho}^{\mathrm{GS}}(\theta,\delta)$, we cannot use the $\epsilon$
perturbation around $\hat{\rho}^{\mathrm{SS}}(\theta,\delta)$ used in Ref.
[44].

##  Appendix C Properties of the Anderson model

This section summarizes the properties of the Anderson model in more detail,
which consists of three subsections. The first subsection summarizes the
evolution matrix. In the second subsection, we provide the explicit forms of
the BSN connection and BSN curvature as well as the expansion coefficients. In
the third subsection, we present some detailed calculations for the Anderson
model. For simplicity, we do not write the $\delta$ dependence of the
variables explicitly in this section.

###  C.1 Evolution matrix and eigenstates in the Anderson model

In this subsection, we summarize the evolution matrix and eigenstates in the
Anderson model. A similar discussion can be found in Ref. [44].

Since $\hat{\rho}$ is a diagonal matrix, $|\hat{\rho}\rangle$ also has only
four components, and the transition matrix $\hat{K}$ in Eq. (4) in the
wideband approximation is given by the $4\times 4$ matrix:

$\hat{K}=-\begin{pmatrix}2f_{-}^{(1)}&-f_{+}^{(1)}&-f_{+}^{(1)}&0\\\
-f_{-}^{(1)}&f_{-}^{(0)}+f_{+}^{(1)}&0&-f_{+}^{(0)}\\\
-f_{-}^{(1)}&0&f_{-}^{(0)}+f_{+}^{(1)}&-f_{+}^{(0)}\\\
0&-f_{-}^{(0)}&-f_{-}^{(0)}&2f_{+}^{(0)}\\\ \end{pmatrix},$ (77)

where we have introduced

$\displaystyle f_{+}^{(j)}:=f_{\mathrm{L}}^{(j)}+f_{\mathrm{R}}^{(j)},\
f_{-}^{(j)}:=2-f_{+}^{(j)},$ (78)

with the Fermi distribution

$f_{\alpha}^{(j)}(\mu_{\alpha}(\theta),U(\theta)):=\frac{1}{1+e^{\beta(\epsilon_{0}+jU(\theta)-\mu_{\alpha}(\theta))}}$
(79)

in the lead $\alpha(=L\ {\rm or}\ R)$ for the single occupancy $j=0$ and
double occupancy $j=1$. The eigenvalues and eigenvectors correspond to the
evolution matrix $\hat{K}(\theta,\delta)$ in Eq. (77) are presented in Ref.
[44].

###  C.2 BSN connection for Anderson model

In this subsection, we present an explicit form of the BSN connection for the
Anderson model. Note that the calculation presented here is derived directly
from the detailed expressions of the eigenfunctions and eigenvalues of Eq. (4)
(see Ref. [44] for details).

For the explicit calculation of the BSN connection, we use Eqs. (5) and (6)
with $\theta=2\pi$. Integrating by parts, we get

$\displaystyle\int_{0}^{2\pi}e^{\int_{\phi}^{2\pi}\varepsilon_{2}(\xi)\frac{d\xi}{\epsilon}}\langle\ell_{2}(\phi)|\frac{d}{d\phi}|r_{0}(\phi)\rangle$
$\displaystyle=\int_{0}^{2\pi}e^{\int_{\phi}^{2\pi}\varepsilon_{2}(\xi)\frac{d\xi}{\epsilon}}\left[\frac{d}{d\phi}\langle\ell_{2}(\phi)|r_{0}(\phi)\rangle-\frac{d\langle\ell_{2}(\phi)|}{d\phi}|r_{0}(\phi)\rangle\right]$
$\displaystyle=0,$ (80)

where $\langle\ell_{2}(\phi)|$ is independent of $\phi$. Thus, the summation
of $i$ in Eq. (5) is reduced to the summation with $i=1$ and 3. The
differentiation of $|r_{0}(\phi)\rangle$ with respect to $\phi$ becomes

$\displaystyle\frac{d}{d\phi}|r_{0}(\phi)\rangle$
$\displaystyle=\left[\frac{d}{d\phi}\frac{1}{2(f_{+}^{(0)}+f_{-}^{(1)})}\right]\begin{pmatrix}f_{+}^{(0)}f_{+}^{(1)}\\\\[5.0pt]
f_{+}^{(0)}f_{-}^{(1)}\\\\[5.0pt] f_{+}^{(0)}f_{-}^{(1)}\\\\[5.0pt]
f_{-}^{(0)}f_{-}^{(1)}\end{pmatrix}$ $\displaystyle\hskip
10.00002pt+\frac{1}{2(f_{+}^{(0)}+f_{-}^{(1)})}\frac{d}{d\phi}\begin{pmatrix}f_{+}^{(0)}f_{+}^{(1)}\\\\[5.0pt]
f_{+}^{(0)}f_{-}^{(1)}\\\\[5.0pt] f_{+}^{(0)}f_{-}^{(1)}\\\\[5.0pt]
f_{-}^{(0)}f_{-}^{(1)}\end{pmatrix}$
$\displaystyle=\alpha(\Lambda)|r_{0}\rangle+\frac{1}{2(f_{+}^{(0)}+f_{-}^{(1)})}\frac{d}{d\phi}\begin{pmatrix}f_{+}^{(0)}f_{+}^{(1)}\\\\[5.0pt]
f_{+}^{(0)}f_{-}^{(1)}\\\\[5.0pt] f_{+}^{(0)}f_{-}^{(1)}\\\\[5.0pt]
f_{-}^{(0)}f_{-}^{(1)}\end{pmatrix},$ (81)

where $\alpha(\Lambda)$ is an unimportant factor because we are only
interested in $\langle\ell_{i}(\phi)|\frac{d}{d\phi}|r_{0}(\phi)\rangle$ when
$i\neq 0$ and $\langle\ell_{i}(\phi)|r_{0}(\phi)\rangle=0$ for $i\neq 0$.

Substituting Eq. (C.2) into Eq. (5) with $\theta=2\pi$ results in

$\displaystyle|\hat{\rho}(2\pi)\rangle\simeq|r_{0}(2\pi)\rangle+\sum_{1,3}\mathcal{C}_{i}|r_{i}(2\pi)\rangle,$
(82)
$\displaystyle\mathcal{C}_{i}=-\int_{0}^{2\pi}d\phi\frac{\langle\tilde{\ell}_{i}(\phi)|}{2(f_{+}^{(0)}+f_{-}^{(1)})}\frac{d}{d\phi}\begin{pmatrix}f_{+}^{(0)}f_{+}^{(1)}\\\\[5.0pt]
f_{+}^{(0)}f_{-}^{(1)}\\\\[5.0pt] f_{+}^{(0)}f_{-}^{(1)}\\\\[5.0pt]
f_{-}^{(0)}f_{-}^{(1)}\end{pmatrix}.$ (83)

By using $\frac{d}{d\theta}(f_{+}^{(j)}+f_{-}^{(j)})=0$, $\mathcal{C}_{i}$
becomes

$\displaystyle\mathcal{C}_{i}=$
$\displaystyle-\int_{0}^{2\pi}d\phi\frac{\langle\tilde{\ell}_{i}(\phi)|}{2(f_{+}^{(0)}+f_{-}^{(1)})}\frac{df_{+}^{(0)}}{d\phi}\begin{pmatrix}f_{+}^{(1)}\\\\[5.0pt]
f_{-}^{(1)}\\\\[5.0pt] f_{-}^{(1)}\\\\[5.0pt] -f_{-}^{(1)}\end{pmatrix}$
$\displaystyle-\int_{0}^{2\pi}d\phi\frac{\langle\tilde{\ell}_{i}(\phi)|}{2(f_{+}^{(0)}+f_{-}^{(1)})}\frac{df_{+}^{(1)}}{d\phi}\begin{pmatrix}f_{+}^{(0)}\\\\[5.0pt]
-f_{+}^{(0)}\\\\[5.0pt] -f_{+}^{(0)}\\\\[5.0pt] -f_{-}^{(0)}\end{pmatrix}.$
(84)

Substituting the expression of $\langle\ell_{1}|$ [44] into Eq. (C.2), we
obtain

$\displaystyle\mathcal{C}_{1}=$ $\displaystyle-\int_{0}^{2\pi}d\phi
e^{\int_{\phi}^{2\pi}\varepsilon_{1}(\xi)\frac{d\xi}{\epsilon}}\frac{2f_{-}^{(1)}}{f_{+}^{(0)}+f_{-}^{(1)}}\frac{df_{+}^{(0)}}{d\phi}$
$\displaystyle-\int_{0}^{2\pi}d\phi
e^{\int_{\phi}^{2\pi}\varepsilon_{1}(\xi)\frac{d\xi}{\epsilon}}\frac{2f_{+}^{(0)}}{f_{+}^{(0)}+f_{-}^{(1)}}\frac{df_{+}^{(1)}}{d\phi}.$
(85)

Similarly, we obtain

$\displaystyle\mathcal{C}_{3}=$ $\displaystyle\int_{0}^{2\pi}d\phi
e^{\int_{\phi}^{2\pi}\varepsilon_{3}(\xi)\frac{d\xi}{\epsilon}}\frac{f_{+}^{(1)}f_{-}^{(1)}}{f_{+}^{(0)}+f_{-}^{(1)}}\frac{df_{+}^{(0)}}{d\phi}$
$\displaystyle-\int_{0}^{2\pi}d\phi
e^{\int_{\phi}^{2\pi}\varepsilon_{3}(\xi)\frac{d\xi}{\epsilon}}\frac{f_{+}^{(0)}f_{-}^{(0)}}{f_{+}^{(0)}+f_{-}^{(1)}}\frac{df_{+}^{(1)}}{d\phi}$
$\displaystyle=$ $\displaystyle\int_{0}^{2\pi}d\phi
e^{4\epsilon^{-1}(\phi-2\pi)}\frac{f_{+}^{(1)}f_{-}^{(1)}}{f_{+}^{(0)}+f_{-}^{(1)}}\frac{df_{+}^{(0)}}{d\phi}$
$\displaystyle-\int_{0}^{2\pi}d\phi
e^{4\epsilon^{-1}(\phi-2\pi)}\frac{f_{+}^{(0)}f_{-}^{(0)}}{f_{+}^{(0)}+f_{-}^{(1)}}\frac{df_{+}^{(1)}}{d\phi}.$
(86)

As shown by Eq. (86), the factor $e^{4\epsilon^{-1}(\phi-2\pi)}$ in the
integrand plays an important role. Because of this factor, it is not necessary
to consider long-term memory in the dynamics. In the case of
$\mathcal{C}_{3}$, the scaling factor does not depend on the choice of the
trajectory, but only depends on $\theta$. Thus, one can estimate the BSN
curvature at $\phi$.

For $\epsilon\ll 1$, the exponential factor $e^{4\epsilon^{-1}(\phi-2\pi)}$
behaves as a cutoff function and thus

$\displaystyle\mathcal{C}_{3}\approx$
$\displaystyle\int_{2\pi-\epsilon}^{2\pi}d\phi\left[\frac{f_{+}^{(1)}f_{-}^{(1)}}{f_{+}^{(0)}+f_{-}^{(1)}}\frac{df_{+}^{(0)}}{d\phi}-\frac{f_{+}^{(0)}f_{-}^{(0)}}{f_{+}^{(0)}+f_{-}^{(1)}}\frac{df_{+}^{(1)}}{d\phi}\right]$
$\displaystyle\sim$
$\displaystyle\epsilon\left[\frac{f_{+}^{(1)}f_{-}^{(1)}}{f_{+}^{(0)}+f_{-}^{(1)}}\frac{df_{+}^{(0)}}{d\phi}-\frac{f_{+}^{(0)}f_{-}^{(0)}}{f_{+}^{(0)}+f_{-}^{(1)}}\frac{df_{+}^{(1)}}{d\phi}\right]_{\phi=2\pi}.$
(87)

Similarly, we can write

$\displaystyle\mathcal{C}_{1}\sim$
$\displaystyle\epsilon\left[\frac{2e^{\varepsilon_{1}(2\pi)}}{f_{+}^{(0)}+f_{-}^{(1)}}\left(f_{-}^{(1)}\frac{df_{+}^{(0)}}{d\phi}-f_{+}^{(0)}\frac{df_{+}^{(1)}}{d\phi}\right)\right]_{\phi=2\pi}.$
(88)

From the above estimates, the leading order contribution of the geometric
phase is the order $O(\epsilon)$.

Let us show that $\mathcal{C}_{3}$ reduces to zero in the non-interacting case
$\beta U_{0}=0$. In this case, $f^{(0)}_{\pm}=f^{(1)}_{\pm}\equiv g_{\pm}$ and
$g_{+}+g_{-}=2$; thus $\mathcal{C}_{3}$ becomes

$\displaystyle\lim_{\beta U_{0}\to 0}\mathcal{C}_{3}$ $\displaystyle=0.$ (89)

In the opposite limit $\beta U_{0}\to\infty$ with $f_{+}^{(1)}=0,\
f_{-}^{(1)}=2$ $\mathcal{C}_{3}$ also becomes zero, because

$\displaystyle\lim_{\beta U_{0}\to\infty}\mathcal{C}_{3}$ $\displaystyle=0.$
(90)

###  C.3 Some detailed results for the Anderson model

In this subsection, we present some detailed results beyond the main text, as
well as a figure showing the control parameters in the parameter space.

In Fig. 12, we plot $\mathcal{C}_{1}(\delta)$ (left figure) and
$\mathcal{C}_{3}(\delta)$ (right figure) versus $\delta$ for various $r$ and
$\beta U_{0}=0.1$. As shown, the coefficients increase with increasing $r$.
The coefficients $\mathcal{C}_{1}$ and $\mathcal{C}_{3}$ become zero at
$\delta=\pi$.

Figure 13 shows the coefficients $\mathcal{C}_{1}(\delta)$ (left figure) and
$\mathcal{C}_{3}(\delta)$ (right figure) versus $\delta$ for various $\beta
U_{0}$ and $r=0.9$. As shown in Fig. 13, $|\mathcal{C}_{i}|$ with $i=1$ and 3
increases with $\beta U_{0}$.

Figure 12: Plots of $\mathcal{C}_{1}$ (Fig. (a)) and $\mathcal{C}_{3}$ (Fig.
(b)) versus $\delta$ for $r=0.5,0.7,0.9$ with $\beta U_{0}=0.1$.  Figure 13:
Plots of $\mathcal{C}_{1}$ (Fig. (a)) and $\mathcal{C}_{3}$ (Fig. (b)) versus
$\delta$ for $\beta U_{0}=0.3,\ 0.5,\ 0.7$ with $r=0.9$.

##  Appendix D Some detailed properties of the density matrix, von Neumann
entropy and heat flow

In this section, we present the detailed results on the density matrix, the
von Neumann entropy, and the heat flow.

First, we verify the validity of the Born-like approximation by comparing its
solution with that obtained by the numerical solution of Eq. (4). In Fig. 14,
we plot the time dependence of the diagonal elements of the density matrix
computed by the Born-like approximation and the numerical solutions. This
shows that the Born-like approximation gives us sufficiently accurate
expressions since the deviation between the numerical solutions and the
expressions based on the Born-like approximation is invisible.

Figure 14: Plots of the time evolution of the density matrix element based on
the Born-like approximation (red-dashed line) and the numerical solutions of
the master equation (black-solid line). The parameters are set to be
$\beta\bar{\mu}=0.1$, $\beta U_{0}=0.1$, and $\beta\epsilon_{0}=0.1$.

Next, we plot $S^{\mathrm{vN}}(\theta,0)$ against $\theta$ in Fig. 15. This
figure contains the comparison of $S^{\mathrm{vN}}(\theta,0)$ obtained by the
Born-like approximation as in Eqs. (5) and (6) with that obtained by the
numerical integration of Eq. (4). As can be seen, the difference between the
two results is invisible. The behavior of $S^{\mathrm{vN}}(\theta,0)$ is
neither non-monotonic nor quasi-periodic.

Figure 15: Plot of von Neumann entropy as a function of $\theta$. The results
obtained by the Born-like approximation (red-dashed line) and the numerical
calculation of the master equation (black-solid line) are plotted.

Third, we also compute the heat flow $\mathcal{Q}(\theta,\delta)$. In Fig. 16,
we plot $\mathcal{Q}(\theta,0)$ as a function of $\theta$. Since the LHS of
Eq. (27) is much smaller than $S_{\rm sys}(\theta,\delta)$ and
$Q(\theta,\delta)$, $Q(\theta,0)$ is almost the same with $S_{\rm
sys}(\theta,0)$ shown in Fig. 9. The deviation between the solution of Eq. (4)
and the result of the Born-like approximation is invisible.

Figure 16:  Plot of the heat $\mathcal{Q}(\theta,0)$ versus $\theta$ for
$\beta U_{0}=0.1$ and $r=0.9$.

## References

  * [1] J. C. Maxwell, Theory of Heat (Appleton, London, 1871).
  * [2] L. Szilard, über die Entropieverminderung in einem thermodynamischen System bei Eingriffen intelligenter Wesen,  Z. Phys. 53, 840 (1929).
  * [3] L. Brillouin, Maxwell’s demon cannot operate: Information and entropy. I, J. Appl. Phys. 22, 334 (1951).
  * [4] L. Brillouin, Physical entropy and information. II, J. Appl. Phys. 22, 338 (1951).
  * [5] L. Brillouin, Science and information theory (Academic Press, New York, 1956).
  * [6] E. T. Jaynes, Information theory and statistical mechanics, Phys. Rev. 106, 620 (1957).
  * [7] R. Landauer, Irreversibility and heat generation in the computing process, IBM J. Res. Dev. 5, 183 (1961).
  * [8] C. H. Bennett, Logical reversibility of computation, IBM J. Res. Dev. 17, 525 (1973).
  * [9] J. V. Koski, V. F. Maisi, T. Sagawa, and J. P. Pekola, Experimental Observation of the Role of Mutual Information in the Nonequilibrium Dynamics of a Maxwell Demon, Phys. Rev. Lett. 113, 030601 (2014).
  * [10] M. D. Vidrighin, O. Dahlsten, M. Barbieri, M. S. Kim, V. Vedral, and I. A. Walmsley, Photonic Maxwell’s demon, Phys. Rev. Lett. 116, 050401 (2016).
  * [11] A. Shu, J. Dai, and V. Scarani, Power of an optical Maxwell’s demon in the presence of photon-number correlations, Phys. Rev. A 95, 022123 (2017).
  * [12] A. Rex, Maxwell’s demon-a historical review, Entropy, 17, 240 (2017).
  * [13] H. Touchette and S. Lloyd, Information-Theoretic Limits of Control.  Phys. Rev. Lett. 84, 1156 (2000).
  * [14] T. Sagawa and M. Ueda, Generalized Jarzynski Equality under Nonequilibrium Feedback Control.  Phys. Rev. Lett. 104, 090602 (2010).
  * [15] T. Sagawa and M. Ueda, Fluctuation Theorem with Information Exchange: Role of Correlations in Stochastic Thermodynamics Phys. Rev. Lett. 109, 180602 (2012).
  * [16] S. Viviana, C.-F. Lee, E. R. Key and D. A. Leigh, A molecular information ratchet, Nature. 445, 523 (2007).
  * [17] S. Toyabe, T. Sagawa, M. Ueda, E. Muneyuki and M. Sano, Experimental demonstration of information-to-energy conversion and validation of the generalized Jarzynski equality, Nature Phys. 6, 988 (2010).
  * [18] J. Parrondo, J. Horowitz and T. Sagawa, Nature Phys. 11, 131 (2015).
  * [19] M. V. Berry, Quantal phase factors accompanying adiabatic changes, Proc. R. Soc. London Ser. A 392, 45 (1984).
  * [20] D. J. Thouless, Quantization of particle transport, Phys. Rev. B 27, 6083 (1983).
  * [21] N. A. Sinitsyn and I. Nemenman, The Berry phase and the pump flux in stochastic chemical kinetics, Europhys. Lett. 77, 58001 (2007) .
  * [22] N. A. Sinitsyn and I. Nemenman, Universal Geometric Theory of Mesoscopic Stochastic Pumps and Reversible Ratchets, Phys. Rev. Lett. 99, 220408 (2007).
  * [23] L. P. Kouwenhoven, A. T. Johnson, N. C. van der Vaart, C. J. P. M. Harmans, and C. T. Foxon, Quantized current in a quantum-dot turnstile using oscillating tunnel barriers, Phys. Rev. Lett. 67, 1626 (1991).
  * [24] H. Pothier, P. Lafarge, C. Urbina, D. Esteve, and M. H. Devoret, Single-Electron Pump Based on Charging Effects, Europhys. Lett. 17, 249 (1992).
  * [25] M. Switkes, C. M. Marcus, K. Campman, and A. C. Gossard, An Adiabatic Quantum Electron Pump, Science 283, 1905 (1999).
  * [26] A. Fuhrer, C. Fasth, and L. Samuelson, Single electron pumping in InAs nanowire double quantum dots, Appl. Phys. Lett. 91, 052109 (2007).
  * [27] S. K. Watson, R. M. Potok, C. M. Marcus, and V. Umansky, Experimental Realization of a Quantum Spin Pump, Phys. Rev. Lett. 91, 258301 (2003).
  * [28] P. W. Brouwer, Scattering approach to parametric pumping, Phys. Rev. B 58, R10135 (1998).
  * [29] J. Ren, P. Hänggi, and B. Li, Berry-Phase-Induced Heat Pumping and Its Impact on the Fluctuation Theorem, Phys. Rev. Lett. 104, 170601 (2010).
  * [30] T. Sagawa and H. Hayakawa, Geometrical expression of excess entropy production, Phys. Rev. E 84, 051110 (2011).
  * [31] T. Yuge, T. Sagawa, A. Sugita, and H. Hayakawa, Geometrical pumping in quantum transport: Quantum master equation approach, Phys. Rev. B 86, 235308 (2012).
  * [32] T. Yuge, T. Sagawa, A. Sugiura, and H. Hayakawa, Geometrical Excess Entropy Production in Nonequilibrium Quantum Systems, J. Stat. Phys. 153, 412 (2013).
  * [33] K. L. Watanabe and H. Hayakawa, Geometric fluctuation theorem for a spin-boson system, Phys. Rev. E 96, 022118 (2017).
  * [34] Y. Hino and H. Hayakawa, Fluctuation relations for adiabatic pumping, Phys. Rev. E 102, 012115 (2020).
  * [35] K. Takahashi, Y. Hino, K. Fujii and H. Hayakawa, Full Counting Statistics and Fluctuation?Dissipation Relation for Periodically Driven Two-State Systems, J. Stat. Phys. 181, 2206 (2020).
  * [36] K Takahashi, K Fujii, Y Hino, and H Hayakawa, Nonadiabatic Control of Geometric Pumping, Phys. Rev. Lett. 124, 150602 (2020).
  * [37] Z. Wang, L. Wang, J. Chen, C. Wang and J. Rie, Geometric heat pump: Controlling thermal transport with time-dependent modulations Frontier, Phys. 17, 13201 (2022).
  * [38] G. E. Crooks, Measuring Thermodynamic Length, Phys. Rev. Lett. 99, 100602 (2007).
  * [39] B. Bhandari, P. T. Alonso, F. Taddei, F. von Oppen, R. Fazio and L. Arranchea, Geometric properties of adiabatic quantum thermal machines, Phys. Rev. B 102, 155407 (2020),
  * [40] P Abiuso, H. J. D. Miller, M. Perarnau-Llobet, and M Scandi, Geometric optimisation of quantum thermodynamic processes, Entropy 22, 1076 (2020).
  * [41] K. Brandner and K. Saito, Thermodynamic Geometry of Microscopic Heat Engines, Phys. Rev. Lett. 124, 040602 (2020).
  * [42] J. Lu, Z. Wang, J. Peng, C. Wang, J.-H. Jiang and J. Ren, Geometric thermodynamic uncertainty relation in a periodically driven thermoelectric heat engine Phys, Rev. B. 105, 115428 (2022).
  * [43] Y. Hino and H. Hayakawa, Geometrical formulation of adiabatic pumping as a heat engine, Phys. Rev. Research 3, 013187 (2021).
  * [44] H. Hayakawa, V. M. M. Paasonen, and R. Yoshii, “Geometrical Quantum Chemical Engine”, arXiv:2112.12370 (2021).
  * [45] S. Kullback, and R. A. Leibler, On information and sufficiency, Annal. Math. Stat., 22, 79 (1951).
  * [46] T. Sagawa, Entropy, Divergence, and Majorization in Classical and Quantum Thermodynamics, to be published in Springer-Brief in Mathematical Physics arXiv:2007.09974.
  * [47] D, Petz, Quasi-entropies for finite quantum systems, Rep. Math. Phys. 23, 57 (1986).
  * [48] D. Petz, Monotonicity of quantum relative entropy revisited, Rev. Math. Phys. 15, 79 (2003).
  * [49] M. B. Ruskai, Inequalities for quantum entropy: A review with conditions for equality, J. Math. Phys. 43, 4358 (2002).
  * [50] F. Hiai, M. Mosonyi, D. Petz, and C. Bény, Quantum f-divergences and error correction, Rev. Math. Phys. 23, 691 (2011).
  * [51] P. Strasberg, Quantum Stochastic Thermodynamics: Foundations and Selected Applications (Oxford Univ. Press, Oxford, 2022).
  * [52] G. Benenti, G. Casati, K. Saito, and R. S. Whitney, Phys. Rep. 694, 1 (2017).
  * [53] G. Kurizki and A. G. Kofman, Thermodynamics and Control of Open Quantum Systems (Cambridge Univ. Press, Cambridge, 2022).
  * [54] H. P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford Univ. Press, Oxford, 2002).
  * [55] A. A. Svidzinsky, K. E. Dorfman, M. O. Scully, Enhancing photocell power by noise-induced coherence, Coherent Opt. Phenomena, 1, 7 (2012).
  * [56] J. Um, K. E. Dorfman and H. Park, Coherence enhanced quantum-dot heat engine, arXiv:2111.09582.
  * [57] R. Yoshii and H. Hayakawa, “Analytical expression of geometrical pumping for a quantum dot based on quantum master equation”, arXiv:1312.3772