# Numerically “exact” simulations of a quantum Carnot cycle: Analysis using
thermodynamic work diagrams

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

(Last updated: )

###### Abstract

We investigate the efficiency of a quantum Carnot engine based on open quantum
dynamics theory. The model includes time-dependent external fields for the
subsystems controlling the isothermal and isentropic processes and for the
system–bath (SB) interactions controlling the transition between these
processes. Numerical simulations are conducted in a nonperturbative and non-
Markovian SB coupling regime using the hierarchical equations of motion under
these fields at different cycle frequencies. The work applied to the total
system and the heat exchanged with the baths are rigorously evaluated. In
addition, by regarding quasi-static work as free energy, we compute the
quantum thermodynamic variables and analyze the simulation results using
thermodynamic work diagrams for the first time. Analysis of these diagrams
indicates that, in the strong SB coupling region, the fields for the SB
interactions are major sources of work, while in other regions, the field for
the subsystem is a source of work. We find that the maximum efficiency is
achieved in the quasi-static case and is determined solely by the bath
temperatures, regardless of the SB coupling strength, which is a numerical
manifestation of Carnot’s theorem.

##  I Introduction

In 1824, Carnot proposed a reversible heat engine (subsequently named the
Carnot engine) that operates between a high-temperature (hot) bath at
temperature $T_{1}$ and a low-temperature (cold) bath at temperature $T_{2}$,
with gas used as a working medium.Carnot (1824) This heat engine performs the
Carnot cycle, which consists of four reversible processes: (i) isothermal
expansion, (ii) isentropic expansion, (iii) isothermal compression, and (iv)
isentropic compression. Here, heat $Q_{1}$ is absorbed from the hot bath, and
work is produced during the isothermal expansion process; heat $Q_{2}$ is
expelled to the cold bath, and some work is lost in the isothermal compression
process. The isentropic processes bridge the two isothermal processes quasi-
statically by lowering or raising the temperature of the system to $T_{2}$ or
$T_{1}$, respectively. The work is done while adiabatically reducing the
temperature, and the net heat gain is equal to the work done, i.e.,
$W=Q_{1}-Q_{2}$. Carnot proved that (a) the thermal efficiency of this engine
is the maximum that is possible and that (b) this maximum efficiency is
determined solely by the temperatures of the heat baths (Carnot’s theorem).
Moreover, in the Carnot cycle, the relationship $Q_{1}/T_{1}-Q_{2}/T_{2}=0$ is
satisfied. To elucidate the characteristic features of the Carnot cycle, a
thermodynamic work diagram (the $P$–$V$ diagram) was later introduced by
Clapeyron.Clapeyron (1834) The work done by the system is represented by the
area enclosed by the curves in the $P$–$V$ diagram.

In 1848, Thomson (later Lord Kelvin) established that the temperatures
appearing in the above formula should be regarded as thermodynamic (or
absolute) temperatures.Thomson (1848) In the 1850s, the foundations of the
second law of thermodynamics were established by Clausius on the basis of
Carnot’s result.Clausius (1850, 1856) This result was generalized by Thomson
and Clausius as $\oint dQ/T=0$, where $T$ is the temperature of the heat
source and $Q$ is the heat reversibly transferred to the system. In 1865,
Clausius introduced the concept of entropy, which is defined as $dS=dQ/T$, and
expressed the second law of thermodynamics for heat $Q$ as $\oint dQ/T<\oint
dS=0$ in any irreversible process.Clausius (1865) Entropy is an extensive
variable that characterizes the reversibility or irreversibility of thermal
processes based on the Carnot cycle, and its conjugate intensive variable is
temperature. The Carnot cycle was a key element in the construction of
thermodynamics.

Recently, the advent of nanotechnology has led to thermodynamic investigations
being extended to the quantum regime,Binder  _et al._ (2018); Kosloff (2013);
Lambert  _et al._ (2019); Härtle  _et al._ (2013); Ludovico  _et al._ (2014);
Esposito, Ochoa, and Galperin (2015a); Cohen and Galperin (2020); Brandão  _et
al._ (2015); Haughian, Esposito, and Schmidt (2018); Perarnau-Llobet  _et al._
(2018); Funo and Quan (2018); Rivas (2020); Cangemi  _et al._ (2020); Talkner
and Hänggi (2020); Zagoskin, Il’ichev, and Nori (2012); Zagoskin  _et al._
(2012); Guthrie  _et al._ (2022); Funo  _et al._ (2019); Nomura, Yoshioka, and
Nori (2021); Funo, Lambert, and Nori (2021); Lu  _et al._ (2022) in which the
work and heat to be manipulated are quantized.An  _et al._ (2015, 2016); I. A
_et al._ (2016); Vitanov  _et al._ (2017); Hoffmann  _et al._ (2018);
Bengtsson  _et al._ (2018); Josefsson  _et al._ (2019); Prete  _et al._
(2019); Ma  _et al._ (2020); de Oliveira  _et al._ (2020); Hernández-Gómez
_et al._ (2021); Ryu  _et al._ (2022) However, such studies involve
fundamental difficulties because the time evolution of the main system
(subsystem) is described by quantum mechanics, while the thermal effects of
the system are described by macroscopic thermodynamics and statistical
mechanics for the equilibrium state. For example, in an investigation of
isothermal processes, we must explicitly treat the system–bath (SB)
interactions so as to maintain the subsystem in a state of thermal equilibrium
under the influence of an external force. Because the energy scale of the SB
interactions is comparable to the energy of the subsystem, we must consider
not only the subsystem, but also the SB interactions in a quantum-mechanically
consistent manner. Moreover, the quantum description of the heat bath is
significant because the subsystem and the bath are entangled (bath
entanglement), particularly in the low-temperature case.Tanimura (2006, 2020,
2014, 2015) In an adiabatic process, a thermodynamic system with a large
number of degrees of freedom is usually assumed to be in an equilibrium state
with a quasi-static change in an extensive thermodynamic variable. An isolated
quantum subsystem, however, cannot reach a thermal equilibrium state of its
own, and thus its temperature cannot be uniquely defined during adiabatic
processes. Additionally, the definitions of thermodynamic variables in a small
subsystem, in particular extensive variables such as magnetization and strain,
are not clear in the quantum case.

In the present paper, to clarify the relationship between quantum mechanics
and statistical thermodynamics, we describe the results of quantum simulations
of a Carnot cycle on the basis of an SB model. Although such investigations
have been conducted in the framework of open quantum dynamics theories, even
for strong SB coupling cases (e.g., for heat transport,Segal, Nitzan, and
Hänggi (2003); Segal and Nitzan (2005); Segal (2006); Boudjada and Segal
(2014); Dhar (2008); Velizhanin, Wang, and Thoss (2008); Dubi and Di Ventra
(2011); Wang, Ren, and Cao (2015); Saito and Kato (2013); Song and Shi (2017);
Aurell and Montana (2019); Kilgour, Agarwalla, and Segal (2019); Xu,
Stockburger, and Ankerhold (2021); Kato and Tanimura (2015); Esposito, Ochoa,
and Galperin (2015b) heat engines and refrigerators,Quan  _et al._ (2007);
Brunner  _et al._ (2014); Chotorlishvili, Toklikishvili, and Berakdar (2011);
Kosloff and Levy (2014); Uzdin, Levy, and Kosloff (2015); Kato and Tanimura
(2016); Xu, Chen, and Liu (2018); Restrepo  _et al._ (2018); Kilgour and Segal
(2018); Das and Mukherjee (2020); Camati, Santos, and Serra (2020); Johal and
Mehta (2021); Ono  _et al._ (2020); Wiedmann, Stockburger, and Ankerhold
(2020, 2021); Xu  _et al._ (2022); Brenes  _et al._ (2020); Pancotti  _et al._
(2020); Lee  _et al._ (2020); Raja  _et al._ (2021); Ivander, Anto-Sztrikacs,
and Segal (2022) entropy production,Goyal, He, and Kawai (2020); Seshadri and
Galperin (2021); Sakamoto and Tanimura (2020) nonequilibrium cases,Schmidt
_et al._ (2015); Cangemi  _et al._ (2020) and quantum information
problems,Goold  _et al._ (2016); Strasberg  _et al._ (2017) in addition to
being applied as a context for the Jarzynski equality,Jarzynski (2004); Yukawa
(2000); Sakamoto and Tanimura (2021) fluctuation theorems,Campisi, Talkner,
and Hänggi (2009); Esposito, Harbola, and Mukamel (2009); Seifert (2012);
Campisi, Hänggi, and Talkner (2011); Aurell, Donvil, and Mallick (2020);
Whitney (2018) and Maxwell’s demonQuan  _et al._ (2006); Maruyama, Nori, and
Vedral (2009); Leff and Rex (2002)), fully quantum investigations face
difficulties because of the lack of a consistent thermodynamic formulation for
a thermal system described by a Hamiltonian with external time-dependent
perturbations.

Recently, it was shown that heat can be defined as the change in the heat-bath
energy.Kato and Tanimura (2015, 2016) The key to investigating this problem is
the hierarchical equations of motion (HEOM) formalism, which enables the
evaluation of the internal energies of not only the subsystem, but also the
bath and the SB interactions, even in low-temperature, non-Markovian, and
nonperturbative conditions.Tanimura and Kubo (1989); Tanimura (1990); Ishizaki
and Tanimura (2005); Tanimura (2006, 2014, 2015, 2020) We show that, by
introducing a time-dependent SB interaction, thermal transitions between
isothermal and adiabatic states can also be investigated. By using HEOM to
calculate the change in work and bath energy applied to the entire system, the
efficiency of the heat engine under an arbitrary time-dependent external field
can then be evaluated without assumptions.Koyanagi and Tanimura (2022)

To enable thermodynamic analysis, the thermodynamic variables can be treated
based on the SB model. For this purpose, we adopt the minimal work principle
for a total isolated system, expressed as $W_{\rm tot}(t)\geq\Delta F_{\rm
tot}(t)$, where $W_{\rm tot}(t)$ is the work done by external fields and
$\Delta F_{\rm tot}(t)$ is the change of free energy, and define the “quasi-
equilibrium” Helmholtz energy (qHE) as $\Delta F_{\rm tot}(t)=W_{\rm tot}^{\rm
qst}(t)$ with the system driven quasi-statically by external fields. Although
$W_{\rm tot}(t)$ cannot be evaluated within the framework of regular open
quantum dynamics theories (because the number of degrees of freedom of the
bath has been reduced), we can evaluate this quantity indirectly using the
hierarchical elements in the HEOM formalism. This is because, in the HEOM
formalism, the higher hierarchical elements store information about the higher
cumulant of the bath coordinates, as previously demonstrated.Kato and Tanimura
(2015, 2016); Song and Shi (2017); Sakamoto and Tanimura (2020, 2021) We then
show that the Kelvin–Planck statement (or heat engine statement) of the second
law of thermodynamics is only validated when we introduce a time-dependent SB
interaction that describes energy conservation for adiabatic
transitions.Koyanagi and Tanimura (2022)

To this end, we extend previous research to verify Carnot’s theorem using the
HEOM approach. The Carnot engine is described as a two-level subsystem (A)
coupled with harmonic heat baths (B). To control the isothermal process, we
consider the external field for the system (the isothermal driving field) and,
to turn the hot and cold baths on and off, we introduce two time-dependent SB
interactions (the adiabatic transition fields). We then investigate the
thermodynamic efficiency of the Carnot cycle under various physical
conditions. In addition, using the qHE formalism, thermodynamic work diagrams
of external forces (such as stresses) and their conjugate variables (such as
strains), similar to Clapeyron’s $P$–$V$ diagram,Clapeyron (1834) are
introduced to analyze the work done in the system. This extension is useful
for analyzing experimental results in the quantum regime, where the quantized
work and heat are to be manipulated.An  _et al._ (2015, 2016); Vitanov  _et
al._ (2017); Hoffmann  _et al._ (2018); Bengtsson  _et al._ (2018); Josefsson
_et al._ (2019); Prete  _et al._ (2019); Ma  _et al._ (2020); de Oliveira  _et
al._ (2020); Hernández-Gómez  _et al._ (2021)

The remainder of this paper is organized as follows. In Sec. II, we introduce
the SB Hamiltonian and the HEOM formalism. We then describe the scheme used to
calculate various thermodynamic variables based on open quantum dynamics
theory. In Sec. III, we explain our model of the quantum Carnot cycle and
present simulation results under the isothermal driving and adiabatic
transition fields. The conjugate properties of these variables are calculated
by regarding quasi-static work as free energy. Thermodynamic work diagrams are
presented as functions of these variables. Finally, Sec. IV presents
concluding remarks.

##  II Theory

###  II.1 Model

We consider a subsystem A coupled to two heat baths at high and low inverse
temperatures $\beta_{1}=1/k_{B}T_{1}$ and $\beta_{2}=1/k_{B}T_{2}$ as the heat
sources, where $k_{B}$ is the Boltzmann constant. The total Hamiltonian is
expressed as

$\displaystyle\hat{H}_{\rm tot}(t)=\hat{H}_{\rm
A}(t)+\sum_{k=1}^{2}\left(\hat{H}_{\rm I}^{k}(t)+\hat{H}_{\rm B}^{k}\right),$
(1)

where $\hat{H}_{\rm A}(t)$, $\hat{H}_{\rm I}^{k}(t)$, and $\hat{H}_{\rm
B}^{k}$ are the Hamiltonians of the system, $k$th SB interaction, and $k$th
bath, respectively. We consider a two-level system (TLS) defined as

$\hat{H}_{\rm A}(t)=-B(t)\hat{\sigma}_{z}+E\hat{\sigma}_{x},$ (2)

where $B(t)$ is the isothermal driving field (IDF), $E$ is the off-diagonal
coupling parameter, and $\hat{\sigma}_{\alpha}$ ($\alpha=x$, $y$, or $z$) are
Pauli matrices. In the case of a spin system, $B(t)$ corresponds to the
longitudinal magnetic field and $E$ is the transverse electric (Stark) field.
The Hamiltonian representing the $k$th SB interaction and the $k$th bath are
given byKoyanagi and Tanimura (2022)

$\displaystyle\hat{H}_{\rm
I}^{k}(t)=A_{k}(t)\hat{V}_{k}\sum_{j}c_{j}^{k}\left[\hat{b}_{j}^{k}+(\hat{b}_{j}^{k})^{\dagger}\right]$
(3)

and

$\displaystyle\hat{H}_{\rm
B}^{k}=\sum_{j}\hbar\omega_{j}^{k}\left[(\hat{b}_{j}^{k})^{\dagger}\hat{b}_{j}^{k}+\frac{1}{2}\right],$
(4)

respectively, where $\hat{V}_{k}$ is the system operator that describes the
coupling to the $k$th bath and $A_{k}(t)$ is the adiabatic transition field
(ATF), which is introduced to describe the operation of an adiabatic wall
between the system and the $k$th heat bath (e.g., the insertion or removal of
the adiabatic wall or attaching or detaching the quantum system to/from the
bath). Here, $\omega_{j}^{k}$, $c_{j}^{k}$, $\hat{b}_{j}^{k}$, and
$(\hat{b}_{j}^{k})^{\dagger}$ are the frequency, coupling strength, and
annihilation and creation operators for the $j$th mode of the $k$th bath,
respectively.

Due to the bosonic nature of the baths, all bath effects on the system are
determined by the $k$th bath correlation function,
$C_{k}(t)\equiv\langle\hat{X}_{k}(t)\hat{X}_{k}(0)\rangle_{\rm B}$, where
$\hat{X}_{k}\equiv\sum_{j}c_{j}^{k}[\hat{b}_{j}^{k}+(\hat{b}_{j}^{k})^{\dagger}]$
is the collective coordinate of the $k$th bath and $\langle\cdots\rangle_{\rm
B}$ represents the average taken with respect to the canonical density
operator of the baths. The bath correlation function is expressed as

$\displaystyle
C_{k}(t)=\int_{0}^{\infty}d\omega\,\frac{J_{k}(\omega)}{\pi}\left[\coth\left(\frac{\beta_{k}\hbar\omega}{2}\right)\cos(\omega
t)-i\sin(\omega t)\right],$ (5)

where
$J_{k}(\omega)\equiv\pi\sum_{j}(c_{j}^{k})^{2}\delta(\omega-\omega_{j}^{k})$
is the bath spectral density. The real part of Eq. (5) is analogous to the
classical correlation function of the bath and corresponds to the
fluctuations, while the imaginary part corresponds to the dissipation. The
fluctuation term is related to the dissipation term through the quantum
version of the fluctuation–dissipation theorem. Tanimura and Kubo (1989);
Tanimura (2006) In this paper, we use the Drude spectral density function,
described as

$J_{k}(\omega)=\frac{\hbar\gamma_{k}^{2}\omega}{\gamma_{k}^{2}+\omega^{2}},$
(6)

where $\gamma_{k}$ is the inverse noise correlation of the $k$th bath.

###  II.2 Hierarchical equations of motion

In the HEOM formalism, the set of equations of motion consists of the
auxiliary density operators (ADOs).Tanimura and Kubo (1989); Tanimura (1990);
Ishizaki and Tanimura (2005); Tanimura (2006, 2014, 2015, 2020); Koyanagi and
Tanimura (2022) Here, we consider the case in which the bath correlation
function, Eq. (5), is written as a linear combination of exponential
functions, $C_{k}(t)=\sum_{l=0}^{K_{k}}\zeta_{l}^{k}e^{-\nu_{l}^{k}|t|}$,
where $\nu_{l}^{k}$, $\zeta_{l}^{k}$, and $K_{k}$ are the frequency, strength,
and cutoff integer value for the $k$th bath obtained from a Padé spectral
decomposition scheme to reduce the hierarchy size.Hu, Xu, and Yan (2010) The
ADOs introduced in the HEOM are defined by $\hat{\rho}_{\vec{n}}(t)$ with a
set of indices
$\vec{n}=(n_{1}^{0},\ldots,n_{1}^{K_{1}},n_{2}^{0},\ldots,n_{2}^{K_{2}})$,
where $n_{l}^{k}$ represents an integer value of zero or above. The zeroth
ADO, $\hat{\rho}_{\vec{0}}(t)$ with $\vec{0}=(0,0,\ldots,0)$, corresponds to
the actual reduced density operator. The HEOM for the IDF and ATFs are then
expressed asKoyanagi and Tanimura (2022)

$\displaystyle\frac{\partial}{\partial t}\hat{\rho}_{\vec{n}}(t)$
$\displaystyle=\left(-\frac{i}{\hbar}\hat{H}_{\rm
A}^{\times}(t)-\sum_{k=1}^{2}\sum_{l=0}^{K_{k}}n_{l}^{k}\nu_{l}^{k}\right)\hat{\rho}_{\vec{n}}(t)$
$\displaystyle-\frac{i}{\hbar}\sum_{k=1}^{2}{A_{k}(t)}\sum_{l=0}^{K_{k}}n_{l}^{k}\hat{\Theta}_{l}^{k}\hat{\rho}_{\vec{n}-\vec{e}_{l}^{k}}(t)$
$\displaystyle-\frac{i}{\hbar}\sum_{k=1}^{2}{A_{k}(t)}\sum_{l=0}^{K_{k}}\hat{V}_{k}^{\times}\hat{\rho}_{\vec{n}+\vec{e}_{l}^{k}}(t),$
(7)

where $\vec{e}_{l}^{~{}k}$ is the ($K_{k}$+1)-dimensional unit vector. We
introduce a set of fluctuation–dissipation operators as

$\hat{\Theta}_{0}^{k}\equiv\left(\frac{\gamma_{k}}{\beta_{k}}+\sum_{m=1}^{K_{m}}\frac{\zeta_{m}\gamma_{k}^{2}}{\beta_{k}}\frac{2\gamma_{k}}{\gamma_{k}^{2}-\nu_{m}^{2}}\right)\hat{V}_{k}^{\times}-\frac{i\hbar\gamma_{k}^{2}}{2}\hat{V}_{k}^{\circ}$
(8)

and

$\displaystyle\hat{\Theta}_{l}^{k}\equiv-\frac{\zeta_{l}\gamma_{k}^{2}}{\beta_{k}}\frac{2\nu_{l}}{\gamma_{k}^{2}-\nu_{l}^{2}}\hat{V}_{k}^{\times},$
(9)

where
$\hat{\mathcal{O}}^{\times}\hat{\mathcal{P}}=[\hat{\mathcal{O}},\hat{\mathcal{P}}]$
and
$\hat{\mathcal{O}}^{\circ}\hat{\mathcal{P}}=\\{\hat{\mathcal{O}},\hat{\mathcal{P}}\\}$
for arbitrary operators $\hat{\mathcal{O}}$ and $\hat{\mathcal{P}}$, and
$\nu_{l}$ and $\zeta_{l}$ are the frequency and strength, respectively.

As the temporal initial conditions, we consider the factorized initial case

$\hat{\rho}_{\rm tot}(0)=\hat{\rho}_{\rm
A}(0)\prod_{k=1}^{2}\frac{e^{-\beta_{k}\hat{H}_{\rm B}^{k}}}{{\rm tr}_{\rm
B}\\{e^{-\beta_{k}\hat{H}_{\rm B}^{k}}\\}},$

where $\hat{\rho}_{\rm A}(t)=\hat{\rho}_{\vec{0}}(t)$ is the reduced density
operator of the subsystem. To obtain the bath-entangled steady-state solution
described as $\hat{\rho}_{\vec{n}}(t)$, we integrate the HEOM under the
periodical external fields until all of the hierarchy elements reach a steady
state.

Using the zeroth member of the hierarchy $\hat{\rho}_{\vec{0}}(t)$, the
nonequilibrium internal energy of the system (or the expectation value of the
system energy) at time $t$ is evaluated as

$\displaystyle U_{\rm A}^{\rm neq}(t)={\rm tr}_{\rm A}\left\\{\hat{H}_{\rm
A}(t)\hat{\rho}_{\vec{0}}(t)\right\\}.$ (10)

The nonequilibrium internal energy of the $k$th SB interaction is then
expressed asKato and Tanimura (2015, 2016); Sakamoto and Tanimura (2020, 2021)

$\displaystyle U_{{\rm I}_{k}}^{\rm neq}(t)=A_{k}(t)\sum_{l=0}^{K_{k}}{\rm
tr}_{\rm A}\left\\{\hat{V}_{k}\hat{\rho}_{\vec{e}_{l}^{k}}(t)\right\\},$ (11)

where $\vec{e}_{l}$ is the index for the first-order hierarchical member.

Using the HEOM formalism, we can evaluate the $k$th bath energy asSakamoto and
Tanimura (2020); Koyanagi and Tanimura (2022)

$\displaystyle\frac{\partial}{\partial t}U_{{\rm B}_{k}}^{\rm
neq}(\beta_{k};t)$ $\displaystyle=A_{k}(t)\sum_{l=0}^{K_{k}}\nu_{l}{\rm
tr}\\{\hat{V}\hat{\rho}_{\vec{e}_{l}}(t)\\}$
$\displaystyle+A_{k}^{2}(t)\gamma^{2}{\rm tr}_{\rm
A}\left\\{\hat{V}^{2}\hat{\rho}_{\vec{0}}(t)\right\\}.$ (12)

Because there is no external force applied to the bath, the $k$th bath heat
can be obtained by integrating the above equation as

$\displaystyle Q_{{\rm B}_{k}}(\beta_{k};t)=\Delta U_{{\rm B}_{k}}^{\rm
neq}(\beta_{k};t)$ (13)

with $\Delta U_{{\rm B}_{k}}^{\rm neq}(\beta_{k};t)\equiv U_{{\rm B}_{k}}^{\rm
neq}(\beta_{k};t+\delta t)-U_{{\rm B}_{k}}^{\rm neq}(\beta_{k};t)$, where
$\delta t$ is the time step of the thermal process.

In this study, work is defined as the change in energy from one state to
another state under a time-dependent perturbation expressed as

$\displaystyle W_{\rm tot}(t)=\int_{t}^{t+\delta
t}dt^{\prime}{\rm{tr}}\left\\{\frac{d\hat{H}_{\rm
tot}(t^{\prime})}{dt^{\prime}}\hat{\rho}_{\rm tot}(t^{\prime})\right\\}.$ (14)

Because there is no external field for the bath, we have $W_{\rm
tot}(t)=W_{{\rm A}+{\rm I}}(t)$, where

$W_{{\rm A}+{\rm I}}(t)\equiv W_{\rm A}(t)+W_{{\rm I}_{1}}(t)+W_{\rm
I_{2}}(t).$ (15)

Using the HEOM, each component is evaluated as

$\displaystyle W_{\rm A}(t)=\int_{t}^{t+\delta t}dt^{\prime}{\rm tr}_{\rm
A}\left\\{\frac{\partial{\hat{H}}_{\rm A}(t^{\prime})}{\partial
t}\hat{\rho}_{\vec{0}}(t^{\prime})\right\\}$ (16)

and

$\displaystyle W_{{\rm I}_{k}}(t)=\int_{t}^{t+\delta
t}dt^{\prime}\frac{dA_{k}(t^{\prime})}{d{t^{\prime}}}\sum_{l=0}^{K_{k}}{\rm
tr}_{\rm A}\\{\hat{V}_{k}\hat{\rho}_{\vec{e}_{l}^{k}}(t^{\prime})\\}.$ (17)

By this definition, when $W_{\rm tot}(t)>0$, work is done from the outside to
the total system and is called “positive work.”

Using the above, we can evaluate the $k$th bath energy for $A_{k}(t)>0$ more
accurately than from Eq. (12) as

$\displaystyle\frac{\partial}{\partial t}U^{\rm neq}_{\rm B_{k}}(\beta_{k};t)$
$\displaystyle=-\frac{\partial}{\partial t}\left[U^{\rm neq}_{\rm A}(t)+U^{\rm
neq}_{\rm I_{k}}(t)\right]$ $\displaystyle+\frac{\partial}{\partial
t}\left[W_{\rm A}(t)+W_{\rm I_{k}}(t)\right].$ (18)

###  II.3 Quasi-equilibrium Helmholtz energy and thermodynamic variables

Using the above definition of work with the HEOM, the efficiency of the system
driven by the IDF and ATFs can be evaluated numerically and rigorously for any
cycle speed. Nevertheless, we concentrate our discussion to the quasi-static
case and attempt to quantify the quantum thermodynamic variables by comparing
the physical quantities calculated in the HEOM with the thermodynamic
quantities evaluated in the qHE, as explained below.

The work for the total system satisfies the minimum work principle, expressed
asLenard (1978); Tasaki (2000); Oono (2017)

$\displaystyle W_{\rm tot}(\beta;t)\geq\Delta F_{{\rm tot}}(\beta;t),$ (19)

which corresponds to the second law of thermodynamics. For the SB model, the
above inequality has been verified numerically by placing the real-time HEOM
for work on the left-hand side and the imaginary-time HEOM for the free energy
on the right-hand side.Tanimura (2014, 2015); Zhang and Tanimura (2022) This
is because, even when ${\hat{H}}_{A}$ and ${\hat{H}}_{I}$ are weakly dependent
on time, the steady-state solution of the HEOM is represented as the reduced
density operator of the total equilibrium state, expressed as ${\rm tr}_{\rm
B}\\{\exp[-\beta({\hat{H}}_{\rm A}+{\hat{H}}_{\rm I}+{\hat{H}}_{\rm
B})]\\}$.Tanimura (2014, 2015) It is known that, for a quasi-static process,
the equality can be expressed within the level of numerical accuracy
asSakamoto and Tanimura (2020, 2021); Koyanagi and Tanimura (2022)

$\displaystyle\Delta F_{{\rm tot}}^{\rm qst}(\beta;t)=W_{\rm tot}^{\rm
qst}(\beta;t).$ (20)

While the above equation holds for the total system, rather than the reduced
system, we can utilize the qHE using the HEOM formalism because, through the
use of the HEOM elements, we can evaluate the total work of the system, as
illustrated in Eqs. (14)–(17). For each component, we also have $\Delta
F_{\alpha}^{\rm qst}(\beta;t)=W_{\alpha}^{\rm qst}(\beta;t)$ for $\alpha={\rm
A}$, ${\rm I}_{1}$, and ${\rm I}_{2}$, where $W_{\alpha}^{\rm qst}(\beta;t)$
is defined by Eqs. (16) and (17) using the quasi-static solution of the HEOM,
which is expressed as $\hat{\rho}_{\vec{n}}^{\rm qst}(t)$. Koyanagi and
Tanimura (2022)

In the Carnot cycle, the subsystem interacts with a single bath at a time. For
external fields, such as the magnetic field $B(t)$ and stress $A_{k}(t)$, we
then have that $\Delta F_{\rm tot}^{\rm qst}(\beta;t)=\Delta F_{\rm A+I}^{\rm
qst}(\beta;t)$ for $\beta=\beta_{k}$ with $k=1$ or 2, where $\Delta F_{\rm
A+I}^{\rm qst}(\beta_{k};t)\equiv\Delta F_{\rm A}^{\rm
qst}(\beta_{k};t)+\Delta F_{\rm I_{k}}^{\rm qst}(\beta_{k};t)$. Then, we can
define the conjugate variables, such as the magnetization and strain,
asKoyanagi and Tanimura (2022)

$M^{\rm qst}(\beta_{k};t)\equiv-\frac{\partial\Delta F_{\rm A}^{\rm
qst}(\beta_{k};t)}{\partial{B(t)}}$ (21)

and

$D_{k}^{\rm qst}(\beta_{k};t)\equiv-\frac{\partial\Delta F_{{\rm I}_{k}}^{\rm
qst}(\beta_{k};t)}{\partial{A_{k}(t)}}.$ (22)

Thus, we can evaluate the above variables in terms of the state described by
the ADOs at time $t$ as

$M^{\rm qst}(\beta_{k};t)={\rm tr}_{\rm
A}\left\\{\hat{\sigma}_{z}\hat{\rho}_{\vec{0}}(t)\right\\}$ (23)

and

$D_{k}^{\rm qst}(\beta_{k};t)=-\sum_{l=0}^{K_{k}}{\rm tr}_{\rm
A}\left\\{\hat{V}_{k}\hat{\rho}_{\vec{e}^{k}_{l}}(t)\right\\}.$ (24)

The change in the quasi-equilibrium Boltzmann entropy from time $t$ to
$t+\delta t$ can also be evaluated by differentiating the qHE with respect to
the $k$th bath inverse temperature $\beta_{k}$ asSakamoto and Tanimura (2020)

$\displaystyle\Delta S_{\rm A+I}^{\rm qst}(\beta_{k};t)\equiv k_{\rm
B}\beta_{k}^{2}\frac{\partial\Delta F^{\rm qst}_{\rm
A+I}(\beta_{k};t)}{\partial\beta_{k}},$ (25)

where we consider the case in which the subsystem interacts with the $k$th
bath at that time. A comparison of this definition with the von Neumann
entropy using the simulation results for the Carnot cycle is given in Appendix
B.

Note that $M^{\rm qst}(\beta;t)$, $D_{k}^{\rm qst}(\beta;t)$, and $\Delta
S_{\rm A+I}^{\rm qst}(\beta;t)$ are state variables in the quasi-static state
because they are uniquely determined by the state specified by the quasi-
equilibrium distribution at $t$ and are independent of the pathway of work.
When the subsystem consists of $n$ noninteracting spins that are independently
coupled to the heat bath, the magnitudes of the above variables are
proportional to $n$. Thus, they are extensive properties, while $B(t)$,
$A_{k}(t)$, and $\beta$ are intensive properties.

The quasi-static isothermal process at temperature $T_{k}=1/k_{B}\beta_{k}$ is
then expressed as

$\displaystyle d\Delta F^{\rm qst}_{\rm A+I}(\beta_{k};t)$
$\displaystyle=-\Delta S_{\rm A+I}^{\rm qst}(\beta_{k};t)dT_{k}-M^{\rm
qst}(\beta_{k};t)dB$ $\displaystyle-D_{k}^{\rm qst}(\beta_{k};t)dA_{k}.$ (26)

##  III Quantum Carnot cycle

The conventional Carnot cycle consists of a system, described by pressure $P$
and volume $V$, that interacts with hot and cold baths. Although isothermal
processes can be easily described by changing an external perturbation $B(t)$,
there are several difficulties in realizing the Carnot cycle, even
theoretically, in a nanoscale quantum system. For example, to conduct a
microscopic investigation of the Carnot cycle, the thermal work done by the
insertion or removal of the adiabatic wall must be included, because otherwise
the energy conservation of the total system would be violated. Accordingly,
because a small isolated system cannot reach thermal equilibrium by itself, an
isentropic process cannot be achieved spontaneously by turning off the heat
bath. Finally, although the Carnot cycle has been commonly characterized by
the $P$–$V$ diagram, in a small quantum system, the definitions of extensive
variables are unclear. As we explained, the qHE provides the necessary means
for analyzing the simulation results of the quantum Carnot cycle obtained from
the HEOM approach.

###  III.1 Case under factorized assumption

Figure 1:  (a) $B$–$M$ diagram of Carnot engine driven by $B(t)$ with
$B_{0}=2$ and $E=0.5$ evaluated under the factorized assumption, i.e., $F_{\rm
A}(t)=-\ln{\rm tr}_{\rm A}\\{\exp[-\beta H_{\rm A}(t)]\\}/\beta$. Here, the
blue and red curves represent the cold ($\beta_{2}=1.0$) isothermal expansion
and hot ($\beta_{1}=0.5$) isothermal compression processes, while the purple
and orange curves represent the hot isentropic expansion and cold isentropic
compression processes, respectively. The cycle starts from the red arrow and
evolves in a counterclockwise fashion over time (heat engine). (b) $T$–$S$
diagram under the same condition as (a).

Before presenting the computational results, we will first illustrate the
$B$–$M$ diagram under the assumption that the free energy of the TLS is
determined by the partition function of the isolated TLS with Markovian and
extremely weak SB interactions (factorized assumption). Thus, from $F_{\rm
A}^{0}(\beta;t)=-\ln{\rm tr}\\{\exp[-\beta\hat{H}_{\rm A}(t)]\\}/\beta$, we
have

$F_{\rm
A}^{0}(\beta;t)=-\frac{1}{\beta}\ln\left[2\cosh\left(\beta\sqrt{B^{2}(t)+E^{2}}\right)\right].$
(27)

The conjugate variable $M_{\rm A}^{0}(\beta_{k};t)\equiv-{\partial F_{\rm
A}^{0}(\beta_{k};t)}/{\partial B}$ at the inverse temperature $\beta_{k}$ is
then evaluated as

$M_{\rm
A}^{0}(\beta_{k};t)=\frac{B(t)}{\sqrt{B^{2}(t)+E^{2}}}\tanh\left(\beta_{k}\sqrt{B^{2}(t)+E^{2}}\right).$
(28)

The entropy can also be evaluated by differentiating $F^{0}_{\rm A}(\beta;t)$
with respect to $\beta$ as

$\begin{split}S^{0}_{\rm A}(\beta;t)&=k_{\rm
B}\left\\{\ln\left[2\cosh\left(\beta\sqrt{B^{2}(t)+E^{2}}\right)\right]\right.\\\
&\quad\quad\left.-\beta\sqrt{B^{2}(t)+E^{2}}\tanh(\beta\sqrt{B^{2}(t)+E^{2}})\right\\}.\end{split}$
(29)

This value agrees with the von Neumann entropy defined by Eq. (40), because
the subsystem is isolated.

Under the factorized assumption, we can depict the $B$–$M$ diagram for the
isothermal process $B(t):$ (i) $B_{1}\rightarrow B_{2}$ and (iii)
$B_{3}\rightarrow B_{4}$ for any form of $B(t)$. To set up the isentropic
processes (ii) $B_{2}\rightarrow B_{3}$ and (iv) $B_{4}\rightarrow B_{1}$, we
must satisfy the condition (see Appendix A)

$B_{b}=\sqrt{\frac{\beta_{a}}{\beta_{b}}\left(B_{a}^{2}+E^{2}\right)-E^{2}},$
(30)

where $\beta_{a}$ and $\beta_{b}$ are the inverse temperatures of the baths
before and after the adiabatic process for $(a,b)=({\rm ii},{\rm iii})$ or
$({\rm iv},{\rm i})$; otherwise, entropy would be generated. In the
conventional Carnot cycle, the isentropic process is regarded as a
temperature-changing process in which the internal energy changes
spontaneously under adiabatic conditions. In the present case, the subsystem
is isolated and the internal energy does not change without an external field.
Moreover, the definition of the temperature during the isentropic process is
not clear, because the system is microscopic and isolated except for the
initial and final states at $\beta_{a}$ and $\beta_{b}$. Thus, instead of
assuming a spontaneous change, we actively control the external field (IDF)
using Eq. (30) to realize the temperature change between $\beta_{a}$ and
$\beta_{b}$ in the isentropic processes.

As in the Carnot cycle, the parameters $B_{1}$ and $B_{2}$ can take any
values, but they must satisfy the relation $B_{1}<B_{2}$ to ensure that the
work done by the cycle is positive. In Fig. 1, we depict the $B$–$M$ diagram
for $B_{1}=\sqrt{15}/2$ and $B_{2}=\sqrt{7}/2$ with $\beta_{\rm 1}/\beta_{\rm
2}=2$. In comparison with the $P$–$V$ diagram for an ideal gas, the $B$–$M$
diagram has the opposite rotational direction. This is because an ideal gas is
described by $dU=TdS-PdV$, whereas here we have $dU=TdS+BdM$ (see Eqs.
(33)-(36)). Thus, this cycle is described as (i) isothermal compression, (ii)
isentropic compression, (iii) isothermal expansion, and (iv) isentropic
expansion, which is the opposite of the $P$–$V$ case. The area enclosed by the
curves corresponds to the work, but a counterclockwise cycle represents
positive work, which is also the opposite of the $P$–$V$ case.

###  III.2 Cases for nonperturbative and non-Markovian conditions

Table 1: Cycle of the IDF [$B(t)$] and the hot and cold ATFs [$A_{1}(t)$ and
$A_{2}(t)$] for an eight-stroke Carnot engine. Because we explicitly treat the
bath removal and bath attachment processes of hot and cold baths denoted by
(i′)–(vi′), the cycle consists of eight strokes. Here, the stroke intervals
are equally spaced and described as $\tau$. Thus, the cycle period is
$T=8\tau$. The parameters $B_{1}$ and $B_{2}$ can take any values, while
$B_{3}$ and $B_{4}$ are determined from Eq. (30). We set $B_{1}=\sqrt{15}/2$,
$B_{2}=\sqrt{7}/2$, $B_{3}=1/2$, and $B_{4}=\sqrt{3}/2$.
| $B(t)$ | $A_{1}(t)/|A_{1}|$ | $A_{2}(t)/|A_{2}|$
---|---|---|---
(i) | $B_{1}+(B_{2}-B_{1})t/\tau$ | $1$ | $0$
(i′) | $B_{2}$ | $1-({t-\tau})/{\tau}$ | $0$
(ii) | $B_{2}+(B_{3}-B_{2})({t-2\tau})/{\tau}$ | $0$ | $0$
(ii′) | $B_{3}$ | $0$ | $({t-3\tau})/{\tau}$
(iii) | $B_{3}+(B_{4}-B_{3})(t-4\tau)/{\tau}$ | $0$ | $1$
(iii′) | $B_{4}$ | $0$ | $1-({t-5\tau})/{\tau}$
(iv) | $B_{4}+(B_{1}-B_{4})({t-6\tau})/{\tau}$ | $0$ | $0$
(iv′) | $B_{1}$ | $({t-7\tau})/{\tau}$ | $0$
Figure 2: Time profiles of $B(t)$, $A_{1}(t)$, and $A_{2}(t)$ for the eight-
stroke Carnot cycle depicted as functions of $\Omega t$, where
$\Omega=\pi/4\tau$. The profile of $B(t)$ consists of the red, purple, blue,
and orange solid lines, corresponding to isothermal and isentropic processes
(i)–(iv), and the red, purple, blue, and orange dashed lines, corresponding to
the bath removal and bath attachment processes of hot and cold baths
(i′)–(vi′), respectively. The profiles of $A_{1}(t)$ and $A_{2}(t)$ consist of
the orange dashed, red solid, and red dashed lines, and the purple dashed,
blue solid, and blue dashed lines, corresponding to the bath attachment,
operation, and removal processes (vi′), (i), and (i′), and (ii′), (iii), and
(iii′), for the hot and cold baths, respectively.

####  III.2.1 Simulation details

We now simulate the Carnot engine using our model, which includes both the IDF
and the ATFs. We compute $W_{\rm tot}(t)=W_{\rm A+I}(t)$ for various cycle
frequencies $\Omega$, and then analyze the characteristics of the engine in
the quasi-static case using thermodynamic variables defined by the qHE,
$\Delta F_{\rm A+I}^{\rm qst}(t)=W_{\rm A+I}^{\rm qst}(t)$. As explained in
the factorized case, to realize the isentropic processes, we must change the
IDF after removing or before attaching the heat bath, as described by (ii) and
(iv). In addition to the regular four Carnot processes (i)–(iv), we must
explicitly treat four additional processes that represent (i′) removing and
(iv′) attaching the hot bath and (ii′) attaching and (iii′) removing the cold
bath.

We conduct numerical simulations for the TLS coupled to two bosonic baths at
different temperatures using the HEOM approach. Here, we consider the equally
spaced stroke period $\tau$. Thus, the cycle period and frequency are
$T=8\tau$ and $\Omega=\pi/4\tau$, respectively. (See Table 1 and Fig. 2.)
Throughout this paper, we fix the system parameter as $E=0.5$ and the bath
parameters as $\beta_{1}=0.5$, $\beta_{2}=1.0$,
$\gamma\equiv\gamma_{1}=\gamma_{2}=1.0$, and $A_{0}=0.1,\sqrt{0.1}$, or
$\sqrt{3.0}$, where $A_{0}\equiv A_{1}=A_{2}$ is the maximum strength of
$A_{1}(t)$ and $A_{2}(t)$. We choose the truncation number of the hierarchy,
representing the depth of the HEOM computation, as $N=8$ for
$A_{0}=\sqrt{3.0}$ and $N=6$ for the other cases. A Padé spectral
decomposition scheme is employed to obtain the expansion coefficients of the
noise correlation functions.Hu, Xu, and Yan (2010) We set the number of Padé
frequencies to $K_{1}=K_{2}=4$. We integrate Eq. (7) from the temporal initial
state until the simulation cycle reaches the steady state with the time-
dependent fields using the fourth-order Runge–Kutta method with a time step of
$\delta t=1.0\times 10^{-2}$. The various thermodynamic variables can then be
evaluated from the HEOM elements.

To reduce the computation time, we set $\hat{\rho}_{\vec{n}}$ to zero for the
elements $N_{1}>1$ and $N_{2}>1$ in $\tau\geq 100$, where
$N_{k}=\sum_{l=0}^{K_{k}}n_{l}^{k}$. This treatment is valid because the ADOs
of the hot bath are negligibly small by the time the bath is restored to the
system; the ADOs decay exponentially in time after the baths are removed from
the system.

####  III.2.2 Thermal efficiency

Figure 3: (i) Efficiencies of the quantum Carnot engine calculated as a
function of $\Omega$ for different SB coupling strengths: (a) $A_{0}=0.1$
(blue curve with square markers), (b) $A_{0}=\sqrt{0.1}$ (green curve with
circular markers), and (c) $A_{0}=\sqrt{3.0}$ (red curve with triangular
markers). (ii) Efficiencies of the quantum Carnot cycle as a function of
$A_{0}$ at the frequencies $\Omega/2\pi=1.25\times 10^{-3}$ (blue dashed
curve), $2.5\times 10^{-4}$ (green dashed curve), and $6.25\times 10^{-5}$
(red dashed curve) marked with black circular data points. The corresponding
$\Omega$ is indicated by the dashed line for each frequency in (i).

We now investigate the efficiency of the engine, which is significant in the
characterization of the Carnot cycle. The total work and heat from the hot
bath per cycle, $W_{\rm tot}$ and $Q_{{\rm B}_{1}}$, are obtained from Eqs.
(15)–(17) and Eqs. (12) and (13). The efficiency is then evaluated as

$\displaystyle\varepsilon=\frac{W_{\rm tot}}{Q_{{\rm B}_{1}}}.$ (31)

Figure 3(i) displays the efficiency $\varepsilon$ as a function of $\Omega$
for (a) the weak ($A_{0}=0.1$), (b) intermediate ($A_{0}=\sqrt{0.1}$), and (c)
strong ($A_{0}=\sqrt{3.0}$) SB coupling cases. For any $A_{0}$, the efficiency
converges to the Carnot limit $\varepsilon_{\rm C}=1-\beta_{\rm H}/\beta_{\rm
C}$ (which is $0.5$ in the present case) for the quasi-static case
$\Omega\rightarrow 0$, providing a manifestation of Carnot’s theorem.

Figure 3(ii) displays the efficiency as a function of the SB coupling strength
$A_{0}$ for the three values of $\Omega$ indicated by the colored dashed lines
in Fig. 3(i). When $A_{0}$ is small, the thermal excitation is weak and a
longer time is required for equilibration. Conversely, when $A_{0}$ is very
large, the efficiency decreases due to the strong relaxation: the heat bath
suppresses the kinetic motion relating to the heat flow.Kato and Tanimura
(2015, 2016); Sakamoto and Tanimura (2020) Thus, for a fixed value of
$\Omega$, the efficiency reaches a maximum in the intermediate region where
$A_{0}$ is neither large nor small. This feature has been observed in various
quantum transport problems, such as exciton transferIshizaki and Fleming
(2009) and chemical-reaction problems,Tanimura and Wolynes (1992) and is known
in the classical case as the Kramers turnover problem. As $\Omega$ becomes
smaller and approaches the quasi-static limit, the $A_{0}$ dependence of the
efficiency is suppressed and the results approach the Carnot limit, as can be
seen from Fig. 3(i).

####  III.2.3 B–M and A–D diagrams

Figure 4:  (1) $B$–$M$, (2) $A_{1}$–$D_{1}$, and (3) $A_{2}$–$D_{2}$ diagrams
for the quasi-static processes with $E=0.5$ in (a) the weak SB coupling case
($A_{0}=0.1$), (b) the intermediate SB coupling case ($A_{0}=\sqrt{0.1}$), and
(c) the strong SB coupling case ($A_{0}=\sqrt{3.0}$). In (1), the $B$–$M$
diagrams, the red and blue curves represent the (i) hot ($\beta_{1}=0.5$)
isothermal compression and (iii) cold ($\beta_{2}=1.0$) isothermal expansion
processes, whereas the purple and orange curves represent the (ii) isentropic
compression and (iv) isentropic expansion processes, respectively. The cycle
starts from the red arrow and evolves in a counterclockwise fashion over time
(heat engine). In (2), the $A_{1}$–$D_{1}$ diagrams, the red and orange dashed
curves represent (i′) removing and (iv′) attaching the hot bath, whereas in
(3), the $A_{2}$–$D_{2}$ diagrams, the purple and blue dashed curves represent
(ii′) attaching and (iii′) removing the cold bath. The cycles in the
$A_{1}$–$D_{1}$ and $A_{2}$–$D_{2}$ diagrams start from the red and blue
arrows and evolve in clockwise and counterclockwise fashions over time
(refrigerator and heat engine), respectively.

Under the quasi-static condition, we can employ a thermodynamic description of
the system because the variables defined from the qHE by Eqs. (21) and (22)
become the state variables.Sakamoto and Tanimura (2020, 2021); Koyanagi and
Tanimura (2022) By comparing these thermodynamic variables with the physical
quantities calculated by HEOM, we can quantitatively verify the qHE
description of quantum thermodynamics.

In Fig. 4, we depict the results for the IDF, hot ATF, and cold ATF as the (1)
$B$–$M$, (2) $A_{1}$–$D_{1}$, and (3) $A_{2}$–$D_{2}$ diagrams in the quasi-
static case for the (a) weak, (b) intermediate, and (c) strong SB coupling
cases. In these diagrams, the bath-removal processes (i′) and (iii′) and bath-
attachment processes (ii′) and (iv′) are represented as dashed curves. The
trajectories of the work diagram are periodic and closed because $M$, $D_{1}$,
and $D_{2}$ are the state variables. The area enclosed by each diagram
corresponds to negative work when evolving in a counterclockwise fashion over
time, and to positive work when evolving in a clockwise fashion. We have
numerically confirmed that $W_{\rm A+I}(\beta;t)$, evaluated as the areas
surrounded by the trajectories of the diagrams, agrees with $W_{\rm A+I}^{\rm
HEOM}(t)$ calculated from the HEOM. (See Table 2.)

In the case of weak and intermediate SB coupling in Figs. 4(2a), 4(2b), 4(3a),
and 4(3b), the positive work done by the hot ATF and the negative work done by
the cold ATF are approximately equivalent and almost cancel each other out.
Thus, the total work is predominantly determined by the work done by the IDF,
as illustrated in the $B$–$M$ diagram. In the cases of weak and intermediate
coupling, the profiles in Figs. 4(1a) and 4(1b) thus become similar to that in
the factorized assumption case in Fig. 2.

In the cases of strong coupling in Figs. 4(1c)–4(3c), because the SB
interaction includes contributions from both the subsystem side and the bath
side, the profiles in the $A$–$D$ and $B$–$M$ diagrams change significantly.
If this TLS is regarded as a spin system, $B(t)$ corresponds to the excitation
energy of the spin. Then, as $B(t)$ increases, the spin becomes aligned with
the ground state, so the magnetization $M(\beta;t)$ increases. Because the SB
interaction with $\hat{V}=\hat{\sigma}_{x}$ excites the spin, an increase or
decrease in $A_{k}(t)$ causes a decrease or increase in $M(\beta;t)$, even if
$B(t)$ remains constant. Thus, in the $B$–$M$ diagram in Fig. 4(1c), we
observe a contribution from the red and orange dashed lines resulting from
changes in $A_{1}(t)$ in processes (i′) and (iv′) and from the purple and blue
dashed lines resulting from changes in $A_{2}(t)$ in processes (ii′) and
(iii′). The areas surrounded by the trajectories of the $B$–$M$ diagrams are
divided vertically into two parts, with the processes evolving in a
counterclockwise fashion over time at the top and those evolving in a
clockwise fashion over time at the bottom. In this case, the work is evaluated
from the area of the top minus the area of the bottom. In Fig. 4(1c), such
effects appear as the vertical shifts of the solid curves corresponding to
(i)–(iv).

In the case of the $A$–$D$ diagrams, an increase in $B(t)$ suppresses the spin
excitation effect of $A_{k}(t)$, so $D_{k}(\beta;t)$ decreases even if
$A_{k}(t)$ does not change, as shown by the red and blue horizontal lines in
the $A_{1}$–$D_{1}$ and $A_{2}$–$D_{2}$ diagrams in Figs. 4(2c) and 4(3c),
respectively. The suppression effect of $B(t)$ is smaller for the hot bath
than for the cold bath, and the absolute value of the work done is larger for
$A_{1}$ than for $A_{2}$. The difference between the work done by $A_{1}$ and
that done by $A_{2}$ is compensated by the work done by $B$, as presented in
Table 2. As we will show below using a $T$–$S$ diagram, the bath removal and
bath attachment processes dominate the heat transfer processes in the strong
SB interaction case.

Although the profiles of the $B$–$M$ diagrams are very different when the SB
coupling strength is large, the total work done by the external fields is
almost the same, regardless of the SB coupling strength, in the quasi-static
limit, providing a manifestation of Carnot’s theorem. (See Fig. 3(ii) and
Table 2.)

Table 2: Work done in one cycle for three values of the SB coupling strength
$A_{0}$. Here, $W_{a}$ with $a=B$–$M$, $A_{1}$–$D_{1}$, and $A_{2}$–$D_{2}$
represents the work evaluated as the areas surrounded by the trajectories of
the work diagrams in Fig. 4, and $W_{b}^{\rm HEOM}$ with $b=A$, $I_{1}$, and
$I_{2}$ represents the work directly evaluated from the HEOM approach with
Eqs. (16) and (17) for one cycle of the external fields $B(t)$, $A_{1}(t)$,
and $A_{2}(t)$, respectively. In both cases, the value of the total work is
identical and, in the HEOM case, it is evaluated as $W_{\rm tot}^{\rm
HEOM}=W_{\rm A}^{\rm HEOM}+W_{\rm I_{1}}^{\rm HEOM}+W_{\rm I_{2}}^{\rm HEOM}$.
Then, as we illustrated in Table 3 in Appendix B, we have $W_{\rm tot}^{\rm
HEOM}+Q_{\rm tot}^{\rm HEOM}=0$
$A_{0}$ | $\quad W_{\rm B-M}\quad$ | $\quad W_{\rm A}^{\rm HEOM}\quad$ | $\qquad W_{\rm A_{1}-D_{1}}\qquad$ | $W_{\rm I_{1}}^{\rm HEOM}$ | $\qquad W_{\rm A_{2}-D_{2}}\qquad$ | $W_{\rm I_{2}}^{\rm HEOM}$ | $\qquad W_{\rm tot}^{\rm HEOM}\qquad$
---|---|---|---|---|---|---|---
$0.1$ | $-0.127$ | $-0.127$ | $4.83\times 10^{-4}$ | $4.83\times 10^{-4}$ | $-4.90\times 10^{-4}$ | $-4.90\times 10^{-4}$ | $-0.127$
$\sqrt{0.1}$ | $-0.128$ | $-0.128$ | $4.81\times 10^{-3}$ | $4.81\times 10^{-3}$ | $-4.85\times 10^{-3}$ | $-4.85\times 10^{-3}$ | $-0.129$
$\sqrt{3.0}$ | $-0.149$ | $-0.149$ | $0.122$ | $0.122$ | $-0.101$ | $-0.101$ | $-0.128$

####  III.2.4 T–S diagram

Figure 5:  $T$–$S$ work diagrams for the quasi-static process
($\Omega/2\pi=1.25\times 10^{-5}$) in the (a) weak, (b) intermediate, and (c)
strong SB coupling cases, corresponding to (a)–(c) in Fig. 4, respectively.
The red and blue lines represent the (i) hot isothermal compression and (iii)
cold isothermal expansion processes, whereas the purple and orange dashed
lines are the (ii) isentropic compression and (iv) isentropic expansion
processes, respectively. In each figure, we set $\Delta S_{\rm A+I}^{\rm
qst}(0)=0$ and the cycle starts from the hot isothermal compression, as
illustrated by the red arrow. In the strong SB coupling case in (c), the bath
removal and attachment processes, represented by the blue and red dashed
lines, determine the amount of heat production.

We now consider the $T$–$S$ diagram. Here, the change in the quasi-equilibrium
Boltzmann entropy (qBE) is calculated from Eq. (25) by numerically
differentiating $\Delta F^{\rm qst}_{\rm A+I}(\beta_{k};t)$ with respect to
$\beta_{k}$ for $k=1$ and 2 using the finite difference expression $[\Delta
F^{\rm qst}_{\rm A+I}(\beta_{k}+\Delta\beta;t)-\Delta F^{\rm qst}_{\rm
A+I}(\beta_{k};t)]/\Delta\beta$, with $\Delta\beta=0.0001$.

Although the $B$–$M$ and $A$–$D$ diagrams are useful for interpreting the work
that is directly evaluated using the HEOM approach, the situation for the
$T$–$S$ diagram is different because the entropy is purely a thermodynamic
variable that is indirectly evaluated using the temperature. Therefore, the
$T$–$S$ diagram is better suited to exploring the thermodynamic features of
the system. The area of the $T$–$S$ diagram is equivalent to the net heat
gain, which can be evaluated from the HEOM formalism as the change in the bath
energy. In Fig. 5, we depict the $T$–$S$ diagram for the (a) weak
[$A_{0}=0.1$], (b) intermediate [$A_{0}=\sqrt{0.1}$], and (c) strong
[$A_{0}=\sqrt{3.0}$] SB interaction cases. After the heat bath is removed, the
system is isolated and the entropy does not change under the isentropic
manipulation described in Sec. III.1. For reference, we also plot the results
based on the von Neumann entropy calculated from the zeroth HEOM element as
$S_{\rm A}^{\rm vN}(t)=-\mathrm{tr}_{\rm
A}\\{\hat{\rho}_{\vec{0}}(t)\ln\hat{\rho}_{\vec{0}}(t)\\}$ in Appendix B. The
trajectories of the $T$–$S$ diagram are periodic and closed because the qBE is
the state variable. The area enclosed by each diagram corresponds to the
positive heat when evolving in a counterclockwise fashion over time. (See
Table 3 in Appendix B.)

For the weak and intermediate SB coupling cases in Figs. 5(a) and 5(b), the
profiles of the $T$–$S$ diagram are similar to the factorized case in Fig.
1(b). This is because the contribution of the entropy from the SB interaction
is small in these cases. Conversely, for the strong SB coupling case in Fig.
5(c), the $T$–$S$ diagram differs significantly due to the entropy change
associated with the bath removal and bath attachment processes. This is
because the system can be excited by the SB interaction through
$\hat{V}=\hat{\sigma}_{x}$, so that entropy increases during the bath
attachment process as the heat flowing from the bath to the system increases,
whereas entropy decreases during the bath removal process due to reverse heat
flow. Thus, although the area enclosed by the lines corresponding to the net
heat gain is the same, as illustrated in Table 3, the heat flow is controlled
by the bath attachment and removal processes in the strong coupling case,
whereas the heat flow is controlled by the isothermal process in the weak and
intermediate coupling cases. However, the area enclosed by the trajectory in
the $T$–$S$ diagram does not depend on the SB coupling strength, as in the
case of work.

##  IV Conclusions

We have conducted accurate numerical simulations of a quantum Carnot engine
based on an SB model, incorporating the ATFs and IDF to describe the
manipulations of the adiabatic wall by external forces. The HEOM formalism
enables the evaluation of the internal energies of not only the subsystem, but
also the bath and the SB interactions, even in low-temperature, non-Markovian,
and nonperturbative conditions. We analyzed the real-time responses of work as
functions of the cycle frequency. As expected, the computational results
approach the Carnot efficiency in the quasi-static limit. Moreover, the
efficiency does not change, regardless of the SB coupling strength, which is
the manifestation of Carnot’s theorem.

If the ATFs are abruptly changed, as in the conventional theory when
introducing an adiabatic process, the cycle is no longer quasi-static and the
efficiency becomes lower than the Carnot limit. The use of reduced equations
of motion with Markovian and rotating wave approximations, such as the
Lindblad equation, should be avoided because this does not accurately treat
the quantum entanglement between the subsystem and the bath or the dynamic
thermal effects described by the fluctuation and dissipation arising from the
heat bath.Tanimura (2014, 2015)

By regarding work as the qHE, we computed the conjugate thermodynamic
variables of the IDF and ATFs and depicted the results as thermodynamic work
diagrams. We further introduced the qBE in an attempt to characterize the
thermal properties of the cycle. We showed that work variables such as
$M(\beta;t)$, $D_{1}(\beta;t)$, and $\Delta S_{{\rm A}+{\rm I}}^{\rm
qst}(\beta;t)$ are state variables. Thus, the areas surrounded by the
trajectories in the $M$–$B$ and $A$–$D$ diagrams correspond to work, whereas
the area within the trajectory of the $T$–$S$ diagram corresponds to the net
heat gain. The results show that ATFs are responsible for the work when the SB
coupling is very strong, whereas the IDF is responsible when the SB coupling
is not strong. As expected, the total work and net heat gain were found to be
independent of the SB coupling strength in the quasi-static limit, as
predicted by Carnot almost 200 years ago. The inclusion of ATFs is the key to
maintaining physical consistency in the study of quantum thermodynamics, and
perhaps for quantum information.

Although here we employed a simple spin-boson system, we can apply the same
method to ideal gas systems characterized by a $P$–$V$ diagram by using the
hierarchical quantum Fokker–Planck equations.Tanimura (2015); Tanimura and
Wolynes (1992); Ikeda and Tanimura (2019); Tanimura (2020) In the case of a
bosonic gas, the results should be similar to those for a classical ideal gas,
except at very low temperatures; conversely, in the case of a fermionic gas,
the results should exhibit quantum effects even at high temperatures. It would
be interesting to investigate heat engines with different mechanisms, although
it seems unlikely that the thermodynamic laws will be violated.

In this investigation, we limited our analysis to work variables in a quasi-
static case. As an extension of this study, it is possible to introduce
nonequilibrium free energy using the nonequilibrium work. Such investigations
are left for future work.

## Acknowledgments

Y.T. was supported by JSPS KAKENHI (Grant No. B21H01884). S.K. acknowledges a
fellowship supported by JST, the establishment of university fellowships
towards the creation of science technology innovation (Grant No. JPMJFS2123).

## 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 Isentropic processes

In the classical Carnot cycle, the isentropic process is regarded as a
temperature-changing process in which the internal energy also changes. In the
present case, we cannot define the temperature during the isentropic process
because the system is microscopic and isolated, aside from the initial and
final states described by $\beta_{1}$ or $\beta_{2}$. Thus, instead of
assuming a spontaneous change in internal energy, we actively control the
external field (IDF) to realize the temperature change between $\beta_{a}$ and
$\beta_{b}$ in the isentropic processes.

Because the heat baths have been removed, we only consider the main system
here. For an isentropic process, we have

$dU_{\rm A}=-M_{\rm A}^{0}dB,$ (32)

because the change in entropy is zero. Alternatively, the internal energy is
expressed as

$dU_{\rm A}=\left(\frac{\partial U_{\rm A}}{\partial\beta}\right)_{\rm
B}d\beta+\left(\frac{\partial U_{\rm A}}{\partial B}\right)_{\beta}dB.$ (33)

Note that, in statistical physics, the free energy defined from the partition
function is ubiquitously referred to as the Helmholtz energy. In the
thermodynamic sense, the Helmholtz energy for $B$ should be called the Gibbs
energy because they appear as the Legendre transformation of $F$ with respect
to $B$ and its conjugate property $M$ ($dG=-SdT-MdB$).Oono (2017) Accordingly,
$dU=TdS+MdB$ should be called enthalpy instead of internal energy. Following
convention, however, we refer to $F$ as Helmholtz energy and $U$ as internal
energy.

Because the internal energy of the TLS is evaluated from Eq. (27) as

$U_{\rm A}^{0}(t)=-\sqrt{B^{2}(t)+E^{2}}\tanh(\beta\sqrt{B^{2}(t)+E^{2}}),$
(34)

we have

$\left(\frac{\partial U_{\rm A}}{\partial\beta}\right)_{\rm
B}=-\frac{{B^{2}(t)+E^{2}}}{\cosh^{2}(\beta\sqrt{B^{2}(t)+E^{2}})},$ (35)
$\displaystyle\left(\frac{\partial U_{\rm A}}{\partial B}\right)_{\beta}=$
$\displaystyle-\frac{B(t)}{\sqrt{B^{2}(t)+E^{2}}}\tanh(\beta\sqrt{B^{2}(t)+E^{2}})$
$\displaystyle-$ $\displaystyle\frac{\beta
B(t)}{\cosh^{2}(\beta\sqrt{B^{2}(t)+E^{2}})}.$ (36)

From Eqs. (32) and (33), we have the relation

$\displaystyle-(B^{2}+E^{2})d\beta=\beta BdB.$ (37)

For isentropic processes from $B(t_{a})=B_{a}$ at $\beta_{a}$ to
$B(t_{b})=B_{b}$ at $\beta_{\rm B}$, we have

$-\ln\left(\frac{\beta_{b}}{\beta_{a}}\right)=\ln\left(\frac{\sqrt{E^{2}+B_{b}^{2}}}{\sqrt{E^{2}+B_{a}^{2}}}\right).$
(38)

Thus, we have the relation

$B_{b}=\sqrt{\frac{\beta_{a}}{\beta_{b}}\left(B_{a}^{2}+E^{2}\right)-E^{2}}.$
(39)

##  Appendix B Quasi-equilibrium Boltzmann entropy and von Neumann entropy

Figure 6:  Time evolution of the qBE (color curves) and vNE (black dashed
curves) in the (a) weak, (b) intermediate, and (c) strong SB coupling cases
for the eight-stroke cycle denoted by (i)–(iv’). (For the correspondence
between each stroke and the line colors, see Fig. 2.) Table 3:  Net heat gain,
evaluated as the area in the $T$–$S$ diagrams based on the qBE [$Q_{qBE}$] and
vNE [$Q_{\rm vNE}$] for different SB coupling strengths $A_{0}$. They are
compared with one cycle of the bath energy [$Q_{tot}^{\rm HEOM}\equiv\Delta
U_{{\rm B}_{1}}^{\rm qst}(\beta_{1};t)+\Delta U_{{\rm B}_{2}}^{\rm
qst}(\beta_{2};t)$], evaluated using Eq. (13) from the HEOM approach.
$A_{0}$ | $\quad Q_{qBE}\quad$ | $\quad Q_{\rm vNE}\quad$ | $Q_{\rm tot}^{\rm HEOM}$
---|---|---|---
$0.1$ | $0.1285$ | $0.1289$ | $0.1272$
$\sqrt{0.1}$ | $0.1286$ | $0.1289$ | $0.1287$
$\sqrt{3.0}$ | $0.1292$ | $0.1290$ | $0.1283$

The von Neumann entropy (vNE) is commonly used in quantum thermodynamics,
particularly in nonequilibrium cases. It is defined as follows:

$\displaystyle S_{A}^{\rm vN}(t)=-{\operatorname{tr}_{\rm
A}}\\{{\hat{\rho}_{\rm A}(t)}\ln{\hat{\rho}_{\rm A}(t)}\\},$ (40)

where ${\hat{\rho}_{A}(t)}$ is the reduced density matrix. In isolated
systems, the qBE and vNE are equivalent, but in a reduced system, this
equivalence does not hold due to the contribution of entropy from the SB
interaction. If we use the reduced density matrix obtained from
nonperturbative approaches, such as the HEOM approach, the description of vNE
becomes reasonably accurate, even in the strong SB coupling case.Sakamoto and
Tanimura (2020) Here, using the zeroth member of the HEOM solution in Eq. (40)
as $\hat{\rho}_{A}(t)=\hat{\rho}_{\vec{0}}$, we introduce the vNE to
characterize the time evolution of the qBE more closely. Because the Carnot
cycle involves two isentropic processes, where the qBE and vNE results must
agree, we can verify the description of vNE in the isotropic processes, as
well as the transitions between the isentropic and isothermal processes.

Using the calculated entropy, we also computed the net heat gain per cycle.
For example, for the vNE, $Q_{\rm vNE}$ is evaluated as

$\displaystyle Q_{\rm
vNE}=\int_{-\frac{\pi}{4\Omega}}^{\frac{\pi}{2\Omega}}dtT_{1}\frac{\partial
S_{\rm A}^{\rm vN}(\beta_{1};t)}{\partial
t}+\int_{\frac{3\pi}{4\Omega}}^{\frac{3\pi}{2\Omega}}dtT_{2}\frac{\partial
S_{\rm A}^{\rm vN}(\beta_{2};t)}{\partial t}.$ (41)

The calculated results are presented in Table 3. Although the net heat gain
calculated directly from HEOM agrees with the total work presented in Table 2,
there is a slight discrepancy between $Q_{\rm tot}^{\rm HEOM}$ and the other
values obtained from the qBE and vNE. This may be due to the insufficient
precision of the $\beta$ derivative in evaluating the entropy.

In Fig. 6, we depict the qBE and vNE as a function of time for different SB
coupling strengths. Here, we set the origin of the two graphs as $S_{\rm
A+I}^{\rm qst}(t)=S_{\rm A}^{\rm vN}(t)$ at $t=\pi/2\Omega$. As illustrated in
Figs. 6(a) and 6(b), the results from the qBE and vNE are similar and almost
overlap. This is because the contribution of the entropy from the SB
interaction is small in these cases. In the strong coupling case in panel (c),
while the time evolution of entropy in the isentropic processes is similar,
that in the isothermal processes controlled by $B(t)$ is different. This is
because the contribution of the entropy from the SB interaction is not
accurately taken into account in the vNE due to its reduced description of the
system. Thus, the vNE underestimates the entropy compared with the qBE in the
isothermal processes. Because their contribution cancels out, the net heat
gain, evaluated as the area in $T$–$S$ diagrams, is similar in both the qBE
and vNE cases, as presented in Table 3.

## References

  * Carnot (1824) S. Carnot,  _Reflexions sur la puissance motrice du feu_ (ChexzBachelier Libraire, Quai des Augustins, Paris, 1824).
  * Clapeyron (1834) B. P. E. Clapeyron, “Memoire sur la puissance motrice de la chaleur,” Journal de l’ Ecole Royale Polytechnique 14, 153–190 (1834).
  * Thomson (1848) W. Thomson, “On an absolute thermometric scale founded on Carnot’s theory of the motive power of heat,” Cambridge Philosophical Society Proceedings for June 5; and Phil. Mag., Oct. , 100–106 (1848).
  * Clausius (1850) R. Clausius, “Ueber die bewegende kraft der wärme und die gesetze, welche sich daraus für die wärmelehre selbst ableiten lassen,” Annalen der Physik 155, 500–524 (1850).
  * Clausius (1856) R. Clausius, “Ueber die anwendung der mechanischen wärmetheorie auf die dampfmaschine,” Annalen der Physik 173, 513–558 (1856).
  * Clausius (1865) R. Clausius, “Ueber verschiedene für die anwendung bequeme formen der hauptgleichungen der mechanischen wärmetheorie,” Annalen der Physik 201, 353–400 (1865).
  * Binder  _et al._ (2018) F. Binder, L. A. Correa, G. C, J. Anders, and G. A. (Eds.),  _Thermodynamics in the Quantum Regime_ (Springer International Publishing, 2018).
  * Kosloff (2013) R. Kosloff, “Quantum thermodynamics: A dynamical viewpoint,” Entropy 15, 2100–2128 (2013).
  * 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).
  * Härtle  _et al._ (2013) R. Härtle, G. Cohen, D. R. Reichman, and A. J. Millis, “Decoherence and lead-induced interdot coupling in nonequilibrium electron transport through interacting quantum dots: A hierarchical quantum master equation approach,” Phys. Rev. B 88, 235426 (2013).
  * Ludovico  _et al._ (2014) M. F. Ludovico, J. S. Lim, M. Moskalets, L. Arrachea, and D. Sánchez, “Dynamical energy transfer in ac-driven quantum systems,” Phys. Rev. B 89, 161306 (2014).
  * Esposito, Ochoa, and Galperin (2015a) M. Esposito, M. A. Ochoa, and M. Galperin, “Quantum thermodynamics: A nonequilibrium Green’s function approach,” Phys. Rev. Lett. 114, 080602 (2015a).
  * Cohen and Galperin (2020) G. Cohen and M. Galperin, “Green’s function methods for single molecule junctions,” The Journal of Chemical Physics 152, 090901 (2020), https://doi.org/10.1063/1.5145210 .
  * Brandão  _et al._ (2015) F. Brandão, M. Horodecki, N. Ng, J. Oppenheim, and S. Wehner, “The second laws of quantum thermodynamics,”  112, 3275–3279 (2015).
  * Haughian, Esposito, and Schmidt (2018) P. Haughian, M. Esposito, and T. L. Schmidt, “Quantum thermodynamics of the resonant-level model with driven system-bath coupling,” Phys. Rev. B 97, 085435 (2018).
  * Perarnau-Llobet  _et al._ (2018) M. Perarnau-Llobet, H. Wilming, A. Riera, R. Gallego, and J. Eisert, “Strong coupling corrections in quantum thermodynamics,” Phys. Rev. Lett. 120, 120602 (2018).
  * Funo and Quan (2018) K. Funo and H. T. Quan, “Path integral approach to heat in quantum thermodynamics,” Phys. Rev. E 98, 012113 (2018).
  * Rivas (2020) A. Rivas, “Strong coupling thermodynamics of open quantum systems,” Phys. Rev. Lett. 124, 160601 (2020).
  * 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).
  * Talkner and Hänggi (2020) P. Talkner and P. Hänggi, “Colloquium: Statistical mechanics and thermodynamics at strong coupling: Quantum and classical,” Rev. Mod. Phys. 92, 041002 (2020).
  * Zagoskin, Il’ichev, and Nori (2012) A. M. Zagoskin, E. Il’ichev, and F. Nori, “Heat cost of parametric generation of microwave squeezed states,” Phys. Rev. A 85, 063811 (2012).
  * Zagoskin  _et al._ (2012) A. M. Zagoskin, S. Savel’ev, F. Nori, and F. V. Kusmartsev, “Squeezing as the source of inefficiency in the quantum Otto cycle,” Phys. Rev. B 86, 014501 (2012).
  * Guthrie  _et al._ (2022) A. Guthrie, C. D. Satrya, Y.-C. Chang, P. Menczel, F. Nori, and J. P. Pekola, “Cooper-pair box coupled to two resonators: An architecture for a quantum refrigerator,” Phys. Rev. Applied 17, 064022 (2022).
  * Funo  _et al._ (2019) K. Funo, N. Lambert, B. Karimi, J. P. Pekola, Y. Masuyama, and F. Nori, “Speeding up a quantum refrigerator via counterdiabatic driving,” Phys. Rev. B 100, 035407 (2019).
  * Nomura, Yoshioka, and Nori (2021) Y. Nomura, N. Yoshioka, and F. Nori, “Purifying deep Boltzmann machines for thermal quantum states,” Phys. Rev. Lett. 127, 060601 (2021).
  * Funo, Lambert, and Nori (2021) K. Funo, N. Lambert, and F. Nori, “General bound on the performance of counter-diabatic driving acting on dissipative spin systems,” Phys. Rev. Lett. 127, 150401 (2021).
  * Lu  _et al._ (2022) Y. Lu, N. Lambert, A. F. Kockum, K. Funo, A. Bengtsson, S. Gasparinetti, F. Nori, and P. Delsing, “Steady-state heat transport and work with a single artificial atom coupled to a waveguide: Emission without external driving,” PRX Quantum 3, 020305 (2022).
  * An  _et al._ (2015) S. An, J.-N. Zhang, M. Um, D. Lv, Y. Lu, J. Zhang, Z.-Q. Yin, H. T. Quan, and K. Kim, “Experimental test of the quantum Jarzynski equality with a trapped-ion system,” Nature Physics 11, 193 (2015).
  * An  _et al._ (2016) S. An, D. Lv, del Campo Adolfo, and K. Kim, “Shortcuts to adiabaticity by counterdiabatic driving for trapped-ion displacement in phase space,” Nature Communications 7, 12999 (2016).
  * I. A  _et al._ (2016) M. I. A, E. Roldan, L. Dinis, D. Petrov, J. M. R. Parrondo, and R. A. Rica, “Brownian Carnot engine,” Nature Phys 12, 67 (2016).
  * Vitanov  _et al._ (2017) N. V. Vitanov, A. A. Rangelov, B. W. Shore, and K. Bergmann, “Stimulated Raman adiabatic passage in physics, chemistry, and beyond,” Rev. Mod. Phys. 89, 015006 (2017).
  * Hoffmann  _et al._ (2018) E. A. Hoffmann, S. Fahlvik, C. Thelander, M. Leijnse, and H. Linke, “A quantum-dot heat engine operating close to the thermodynamic efficiency limits,” Nature Nanotechnology 13, 920–924 (2018).
  * Bengtsson  _et al._ (2018) J. Bengtsson, M. N. Tengstrand, A. Wacker, P. Samuelsson, M. Ueda, H. Linke, and S. M. Reimann, “Quantum Szilard engine with attractively interacting Bosons,” Phys. Rev. Lett. 120, 100601 (2018).
  * Josefsson  _et al._ (2019) M. Josefsson, A. Svilans, H. Linke, and M. Leijnse, “Optimal power and efficiency of single quantum dot heat engines: Theory and experiment,” Phys. Rev. B 99, 235432 (2019).
  * Prete  _et al._ (2019) D. Prete, P. A. Erdman, V. Demontis, V. Zannier, D. Ercolani, L. Sorba, F. Beltram, F. Rossella, F. Taddei, and S. Roddaro, “Thermoelectric conversion at 30 K in InAs/InP nanowire quantum dots,” Nano Letters 19, 3033–3039 (2019), pMID: 30935206.
  * Ma  _et al._ (2020) Y.-H. Ma, R.-X. Zhai, J. Chen, C. P. Sun, and H. Dong, “Experimental test of the $1/\tau$-scaling entropy generation in finite-time thermodynamics,” Phys. Rev. Lett. 125, 210601 (2020).
  * de Oliveira  _et al._ (2020) A. de Oliveira, R. Gomes, V. Brasil, N. Rubiano da Silva, L. Céleri, and P. Souto Ribeiro, “Full thermalization of a photonic qubit,” Physics Letters A 384, 126933 (2020).
  * Hernández-Gómez  _et al._ (2021) S. Hernández-Gómez, N. Staudenmaier, M. Campisi, and N. Fabbri, “Experimental test of fluctuation relations for driven open quantum systems with an NV center,” New Journal of Physics 23, 065004 (2021).
  * Ryu  _et al._ (2022) S. Ryu, R. López, L. Serra, and D. Sánchez, “Beating Carnot efficiency with periodically driven chiral conductorse,” Nature Communications 13, 2512 (2022).
  * 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 (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 (2015) Y. Tanimura, “Real-time and imaginary-time quantum hierarchal Fokker-Planck equations,” The Journal of Chemical Physics 142, 144110 (2015).
  * Segal, Nitzan, and Hänggi (2003) D. Segal, A. Nitzan, and P. Hänggi, “Thermal conductance through molecular wires,” The Journal of Chemical Physics 119, 6840–6855 (2003), https://doi.org/10.1063/1.1603211 .
  * Segal and Nitzan (2005) D. Segal and A. Nitzan, “Spin-Boson thermal rectifier,” Phys. Rev. Lett. 94, 034301 (2005).
  * Segal (2006) D. Segal, “Heat flow in nonlinear molecular junctions: Master equation analysis,” Phys. Rev. B 73, 205415 (2006).
  * Boudjada and Segal (2014) N. Boudjada and D. Segal, “From dissipative dynamics to studies of heat transfer at the nanoscale: Analysis of the spin-Boson model,” The Journal of Physical Chemistry A 118, 11323–11336 (2014), pMID: 25396751.
  * Dhar (2008) A. Dhar, “Heat transport in low-dimensional systems,” Advances in Physics 57, 457–537 (2008).
  * 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).
  * Dubi and Di Ventra (2011) Y. Dubi and M. Di Ventra, “Colloquium: Heat flow and thermoelectricity in atomic and molecular junctions,” Rev. Mod. Phys. 83, 131–155 (2011).
  * Wang, Ren, and Cao (2015) C. Wang, J. Ren, and J. Cao, “Nonequilibrium energy transfer at nanoscale: A unified theory from weak to strong coupling,” Scientific Reports 5 (2015), 10.1038/srep11787.
  * Saito and Kato (2013) K. Saito and T. Kato, “Kondo signature in heat transfer via a local two-state system,” Phys. Rev. Lett. 111, 214301 (2013).
  * 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).
  * Aurell and Montana (2019) E. Aurell and F. Montana, “Thermal power of heat flow through a qubit,” Phys. Rev. E 99, 042130 (2019).
  * Kilgour, Agarwalla, and Segal (2019) M. Kilgour, B. K. Agarwalla, and D. Segal, “Path-integral methodology and simulations of quantum thermal transport: Full counting statistics approach,” The Journal of Chemical Physics 150, 084111 (2019).
  * Xu, Stockburger, and Ankerhold (2021) M. Xu, J. T. Stockburger, and J. Ankerhold, “Heat transport through a superconducting artificial atom,” Phys. Rev. B 103, 104304 (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).
  * Esposito, Ochoa, and Galperin (2015b) M. Esposito, M. A. Ochoa, and M. Galperin, “Nature of heat in strongly coupled open quantum systems,” Phys. Rev. B 92, 235440 (2015b).
  * 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).
  * Brunner  _et al._ (2014) N. Brunner, M. Huber, N. Linden, S. Popescu, R. Silva, and P. Skrzypczyk, “Entanglement enhances cooling in microscopic quantum refrigerators,” Phys. Rev. E 89, 032115 (2014).
  * Chotorlishvili, Toklikishvili, and Berakdar (2011) L. Chotorlishvili, Z. Toklikishvili, and J. Berakdar, “Thermal entanglement and efficiency of the quantum Otto cycle for the su(1,1) Tavis–Cummings system,”  44, 165303 (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.
  * 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).
  * 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).
  * Xu, Chen, and Liu (2018) Y. Y. Xu, B. Chen, and J. Liu, “Achieving the classical Carnot efficiency in a strongly coupled quantum heat engine,” Phys. Rev. E 97, 022130 (2018).
  * Restrepo  _et al._ (2018) S. Restrepo, J. Cerrillo, P. Strasberg, and G. Schaller, “From quantum heat engines to laser cooling: Floquet theory beyond the Born–Markov approximation,”  20, 053063 (2018).
  * Kilgour and Segal (2018) M. Kilgour and D. Segal, “Coherence and decoherence in quantum absorption refrigerators,” Phys. Rev. E 98, 012117 (2018).
  * Das and Mukherjee (2020) A. Das and V. Mukherjee, “Quantum-enhanced finite-time Otto cycle,” Physical Review Research 2 (2020), 10.1103/physrevresearch.2.033083.
  * Camati, Santos, and Serra (2020) P. A. Camati, J. F. G. Santos, and R. M. Serra, “Employing non-Markovian effects to improve the performance of a quantum Otto refrigerator,” Phys. Rev. A 102, 012217 (2020).
  * Johal and Mehta (2021) R. S. Johal and V. Mehta, “Quantum heat engines with complex working media, complete Otto cycles and heuristics,” Entropy 23 (2021), 10.3390/e23091149.
  * 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).
  * Wiedmann, Stockburger, and Ankerhold (2020) M. Wiedmann, J. T. Stockburger, and J. Ankerhold, “Non-Markovian dynamics of a quantum heat engine: out-of-equilibrium operation and thermal coupling control,” New Journal of Physics 22, 033007 (2020).
  * 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).
  * Xu  _et al._ (2022) M. Xu, J. T. Stockburger, G. Kurizki, and J. Ankerhold, “Minimal quantum thermal machine in a bandgap environment: non-Markovian features and anti-Zeno advantage,” New Journal of Physics 24, 035003 (2022).
  * Brenes  _et al._ (2020) M. Brenes, J. J. Mendoza-Arenas, A. Purkayastha, M. T. Mitchison, S. R. Clark, and J. Goold, “Tensor-network method to simulate strongly interacting quantum thermal machines,” Phys. Rev. X 10, 031040 (2020).
  * Pancotti  _et al._ (2020) N. Pancotti, M. Scandi, M. T. Mitchison, and M. Perarnau-Llobet, “Speed-ups to isothermality: Enhanced quantum thermal machines through control of the system-bath coupling,” Phys. Rev. X 10, 031015 (2020).
  * Lee  _et al._ (2020) S. Lee, M. Ha, J.-M. Park, and H. Jeong, “Finite-time quantum Otto engine: Surpassing the quasistatic efficiency due to friction,” Phys. Rev. E 101, 022127 (2020).
  * Raja  _et al._ (2021) S. H. Raja, S. Maniscalco, G. S. Paraoanu, J. P. Pekola, and N. L. Gullo, “Finite-time quantum Stirling heat engine,” New Journal of Physics 23, 033034 (2021).
  * Ivander, Anto-Sztrikacs, and Segal (2022) F. Ivander, N. Anto-Sztrikacs, and D. Segal, “Strong system-bath coupling effects in quantum absorption refrigerators,” Phys. Rev. E 105, 034112 (2022).
  * Goyal, He, and Kawai (2020) K. Goyal, X. He, and R. Kawai, “Entropy production of a small quantum system under strong coupling with an environment: A computational experiment,” Physica A: Statistical Mechanics and its Applications 552, 122627 (2020), tributes of Non-equilibrium Statistical Physics.
  * Seshadri and Galperin (2021) N. Seshadri and M. Galperin, “Entropy and information flow in quantum systems strongly coupled to baths,” Phys. Rev. B 103, 085415 (2021).
  * 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).
  * 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).
  * Goold  _et al._ (2016) J. Goold, M. Huber, A. Riera, L. del Rio, and P. Skrzypczyk, “The role of quantum information in thermodynamics—a topical review,”  49, 143001 (2016).
  * Strasberg  _et al._ (2017) P. Strasberg, G. Schaller, T. Brandes, and M. Esposito, “Quantum and information thermodynamics: A unifying framework based on repeated interactions,” Phys. Rev. X 7, 021003 (2017).
  * Jarzynski (2004) C. Jarzynski, “Nonequilibrium work theorem for a system strongly coupled to a thermal environment,” Journal of Statistical Mechanics: Theory and Experiment 2004, P09005 (2004).
  * Yukawa (2000) S. Yukawa, “A quantum analogue of the Jarzynski equality,” Journal of the Physical Society of Japan 69, 2367–2370 (2000).
  * 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).
  * Campisi, Talkner, and Hänggi (2009) M. Campisi, P. Talkner, and P. Hänggi, “Fluctuation theorem for arbitrary open quantum systems,” Phys. Rev. Lett. 102, 210401 (2009).
  * 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).
  * Seifert (2012) U. Seifert, “Stochastic thermodynamics, fluctuation theorems and molecular machines,”  75, 126001 (2012).
  * 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).
  * Aurell, Donvil, and Mallick (2020) E. Aurell, B. Donvil, and K. Mallick, “Large deviations and fluctuation theorem for the quantum heat current in the spin-boson model,” Phys. Rev. E 101, 052116 (2020).
  * Whitney (2018) R. S. Whitney, “Non-Markovian quantum thermodynamics: Laws and fluctuation theorems,” Phys. Rev. B 98, 085415 (2018).
  * 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).
  * 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).
  * Leff and Rex (2002) H. Leff and A. F. Rex, “Maxwell’s demon 2 entropy, classical and quantum information, computing,” (2002).
  * 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).
  * 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).
  * Koyanagi and Tanimura (2022) 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 (2022), https://doi.org/10.1063/5.0093666 .
  * 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).
  * Lenard (1978) A. Lenard, “Thermodynamical proof of the Gibbs formula for elementary quantum systems,” Journal of Statistical Physics 19, 575 (1978).
  * Tasaki (2000) H. Tasaki, “Jarzynski relations for quantum systems and some applications,” (2000), arXiv:cond-mat/0009244 [cond-mat.stat-mech] .
  * Oono (2017) Y. Oono,  _Perspectives on Statistical Thermodynamics_ (Cambridge University Press, 2017).
  * Zhang and Tanimura (2022) J. Zhang and Y. Tanimura, “Imaginary-time hierarchical equations of motion for thermodynamic variables,” The Journal of Chemical Physics 156, 174112 (2022).
  * Ishizaki and Fleming (2009) A. Ishizaki and G. R. Fleming, “On the adequacy of the Redfield equation and related approaches to the study of quantum dynamics in electronic energy transfer,” The Journal of Chemical Physics 130, 234110 (2009).
  * 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).