# Classical and quantum thermodynamics described as a system–bath model: The
dimensionless minimum work principle

Shoki Koyanagi  0000-0002-8607-1699<EMAIL_ADDRESS>and
<EMAIL_ADDRESS>Department of Chemistry, Graduate School of
Science, Kyoto University, Kyoto 606-8502, Japan  Yoshitaka Tanimura
0000-0002-7913-054X<EMAIL_ADDRESS>and
<EMAIL_ADDRESS>Department of Chemistry, Graduate School of
Science, Kyoto University, Kyoto 606-8502, Japan

(Last updated: )

###### Abstract

We formulate a thermodynamic theory applicable to both classical and quantum
systems. These systems are depicted as thermodynamic system–bath models
capable of handling isothermal, isentropic, thermostatic, and entropic
processes. Our approach is based on the use of a dimensionless thermodynamic
potential expressed as a function of the intensive and extensive thermodynamic
variables. Using the principles of dimensionless minimum work and
dimensionless maximum entropy derived from quasi-static changes of external
perturbations and temperature, we obtain the Massieu-Planck potentials as
entropic potentials and the Helmholtz-Gibbs potentials as free energy. These
potentials can be interconverted through time-dependent Legendre
transformations. Our results are verified numerically for an anharmonic
Brownian system described in phase space using the low-temperature quantum
Fokker–Planck equations in the quantum case and the Kramers equation in the
classical case, both developed for the thermodynamic system–bath model. Thus,
we clarify the conditions for thermodynamics to be valid even for small
systems described by Hamiltonians and establish a basis for extending
thermodynamics to non–equilibrium conditions.

##  I Introduction

Exactly 200 years have passed since Carnot published his work on the
efficiency of heat engines.[1] Thermodynamics describes macroscopic thermal
phenomena in equilibrium and quasi-static processes, independently of the
system dynamics. It is widely applied and has shown great success in
describing thermal phenomena characterized by intensive and extensive
thermodynamic properties. Massieu [2] and later Planck[3] combined the total
energy with temperature and entropy to derive the entropic potentials.[4]
Subsequently, the free energy was introduced by Gibbs and Helmholtz,
completing the foundations of thermodynamics.[5]

Statistical mechanics, on the other hand, is a system–specific theory that
describes many–body phenomena in an equilibrium state based on its microscopic
statistical properties.[6] It was initiated by Boltzmann’s establishment of
the equality $S=k_{\rm B}\ln{\mathcal{W}}$, where $S$ is the entropy,
${\mathcal{W}}$ is the number of possible microscopic states, and $k_{\rm B}$
is the Boltzmann constant, followed by the introduction of probability
distribution functions in phase space by Gibbs and Einstein, and then by the
1922 definition of partition functions as quantum discretized eigenstates.[7]
It has thus been a century since statistical mechanics took its current form.

Attempts to find a relationship between thermodynamics and statistical
mechanics have been successful under limited conditions for specific systems,
especially quantum systems, which are the subject of quantum thermodynamics.
[8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 17, 19, 20, 21, 22, 23, 24, 25, 26,
27, 28] However, a fundamental difference exists between the theoretical
foundations of statistical physics, which is based on a first-principles
description of kinetic systems from a microscopic perspective, and those of
thermodynamics, which relies on phenomenological description of thermal
systems from a macroscopic viewpoint. For example, in quantum mechanics,
observables are defined as expectation values, whereas in thermodynamics, they
are described by macroscopic intensive and extensive variables. Consequently,
a systematic theory based on a system–bath (SB) model that describes the
relationship between thermodynamic potentials characterized as intensive and
extensive variables has yet to be firmly established.

Recently, it has been shown that it is possible to evaluate the free energy on
the basis of an SB model from dynamically calculated work as $W\geq\Delta G$
(i.e., the minimum work principle), where $W$ is the work defined as being
done from the outside to the subsystem by time-dependent external fields under
an isothermal time-irreversible process, and $\Delta G$ is the change in free
energy. Thus, the free energy is evaluated as $\Delta G=W^{\rm qst}$, with the
system being driven quasi-statically by external fields.[26, 29] Using this
definition of thermodynamic potentials, the first and second laws of
thermodynamics have been verified under fully quantum conditions.[30, 31]

Note that the free energy defined by the partition function is ubiquitously
referred to as the Helmholtz energy $\Delta F$. However, in the thermodynamic
sense, the free energy as a function of the external field, which is an
intensive variable, should be called the “Gibbs energy,” while the minimum
work principle or the Kelvin–Planck statement is typically expressed as
$W\geq\Delta F$. Since both free energies play an essential role in the
present study, hereinafter we shall use the term “Gibbs energy” to refer to
what we have previously called the Helmholtz energy with reference both to our
own work and to that of others.

The SB model exhibits time-irreversible dynamics toward the thermal
equilibrium state of the total system owing to the infinite bath degrees of
freedom, while the total equilibrium state obeys a microcanonical
ensemble.[32, 33, 34, 35, 36, 37] The effect of a heat bath is described by
fluctuations and dissipation satisfying the fluctuation–dissipation theorem,
and the system reaches a thermal equilibrium state in which the energy
supplied by fluctuation and the energy lost by dissipation are balanced.[37,
38, 39, 40, 41] Note that in the theory of open quantum dynamics, fluctuations
are essentially non-Markovian noise expressed in terms of Matsubara
frequencies, even in the case of an Ohmic spectral distribution. Therefore, it
is important to adopt a non-Markovian treatment of fluctuations to describe
the correct thermal equilibrium state in which the subsystem and bath are
entangled (“bathentanglement”).[38, 39] Moreover, to calculate thermodynamic
variables, it is necessary to accurately evaluate the heat $Q(t)$, which is
the energy change in the heat bath, from the Hamiltonian dynamics.[27, 28] In
microscopic quantum systems, it is also essential to include contributions
from the SB interactions to maintain the first law of thermodynamics, which
corresponds to the energy conservation law for the total system.[29, 30, 31]

Although the SB interaction part of the internal energy and heat, whether in
quasi-static or non–equilibrium systems, is difficult to evaluate within the
conventional framework of open quantum dynamics theory because the bath
degrees of freedom have been reduced, we can evaluate these quantities
numerically by appropriate treatment of bathentanglement, such as by using
multiconfigurational time-dependent Hartree (MCTDH),[42, 43] polaron
transformation,[44] and other approaches.[45, 46, 47, 48, 49, 50, 51, 52]
Among these, the hierarchical equations of motion (HEOM) formalism provides a
stable and versatile scheme for performing numerical simulations of a wide
range of problems,[41, 37, 38, 53, 54, 55, 56, 39, 40, 57, 58, 59, 60, 61]
such as for spin–lattice systems,[62] the Hornstein–Hubbard model,[63] and
vibrational modes of liquid water.[64, 65]

Although investigations of quantum thermodynamics have progressed with the
help of experimental advances, [66, 67, 68, 69] several issues remain unclear.
For example, extensive variables such as entropy and susceptibility have been
obtained from the free energy, since they satisfy Legendre transformations,
but how are these defined in the SB model? What is meant by the temperature
derivative of the free energy evaluated from the isothermal process described
by the SB model? Are the arguments discussed in the theory of open quantum
dynamics consistent with ordinary thermodynamics in the classical limit? The
answers to these questions should provide a basis for extending
thermodynamics, defined by quasi-static processes, to non–equilibrium
processes.

In this paper, we develop a thermodynamic theory applicable to classical and
quantum systems, expressed as time-dependent thermodynamic potentials as
functions of extensive and intensive variables, based on the
Ullersma–Caldeira–Leggett (or Brownian) model[32, 33, 34, 35, 36] extended for
thermodynamic studies. Our argument is based on the dimensionless minimum work
principle, which is expressed using the dimensionless (or entropic) work and
entropic potentials. The results are verified by numerical simulations based
on the low-temperature quantum Fokker–Planck equations (LT-QFPE)[57] and the
Kramers equation in the classical limit.

The remainder of this paper is organized as follows. In Sec. II, we introduce
the SB Hamiltonian for thermodynamic processes. We then define non–equilibrium
enthalpy and internal energies. In Sec. III, we introduce dimensionless
thermodynamic variables such as dimensionless entropy and polarization. Then,
in Sec. IV, we discuss the laws of thermodynamics within the framework of open
quantum dynamics theory. Using the quasi-static values of work and heat, we
introduce the dimensionless minimum work principle and maximum entropy
principle. We then express the thermodynamic potentials in terms of intensive
and extensive variables in total differential form. In Sec. V, to verify our
results in a numerically rigorous manner, we introduce quantum and classical
reduced equations of motion for thermodynamic processes and perform
simulations. Finally, in Sec. VI, we present concluding remarks.

##  II Thermodynamic system–bath model

###  II.1 Model Hamiltonian

We develop a thermodynamic theory applicable to classical and quantum systems
that can describe isothermal, isentropic, thermostatic, and entropic
processes. To achieve this, we extend the Ullersma–Caldeira–Leggett model (or
Brownian model),[32, 33, 34, 35, 36, 55, 56, 39, 37, 38, 40, 57] which
consists of a subsystem A defined in phase space coupled to a heat bath B. By
using this model instead of the spin–boson model,[45, 53, 54] we can obtain
the classical results by taking $\hbar\rightarrow 0$. The total Hamiltonian is
then expressed as

$\displaystyle\hat{H}_{\rm tot}(t)=\hat{H}_{\rm A}(t)+\hat{H}_{\rm I+B}(t),$
(1)

where

$\displaystyle\hat{H}_{\rm A}(t)=\hat{H}_{\rm A}^{0}+\hat{H}^{\prime}_{\rm
A}(t),$ (2)

with

$\displaystyle\hat{H}_{\rm A}^{0}\equiv\frac{\hat{p}^{2}}{2m}+U(\hat{q})$ (3)

being the Hamiltonian for the subsystem with mass $m$ and potential
$U(\hat{q})$ described by the momentum $\hat{p}$ and position $\hat{q}$, and

$\displaystyle\hat{H}^{\prime}_{\rm A}(t)\equiv-E(t){\mu}(\hat{q})$ (4)

is the perturbation described as ${\mu}(\hat{q})$ as a function of the
subsystem coordinate and the external field $E(t)$, which is a
thermodynamic–intensive variable.[27, 28] In a typical example, $E(t)$ denotes
the laser field and ${\mu}(\hat{q})$ the electric dipole moment. In laser
spectroscopy, the expectation value of the dipole is an extensive variable and
is observed as the polarization $P_{\rm A}(t)=\langle\mu(\hat{q})\rangle$,
where $\langle\cdots\rangle$ denotes the thermal average.[37]

Although the conventional SB model has been limited to the investigation of
isothermal processes at a constant temperature, here we extend it to describe
thermostatic processes in which temperature varies with time by introducing
multiple heat baths. Thus, we consider a situation in which $N$ independent
heat baths, each in the thermal equilibrium state at the inverse temperature
$\beta_{k}\equiv 1/k_{\rm B}T_{k}$ are connected to or disconnected from the
subsystem A using the window function $\xi_{k}(t)$. The bath part of the
Hamiltonian is expressed as follows:

$\displaystyle\hat{H}_{\rm I+B}(t)=\sum_{k=1}^{N}\left(\hat{H}_{\rm
B}^{k}+\xi_{k}(t)\hat{H}_{\rm I}^{k}\right).$ (5)

Here, the $k$th bath Hamiltonian is expressed as an ensemble of harmonic
oscillators and is given by

$\displaystyle\hat{H}_{\rm
B}^{k}\equiv\sum_{j}\left\\{\frac{(\hat{p}^{k}_{j})^{2}}{2m_{j}}+\frac{1}{2}m_{j}^{k}(\omega_{j}^{k})^{2}(\hat{x}_{j}^{k})^{2}\right\\},$
(6)

where the momentum, position, mass, and frequency of the $j$th bath oscillator
are given by $\hat{p}_{j}^{k}$, $\hat{x}_{j}^{k}$, $m_{j}^{k}$, and
$\omega_{j}^{k}$, respectively. Here, we consider the situation where the
$k$th bath is always in the thermal equilibrium state
$\exp(-\beta_{k}\hat{H}_{\rm B}^{k})$ and heat is transferred to the bath only
when the bath is connected to the subsystem.

The SB interaction, including the counter term, is expressed as

$\displaystyle\hat{H}_{\rm
I}^{k}\equiv\sum_{j}\left\\{-A_{k}V(\hat{q})c_{j}^{k}\hat{x}_{j}^{k}+\frac{(c_{j}^{k})^{2}A_{k}^{2}V^{2}(\hat{q})}{2m_{j}^{k}(\omega_{j}^{k})^{2}}\right\\},$
(7)

where $c_{j}^{k}$ is the SB coupling coefficient of the $j$th bath oscillator,
and $A_{k}$ and $V(\hat{q})$ are the coupling strength and the subsystem part
of the $k$th SB interaction, respectively. In the study of optical
spectroscopy, $V(\hat{q})$ describes the effects of vibrational relaxation and
dephasing.[70, 64, 65] The counterterm is introduced to preserve the
translational symmetry of the total Hamiltonian at $U({\hat{q}})=0$ and to
remove the undesired self-energy divergence that occurs in the Markovian
case.[32, 33, 34, 35, 36, 55, 56, 39, 37, 38, 40, 57]

The major difference when evaluating thermodynamic properties is that in the
present model, the interaction term $\hat{H}_{I}^{k}$ is reduced with the
bath, whereas in the spin–boson model,[45, 54] it is treated separately from
the bath.[28, 29, 30, 31] This makes the classical limit of the thermal
equilibrium state of the time-independent subsystem [$E(t)=0$] at $\beta_{k}$
identical to that of the isolated subsystem as ${\hat{\rho}}_{\rm A}^{\rm
eq}=\exp(-\beta_{k}\hat{H}_{\rm A}^{0})$, not as ${\hat{\rho}}_{\rm A}^{\rm
eq}={\rm tr}_{\rm B}\\{\exp[-\beta_{k}(\hat{H}_{\rm A}^{0}+\hat{H}_{\rm
I+B}^{k})]\\}$, while the quantum equilibrium state is still different from
the isolated one, owing to the bathentanglement.[38, 39] This is a desirable
feature for comparison with conventional (classical) thermodynamic results,
where the interaction is not explicitly considered.

The $k$th bath at $\beta_{k}$ is characterized by a spectral distribution
function (SDF) defined as

$\displaystyle
J^{k}(\omega)\equiv\sum_{j}\frac{\hbar(c_{j}^{k}A_{k})^{2}}{2m_{j}^{k}\omega_{j}^{k}}\delta(\omega-\omega_{j}^{k}).$
(8)

The window function, which we call the “thermostatic field,” is, for example,
defined as

$\displaystyle\xi_{k}(t)=\theta(t-t_{k})\theta(t_{k}+\Delta t-t),$ (9)

where $\theta(t)$ is the step function and the time $t_{k}$ is defined as
$t_{k}=t_{0}+(k-1)\Delta t$, with initial time $t_{0}$ and duration $\Delta
t$.

For a fixed temperature $\beta_{k}=\beta$, $A_{k}(t)\equiv A_{k}\xi_{k}(t)$
[i.e., $A_{k}^{2}(t)=A_{k}^{2}\xi_{k}(t)$] is the “adiabatic transition field”
introduced to describe isothermal–adiabatic manipulations (e.g., the
insertion/removal of an adiabatic wall or the connection/disconnection of the
subsystem and bath[30, 31]).

For a fixed coupling strength $A_{k}=A$, the thermostatic processes can be
described by utilizing $\xi_{k}(t)$ for multiple heat baths with different
temperatures $T_{k}$ to set only one bath as being connected to a subsystem at
a time. In this case, the bath temperature is effectively expressed as

$\displaystyle T(t)=\sum_{k=1}^{N}T_{k}\xi_{k}(t),$ (10)

or the inverse temperature as $\beta(t)=[k_{\rm B}T(t)]^{-1}$, except in the
adiabatic case.

###  II.2 Quantum fluctuation–dissipation theorem

For the $k$th bath, if we consider the interaction coordinate of the bath
modes as $\hat{X}^{k}\equiv\sum_{j}c_{j}^{k}\hat{x}_{j}^{k}$, the subsystem A
is driven by the external force $\hat{X}^{k}(t)$ through the interaction
$-V(\hat{q})\hat{X}^{k}$, where $\hat{X}^{k}(t)$ is the Heisenberg
representation of $\hat{X}^{k}$ for $\hat{H}_{\rm B}^{k}$. Because each bath
is a harmonic that is Gaussian in nature, the character of $\hat{X}^{k}(t)$ is
specified by its two-time correlation functions, such as the symmetrized and
canonical correlation functions defined by[41, 37, 39, 38, 55, 56, 39, 37, 38,
40, 57]

$\displaystyle
C^{k}(t)=\dfrac{1}{2}\langle\hat{X}^{k}(t)\hat{X}^{k}(0)+\hat{X}^{k}(0)\hat{X}^{k}(t)\rangle_{\mathrm{B}}$
(11)

and

$\displaystyle\begin{split}R^{k}(t)&=\frac{\beta_{k}}{2}\langle\hat{X}^{k};\hat{X}^{k}(t)\rangle_{\mathrm{B}}\\\
&\equiv\frac{1}{2}\int_{0}^{\beta_{k}}{d\lambda}\,\left\langle{\hat{X}}^{k}(-{\rm
i}\hbar\lambda)\hat{X}^{k}(t)\right\rangle_{\mathrm{B}},\end{split}$ (12)

where $\langle\cdots\rangle_{\mathrm{B}}$ represents the thermal average of
the $k$th bath degree of freedom.

The function $C^{k}(t)$ corresponds to fluctuation and is analogous to the
classical correlation function $X^{k}(t)$, whereas $R^{k}(t)$ corresponds to
dissipation. Both are induced by the bath, and they are related through the
fluctuation–dissipation theorem in Fourier form as
$C^{k}[\omega]=\hbar\omega\coth(\beta_{k}\hbar\omega/2)R^{k}[\omega]$.[41]
Note that in the SB system, the fluctuation–dissipation relation is the
condition to have the thermal equilibrium state, while the detailed balance
condition is not satisfied under the strong or non-Markovian SB coupling
conditions.[37, 38, 39, 40]

With the SDF, these are expressed as

$\displaystyle C^{k}(t)=\int_{0}^{\infty}d\omega
J^{k}(\omega)\coth\left(\frac{\beta_{k}\hbar\omega}{2}\right)\cos(\omega t)$
(13)

and

$\displaystyle
R^{k}(t)=\int_{0}^{\infty}\\!d\omega\,\frac{J^{k}(\omega)}{\hbar\omega}\cos(\omega
t).$ (14)

In any situation, by choosing the SDF in a continuous form, the bath degrees
of freedom become effectively infinite. The energy eigenstates of the total
system then become continuous and exhibit time-irreversible dynamics toward
the thermal steady state, with or without time-dependent external fields,
while the total system remains an isolated Hamiltonian system that conserves
energy. When $E(t)=0$, if the change in $T_{k}$ is small and $\Delta t$ is
large, the total system enters a canonical distribution described as
$\hat{\rho}^{\rm eq}_{\rm tot}(t_{k})\approx\exp[-\beta_{k}(\hat{H}_{\rm
A}^{0}+\hat{H}_{\rm I+B}^{k})]$, which can be verified from the steady-state
solution of the HEOM and from the solution of the imaginary HEOM[39] (see also
Appendix A).

In this study, we focus on the quasi-static case and assume that the
correlation time of the noise is much shorter than the characteristic time
scale of the system dynamics. Thus, we consider the Ohmic case described by

$\displaystyle J^{k}(\omega)=\frac{\hbar A_{k}^{2}\omega}{\pi}.$ (15)

The correlation functions in Eqs. (13) and (14) are then evaluated as[57]

$\displaystyle\begin{split}C^{k}(t)&\simeq\frac{2A_{k}^{2}}{\beta_{k}}\left(1+\sum_{l=1}^{K}2\right)\delta(t)-\sum_{l=1}^{K}\frac{2A_{k}^{2}\nu_{l}^{k}}{\beta_{k}}e^{-\nu_{l}^{k}|t|}\end{split}$
(16)

and

$\displaystyle R^{k}(t)=A_{k}^{2}\delta(t).$ (17)

Note that, in this Ohmic SDF, some of the physical observables, including the
mean square of the momentum, $\langle p^{2}\rangle$, diverge without the
cutoff $K$ because of the divergence of the first and second terms in Eq. (16)
under the infinite summation over $l$, which is often referred to as
ultraviolet divergence.[35, 36, 71, 72] However, because the contributions of
random forces from such terms are averaged over a sufficiently short time,
their effect on the dynamics of interest can be ignored by choosing a finite
$K$.[57]

We choose the coefficients $\nu_{l}^{k}$ and $\eta_{l}^{k}$ to realize the
relation for finite $K$, where the first term on the right-hand side of Eq.
(16) is the classical contribution from the temperature, and the remaining
terms are the quantum low-temperature (QLT) corrections.[57] Note that the
fluctuation term is always non-Markovian because of the quantum nature of the
noise; it can be regarded as Markovian only in the high-temperature limit,
$\beta_{k}\hbar\omega_{\mathrm{0}}\ll 1$, in which the heat bath exhibits
classical behavior.[37, 36] This is an important conclusion obtained from the
quantum fluctuation–dissipation theorem. That is, negative non-Markov terms
always appear in the Ohmic case unless a time coarse-grained Markov assumption
is adopted, which is often unphysical as a description of quantum
thermodynamics.[37, 38, 39, 40]

###  II.3 Non–equilibrium enthalpy and internal energy

In a typical thermodynamic theory, a system is characterized by thermodynamic
potentials described in terms of intensive variables such as electric field
$E$ and temperature $T$ and extensive variables such as polarization $P$ and
entropy $S$. To construct a quantum thermodynamic theory similar to
thermodynamics, isothermal and isentropic (d$E$ and d$P$) processes, and
thermostatic and entropic (d$T$ and d$S$) processes must be investigated.
However, most investigations including traditional thermodynamics to date have
been limited to isothermal processes in which only the external field
corresponding to the intensive variable is manipulated.

Although various quantum thermodynamic studies have been conducted, it has
been difficult to derive thermodynamic laws involving the total derivative of
temperature (d$T$), owing to the difficulty of including thermostatic
processes. Moreover, these studies have not introduced extensive variables
such as $P$ as observables in open–quantum systems, nor have they provided
their Legendre transformations, which play an essential role in
thermodynamics.

Here, we show that we can construct a complete description of a thermodynamic
theory from work in isothermal processes ($T_{\rm fix}$) and heat in constant-
external-field (or thermostatic) processes ($E_{\rm fix}$). For this purpose,
we introduce the expectation values of energy for each component of the
Hamiltonian. From the Hamiltonian [Eqs. (2)–(4)], the total energy of the
subsystem is expressed as

$\displaystyle H_{\rm A}(t)=U_{\rm A}(t)+H^{\prime}(t),$ (18)

where

$\displaystyle U_{\rm A}(t)\equiv\mathrm{tr}_{\rm A}\\{\hat{H}_{\rm
A}^{0}\hat{\rho}_{\rm A}(t)\\}$ (19)

is the self-energy of the subsystem, which is considered “internal energy” in
a quasi-static case, and

$\displaystyle H^{\prime}(t)=-E(t)P_{\rm A}(t)$ (20)

is the interaction energy described with the optical polarization defined as

$\displaystyle P_{\rm A}(t)\equiv\mathrm{tr}_{\rm
A}\left\\{{\mu}(\hat{q})\hat{\rho}_{\rm A}(t)\right\\},$ (21)

where $\hat{\rho}_{\rm A}(t)=\mathrm{tr}_{\rm B}\\{\hat{\rho}_{\rm tot}(t)\\}$
is the reduced density operator for the total density operator
$\hat{\rho}_{\rm tot}(t)$.

While $E(t)$ is an intensive variable, $P_{\rm A}(t)$ is an extensive
variable. This is because for distinguishable $M$ subsystems without mutual
interaction, the total density operator is expressed as $\hat{\rho}_{\rm
A}^{(M)}=\hat{\rho}_{\rm A}\otimes\cdots\otimes\hat{\rho}_{\rm A}$, and thus
$P_{\rm A}(t)$ becomes proportional to the size of the system. Moreover,
because the evaluation of $P_{\rm A}(t)$ requires a statistical average of the
distribution function $\hat{\rho}_{\rm A}(t)$ even for a single-particle
system, we can evaluate it as a thermal variable.

Because $E(t)$ and $P_{\rm A}(t)$ are conjugate to each other, the relation

$\displaystyle H_{\rm A}(t)=U_{\rm A}(t)-E(t)P_{\rm A}(t)$ (22)

is regarded as a time-dependent Legendre transformation between
non–equilibrium enthalpy $H_{\rm A}(t)$ and internal energy $U_{\rm A}(t)$.

The total non–equilibrium enthalpy is now expressed as

$\displaystyle H_{\rm tot}(t)=H_{\rm A}(t)+\sum_{k=1}^{N}H_{\rm I+B}^{k}(t),$
(23)

where $H_{\rm I+B}^{k}(t)$ is the enthalpy of the $k$th bath, defined as

$\displaystyle H_{\rm I+B}^{k}(t)=\mathrm{tr}_{\rm tot}\left\\{\hat{H}_{\rm
I+B}^{k}(t)\hat{\rho}_{\rm tot}(t)\right\\}.$ (24)

The bath part of the energy (non–equilibrium bath enthalpy) is then expressed
as

$\displaystyle H_{\rm I+B}(t)\equiv\sum_{k=1}^{N}H_{\rm I+B}^{k}(t).$ (25)

##  III Dimensionless thermodynamic variables

Although the minimum work principle for an intensive variable, expressed as
$W_{\rm A}^{int}(t)\geq\Delta G_{\rm A}^{\rm qst}$, provides the condition for
determining the thermodynamic potential from work,[29] it does not do so for
heat, which is important in determining entropy. We find that this difficulty
can be overcome by introducing the dimensionless enthalpy of the subsystem,
defined as $\tilde{H}_{\rm A}(t)\equiv\beta(t)H_{\rm A}(t)$. We use this to
define the changes in dimensionless (or entropic) “intensive work” and
“extensive heat” in time corresponding to power and heat flow, respectively,
as follows:

$\displaystyle\frac{d\tilde{W}_{\rm A}^{int}(t)}{dt}\equiv\mathrm{tr}_{\rm
A}\left\\{\frac{\partial}{\partial t}\left[\beta(t)\hat{H}_{\rm
A}(t)\right]\hat{\rho}_{\rm A}(t)\right\\}$ (26)

and

$\displaystyle\frac{d\tilde{Q}_{\rm A}^{ext}(t)}{dt}\equiv\mathrm{tr}_{\rm
A}\left\\{\left[\beta(t)\hat{H}_{\rm A}(t)\right]\frac{\partial\hat{\rho}_{\rm
A}(t)}{\partial t}\right\\}.$ (27)

These are interrelated by the following time-dependent Legendre
transformation:

$\displaystyle\begin{split}\frac{d\tilde{W}^{int}_{\rm
A}(t)}{dt}&=-\frac{d\tilde{Q}_{\rm
A}^{ext}(t)}{dt}+\frac{d}{dt}\left[\beta(t)H_{\rm A}(t)\right].\end{split}$
(28)

They are then evaluated as

$\displaystyle\frac{d\tilde{W}^{int}_{\rm A}(t)}{dt}=H_{\rm
A}(t)\frac{d\beta(t)}{dt}-\tilde{P}_{\rm A}(t)\frac{dE(t)}{dt}$ (29)

and

$\displaystyle\frac{d\tilde{Q}_{\rm A}^{ext}(t)}{dt}=\beta(t)\frac{dH_{\rm
A}(t)}{dt}+\tilde{P}_{\rm A}(t)\frac{dE(t)}{dt},$ (30)

where $\tilde{P}_{\rm A}(t)$ is the dimensionless polarization, defined as

$\displaystyle\tilde{P}_{\rm A}(t)=\beta(t)P_{\rm A}(t).$ (31)

In quantum thermodynamics, work has been defined using the change in an
intensive variable as $P_{\rm A}(t)dE(t)$. In typical thermodynamics, however,
work can also be defined as the change in an extensive variable as
$-E(t)dP_{\rm A}(t)$. To treat work defined as such, we further introduce the
dimensionless (or entropic) “extensive work” using a time-dependent Legendre
transformation of $\tilde{W}^{int}_{\rm A}(t)$ as follows:

$\displaystyle\begin{split}\frac{d\tilde{W}_{\rm
A}^{ext}(t)}{dt}&=\frac{d\tilde{W}^{int}_{\rm
A}(t)}{dt}+\frac{d}{dt}\left[\tilde{P}_{\rm A}(t)E(t)\right]\\\ &=H_{\rm
A}(t)\frac{d\beta(t)}{dt}+E(t)\frac{d\tilde{P}_{\rm A}(t)}{dt}.\end{split}$
(32)

The fundamental difference between the present treatment and the conventional
treatment of thermodynamics is that the total energy, including the heat bath,
is explicitly described by the Hamiltonian. Thus, we consider the
dimensionless total enthalpy expressed as

$\displaystyle\tilde{H}_{\rm tot}(t)=\mathrm{tr}_{\rm
tot}\left\\{\hat{\tilde{H}}_{\rm tot}(t)\hat{\rho}_{\rm tot}(t)\right\\},$
(33)

where $\hat{\tilde{H}}_{\rm tot}(t)$ is the dimensionless total Hamiltonian
defined as

$\displaystyle\hat{\tilde{H}}_{\rm tot}(t)\equiv\beta(t)\hat{H}_{\rm
A}(t)+\sum_{k=1}^{N}\beta_{k}\hat{H}_{\rm I+B}^{k}(t).$ (34)

The relationship between dimensionless (or entropic) total work and heat is
then expressed as

$\displaystyle\frac{d\tilde{W}_{\rm tot}^{int}(t)}{dt}=-\frac{d\tilde{Q}_{\rm
tot}^{ext}(t)}{dt}+\frac{d}{dt}\tilde{H}_{\rm tot}(t),$ (35)

where

$\displaystyle\quad\frac{d\tilde{W}_{\rm tot}^{int}(t)}{dt}=\mathrm{tr}_{\rm
tot}\left\\{\frac{\partial\hat{\tilde{H}}_{\rm tot}(t)}{\partial
t}\hat{\rho}_{\rm tot}(t)\right\\}$ (36)

and

$\displaystyle\frac{d\tilde{Q}_{\rm tot}^{ext}(t)}{dt}=0.$ (37)

Here, to obtain Eq. (37), we consider the case in which the subsystem is
connected to only one bath at the same time and use the following identity:

$\displaystyle\mathrm{tr}_{\rm tot}\left\\{\hat{H}_{\rm
tot}(t){\frac{\partial\hat{\rho}_{\rm tot}(t)}{\partial t}}\right\\}=0.$ (38)

From Eq. (35), we have

$\displaystyle\frac{d\tilde{W}_{\rm tot}^{int}(t)}{dt}=\frac{d\tilde{H}_{\rm
tot}(t)}{dt}.$ (39)

Thus, the thermal current between the subsystem and bath is conserved as

$\displaystyle\frac{d\tilde{Q}_{\rm A}^{ext}(t)}{dt}=-\frac{d\tilde{Q}_{\rm
I+B}^{ext}(t)}{dt},$ (40)

where

$\displaystyle\frac{d\tilde{Q}_{\rm
I+B}^{ext}(t)}{dt}=\sum_{k=1}^{N}{\mathrm{tr}_{\rm
tot}\left\\{\beta_{k}\hat{H}_{\rm I+B}^{k}(t)\frac{\partial\hat{\rho}_{\rm
tot}(t)}{\partial t}\right\\}}.$ (41)

We then identify the time derivative of the dimensionless bath entropy, i.e.,
$\tilde{S}_{\rm I+B}(t)={S}_{\rm I+B}(t)/k_{\rm B}$, with the dimensionless
bath heat current as follows:

$\displaystyle{\frac{d\tilde{S}_{\rm I+B}(t)}{dt}=\frac{d\tilde{Q}_{\rm
I+B}^{ext}(t)}{dt}.}$ (42)

##  IV The laws of thermodynamics

Quantum mechanics is a first-principles theory, and its logical structure is
reasonably simple. By contrast, thermodynamics is a phenomenological theory
developed on the basis of several principles and statements. In this section,
we deduce thermodynamics from quantum mechanics using observables, defined as
quantum mechanical expectation values, to clarify the central dogma of
thermodynamics.

###  IV.1 First to third laws of thermodynamics

####  IV.1.1 First law: Energy conservation law

The key principle for energy in quantum mechanics is the energy conservation
law for work, whereas in thermodynamics, energy is described not only as work
but also as heat. In fact, by adding Eqs. (30) and (32), we have the first law
of thermodynamics expressed as

$\displaystyle\frac{d\tilde{U}_{\rm A}(t)}{dt}=\frac{d\tilde{W}_{\rm
A}^{ext}(t)}{dt}+\frac{d\tilde{Q}_{\rm A}^{ext}(t)}{dt},$ (43)

where

$\displaystyle\tilde{U}_{\rm A}(t)=\beta(t)H_{\rm A}(t)+E(t)\tilde{P}_{\rm
A}(t)$ (44)

is the dimensionless non–equilibrium internal energy expressed also as
$\tilde{U}_{\rm A}(t)=\beta(t)U_{\rm A}(t)$. The above equality is the
consequence of the energy conservation law expressed as the time derivative of
Eq. (23) as follows:

$\displaystyle\frac{dH_{\rm tot}(t)}{dt}=\frac{dH_{\rm
A}(t)}{dt}+\frac{dH_{\rm I+B}(t)}{dt}.$ (45)

Then, using Eqs. (38) and (40), we obtain Eq. (43).

####  IV.1.2 Second law: Increasing internal energy under time-irreversible
process

As an extension of the minimum work principle for work as being done from the
outside to the subsystem, we consider “the dimensionless (or entropic) minimum
work principle” for the total system from one equilibrium state to another,
expressed as

$\displaystyle\tilde{W}_{\rm tot}^{int}\geq\left(\tilde{W}^{int}_{\rm
tot}\right)^{\rm qst},$ (46)

where $\tilde{W}_{\rm tot}^{int}$ is defined by Eq. (36), and
$(\tilde{W}^{int}_{\rm tot})^{\rm qst}$ represents a transition that occurs
quasi-statically. The proof of the above inequality is given in Appendixes A
and B. Unlike the existing minimum work principle, this definition can treat a
thermostatic process. From Eq. (39), the above inequality can also be
expressed as

$\displaystyle\Delta\tilde{H}_{\rm tot}\geq\Delta\tilde{H}_{\rm tot}^{\rm
qst},$ (47)

which corresponds to “the law of increasing enthalpy” under a time-
irreversible process for a closed system.

When the contribution of work from the SB interaction is incorporated into the
bath, which is the case in the present Brownian-based model, we have the
equality

$\displaystyle\frac{d\tilde{W}_{\rm A}^{int}(t)}{dt}=\frac{d\tilde{W}_{\rm
tot}^{int}(t)}{dt},$ (48)

and Eq. (46) reduces to

$\displaystyle\tilde{W}_{\rm A}^{int}\geq\left(\tilde{W}^{int}_{\rm
A}\right)^{\rm qst}.$ (49)

This relation is an extension to thermostatic processes of the minimum work
principle for isothermal processes, which states that thermodynamic weight is
maximized in a quasi-static equilibrium state. Because not only the intensive
variables $\beta(t)$ and $E(t)$ but also the extensive variables
$\tilde{P}_{\rm A}(t)$ and $H_{\rm A}(t)$ may change independently during a
non–equilibrium transition between the two equilibrium states, we can obtain
the relationship for heat with the use of Eq. (28) as

$\displaystyle\tilde{Q}_{\rm A}^{ext}\leq\left(\tilde{Q}_{\rm
A}^{ext}\right)^{\rm qst},$ (50)

which is equivalent to “the principle of maximum entropy.” This is because,
for the transition from equilibrium state 1 at time $t_{1}$ to equilibrium
state 2 at $t_{2}$, we have

$\displaystyle\int^{t_{2}}_{t_{1}}\frac{1}{T(t^{\prime})}\frac{dQ_{\rm
A}(t^{\prime})}{dt^{\prime}}dt^{\prime}\leq\Delta S_{\rm A}^{\rm qst},$ (51)

where

$\displaystyle{\frac{dQ_{\rm A}(t)}{dt}=\mathrm{tr}_{\rm
A}\left\\{\hat{H}_{\rm A}(t)\frac{\partial\hat{\rho}_{\rm A}(t)}{\partial
t}\right\\}}$ (52)

is the system heat current,[28] which satisfies $\beta(t)dQ_{\rm
A}(t)/dt=d\tilde{Q}_{\rm A}^{ext}(t)/dt$.

From Eqs. (40) and (51), the law of total entropy production is expressed as

$\displaystyle\Delta S_{\rm tot}\geq 0,$ (53)

where $\Delta S_{\rm tot}$ is a total entropy difference between the
equilibrium states at time $t_{1}$ and $t_{2}$, defined as

$\displaystyle\Delta S_{\rm tot}={\Delta S_{\rm A}^{\rm
qst}}+\int^{t_{2}}_{t_{1}}\frac{dS_{\rm I+B}(t)}{dt}dt.$ (54)

Here, $S_{\rm I+B}(t)=k_{\rm B}\tilde{S}_{\rm I+B}(t)$ is the bath entropy
[see Eq. (42)].

We can also obtain the inequality for dimensionless (or entropic) extensive
work using the Legendre transformation (32) for Eq. (49) as

$\displaystyle\tilde{W}_{\rm A}^{ext}\geq\left(\tilde{W}_{\rm
A}^{ext}\right)^{\rm qst}.$ (55)

####  IV.1.3 Third law: Zero-temperature limit of SB model

In the framework of open quantum dynamics theory, the third law of
thermodynamics should describe the state of the subsystem in the limit of zero
temperature, but anomalous behavior has been found in the spin–boson system at
zero temperature.[45, 46, 43, 44, 47] Thus, entropy may not be zero even at
zero temperature, because of the essential role played by bathentanglement:
the entropy may diverge as $\beta\rightarrow\infty$ as the classical
approximation ($\beta\hbar\omega_{\rm A}\ll 1$, where $\omega_{\rm A}$ is the
characteristic frequency of the subsystem) collapses in the low-temperature
regime. This implies that some statements of the third law of thermodynamics,
such as “the entropy of a system at absolute zero is a well-defined constant,”
do not hold for quantum systems described by an SB model, whereas Nernst’s
statement “it is impossible for any procedure to reach $T=0$ in a finite
number of steps,” which is equivalent to the assertion that “thermodynamic
states with $T=0$ do not exist” seems reasonable—although theoretically
anything is possible.

###  IV.2 Thermodynamic potentials

####  IV.2.1 Massieu–Planck potentials: Entropic state representation

Quasi-static dimensionless (or entropic) work and heat are related to the
(dimensionless) Massieu–Planck potentials, which are entropic thermodynamic
potentials defined as Legendre transforms of entropy.[5, 2, 3, 4] Two of these
functions are the Massieu and Planck potentials, introduced as dimensionless
forms of the Helmholtz and Gibbs energies and defined as $\Phi_{\rm
A}\equiv\beta F_{\rm A}$ and $\Xi_{\rm A}\equiv\beta G_{\rm A}$, respectively.
The other two are entropy–derived potentials: the dimensionless entropy
$\Pi_{\rm A}$ and its conjugate entropy $\Theta_{\rm A}$. Because these two do
not seem to have previously been given names, we refer here to $\Pi_{\rm A}$
as the Clausius entropy (C-entropy) after Clausius, who first introduced
entropy,[6] and to $\Theta_{\rm A}$, associated with the partition function
that manifestly includes an intensive variable, as the Boltzmann entropy
(B-entropy). In the present case, we can introduce the Massieu and Planck
potentials and the B-entropy as $\Delta\Phi_{\rm A}^{\rm qst}=-(W_{\rm
A}^{ext})^{\rm qst}$, $\Delta\Xi_{\rm A}^{\rm qst}=-(W_{\rm A}^{int})^{\rm
qst},$ and $\Delta\Theta_{\rm A}^{\rm qst}=(\tilde{Q}_{\rm A}^{ext})^{\rm
qst}$, respectively.

Thus, from Eqs. (49) and (55), we obtain the dimensionless (or entropic)
minimum work principle for the intensive and expensive work defined as being
done from the outside as

$\displaystyle\tilde{W}_{\rm A}^{int}\geq-\Delta\Xi_{\rm A}^{\rm qst}$ (56)

and

$\displaystyle\tilde{W}_{\rm A}^{ext}\geq-\Delta\Phi_{\rm A}^{\rm qst}.$ (57)

The latter inequality is a generalization of the Kelvin–Planck statement,
often used as a definition of the second law of thermodynamics, for isothermal
processes. From Eq. (50), we obtain

$\displaystyle\tilde{Q}_{\rm A}^{ext}\leq\Delta\Theta_{\rm A}^{\rm qst},$ (58)

which we call “the principle of maximum dimensionless heat generation” and
which is equivalent to “the principle of maximum B-entropy.”

Table 1:  Total differential expressions for the quasi-static (qst.) entropic potentials as functions of the intensive variables $\beta^{\rm qst}(t)$ and $E^{\rm qst}(t)$ and the extensive variables $H_{\rm A}^{\rm qst}(t)$ and $\tilde{P}_{\rm A}^{\rm qst}(t)$. Entropy has two definitions, depending on whether the work variable is intensive or extensive. Of these dimensionless entropies, the commonly used one, which we call Clausius entropy (C-entropy) and involves only extensive variables, is expressed as $\Pi_{\rm A}^{\rm qst}[H_{\rm A}^{\rm qst},P_{\rm A}^{\rm qst}]$, whereas the less widely used one, which we call Boltzmann entropy (B-entropy), is expressed as $\Theta_{\rm A}^{\rm qst}[H_{\rm A}^{\rm qst},E^{\rm qst}]$. The potentials are related by the Legendre transformations as shown. Qst. Potential | Differential form | Natural var. | Legendre transformation
---|---|---|---
Massieu | $d\Phi_{\rm A}^{\rm qst}=-H_{\rm A}^{\rm qst}d\beta^{\rm qst}-E^{\rm qst}d\tilde{P}_{\rm A}^{\rm qst}$ | $\beta^{\rm qst},\tilde{P}_{\rm A}^{\rm qst}$ | $\cdots$
Planck | $d\Xi_{\rm A}^{\rm qst}=-H_{\rm A}^{\rm qst}d\beta^{\rm qst}+\tilde{P}_{\rm A}^{\rm qst}dE^{\rm qst}$ | $\beta^{\rm qst},E^{\rm qst}$ | $\Xi_{\rm A}^{\rm qst}=\Phi_{\rm A}^{\rm qst}+E^{\rm qst}\tilde{P}_{\rm A}^{\rm qst}$
B-Entropy | $d\Theta_{\rm A}^{\rm qst}=\beta^{\rm qst}dH_{\rm A}^{\rm qst}+\tilde{P}_{\rm A}^{\rm qst}dE^{\rm qst}$ | $H_{\rm A}^{\rm qst},E^{\rm qst}$ | $\Theta_{\rm A}^{\rm qst}=\Xi_{\rm A}^{\rm qst}+\beta^{\rm qst}H^{\rm qst}_{\rm A}$
C-Entropy | $d\Pi_{\rm A}^{\rm qst}=\beta^{\rm qst}dH_{\rm A}^{\rm qst}-E^{\rm qst}d\tilde{P}_{\rm A}^{\rm qst}$ | $H_{\rm A}^{\rm qst},\tilde{P}_{\rm A}^{\rm qst}$ | $\Pi_{\rm A}^{\rm qst}=\Phi_{\rm A}^{\rm qst}+\beta^{\rm qst}H^{\rm qst}_{\rm A}$

From Eqs. (28) and (32), we find that these entropic potentials are related by
the following Legendre transformations:

$\displaystyle\Phi_{\rm A}^{\rm qst}(t)=\Xi_{\rm A}^{\rm
qst}(t)-\tilde{P}_{\rm A}^{\rm qst}(t)E^{\rm qst}(t)$ (59)

and

$\displaystyle\Theta_{\rm A}^{\rm qst}(t)=\Xi_{\rm A}^{\rm qst}(t)+\beta^{\rm
qst}(t)H_{\rm A}^{\rm qst}(t),$ (60)

respectively, where the quasi-equilibrium values of the potentials are
evaluated from Eqs. (49), (50), and (55). From these, the C-entropy can be
expressed as

$\displaystyle\Pi_{\rm A}^{\rm qst}(t)=\Phi_{\rm A}^{\rm qst}(t)+\beta^{\rm
qst}(t)H_{\rm A}^{\rm qst}(t).$ (61)

The Planck potential is convex for the inverse temperature $\beta^{\rm
qst}(t)$ and the external field $E^{\rm qst}(t)$ (see Appendix C). Thus, the
natural variables of the Planck potential are $\beta^{\rm qst}(t)$ and $E^{\rm
qst}(t)$, described as $\Xi_{\rm A}^{\rm qst}[\beta^{\rm qst}(t),\,E^{\rm
qst}(t)]$. From Eqs. (59)–(61), we can express the Massieu potential and the
dimensionless entropies as $\Phi_{\rm A}^{\rm qst}[\beta^{\rm
qst}(t),\,\tilde{P}_{\rm A}^{\rm qst}(t)]$, $\Theta_{\rm A}^{\rm qst}[H_{\rm
A}^{\rm qst}(t),\,E^{\rm qst}(t)]$, and $\Pi_{\rm A}^{\rm qst}[H_{\rm A}^{\rm
qst}(t),\,\tilde{P}_{\rm A}^{\rm qst}(t)]$. The various relationships between
the entropic potentials and intensive and extensive variables can be derived
for these potentials. [5, 4] Note that, although we chose enthalpy as the
natural variable, it is possible to choose internal energy for the description
of the dimensionless thermodynamic potentials (see Appendix D).

We summarize the dimensionless thermodynamic (or entropic) potentials as
functions of the natural variables in total differential form in Table 1.

####  IV.2.2 Helmholtz–Gibbs potentials: Energy state representation

From dimensionless thermodynamic variables, potentials are naturally expressed
in the entropic representation, while the energy state representation of
potentials is commonly used in thermodynamics. For convenience in numerical
simulations and to facilitate extension to non–equilibrium states,[29, 30, 31]
we introduce these potentials by evaluating them for the isothermal case and
constant-external field (thermostatic) cases.

For fixed $\beta(t)=\beta$ (isothermal case), the inequalities (49) and (55)
for $F_{\rm A}^{\rm qst}(t)=-(\tilde{W}_{\rm A}^{ext})^{\rm qst}/\beta$ and
$G_{\rm A}^{\rm qst}(t)=-(\tilde{W}_{\rm A}^{int})^{\rm qst}/\beta$ reduce to

$\displaystyle W_{\rm A}^{ext}\geq\Delta F_{\rm A}^{\rm qst}$ (62)

and

$\displaystyle W_{\rm A}^{int}\geq\Delta G_{\rm A}^{\rm qst},$ (63)

where we define the extensive work $W_{\rm A}^{ext}$ and intensive work
$W_{\rm A}^{int}$ by

$\displaystyle\frac{dW_{\rm A}^{ext}(t)}{dt}=-\mathrm{tr}_{\rm
A}\left\\{\hat{H}^{\prime}_{\rm A}(t)\frac{\partial\hat{\rho}_{\rm
A}(t)}{\partial t}\right\\}$ (64)

and

$\displaystyle\frac{dW_{\rm A}^{int}(t)}{dt}=\mathrm{tr}_{\rm
A}\left\\{\frac{\partial\hat{H}^{\prime}_{\rm A}(t)}{\partial
t}\hat{\rho}_{\rm A}(t)\right\\}.$ (65)

We can then evaluate the two free energies as $\Delta F_{\rm A}^{\rm
qst}=(W_{\rm A}^{ext})^{\rm qst}$ and $\Delta G_{\rm A}^{\rm qst}=(W_{\rm
A}^{int})^{\rm qst}$. The inequality (62) corresponds to the Kelvin–Planck
statement. From Eq. (32), these intensive and extensive works satisfy the
time-dependent Legendre transformation, expressed as

$\displaystyle\frac{dW_{\rm A}^{ext}(t)}{dt}=\frac{dW_{\rm
A}^{int}(t)}{dt}+P_{\rm A}(t)E(t).$ (66)

Whereas the Planck potential is convex for the external field $E^{\rm
qst}(t)$, the Gibbs energy is concave for $E^{\rm qst}(t)$, because $G_{\rm
A}^{\rm qst}(t)=-\Xi_{\rm A}^{\rm qst}(t)/\beta$. From Eq. (59), the Helmholtz
and Gibbs energies satisfy the Legendre transformation expressed as

$\displaystyle F_{\rm A}^{\rm qst}(t)=G_{\rm A}^{\rm qst}(t)+P_{\rm A}^{\rm
qst}(t)E^{\rm qst}(t).$ (67)

Table 2:  Total differential expressions for the quasi-static (qst.) thermodynamic potentials as functions of the intensive variables $T^{\rm qst}(t)$ and $E^{\rm qst}(t)$ and the extensive variables $S_{\rm A}^{\rm qst}(t)$ and ${P}_{\rm A}^{\rm qst}(t)$, which are related through Legendre transformations. Qst. potential | Differential form | Natural var. | Legendre transformation
---|---|---|---
Helmholtz | $dF_{\rm A}^{\rm qst}=-S_{\rm A}^{\rm qst}dT^{\rm qst}+E^{\rm qst}dP_{\rm A}^{\rm qst}$ | $T^{\rm qst},P_{\rm A}^{\rm qst}$ | $\cdots$
Gibbs | $dG_{\rm A}^{\rm qst}=-S_{\rm A}^{\rm qst}dT^{\rm qst}-P_{\rm A}^{\rm qst}dE^{\rm qst}$ | $T^{\rm qst},E^{\rm qst}$ | $G_{\rm A}^{\rm qst}=F_{\rm A}^{\rm qst}-E^{\rm qst}P_{\rm A}^{\rm qst}$
Internal | $dU_{\rm A}^{\rm qst}=T^{\rm qst}dS_{\rm A}^{\rm qst}+E^{\rm qst}dP_{\rm A}^{\rm qst}$ | $S_{\rm A}^{\rm qst},P_{\rm A}^{\rm qst}$ | $U_{\rm A}^{\rm qst}=F_{\rm A}^{\rm qst}+T^{\rm qst}S_{\rm A}^{\rm qst}$
Enthalpy | $dH_{\rm A}^{\rm qst}=T^{\rm qst}dS_{\rm A}^{\rm qst}-P_{\rm A}^{\rm qst}dE^{\rm qst}$ | $S_{\rm A}^{\rm qst},E^{\rm qst}$ | $H^{\rm qst}_{\rm A}=G_{\rm A}^{\rm qst}+T^{\rm qst}S_{\rm A}^{\rm qst}$

When we fix the external field $E(t)=E_{\rm fix}$, the inequality for the
dimensionless entropy reduces to

$\displaystyle\int^{t_{2}}_{t_{1}}\beta(t)\frac{dH_{\rm
A}(t)}{dt}dt\leq\Delta\Theta_{\rm A}^{\rm qst},$ (68)

where the system is in equilibrium at times $t_{1}$ and $t_{2}$ and we have
introduced the time-integral form so that $H_{\rm A}(t)$ can be treated even
when it has a singular point as a function of $t$. From Appendix C, we can
prove that the dimensionless entropy is concave as a function of enthalpy,
expressed as

$\displaystyle\frac{\partial^{2}\Theta_{\rm A}^{\rm
qst}(t)}{\partial\left(H_{\rm A}^{\rm qst}(t)\right)^{2}}<0.$ (69)

From the total differential form of the dimensionless entropy presented in
Table 1, we obtain $\partial\Theta_{\rm A}^{\rm qst}(t)/\partial H_{\rm
A}^{\rm qst}(t)=\beta^{\rm qst}(t)$. Thus, Eq. (69) reduces to

$\displaystyle\frac{\partial\beta^{\rm qst}(t)}{\partial H_{\rm A}^{\rm
qst}(t)}<0,$ (70)

which means that the heat capacity is positive, i.e., $\partial H_{\rm A}^{\rm
qst}(t)/\partial T^{\rm qst}(t)>0$, with the use of $d\beta^{\rm
qst}(t)=-{d}T^{\rm qst}(t)/k_{\rm B}[T^{\rm qst}(t)]^{2}$.

For the constant-external-field process, we also have the relation (see
Appendix E), expressed as

$\displaystyle S^{\rm qst}_{\rm A}(t)=-\frac{\partial G_{\rm A}^{\rm
qst}(t)}{\partial T^{\rm qst}(t)},$ (71)

where $S_{\rm A}^{\rm qst}(t)=k_{\rm B}\Theta_{\rm A}^{\rm qst}(t)$. Using Eq.
(71) and the principle of minimum dimensionless work, we can prove that the
Gibbs energy is convex for the temperature $T^{\rm qst}(t)$ (see Appendix F).
Thus, the natural variables of the Gibbs energy are the temperature $T^{\rm
qst}(t)$ and the external field $E^{\rm qst}(t)$, expressed as $G_{\rm A}^{\rm
qst}[T^{\rm qst}(t),E^{\rm qst}(t)]$.

From Eq. (60), we obtain the Legendre transformation for the enthalpy as

$\displaystyle H_{\rm A}^{\rm qst}(t)=G_{\rm A}^{\rm qst}(t)+T^{\rm
qst}(t)S_{\rm A}^{\rm qst}(t).$ (72)

Because the Gibbs and Helmholtz energies are related by the Legendre
transformation for $E^{\rm qst}(t)$ and $P_{\rm A}^{\rm qst}(t)$, the
Helmholtz energy is also convex for the temperature $T^{\rm qst}(t)$. Thus,
from Eq. (18), we obtain the Legendre transformation between the Helmholtz
energy and internal energy as follows:

$\displaystyle U_{\rm A}^{\rm qst}(t)=F_{\rm A}^{\rm qst}(t)+T^{\rm
qst}(t)S_{\rm A}^{\rm qst}(t).$ (73)

From Eqs. (67), (72), and (73), the thermodynamic potentials are expressed in
terms of the natural variables as $F_{\rm A}^{\rm qst}[T^{\rm qst}(t),P_{\rm
A}^{\rm qst}(t)]$, $H_{\rm A}^{\rm qst}[S_{\rm A}^{\rm qst}(t),E^{\rm
qst}(t)]$, and $U^{\rm qst}_{\rm A}[S_{\rm A}^{\rm qst}(t),P_{\rm A}^{\rm
qst}(t)]$, respectively.

We summarize the thermodynamic potentials as functions of the natural
variables in total differential form in Table 2.

Although Table 2 is similar to the well-known table of thermodynamics for the
equilibrium state, we have obtained it from a dynamical approach using a
quasi-static process and have not relied on calculation of the partition
function. This indicates that the same results can be obtained not only
theoretically but also experimentally by measuring heat and polarization in
quasi-static processes.

##  V Numerical Demonstrations

###  V.1 Reduced equations of motion for thermodynamic processes

We consider the standard Brownian model described as
$V(\hat{q})\equiv\hat{q}$, although, if necessary, we can treat nonlinear SB
coupling, which is critical for the description of molecular motion.[70, 64,
65] We then set $\mu(\hat{q})\equiv\mu\hat{q}$ and introduce the time-
dependent potential defined as

$\displaystyle U^{\prime}(\hat{q};t)=U(\hat{q})-\mu\hat{q}E(t).$ (74)

For the reduced density matrix element $\rho_{A}(q,q^{\prime};t)=\langle
q|\hat{\rho}_{A}(t)|q^{\prime}\rangle$, we introduce the Wigner distribution
function, which is the quantum analog of the classical distribution function
in phase space, defined as

$\displaystyle W_{\rm
A}(p,q)\equiv\frac{1}{{2\pi\hbar}}\int_{-\infty}^{\infty}{dx}\,e^{ipx/\hbar}\rho_{A}\left({q-\frac{x}{2},q+\frac{x}{2}}\right).$
(75)

The Wigner distribution function is a real function, in contrast to the
complex density matrix: it reduces to the classical distribution function in
the classical limit.

We choose the coefficients $\nu_{l}^{k}$ and $\eta_{l}^{k}$ in Eq. (16) for
finite $K$. Then, we incorporate this contribution using the HEOM
formalism.[38, 57] To reduce computational cost, we employ the Padé spectral
decomposition (PSD) scheme for $\eta_{l}^{k}$ and $\nu_{l}^{k}$ to enhance
computational efficiency while maintaining accuracy.[73] For time-dependent
$\beta(t)$, the decomposition constants become time-dependent functions as
$\eta_{l}(t)$ and $\nu_{l}(t)$.

Under quasi-static conditions, we assume that the time scale of the quantum
thermal fluctuations $\beta(t)\hbar/2\pi$ is shorter than that of the
subsystem. Thus, the SB coherence among different heat baths [e.g., the $k$th
bath and the $(k+1)$th bath] is negligible. Under this condition, we can
describe the system dynamics using a $K$-dimensional hierarchy instead of a
$(K\times N$)-dimensional hierarchy, where $N$ is the number of heat baths.

####  V.1.1 Low-temperature quantum Fokker–Planck equations (LT-QFPE) for
thermodynamic processes

Under the PSD scheme, the equations of motion for the Wigner function are
expressed as[57]

$\displaystyle\frac{\partial}{\partial t}{W}_{\vec{n}}(p,q;t)$
$\displaystyle=-\left(\mathcal{\hat{L}}_{qm}(t)+\sum_{l}^{K}n_{l}\nu_{l}(t)+\hat{\Xi}_{K}(p,q;t)\right){W}_{\vec{n}}(p,q;t)$
$\displaystyle\quad-\sum_{l=1}^{K}\hat{\Phi}_{p}(t){W}_{\vec{n}+\vec{e}_{l}}(p,q;t)$
$\displaystyle\quad-\sum_{l=1}^{K}n_{l}\nu_{l}(t)\hat{\Theta}_{l}(p,q;t){W}_{\vec{n}-\vec{e}_{l}}(p,q;t),$
(76)

where $\vec{n}\equiv(\dots,n_{k},\dots)$ is a $K$-dimensional multi-index
whose components are all non-negative integers and
$\vec{e}_{k}\equiv(0,\dots,1,0,\dots)$ is the $k$th unit vector. The multi-
index $\vec{n}$ represents the index of the hierarchy. Physically, the first
hierarchical element, ${W}_{\vec{0}}(p,q,t)$, corresponds to ${W}_{A}(p,q,t)$,
and the rest of the hierarchical elements serve only to facilitate the
treatment of the non-Markovian system–bath interaction that arises from the
hierarchical low-temperature expansion of the noise correlation functions in
terms of $e^{-\nu(t)t}$. In the Wigner representation, the quantum Liouvillian
takes the form[57, 74]

$\displaystyle-$
$\displaystyle\mathcal{\hat{L}}_{qm}(t)W(p,\,q)\equiv-\frac{p}{m}\frac{\partial}{{\partial
q}}W(p,\,q)$ $\displaystyle+$
$\displaystyle\sum_{n=0}^{\infty}\frac{1}{(2n+1)!}\frac{\partial^{2n+1}U^{\prime}(q;t)}{\partial
q^{2n+1}}\left(-\frac{\hbar^{2}}{4}\frac{\partial^{2}}{\partial
p^{2}}\right)^{n}\frac{\partial}{\partial p}W(p,q;t).$

The operators appearing in Eq. (76) are defined as

$\displaystyle\hat{\Phi}_{p}(t)$
$\displaystyle\equiv-\frac{A}{\beta(t)}\frac{\partial}{\partial p},$ (78)
$\displaystyle\hat{\Theta}_{0}(p,q;t)=\frac{A\beta(t)}{m}\biggl{(}p+\frac{m}{\beta(t)}\frac{\partial}{\partial
p}\biggr{)},$ (79) $\displaystyle\hat{\Theta}_{l}(p,q;t)$ $\displaystyle\equiv
2A\eta_{l}\frac{\partial}{\partial p}~{}~{}({\rm for}~{}1\leq l\leq K),$ (80)

and

$\displaystyle\hat{\Xi}_{K}(p,q;t)$
$\displaystyle\equiv\hat{\Phi}_{p}(t)\sum_{l=0}^{K}\hat{\Theta}_{l}(p,q;t).$
(81)

Owing to the presence of low-temperature correction terms, the system and bath
are entangled,[38] i.e., $\hat{\rho}_{\rm
tot}(t)\neq\hat{\rho}_{A}^{0}(t)\hat{\rho}_{\rm B}(t)$, where
$\hat{\rho}_{A}^{0}(t)={\rm tr}_{\rm B}\\{\hat{\rho}_{\rm tot}(t)\\}$.

Because we consider the case in which the time scale of the quantum thermal
fluctuations $\beta(t)\hbar/2\pi$ is shorter than the time scale of the
subsystem $1/\omega_{\rm A}$, the coherence between the subsystem and
different heat baths [e.g., the $k$th and $(k+1)$th baths] is taken into
account by the time-dependent Matsubara frequencies expressed as $\nu_{l}(t)$.

####  V.1.2 Kramers equation for thermodynamic processes

The Wigner distribution function reduces to the classical distribution
function in the limit $\hbar\to 0$, and hence we can directly compare the
quantum results with the classical results.[39] The classical limit of LT-QFPE
is the Kramers equation expressed as[55, 56, 57, 75, 76]

$\displaystyle\frac{\partial}{\partial t}{W}(p,q;t)$ $\displaystyle=$
$\displaystyle-\mathcal{\hat{L}}_{cl}(t)W(p,q;t)$ $\displaystyle+$
$\displaystyle\frac{A^{2}}{m}\frac{\partial}{\partial
p}\left(p+\frac{m}{\beta(t)}\frac{\partial}{\partial p}\right)W(p,q;t),$

where the classical Liouvillian is defined as

$\displaystyle-\mathcal{\hat{L}}_{cl}(t)W(p,q)=-\frac{p}{m}\frac{\partial}{\partial
q}W(p,q)+\frac{\partial U^{\prime}(q;t)}{\partial q}\frac{\partial}{\partial
p}W(p,q).$ (83)

The description of the Kramers equation is equivalent to that of the Langevin
equation[55]

$\displaystyle
m\ddot{q}+A^{2}\dot{q}+\frac{dU^{\prime}(q;t)}{dq}+\Omega(t)=0,$ (84)

with the Gaussian white noise defined as

$\displaystyle\left\langle{\Omega(t)}\right\rangle=0,\quad\left\langle{\Omega(t)\Omega(0)}\right\rangle=\frac{A^{2}}{\beta(t)}\delta(t).$
(85)

Note that $W(p,q;t)$ in the Kramers equation is a probability distribution
function, whereas $q$ in the Langevin equation is a sampling trajectory and
cannot be a distribution function unless the trajectories are averaged over
for noise sampling.

Although it is physically impossible to change the temperature of a heat bath
with infinite specific heat, Eqs. (76), (LABEL:heom_cl), and (84) take forms
in which $\beta$ in the equations of motion derived under the thermostatic
process has been replaced with $\beta(t)$. It should be noted, however, that
in the non-Ohmic case, the temperature changes even on the time scale on which
the noise correlations are defined. This makes it difficult to apply the
fluctuation–dissipation theorem to characterize the noise, and thus a simple
replacement $\beta\rightarrow\beta(t)$ is not allowed.

####  V.1.3 Isothermal work and constant-external-field heat

We denote the solution of the reduced density elements obtained from Eq. (76)
under any $E(t)$ by ${W}_{\vec{n}}(p,q,t)$, whereas that of the classical
distribution function obtained from Eq. (LABEL:heom_cl) is denoted by
${W}(q,p;t)$. The polarization and enthalpy at time $t$ are evaluated as

$\displaystyle P_{\rm A}(t)=\mu\left.\mathrm{tr}_{\rm
A}\\{q{W}_{\vec{0}}(p,q;t)\\}\right|_{T(t)=T}$ (86)

and

$\displaystyle H_{\rm A}(t)=\mathrm{tr}_{\rm
A}\left\\{\left[\frac{p^{2}}{2m}+U^{\prime}(q;t)\right]W_{\vec{0}}(p,q;t)\right\\},$
(87)

respectively. By replacing ${W}_{\vec{0}}(p,q,t)$ with ${W}(p,q,t)$ in the
above, we can evaluate the classical results. The dimensionless heat current
in Eq. (27) is expressed as

$\displaystyle\begin{split}\frac{dQ_{\rm A}(t)}{dt}&=-\mathrm{tr}_{\rm
A}\left\\{\left(\frac{Ap}{m}\right)^{2}W_{\vec{0}}(p,q;t)\right\\}\\\
&\quad+\frac{A^{2}}{m\beta(t)}\left(1+2\sum_{l=1}^{K}\eta_{l}\right)\\\
&\quad-\sum_{l=1}^{K}\mathrm{tr}_{\rm
A}\left\\{\frac{Ap}{m}W_{\vec{e}_{l}}(p,q;t)\right\\}\end{split}$ (88)

for the quantum case and as

$\displaystyle\frac{dQ_{\rm A}(t)}{dt}=-\mathrm{tr}_{\rm
A}\left\\{\left(\frac{Ap}{m}\right)^{2}W(p,q;t)\right\\}+\frac{A^{2}}{m\beta(t)}$
(89)

for the classical case. However, we find that the accuracy of the numerical
result obtained from Eqs. (88) and (89) is not sufficient, because the heat
current $dQ_{\rm A}(t)/dt$ does not become zero even when the system
approaches equilibrium, and the errors accumulate over a long simulation time.
Thus, we calculate the dimensionless heat current using Eqs. (86) and (87) as

$\displaystyle\frac{d\tilde{Q}_{\rm A}(t)}{dt}=\beta(t)\frac{dH_{\rm
A}(t)}{dt}+\beta(t)P_{\rm A}(t)\frac{dE(t)}{dt}.$ (90)

From Eq. (52), the bath plus SB interaction heat is defined as

$\displaystyle{\frac{dQ_{I+B}(t)}{dt}=\mathrm{tr}_{\rm
tot}\left\\{\hat{H}_{\rm I+B}(t)\frac{\partial\hat{\rho}_{\rm
tot}(t)}{\partial t}\right\\}}$ (91)

and satisfies $d(Q_{\rm A}(t)+Q_{\rm I+B}(t))/dt=0$.

###  V.2 Numerical Results

####  V.2.1 Simulation details

Figure 1:  Potential surfaces with (a) $E(t=0)=0.2$ and (b) $E(t=\tau^{\rm
qst})=0.5$ are shown as black curves. The blue and red curves represent the
ground and excitation states at their respective excitation energies. From
this figure, the bath temperature is considered to be low when $T=1/(k_{\rm
B}\beta)\approx 0.3$ and high when $T=1/(k_{\rm B}\beta)\approx 5.0$.

We perform simulations for an anharmonic potential system to demonstrate that
our theory provides a practical means of computing thermodynamic variables and
thermodynamic potentials for arbitrary systems. Thus, we consider a quartic
anharmonic potential with external interaction described by
$\mu(\hat{q})=\hat{q}$. The potential function is expressed as

$\displaystyle
U^{\prime}(\hat{q})=U_{2}\hat{q}^{2}+U_{3}\hat{q}^{3}+U_{4}\hat{q}^{4}-E(t)q,$
(92)

where the harmonic and anharmonic constants are $U_{2}=0.1$, $U_{3}=0.02$, and
$U_{4}=0.05$.

Snapshots of the potential surface with the eigenstates and eigenenergies are
depicted in Fig. 1. In the isothermal process, the first excitation energy is
$\Delta E_{\rm g\rightarrow e}(t)\sim 0.8$, and so the bath temperature is
considered low when $T=1/(k_{\rm B}\beta)\approx 0.3$ and high when
$T=1/(k_{\rm B}\beta)\approx 5.0$.

To obtain a comprehensive theory for thermodynamic potentials, not only the
external field $E(t)$ but also the temperature $T(t)$ must be time-dependent
as intensive variables. Thus, we simulate both isothermal and thermostatic
processes.

We consider an isothermal process driven by a quasi-static change in $E(t)$
that consists of (i) equilibrium, (ii) quasi-static, and (iii) equilibrium
steps, defined as

$E^{\rm qst}(t)=\left\\{\begin{array}[]{rll}{\rm(i)}&0.2&(t<0),\\\
{\rm(ii)}&0.2+0.3\;t/\tau^{\rm qst}&(0\leq t<\tau^{\rm qst}),\\\
{\rm(iii)}&0.5&(\tau^{\rm qst}\leq t).\end{array}\right.$ (93)

Here, the constant $\tau^{\rm qst}$ is the time duration parameter for the
quasi-static process, and we set $\tau^{\rm qst}=1.0\times 10^{4}$.

To calculate $S_{\rm A}^{\rm qst}(t)$, we consider the thermostatic transition
with fixed external field $E(t)=0.2$ after isothermal evolution until time
$t=0$. The time profile of $T(t)$ is set by

$T^{\rm qst}(t)=\left\\{\begin{array}[]{rll}{\rm(i)}&T_{0}&(t<0),\\\
{\rm(ii)}&(1.0+t/\tau^{\rm qst})T_{0}&(0\leq t<\tau^{\rm qst}),\\\
{\rm(iii)}&2T_{0}&(\tau^{\rm qst}\leq t),\end{array}\right.$ (94)

where $T_{0}$ is the initial temperature. We perform the simulation for three
different cases: $T_{0}=5.0$ (hot), $1.0$ (intermediate), and $0.3$ (cold).

We summarize the conditions of the numerical simulation for the above two
cases in Appendix G.

####  V.2.2 Results

Figure 2: (a) Time profiles of $E^{\rm qst}(t)$ for (i) equilibrium, (ii)
quasi-static, and (iii) equilibrium steps. (b) Calculated values of $P_{\rm
A}(t)$ in the classical case (dashed line) and the quantum case (solid line)
for different SB coupling strengths. The red, green, and blue curves represent
the strong ($A=1.5$), intermediate ($A=1.0$), and weak ($A=0.5$) SB
interaction cases, respectively. The classical results are independent of the
SB bond strength, and all overlap.

Figure 2 depicts the time profiles of $E^{\rm qst}(t)$ and the calculated
$P_{\rm A}^{\rm qst}(t)$ in the isothermal case ($T(t)=0.3$). As $E^{\rm
qst}(t)$ increases, $P_{\rm A}^{\rm qst}(t)$ also increases, because of the
conjugate relationship ${\partial G^{\rm qst}_{\rm A}(t)}/{\partial E^{\rm
qst}(t)}=P_{\rm A}^{\rm qst}(t)$. In the classical case, as can be seen from
the steady-state solution of the Kramers equation (LABEL:heom_cl), we have
$Z_{\rm A}^{\rm cl}=\mathrm{tr}_{\rm A}\\{\exp[-\beta(t)\hat{H}_{\rm
A}(t)]\\}$ and, as in the conventional thermodynamics case, the results are
independent of the SB coupling strength. In the quantum case, however, $P_{\rm
A}^{\rm qst}(t)$ becomes smaller with a smaller SB coupling strength. This
difference arises from the bathentanglement, which is described by the low-
temperature correction term in the LT-QFPE. (See also Refs. 77, 78) However,
as the SB coupling becomes stronger, the system approaches the Smoluchowski
limit, where motion is suppressed, and the difference from the classical
result becomes smaller.[57]

The change in the Gibbs energy $\Delta G^{\rm qst}_{\rm A}$ in the isothermal
process is evaluated from $(W_{\rm A}^{int})^{\rm qst}$. Then, to examine the
description of the thermodynamic relation, we compute $\Delta H_{\rm A}^{\rm
G}=-T^{2}\partial(\Delta G_{\rm A}^{\rm qst}/T)/\partial T$ and $\Delta S_{\rm
A}^{\rm G}=-\partial\Delta G_{\rm A}^{\rm qst}/\partial T$ and compare these
values with separately calculated values of $\Delta H_{\rm A}^{\rm qst}(t)$
and $\Delta S_{\rm A}^{\rm qst}(t)$ from Eqs. (18) and (58), respectively. The
results are summarized in Table 3, from which it can be seen that Eqs. (125)
and (127) are valid

Table 3:  Enthalpy and entropy changes and temperature derivative of the Gibbs energy for each interaction strength $A$. Here, $\Delta H_{\rm A}^{\rm G}$ and $\Delta S_{\rm A}^{\rm G}$ are calculated from Eqs. (125) and (127), respectively. The results in this table indicate that the relation $\Delta G_{\rm A}^{\rm qst}=\Delta H_{\rm A}^{\rm G}+T^{\rm qst}\Delta S_{\rm A}^{\rm G}$ is satisfied. $A$ | $\Delta H_{\rm A}^{\rm qst}$ | $\Delta H_{\rm A}^{\rm G}$ | $\Delta S_{\rm A}^{\rm qst}$ | $\Delta S_{\rm A}^{\rm G}$ | $\Delta G_{\rm A}^{\rm qst}$
---|---|---|---|---|---
$0.5$ | $-0.170$ | $-0.170$ | $-7.63\times 10^{-2}$ | $-7.62\times 10^{-2}$ | $-0.147$
$1.0$ | $-0.181$ | $-0.180$ | $-9.81\times 10^{-2}$ | $-9.50\times 10^{-2}$ | $-0.152$
$1.5$ | $-0.192$ | $-0.192$ | $-0.117$ | $-0.116$ | $-0.157$
Figure 3:  (a) Time profiles of $\beta^{\rm qst}(t)$ for (i) equilibrium, (ii)
quasi-static, and (iii) equilibrium steps. (b) Time profile of $S_{\rm A}^{\rm
qst}(t)$ for different temperatures in the classical case (dashed line) and
the quantum case (solid line) for different SB coupling strengths. The red,
green, and blue curves represent the cases of high ($T_{0}=5.0$), intermediate
($T_{0}=1.0$), and low ($T_{0}=0.3$) bath temperature, respectively. The
classical and quantum results approach each other at higher temperatures and
almost overlap for the intermediate- and high-temperature cases.

Figure 3 presents the time profiles of $T(t)$ and the change in entropy
$\Delta S_{\rm A}^{\rm qst}(t)$ as a function of $t$ in the constant-external-
field case $A=1.0$. Because $S_{\rm A}^{\rm qst}(t)$ satisfies the conjugate
relation ${\partial G^{\rm qst}(t)}/{\partial T^{\rm qst}(t)}=-S_{\rm A}^{\rm
qst}(t)$, the variation of $S_{\rm A}^{\rm qst}(t)$ in time is similar to that
of $P_{\rm A}^{\rm qst}(t)$, while $S_{\rm A}^{\rm qst}(t)$ does not increase
linearly with $t$ because the role of $T(t)$ in ${G^{\rm qst}_{\rm A}(t)}$ is
completely different from $E(t)$. At low temperatures, quantum
bathentanglement becomes important, and the differences between classical and
quantum results are pronounced. At high and intermediate temperatures
$T_{0}=5.0$ and $1.0$, on the other hand, the contribution from the Matsubara
frequency terms becomes smaller, and so the quantum result almost overlaps
with the classical result.

##  VI Conclusions

The virtue of thermodynamics lies in its ability to describe macroscopic
thermal phenomena resulting from complex microscopic interactions in a system-
independent manner as changes in thermodynamic potentials, which are described
as interrelated intensive and extensive variables using Legendre
transformations. This virtue should be preserved when we develop a quantum
thermodynamic theory rather than an open quantum dynamical theory, although in
either case, the theory must be specific to the SB model. Moreover, as a
first-principles argument based on the SB model, the thermodynamic laws
themselves should be derived within the framework of open quantum dynamical
theory. Thus, in the present study, we have demonstrated the following:

  1. 1.

A model consisting of multiple heat baths at different temperatures
(thermodynamic SB model) has been introduced to describe the thermostatic
process. Work and heat have then been defined as the expectation values of the
Hamiltonian system.

  2. 2.

Extensive variables, such as polarization and enthalpy, have been defined as
physical variables that can describe not only equilibrium but also
non–equilibrium processes.

  3. 3.

Dimensionless (or entropic) work and heat satisfy the time-dependent
(non–equilibrium) Legendre transformations as extensive variables and their
conjugate intensive variables.

  4. 4.

We evaluate the Massieu and Planck potentials using the minimum value of the
dimensionless extensive and intensive work, respectively, and we evaluate the
Boltzmann entropy using the maximum value of the dimensionless heat (the
principles of dimensionless minimum work and dimensionless maximum entropy).
Expressions for the entropic potentials in total differential form as
functions of natural variables are presented in Table 1.

  5. 5.

These principles are a consequence of the fact that the number of energy
states (thermodynamic weights) is maximized (or entropy is maximized) when the
distribution of states becomes a canonical ensemble for the energy given at
that instant, even if the energy of the subsystem changes with time due to
temperature changes in the heat bath.

  6. 6.

The Gibbs and Helmholtz energies, enthalpy, and internal energy can be
evaluated from entropic potentials. The expressions for these functions in
total differential form as functions of natural variables are presented in
Table 2.

  7. 7.

The first and second laws of thermodynamics also follow naturally from these
arguments.

  8. 8.

Numerically accurate simulations in both classical and quantum cases have been
performed using the reduced equations of motion developed for the
thermodynamic SB model.

  9. 9.

The differences between the classical and quantum results can be attributed to
quantum entanglement between system and bath (i.e. bathentanglement).

The above statements indicate that the model and formalism introduced in the
present study are powerful means for accurate simulation and analysis of
thermodynamic problems in both classical and quantum regimes. If the
limitations due to computational cost can be removed, it will be possible to
investigate the thermodynamic properties of any system as a subsystem, such as
molecular liquids and proteins using the Langevin formalism as classical
examples,[79] and spin–lattice[62] and superconductor systems.[63]

Finally, using the present model, the definitions of non–equilibrium intensive
and extensive variables, and the reduced equations of motion, it is possible
to study thermodynamic processes in a fully non–equilibrium regime. Thus, we
will show in a subsequent paper that the thermodynamic relationships presented
in Tables 1 and 2 can be systematically extended to the non–equilibrium regime
in terms of the time-dependent intensive variables $E(t)$ and $T(t)$ and
extensive variables $P_{\rm A}(t)$ and $S_{\rm A}(t)$ by introducing waste
heat.[80]

## Acknowledgments

Y. T. thanks to Peter Wolynes for pointing out the relationship between the
dimensionless thermodynamic potentials and the Massiue–Planck potentials. Y.T.
was supported by JSPS KAKENHI (Grant No. B21H01884). S.K. acknowledges a
fellowship supported by JST, the establishment of university fellowships
toward the creation of science technology innovation (Grant No. JPMJFS2123).
S. K. was also supported by Grant-in-Aid for JSPS Fellows (Grant No.
24KJ1373).

## Author declarations

### Conflict of Interest

The authors have no conflicts to disclose.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Dimensionless minimum work principle

To obtain the dimensionless minimum work principle for multiple-bath systems,
the derivation of the Jarzynski equality is modified.[81, 29]

We assume that the subsystem is initially connected only to the first bath
($k=1$), and the total density operator at time $t_{0}$ is expressed as

$\displaystyle\hat{\rho}_{\rm tot}^{\rm init}=\frac{1}{Z_{\rm tot}^{\rm
init}}e^{-\hat{\tilde{H}}_{\rm tot}(t_{0})},$ (95)

where $\hat{\tilde{H}}_{\rm tot}(t)$ is defined in Eq. (34) and $Z_{\rm
tot}^{\rm init}$ is the total partition function expressed as

$\displaystyle Z_{\rm tot}^{\rm init}(t_{0})=Z_{\rm
A+IB}^{1}(t_{0})\prod_{k=2}^{N}Z_{\rm B}^{k},$ (96)

with

$\displaystyle Z_{\rm
A+IB}^{k}(t)=\mathrm{tr}_{\mathrm{A+I}\mathrm{B}^{k}}\left\\{e^{-\beta_{k}\left(\hat{H}_{\rm
A}(t)+\hat{H}_{\rm I+B}^{k}(t)\right)}\right\\}$ (97)

and

$\displaystyle Z_{\rm
B}^{k}=\mathrm{tr}_{\mathrm{B}^{k}}\left\\{e^{-\beta_{k}\hat{H}_{\rm
B}^{k}}\right\\}.$ (98)

Under the assumption that only the $N$th bath as being connected to the
subsystem at a time $t$, the partition function of the final state is given by

$\displaystyle Z^{\rm fin}_{\rm tot}=Z_{\rm
A+IB}^{N}(t)\prod_{k=1}^{N-1}Z_{\rm B}^{k}.$ (99)

The time evolution of the total system from $t_{0}$ to $t$ is described by the
operator $\hat{\mathcal{G}}(t,t_{0})$. We then have

$\frac{Z_{\rm A}(t)}{Z_{\rm A}(t_{0})}=\\\ \mathrm{tr}_{\rm
tot}\left\\{\hat{\mathcal{G}}^{\dagger}_{\rm
tot}(t,t_{0})e^{-\hat{\tilde{H}}_{\rm tot}(t)}\hat{\mathcal{G}}_{\rm
tot}(t,t_{0})e^{\hat{\tilde{H}}_{\rm tot}(t_{0})}\hat{\rho}_{\rm tot}^{\rm
init}\right\\},$ (100)

where $Z_{\rm A}(t)=Z_{\rm A+IB}^{k}/Z_{\rm B}^{k}$ [for $\xi_{k}(t)\neq 0$]
is the partition function of the subsystem.

Let $|\tilde{D}_{i}\rangle$ and $|\tilde{E}_{j}\rangle$ be the eigenkets of
the operators $\hat{\tilde{H}}_{\rm tot}(t_{0})$ and $\hat{\tilde{H}}_{\rm
tot}(t)$. The right-hand side of Eq. (100) is then expressed as

$\displaystyle\sum_{i,j}e^{-(\tilde{E}_{j}-\tilde{D}_{i})}|\langle\tilde{E}_{j}|\hat{\mathcal{G}}_{\rm
tot}(t,t_{0})|\tilde{D}_{i}\rangle|^{2}p(\tilde{D}_{i}),$ (101)

where $p(\tilde{D}_{i})=\langle\tilde{D}_{i}|\exp(-\hat{\tilde{H}}_{\rm
tot}(t_{0}))|\tilde{D}_{i}\rangle/Z_{\rm tot}^{\rm init}$ is the population of
the $i$th state. Applying Jensen’s inequality to Eq. (101), we obtain

$\displaystyle\frac{e^{-\langle\hat{\tilde{H}}_{\rm
tot}(t)\rangle}}{e^{-\langle\hat{\tilde{H}}_{\rm
tot}(t_{0})\rangle}}\leq\frac{Z_{\rm A}(t)}{Z_{\rm A}(t_{0})}.$ (102)

Taking the logarithm of both sides of Eq. (102), we have the inequality

$\displaystyle\int^{t}_{t_{0}}\frac{\partial}{\partial
t^{\prime}}\mathrm{tr}_{\rm tot}\left\\{\hat{\tilde{H}}_{\rm
tot}(t^{\prime})\hat{\rho}_{\rm
tot}(t^{\prime})\right\\}dt^{\prime}\geq-\Delta\Xi_{\rm A},$ (103)

where $\Xi_{\rm A}(t)=\ln Z_{\rm A}(t)$ is the Planck potential.

When the $k$th bath is connected to and disconnected from the subsystem, we
have

$\displaystyle\mathrm{tr}_{\rm tot}\left\\{\left[\hat{H}_{\rm
A}(t)+\hat{H}_{\rm I+B}(t)\right]\frac{\partial\hat{\rho}_{\rm
tot}(t)}{\partial t}\right\\}=0$ (104)

and

$\displaystyle\mathrm{tr}_{\rm tot}\left\\{\hat{H}_{\rm
B}^{k}\frac{\partial\hat{\rho}_{\rm tot}(t)}{\partial t}\right\\}=0.$ (105)

Thus, the left-hand side of Eq. (103) reduces to

$\displaystyle\int^{t}_{t_{0}}\mathrm{tr}_{\rm
tot}\left\\{\frac{\partial\hat{\tilde{H}}_{\rm tot}(t^{\prime})}{\partial
t^{\prime}}\hat{\rho}_{\rm tot}(t^{\prime})\right\\}dt^{\prime}.$ (106)

In the case where the work due to the SB interaction does not contribute, we
can omit $\partial\hat{H}_{\rm I+B}^{k}(t)/\partial t$. Thus, we obtain the
dimensionless (or entropic) minimum work principle expressed as

$\displaystyle\int^{t}_{t_{0}}\mathrm{tr}_{\rm
A}\left\\{\frac{\partial}{\partial
t^{\prime}}\left[\beta(t^{\prime})\hat{H}_{\rm
A}(t^{\prime})\right]\hat{\rho}_{\rm
A}(t^{\prime})\right\\}\geq-\Delta\Xi_{\rm A}.$ (107)

##  Appendix B Equality of the minimum work principle

In this appendix, we consider only the quasi-static case and omit the
superscript $k$ in Eqs. (5)–(7). When $J(\omega)$ is Ohmic, the spectral
density is invariant under the transformation $\hat{H}_{\rm I+B}\rightarrow
a\hat{H}_{\rm I+B}$ with $m_{j}\rightarrow m_{j}/a$ , $\omega_{j}\rightarrow
a\omega_{j}$, and $c_{j}\rightarrow ac_{j}$, where $a>0$ is a dimensionless
scaling parameter.

Because the dynamics of the subsystem depend only on temperature, $J(\omega)$,
$\hat{H}_{\rm A}(t)$, and $\hat{\rho}_{\rm A}(t)$ in the quasi-static process
are independent of the choice of $a$. Therefore, we can set the parameter
$a=\beta_{0}/\beta(t)$ for later convenience without loss of generality, where
$\beta_{0}>0$ is the arbitrary inverse temperature.

Using Kubo’s identity, we have[29]

$\begin{split}&e^{-\beta(t+\Delta t)\hat{H}_{\rm tot}(t+\Delta
t)}-e^{-\beta(t)\hat{H}_{\rm tot}(t)}\\\ &=-\int^{1}_{0}d\lambda
e^{-(1-\lambda)\beta(t)\hat{H}_{\rm tot}(t)}\\\
&\quad\times\left\\{\beta(t+\Delta t)\hat{H}_{\rm tot}(t+\Delta
t)-\beta(t)\hat{H}_{\rm tot}(t)\right\\}\\\ &\quad\times
e^{-\lambda\beta(t+\Delta t)\hat{H}_{\rm tot}(t+\Delta t)}.\end{split}$ (108)

Expanding Eq. (108) in $\Delta t$ and taking the limit $\Delta t\rightarrow
0$, we obtain

$\begin{split}\frac{\partial}{\partial t}e^{-\beta(t)\hat{H}_{\rm
tot}(t)}={}&-\int^{1}_{0}d\lambda e^{-(1-\lambda)\beta(t)\hat{H}_{\rm
tot}(t)}\\\ &\times\frac{\partial}{\partial t}\left[\beta(t)\hat{H}_{\rm
A}(t)\right]e^{-\lambda\beta(t)\hat{H}_{\rm tot}(t)}.\end{split}$ (109)

Here, we have employed th equality

$\displaystyle\frac{\partial}{\partial t}\left[\beta(t)\hat{H}_{\rm
tot}(t)\right]=\frac{\partial}{\partial t}\left[\beta(t)\hat{H}_{\rm
A}(t)\right],$ (110)

using the fact that $\beta(t)\hat{H}_{\rm I+B}(t)$ is time-independent because
the time-dependent terms in $\beta(t)$ and $\hat{H}_{\rm I+B}(t)$ cancel out.
Taking the trace of the total system on both sides of Eq. (109), we obtain

$\frac{\partial}{\partial t}Z_{\rm A}(t)=-\mathrm{tr}_{\rm
tot}\left\\{\frac{\partial}{\partial t}\left[\beta(t)\hat{H}_{\rm
A}(t)\right]\frac{e^{-\beta(t)\hat{H}_{\rm tot}(t)}}{Z_{\rm B}}\right\\},$
(111)

where $Z_{\rm A}(t)=\mathrm{tr}_{\rm tot}\\{e^{-\beta(t)\hat{H}_{\rm
tot}(t)}\\}/Z_{\rm B}$ is the partition function of the subsystem and $Z_{\rm
B}=\mathrm{tr}_{\rm B}\\{e^{-\beta(t)\hat{H}_{\rm B}(t)}\\}$ is the bath
partition function, which is time-independent. Dividing both sides of Eq.
(111) by $Z_{\rm A}(t)$, we obtain the equality

$\frac{\partial}{\partial t}\left[\beta(t)G_{\rm A}(t)\right]=\mathrm{tr}_{\rm
A}\left\\{\frac{\partial}{\partial t}\left[\beta(t)\hat{H}_{\rm
A}(t)\right]\hat{\rho}_{\rm A}^{\rm qst}(t)\right\\},$ (112)

where we have introduced the quasi-static Gibbs energy $G_{\rm A}(t)=-\ln
Z_{\rm A}(t)/\beta(t)$ and the reduced density operator in the quasi-static
process $\hat{\rho}_{\rm A}^{\rm qst}(t)$. Solving Eq. (112) for $dG_{\rm
A}(t)/dt$, we obtain

$\displaystyle\frac{dG_{\rm A}^{\rm qst}(t)}{dt}={}$
$\displaystyle-\frac{\mathrm{tr}_{\rm A}\\{\hat{H}_{\rm A}(t)\hat{\rho}_{\rm
A}^{\rm qst}(t)\\}-G_{\rm A}^{\rm qst}(t)}{T^{\rm qst}(t)}\frac{dT^{\rm
qst}(t)}{dt}$ $\displaystyle+\mathrm{tr}_{\rm A}\\{\mu(\hat{q})\hat{\rho}_{\rm
A}^{\rm qst}(t)\\}\frac{dE^{\rm qst}(t)}{dt}.$ (113)

By comparing Eq. (B) with $dG_{\rm A}^{\rm qst}=-S_{\rm A}^{\rm qst}dT^{\rm
qst}+P^{\rm qst}_{\rm A}dE^{\rm qst}$, the entropy and polarization are
evaluated in terms of the partition function as

$S_{\rm A}^{\rm qst}(t)=k_{\rm B}\left[\beta^{\rm qst}(t)\mathrm{tr}_{\rm
A}\\{\hat{H}_{\rm A}(t)\hat{\rho}_{\rm A}^{\rm qst}(t)\\}-\ln Z_{\rm
A}(t)\right]$ (114)

and

$\displaystyle P^{\rm qst}_{\rm A}(t)=\mathrm{tr}_{\rm
A}\\{\mu(\hat{q})\hat{\rho}_{\rm A}^{\rm qst}(t)\\}.$ (115)

##  Appendix C Convexity and concavity of thermodynamic functions

Figure 4:  Schematic illustration of the convexity of $\Xi_{\rm A}^{\rm qst}$.
The blue curve represents the Planck potential as a function of the inverse
temperature $\beta$ for fixed $E=E_{\rm fix}$. The red line represents the
tangent line of the blue curve at the inverse temperature $\beta_{1}$,
expressed as $\Xi_{\rm A}^{\rm qst})_{1}-(H_{\rm A}^{\rm
qst})_{1}(\beta-\beta_{1})$. The red line is below the blue curve.

We show the convexity of the Planck potential for $\beta^{\rm qst}(t)$. For a
fixed external field ${E(t)=E_{\rm fix}}$, consider the two equilibrium states
$\alpha=1$ at time $t_{1}$ and $\alpha=2$ at time $t_{2}$. Under this
condition, the dimensionless minimum work principle is expressed as

$\displaystyle-\int^{t_{2}}_{t_{1}}H_{\rm
A}(t^{\prime})\frac{d\beta(t^{\prime})}{dt^{\prime}}dt^{\prime}\leq\left(\Xi_{\rm
A}^{\rm qst}\right)_{2}-\left(\Xi_{\rm A}^{\rm qst}\right)_{1},$ (116)

where $\left(\Xi_{\rm A}^{\rm qst}\right)_{\alpha}$ is the Planck potential
for $\alpha=1$ or 2. We then apply the inequality (116) for the thermostatic
process as follows:

$\beta(t)=\begin{cases}\beta_{1}&(t\leq t_{1}),\\\
\beta_{2}&(t>t_{1}).\end{cases}$ (117)

The left-hand side of Eq. (116) can be evaluated as $-H_{\rm
A}(t_{1})(\beta_{2}-\beta_{1})$. Because the system is in equilibrium at
$t_{1}$, the enthalpy is expressed in terms of its equilibrium value as
$H_{\rm A}(t_{1})=(H_{\rm A}^{\rm qst})_{1}$. Using $H_{\rm A}^{\rm
qst}(t)=-\partial\Xi_{\rm A}^{\rm qst}(t)/\partial\beta^{\rm qst}(t)$ for Eq.
(116), we obtain

$\displaystyle\left(\Xi_{\rm A}^{\rm
qst}\right)_{1}+\left.\frac{\partial\Xi_{\rm A}^{\rm qst}[\beta,E_{\rm
fix}]}{\partial\beta}\right|_{\beta=\beta_{1}}(\beta_{2}-\beta_{1})\leq\left(\Xi_{\rm
A}^{\rm qst}\right)_{2}.$ (118)

As shown in Fig. 4, the left-hand side of Eq. (118) corresponds to the tangent
line of the Planck potential. The inverse temperature $\beta_{1}$ and
$\beta_{2}$ are arbitrary, indicating that the Planck potential is convex with
respect to $\beta^{\rm qst}(t)$.

In the same manner, we can show the convexity and concavity of the
dimensionless entropic potentials. We summarize the results in Table 4.

Table 4:  Convexity and concavity of the entropic potentials for the natural variables. Potential | Symbol | Convexity and concavity
---|---|---
Massieu | $\Phi_{\rm A}^{\rm qst}$ |  $\beta^{\rm qst}$: convex |  $\tilde{P}_{\rm A}^{\rm qst}$: concave
Planck | $\Xi_{\rm A}^{\rm qst}$ |  $\beta^{\rm qst}$: convex |  $E^{\rm qst}$: convex
B-entropy | $\Theta_{\rm A}^{\rm qst}$ |  $H_{\rm A}^{\rm qst}$: concave |  $E^{\rm qst}$: convex
C-entropy | $\Pi_{\rm A}^{\rm qst}$ |  $H_{\rm A}^{\rm qst}$: concave |  $\tilde{P}_{\rm A}^{\rm qst}$: concave

##  Appendix D Internal energy representation

For the description of the entropic potentials, here we employ the internal
energy instead of the enthalpy. From Table 1, we have

$\displaystyle d\Xi_{\rm A}^{\rm qst}=-H_{\rm A}^{\rm qst}d\beta^{\rm
qst}+\tilde{P}_{\rm A}^{\rm qst}dE^{\rm qst}.$ (119)

Substituting Eq. (22) into this, the Planck potential is expressed as

$\displaystyle d\Xi_{\rm A}^{\rm qst}(t)=-U_{\rm A}^{\rm qst}(t)d\beta^{\rm
qst}(t)+P_{\rm A}^{\rm qst}(t)d\tilde{E}^{\rm qst}(t),$ (120)

where $\tilde{E}^{\rm qst}(t)=\beta^{\rm qst}(t)E^{\rm qst}(t)$. Applying the
Legendre transformation to this, we obtain the Massieu potential and B-entropy
in the internal energy representation as

$\displaystyle d\Phi_{\rm A}^{\rm qst}(t)=-U_{\rm A}^{\rm qst}(t)d\beta^{\rm
qst}(t)-\tilde{E}^{\rm qst}(t)dP_{\rm A}^{\rm qst}(t)$ (121)

and

$\displaystyle d\Theta_{\rm A}^{\rm qst}(t)=\beta^{\rm qst}(t)dU_{\rm A}^{\rm
qst}(t)-\tilde{E}^{\rm qst}(t)dP_{\rm A}^{\rm qst}(t),$ (122)

where we have used $E^{\rm qst}(t)\tilde{P}_{\rm A}^{\rm
qst}(t)=\tilde{E}^{\rm qst}(t)P_{\rm A}^{\rm qst}(t)$ and $\beta^{\rm
qst}(t)H_{\rm A}^{\rm qst}(t)=\beta^{\rm qst}(t)U_{\rm A}^{\rm
qst}(t)-\tilde{E}^{\rm qst}(t)P_{\rm A}^{\rm qst}(t)$. In Eq. (122), the
B-entropy is expressed in terms of extensive variables, while the value of
B-entropy is unchanged. However, in the internal energy representation, we
cannot obtain the C-entropy from the Legendre transformation of the B-entropy,
because the sign of the second term in Eq. (122) is opposite to that in Table
1.

##  Appendix E Enthalpy and Gibbs energy as functions of $T^{\rm qst}(t)$ and
$S^{\rm qst}(t)$

Here, we consider the fixed-external-field process [$E(t)=E_{\rm fix}$]. From
the total differential form of the Planck potential presented in Table 1, we
have

$\displaystyle\frac{d\Xi_{\rm A}^{\rm qst}(t)}{dt}=-H_{\rm A}^{\rm
qst}(t)\frac{d\beta^{\rm qst}(t)}{dt}.$ (123)

Because the Gibbs energy is defined as $G_{\rm A}^{\rm qst}(t)=-\Xi_{\rm
A}^{\rm qst}(t)/\beta$ for arbitrary $\beta$, we can replace $\beta$ with
$\beta^{\rm qst}(t)$. In this way, we can extend this definition of Gibbs
energy to the thermostatic process. Then, substituting $\Xi_{\rm A}^{\rm
qst}(t)={-\beta^{\rm qst}(t)}G_{\rm A}^{\rm qst}(t)$ in Eq. (123), we have

$\displaystyle\frac{d}{dt}\left[\beta^{\rm qst}(t)G_{\rm A}^{\rm
qst}(t)\right]=H_{\rm A}^{\rm qst}(t)\frac{d\beta^{\rm qst}(t)}{dt},$ (124)

which leads to

$\displaystyle H_{\rm A}^{\rm qst}(t)=\frac{\partial}{\partial\beta^{\rm
qst}(t)}\left[\beta^{\rm qst}(t)G_{\rm A}^{\rm qst}(t)\right].$ (125)

Equation (124) can then be rewritten as

$\displaystyle\frac{dG_{\rm A}(t)}{dt}=\frac{\Xi_{\rm A}^{\rm
qst}(t)+\beta^{\rm qst}(t)H_{\rm A}^{\rm qst}(t)}{(\beta^{\rm
qst}(t))^{2}}\frac{d\beta^{\rm qst}(t)}{dt},$ (126)

and thus we have

$\displaystyle S_{\rm A}(t)=k_{\rm B}(\beta^{\rm qst}(t))^{2}\frac{\partial
G_{\rm A}(t)}{\partial\beta^{\rm qst}(t)},$ (127)

where we have applied the Legendre transformation in Eq. (60) to the right-
hand side of Eq. (126).

##  Appendix F Concavity of Gibbs energy

Here, we prove the concavity of the Gibbs energy for $T^{\rm qst}(t)$ for a
fixed external field $E^{\rm qst}(t)=E_{\rm fix}$. We then consider the
transition from equilibrium state 1 to another equilibrium state 2, as in
Appendix C, where the inequality (116) holds. For the thermostatic process
defined in Eq. (117), we obtain the inequality expressed as

$\displaystyle-\left(H_{\rm A}^{\rm
qst}\right)_{1}(\beta_{2}-\beta_{1})\leq\left(\Xi_{\rm A}^{\rm
qst}\right)_{2}-\left(\Xi_{\rm A}^{\rm qst}\right)_{1}.$ (128)

Then, using the equality $(H_{\rm A}^{\rm qst})_{1}=((\Theta_{\rm A}^{\rm
qst})_{1}-(\Xi_{\rm A}^{\rm qst})_{1})/\beta_{1}$, which we obtain from the
Legendre transformation for the dimensionless B-entropy in the quasi-static
(equilibrium) state 1 (i.e., $(\Theta_{\rm A}^{\rm qst})_{1}=\beta_{1}(H_{\rm
A}^{\rm qst})_{1}+(\Xi_{\rm A}^{\rm qst})_{1}$ ), we obtain

$\displaystyle\left(\Theta_{\rm A}^{\rm
qst}\right)_{1}\left(\frac{\beta_{2}}{\beta_{1}}-1\right)\leq\left(\Xi_{\rm
A}^{\rm qst}\right)_{2}-\frac{\beta_{2}}{\beta_{1}}\left(\Xi_{\rm A}^{\rm
qst}\right)_{1}.$ (129)

Dividing both sides of Eq. (129) by $-\beta_{2}$, we have the inequality for
the Gibbs energy as

$\displaystyle-\left(S_{\rm A}^{\rm
qst}\right)_{1}(T_{2}-T_{1})\geq\left(G_{\rm A}^{\rm
qst}\right)_{2}-\left(G_{\rm A}^{\rm qst}\right)_{1},$ (130)

where we have used the definitions of the quasi-static entropy $(S_{\rm
A}^{\rm qst})_{1}=k_{\rm B}(\Theta_{\rm A}^{\rm qst})_{1}$ and of the quasi-
static Gibbs energy $(G_{\rm A}^{\rm qst})_{\alpha}=-(\Xi_{\rm A}^{\rm
qst})_{\alpha}/\beta_{\alpha}\;(\alpha=1,2)$. Because we have Eq. (71), Eq.
(130) leads to the concavity of the Gibbs energy for the temperature.

We summarize the convexity and concavity of the thermodynamic potentials in
Table 5.

Table 5:  Convexity and concavity of the thermodynamic potentials for the natural variables. Potential | Symbol | Convexity and concavity
---|---|---
Helmholtz | $F_{\rm A}^{\rm qst}$ |  $T^{\rm qst}$: concave |  $\tilde{P}_{\rm A}^{\rm qst}$: convex
Gibbs | $G_{\rm A}^{\rm qst}$ |  $T^{\rm qst}$: concave |  $E^{\rm qst}$: concave
Internal energy | $U_{\rm A}^{\rm qst}$ |  $S_{\rm A}^{\rm qst}$: convex |  $\tilde{P}_{\rm A}^{\rm qst}$: convex
Enthalpy | $H_{\rm A}^{\rm qst}$ |  $S_{\rm A}^{\rm qst}$: convex |  $E^{\rm qst}$: concave

##  Appendix G Numerical simulation details

The numerical calculations in Sec. V.2 were carried out to integrate Eq. (76)
with Eqs. (V.1.1)–(81) in the quantum cases and Eq. (LABEL:heom_cl) with Eq.
(83) in the classical cases, using a fourth-order Runge–Kutta algorithm with
the predictor–corrector method incorporated into the fourth-order
Adams–Bashforth method and the fourth-order Adams–Moulton method.[82] We set
the time step $\delta t=0.001$. Uniform meshes were employed to discretize the
Wigner function with mesh sizes of $N_{q}=64$ and $N_{p}=64$ in the $q$ and
$p$ directions, respectively.

Other parameters for the isothermal and thermostatic processes are listed in
Tables 6 and 7, respectively.

Table 6:  Parameter values used for the simulations of the isothermal process. Here, $dx$ and $dp$ are the mesh sizes for position and momentum, respectively, in Wigner space. The integers $N$ and $K$ are the cutoff numbers used in the LT-QFPE. | $A$ | $N$ | $K$ | $dx$ | $dp$
---|---|---|---|---|---
Classical | $0.5$ | $\cdots$ | $\cdots$ | $0.2$ | $0.2$
$1.0$ | $\cdots$ | $\cdots$ | $0.2$ | $0.2$
$1.5$ | $\cdots$ | $\cdots$ | $0.2$ | $0.2$
Quantum | $0.5$ | $6$ | $3$ | $0.2$ | $0.2$
$1.0$ | $7$ | $3$ | $0.2$ | $0.2$
$1.5$ | $8$ | $3$ | $0.25$ | $0.35$
Table 7:  Parameter values used for the simulations of the thermostatic process. Here, $dx$ and $dp$ are the mesh sizes for position and momentum, respectively, in Wigner space. The integers $N$ and $K$ are the cutoff numbers used in the LT-QFPE. | $T_{0}$ | $N$ | $K$ | $dx$ | $dp$
---|---|---|---|---|---
Classical | $0.3$ | $\cdots$ | $\cdots$ | $0.2$ | $0.2$
$1.0$ | $\cdots$ | $\cdots$ | $0.25$ | $0.25$
$5.0$ | $\cdots$ | $\cdots$ | $0.3$ | $0.4$
Quantum | $0.3$ | $7$ | $3$ | $0.25$ | $0.3$
$1.0$ | $7$ | $2$ | $0.3$ | $0.45$
${5.0}$ | $7$ | $1$ | $0.3$ | $0.6$

## References

  * Carnot [1824] S. Carnot,  _Reflexions sur la puissance motrice du feu_ (ChexzBachelier Libraire, Quai des Augustins, Paris, 1824).
  * Massieu [1869] F. Massieu, “Sur les fonctions caractéristiques des divers fluides,” CR Acad. Sci. Paris 69, 858–862 (1869).
  * Planck [1922] M. Planck,  _Vorlesungen uber Thermodynamik_ (De Gruyter, 1922).
  * Planes and Vives [2002] A. Planes and E. Vives, “Entropic Formulation of Statistical Mechanics,” Journal of Statistical Physics 106, 827–850 (2002).
  * Guggenheim [1986] E. A. Guggenheim,  _Thermodynamics: An Advanced Treatment for Chemists and Physicists_ (North Holland, 1986).
  * Oono [2017] Y. Oono,  _Perspectives on Statistical Thermodynamics_ (Cambridge University Press, 2017).
  * Darwin and Fowler [1922] C. Darwin and R. Fowler, “Xliv. on the partition of energy,” The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 44, 450–479 (1922).
  * Quan  _et al._ [2006] H. T. Quan, Y. D. Wang, Y.-x. Liu, C. P. Sun, and F. Nori, “Maxwell’s demon assisted thermodynamic cycle in superconducting quantum circuits,” Phys. Rev. Lett. 97, 180402 (2006).
  * Esposito, Harbola, and Mukamel [2009] M. Esposito, U. Harbola, and S. Mukamel, “Nonequilibrium fluctuations, fluctuation theorems, and counting statistics in quantum systems,” Rev. Mod. Phys. 81, 1665–1702 (2009).
  * Campisi, Hänggi, and Talkner [2011] M. Campisi, P. Hänggi, and P. Talkner, “Colloquium: Quantum fluctuation relations: Foundations and applications,” Rev. Mod. Phys. 83, 771–791 (2011).
  * Kosloff and Levy [2014] R. Kosloff and A. Levy, “Quantum heat engines and refrigerators: Continuous devices,” Annual Review of Physical Chemistry 65, 365–393 (2014), pMID: 24689798.
  * Esposito, Ochoa, and Galperin [2015a] M. Esposito, M. A. Ochoa, and M. Galperin, “Nature of heat in strongly coupled open quantum systems,” Phys. Rev. B 92, 235440 (2015a).
  * Schmidt  _et al._ [2015] R. Schmidt, M. F. Carusela, J. P. Pekola, S. Suomela, and J. Ankerhold, “Work and heat for two-level systems in dissipative environments: Strong driving and non-Markovian dynamics,” Phys. Rev. B 91, 224303 (2015).
  * Uzdin, Levy, and Kosloff [2015] R. Uzdin, A. Levy, and R. Kosloff, “Equivalence of quantum heat machines, and quantum-thermodynamic signatures,” Phys. Rev. X 5, 031044 (2015).
  * Whitney [2018] R. S. Whitney, “Non-Markovian quantum thermodynamics: Laws and fluctuation theorems,” Phys. Rev. B 98, 085415 (2018).
  * Quan  _et al._ [2007] H. T. Quan, Y.-x. Liu, C. P. Sun, and F. Nori, “Quantum thermodynamic cycles and quantum heat engines,” Phys. Rev. E 76, 031105 (2007).
  * Maruyama, Nori, and Vedral [2009] K. Maruyama, F. Nori, and V. Vedral, “Colloquium: The physics of maxwell’s demon and information,” Rev. Mod. Phys. 81, 1–23 (2009).
  * Ono  _et al._ [2020] K. Ono, S. N. Shevchenko, T. Mori, S. Moriyama, and F. Nori, “Analog of a quantum heat engine using a single-spin qubit,” Phys. Rev. Lett. 125, 166802 (2020).
  * Talkner and Hänggi [2020a] P. Talkner and P. Hänggi, “Colloquium: Statistical mechanics and thermodynamics at strong coupling: Quantum and classical,” Rev. Mod. Phys. 92, 041002 (2020a).
  * Cangemi  _et al._ [2020] L. M. Cangemi, V. Cataudella, G. Benenti, M. Sassetti, and G. De Filippis, “Violation of thermodynamics uncertainty relations in a periodically driven work-to-work converter from weak to strong dissipation,” Phys. Rev. B 102, 165418 (2020).
  * Strasberg and Winter [2021] P. Strasberg and A. Winter, “First and second law of quantum thermodynamics: A consistent derivation based on a microscopic definition of entropy,” PRX Quantum 2, 030202 (2021).
  * Cockrell and Ford [2022] C. Cockrell and I. J. Ford, “Stochastic thermodynamics in a non-markovian dynamical system,” Phys. Rev. E 105, 064124 (2022).
  * Talkner and Hänggi [2020b] P. Talkner and P. Hänggi, “Colloquium: Statistical mechanics and thermodynamics at strong coupling: Quantum and classical,” Rev. Mod. Phys. 92, 041002 (2020b).
  * Dann and Kosloff [2023] R. Dann and R. Kosloff, “Unification of the first law of quantum thermodynamics,” New Journal of Physics 25, 043019 (2023).
  * Ferreri  _et al._ [2023] A. Ferreri, V. Macrì, F. K. Wilhelm, F. Nori, and D. E. Bruschi, “Quantum field heat engine powered by phonon-photon interactions,” Phys. Rev. Res. 5, 043274 (2023).
  * Sakamoto and Tanimura [2021] S. Sakamoto and Y. Tanimura, “Open quantum dynamics theory for non-equilibrium work: Hierarchical equations of motion approach,” Journal of the Physical Society of Japan 90, 033001 (2021).
  * Kato and Tanimura [2015] A. Kato and Y. Tanimura, “Quantum heat transport of a two-qubit system: Interplay between system-bath coherence and qubit-qubit coherence,” The Journal of Chemical Physics 143, 064107 (2015).
  * Kato and Tanimura [2016] A. Kato and Y. Tanimura, “Quantum heat current under non-perturbative and non-Markovian conditions: Applications to heat machines,” The Journal of Chemical Physics 145, 224105 (2016).
  * Sakamoto and Tanimura [2020] S. Sakamoto and Y. Tanimura, “Numerically ”exact” simulations of entropy production in the fully quantum regime: Boltzmann entropy vs von Neumann entropy,” The Journal of Chemical Physics 153, 234107 (2020).
  * Koyanagi and Tanimura [2022a] S. Koyanagi and Y. Tanimura, “The laws of thermodynamics for quantum dissipative systems: A quasi-equilibrium Helmholtz energy approach,” The Journal of Chemical Physics 157, 014104 (2022a).
  * Koyanagi and Tanimura [2022b] S. Koyanagi and Y. Tanimura, “Numerically “exact” simulations of a quantum carnot cycle: Analysis using thermodynamic work diagrams,” The Journal of Chemical Physics 157, 084110 (2022b).
  * Ullersma [1966a] P. Ullersma, “An exactly solvable model for brownian motion: I. derivation of the langevin equation,” Physica 32, 27–55 (1966a).
  * Ullersma [1966b] P. Ullersma, “An exactly solvable model for brownian motion: Ii. derivation of the fokker-planck equation and the master equation,” Physica 32, 56–73 (1966b).
  * Caldeira and Leggett [1983] A. Caldeira and A. Leggett, “Path integral approach to quantum brownian motion,” Physica A: Statistical Mechanics and its Applications 121, 587–616 (1983).
  * Grabert, Schramm, and Ingold [1988] H. Grabert, P. Schramm, and G.-L. Ingold, “Quantum brownian motion: The functional integral approach,” Physics Reports 168, 115–207 (1988).
  * Weiss [2012] U. Weiss,  _Quantum Dissipative Systems_, 4th ed. (WORLD SCIENTIFIC, 2012).
  * Tanimura [2006] Y. Tanimura, “Stochastic Liouville, Langevin, Fokker-Planck, and master equation qpproaches to quantum dissipative systems,” Journal of the Physical Society of Japan 75, 082001 (2006).
  * Tanimura [2020] Y. Tanimura, “Numerically ”exact” approach to open quantum dynamics: The hierarchical equations of motion (HEOM),” The Journal of Chemical Physics 153, 020901 (2020).
  * Tanimura [2015] Y. Tanimura, “Real-time and imaginary-time quantum hierarchal Fokker-Planck equations,” The Journal of Chemical Physics 142, 144110 (2015).
  * Kato and Tanimura [2013] A. Kato and Y. Tanimura, “Quantum suppression of ratchet rectification in a Brownian system driven by a biharmonic force,” The Journal of Physical Chemistry B 117, 13132–13144 (2013).
  * Tanimura and Kubo [1989] Y. Tanimura and R. Kubo, “Time evolution of a quantum system in contact with a nearly Gaussian-Markoffian noise bath,” Journal of the Physical Society of Japan 58, 101–114 (1989).
  * Velizhanin, Wang, and Thoss [2008] K. A. Velizhanin, H. Wang, and M. Thoss, “Heat transport through model molecular junctions: A multilayer multiconfiguration time-dependent Hartree approach,” Chemical Physics Letters 460, 325–330 (2008).
  * Wang and Thoss [2008] H. Wang and M. Thoss, “From coherent motion to localization: dynamics of the spin-boson model at zero temperature,” New Journal of Physics 10, 115005 (2008).
  * Chin  _et al._ [2011] A. W. Chin, J. Prior, S. F. Huelga, and M. B. Plenio, “Generalized polaron ansatz for the ground state of the sub-ohmic spin-boson model: An analytic theory of the localization transition,” Phys. Rev. Lett. 107, 160601 (2011).
  * Leggett  _et al._ [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, “Dynamics of the dissipative two-state system,” Rev. Mod. Phys. 59, 1–85 (1987).
  * Bulla, Tong, and Vojta [2003] R. Bulla, N.-H. Tong, and M. Vojta, “Numerical renormalization group for bosonic systems and application to the sub-ohmic spin-boson model,” Phys. Rev. Lett. 91, 170601 (2003).
  * Duan  _et al._ [2017] C. Duan, Z. Tang, J. Cao, and J. Wu, “Zero-temperature localization in a sub-ohmic spin-boson model investigated by an extended hierarchy equation of motion,” Phys. Rev. B 95, 214308 (2017).
  * Anto-Sztrikacs, Nazir, and Segal [2023] N. Anto-Sztrikacs, A. Nazir, and D. Segal, “Effective-hamiltonian theory of open quantum systems at strong coupling,” PRX Quantum 4, 020307 (2023).
  * Ma  _et al._ [2012] J. Ma, Z. Sun, X. Wang, and F. Nori, “Entanglement dynamics of two qubits in a common bath,” Phys. Rev. A 85, 062323 (2012).
  * Lambert  _et al._ [2019] N. Lambert, S. Ahmed, M. Cirio, and F. Nori, “Modelling the ultra-strongly coupled spin-boson model with unphysical modes,” Nature communications 10, 1–9 (2019).
  * Wiedmann, Stockburger, and Ankerhold [2021] M. Wiedmann, J. T. Stockburger, and J. Ankerhold, “Non-Markovian quantum Otto refrigerator,” Eur. Phys. J. Spec. Tops 230, 851–857 (2021).
  * Esposito, Ochoa, and Galperin [2015b] M. Esposito, M. A. Ochoa, and M. Galperin, “Quantum thermodynamics: A nonequilibrium Green’s function approach,” Phys. Rev. Lett. 114, 080602 (2015b).
  * Ishizaki and Tanimura [2005] A. Ishizaki and Y. Tanimura, “Quantum dynamics of system strongly coupled to low-temperature colored noise bath: Reduced hierarchy equations approach,” Journal of the Physical Society of Japan 74, 3131–3134 (2005).
  * Tanimura [2014] Y. Tanimura, “Reduced hierarchical equations of motion in real and imaginary time: Correlated initial states and thermodynamic quantities,” The Journal of Chemical Physics 141, 044114 (2014).
  * Tanimura and Wolynes [1991] Y. Tanimura and P. G. Wolynes, “Quantum and classical Fokker-Planck equations for a Gaussian-Markovian noise bath,” Phys. Rev. A 43, 4131–4142 (1991).
  * Tanimura and Wolynes [1992] Y. Tanimura and P. G. Wolynes, “The interplay of tunneling, resonance, and dissipation in quantum barrier crossing: A numerical study,” The Journal of Chemical Physics 96, 8485–8496 (1992).
  * Ikeda and Tanimura [2019] T. Ikeda and Y. Tanimura, “Low-temperature quantum Fokker-Planck and Smoluchowski equations and their extension to multistate systems,” Journal of Chemical Theory and Computation 15, 2517–2534 (2019).
  * Song and Shi [2017] L. Song and Q. Shi, “Hierarchical equations of motion method applied to nonequilibrium heat transport in model molecular junctions: Transient heat current and high-order moments of the current operator,” Phys. Rev. B 95, 064308 (2017).
  * Gong  _et al._ [2022] H. Gong, Y. Wang, X. Zheng, R. Xu, and Y. Yan, “Nonequilibrium work distributions in quantum impurity system–bath mixing processes,” The Journal of Chemical Physics 157, 054109 (2022).
  * Latune, Pleasance, and Petruccione [2023] C. L. Latune, G. Pleasance, and F. Petruccione, “Cyclic quantum engines enhanced by strong bath coupling,” Phys. Rev. Appl. 20, 024038 (2023).
  * Boettcher  _et al._ [2024] V. Boettcher, R. Hartmann, K. Beyer, and W. T. Strunz, “Dynamics of a strongly coupled quantum heat engine—Computing bath observables from the hierarchy of pure states,” The Journal of Chemical Physics 160, 094108 (2024).
  * Nakamura and Tanimura [2018] K. Nakamura and Y. Tanimura, “Hierarchical Schrödinger equations of motion for open quantum dynamics,” Phys. Rev. A 98, 012109 (2018).
  * Nakamura and Tanimura [2021] K. Nakamura and Y. Tanimura, “Optical response of laser-driven charge-transfer complex described by Holstein-Hubbard model coupled to heat baths: Hierarchical equations of motion approach,” The Journal of Chemical Physics 155, 064106 (2021).
  * Takahashi and Tanimura [2023a] H. Takahashi and Y. Tanimura, “Discretized hierarchal equations of motion in mixed liouville-wigner space for two-dimensional vibrational spectroscopies of water,” The Journal of Chemical Physics 158, 044115 (2023a).
  * Takahashi and Tanimura [2023b] H. Takahashi and Y. Tanimura, “Simulating two-dimensional correlation spectroscopies with third-order infrared and fifth-order infrared–Raman processes of liquid water,” The Journal of Chemical Physics 158, 124108 (2023b).
  * Chang  _et al._ [2006] C. W. Chang, A. M. Fennimore, A. Afanasiev, D. Okawa, T. Ikuno, H. Garcia, D. Li, A. Majumdar, and A. Zettl, “Isotope effect on the thermal conductivity of boron nitride nanotubes,” Phys. Rev. Lett. 97, 085901 (2006).
  * Knowles, Koch, and Shelton [2018] K. E. Knowles, M. D. Koch, and J. L. Shelton, “Three applications of ultrafast transient absorption spectroscopy of semiconductor thin films: spectroelectrochemistry, microscopy, and identification of thermal contributions,” J. Mater. Chem. C 6, 11853–11867 (2018).
  * Albert-Seifried and Friend [2011] S. Albert-Seifried and R. H. Friend, “Measurement of thermal modulation of optical absorption in pump-probe spectroscopy of semiconducting polymers,” Applied Physics Letters 98, 223304 (2011).
  * Rao  _et al._ [2011] A. Rao, M. W. B. Wilson, S. Albert-Seifried, R. Di Pietro, and R. H. Friend, “Photophysics of pentacene thin films: The role of exciton fission and heating effects,” Phys. Rev. B 84, 195411 (2011).
  * Okumura and Tanimura [1997] K. Okumura and Y. Tanimura, “Two-time correlation functions of a harmonic system nonbilinearly coupled to a heat bath: Spontaneous Raman spectroscopy,” Phys. Rev. E 56, 2747–2750 (1997).
  * Pechukas, Ankerhold, and Grabert [2000] P. Pechukas, J. Ankerhold, and H. Grabert, “Quantum smoluchowski equation,” Annalen der Physik 512, 794–803 (2000).
  * Ankerhold, Pechukas, and Grabert [2001] J. Ankerhold, P. Pechukas, and H. Grabert, “Strong friction limit in quantum mechanics: The quantum smoluchowski equation,” Phys. Rev. Lett. 87, 086802 (2001).
  * Hu, Xu, and Yan [2010] J. Hu, R.-X. Xu, and Y. Yan, “Communication: Padé spectrum decomposition of Fermi function and Bose function,” The Journal of Chemical Physics 133, 101106 (2010).
  * Frensley [1990] W. R. Frensley, “Boundary conditions for open quantum systems driven far from equilibrium,” Rev. Mod. Phys. 62, 745–791 (1990).
  * Kramers [1940] H. Kramers, “Brownian motion in a field of force and the diffusion model of chemical reactions,” Physica 7, 284–304 (1940).
  * Risken and Risken [1996] H. Risken and H. Risken,  _Fokker-planck equation_ (Springer, 1996).
  * Lee, Cao, and Gong [2012] C. K. Lee, J. Cao, and J. Gong, “Noncanonical statistics of a spin-boson model: Theory and exact monte carlo simulations,” Phys. Rev. E 86, 021109 (2012).
  * Xu and Jianshu [2016] D. Xu and C. Jianshu, “Non-canonical distribution and non-equilibrium transport beyond weak system-bath coupling regime: A polaron transformation approach,” Frontiers of Physics 11, 110308 (2016).
  * Ikeda, Ito, and Tanimura [2015] T. Ikeda, H. Ito, and Y. Tanimura, “Analysis of 2D THz-Raman spectroscopy using a non-Markovian Brownian oscillator model with nonlinear system-bath interactions,” The Journal of Chemical Physics 142, 212421 (2015).
  * Koyanagi and Tanimura [2024] S. Koyanagi and Y. Tanimura, “Classical and quantum thermodynamics in non–equilibrium regime: Application to stirling engine,” (2024), xxxxxx:xxxxxx .
  * Tasaki [2000] H. Tasaki, “Jarzynski relations for quantum systems and some applications,” (2000), arXiv:cond-mat/0009244 [cond-mat.stat-mech] .
  * Press  _et al._ [1988] W. H. Press, W. T. Vetterling, S. A. Teukolsky, and B. P. Flannery,  _Numerical recipes_ (Cambridge University Press, 1988).