# Classical and quantum thermodynamics in a non–equilibrium regime:
Application to Stirling engine

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 develop a thermodynamic theory in the non–equilibrium regime by extending a
theory based on the dimensionless (DL) minimum work principle previously
developed for a thermodynamic system–bath model [S. Koyanagi and Y. Tanimura,
J. Chem. Phys. in press (http://arxiv.org/abs/2405.16787)]. Our results are
described by non–equilibrium thermodynamic potentials expressed in time-
derivative form in terms of extensive and intensive variables. This is made
possible through the incorporation of waste heat, which is equivalent to the
loss work consumed by the bath, into the definitions of thermodynamic
potentials in a non–equilibrium regime. These potentials can be evaluated from
the DL non–equilibrium-to-equilibrium minimum work principle, which is derived
from the principle of DL minimum work and is equivalent to the second law of
thermodynamics. We thus obtain the non–equilibrium Massieu–Planck potentials
as entropic potentials and the non–equilibrium Helmholtz–Gibbs potentials as
free energies. Our results are verified numerically by simulating a Stirling
engine consisting of two isothermal and two thermostatic processes using the
quantum Fokker–Planck equations and the classical Kramers equation derived
from the thermodynamic system–bath model. We then show that any thermodynamic
process can be analyzed using a non–equilibrium work diagram analogous to the
equilibrium one for given time-dependent intensive variables. The results can
be used to develop efficient heat machines in non–equilibrium regimes.

##  I Introduction

Ever since Carnot explored the efficiency of heat engines 200 years ago,Carnot
(1824) there have been longstanding attempts to study thermodynamics in
non–equilibrium regimes, driven by both academic curiosity and practical
interest. In particular, recent advances in nanotechnology have led to
increased interest in the study of the thermodynamics of microscopic quantum
systems. Quan  _et al._ (2006); Esposito, Harbola, and Mukamel (2009);
Campisi, Hänggi, and Talkner (2011); Kosloff and Levy (2014); Esposito, Ochoa,
and Galperin (2015); Schmidt  _et al._ (2015); Uzdin, Levy, and Kosloff
(2015); Whitney (2018); Quan  _et al._ (2007); Maruyama, Nori, and Vedral
(2009); Ono  _et al._ (2020); Maruyama, Nori, and Vedral (2009); Talkner and
Hänggi (2020a); Cangemi  _et al._ (2020); Strasberg and Winter (2021);
Cockrell and Ford (2022); Talkner and Hänggi (2020b); Dann and Kosloff (2023);
Ferreri  _et al._ (2023); Sakamoto and Tanimura (2021); Kato and Tanimura
(2015, 2016)

As a kinetic theory, such phenomena are typically explained using open quantum
dynamics theories based on the system–bath (SB) model.Leggett  _et al._
(1987); Caldeira and Leggett (1983); Grabert, Schramm, and Ingold (1988);
Weiss (2012) Although thermodynamics is a system-independent theory, it is
derived from the presence of a heat bath and is consistent with theories based
on the SB model.Talkner and Hänggi (2020b) The main (subsystem) system of the
SB model can describe a variety of systems, from a simple two-level system to
complex molecular systems with many degrees of freedom, and it is also
possible to take the classical limit in the case where the system can be
described in phase space.Caldeira and Leggett (1983); Grabert, Schramm, and
Ingold (1988); Weiss (2012); Tanimura and Kubo (1989); Tanimura (1990, 2014,
2006, 2020) Under the condition of a heat bath with infinite degrees of
freedom, the total system irreversibly relaxes toward its equilibrium state in
time, and so the second law of thermodynamics is naturally obeyed without the
need for the assumption of ergodicity as a dynamical system.Tanimura (2006,
2020, 2014, 2015)

However, to investigate whether the thermodynamic laws are satisfied even in
quantum cases where the subsystem and bath are quantum mechanically entangled,
it is necessary to treat the bath in a non-Markovian and nonperturbative
manner so that the thermal equilibrium state of the total systems satisfies
the energy conservation law, including the SB interaction. Therefore,
equations of motion derived using the Markovian or rotational wave (or
secular) approximation, such as the Lindblad equation and the quantum master
equation, can only be applied to high-temperature regions where the subsystem
exhibits semiclassical dynamics.Tanimura (2006, 2020, 2014, 2015) Over the
past 30 years, several methodologies have been developed to accurately
describe the effects of quantum entanglement for condensed systems while
maintaining strict energy conservation to satisfy the first law of
thermodynamics. Such methods include the hierarchical equations of motion
(HEOM), Tanimura and Kubo (1989); Tanimura (1990); Ishizaki and Tanimura
(2005); Tanimura (2006, 2020, 2014, 2015) the quasi-adiabatic path integral
(QUAPI),Makri (1995); Makri and Makarov (1995a, b); Thorwart, Reimann, and
Hänggi (2000); Jadhao and Makri (2008); Makri (2014); Segal, Millis, and
Reichman (2010) and multiconfigurational time-dependent Hartree (MCTDH).Meyer,
Manthe, and Cederbaum (1990); Manthe, Meyer, and Cederbaum (1992); Wang and
Thoss (2007) Among these, the HEOM approach is ideal for thermodynamic
investigations, not only because it can perform numerically “exact” dynamic
simulations, but also because it can evaluate the change in energy of the heat
bath and the SB interaction separately, even for processes far from
equilibrium.Kato and Tanimura (2015, 2016); Sakamoto and Tanimura (2020,
2021); Song and Shi (2017); Gong  _et al._ (2022); Latune, Pleasance, and
Petruccione (2023); Boettcher  _et al._ (2024); Koyanagi and Tanimura (2022a,
b, 2024a)

While the difficulty of numerical simulation has been solved, there is a
fundamental difference between open quantum dynamics theory, which is based on
a first-principles description of kinetic systems from a microscopic
perspective, and thermodynamic theory, which is based on a phenomenological
description of thermal systems from a macroscopic perspective. For example, in
quantum mechanics, observables are defined as expectation values, whereas in
thermodynamics, they are described by macroscopic intensive and extensive
variables. Furthermore, quantum mechanics is a formalism for time evolution,
whereas thermodynamics deals primarily with static and quasi-static states
near thermal equilibrium. In this respect, most studies of quantum
thermodynamics are merely kinetic simulations of open quantum dynamics, and
their relation to thermodynamics has not been studied in depth. Thus, the
fundamental difference between microscopic quantum mechanics and macroscopic
thermodynamics raises many open questions, such as whether the Carnot limit
can be violated in a quantum case or the existence of Maxwell’s demon.

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 described as
interrelated intensive and extensive variables through 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.

It has been found that the Gibbs energy can be obtained directly from kinetic
simulations or experiments based on the minimum work principle (or
Kelvin–Planck statement), expressed as $W(t)\geq\Delta G(t)$, where $W(t)$ is
the work done from the outside to the subsystem by external fields and $\Delta
G(t)$ is the change in free energy, by evaluating the work in a quasi-static
process,Sakamoto and Tanimura (2020, 2021) which allows us to draw a work
diagram corresponding to the $P$–$V$ diagram.Koyanagi and Tanimura (2022a, b)
However, with this approach, the contributions of temperature $T$ and entropy
$S$ in the thermodynamic potential expressed as $TdS$ and $SdT$ cannot be
evaluated, because the minimum work principle is defined as an isothermal
process $dT=0$.

To overcome this limitation, we developed a thermostatic SB model that is
defined by a system coupled to multiple heat baths with different
temperatures.Koyanagi and Tanimura (2024a) We then extended the minimum work
principle to thermostatic processes in a dimensionless (DL) form (the DL
minimum work principle) as ${\tilde{W}}(t)\geq\Delta\Xi(t)$, where
${\tilde{W}}(t)\equiv\beta(t)W(t)$ is the DL (entropic) work, with
$\beta(t)\equiv 1/k_{\rm B}T(t)$ (where $k_{\rm B}$ is the Boltzmann constant)
being the time-dependent inverse temperature, and $\Delta\Xi$ is the change in
the DL Planck potential.Koyanagi and Tanimura (2024a) Not only intensive
variables but also extensive variables, which are related by the time-
dependent Legendre transformations, were introduced as quantum expectation
values of the SB system. The validity of these results has been numerically
verified using the low-temperature quantum Fokker–Planck equations in the
quantum case and the Kramers equation in the classical case,Ikeda and Tanimura
(2019) both developed for the thermodynamic SB model.Koyanagi and Tanimura
(2024a) Although these results were restricted to quasi-static cases, the
definitions of the thermodynamic potentials, intensive and extensive
variables, and Legendre transformations were defined in such a way that they
hold for any non–equilibrium process. Taking advantage of this, we extend our
thermodynamic theory here to the non–equilibrium regime.

The rest of this paper is organized as follows. In Sec. II, our previous
results are summarized with regard to thermodynamics as applied to work in a
system–independent manner. We then derive the principle of non–equilibrium DL
minimum work to obtain the DL Massieu–Planck potentials and Helmholtz–Gibbs
potentials in the non–equilibrium regime. Results are verified in Sec. III by
numerical simulations using the SB model. Finally, Sec. IV presents concluding
remarks.

##  II Reflections on work of heat

In our previous paper,Koyanagi and Tanimura (2024a) we presented
system–specific thermodynamic laws based on open quantum dynamics theory.
Here, we develop the same laws by introducing several thermodynamic
“statements” without going into the details of the system, as in traditional
thermodynamic theory. By doing so, we clarify the distinctive features of
thermodynamic theory that allow it to treat systems in non–equilibrium
regimes.

###  II.1 Laws of thermodynamics as applied to work

We consider a thermodynamic system consisting of a subsystem A and a heat bath
B at an inverse temperature $\beta(t)$. The presence of a heat bath
temperature can be regarded as a consequence of the zeroth law of
thermodynamics.

An important statement of thermodynamics is that thermodynamic systems are
described in terms of intensive and extensive variables. In particular,
extensive variables that are proportional to the size are essential. This
statement, controversial, is sometimes called the fourth law.Oono (2017)
However, this is the premise of thermodynamics and we should call the minus
first law of thermodynamics.

The energy of the subsystem corresponds to internal energy and is expressed as
$U_{\rm A}(t)$, which changes with time owing to changes in the bath
temperature. Internal energy is an extensive variable, whereas inverse
temperature is an intensive variable. The external perturbation we consider
here is expressed as $-x(t)X_{\rm A}(t)$, where $x(t)$ and $X_{\rm A}(t)$ are
intensive and extensive variables, respectively. In traditional
thermodynamics, $x(t)$ is derived from the Euler relation. We then introduce
the total energy, which is related to enthalpy by ${H}_{\rm A}(t)=U_{\rm
A}(t)-x(t)X_{\rm A}(t)$, which can also be regarded as the Legendre
transformation between internal energy and enthalpy.

Because $U_{\rm A}(t)$ is the conjugate variable of $\beta(t)$, we introduce
the DL (or entropic) internal energy denoted by $\tilde{U}_{\rm A}(t)$, where
$\tilde{B}(t)\equiv\beta(t)B(t)$ for any variable $B(t)$. Because $\beta(t)$
diverges for $T\rightarrow 0$, states with $T=0$ do not exist, which
corresponds to the third law of thermodynamics.

From the time derivative of the DL internal energy $\tilde{U}_{\rm A}(t)$, we
have

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

where

$\displaystyle\frac{d\tilde{W}^{ext}_{\rm A}(t)}{dt}=U_{\rm
A}(t)\frac{d\beta(t)}{dt}+\tilde{x}(t)\frac{dX_{\rm A}(t)}{dt}$ (2)

and

$\displaystyle\frac{d\tilde{Q}_{\rm
A}^{ext}(t)}{dt}\equiv\beta(t)\frac{dU_{\rm
A}(t)}{dt}-\tilde{x}(t)\frac{dX_{\rm A}(t)}{dt}$ (3)

are the DL (or entropic) extensive work defined as being done from the outside
and DL (or entropic) extensive heat, respectively. For the enthalpy, we have

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

where

$\displaystyle\frac{d\tilde{W}_{\rm A}^{int}(t)}{dt}\equiv U_{\rm
A}(t)\frac{d\beta(t)}{dt}-X_{\rm A}(t)\frac{d\tilde{x}(t)}{dt}$ (5)

is the intensive work. Equations (1) and (4) correspond to the first law of
thermodynamics.

From the definitions, these intensive and extensive variables satisfy the
time-dependent Legendre (TDL) transformation expressed as

$\displaystyle\frac{d\tilde{W}_{\rm
A}^{int}(t)}{dt}=\frac{d\tilde{W}^{ext}_{\rm
A}(t)}{dt}-\frac{d}{dt}\left[\tilde{x}(t)X_{\rm A}(t)\right].$ (6)

Equation (1) is also expressed in the form of the TDL transformation as

$\displaystyle\frac{d\tilde{W}_{\rm
A}^{ext}(t)}{dt}=-\frac{d\tilde{Q}^{ext}_{\rm
A}(t)}{dt}+\frac{d}{dt}\left[\beta(t)U_{\rm A}(t)\right].$ (7)

The intensive heat $\tilde{Q}_{\rm A}^{int}(t)$ can also be defined from
$\tilde{Q}_{\rm A}^{ext}(t)$ using the TDL [Eq. (46) in Appendix A]. While the
above equations, including the first law of thermodynamics, hold for any
non–equilibrium processes described by intensive variables $\beta(t)$ and
$\tilde{x}(t)$ and extensive variables $U_{\rm A}(t)$ and $X_{\rm A}(t)$ which
are the state variables, as they are kinematic observables at time $t$, work
and heat are not state variables.

As the second law of thermodynamics, we adopt the DL minimum work principle
for the subsystem from one equilibrium state to another, which extends the
Kelvin-Planck statement for isothermal processes to thermostatic processes.
For work defined as being done from the outside to the subsystem, we have

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

where $\Phi_{\rm A}^{\rm qst}$ is the DL Massieu potential Massieu (1869) and
equality holds under quasi-static changes of natural variables $\beta(t)$ and
$X_{\rm A}(t)$ as $\Phi_{\rm A}^{\rm qst}[\beta^{\rm qst},X_{\rm A}^{\rm
qst}]=\tilde{W}_{\rm A}^{ext}[\beta^{\rm qst},X_{\rm A}^{\rm qst}]$. The above
inequality states that for any process occurring between two states, the work
done from the outside is minimized if the process is quasistatic (or
reversible), because no energy is dissipated in the heat bath. Note that if we
define work as being done from the subsystem to outside, the inequality in
Eq.(8) is reversed and the relationship is called the maximum work principle
(to the outside). Using the TDL transformation (6), we also have the Planck
potential Planck (1922) in terms of natural variables expressed as $\Xi_{\rm
A}^{\rm qst}[\beta^{\rm qst},\tilde{x}^{\rm qst}]$ (see Appendix A).

For the DL extensive heat $\tilde{Q}_{\rm A}^{ext}(t)$, we have

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

which corresponds to the principle of maximum entropy, where $\Lambda_{\rm
A}^{\rm qst}[U_{\rm A}^{\rm qst},{X}_{\rm A}^{\rm qst}]$ is the entropic
potential. Because both natural variables are extensive, this potential is
fundamental and we call it Massieu entropy (M-entropy). Using the TDL
tranformation Eq. (46), we also have Planck entropy (P-entropy) defined by
$\Gamma_{\rm A}^{\rm qst}[U_{\rm A}^{\rm qst},\tilde{x}^{\rm qst}]$ from
$\tilde{Q}_{\rm A}^{int}(t)$ (see Appendix A). Similar to the entropy obtained
from the partition function in statistical physics, this entropy is a function
of the extensive variable $\tilde{x}^{\rm qst}$.

Here, as the conjugate variable for $\beta(t)$, we chose the internal energy
$U_{\rm A}^{\rm qst}$, which is considered a fundamental variable in
thermodynamics, but enthalpy $H_{\rm A}^{\rm qst}$ could also be
chosen.Koyanagi and Tanimura (2024a) In such case, there are also two
entropies, one involving both extensive variables (Clasius entropy or
C-entropy) and one for intensive variable with respect to external forces
(Boltzmann entropy or B-entropy). The B-entropy and M-entopy are related to
the others by the Legendre transformations between $U_{\rm A}^{\rm qst}$ and
$H_{\rm A}^{\rm qst}$ and the values do not change, although the former
inculdes one intensive variable whereas the later includes extensive variables
only. However, there is no explicit relationship between C-entropy and
M-entropy, because the Legendre transformation between two entropy in the
$U_{\rm A}^{\rm qst}$ representation and the $H_{\rm A}^{\rm qst}$
representation are different. (See Appendix D in Ref. Koyanagi and Tanimura,
2024a). In any case, the Massieu–Planck potentials are entropic
potentials.Guggenheim (1986)

The total differential form of the Helmholtz–Gibbs potentials can be derived
from the Massieu–Planck potentials using the definitions $F^{\rm qst}_{\rm
A}(t)=-\Phi_{\rm A}^{\rm qst}(t)/\beta^{\rm qst}(t)$ and $G^{\rm qst}_{\rm
A}(t)=-\Xi_{\rm A}^{\rm qst}(t)/\beta^{\rm qst}(t)$. The results are
summarized and presented in Table 6 in Appendix A.

###  II.2 Thermodynamic potentials in a non–equilibrium regime

####  II.2.1 The DL non–equilibrium minimum work principle

On the basis of Eq. (8), we can define thermodynamic potentials applicable to
the non–equilibrium regime. Consider a non–equilibrium state $A$ and
equilibrium state $n$. The DL non–equilibrium-to-equilibrium minimum work
principle for $A\rightarrow n$ is expressed as (see Appendix B)

$\displaystyle(\tilde{W}_{\rm A}^{ext})_{{A}\rightarrow
n}\geq-(\Delta\Phi_{A}^{\rm neq})_{A\rightarrow n},$ (10)

where $\Phi_{\rm A}^{\rm neq}(t)$ is the non–equilibrium Massieu potential.
This inequality indicates that the path from non–equilibrium to equilibrium
has a lower bound on work. In engineering, the effective energy from non-
equilibrium to equilibrium is referred to as exergy.Rant (1956); Rosen and
Dincer (2001) The non–equilibrium thermodynamic potential we introduced here
can be regarded as a generalization of it.

The difference between $(\tilde{W}_{\rm A}^{ext})_{{A}\rightarrow n}$ and
$-(\Delta\Phi_{A}^{\rm neq})_{{A}\rightarrow n}$ corresponds to the DL heat
generated as loss work, denoted by $-(Q_{\rm A}^{\rm wst})_{{A}\rightarrow
n}$, which we call DL waste heat:

$\displaystyle-(\tilde{Q}_{\rm A}^{\rm wst})_{{A}\rightarrow
n}=(\tilde{W}_{\rm A}^{ext})_{{A}\rightarrow n}+(\Delta\Phi_{\rm A}^{\rm
neq})_{{A}\rightarrow n}\geq 0.$ (11)

Because the DL waste heat between two non–equilibrium states $A$ and $B$ is
evaluated as $(\tilde{Q}_{\rm A}^{\rm
wst})_{{A}\rightarrow{B}}=(\tilde{Q}_{\rm A}^{\rm wst})_{{A}\rightarrow
n}-(\tilde{Q}_{\rm A}^{\rm wst})_{B\rightarrow n}$, we the waste heat for any
non–equilibrium process at times $t$ and $t+dt$ (where $dt$ is an
infinitesimal time) in time-derivative form:

$\displaystyle-\frac{d\tilde{Q}_{\rm A}^{\rm
wst}(t)}{dt}=\frac{d\tilde{W}_{\rm A}^{ext}(t)}{dt}+\frac{d\Phi_{\rm A}^{\rm
neq}(t)}{dt}.$ (12)

As shown in Appendix B, the waste heat is nonpositive, i.e., $d\tilde{Q}_{\rm
A}^{\rm wst}(t)/dt\leq 0$. Thus, we have the inequality $(\tilde{W}_{\rm
A}^{ext})_{A\rightarrow B}^{\rm min}\geq-(\Delta\Phi_{\rm A}^{\rm neq})_{\rm
A\rightarrow B}$, or

$\displaystyle\frac{d\tilde{W}_{\rm A}^{ext}(t)}{dt}\geq-\frac{d\Phi_{\rm
A}^{\rm neq}(t)}{dt},$ (13)

which we call the DL non–equilibrium minimum work principle. While
$\tilde{W}_{\rm A}^{ext}(t)$ and $\tilde{Q}_{\rm A}^{\rm wst}(t)$ are not
state variables because they depend on a path, $\Phi_{\rm A}^{\rm neq}(t)$ is
a state variable defined by the non–equilibrium-to-equilibrium minimum work
path. When state $B$ is on the non–equilibrium-to-equilibrium minimum work
path from $A$ to $n$, equality holds in (13).

Using Eqs. (2) and (12), we obtain the time-derivative form of the
non–equilibrium Massieu potential as

$\displaystyle\frac{d\Phi_{\rm A}^{\rm neq}(t)}{dt}=-U_{\rm A}^{\rm
neq}\frac{d\beta(t)}{dt}-\tilde{x}(t)\frac{d{X}_{\rm A}^{\rm
neq}(t)}{dt}-\frac{d\tilde{Q}_{\rm A}^{\rm wst}(t)}{dt}.$ (14)

The convexity of DL potentials is discussed on the basis of the SB model in
Appendix C. Thus, we can regard these as thermodynamic potentials.

Using Eqs. (6), we have the non–equilibrium Planck potential expressed as

$\displaystyle\frac{d\Xi_{\rm A}^{\rm neq}(t)}{dt}=-U_{\rm A}^{\rm
neq}\frac{d\beta(t)}{dt}+{X}_{\rm A}^{\rm
neq}(t)\frac{d\tilde{x}(t)}{dt}-\frac{d\tilde{Q}_{\rm A}^{\rm wst}(t)}{dt}.$
(15)

The M-entropy $\Lambda_{\rm A}^{\rm neq}(t)$ and P-entropy $\Gamma_{\rm
A}^{\rm neq}(t)$ can be evaluated from the TDL transformations (6) and (7).
The non–equilibrium potentials satisfy the Legendre transformations:

$\displaystyle\Xi_{\rm A}^{\rm neq}(t)=\Phi_{\rm A}^{\rm
neq}(t)+\tilde{x}(t){X}_{\rm A}^{\rm neq}(t)$ (16)

and

$\displaystyle\Lambda_{\rm A}^{\rm neq}(t)=\Phi_{\rm A}^{\rm
neq}(t)+\beta(t)U_{\rm A}^{\rm neq}(t).$ (17)

The time-derivative forms of the DL non–equilibrium entropic potentials are
summarized in Table 1. Note that while the entropic potentials, the intensive
variables and non–equilbrium extensive variables are state variables, the
time–deribative of them including waste heat are not state variables.

As shown in Appendix D, the non–equilibrium Massieu potential is always
smaller than the quasi-static Massieu potential. This indicates that the
non–equilibrium Massieu potential is minimum when the state is equilibrium. In
the quasi-static case, Table 1 is reduced to Table 5 (Appendix A). Thus, we
can regard these potentials as extensions of the thermodynamic potentials to a
non–equilibrium regime.

From Eqs. (12) and (17), the waste heat current is expressed using the
non–equilibrium M-entropy as

$\displaystyle\frac{d\tilde{Q}_{\rm A}^{\rm wst}(t)}{dt}=\frac{d\tilde{Q}_{\rm
A}^{ext}(t)}{dt}-\frac{d\Lambda_{\rm A}^{\rm neq}(t)}{dt}.$ (18)

Table 1:  Time-derivative forms of the non–equilibrium (neq) Massieu–Planck potentials as functions of the intensive variables $\beta(t)$ and $\tilde{x}(t)$ and the extensive variables $U_{\rm A}^{\rm neq}(t)$ and $X_{\rm A}^{\rm neq}(t)$. Of the DL entropies, the commonly used one, which we call the Massieu entropy (M-entropy), involves only extensive variables and is denoted by $\Lambda_{\rm A}^{\rm neq}[U_{\rm A}^{\rm neq}(t),{X}_{\rm A}^{\rm neq}(t)]$, whereas the less widely used one, which we call the Planck entropy (P-entropy), is denoted by $\Gamma_{\rm A}^{\rm qst}[U_{\rm A}^{\rm qst},\tilde{x}]$. Each potential is related to the others by the Legendre transformations shown in the final column.  neq DL Thermodynamic pot. | Differential Form | Natural var. | Legendre Transformation
---|---|---|---
Massieu | $\frac{d}{dt}\Phi_{\rm A}^{\rm neq}=-U_{\rm A}^{\rm neq}\frac{d}{dt}\beta-\tilde{x}\frac{d}{dt}{X}_{\rm A}^{\rm neq}-\frac{d}{dt}\tilde{Q}^{\rm wst}_{\rm A}$ | $\beta,{X}_{\rm A}^{\rm neq}$ | $\cdots$
Planck | $\frac{d}{dt}\Xi^{\rm neq}_{\rm A}=-U_{\rm A}^{\rm neq}\frac{d}{dt}\beta+{X}_{\rm A}^{\rm neq}\frac{d}{dt}\tilde{x}-\frac{d}{dt}\tilde{Q}^{\rm wst}_{\rm A}$ | $\beta,\tilde{x}$ | $\Xi_{\rm A}^{\rm neq}=\Phi_{\rm A}^{\rm neq}+\tilde{x}{X}_{\rm A}^{\rm neq}$
M-Entropy | $\frac{d}{dt}\Lambda_{\rm A}^{\rm neq}=\beta\frac{d}{dt}U_{\rm A}^{\rm neq}-\tilde{x}\frac{d}{dt}X_{\rm A}^{\rm neq}-\frac{d}{dt}\tilde{Q}_{\rm A}^{\rm wst}$ | $U_{\rm A}^{\rm neq},X_{\rm A}^{\rm neq}$ | $\Lambda_{\rm A}^{\rm neq}=\Phi_{\rm A}^{\rm neq}+\beta U_{\rm A}^{\rm neq}$
P-Entropy | $\frac{d}{dt}\Gamma_{\rm A}^{\rm neq}=\beta\frac{d}{dt}U_{\rm A}^{\rm neq}+X_{\rm A}^{\rm neq}\frac{d}{dt}\tilde{x}-\frac{d}{dt}\tilde{Q}_{\rm A}^{\rm wst}$ | $U_{\rm A}^{\rm neq},\tilde{x}$ | $\Gamma_{\rm A}^{\rm neq}=\Lambda_{\rm A}^{\rm neq}+\tilde{x}{X}_{\rm A}^{\rm neq}$

####  II.2.2 Nonequilibrium Helmholtz–Gibbs potentials

Table 2:  Time-derivative forms of the non–equilibrium (neq) Helmholtz–Gibbs potentials as functions of the intensive variables $T(t)$ and $x(t)$ and the extensive variables $S_{\rm A}^{\rm neq}(t)$ and $X_{\rm A}^{\rm neq}(t)$, which are interrelated through the Legendre transformations shown in the final column. Because heat is always lost in non–equilibrium processes, the waste heat $Q^{\rm wst}(t)$ appears in the equations. neq Thermodynamic pot. | Differential Form | Natural var. | Legendre Transformation
---|---|---|---
Helmholtz Energy | $\frac{d}{dt}F^{\rm neq}_{\rm A}=-S^{\rm neq}_{\rm A}\frac{d}{dt}T+x\frac{d}{dt}X_{\rm A}^{\rm neq}+\frac{d}{dt}Q^{\rm wst}$ | $T,x$ | -
Gibbs Energy | $\frac{d}{dt}G^{\rm neq}_{\rm A}=-S^{\rm neq}_{\rm A}\frac{d}{dt}T-X_{\rm A}^{\rm neq}\frac{d}{dt}x+\frac{d}{dt}Q^{\rm wst}$ | $T,X_{\rm A}^{\rm neq}$ | $G^{\rm neq}_{\rm A}=F^{\rm neq}_{\rm A}-xX_{\rm A}^{\rm neq}$
Internal Energy | $\frac{d}{dt}U^{\rm neq}_{\rm A}=T\frac{d}{dt}S^{\rm neq}_{\rm A}+x\frac{d}{dt}X_{\rm A}^{\rm neq}+\frac{d}{dt}Q^{\rm wst}$ | $S_{\rm A}^{\rm neq},x$ | $U^{\rm neq}_{\rm A}=F^{\rm neq}_{\rm A}+TS^{\rm neq}_{\rm A}$
Enthalpy | $\frac{d}{dt}H^{\rm neq}_{\rm A}=T\frac{d}{dt}S^{\rm neq}_{\rm A}-X_{\rm A}^{\rm neq}\frac{d}{dt}x+\frac{d}{dt}Q^{\rm wst}$ | $S_{\rm A}^{\rm neq},X_{\rm A}^{\rm neq}$ | $H^{\rm neq}_{\rm A}=U^{\rm neq}_{\rm A}-xX_{\rm A}^{\rm neq}$

We introduce non–equilibrium Helmholtz and Gibbs energies defined as $F^{\rm
neq}_{\rm A}(t)=-\Phi_{\rm A}^{\rm neq}(t)/\beta(t)$ and $G^{\rm neq}_{\rm
A}(t)=-\Xi_{\rm A}^{\rm neq}(t)/\beta(t)$. Because $d\beta(t)/dt=-(1/k_{\rm
B}T^{2}(t))dT(t)/dt$, we obtain the time-derivative forms of these from Eqs.
(14) and (15) as

$\frac{dF_{\rm A}^{\rm neq}(t)}{dt}=-S_{\rm A}^{\rm
neq}(t)\frac{dT(t)}{dt}+x(t)\frac{dX_{\rm A}^{\rm neq}(t)}{dt}+\frac{dQ_{\rm
A}^{\rm wst}(t)}{dt}$ (19)

and

$\frac{dG_{\rm A}^{\rm neq}(t)}{dt}=-S_{\rm A}^{\rm
neq}(t)\frac{dT(t)}{dt}+X_{\rm A}^{\rm neq}(t)\frac{dx(t)}{dt}+\frac{dQ_{\rm
A}^{\rm wst}(t)}{dt},$ (20)

where $S_{\rm A}^{\rm neq}(t)=k_{\rm B}\Lambda_{\rm A}^{\rm neq}(t)$ is a
non–equilibrium entropy and we have introduced the waste heat current as

$\displaystyle\frac{dQ_{\rm A}^{\rm
wst}(t)}{dt}=\frac{1}{\beta(t)}\frac{d\tilde{Q}_{\rm A}^{\rm wst}(t)}{dt}.$
(21)

From Eq. (16), we obtain the TDL transformation between the non–equilibrium
Gibbs and Helmholtz potentials as

$F_{\rm A}^{\rm neq}(t)=G_{\rm A}^{\rm neq}(t)+x(t)X_{\rm A}^{\rm neq}(t).$
(22)

Solving Eq. (17) for $U_{\rm A}^{\rm neq}(t)$, we obtain the following TDL
transformations:

$U_{\rm A}^{\rm neq}(t)=F_{\rm A}^{\rm neq}(t)+T(t)S_{\rm A}^{\rm neq}(t).$
(23)

From the above, we obtain the time-derivative expressions for the enthalpy and
internal energy. The results are summarized in Table 2. These results reduce
to those in Table 6 (Appendix A) in the quasi-static case. Thus, we can regard
these potentials as extensions of the thermodynamic potentials to a
non–equilibrium regime.

Dividing both sides of Eq. (3) by $\beta(t)$ yields the first law of
thermodynamics expressed for the internal energy as

$\displaystyle\frac{dU_{\rm A}^{\rm neq}(t)}{dt}=\frac{dW_{\rm
A}^{ext}(t)}{dt}+\frac{dQ^{ext}_{\rm A}(t)}{dt},$ (24)

where

$\displaystyle\frac{dQ_{\rm
A}^{ext}(t)}{dt}=\frac{1}{\beta(t)}\frac{d\tilde{Q}_{\rm A}^{ext}(t)}{dt}$
(25)

is the extensive heat current. From Eqs. (18) and (25), we can evaluate the DL
waste heat current as

$\displaystyle-\frac{d\tilde{Q}_{\rm A}^{\rm wst}(t)}{dt}=\frac{1}{k_{\rm
B}}\left\\{\frac{dS_{\rm A}^{\rm neq}(t)}{dt}-\frac{1}{T(t)}\frac{dQ_{\rm
A}^{ext}(t)}{dt}\right\\}.$ (26)

In the isothermal case, by integrating both sides between the two equilibrium
states over time $t$, we obtain

$\displaystyle-\tilde{Q}_{\rm A}^{\rm wst}=\frac{1}{k_{\rm B}}\left\\{\Delta
S_{\rm A}^{\rm qst}-\frac{Q_{\rm A}^{ext}}{T}\right\\},$ (27)

which indicates that the DL waste heat is equivalent to the entropy production
evaluated for a heat engine.Sakamoto and Tanimura (2020); Koyanagi and
Tanimura (2022a)

The non–equilibrium Gibbs energy satisfies

$S^{\rm neq}_{\rm A}(t)=-\left(\frac{\partial G^{\rm neq}_{\rm A}}{\partial
T(t)}\right)_{x(t),\tilde{Q}^{\rm wst}},$ (28) $X_{\rm
A}(t)=-\left(\frac{\partial G_{\rm A}^{\rm neq}}{\partial
x(t)}\right)_{T(t),\tilde{Q}^{\rm wst}},$ (29)

and

$\displaystyle H_{\rm A}(t)=-T^{2}(t)\frac{\partial}{\partial
T(t)}\left(\frac{G^{\rm neq}_{\rm A}(t)}{T(t)}\right)_{x(t),\tilde{Q}^{\rm
wst}}$ (30)

Equation (30) is an extension of the Gibbs–Helmholtz relation to a
non–equilibrium regime.

##  III Numerical Demonstration

Although the results presented in Tables 1 and 2 hold for any non–equilibrium
system consisting of subsystem and bath, the extensive variables and waste
heat that appear in these relationships are system-specific, and there is no
general theory on the basis of which non–equilibrium thermodynamic potentials
can be obtained. However, it is possible to evaluate them as functions of the
intensive and extensive variables using an optimization algorithm. Here, as a
demonstration, we evaluate non–equilibrium thermodynamic potentials
numerically using the thermodynamic SB model.Koyanagi and Tanimura (2024a)

###  III.1 Thermodynamic system–bath model

We employ the Ullersma–Caldeira–Leggett (or Brownian) model,Ullersma (1966a,
b); Caldeira and Leggett (1983); Grabert, Schramm, and Ingold (1988); Weiss
(2012); Tanimura and Wolynes (1991, 1992) which is ideal for thermodynamic
simulations because the subsystem and bath are well defined, and rigorous
numerical solutions can be obtained in both classical and quantum cases under
any time-dependent external perturbation. Many of the favorable features for
thermodynamic investigations arise from the presence of a counter term, which
allows us to include the contribution of the SB interaction in the
bath.Tanimura (2015, 2006, 2020); Tanimura and Wolynes (1991, 1992); Kato and
Tanimura (2013) By introducing multiple heat baths at different temperatures
controlled by time-dependent SB coupling functions, we can investigate
isothermal, isentropic, thermostatic, and entropic processes. The total
Hamiltonian for isothermal and thermostatic processes is written asKoyanagi
and Tanimura (2024a)

$\hat{H}_{\rm tot}(t)=\hat{H}_{A}^{0}+\hat{H}_{\rm
A}^{\prime}(t)+\sum_{k=0}^{N}\hat{H}_{\rm IB}^{k}(t),$ (31)

where

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

is the unperturbed Hamiltonian of the subsystem with mass $m$ and potential
$U(\hat{q})$ described by momentum $\hat{p}$ and position $\hat{q}$. The
internal energy is then evaluated as $U_{\rm A}^{\rm
neq}(t)=\mathrm{tr}\\{\hat{H}_{\rm A}^{0}\hat{\rho}_{\rm A}(t)\\}$, where
$\hat{\rho}_{\rm A}(t)$ is the reduced density operator of the subsystem. The
external perturbation is expressed as $\hat{H}_{\rm
A}^{\prime}(t)\equiv-x(t)\hat{X}_{\rm A}$, where $\hat{X}_{\rm A}$ is an
operator of the subsystem coordinate [i.e., $\hat{X}_{\rm A}(\hat{q})$] and
$x(t)$ is regarded as the thermodynamic intensive variable. The extensive
variable is evaluated as $X_{\rm A}^{\rm neq}(t)=\mathrm{tr}\\{\hat{X}_{\rm
A}\hat{\rho}_{\rm A}(t)\\}$.

Although the conventional SB model has been limited to the investigation of
isothermal processes at constant temperature, we can extend it to describe
thermostatic processes in which temperature varies with time by introducing
$N$ independent heat baths, each in the thermal equilibrium state at the
inverse temperature $\beta_{k}\equiv 1/k_{\rm B}T_{k}$ connected to or
disconnected from subsystem A using the window function $\xi_{k}(t)$.Koyanagi
and Tanimura (2024a) The $k$th bath Hamiltonian is expressed as an ensemble of
harmonic oscillators and is given by

$\hat{H}_{\rm
IB}^{k}(t)\equiv\sum_{j}\left\\{\frac{(\hat{p}^{k}_{j})^{2}}{2m_{j}}+\frac{m_{j}^{k}(\omega_{j}^{k})^{2}}{2}\left[\hat{x}_{j}^{k}-\frac{c_{j}^{k}A_{k}\xi_{k}(t)\hat{q}}{m_{j}^{k}(\omega_{j}^{k})^{2}}\right]^{2}\right\\},$
(33)

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 $N$
independent heat baths, each in the thermal equilibrium state
$\exp(-\beta_{k}\hat{H}_{\rm I+B}^{k})$ at the inverse temperature
$\beta_{k}\equiv 1/k_{\rm B}T_{k}$, are connected or disconnected to subsystem
A according to a control function $\xi_{k}(t)$. The bath temperature is
effectively expressed as

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

or the inverse temperature as $\beta(t)=[k_{\rm B}T(t)]^{-1}$. As a spectral
distribution function for the $k$th bath
$J^{k}(\omega)\equiv\sum_{j}{\hbar(c_{j}^{k}A_{k})^{2}}/({2m_{j}^{k})\omega_{j}^{k}}\delta(\omega-\omega_{j}^{k})$,
we consider the Ohmic case described by

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

and assume that the time scale of the quantum thermal fluctuations
$\beta(t)\hbar/2\pi$ is shorter than the time scale of the external
perturbations $\beta(t)$ and $x(t)$. This allows us to extend the low-
temperature quantum Fokker–Planck equations (LT-QFPE)Ikeda and Tanimura (2019)
for the set of Wigner distribution functions $W_{\vec{n}}(p,q;t)$, where
${\vec{n}}$ is the index of hierarchy member, and the Kramers equation for
phase space distribution function $W(p,q;t)$ to apply the thermostatic case by
introducing time-dependent Matsubara frequencies $\nu(t)=1/\hbar\beta(t)$.
Expressions for these equations are presented in Ref. Koyanagi and Tanimura,
2024a. The source code for thermodynamic HEOM in phase space representation
and thermodynamic Kramers equation are provided as a separated paper.Koyanagi
and Tanimura (2024b).

It should be noted that in the non-Ohmic case, the temperature may change 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.

###  III.2 DL work and extensive variables

We express the solution for the reduced density elements under any $x(t)$ in
the Wigner representation using the zeroth member of the hierarchical Wigner
functions as $W(p,q,t)\equiv W_{\vec{0}}(p,q;t)$. In the classical limit
$\hbar\rightarrow 0$, ${W}(q,p;t)$ corresponds to the classical distribution
function. We use this to define the change in DL intensive work over time
corresponding to power and heat flow as follows:

$\frac{d\tilde{W}_{\rm A}^{int}(t)}{dt}=U_{\rm A}^{\rm
neq}(t)\frac{d\beta(t)}{dt}-X_{\rm A}^{\rm neq}(t)\frac{d\tilde{x}(t)}{dt},$
(36)

where the extensive variables in the Wigner representation at time $t$ are
expressed as

$\displaystyle X_{\rm A}^{\rm neq}(t)=\mathrm{tr}_{\rm A}\\{X_{\rm
A}(q){W}(p,q;t)\\}$ (37)

and

$\displaystyle U_{\rm A}^{\rm neq}(t)=\mathrm{tr}_{\rm
A}\left\\{\left[\frac{p^{2}}{2m}+U(q)\right]W(p,q;t)\right\\}.$ (38)

From the definition (2) of extensive work and the TDL transformations (6) and
(7), the extensive heat current is evaluated as

$\frac{dQ_{\rm A}^{ext}(t)}{dt}=\\\ \mathrm{tr}_{\rm
A}\left\\{\left[\frac{p^{2}}{2m}+U(q)-x(t)X_{\rm A}(q)\right]\frac{\partial
W(p,q;t)}{\partial t}\right\\}.$ (39)

###  III.3 Evaluation of the non–equilibrium Planck potential

Figure 1:  Schematic of the evaluation of the non–equilibrium Planck potential
$\Xi_{\rm A}^{\rm neq}(t)$ at $t_{k}$ and $t_{k+1}$. A non–equilibrium process
driven by the intensive variables $\beta(t)$ and $\tilde{x}(t)$ is described
as the non–equilibrium distribution $W(p,q;t)$ expressed by the black arrows.
The equilibrium states for different sets of $\beta^{\rm qst}_{n_{k}}$ and
$\tilde{x}_{n_{k}}^{\rm qst}$ are expressed by the red dashed line. For each
state $k$ to the equilibrium state $n_{k}$, we minimize $(\tilde{W}_{\rm
A}^{int})_{k\rightarrow n_{k}}+(\Xi_{\rm A}^{\rm qst})_{n_{k}}$ by choosing
the trial functions $\beta_{k}^{\rm trial}(t^{\prime})$ and
$\tilde{x}_{k}^{\rm trial}(t^{\prime})$, where $(\tilde{W}_{\rm
A}^{int})_{k\rightarrow n_{k}}^{\rm min}$ is the DL minimum intensive work for
$k\rightarrow n_{k}$ and $(\Xi_{\rm A}^{\rm qst})_{n_{k}}$ is the quasi-static
Planck potential at the state $n_{k}$. We can then evaluate the
non–equilibrium Planck potential at the state $k$ as $(\Xi_{\rm A}^{\rm
neq})_{k}=(\tilde{W}_{\rm A}^{int})_{k\rightarrow n_{k}}^{\rm min}+(\Xi_{\rm
A}^{\rm neq})_{n_{k}}$.

From Eq. (51), we obtain the expression for the non–equilibrium Planck
potential $(\Xi_{\rm A}^{\rm neq})_{k}=\Xi_{\rm A}^{\rm
neq}[\beta(t_{k}),\tilde{x}(t_{k})]$ of the non–equilibrium-to-equilibrium
process $k\rightarrow n_{k}$ as

$\displaystyle(\Xi_{\rm A}^{\rm neq})_{k}=(\tilde{W}_{\rm
A}^{int})_{k\rightarrow n_{k}}^{\rm min}+(\Xi_{\rm A}^{\rm qst})_{n_{k}},$
(40)

where $n_{k}$ represents an equilibrium state, $(\tilde{W}_{\rm
A}^{int})_{k\rightarrow n_{k}}^{\rm min}$ is the DL minimum intensive work for
$k\rightarrow n_{k}$, and $(\Xi_{\rm A}^{\rm qst})_{n_{k}}$ is the quasi-
static Planck potential at state $n_{k}$. As shown in Eq. (55), because the
right-hand side of Eq. (40) is independent of the choice of $n_{k}$, to
evaluate $(\Xi_{\rm A}^{\rm neq})_{k}$, we do not specify $n_{k}$ to minimize
$(\tilde{W}_{\rm A}^{int})_{k\rightarrow n_{k}}+(\Xi_{\rm A}^{\rm
qst})_{n_{k}}$.

To perform numerical calculations, we express Eq. (40) with Eq. (15) in terms
of trial functions $\beta^{\rm trial}(t_{k},t^{\prime})$ and $\tilde{x}^{\rm
trial}(t_{k},t^{\prime})$ as

$\displaystyle Xi_{\rm A}^{\rm target}[\beta^{\rm trial}(t),x^{\rm
trial}(t)]=\Xi^{\rm qst}_{\rm A}(\beta^{\rm qst}_{n_{k}},x^{\rm qst}_{n_{k}})$
$\displaystyle+\int_{t}^{t+\Delta t}\left[U_{\rm
A}(t^{\prime})\frac{d\beta^{\rm trial}(t^{\prime})}{dt^{\prime}}-X_{\rm
A}(t^{\prime})\frac{d\tilde{x}^{\rm
trial}(t^{\prime})}{dt^{\prime}}\right]dt^{\prime},$ (41)

where $\Xi_{\rm A}^{\rm target}[\beta^{\rm trial}(t),\tilde{x}^{\rm
trial}(t)]$ is the target function to be minimized and $\Xi_{\rm A}^{\rm
qst}(\beta_{n_{k}}^{\rm qst},\tilde{x}_{n_{k}}^{\rm qst})$ is the quasi-static
Planck potential at $n_{k}$, which is also evaluated from the LT-QFPE and
Kramers equation. The DL trail function $\tilde{x}_{k}^{\rm trial}(t)$ in
Eq.(III.3) is evaluated from $x_{k}^{\rm trial}(t)$ as $\tilde{x}_{k}^{\rm
trial}(t)=\beta_{k}^{\rm trial}(t)x_{k}^{\rm trial}(t)$. To reduce the
computational cost, we assume that the optimized trial functions become
constant after the characteristic time for the equilibration $\Delta t$, as
$\beta^{\rm qst}_{n_{k}}=\beta_{k}^{\rm trial}(t_{k}+\Delta t)$ and
$x_{n_{k}}^{\rm qst}=x_{k}^{\rm trial}(t_{k}+\Delta t)$.

As the trial functions for time $t^{\prime}>t_{k}$, we chose the
${N_{\beta}}$th and ${N_{x}}$th Taylor expansion forms expressed as

$\displaystyle\beta_{k}^{\rm
trial}(t^{\prime})=\sum_{n=0}^{N_{\beta}}\beta_{k}^{(n)}(t^{\prime}-t_{k})^{n},$
(42)

and

$\displaystyle x_{k}^{\rm
trial}(t^{\prime})=\sum_{n=0}^{N_{x}}x_{k}^{(n)}(t^{\prime}-t_{k})^{n},$ (43)

where $\beta_{k}^{(n)}$ and $x_{k}^{(n)}$ are the $n$th-order Taylor
coefficients. Thus, the functional minimization of $\Xi_{\rm A}^{\rm
target}[\beta^{\rm trial}(t),\tilde{x}^{\rm trial}(t)]$ becomes a
multivariable function minimization for $\beta_{k}^{(n)}$ and $x_{k}^{(n)}$.

Then, $\Xi_{\rm A}^{\rm neq}(t)$ is evaluated as follows (see Fig. 1). First,
we perform the non–equilibrium simulations for given $\beta(t)$ and
$\tilde{x}(t)$ to obtain $W_{\vec{n}}(p,q;t)$ for each $t_{k}$. The minimum
work from the non–equilibrium state $k$ to the equilibrium state $n_{k}$
expressed as $(\tilde{W}_{\rm A}^{int})_{k\rightarrow n_{k}}^{\rm min}$ is
then evaluated using an optimization algorithm for $\beta_{k}^{\rm
trial}(t^{\prime})$ and $\tilde{x}_{k}^{\rm trial}(t^{\prime})$. From Eq.
(40), we set this value as $(\Xi_{\rm A}^{\rm neq})_{k}$. By repeating this
procedure for different values of $k$, we obtain $\Xi_{\rm A}^{\rm
neq}(t_{k})$ at each step. The other potentials can be evaluated from the TDL
transformations.

###  III.4 Quantum and classical Stirling engine

The non–equilibrium thermodynamic potentials are state variables as functions
of the non–equilibrium intensive and extensive variables, which are also state
variables, while $Q_{\rm A}^{\rm wst}$ is not, and, in the quasi-static limit,
they reduce to conventional thermodynamic potentials. These non–equilibrium
potentials are useful because they can be used to analyze thermodynamic
processes using work diagrams, as in the case of equilibrium thermodynamics.
By finding a thermodynamic process to minimize waste heat, we can construct a
heat machine with maximum efficiency under non–equilibrium conditions,
although the efficiency is lower than the Carnot limit.Koyanagi and Tanimura
(2022a, b)

Here, we demonstrate how to evaluate non–equilibrium potentials for a case in
which $\beta(t)$ and $x(t)$ are specified. For this purpose, we consider the
classical/quantum Stirling engine (CSE/QSE) consisting of four steps: (i) a
hot isothermal process, (ii) a transition from a hot to cold thermostatic
process, (iii) a cold isothermal process, and (iv) a transition from a cold to
hot thermostatic process for an anharmonic potential system. Since the purpose
of the calculation is to illustrate how the potentials are calculated, we fix
the changes in $\beta(t)$ and $x(t)$ as shown in Fig. 2 with $\tau=20$ and do
not optimize them to improve their efficiency.

####  III.4.1 Simulation details

Figure 2:  Time profiles of (a) the electric field $x(t)$ and (b) the inverse
temperature $\beta(t)$ for the Stirling engine. The red, green, blue, and
orange lines represent (i) hot isothermal, (ii) isoelectric (from hot to
cold), (iii) cold isothermal, and (iv) isoelectric (from cold to hot)
processes, respectively. Time $t$ is normalized by the period of the cycle
$\tau$.

We performed simulations using the anharmonic potential system employed in our
previous study.Koyanagi and Tanimura (2024a) Thus, we considered a quartic
anharmonic potential with the external interaction described as $X_{\rm
A}(\hat{q})=\hat{q}$. The potential function in Eq. (32) is expressed as

$\displaystyle U(\hat{q})=U_{2}\hat{q}^{2}+U_{3}\hat{q}^{3}+U_{4}\hat{q}^{4},$
(44)

where the constants are given by $U_{2}=0.1$, $U_{3}=0.02$, and $U_{4}=0.05$.
Numerical calculations were carried out to integrate the thermodynamic HEOM in
the quantum cases and the thermodynamic Kramers equation in the classical
cases. The detailed conditions for the numerical calculations, including the
working parameters, are presented in Table 3 and in Ref. Koyanagi and
Tanimura, 2024a

Table 3:  Parameter values used for the simulations of the Stirling engine. 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-FPE. | $A$ | $N$ | $K$ | $dx$ | $dp$ |
---|---|---|---|---|---|---
Classical | $0.5$–$1.5$ | $\cdots$ | $\cdots$ | $0.25$ | $0.25$ |
Quantum | $0.5$ | $6$ | $2$ | $0.3$ | $0.4$ |
$1.0$ | $7$ | $2$ | $0.3$ | $0.5$ |
$1.5$ | $8$ | $2$ | $0.3$ | $0.6$ |

To set the trial functions defined in Eqs. (42) and (43), we chose
$N_{\beta}=N_{x}=5$. Then, using the Nelder–Mead method, we minimized
$\Xi_{\rm A}^{\rm target}[\beta^{\rm trial}(t),x^{\rm trial}(t)]$ in Eq.
(III.3) with a cutoff time $\Delta t=1.0$. We evaluated $\Xi_{\rm A}^{\rm
neq}(t)$ every 20 steps. The value of the quasi-static Planck potential
$\Xi_{\rm A}^{\rm qst}(\beta^{\rm qst},\tilde{x}^{\rm qst})$ in Eq. (III.3)
was obtained by integrating the thermodynamic LT-QFPE and the Kramers
equation.

####  III.4.2 Results

Figure 3:  $x(t)$–$X_{\rm A}^{\rm neq}(t)$ diagrams for the Stirling engine in the classical case (1, left column) and quantum case (2, right column) for (a) $A=0.5$ (weak), (b) $1.0$ (intermediate), and (c) $1.5$ (strong) SB coupling strengths. In each plot, the four curves (or lines) represent (i) hot isothermal (red), (ii) from hot to cold thermostatic (green), (iii) cold isothermal (blue), and (iv) from cold to hot thermostatic (orange) processes, respectively. The processes evolve in a counterclockwise fashion over time in a heat engine, whereas they evolve in a clockwise fashion over time in a refrigerator. Figure 4:  $T(t)$–$S_{\rm A}^{\rm neq}(t)$ diagrams for the Stirling engine in the classical case (1, left column) and quantum case (2, right column) for (a) $A=0.5$ (weak), (b) $1.0$ (intermediate), and (c) $1.5$ (strong) SB coupling strength, respectively. Each cycle starts with the red arrow, and the four curves represent (i) hot isothermal (red), (ii) from hot to cold thermostatic (green), (iii) cold isothermal (blue), and (iv) from cold to hot thermostatic (orange) processes, respectively. Table 4:  Work performed in one cycle for the classical and quantum cases for different SB coupling strengths. $A$ |  $W$ (classical) |  $W$ (quantum)
---|---|---
$0.5$ (weak) | $-1.44\times 10^{-2}$ | $-1.07\times 10^{-2}$
$1.0$ (intermediate) | $-5.77\times 10^{-3}$ | $-3.09\times 10^{-3}$
$1.5$ (strong) | $1.01\times 10^{-2}$ | $1.10\times 10^{-2}$

The extensive variables $X_{\rm A}^{\rm neq}(t)$ and $U_{\rm A}^{\rm neq}(t)$
were obtained from Eqs. (37) and (38), respectively. The non–equilibrium
entropy expressed as $S_{\rm A}^{\rm neq}(t)=k_{\rm B}\Lambda_{\rm A}^{\rm
neq}(t)$ was then obtained from $X_{\rm A}^{\rm neq}(t)$, and $U_{\rm A}^{\rm
neq}(t)$ and $\Xi_{\rm A}^{\rm neq}(t)$ were evaluated according to the
procedure described in Sec. III.3 using the TDL transformations Eqs. (16) and
(17). The results are depicted as non–equilibrium $x(t)$–$X_{\rm A}^{\rm
neq}(t)$ and $T(t)$–$S_{\rm A}^{\rm neq}(t)$ diagrams, which are closed
because the cycle is steady and because $T(t)$, $S_{\rm A}^{\rm neq}(t)$,
$x(t)$, and $X_{\rm A}^{\rm neq}(t)$ are state variables.

We first present the $x(t)$–$X_{\rm A}^{\rm neq}(t)$ diagrams for weak
($A=0.5$), intermediate ($A=1.0$), and strong ($A=1.5$) system–bath coupling
strengths in Fig. 3. The area enclosed by each diagram corresponds to positive
work when evolving in the counterclockwise direction over time, whereas it
corresponds to negative work when evolving clockwise. The work performed in
one cycle is presented in Table 4. In each diagram in Fig. 3, the red and blue
horizontal lines appear because there is a time delay in the change of $X_{\rm
A}^{\rm neq}(t)$ with respect to $x(t)$ for the heat bath to take effect.

As shown in Fig. 3 and Table 4, the larger the SB coupling, the smaller the
quantum effect, because it is suppressed by relaxation.Ikeda and Tanimura
(2019) Thus, when the coupling is large, dissipation becomes dominant, and the
system becomes a refrigerator (or damper). This is because the time scales of
fluctuations and dissipation are different: in the non–equilibrium case,
dissipation dominates when the SB coupling is large, whereas in the quasi-
static case, fluctuations and dissipation are balanced regardless of the
coupling strength. This is a distinct difference from the Carnot case; the
efficiency reaches a maximum in the intermediate SB coupling region where $A$
is neither large nor small.KT20JCP2

Next, we present the $T(t)$–$S_{\rm A}^{\rm neq}(t)$ diagrams in Fig. 4. The
area enclosed by a clockwise curve represents the difference between extensive
heat and waste heat, $Q_{\rm A}^{ext}-Q_{\rm A}^{\rm wst}$, per cycle. In the
quasi-static case, the area agrees with the extensive heat per cycle because
$Q^{\rm wst}$ becomes zero. In the thermostatic processes (green and orange
curves), when the SB coupling is small, there is a time delay in the change of
$S_{\rm A}^{\rm neq}(t)$ with respect to $T(t)$. This occurs because it takes
time for the system to be excited when thermal fluctuations are small, whereas
the delay is almost negligible when the SB coupling is strong.

Figure 5:  (a) Waste heat current $dQ^{\rm wst}(t)/dt$, (b) enthalpy, and (c)
internal energy as functions of $t$. The blue, green, and red curves represent
the weak ($A=0.5$), intermediate ($A=1.0$), and strong ($A=1.5$) SB coupling
cases, respectively. The solid and dashed curves represent the quantum and
classical results, respectively. The plots of waste heat currents have been
time-averaged over 10 steps each using the raw data.

In Fig. 5, we show the (time-averaged) waste heat current, enthalpy, and
internal energy as functions of $t$. As we show in Sec. II, the waste heat
current is always negative, and it becomes large in the isothermal case for
larger SB coupling because of the strong dissipation, whereas it becomes small
in the thermostatic processes for the case of weak SB coupling because
fluctuations are suppressed.

Using the graph of the waste heat current, we can improve the cycle
efficiency. For instance, in the case of strong coupling, because a large
waste heat current occurs in the isothermal processes, by increasing the
duration of the isothermal processes, we can reduce the waste heat in the
cycle.

For weak SB coupling, both the enthalpy and internal energy in the classical
case agree with those in the quantum case, whereas for strong SB coupling,
they disagree, owing to the quantum coherence between the system and bath
(bathentanglement).Tanimura (2020)

##  IV Conclusions

The virtue of thermodynamics lies in its ability to describe, in a system-
independent manner, macroscopic thermal phenomena resulting from complex
microscopic interactions as changes in thermodynamic potentials described as
interrelated intensive and extended variables through Legendre
transformations.

We have developed the laws of thermodynamics as applied to work in a system-
independent manner, based on the principle of DL minimum work. Subsequently,
we have developed a non–equilibrium thermodynamic theory that describes
non–equilibrium thermodynamic potential in terms of time-dependent intensive
and extensive variables, along with waste heat.

Thus, the non–equilibrium Massieu–Planck and Helmholtz–Gibbs potentials are
obtained in time-derivative form in terms of non–equilibrium extensive and
intensive variables, which are state variables along with waste heat, which is
not a state variable. Our results are summarized in Tables 1 and 2 and are
consistent with traditional thermodynamics in a quasi-static case, as
presented in Tables 5 and 6.

These results have been validated for quantum and classical Stirling engines
through numerical simulations based on the thermodynamic SB model, which can
describe both isothermal and thermostatic processes. An optimization algorithm
has been used to evaluate non–equilibrium thermodynamic potentials. Work
diagrams in non–equilibrium regimes have been presented and analyzed.

Although the time evolution of the system in the non–equilibrium regime is
model-specific, it is possible to test our theory even in real systems through
the use of optimization algorithms. By designing processes to reduce waste
heat, it may be possible to develop efficient heat engines even for
non–equilibrium processes.

In engineering, the effective energy from non-equilibrium to equilibrium is
referred to as exergy.Rant (1956); Rosen and Dincer (2001) The non–equilibrium
thermodynamic potentials we have introduced can be regarded as a
generalization of it. The present theory provides a methodology to
systematically evaluate and improve it, which should be useful from the
perspective of Sustainable Development Goal (SDGs).

## Acknowledgments

The authors thank Yoshi Oono for his critical comments on the definition and
terminology of thermodynamic variables, particularly for extensive variables.
Y.T. was supported by JSPS KAKENHI (Grant No. B21H01884). S.K. was supported
by JST fellowship, the establishment of university fellowships towards the
creation of science technology innovation (Grant No. JPMJFS2123), and 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 Quasi-static thermodynamic potentials

Table 5:  Total differential expressions for the quasi-static (qst.) entropic potentials as functions of the intensive variables $\beta^{\rm qst}(t)$ and $\tilde{x}^{\rm qst}(t)$ and the extensive variables $U_{\rm A}^{\rm qst}(t)$ and ${X}_{\rm A}^{\rm qst}(t)$. Entropy has two definitions, depending on whether the work variable is intensive or extensive. Of these DL entropies, the commonly used one, which we call Massieu entropy (M-entropy), involves only extensive variables and is denoted by $\Lambda_{\rm A}^{\rm qst}[U_{\rm A}^{\rm qst},{X}_{\rm A}^{\rm qst}]$, whereas the less widely used one, which we call Planck entropy (P-entropy), is denoted by $\Gamma_{\rm A}^{\rm qst}[U_{\rm A}^{\rm qst},\tilde{x}^{\rm qst}]$. Whereas the enthalpy $H_{\rm A}^{\rm qst}(t)$ was chosen as the natural variable in Ref. Koyanagi and Tanimura, 2024a, here we chose the internal energy $U_{\rm A}^{\rm qst}$ instead. Each potential is related to the others by the Legendre transformations shown in the final column. Qst. Potential | Differential form | Natural var. | Legendre transformation
---|---|---|---
Massieu | $d\Phi_{\rm A}^{\rm qst}=-U_{\rm A}^{\rm qst}d\beta^{\rm qst}-\tilde{x}^{\rm qst}d{X}_{\rm A}^{\rm qst}$ | $\beta^{\rm qst},{X}_{\rm A}^{\rm qst}$ | $\cdots$
Planck | $d\Xi_{\rm A}^{\rm qst}=-U_{\rm A}^{\rm qst}d\beta^{\rm qst}+{X}_{\rm A}^{\rm qst}d\tilde{x}^{\rm qst}$ | $\beta^{\rm qst},\tilde{x}^{\rm qst}$ | $\Xi_{\rm A}^{\rm qst}=\Phi_{\rm A}^{\rm qst}+\tilde{x}^{\rm qst}{X}_{\rm A}^{\rm qst}$
M-Entropy | $d\Lambda_{\rm A}^{\rm qst}=\beta^{\rm qst}dU_{\rm A}^{\rm qst}-\tilde{x}^{\rm qst}dX_{\rm A}^{\rm qst}$ | $U_{\rm A}^{\rm qst},X_{\rm A}^{\rm qst}$ | $\Lambda_{\rm A}^{\rm qst}=\Xi_{\rm A}^{\rm qst}+\beta^{\rm qst}U^{\rm qst}_{\rm A}$
P-Entropy | $d\Gamma_{\rm A}^{\rm qst}=\beta^{\rm qst}dU_{\rm A}^{\rm qst}+X_{\rm A}^{\rm qst}d\tilde{x}^{\rm qst}$ | $U_{\rm A}^{\rm qst},\tilde{x}^{\rm qst}$ | $\Gamma_{\rm A}^{\rm qst}=\Phi_{\rm A}^{\rm qst}+\beta^{\rm qst}U^{\rm qst}_{\rm A}$

We introduce the DL intensive heat defined as

$\displaystyle\frac{d\tilde{Q}_{\rm
A}^{int}(t)}{dt}\equiv\beta(t)\frac{dU_{\rm A}(t)}{dt}+X_{\rm
A}(t)\frac{d\tilde{x}(t)}{dt}$ (45)

which satisfies the Legendre transformation

$\displaystyle\frac{d\tilde{Q}_{\rm
A}^{int}(t)}{dt}=\frac{d\tilde{Q}^{ext}_{\rm
A}(t)}{dt}+\frac{d}{dt}\left[\tilde{x}(t)X_{\rm A}(t)\right].$ (46)

The extensive and intensive work and heat $\tilde{W}_{\rm A}^{ext}(t)$,
$\tilde{W}_{\rm A}^{int}(t)$, $\tilde{Q}_{\rm A}^{ext}(t)$, and
$\tilde{Q}_{\rm A}^{int}(t)$ are interrelated through the Legendre
transformations (6), (7), and (46). Therefore, the inequalities (8), (9),

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

and

$\displaystyle\tilde{Q}_{\rm A}^{int}\leq\Delta\Gamma_{\rm A}^{\rm qst},$ (48)

where $\Gamma_{\rm A}^{\rm qst}$ is the Planck entropy (P-entropy), are all
expressions of the second law of thermodynamics.

The total differential forms of the Massieu–Planck potentials are presented in
Table 5.

The Helmholtz–Gibbs potentials can be obtained form the Massieu–Planck
potentials using the definitions $F^{\rm qst}_{\rm A}(t)=-\Phi_{\rm A}^{\rm
qst}(t)/\beta^{\rm qst}(t)$ and $G^{\rm qst}_{\rm A}(t)=-\Xi_{\rm A}^{\rm
qst}(t)/\beta^{\rm qst}(t)$. From Eq. (6), we have

$F_{\rm A}^{\rm qst}(t)=G_{\rm A}^{\rm qst}(t)+x^{\rm qst}(t)X_{\rm A}^{\rm
qst}(t).$ (49)

Accordingly, from Eq. (7), we have

$U_{\rm A}^{\rm qst}(t)=F_{\rm A}^{\rm qst}(t)+T^{\rm qst}(t)S_{\rm A}^{\rm
qst}(t),$ (50)

where we have used $d\beta^{\rm qst}(t)/dt=-(1/k_{\rm B}[T^{\rm
qst}(t)]^{2})dT^{\rm qst}(t)/dt$ and $S_{\rm A}^{\rm qst}(t)=k_{\rm
B}\Lambda_{\rm A}^{\rm qst}(t)$.

The total differential forms of the Helmholtz–Gibbs potentials are presented
in Table 6.

Table 6:  Total differential expressions for the quasi-static (qst.) thermodynamic potentials as functions of the intensive variables $T^{\rm qst}(t)$ and $x^{\rm qst}(t)$ and the extensive variables $S_{\rm A}^{\rm qst}(t)$ and ${X}_{\rm A}^{\rm qst}(t)$, which are related through the Legendre transformations shown in the final column.Koyanagi and Tanimura (2024a) Qst. potential | Differential form | Natural var. | Legendre transformation
---|---|---|---
Helmholtz | $dF_{\rm A}^{\rm qst}=-S_{\rm A}^{\rm qst}dT^{\rm qst}+x^{\rm qst}dX_{\rm A}^{\rm qst}$ | $T^{\rm qst},X_{\rm A}^{\rm qst}$ | $\cdots$
Gibbs | $dG_{\rm A}^{\rm qst}=-S_{\rm A}^{\rm qst}dT^{\rm qst}-X_{\rm A}^{\rm qst}dx^{\rm qst}$ | $T^{\rm qst},x^{\rm qst}$ | $G_{\rm A}^{\rm qst}=F_{\rm A}^{\rm qst}-x^{\rm qst}X_{\rm A}^{\rm qst}$
Internal | $dU_{\rm A}^{\rm qst}=T^{\rm qst}dS_{\rm A}^{\rm qst}+x^{\rm qst}dX_{\rm A}^{\rm qst}$ | $S_{\rm A}^{\rm qst},X_{\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}-X_{\rm A}^{\rm qst}dx^{\rm qst}$ | $S_{\rm A}^{\rm qst},x^{\rm qst}$ | $H^{\rm qst}_{\rm A}=G_{\rm A}^{\rm qst}+T^{\rm qst}S_{\rm A}^{\rm qst}$

##  Appendix B DL non–equilibrium minimum work principle

Figure 6:  Schematic to derive Eq. (51). The blue circles 1 and 2 represent
the equilibrium states, and the red square A represents the non–equilibrium
state. The quasi-static work performed for $1\rightarrow 2$ (red dashed line)
is denoted by $(\tilde{W}_{\rm A}^{ext})_{1\rightarrow 2}^{\rm min}$, while
the non–equilibrium work from 1 to A (dark blue curve) is denoted by
$(\tilde{W}_{\rm A}^{ext})_{1\rightarrow A}$. When $(\tilde{W}_{\rm
A}^{ext})_{1\rightarrow A}$ is fixed, the value of $(\tilde{W}_{\rm
A}^{ext})_{A\rightarrow 2}$ depends on the path (green and light blue curves).
From the second law of thermodynamics, we have $(W_{\rm
A}^{ext})_{1\rightarrow A}+(\tilde{W}_{\rm A}^{ext})_{\rm A\rightarrow
2}\geq-(\Delta\Phi_{\rm A}^{\rm qst})_{1\rightarrow 2}$. This indicates a
lower bound of work for ${A\rightarrow 2}$ denoted by $(\tilde{W}_{\rm
A}^{ext})_{A\rightarrow 2}^{\rm min}$ (light blue curve) and given by Eq.
(51).

Here, we derive fundamental equations to develop the DL non–equilibrium
minimum work principle. Consider two equilibrium states 1 and 2 that are
connected by a quasi-static process (see Fig. 6). From the principle of DL
minimum work [Eq. (8)], the work done in this process is equivalent to the
difference in Massieu potentials: $(\tilde{W}_{\rm A}^{ext})_{1\rightarrow
2}^{\rm min}=-(\Delta\Phi_{\rm A}^{\rm qst})_{1\rightarrow 2}\equiv(\Phi_{\rm
A}^{\rm qst})_{1}-(\Phi_{\rm A}^{\rm qst})_{2}$, where ${n\rightarrow
n^{\prime}}$ represents the transition from any state $n$ to $n^{\prime}$.
Separately, we consider the non–equilibrium state $A$ and introduce the
non–equilibrium process $1\rightarrow A$ (dark blue curve in Fig. 6), and any
process from $A$ to 2 (green and light blue curves). For each of these
processes, the work done is denoted by $(\tilde{W}_{\rm
A}^{ext})_{1\rightarrow A}$ and $(\tilde{W}_{\rm A}^{ext})_{A\rightarrow 2}$,
respectively. For fixed $(\tilde{W}_{\rm A}^{ext})_{1\rightarrow A}$, the
value of $(\tilde{W}_{\rm A}^{ext})_{A\rightarrow 2}$ changes depending on the
path.

From the second law of thermodynamics [Eq. (8)], we then have $(\tilde{W}_{\rm
A}^{ext})_{1\rightarrow A}+(\tilde{W}_{\rm A}^{ext})_{A\rightarrow
2}>-(\Delta\Phi_{\rm A}^{\rm qst})_{1\rightarrow 2}.$ This indicates that
$(\tilde{W}_{\rm A}^{ext})_{A\rightarrow 2}$ has a lower bound for fixed
$(\tilde{W}_{\rm A}^{ext})_{1\rightarrow A}$; otherwise the inequality is
violated, for example, as $-\infty>-(\Delta\Phi_{\rm A}^{\rm
qst})_{1\rightarrow 2}-(\tilde{W}_{\rm A}^{ext})_{1\rightarrow A}$. Using the
lower bound of the work, we introduce the non–equilibrium Massieu potential
$\Phi_{\rm A}^{\rm neq}$ defined as

$\displaystyle(\tilde{W}_{\rm A}^{ext})_{A\rightarrow 2}^{\rm
min}=-[(\Phi_{\rm A}^{\rm qst})_{2}-(\Phi_{\rm A}^{\rm neq})_{A}].$ (51)

Figure 7: Schematic to derive the principle of DL non–equilibrium minimum work
expressed by Eq. (53). We consider two non–equilibrium states $A$ and $B$ (red
square) that are connected by a non–equilibrium process described by the red
arrow for $A\rightarrow B$. The light blue curves $A\rightarrow 2$ and
$B\rightarrow 2$ represent non–equilibrium-to-equilibrium minimum work paths
whose work is denoted by $(\tilde{W}_{\rm A}^{ext})_{A\rightarrow 2}^{\rm
min}$ and $(\tilde{W}_{\rm A}^{ext})_{B\rightarrow 2}^{\rm min}$,
respectively.

Introducing a second non–equilibrium state $B$ on the pathway $A\rightarrow 2$
as depicted in Fig. 7, we now discuss the non–equilibrium transition
${A\rightarrow B}$ (red arrow). We consider the non–equilibrium-to-equilibrium
minimum work from $B$ to 2 (light blue curve), expressed as

$\displaystyle(\tilde{W}_{\rm A}^{ext})_{B\rightarrow 2}^{\rm
min}=-[(\Phi_{\rm A}^{\rm qst})_{2}-(\Phi_{\rm A}^{\rm neq})_{B}].$ (52)

From the inequality $(\tilde{W}_{\rm A}^{ext})_{A\rightarrow
2}\geq-[(\Phi_{\rm A}^{\rm qst})_{2}-(\Phi_{\rm A}^{\rm neq})_{A}]$ and Eq.
(52), we have

$\displaystyle(\tilde{W}_{\rm A}^{ext})_{A\rightarrow B}\geq-(\Delta\Phi_{\rm
A}^{\rm neq})_{A\rightarrow B},$ (53)

where $(\tilde{W}_{\rm A}^{ext})_{A\rightarrow B}\equiv(\tilde{W}_{\rm
A}^{ext})_{A\rightarrow 2}-(\tilde{W}_{\rm A}^{ext})_{B\rightarrow 2}^{\rm
min}$. Thus, equality holds in (53) when $B$ in on the minimal pathway
$(A\rightarrow 2)^{\rm min}$.

Because a non–equilibrium process changes as a function of time, it is more
convenient to use the time-derivative form of the above. Thus, for the process
from $A$ to $B$ at times $t$ and $t+dt$, where $dt$ is an infinitesimal time,
we have the following inequality:

$\displaystyle\frac{d\tilde{W}_{\rm A}^{ext}(t)}{dt}\geq-\frac{d\Phi_{\rm
A}^{\rm neq}(t)}{dt}.$ (54)

This inequality is an extension of the principle of minimum work to a
non–equilibrium regime.

Figure 8:  Schematic showing that the non–equilibrium Massieu potential can be
evaluated from any equilibrium state 3 instead of 2. To illustrate this, we
consider the non–equilibrium-to-equilibrium path $A\rightarrow 3$, where we
have introduced a third equilibrium state denoted by 3 (blue circle). The work
between the two equilibrium states is denoted by $(W_{\rm
A}^{ext})_{2\rightarrow 3}^{\rm qst}$ and $(W_{\rm A}^{ext})_{3\rightarrow
2}^{\rm qst}$, respectively.

Although we have derived Eqs. (53) and (54) for the specific equilibrium state
2, actually the value $(\Phi_{\rm A}^{\rm neq})_{A}$ is the same for any
quasi-equilibrium state along $\beta^{\rm qst}(t)$ and $x^{\rm qst}(t)$. To
illustrate this, we introduce a third equilibrium state 3 and consider the
transition $A\rightarrow 3$ (see Fig. 8). Because $(\tilde{W}_{\rm
A}^{ext})_{A\rightarrow 2}^{\rm min}+(\tilde{W}_{\rm A}^{ext})_{2\rightarrow
3}^{\rm qst}\geq(\tilde{W}_{\rm A}^{ext})_{A\rightarrow 3}^{\rm min}$ for
$A\rightarrow 2\rightarrow 3$ and $(\tilde{W}_{\rm A}^{ext})_{A\rightarrow
3}^{\rm min}+(\tilde{W}_{\rm A}^{ext})_{3\rightarrow 2}^{\rm
qst}\geq(\tilde{W}_{\rm A}^{ext})_{A\rightarrow 2}^{\rm min}$ for
$A\rightarrow 2\rightarrow 3$, we have

$\displaystyle(W_{\rm A}^{ext})_{A\rightarrow 2}^{\rm min}+(\Phi_{\rm A}^{\rm
neq})_{2}=(W_{\rm A}^{ext})_{A\rightarrow 3}^{\rm min}+(\Phi_{\rm A}^{\rm
neq})_{3},$ (55)

where we have used $(W_{\rm A}^{ext})_{2\rightarrow 3}^{\rm qst}=(\Phi_{\rm
A}^{\rm qst})_{3}-(\Phi_{\rm A}^{\rm qst})_{2}$. From Eq. (51), this indicates
that the non–equilibrium Massieu potential is independent of the choice of the
equilibrium state 2.

Because the other potentials are interrelated through TDL transformations, we
have

$\displaystyle\frac{d\tilde{W}_{\rm A}^{int}(t)}{dt}\geq-\frac{d\Xi_{\rm
A}^{\rm neq}(t)}{dt},$ (56) $\displaystyle\frac{d\tilde{Q}_{\rm
A}^{ext}(t)}{dt}\leq\frac{d\Lambda_{\rm A}^{\rm neq}(t)}{dt},$ (57)

and

$\displaystyle\frac{d\tilde{Q}_{\rm A}^{int}(t)}{dt}\leq\frac{d\Gamma_{\rm
A}^{\rm neq}(t)}{dt}.$ (58)

For an isothermal case, from the definitions $F^{\rm neq}_{\rm
A}(t)=-\Phi_{\rm A}^{\rm neq}(t)/\beta$ and $G^{\rm neq}_{\rm A}(t)=-\Xi_{\rm
A}^{\rm neq}(t)/\beta$, we have

$\displaystyle\frac{dW_{\rm A}^{ext}(t)}{dt}\geq\frac{dF_{\rm A}^{\rm
neq}(t)}{dt}$ (59)

and

$\displaystyle\frac{dW_{\rm A}^{int}(t)}{dt}\geq\frac{dG_{\rm A}^{\rm
neq}(t)}{dt}.$ (60)

However, unlike in the Massieu–Planck case, the equality in Eqs. (59) and (60)
cannot be obtained in general because the minimum DL work path in Fig. 6 is
evaulaed by optimizing both $\beta(t)$ and $x(t)$, while $\beta$ has been
fixed to introduce the Helmholtz energy as $F^{\rm neq}_{\rm A}(t)=-\Phi_{\rm
A}^{\rm neq}(t)/\beta$.

The above discussion indicates that by including the thermostatic process, we
may develop an efficient heat engine under non–equilibrium conditions.

##  Appendix C Convexity of non–equilibrium thermodynamic potentials as
functions of extensive variables

Table 7:  Inequalities holding between the non–equilibrium state at time $t_{1}$ and the equilibrium state at time $t_{2}$, together with the conditions under which these inequalities hold. From these inequalities, we find that, for example, for a given inverse temperature and electric field, the non–equilibrium Planck potential takes its maximum value when the state is in equilibrium. neq Potentials | Convexity and Concavity | Condition
---|---|---
Massieu | $(1-\lambda)(\Phi_{\rm A}^{\rm neq})_{A}+\lambda(\Phi_{\rm A}^{\rm neq})_{B}\leq(\Phi_{\rm A}^{\rm neq})_{\lambda}$ | $\beta_{A}=\beta_{B}=\beta_{\lambda}$
Planck | $(1-\lambda)(\Xi_{\rm A}^{\rm neq})_{A}+\lambda(\Xi_{\rm A}^{\rm neq})_{B}\leq(\Xi_{\rm A}^{\rm neq})_{\lambda}$ |  $\beta_{A}=\beta_{B}=\beta_{\lambda}$ and $\tilde{x}_{A}=\tilde{x}_{B}=\tilde{x}_{\lambda}$
Helmholtz | $(1-\lambda)(F_{\rm A}^{\rm neq})_{A}+\lambda(F_{\rm A}^{\rm neq})_{B}\geq(F_{\rm A}^{\rm neq})_{\lambda}$ | $\beta_{A}=\beta_{B}=\beta_{\lambda}$
Gibbs | $(1-\lambda)(G_{\rm A}^{\rm neq})_{A}+\lambda(G_{\rm A}^{\rm neq})_{B}\geq(G_{\rm A}^{\rm neq})_{\lambda}$ |  $\beta_{A}=\beta_{B}=\beta_{\lambda}$ and $x_{A}=x_{B}=x_{\lambda}$
Figure 9:  Schematic for deriving convexity relations. The red squares
represent the non–equilibrium states. The arrows are the relaxation paths to
equilibrium state 2 described by the time-evolution operator $\hat{U}_{\rm
A}(t,t_{0})$: The blue arrow is the DL minimum work path from state $\lambda$.

To obtain stable thermodynamic properties, the thermodynamic potential must be
a convex function of the work variables. Even in the non–equilibrium case, we
find that convexity holds under certain conditions. Here, we show this in the
case of the SB model. Because other potentials can be discussed similarly, we
limit our discussion here to the case of a non–equilibrium Massieu potential.

We consider three non–equilibrium states $A$, $B$, and $\lambda$
($\lambda\in[0,1]$) with the same initial inverse temperature $\beta(t_{0})$
at time $t_{0}$. We then assume that the state $\lambda$ satisfies the
relation

$\displaystyle(W_{\rm A})_{\lambda}=(1-\lambda)(W_{\rm A})_{A}+\lambda(W_{\rm
A})_{B},$ (61)

where $(W_{\rm A})_{A}=\left(W_{\rm A}(p,q;t)\right)_{A}$ and $(W_{\rm
A})_{B}=\left(W_{\rm A}(p,q;t)\right)_{B}$ are Wigner distributions in the
quantum case (or phase space distributions in the classical case) in the $A$
and $B$ states. The time-evolution operator for a subsystem is expressed as
$\hat{U}_{\rm
A}(t,t_{0})=\exp\\!\big{[}-\int\hat{\mathscr{L}}(p,q)\,dt\big{]}$, where
$\hat{\mathscr{L}}(p,q)$ is the Liouville operator, which includes fluctuation
and dissipation operators for equilibration. In the quantum case,
$\hat{\mathscr{L}}(p,q)$ acts on $W_{\rm A}$, which is represented by
hierarchical elements.Tanimura (2006, 2020); Ikeda and Tanimura (2019) We then
obtain the relations for the extensive variable and internal energies as

$\displaystyle({X}_{\mathrm{A}})_{\lambda}(t)=(1-\lambda)({X}_{\mathrm{A}})_{A}(t)+\lambda({X}_{\mathrm{A}})_{B}(t)$
(62)

and

$\displaystyle(U_{\mathrm{A}})_{\lambda}(t)=(1-\lambda)(U_{\mathrm{A}})_{A}(t)+\lambda(U_{\mathrm{A}})_{B}(t),$
(63)

where we have defined $(X_{\mathrm{A}})_{\alpha}(t)=\mathrm{tr}_{\rm
tot}\\{X_{\rm A}(q)\hat{U}_{\rm tot}(t,t_{0})(W_{\rm tot})_{\alpha}\\}$ and
$(U_{\mathrm{A}})_{\alpha}(t)=\mathrm{tr}_{\rm
tot}\\{[p^{2}/2m+U(q)]\hat{U}_{\rm tot}(t,t_{0})(W_{\rm tot})_{\alpha}\\}$ for
$\alpha=A,B,$ and $\lambda$. Thus, from Eq. (2), the DL extensive work
$(W_{\rm A}^{ext})_{\alpha}(t,t_{0})$ ($\alpha=A,B,$ and $\lambda$), performed
between times $t_{0}$ and $t$ also satisfies the equality

$\displaystyle(W_{\rm A}^{ext})_{\lambda}(t,t_{0})=(1-\lambda)(W_{\rm
A}^{ext})_{A}(t,t_{0})+\lambda(W_{\rm A}^{ext})_{B}(t,t_{0}).$ (64)

The DL minimum work path from the state $\lambda$ to the equilibrium state 2
is described by $\hat{U}_{\rm A}(t,t_{0})$ (see Fig. 9). This indicates that
when $t$ is large, the intensive variables must be
$\beta(t)\rightarrow\beta_{2}$ and ${x}(t)\rightarrow x_{2}$, where
$\beta_{2}$ and $x_{2}$ are the values at the equilibrium state 2. Therefore,
the processes starting from $A$ and $B$ also relax to the equilibrium state 2
as $t\rightarrow\infty$. Then, from Eq. (64), we obtain

$\displaystyle(W_{\rm A}^{ext})_{\lambda\rightarrow 2}^{\rm
min}=(1-\lambda)(W_{\rm A}^{ext})_{A\rightarrow 2}+\lambda(W_{\rm
A}^{ext})_{B\rightarrow 2}.$ (65)

Using the inequality $(W_{\rm A}^{ext})_{\alpha\rightarrow 2}\geq(W_{\rm
A}^{ext})_{\alpha\rightarrow 2}^{\rm min}$ for the states $\alpha=A$ and $B$,
and Eq. (51), we obtain the inequality to prove convexity as

$\displaystyle(1-\lambda)(\Phi_{\rm A}^{\rm neq})_{A}+\lambda(\Phi_{\rm
A}^{\rm neq})_{B}\leq(\Phi_{\rm A}^{\rm neq})_{\lambda},$ (66)

where $(\Phi_{\rm A}^{\rm neq})_{\alpha}$ is the non–equilibrium Massieu
potential in the state $\alpha$ ($\alpha=A$, $B$, and $\lambda$). When $A$ and
$B$ are in equilibrium states, we obtain the convexity relation in the quasi-
static case for fixed inverse temperature using the inequality for the Massieu
potential in Table 8 as Koyanagi and Tanimura (2024a)

$\displaystyle(1-\lambda)(\Phi_{\rm A}^{\rm qst})_{A}+\lambda(\Phi_{\rm
A}^{\rm qst})_{B}\leq(\Phi_{\rm A}^{\rm qst})_{\lambda}.$ (67)

We summarize the convexity and concavity properties for the other
non–equilibrium thermodynamic potential in Table 7.

##  Appendix D Inequality between non–equilibrium and quasi-static
thermodynamic potentials

Table 8:  Inequalities for Massieu–Planck potentials holding between the non–equilibrium state $A$ at time $t_{1}$ and the equilibrium state 2 at time $t_{2}$, together with the conditions under which these inequalities hold. From these inequalities, we find that, for example, for a given inverse temperature and intensive variable, the non–equilibrium Planck potential takes its maximum value when $A$ is in the equilibrium state. neq Potentials | Inequalities | Condition
---|---|---
Massieu | $\Phi_{\rm A}^{\rm qst}(t_{2})\geq\Phi_{\rm A}^{\rm neq}(t_{1})$ |  $\beta(t_{1})=\beta(t_{2})$ and $X_{\rm A}(t_{1})=X_{\rm A}(t_{2})$
Planck | $\quad\Xi^{\rm qst}_{\rm A}(t_{2})\geq\Xi_{\rm A}^{\rm neq}(t_{1})\quad$ |  $\beta(t_{1})=\beta(t_{2})$ and $\tilde{x}(t_{1})=\tilde{x}(t_{2})$
C–Entropy | $\Lambda_{\rm A}^{\rm qst}(t_{2})\geq\Lambda_{\rm A}^{\rm neq}(t_{1})$ |  $U_{\rm A}(t_{1})=U_{\rm A}(t_{2})$ and $X_{\rm A}(t_{1})=X_{\rm A}(t_{2})$
B–Entropy | $\Gamma_{\rm A}^{\rm qst}(t_{2})\geq\Gamma_{\rm A}^{\rm neq}(t_{1})$ |  $U_{\rm A}(t_{1})=U_{\rm A}(t_{2})$ and $\tilde{x}(t_{1})=\tilde{x}(t_{2})$
Table 9: Inequalities for Helmholtz–Gibbs potentials holding between the non–equilibrium state at time $t_{1}$ and the equilibrium state at time $t_{2}$, together with the conditions under which these inequalities hold.  neq Potentials | Inequalities | Condition
---|---|---
Helmholtz | $F_{\rm A}^{\rm qst}(t_{2})\leq F_{\rm A}^{\rm neq}(t_{1})$ |  $T(t_{1})=T(t_{2})$ and $X_{\rm A}(t_{1})=X_{\rm A}(t_{2})$
Gibbs | $G_{\rm A}^{\rm qst}(t_{2})\leq G_{\rm A}^{\rm neq}(t_{1})$ |  $T(t_{1})=T(t_{2})$ and $x(t_{1})=x(t_{2})$
Internal Energy | $U_{\rm A}(t_{2})\leq U_{\rm A}(t_{1})$ |  $S_{\rm A}^{\rm neq}(t_{1})=S_{\rm A}^{\rm qst}(t_{2})$ and $X_{\rm A}(t_{1})=X_{\rm A}(t_{2})$
Enthalpy | $H_{\rm A}(t_{2})\leq H_{\rm A}(t_{1})$ |  $S_{\rm A}^{\rm neq}(t_{1})=S_{\rm A}^{\rm qst}(t_{2})$ and $x(t_{1})=x(t_{2})$

Under given conditions, the non–equilibrium DL thermodynamic potentials are
smaller than the quasi-static ones. Here, we derive the inequalities between
the quasi-static and non–equilibrium thermodynamic potentials.

We consider a non–equilibrium-to-equilibrium transition $A\rightarrow 2$ from
time $t_{1}$ to $t_{2}$ described by the inverse temperature $\beta(t)$ and
the extensive variable $X_{\rm A}^{\rm neq}(t)$. The DL extensive work in this
process is evaluated as

$\displaystyle\int^{t_{2}}_{t_{1}}\left[U_{\rm
A}(t^{\prime})\frac{d\beta(t^{\prime})}{dt^{\prime}}-\tilde{x}(t^{\prime})\frac{dX_{\rm
A}^{\rm neq}(t^{\prime})}{dt^{\prime}}\right]dt^{\prime}$
$\displaystyle=\tilde{x}_{2}(X_{\rm A}^{\rm neq}(t_{2})-X_{\rm A}^{\rm
neq}(t_{1})),$ (68)

where we have assumed that $\beta(t_{1})=\beta_{2}$ and $X_{\rm A}^{\rm
neq}(t_{1})=(X_{\rm A}^{\rm qst})_{2}$, and that $\beta_{2}=\beta(t)$ and
$\tilde{x}_{2}=\tilde{x}(t)$ for time $t>t_{1}$, to make the subsystem relax
to the equilibrium state 2 at time $t_{2}$. Because $X_{\rm A}^{\rm
neq}(t_{2})=(X_{\rm A}^{\rm qst})_{2}=X_{\rm A}^{\rm neq}(t_{1})$, the right-
hand side of Eq. (D) vanishes. Thus, the net DL extensive work is zero, and so
the DL minimum work principle reduces to $0\geq(\Phi_{\rm A}^{\rm
neq})_{A\rightarrow 2}$, from which it follows that $\Phi_{\rm A}^{\rm
qst}(t_{2})\geq\Phi_{\rm A}^{\rm neq}(t_{1})$.

The inequalities for the other DL thermodynamic potentials are obtained in the
same manner. We summarize the results in Tables 8 and 9

Figure 10: Difference between the quasi-static and non–equilibrium Planck
potential in the quantum and classical cases (solid and dashed curves,
respectively) for one cycle of the Stirling engine. The blue, green, and red
curves represent the cases of weak ($A=0.5$), intermediate ($A=1.0$), and
strong ($A=1.5$) SB coupling, respectively.

In Fig. 10, to illustrate the inequalities presented in Table 8, we plot the
difference between the quasi-static and non–equilibrium Planck potential,
$\Delta\Xi_{\rm A}(t)=\Xi_{\rm A}^{\rm qst}[\beta(t),\tilde{x}(t)]-\Xi_{\rm
A}^{\rm neq}(t)$, for one cycle of the Stirling engine described in Sec.
III.4. As can be seen, $\Delta\Xi_{\rm A}(t)$ is always positive and satisfies
the inequalities in Table 8.

## References

  * Carnot (1824) S. Carnot,  _Reflexions sur la puissance motrice du feu_ (ChexzBachelier Libraire, Quai des Augustins, Paris, 1824).
  * 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 (2015) M. Esposito, M. A. Ochoa, and M. Galperin, “Nature of heat in strongly coupled open quantum systems,” Phys. Rev. B 92, 235440 (2015).
  * 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).
  * 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).
  * 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 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).
  * Tanimura (1990) Y. Tanimura, “Nonperturbative expansion method for a quantum system coupled to a harmonic-oscillator bath,” Phys. Rev. A 41, 6676–6687 (1990).
  * 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 (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).
  * 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).
  * Makri (1995) N. Makri, “Numerical path integral techniques for long time dynamics of quantum dissipative systems,” Journal of Mathematical Physics 36, 2430–2457 (1995).
  * Makri and Makarov (1995a) N. Makri and D. E. Makarov, “Tensor propagator for iterative quantum time evolution of reduced density matrices. I. Theory,” The Journal of Chemical Physics 102, 4600–4610 (1995a).
  * Makri and Makarov (1995b) N. Makri and D. E. Makarov, “Tensor propagator for iterative quantum time evolution of reduced density matrices. II. Numerical methodology,” The Journal of Chemical Physics 102, 4611–4618 (1995b).
  * Thorwart, Reimann, and Hänggi (2000) M. Thorwart, P. Reimann, and P. Hänggi, “Iterative algorithm versus analytic solutions of the parametrically driven dissipative quantum harmonic oscillator,” Phys. Rev. E 62, 5808–5817 (2000).
  * Jadhao and Makri (2008) V. Jadhao and N. Makri, “Iterative monte carlo for quantum dynamics,” The Journal of Chemical Physics 129, 161102 (2008).
  * Makri (2014) N. Makri, “Blip decomposition of the path integral: Exponential acceleration of real-time calculations on quantum dissipative systems,” The Journal of Chemical Physics 141, 134117 (2014).
  * Segal, Millis, and Reichman (2010) D. Segal, A. J. Millis, and D. R. Reichman, “Numerically exact path-integral simulation of nonequilibrium quantum transport and dissipation,” Phys. Rev. B 82, 205323 (2010).
  * Meyer, Manthe, and Cederbaum (1990) H.-D. Meyer, U. Manthe, and L. Cederbaum, “The multi-configurational time-dependent hartree approach,” Chemical Physics Letters 165, 73–78 (1990).
  * Manthe, Meyer, and Cederbaum (1992) U. Manthe, H. Meyer, and L. S. Cederbaum, “Wave‐packet dynamics within the multiconfiguration Hartree framework: General aspects and application to NOCl,” The Journal of Chemical Physics 97, 3199–3213 (1992).
  * Wang and Thoss (2007) H. Wang and M. Thoss, “Quantum dynamical simulation of electron-transfer reactions in an anharmonic environment,” The Journal of Physical Chemistry A 111, 10369–10375 (2007).
  * 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).
  * 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).
  * 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).
  * Koyanagi and Tanimura (2024a) S. Koyanagi and Y. Tanimura, “Classical and quantum thermodynamic theory described by a system–baths model: The dimensionless minimum work principle,” The Journal of Chemical Physics xxx, xxxxx (2024a), 2405.16787 .
  * 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).
  * Oono (2017) Y. Oono,  _Perspectives on Statistical Thermodynamics_ (Cambridge University Press, 2017).
  * 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).
  * Guggenheim (1986) E. A. Guggenheim,  _Thermodynamics: An Advanced Treatment for Chemists and Physicists_ (North Holland, 1986).
  * Rant (1956) Z. Rant, “Exergie, ein neues wortfur ’technishe arbeitsfahigkeit’,” (1956) p. 33.
  * Rosen and Dincer (2001) M. A. Rosen and I. Dincer, “Exergy as the confluence of energy, environment and sustainable development,” Exergy, An International Journal 1, 3–13 (2001).
  * 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).
  * 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).
  * 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).
  * Koyanagi and Tanimura (2024b) S. Koyanagi and Y. Tanimura, “Thermodynamic hierarchical equations of motion in phase space representation,” (2024b), xxxxxx:xxxxxx .