# Cooling and work extraction under memory-assisted Markovian thermal
processes

Yuancheng Liu  Xueyuan Hu<EMAIL_ADDRESS>School of Information Science and
Engineering, Shandong University, Qingdao 266237, China

###### Abstract

We investigate the limits on cooling and work extraction via Markovian thermal
processes assisted by a finite-dimensional memory. Here the memory is a
$d$-dimensional quantum system with trivial Hamiltonian and initially in a
maximally mixed state. For cooling a qubit system, we consider two paradigms,
cooling under coherent control and cooling under incoherent control. For both
paradigms, we derive the optimal ground-state populations under the set of
general thermal processes (TP) and the set of Markovian thermal processes
(MTP), and we further propose memory-assisted protocols, which bridge the gap
between the performances of TP and MTP. For the task of work extraction, we
prove that when the target system is a qubit in the excited state, the minimum
extraction error achieved by TP can be approximated by Markovian thermal
processes assisted by a large enough memory. Our results can bridge the
performances of TP and MTP in thermodynamic tasks including cooling and work
extraction.

##  I Introduction

Laws of thermodynamics pose essential restrictions on processes such as
cooling and work extraction. Ever since Szilard’s discussion of Maxwell’s
demon, where the amount of extracted work is related to the knowledge of the
target system, the information-theoretic approach has been employed to
investigate thermodynamic processes [1]. In particular, proposals for quantum
engines have been proposed [2, 3, 4, 5] and experimentally demonstrated [6, 7,
8, 9]. Further, quantum gates are employed to design schemes for cooling,
where the entropy of the target system is transferred to an auxiliary system,
which can then release the entropy to a heat bath. Such scheme of cooling is
called heat bath algorithmic cooling [10, 11, 12, 13].

Recently, great progresses have been made in resource-theoretic approaches to
quantum thermodynamics [14, 15]. In particular, conditions on state
transformations under the set of thermal processes (TP) are established [16].
Based on these conditions, bounds on work extraction are derived [16]. Also,
limitations on cooling in different regimes are derived [17, 18, 19, 20].
Further, the heat-bath algorithmic cooling is extended to involve general
thermal processes, such that the ground-state population of the target system
goes to 1 exponentially fast in the number of rounds [21].

When the set of thermal processes is restricted by the Markovian condition,
which means that the evolution of states under such processes can be described
by a master equation with the thermal state a fixed point, state
transformations are restricted by stronger limitations. Such processes are
defined as the set of Markovian thermal processes (MTP) in Ref. [22, 23],
where it is proved that any operation in MTP can be realized by a sequence of
elementary thermalizations. A direct consequence of this result is that under
MTP and without any auxiliary system, a two-level system cannot be cooled to a
temperature lower than that of the reservoir.

By noticing that memory effect of the environment sets the difference between
Markovian and non-Markovian processes, we introduce a $d$-dimensional
auxiliary system, with trivial Hamiltonian and initially in a maximally mixed
state, to assist state transformations under MTP. Because this auxiliary
system does not provide energy or non-equilibrium resources, it acts as a
“memory” in state transformations.

In this paper, we study the limits on cooling and work extraction via memory-
assisted Markovian thermal processes, and compare the performance with TP and
MTP. For cooling a qubit system, we consider two paradigms, cooling under
coherent control and cooling under incoherent control. For both paradigms, we
derive the optimal ground-state populations under TP and MTP, and find that
for the paradigm with coherence control, there is a gap between asymptotic
ground-state populations (in the limit of infinite rounds) achieved by TP and
MTP, while for the paradigm with incoherent control, both TP and MTP can reach
the same asymptotic ground-state population but the convergence rate differs.
We further propose memory-assisted protocols for cooling under coherent and
incoherent control, which bridges the gap between the performances of TP and
MTP. For the problem of extracting work from an out-of-equilibrium qubit
state, we prove that if the qubit is initially in the excited state, the
optimal strategy by TP can be approximated by MTP assisted by a large enough
memory, indicating that memory-assisted MTP can outperform elementary thermal
process when the dimension of the target is larger than 2. Our results may
shed light on the study of memory effect of non-Markovian thermodynamic
processes.

##  II Prelimitaries

###  II.1 Markovian thermal processes and related concepts

Consider a quantum system $S$ with Hamiltonian $H_{S}$ and surrounded by a
reservoir $R$ at inverse temperature $\beta$. In equilibrium, the state of $S$
reads $\tau\equiv e^{-\beta H_{S}}/\mathrm{Tr}(e^{-\beta H_{S}})$, which is
also called a Gibbs state.

Here we briefly review three sets of completely positive and trace-preserving
(CPTP) maps: the TP [24, 25],the MTP [22, 23], and ETP [26].

A CPTP map $\Lambda$ belongs to TP if and only if it satisfies the following
two conditions [24, 25].
(P1) $\Lambda$ is time-translationally symmetric

$\Lambda(e^{-iH_{S}t}\rho
e^{iH_{S}t})=e^{-iH_{S}t}\Lambda(\rho)e^{iH_{S}t},\forall\rho,t.$ (1)

(P2) $\Lambda$ preserves the Gibbs state

$\Lambda(\tau)=\tau.$ (2)

Here and following, we focus on quantum states which are diagonal on the
eigenbasis of Hamiltonian. Such states can be fully characterized as vectors
$\bm{p}$ of occupation probabilities $p_{k}=\langle k|\rho|k\rangle$, where
$|k\rangle$ are eigenstates of $H_{S}$ corresponding to energy levels $E_{k}$.
Then the action of a CPTP map on $\rho$ can be described as a stochastic
matrix acting on $\bm{p}$

$\rho^{\prime}=\Lambda(\rho)\Rightarrow\bm{p^{\prime}}=G\bm{p},$ (3)

Here $G$ is a matrix of transition probabilities
$G_{k^{\prime}k}=p_{k^{\prime}|k}\equiv\langle
k^{\prime}|\Lambda(|k\rangle\langle k|)|k^{\prime}\rangle$ from state
$|k\rangle$ to state $|k^{\prime}\rangle$, and $\bm{p^{\prime}}$ is a vector
of occupation probabilities $p^{\prime}_{k}=\langle k|\rho^{\prime}|k\rangle$
for the output state $\rho^{\prime}$. From (P2), the population dynamics $G$
induced by TP is a stochastic matrix which preserves the Gibbs distribution.
According to Ref. [16, 27],

$\bm{p}\stackrel{{\scriptstyle\mathrm{TP}}}{{\longrightarrow}}\bm{p^{\prime}}\Leftrightarrow\bm{p}\succ_{{}_{T}}\bm{p^{\prime}},$
(4)

where $\succ_{{}_{T}}$ stands for thermo-majorization.

An elementary thermal process [26] is realized as a sequence of thermal
processes which involve only two energy levels at a time. The set of ETP is
equivalent to TP for the qubit case, while for high dimensional systems, ETP
are strict subsets of TP [26]. When a thermal process involves only two energy
levels $|i\rangle$ and $|j\rangle$ with energy gap $E_{ij}$, the reduced
transition matrix in the subspace spanned by $\\{|i\rangle,|j\rangle\\}$ can
be written as

$G^{ij}=(1-\lambda)\mathbb{1}^{ij}+\lambda\beta^{ij},$ (5)

where $\lambda\in[0,1]$, $\mathbb{1}^{ij}$ is the identity matrix, and
$\beta^{ij}$ is called $\beta$-swap between $i$ and $j$, expressed as

$\beta^{ij}=\begin{pmatrix}1-e^{-\beta E_{ij}}&1\\\ e^{-\beta
E_{ij}}&0\end{pmatrix}.$ (6)

The following lemma follows directly from Eq. (5).

######  Lemma 1.

Consider a qubit system with Hamiltonian $H_{S}=E|e\rangle\langle e|$, and two
states which are incoherent on energy eigenbasis with occupation vectors
$\bm{p}=[p,1-p]$ and $\bm{p^{\prime}}=[p^{\prime},1-p^{\prime}]$. Then
$\bm{p}$ can be transformed to $\bm{p^{\prime}}$ via TP, if and only if

$\left\\{\begin{matrix}p\leq p^{\prime}\leq p_{\beta},&p\leq\gamma,\\\
p{{}_{\beta}}\leq p^{\prime}\leq p,&p\geq\gamma,\end{matrix}\right.$ (7)

where $p_{\beta}=[\beta^{g,e}\bm{p}]_{g}=1-pe^{-\beta E}$ is the ground-state
population of the output state obtained from $\beta$-swap, and $\gamma\equiv
1/(1+e^{-\beta E})$ is the ground-state population of the Gibbs state.

A Markovian thermal process is defined as a thermal process generated by a
Markovian master equation. A Markovian process is memoryless, meaning that the
evolution rate of a system at time $t$ solely depends on its state at time
$t$. As recently proved in Ref. [22], any state transformation induced by the
set of MTP can be achieved by a sequence of elementary partial
thermalizations, which involves two energy levels at a time. It follows that
MTP are subsets of ETP. A partial thermalization between energy levels $E_{i}$
and $E_{j}$ ($E\equiv E_{j}-E_{i}>0$) is expressed as

$T_{\lambda}^{ij}=\begin{pmatrix}1-\lambda(1-\gamma)&\lambda\gamma\\\
\lambda(1-\gamma)&1-\lambda\gamma\end{pmatrix},$ (8)

where $\lambda\in[0,1]$. For $\lambda=1$, $T_{\lambda}^{ij}$ transforms any
state to a balanced distribution on $i$ and $j$. It means that if
$\bm{p^{\prime}}=T_{1}^{ij}\bm{p}$, then
$p_{j}^{\prime}=p_{i}^{\prime}e^{-\beta E}$, $\forall\bm{p}$. This is called a
full thermalization between $i$ and $j$, which we will label $T^{ij}$ in the
following.

When a catalytic system is employed, the state transformations are enhanced
[28, 29, 30, 31, 32, 33, 34, 35, 36]. A catalyst is finite-dimensional
ancilla, which interacts with the system via free operations and then returns
to the exact original state. When the set of free operations is TP, the
initial state of a catalyst is non-equilibrium. Interestingly, as recently
discovered in Ref. [36], with Gibbs state catalysts, elementary thermal
processes can emulate any operation in TP.

###  II.2 Memory-assisted Markovian thermal process

Consider a qubit system with Hamiltonian $H_{S}=E|e\rangle\langle e|$ and
initially in state $\bm{p}=[p_{0},1-p_{0}]$. According to the results in Refs.
[22, 23], $\rho$ can be transformed to state
$\bm{p}^{\mathrm{MTP}}=[p^{\mathrm{MTP}},1-p^{\mathrm{MTP}}]$ by MTP, if and
only if, $\gamma\leq p^{\mathrm{MTP}}\leq p_{0}$ for $p_{0}>\gamma$, and
$p_{0}\leq p^{\mathrm{MTP}}\leq\gamma$ for $p_{0}\leq\gamma$. In other words,
any state which can be reached by MTP from a qubit state $\rho$, is a mixture
of $\rho$ and $\gamma$.

By employing a qubit ancillary with trivial Hamiltonian $H_{A}=0$ and
initially in a maximally mixed state
$\tau^{(2)}_{A}=\frac{1}{2}(|1\rangle\langle 1|+|2\rangle\langle 2|)$, a
global MTP can transform the state of $S$ to
$\rho^{(2)}=p^{(2)}|g\rangle\langle g|+(1-p^{(2)})|e\rangle\langle e|$ with

$p^{(2)}=\gamma+\gamma(1-\gamma)(\gamma-p_{0}).$ (9)

See Appendix B for detailed calculations. Because $p^{(2)}>\gamma$ for
$p_{0}<\gamma$ and $p^{(2)}<\gamma$ for $p_{0}>\gamma$, $\rho^{(2)}$ cannot be
reached from $\rho$ via MTP.

It is worth mentioning that the auxiliary system considered here is similar to
the catalyst system proposed in Ref. [22, 23], in the sense that the initial
and final states of the auxiliary system are the same. Nevertheless, here we
impose stronger requirements on the auxiliary system, i.e., trivial
Hamiltonian and maximally mixed initial state. In other words, it does not
provide energy or non-equilibrium resource for the state transformation of the
target system. The reason why it can enable state transformations which cannot
be realized by MTP is that, instead of fully thermalizing after each
elementary thermalization, the assisted system memorizes its state as well as
its correlation to the target system. Therefore, we will call state
transformations which are assisted by such auxiliary systems the memory-
assisted Markovian thermal processes.

######  Definition 1.

(Memory-assisted $\mathrm{MTP}$, $\mathrm{MMTP}^{(d)}$) A quantum channel
$\mathcal{E}$ is a memory-assisted Markovian thermal process if there are an
assisted system $A$ with trivial Hamiltonian $H_{A}=0$ and initially in a
maximally mixed state $\tau_{A}=\mathbb{1}_{A}/d$ and a joint Markovian
thermal process $\mathcal{E}^{\mathrm{MTP}}$ such that

$\mathcal{E}(\cdot)=\mathrm{Tr}_{A}[\mathcal{E}^{\mathrm{MTP}}(\cdot\otimes\tau_{A})].$
(10)

When the dimension of $A$ is restricted to $d$, the set of memory-assisted
Markovian thermal processes is labeled as $\mathrm{MMTP}^{(d)}$.

The above definition is similar to that proposed in the recent paper [37], but
we impose the trivial Hamiltonian condition to avoid energy exchange between
the memory and the other systems. The following theorem is also proved in Ref.
[37]. For consistency in this paper, we provide our proof in Appendix B.

######  Theorem 1.

Let $\bm{p}$ be an initial state of the target system, which undergoes
$\Lambda^{ij}\in\mathrm{MMTP}^{(d)}$ involving only two energy levels $E_{i}$
and $E_{j}$ ($E=E_{j}-E_{i}>0$). Then the state $\bm{p^{(d)}}$ with
$p_{i}^{(d)}$ in the following form can be reached,

$p_{i}^{(d)}=(1-e^{-\beta E})p_{i}+p_{j}+[(1-\gamma)p_{i}-\gamma
p_{j}]\delta_{d}(\gamma),$ (11)

where $\delta_{d}(\gamma)<o[(4\gamma(1-\gamma))^{d}d^{-3/2}]$.

Notice that the action of $\beta$-swap $\beta^{ij}$ leads to
$[\beta^{i,j}\bm{p}]_{i}=(1-e^{-\beta E})p_{i}+p_{j}$. Therefore, $\beta$-swap
can be simulated by $\mathrm{MMTP}^{(d)}$ with error exponentially decreasing
in the dimension of the assisted memory.

##  III Cooling

###  III.1 Settings and main results

Our target system is a qubit system with Hamiltonian $H_{S}=E|e\rangle\langle
e|$ and is surrounded by a reservoir $R$ at inverse temperature $\beta$.
Initially, the target system is thermalized to an equilibrium state
$\tau=\gamma|g\rangle\langle g|+(1-\gamma)|e\rangle\langle e|$ with
$\gamma=1/(1+e^{-\beta E})$.

Figure 1: Cooling Paradigms. (a) Cooling under coherence control. The target
qubit is first initialized as a thermal state in $R$. In step 1 of each round,
the target qubit transforms under the unitary $\sigma_{x}$; in step 2, it
undergoes a process $\Lambda^{\mathrm{F}}\in\mathrm{F}$, where
$\mathrm{F}=\mathrm{TP},\ \mathrm{MTP},$ or $\mathrm{MMTP}^{(d)}$. (b) Cooling
under incoherent control. To initialize the composed system, the target qubit
is thermalized in $R$ while the auxiliary is thermalized in $H$. In step 1 of
each round, one preserves the state of $S$ and thermalizes $A$ in $H$; in step
2 , the composed system of $S$ and $A$ undergoes a joint process
$\Lambda^{\mathrm{F}}\in\mathrm{F}$.

In order to cool the target system, we consider two paradigms: cooling under
coherent control and cooling under incoherent control, which are schematically
depicted in Fig. 1. Both paradigms consist of $n$ rounds. In each round, the
target system absorbs the energy for cooling and then undergoes a process
which belongs to TP, MTP, or $\mathrm{MMTP}^{(d)}$. The difference between the
two paradigms lies in the source of energy. In the paradigm of coherent
control, the energy for cooling comes from a unitary operation, while the
energy in incoherent control is provided by a hot bath $H$ at inverse
temperature $\beta_{H}<\beta$.

In the paradigm of coherent control, each round consists of two steps. In the
first step, a unitary operation is applied to invert the populations of the
target system, while the target system is decoupled from the heat bath. The
second step is to maximize the ground-state population via TP, MTP or
$\mathrm{MMTP}^{(d)}$. We denote the state of the target qubit after the $n$th
cycle as $\rho_{n}^{\mathrm{F},coh}=p_{n}^{\mathrm{F},coh}|g\rangle\langle
g|+(1-p_{n}^{\mathrm{F},coh})|e\rangle\langle e|$, where
$\mathrm{F}=\mathrm{TP},\ \mathrm{MTP}$, or $\mathrm{MMTP}^{(d)}$ denotes the
allowed operations in the cooling process, and the superscript $coh$ means
that the energy comes from coherent control. Our main result for coherent
control is summarized in the following theorem.

######  Theorem 2.

Under coherent control, the ground-state populations after the $n$th round,
under cooling processes involving $\mathrm{TP}$ and $\mathrm{MTP}$, are
respectively upper bounded by

$\displaystyle p_{n,*}^{\mathrm{TP},coh}$ $\displaystyle=$ $\displaystyle
1-(1-\gamma)e^{-n\beta E},$ (12) $\displaystyle p_{n,*}^{\mathrm{MTP},coh}$
$\displaystyle=$ $\displaystyle\gamma,$ (13)

where the bounds are reachable. Moreover, under $\mathrm{MMTP}^{(d)}$, the
ground-state population can reach

$p_{n,*}^{(d),coh}=p_{\max}^{(d),coh}-(e^{-\beta
E}-\delta_{d}(\gamma))^{n}(p_{\max}^{(d),coh}-\gamma),$ (14)

where

$p_{\max}^{(d),coh}=1-\frac{\gamma}{1+(1-e^{-\beta E})/\delta_{d}(\gamma)}.$
(15)

Detailed proof and discussion of the above theorem are in Sec. III.2.

In the paradigm of incoherent control, we employ an auxiliary qubit with
Hamiltonian $H_{A}=(\mathcal{E}-E)|1_{A}\rangle\langle 1_{A}|$, where we set
$\mathcal{E}>E$. Initially, the target system $S$ is thermalized in $R$ while
the auxiliary is thermalized in $H$. Each round of the protocol consists of
two steps. In the first step, the auxiliary is fully thermalized in $H$,
absorbing the energy for cooling, while the state of target system is
preserved. In the second step, the target system and the auxiliary are brought
in contact with $R$ and undergo a thermal process. After the $n$th round, the
state of the target system becomes
$\rho_{n}^{\mathrm{F},inc}=p_{n}^{\mathrm{F},inc}|g\rangle\langle
g|+(1-p_{n}^{\mathrm{F},inc})|e\rangle\langle e|$, where the superscript $inc$
indicates incoherent control. We prove the following.

######  Theorem 3.

After the $n$th round, ground-state population of the target system under
incoherent control can reach

$p_{n}^{\mathrm{F},inc}=p_{*}^{inc}-v_{\mathrm{F}}^{n}(p_{*}^{inc}-\gamma),$
(16)

$\forall\mathrm{F}=\mathrm{TP},\mathrm{MTP},\mathrm{MMTP}^{(d)}$. the
asymptotic ground-state population reads

$p_{*}^{inc}=\frac{1}{1+e^{-\beta\mathcal{E}}e^{\beta_{H}(\mathcal{E}-E)}}.$
(17)

The convergence rates for $\mathrm{F}=\mathrm{TP},\mathrm{MTP}$, and
$\mathrm{MMTP}^{(d)}$ are respectively

$\displaystyle v_{\mathrm{TP}}$ $\displaystyle=$
$\displaystyle\eta(1-e^{-\beta\mathcal{E}}),$ (18) $\displaystyle
v_{\mathrm{MTP}}$ $\displaystyle=$ $\displaystyle
v_{\mathrm{TP}}+\frac{1-\eta+\eta
e^{-\beta\mathcal{E}}}{1+e^{\beta\mathcal{E}}},$ (19) $\displaystyle
v_{\mathrm{MMTP}^{(d)}}$ $\displaystyle=$ $\displaystyle
v_{\mathrm{TP}}+\frac{\eta-(2\eta-1)}{1+e^{-\beta\mathcal{E}}}\delta_{d}(\gamma_{\mathcal{E}}),$
(20)

where $\gamma=1/(1+e^{-\beta E})$, $\eta=1/(1+e^{-\beta_{H}(\mathcal{E}-E)})$,
and $\gamma_{\mathcal{E}}=1/(1+e^{-\beta\mathcal{E}})$. Further, for
$\mathrm{F}=\mathrm{TP},\mathrm{MTP}$, Eq. (16) is also the upper bound for
ground-state population.

It can be inferred from the theorem that, cooling processes under
$\mathrm{TP}$, $\mathrm{MTP}$, and $\mathrm{MMTP}^{(d)}$ can all approach
$p_{*}^{inc}$ in the asymptotic limit, while the rates of convergence are
different. Further, for $H$ in the high temperature limit
$\beta_{H}\rightarrow 0$, the asymptotic ground-state population $p_{*}^{inc}$
reaches its maximum. Then we arrive at the following corollary.

######  Corollary 1.

Under incoherent control, in the limit of infinity cycles, the ground-state
population of a qubit target system is upper bounded by

$p_{*}=\frac{1}{1+e^{-\beta\mathcal{E}}}.$ (21)

where $\mathcal{E}$ stands for the largest energy gap of the composed system
consisting of the target qubit $S$ and the auxiliary qubit $A$.

###  III.2 Cooling under coherence control

In the $n$th round, the state of the system first transforms under the unitary
$\sigma_{x}:=|g\rangle\langle e|+|e\rangle\langle g|$, and then undergoes a
process $\Lambda^{\mathrm{F}}\in\mathrm{F}$. Precisely, the ground-state
population after the $n$th round can be expressed as

$\displaystyle p_{n}^{\mathrm{F},coh}$ $\displaystyle=$ $\displaystyle\langle
g|\Lambda^{\mathrm{F}}\left((1-p_{n-1}^{\mathrm{F},coh})|g\rangle\langle
g|+p_{n-1}^{\mathrm{F},coh}|e\rangle\langle e|\right)|g\rangle,$ (22)
$\displaystyle=$ $\displaystyle
p^{\mathrm{F}}_{g|g}+(p^{\mathrm{F}}_{g|e}-p^{\mathrm{F}}_{g|g})p_{n-1}^{\mathrm{F},coh}.$

Because $p_{0}^{\mathrm{F},coh}=\gamma$, it follows that

$p_{n}^{\mathrm{F},coh}=p_{*}^{\mathrm{F},coh}-v^{n}_{\mathrm{F},coh}(p_{*}^{\mathrm{F},coh}-\gamma),$
(23)

where
$p_{*}^{\mathrm{F},coh}=p^{\mathrm{F}}_{g|g}/(p^{\mathrm{F}}_{g|g}+p^{\mathrm{F}}_{e|e})$
and $v_{\mathrm{F},coh}=p^{\mathrm{F}}_{g|e}-p^{\mathrm{F}}_{g|g}$. Because
$|v_{\mathrm{F},coh}|<1$ as long as $\Lambda^{\mathrm{F}}$ is not an identical
operation, the ground-state population converges to $p_{*}^{\mathrm{F},coh}$
at a speed exponentially fast in the number of rounds.

For $\mathrm{F}=\mathrm{TP}$, Eq. (12) has been proved in [21]. For
consistency, we give the proof as follows. From Eq. (5), we have

$p_{*}^{\mathrm{TP},coh}=(1+\frac{1-\lambda}{1-\lambda e^{-\beta E}})^{-1},$
(24)

and $v_{\mathrm{TP},coh}=\lambda(1+e^{-\beta E})-1$ with $\lambda\in[0,1]$. It
follows that if one requires to cool the system to a temperature below that of
the reservoir, i.e., $p_{*}^{\mathrm{TP},coh}\geq\gamma$, the thermal process
employed should satisfy $\lambda\geq\gamma$. Further, from Lemma 1, in each
round, $p_{n}^{\mathrm{TP},coh}$ is maximized at $\lambda=1$, as long as
$1-p_{n-1}^{\mathrm{TP},coh}\leq\gamma$, which naturally holds. For
$\lambda=1$, we have $p_{*}^{\mathrm{TP},coh}=1$ and
$v_{\mathrm{TP},coh}=e^{-\beta E}$. It then follows from Eq. (23) that the
maximum of $p_{n}^{\mathrm{TP},coh}$ is in the form of Eq. (12). A direct
consequence is that, under coherence control and with
$\mathrm{F}=\mathrm{TP}$, the target qubit can be cooled to absolute zero in
the asymptotic limit .

In contrast, when $\mathrm{F}=\mathrm{MTP}$, the system can only be cooled to
the temperature of the reservoir. This is because MTP acting on a qubit system
is equivalent to partial thermalizations.

When $\mathrm{F}=\mathrm{MMTP}^{(d)}$, we apply simulated $\beta$-swap on the
target qubit and thermalize the memory before each round. The ground-state
population after the $n$th round can be calculated by using Eq. (11) and
setting $p_{i}=1-p_{n-1}^{(d),coh}$ and $p_{j}=1-p_{i}$. Precisely, We obtain
the following

$\displaystyle p^{(d),coh}_{n}$ $\displaystyle=$ $\displaystyle[e^{-\beta
E}-\delta_{d}(\gamma)]p_{n-1}^{(d),coh}$ (25) $\displaystyle+1-e^{-\beta
E}+(1-\gamma)\delta_{d}(\gamma)$ $\displaystyle=$ $\displaystyle
p_{\max}^{(d),coh}-[e^{-\beta
E}-\delta_{d}(\gamma)](p_{\max}^{(d),coh}-p_{n-1}^{(d),coh})$ $\displaystyle=$
$\displaystyle p_{\max}^{(d),coh}-[e^{-\beta
E}-\delta_{d}(\gamma)]^{n}(p_{\max}^{(d),coh}-\gamma),$

where $p_{\max}^{(d),coh}$ is in the form of Eq. (15). This completes the
proof of Theorem 2.

Careful analysis shows that $p_{\max}^{(d),coh}>\gamma$, $\forall d\geq 2$,
indicating that cooling is enhanced by the memory. Further,
$p_{\max}^{(d),coh}$ is monotonically increasing in $d$, which means that the
enhancement grows with the size of the memory. Last but not least,
$p_{\max}^{(1),coh}=p_{*}^{\mathrm{MTP},coh}=\gamma$ and
$p_{\max}^{(\infty),coh}=p_{*}^{\mathrm{TP},coh}=1$, and this bridges the gap
between the performances of TP and MTP.

###  III.3 Cooling under incoherent control

Initially, the target qubit is fully thermalized in $R$ while the auxiliary is
fully thermalized in $H$. On basis
$\\{|g0_{A}\rangle,|g1_{A}\rangle,|e0_{A}\rangle,|e1_{A}\rangle\\}$, the
initial state of $SA$ can be written as

$\displaystyle\bm{p^{(0)}}$ $\displaystyle=$
$\displaystyle[p^{(0)}_{g0},p^{(0)}_{g1};p^{(0)}_{e0},p^{(0)}_{e1}]$ (26)
$\displaystyle=$
$\displaystyle[\gamma\eta,\gamma(1-\eta);(1-\gamma)\eta,(1-\gamma)(1-\eta)],$

where $\eta=1/(1+e^{-\beta_{H}(\mathcal{E}-E)})$. The target in the first
round is then to maximize the sum of the occupations on $|g0_{A}\rangle$ and
$|g1_{A}\rangle$ by a quantum operation in $\mathrm{F}$.

We first consider the case $\mathrm{F}=\mathrm{TP}$. Let
$\bm{\tau}=[\gamma,1-\gamma]\otimes[\gamma_{A},1-\gamma_{A}]$ with
$\gamma_{A}=1/(1+e^{-\beta(\mathcal{E}-E)})$ be the Gibbs state of $SA$ in the
reservoir $R$. Because $\beta_{H}<\beta$, we have

$\frac{p^{(0)}_{g1}}{\tau_{g1}}=\frac{p^{(0)}_{e1}}{\tau_{e1}}>\frac{p^{(0)}_{g0}}{\tau_{g0}}=\frac{p^{(0)}_{e0}}{\tau_{e0}}.$
(27)

It means that one can neither increase the population on $|g1_{A}\rangle$, nor
decrease the population on $|e0_{A}\rangle$ via any operation in TP. The
optimal strategy is then the $\beta$-swap between the $|g0_{A}\rangle$ and
$|e1_{A}\rangle$. The ground-state population of $S$ then becomes

$p_{1}^{\mathrm{TP},inc}=1-\eta+\gamma\eta(1-e^{-\beta\mathcal{E}}).$ (28)

Equivalently,

$(p_{*}^{inc}-p_{1}^{\mathrm{TP},inc})=v_{\mathrm{TP}}(p_{*}^{inc}-\gamma),$
(29)

where $p_{*}^{inc}$ and $v_{\mathrm{TP}}$ are defined in Eqs. (17) and (18)
respectively. Notice that $p_{*}^{inc}>p_{1}^{\mathrm{TP},inc}$, because
$p_{*}^{inc}>\gamma$ and $v_{\mathrm{TP}}>0$.

In the second round, the auxiliary is first fully thermalized in $H$,
resulting in state $[\eta,1-\eta]$. The joint state of $SA$ becomes

$\displaystyle\bm{p^{(1)}}$ $\displaystyle=$
$\displaystyle[p^{(1)}_{g0},p^{(1)}_{g1};p^{(1)}_{e0},p^{(1)}_{e1}]$ (30)
$\displaystyle=$
$\displaystyle[p_{1}^{\mathrm{TP},inc}\eta,p_{1}^{\mathrm{TP},inc}(1-\eta);$
$\displaystyle(1-p_{1}^{\mathrm{TP},inc})\eta,(1-p_{1}^{\mathrm{TP},inc})(1-\eta)].$

Then we check that

$\frac{p^{(1)}_{g1}}{\tau_{g1}}>\frac{p^{(1)}_{e1}}{\tau_{e1}}>\frac{p^{(1)}_{g0}}{\tau_{g0}}>\frac{p^{(1)}_{e0}}{\tau_{e0}}.$
(31)

Therefore, the best strategy is still the $\beta$-swap between the
$|g0_{A}\rangle$ and $|e1_{A}\rangle$. The resulted ground-state population
$p_{2}^{\mathrm{TP},inc}$ then satisfies

$(p_{*}^{inc}-p_{2}^{\mathrm{TP},inc})=v_{\mathrm{TP}}(p_{*}^{inc}-p_{1}^{\mathrm{TP},inc}).$
(32)

We apply the above discussion $n$ times and obtain the recurrence relation

$(p_{*}^{inc}-p_{n}^{\mathrm{TP},inc})=v_{\mathrm{TP}}(p_{*}^{inc}-p_{n-1}^{\mathrm{TP},inc}).$
(33)

Importantly, this recurrence relation implies
$p_{*}^{inc}>p_{n}^{\mathrm{TP},inc}$ as long as
$p_{*}^{inc}>p_{n-1}^{\mathrm{TP},inc}$, which in turn ensures that the
ordering in Eq. (31) holds for $p^{(n)}$. This is the reason for employing the
$\beta$-swap between $|g0_{A}\rangle$ and $|e1_{A}\rangle$ as the optimal
cooling strategy in every round. Combining Eqs. (29) and (33), we arrive at
Eq. (16) in Theorem 3 for $\mathrm{F}=\mathrm{TP}$.

As for $\mathrm{F}=\mathrm{MTP}$, the discussions are similar. The only
difference is that $\beta$-swap is not in $\mathrm{MTP}$. The optimal strategy
in each cycle then becomes full thermalization between $|g0_{A}\rangle$ and
$|e1_{A}\rangle$. Therefore, we have the following recurrence relation

$\displaystyle p_{n}^{\mathrm{MTP},inc}$ $\displaystyle=$
$\displaystyle\gamma[p_{n-1}^{\mathrm{MTP},inc}\eta+(1-p_{n-1}^{\mathrm{MTP},inc})(1-\eta)]$
(34) $\displaystyle+(1-\eta)p_{n-1}^{\mathrm{MTP},inc}.$

where $p_{0}^{\mathrm{MTP},inc}=\gamma$. This is equivalent to

$(p_{*}^{inc}-p_{n}^{\mathrm{MTP},inc})=v_{\mathrm{MTP}}(p_{*}^{inc}-p_{n-1}^{\mathrm{MTP},inc}).$
(35)

Equation (16) then follows directly for $\mathrm{F}=\mathrm{MTP}$.
Interestingly, even though the allowed thermal process is restricted to be
Markovian, one can still approach $p_{*}^{inc}$ in the asymptotic limit. The
difference between the performances of $\mathrm{TP}$ and $\mathrm{MTP}$ is the
convergence rate. Because $v_{\mathrm{MTP}}>v_{\mathrm{TP}}$,
$p_{n}^{\mathrm{TP},inc}$ converges faster than $p_{n}^{\mathrm{MTP},inc}$.

For $\mathrm{F}=\mathrm{MMTP}^{(d)}$, our strategy is to implement simulated
$\beta$-swap between $|g0_{A}\rangle$ and $|e1_{A}\rangle$ in each cycle. From
Theorem 1, where we set $p_{i}=\eta p_{n-1}^{(d),inc}$ and
$p_{j}=(1-\eta)(1-p_{n-1}^{(d),inc})$, the sum of populations on
$|g0_{A}\rangle$ and $|g1_{A}\rangle$ after the $n$th round then reads

$\displaystyle p_{n}^{(d),inc}$ $\displaystyle=$
$\displaystyle(1-\eta)p_{n-1}^{(d),inc}+(1-e^{-\beta\mathcal{E}})\eta
p_{n-1}^{(d),inc}$ (36)
$\displaystyle+(1-\eta)(1-p_{n-1}^{(d),inc})[(1-\gamma_{\mathcal{E}})\eta
p_{n-1}^{(d),inc}$
$\displaystyle-\gamma_{\mathcal{E}}(1-\eta)(1-p_{n-1}^{(d),inc})]\delta_{d}(\gamma_{\mathcal{E}}),$
$\displaystyle=$ $\displaystyle
p_{*}^{inc}-v_{\mathrm{MMTP}^{(d)}}(p_{*}^{inc}-p_{n-1}^{(d),inc}),$
$\displaystyle=$ $\displaystyle
p_{*}^{inc}-v_{\mathrm{MMTP}^{(d)}}^{n}(p_{*}^{inc}-\gamma),$

with $v_{\mathrm{MMTP}^{(d)}}$ in the form of Eq. (20). This completes the
proof for $\mathrm{F}=\mathrm{MMTP}^{(d)}$. Still, the asymptotic limit of
$p_{n}^{(d),inc}$ is the same as that of $p_{n}^{\mathrm{TP},inc}$ and
$p_{n}^{\mathrm{MTP},inc}$. Further, $v_{\mathrm{MMTP}^{(d)}}$ is
monotonically decreasing in $d$. It bridges the gap between the convergence
rates for $\mathrm{TP}$ and $\mathrm{MTP}$, in the sense that
$v_{\mathrm{MMTP}^{(1)}}=v_{\mathrm{MTP}}$ and
$v_{\mathrm{MMTP}^{(\infty)}}=v_{\mathrm{TP}}$.

##  IV Work extraction

Consider a qubit $S$ initially in a non-equilibrium state $[p_{0},1-p_{0}]$,
and a work bit with Hamiltonian $H_{{}_{\tiny W}}={\tiny W}|1\rangle\langle
1|$ initially in the ground state $[1,0]$. The task of work extraction is to
transform the state of the work bit to its exited state $|1\rangle$ by joint
thermal processes on $SW$. Precisely, in a single-shot work extraction,

$[p_{0},1-p_{0}]_{S}\otimes[1,0]_{{\tiny
W}}\stackrel{{\scriptstyle\mathrm{F}}}{{\longmapsto}}[\gamma,1-\gamma]_{S}\otimes[\epsilon,1-\epsilon]_{{\tiny
W}},$ (37)

where $\mathrm{F}=\mathrm{TP},\mathrm{MTP},\mathrm{ETP}$, or
$\mathrm{MMTP}^{(d)}$, and $\epsilon$ is called the error of work extraction.
The task of work extraction is to minimize $\epsilon$ over $\mathrm{F}$.

Here we analytically solve the above problem for the case where the target
qubit is in the excited state, i.e., $p_{0}=0$, and prove the following
theorem.

######  Theorem 4.

Consider a qubit $S$ with Hamiltonian $H_{S}=E|e\rangle\langle e|$ and
initially in state $[0,1]_{S}$, and a work bit with Hamiltonian
$H_{{}_{W}}=W|1\rangle\langle 1|$ and initially in state $[1,0]_{W}$. By
$\mathrm{MMTP}^{(d)}$, the error $\epsilon^{(d)}$ of work extraction can reach

$\epsilon^{(d)}(W)=I_{d}(\frac{1}{1+e^{-\beta(E-W)}},\frac{1}{1+e^{\beta
W}}).$ (38)

Moreover, for $d\rightarrow\infty$, the error can approach

$\displaystyle\epsilon^{(\infty)}(W)=\epsilon^{\mathrm{TP}}(W)$ (39)
$\displaystyle=$ $\displaystyle\left\\{\begin{matrix}0,&{\tiny W}\leq
W_{0},\\\ 1-e^{\beta(E-W)}-e^{-\beta{\tiny W}},&{\tiny
W}>W_{0}.\end{matrix}\right.$

where $Z\equiv 1+e^{-\beta E}$ is the partition function of the target qubit,
$W_{0}\equiv E+k_{B}T\ln Z$ is the maximal work which can be extracted from
$S$ via $\mathrm{TP}$ with vanishing error, and $\epsilon^{\mathrm{TP}}(W)$ is
the minimum error that can be reached via $\mathrm{TP}$ for given $W$.

Here the function $I_{d}(x,y)$ is defined in Eq. (88) in Appendix A.

This theorem provides an evidence in the regime of high-dimensional systems
for the conjecture that $\mathrm{MMTP}^{(d)}$ can approach TP for
$d\rightarrow\infty$. This conjecture is not trivial, because if it is proved,
then the thermal bath can be divided into two parts. The first one can
exchange energy with the system but does not have memory effect, while the
second one serves as a memory and does not exchange energy with the system.

###  IV.1 Error of work extraction via TP, MTP and ETP

In the following, we derive the minimum error of work extraction as a function
of $W$, when the allowed operations are TP, MTP and ETP. In particular, for
comparison with the performance of $\mathrm{MMTP}^{(d)}$, we focus on the case
with $p_{0}=0$. On basis $\\{|g0\rangle,|g1\rangle,|e0\rangle,|e1\rangle\\}$,
the initial state of $SW$ then reads $\bm{p}_{0}=[0,0,1,0]$. The task of work
extraction is then to minimize the sum of populations on $|g0\rangle$ and
$|e0\rangle$ in the output by the allowed thermal processes.

First, we consider the case with $\mathrm{F}=\mathrm{TP}$. The optimal
extraction error $\epsilon^{\mathrm{TP}}(W)$ can be computed via the thermo-
majorization condition

$[0,1]_{S}\otimes[1,0]_{{\tiny
W}}\succ_{{}_{T}}[\gamma,1-\gamma]_{S}\otimes[\epsilon,1-\epsilon]_{{\tiny
W}}.$ (40)

This condition leads to

$\epsilon_{\mathrm{TP}}(W)=\left\\{\begin{matrix}0,&{\tiny W}\leq W_{0},\\\
1-e^{\beta(E-W)}-e^{-\beta{\tiny W}},&{\tiny W}>W_{0}.\end{matrix}\right.$
(41)

where $W_{0}$ is defined in theorem 4. The optimal strategy is to implement a
joint thermal operation described by the transition matrix $G^{\mathrm{TP}}$,
followed by a local full thermalization on $S$. The form of $G^{\mathrm{TP}}$
reads as follows.
For $W\leq E$,

$G^{\mathrm{TP}}_{1}=\left(\begin{matrix}1&0&0&0\\\ 0&1-e^{-\beta(E-W)}&1&0\\\
0&e^{-\beta(E-W)}&0&0\\\ 0&0&0&1\end{matrix}\right).$ (42)

For $E<W\leq W_{0}$, which implies $e^{\beta E}<e^{\beta W}\leq 1+e^{\beta
E}$,

$G^{\mathrm{TP}}_{2}=\left(\begin{matrix}1&0&0&0\\\ 0&0&e^{-\beta(W-E)}&0\\\
0&1&0&e^{\beta W}-e^{\beta E}\\\ 0&0&1-e^{-\beta(W-E)}&1-e^{\beta W}+e^{\beta
E}\end{matrix}\right).$ (43)

For $W>W_{0}$, which implies $e^{\beta W}>1+e^{\beta E}$,

$G^{\mathrm{TP}}_{3}=\left(\begin{matrix}1&0&0&0\\\ 0&0&e^{-\beta(W-E)}&0\\\
0&1&1-e^{-\beta W}(1+e^{\beta E})&1\\\ 0&0&e^{-\beta W}&0\end{matrix}\right).$
(44)

The above form of $G^{\mathrm{TP}}$ minimizes the sum of $G_{g0|e0}$ and
$G_{e0|e0}$, and thus minimizes $\epsilon$ in the output.

Next we study the performance of MTP. The optimal strategy reads

$G^{\mathrm{MTP}}=T_{S}T^{g1,e0}T^{g1,e1},$ (45)

or

$\tilde{G}^{\mathrm{MTP}}=T_{S}T^{g1,e1}T^{g1,e0},$ (46)

where $T_{S}$ denotes full thermalization of system $S$. Both
$G^{\mathrm{MTP}}$ and $\tilde{G}^{\mathrm{MTP}}$ give

$\epsilon_{\mathrm{MTP}}(W)=\frac{1}{(1+e^{\beta(E-{\tiny
W}})(1+e^{-\beta{\tiny W}})}.$ (47)

Notice that, for ${\tiny W}\to 0$,
$\epsilon_{\mathrm{MTP}}(W)\to\frac{1}{2}(1-\gamma)\neq 0$.

For $\mathrm{F}=\mathrm{ETP}$, we find that the optimal transition matrix is

$G^{\mathrm{ETP}}=T_{S}\beta^{g1,e0}\beta^{g1,e1},$ (48)

or

$\tilde{G}^{\mathrm{ETP}}=T_{S}\beta^{g1,e1}\beta^{g1,e0}.$ (49)

It then follows that

$\epsilon_{\mathrm{ETP}}(W)=\left\\{\begin{matrix}0,&{\tiny W}\leq E,\\\
(1-e^{\beta\delta})(1-e^{-\beta{\tiny W}}),&{\tiny W}>E.\end{matrix}\right.$
(50)

Figure 2:  The error $\bm{\epsilon}$ of work extraction as a function of
extracted work $\bm{W}$. Work is extracted from a qubit $S$ with energy gap
$E$ and initially in the excited state $[0,1]_{S}$.

Here We briefly mention that the Gibbs-preserving stochastic matrices
$G_{2}^{\mathrm{TP}}$ and $G_{3}^{\mathrm{TP}}$ in Eqs. (43) and (44) satisfy
the detailed balanced condition, i.e.,
$G_{j|i}=e^{-\beta(E_{j}-E_{i})}G_{i|j}$. However, these transformations are
impossible via ETP, because as shown in Fig. 2, the error of work extraction
achieved by these two stochastic matrices cannot be reached by ETP. This
observation shows that detailed balance is a necessary but not sufficient
condition for a thermal process to be realized as a sequence of thermal
processes which involve only two energy levels at a time.

###  IV.2 Memory-assisted protocol for work extraction

Here we propose a protocol for work extraction from a qubit system in the
excited state with the assistance of a $d$-dimensional memory. The basis of
the Hilbert space of the composed system $SWM$ reads
$\\{|g\rangle,|e\rangle\\}_{S}\otimes\\{|0\rangle,|1\rangle\\}_{W}\otimes\\{|1\rangle,\dots,|d\rangle\\}_{M}=\\{|\xi\zeta
k\rangle\\}_{\xi=g,e;\zeta=0,1;k=1,\dots,d}$. Initially, the system is in the
excited state $|e\rangle$, the work bit is in the ground state $|0\rangle$,
and the memory is in a maximally mixed state $\frac{\mathbb{1}}{d}$. Our
protocol is schematically depicted in Fig. 3.

Figure 3: Scheme for work extraction under $\bm{\mathrm{MMTP}^{(d)}}$. The
color of each dot stands for the population distribution. Initially, the qubit
$S$ is in the excited state and thus the composed state of $SW$ reads
$[0,0,1,0]_{SW}$, while the $d$-dimensional memory is in a maximally mixed
state. In step (I), MTP is applied to the subspace spanned by
$\\{|g11\rangle,\dots,|g1d\rangle;|e01\rangle,\dots,|e0d\rangle\\}$. In step
(II), MTP is applied to the subspace spanned by
$\\{|e01\rangle,\dots,|e0d\rangle;|e11\rangle,\dots,|e1d\rangle\\}$. The error
$\epsilon^{(d)}$ of extraction equals to the sum of populations on
$|e01\rangle,\dots,|e0d\rangle$.

Precisely, our protocol consist of two steps, which are summarized as

$\displaystyle[0,0,1,0]\otimes[\overbrace{\frac{1}{d},...,\frac{1}{d}}^{d\
times}]$
$\displaystyle=[{0,...,0};0,...,0;\frac{1}{d},...,\frac{1}{d};0,...,0]$
$\displaystyle\stackrel{{\scriptstyle\textbf{(\@slowromancap
i@)}}}{{\longmapsto}}[0,...,0;a_{d}^{(1)},...,a_{d}^{(d)};b_{1}^{(d)},...,b_{d}^{(d)};0,...,0]$
$\displaystyle\stackrel{{\scriptstyle\textbf{(\@slowromancap
ii@)}}}{{\longmapsto}}[0,...,0;a_{d}^{(1)},...,a_{d}^{(d)};\epsilon_{1},...,\epsilon_{d};c_{1}^{(d)},...,c_{d}^{(d)}].$

In step (I), simulated $\beta$-swap is implemented between
$|e\rangle_{S}\otimes|0\rangle_{W}$ and $|g\rangle_{S}\otimes|1\rangle_{W}$.
Notice that instead of fully thermalizing $M$ after the above operation, we
preserve the total state of $SWM$ for the second step. By combining Eqs.
(103), (105) and (107) together in Appendix B.2, we obtain

$\displaystyle b_{j}^{(d)}$ $\displaystyle=$
$\displaystyle\frac{(1-\gamma_{\delta})^{d}}{d}\sum^{j-1}_{j^{\prime}=0}f^{(j^{\prime})}_{d}\gamma_{\delta}^{j^{\prime}},$
(51)

where $\delta=W-E$, and $\gamma_{\delta}=1/(1+e^{-\beta\delta})$. See Appendix
D for detailed derivation. Notice that $b^{(d)}_{j}$ is monotonically
increasing with $j$.

Step (II) consists of $d$ subroutines. In the $k$th subroutine, full
elementary thermalizations are sequentially implemented between $|e0k\rangle$
and $|e1j\rangle$ with $j=1,\dots,d$. This results in

$\epsilon_{k}=\sum_{j=1}^{k}b_{k+j-1}^{(d)}\gamma_{W}^{d}(1-\gamma_{W})^{j-1}f_{d}^{(j-1)}.$
(52)

Consequently, we have

$\displaystyle\epsilon^{(d)}=\sum^{d}_{k=1}\epsilon_{k}=I_{d}(\gamma_{\delta},1-\gamma_{W}).$
(53)

The detailed derivation of the above formulas can be found in Appendix D.
Further, in Appendix A, we prove that

$\displaystyle\lim_{d\to\infty}I_{d}(x,y)$ (54) $\displaystyle=$
$\displaystyle\left\\{\begin{matrix}1-\frac{x}{1-x}-\frac{y}{1-y},&\frac{x}{1-x}+\frac{y}{1-y}<1,\\\
0,&\frac{x}{1-x}+\frac{y}{1-y}\geq 1.\end{matrix}\right.$

This directly leads to Eq. (39) in Theorem 4. The comparison between the
performances of TP, MTP, ETP and $\mathrm{MMTP}^{(d)}$ with finite $d$ is
plotted in Fig. 2. Interestingly, although $\epsilon^{(d)}(W)$ with finite $d$
is an analytic function, in the limit of $d\rightarrow\infty$, it approaches
$\epsilon^{\mathrm{TP}}(W)$, whose first derivative on $W$ is not continuous
at $W=W_{0}$.

##  V conclusion

We have introduced a finite-dimensional memory system with trivial Hamiltonian
and initially in a maximally mixed state to bridge the gap between the
performances of Markovian thermal processes and general thermal processes in
thermodynamic tasks, including cooling and work extraction. We study two
cooling paradigms, under coherent control and under incoherent control. In
both paradigms, the ground-state population of the target qubit approaches a
limit exponentially fast in the number of rounds, whether the set of allowed
processes is TP, MTP or $\mathrm{MMTP}^{(d)}$. For the paradigm under
coherence control, the limit of achievable ground-state population via
$\mathrm{MMTP}^{(d)}$ is monotonically increasing in $d$, which reduces to the
limit for MTP when $d=1$, and approaches the limit for TP when
$d\rightarrow\infty$. For the paradigm under incoherent control, we show that
the limits of achievable ground-state population are the same for TP, MTP, and
$\mathrm{MMTP}^{(d)}$, while the convergence rates are different. Still, the
rate for $\mathrm{MMTP}^{(d)}$ bridges the rates for TP and MTP.

In the problem of work extraction, the system which undergoes
$\mathrm{MMTP}^{(d)}$ consists of an out-of-equilibrium qubit $S$ and a work
qubit $W$. We derive the analytic expression for the error of work extraction
achieved by $\mathrm{MMTP}^{(d)}$ when $S$ is initially in the excited state.
Interestingly, the extraction error achieved by $\mathrm{MMTP}^{(d)}$, which
is an analytic function of the energy gap $W$ of the work bit, approaches in
the limit $d\rightarrow\infty$ to the extraction error achieved by TP, whose
first-order derivation is not continuous.

Our results provide evidences for the conjecture that $\mathrm{MMTP}^{(d)}$
can simulate transformations under TP with vanishing error, when the size $d$
of the memory is large enough.

###### Acknowledgements.

This work was supported by NSFC under Grant No. 11774205.

##  Appendix A SEVERALe FUNCTIONS

In this appendix, we give the definitions to several functions, and derive
some equivalent expressions for these functions.

######  Definition 2.

We define functions $L_{n}^{(m)}(x)$, $K_{n}^{(m)}(x)$, and $I_{n}^{(m)}(x)$
as follows:

$\displaystyle L_{n}^{(m)}(x)$ $\displaystyle:=$
$\displaystyle(1-x)^{n}\sum_{j=0}^{m}f_{n}^{(j)}x^{j},$ (55) $\displaystyle
K_{n}^{(m)}(x)$ $\displaystyle:=$
$\displaystyle\frac{(1-x)^{n}}{n}\sum_{j=0}^{m}jf_{n}^{(j)}x^{j},$ (56)
$\displaystyle I_{n}^{(m)}(x)$ $\displaystyle:=$
$\displaystyle\frac{(1-x)^{n}}{n}\sum_{j=0}^{m}(n-j)f_{n}^{(j)}x^{j},$ (57)

where $x\in[0,1]$, $m$ and $n$ are positive integers satisfying $m\leq n-1$,
and $f_{n}^{(j)}$ are defined by the following recurrence relation

$f_{j}^{(0)}=1,\ f_{j}^{(k+1)}=\sum_{j^{\prime}=1}^{j}f_{j^{\prime}}^{(k)},$
(58)

or equivalently,

$f_{j}{(k)}=C_{j-1+k}^{k}.$ (59)

By definition, $K_{n}^{(m)}(x)$ and $I_{n}^{(m)}(x)$ are related to
$L_{n}^{(m)}(x)$ as

$\displaystyle K_{n}^{(m)}(x)$ $\displaystyle=$
$\displaystyle\frac{x}{1-x}L_{n}^{(m)}(x)+\frac{x}{n}\frac{d}{dx}L_{n}^{(m)}(x),$
(60) $\displaystyle I_{n}^{(m)}(x)$ $\displaystyle=$ $\displaystyle
L_{n}^{(m)}(x)-K_{n}^{(m)}(x).$ (61)

######  Lemma 2.

The function $L_{n}^{(m)}(x)$ has the following equivalent expressions

$\displaystyle L_{n}^{(m)}(x)$ $\displaystyle=$ $\displaystyle
1-nC_{n+m}^{m}x^{m+1}\sum_{l=0}^{n-1}C_{n-1}^{l}\frac{(-x)^{l}}{m+l+1},$ (62)
$\displaystyle=$ $\displaystyle 1-nC_{n+m}^{m}\int_{0}^{x}x^{m}(1-x)^{n-1}dx.$
(63)

###### Proof.

In order to prove Eq. (62), we decompose $L_{n}^{(m)}(x)$ into a polynomial
$\sum_{k=0}^{n+m}c_{k}x^{k}$, and calculate the coefficient $c_{k}$.
(i) For $k=0$, we have $c_{0}=f_{n}^{(0)}=1$.
(ii) For $1\leq k\leq m$, direct calculation leads to
$c_{k}=\sum_{k^{\prime}=0}^{k}(-1)^{k^{\prime}}M_{k^{\prime}}^{(k)}$, where
$M_{k^{\prime}}^{(k)}=\sum_{k^{\prime}=0}^{k}C_{n}^{k^{\prime}}f_{n}^{(k-k^{\prime})}=\sum_{k^{\prime}=0}^{k}C_{n}^{k^{\prime}}C_{n+k-k^{\prime}-1}^{k-k^{\prime}}$.
Further, by the definition of $M_{k^{\prime}}^{(k)}$, the following equation
holds for $k\geq 1$ and $k^{\prime}=1,\dots,k$,

$\displaystyle
M_{k^{\prime}}^{(k)}=\Delta_{k^{\prime}-1}^{(k)}+\Delta_{k^{\prime}}^{(k)},$
(64)

where $\Delta_{0}^{(k)}=C_{n}^{0}f_{n}^{(k)}=f_{n}^{(k)}$, and
$\Delta_{k^{\prime}}^{(k)}=M_{k^{\prime}}^{(k)}\frac{(n-k^{\prime})(k-k^{\prime})}{nk}$.

Therefore, we get

$\displaystyle c_{k}$ $\displaystyle=$ $\displaystyle
M_{0}+\sum_{k^{\prime}=1}^{k}(-1)^{k^{\prime}}(\Delta_{k^{\prime}-1}+\Delta_{k^{\prime}})$
(65) $\displaystyle=$
$\displaystyle\sum_{k^{\prime}=0}^{k}(-1)^{k^{\prime}}\Delta_{k^{\prime}}+\sum_{k^{\prime}=0}^{k-1}(-1)^{k^{\prime}+1}\Delta_{k^{\prime}}$
(66) $\displaystyle=$ $\displaystyle\Delta_{k}=0.$ (67)

(iii) For $m+1\leq k\leq m+n$, let $k=m+l$ and thus, $0\leq l\leq n-1$. The
coefficient $c_{k}$ is calculated as

$\displaystyle c_{k}$ $\displaystyle=$
$\displaystyle\sum_{k^{\prime}=l+1}^{n}(-1)^{k^{\prime}}M_{k^{\prime}}$ (68)
$\displaystyle=$
$\displaystyle\sum_{k^{\prime}=0}^{n}(-1)^{k^{\prime}}M_{k^{\prime}}-\sum_{k^{\prime}=0}^{l}(-1)^{k^{\prime}}M_{k^{\prime}}$
(69) $\displaystyle=$ $\displaystyle 0-(-1)^{l}\Delta_{l}$ (70)
$\displaystyle=$
$\displaystyle(-1)^{l+1}\frac{n-l}{n}C_{n}^{l}\frac{n(m+1)}{n(m+l+1)}C_{n+m}^{m+1}$
(71) $\displaystyle=$
$\displaystyle(-1)^{l+1}nC_{n+m}^{m}\frac{C_{n-1}^{l}}{m+l+1}.$ (72)

This leads to the expression as Eq. (62).

Further, from Eq. (62) we have

$\displaystyle L_{n}^{(m)}(0)$ $\displaystyle=$ $\displaystyle 1,$
$\displaystyle\frac{d}{dx}L_{n}^{(m)}(x)$ $\displaystyle=$ $\displaystyle-
nC_{n+m}^{m}\sum_{l=0}^{n-1}C_{n-1}^{l}(-1)^{l}x^{m+l}$ (73) $\displaystyle=$
$\displaystyle-nC_{n+m}^{m}x^{m}(1-x)^{n-1}.$

This leads to the expression as Eq. (63).

∎

######  Lemma 3.

For $x\in(0,\frac{1}{2})$,

$I_{n}^{(n-1)}(x)=\frac{1-2x}{1-x}+x\delta_{n}(1-x),$ (74)

with
$\delta_{n}(1-x):=\sum_{n^{\prime}=n}^{\infty}\frac{C_{2n^{\prime}}^{n^{\prime}}}{n^{\prime}+1}[x(1-x)]^{n^{\prime}}<o[(4x(1-x))^{n}n^{-3/2}]$.

###### Proof.

From Eq.(57), $I^{(0)}_{1}(x)=1-x$, and

$\displaystyle I^{(n-1)}_{n}(x)$ $\displaystyle=$
$\displaystyle(1-x)^{n}[1+C_{n-1}^{1}x+\frac{1}{n}\sum_{j=2}^{n-1}(n-j)C_{n+j-1}^{j}x^{j}],$
(75) $\displaystyle=$
$\displaystyle(1-x)^{n}[(1+x)^{n-1}+\frac{1}{n}\sum_{j=2}^{n-1}(n-j)C_{n+j-1}^{j}x^{j}-\sum_{j=2}^{n-1}C_{n-1}^{j}x^{j}],$
$\displaystyle=$
$\displaystyle(1-x)^{n}[(1+x)^{n-1}+\sum_{j=2}^{n-1}(\frac{n-j}{n}C_{n+j-1}^{j}-C_{n-1}^{j})x^{j}],$
$\displaystyle I^{(n)}_{n+1}(x)$ $\displaystyle=$
$\displaystyle(1-x)^{n}[(1-x)(1+x)(1+x)^{n-1}+(1-x)\sum_{j=2}^{n}(\frac{n+1-j}{n+1}C_{n+j}^{j}-C_{n}^{j})x^{j}].$
(76)

It follows that

$\displaystyle I^{(n)}_{n+1}(x)-I^{(n-1)}_{n}(x)$ $\displaystyle=$
$\displaystyle(1-x)^{n}\\{-x^{2}(1+x)^{n-1}+\sum_{j=2}^{n-1}[(1-x)(\frac{n+1-j}{n+1}C_{n+j}^{j}-C_{n}^{j})$
$\displaystyle-(\frac{n-j}{n}C_{n+j-1}^{j}-C_{n-1}^{j})]x^{j}+(1-x)x^{n}(\frac{C^{n}_{2n}}{n+1}-1)\\},$
$\displaystyle=$ $\displaystyle
x^{2}(1-x)^{n}\\{-\sum^{n-1}_{j=0}C^{j}_{n-1}x^{j}+\sum_{j=0}^{n-3}[(1-x)(\frac{n-1-j}{n+1}C_{n+j+2}^{j+2}-C_{n}^{j+2})$
$\displaystyle-(\frac{n-2-j}{n}C_{n+j+1}^{j+2}-C_{n-1}^{j+2})]x^{j}+(1-x)x^{n-2}(\frac{C^{n}_{2n}}{n+1}-1)\\},$
$\displaystyle=$ $\displaystyle
x^{2}(1-x)^{n}\\{\sum_{j=0}^{n-3}[-x(\frac{n-(j+1)}{n+1}C_{n+j+2}^{j+2}-C_{n}^{j+2})$
$\displaystyle-(-\frac{n-j}{n+1}C_{n+j+1}^{j+1}+C^{j+1}_{n})]x^{j}-(n-1)x^{n-2}-x^{n-1}+(1-x)x^{n-2}(\frac{C^{n}_{2n}}{n+1}-1)\\},$
$\displaystyle=$ $\displaystyle
x^{2}(1-x)^{n}\\{\sum_{j=0}^{n-3}(C_{n}^{j+2}-\frac{n-(j+1)}{n+1}C_{n+j+2}^{j+2})x^{j+1}+\sum_{j=0}^{n-3}(\frac{n-j}{n+1}C_{n+j+1}^{j+1}-C^{j+1}_{n})x^{j}$
$\displaystyle-(n-1)x^{n-2}-x^{n-1}+(1-x)x^{n-2}(\frac{C^{n}_{2n}}{n+1}-1)\\},$
$\displaystyle=$ $\displaystyle
x^{2}(1-x)^{n}\\{(n-\frac{2}{n+1}C^{n-1}_{2n-1})x^{n-2}-(n-1)x^{n-2}-x^{n-1}+(1-x)x^{n-2}(\frac{C^{n}_{2n}}{n+1}-1)\\},$
$\displaystyle=$ $\displaystyle
x^{2}(1-x)^{n}\\{-x^{n-1}-x^{n-1}(\frac{C^{n}_{2n}}{n+1}-1)\\},$
$\displaystyle=$ $\displaystyle-x^{n+1}(1-x)^{n}\frac{C^{n}_{2n}}{n+1}.$ (78)

In obtaining Eq.(3), we used the relations
$\frac{n-2-j}{n}C_{n+j+1}^{j+2}=\frac{n-1-j}{n+1}C_{n+j+2}^{j+2}-\frac{n-j}{n+1}C_{n+j+1}^{j+1}$
and
$C_{n-1}^{j}-C_{n-1}^{j+2}=C_{n-1}^{j+1}+C_{n-1}^{j}-(C_{n-1}^{j+2}+C_{n-1}^{j+1})=C_{n}^{j+1}-C_{n}^{j+2}$.
Furthermore, for $x\in(0,\frac{1}{2})$, the following equation holds

$\sum^{\infty}_{n=1}\frac{C^{n}_{2n}}{n+1}[x(1-x)]^{n}=\frac{x}{1-x}.$ (79)

The proof is as follows. Let $x(1-x)=t$, and for $x\in(0,\frac{1}{2})$, we
have $x=\frac{1-\sqrt{1-4t}}{2}$. Consequently, Eq. (79) becomes

$\sum^{\infty}_{n=1}\frac{C^{n}_{2n}}{n+1}t^{n}=\frac{1-2t-\sqrt{1-4t}}{2t}.$
(80)

Taylor expansion of $\sqrt{1-4t}$ at $t=0$ gives

$\displaystyle\sqrt{1-4t}$ $\displaystyle=$ $\displaystyle
1-2t-2\sum^{\infty}_{n=2}\frac{(2n-2)!}{n!(n-1)!}t^{n},$ (81) $\displaystyle=$
$\displaystyle 1-2t-2\sum^{\infty}_{n=1}\frac{(2n)!}{n!(n+1)!}t^{n+1},$
$\displaystyle=$ $\displaystyle
1-2t-2\sum^{\infty}_{n=1}\frac{C^{n}_{2n}}{n+1}t^{n+1},$

which is equivalent to Eq. (80), and thus we arrive at Eq. (79).

Therefore, for $x\in(0,\frac{1}{2})$,

$\displaystyle I^{(n-1)}_{n}(x)$ $\displaystyle=$ $\displaystyle
I^{0}_{1}(x)+\sum^{n-1}_{n^{\prime}=1}I^{(n)}_{n+1}(x)-I^{(n-1)}_{n}(x),$
$\displaystyle=$ $\displaystyle
1-x-x\sum^{n-1}_{n^{\prime}=1}\frac{C^{n^{\prime}}_{2n^{\prime}}}{n^{\prime}+1}[x(1-x)]^{n^{\prime}},$
$\displaystyle=$ $\displaystyle
1-x-x(\frac{x}{1-x}-\sum^{\infty}_{n^{\prime}=n}\frac{C^{n^{\prime}}_{2n^{\prime}}}{n^{\prime}+1}[x(1-x)]^{n^{\prime}}),$
$\displaystyle=$ $\displaystyle\frac{1-2x}{1-x}+x\delta_{n}(1-x).$

For $n\gg 0$, by Stirling’s approximation $n!\approx\sqrt{2\pi
n}(\frac{n}{e})^{n}$, we have

$\displaystyle\delta_{n}(1-x)$ $\displaystyle=$
$\displaystyle\sum_{n^{\prime}=n}^{\infty}\frac{C_{2n^{\prime}}^{n^{\prime}}}{n^{\prime}+1}[x(1-x)]^{n^{\prime}},$
$\displaystyle\approx$
$\displaystyle\sum_{n^{\prime}=n}^{\infty}\frac{\sqrt{2\pi(2n^{\prime})}}{2\pi
n^{\prime}}\frac{(\frac{2n^{\prime}}{e})^{2n^{\prime}}}{(\frac{n^{\prime}}{e})^{2n^{\prime}}}\frac{1}{n^{\prime}+1}[x(1-x)]^{n^{\prime}},$
$\displaystyle=$ $\displaystyle\sum_{n^{\prime}=n}^{\infty}\frac{1}{\sqrt{\pi
n^{\prime}}(n^{\prime}+1)}[4x(1-x)]^{n^{\prime}},$ $\displaystyle<$
$\displaystyle\sum_{n^{\prime}=n}^{\infty}\frac{1}{\sqrt{\pi
n}(n+1)}[4x(1-x)]^{n^{\prime}}$ $\displaystyle<$
$\displaystyle\frac{1}{\sqrt{\pi}n^{3/2}}\frac{[4x(1-x)]^{n}}{(2x-1)^{2}}=o[(4x(1-x))^{n}n^{-3/2}].$

∎

######  Lemma 4.

In the limit $n\to\infty$, the function $L_{n}^{(m)}(x)$ approaches a
discontinuous function as

$\lim_{n\to\infty}L_{n}^{(m)}(x)=\left\\{\begin{matrix}1,&x<\frac{m}{m+n-1},\\\
0,&x>\frac{m}{m+n-1}.\end{matrix}\right.$ (83)

###### Proof.

By definition, we have $L_{n}^{(m)}(0)=1$ and $L_{n}^{(m)}(1)=0$, $\forall
m,n$.

Now let
$g^{(m)}_{n}(x):=-\frac{d}{dx}L_{n}^{(m)}(x)=nC_{n+m}^{m}x^{m}(1-x)^{n-1}$.

For $n\gg 0$ and $m$ being finite, or equivalently,
$\frac{m}{m+n-1}\rightarrow 0_{+}$, we have

$\displaystyle g^{(m)}_{n}(x)$ $\displaystyle=$
$\displaystyle\frac{n}{m!}\frac{(n+m)!}{n!}x^{m}(1-x)^{n-1},$
$\displaystyle\approx$
$\displaystyle\frac{n}{m!}\sqrt{\frac{n+m}{n}}\frac{(\frac{n+m}{e})^{m+n}}{(\frac{n}{e})^{n}}x^{m}(1-x)^{n-1},$
$\displaystyle=$
$\displaystyle\frac{n}{m!}\sqrt{\frac{n+m}{n}}(\frac{n+m}{e})^{m}(\frac{n+m}{n})^{n}x^{m}(1-x)^{n-1},$
$\displaystyle\approx$
$\displaystyle\frac{n}{m!}\sqrt{\frac{n+m}{n}}(n+m)^{m}x^{m}(1-x)^{n-1}.$

In the second line, we used Stirling’s formula, and in the fourth line, we
used $(\frac{n+m}{n})^{n}\approx e^{m}$ for $n\gg 0$. The above equation leads
to $\lim_{n\rightarrow\infty}g^{(m)}_{n}(x)=0$, $\forall x\in(0,1]$.
Therefore,

$\lim_{n\to\infty}L_{n}^{(m)}(x)=\left\\{\begin{matrix}1,&x=0,\\\ 0,&0<x\leq
1.\end{matrix}\right.$ (84)

It means that Eq. (83) holds for $\frac{m}{m+n-1}\rightarrow 0_{+}$.

For $n\gg 0$ and $\frac{m}{m+n-1}\neq 0$, which means that $m\gg 0$,

$\displaystyle g^{(m)}_{n}(x)$ $\displaystyle=$
$\displaystyle(n+m)\frac{(n-1+m)!}{m!(n-1)!}x^{m}(1-x)^{n-1},$
$\displaystyle\approx$
$\displaystyle\frac{n+m}{\sqrt{2\pi}}\sqrt{\frac{n-1+m}{m(n-1)}}\frac{(n-1+m)^{n-1+m}}{m^{m}(n-1)^{n-1}}x^{m}(1-x)^{n-1},$
$\displaystyle=$
$\displaystyle\frac{n+m}{\sqrt{2\pi}}\sqrt{\frac{n-1+m}{m(n-1)}}[\frac{n-1+m}{m}x]^{m}[\frac{n-1+m}{n-1}(1-x)]^{n-1},$
$\displaystyle=$
$\displaystyle\frac{n+m}{\sqrt{2\pi}}\sqrt{\frac{n-1+m}{m(n-1)}}[(1+\frac{1}{\mu})^{\mu}x^{\mu}(1+\mu)(1-x)]^{n-1}.$

where $\mu=\frac{m}{n-1}\leq 1$. Notice that
$x^{\mu}(1-x)\leq[(1+\frac{1}{\mu})^{\mu}(1+\mu)]^{-1}$, and the equation
holds if and only if $x=\frac{\mu}{1+\mu}=\frac{m}{m+n-1}$. Hence,

$\lim_{n\to\infty}-\frac{d}{dx}L_{n}^{(m)}(x)=nC_{n+m}^{m}x^{m}(1-x)^{n-1}=\left\\{\begin{matrix}\infty,&x=\frac{m}{m+n-1},\\\
0,&otherwise.\end{matrix}\right.$ (86)

Combining it with the fact that $L_{n}^{(m)}(0)=1$ and $L_{n}^{(m)}(1)=0$, we
arrive at Eq. (83). ∎

######  Definition 3.

The function $I_{d}(x,y)$ with $d$ being a positive integer and $x,y\in[0,1]$
is defined as

$I_{d}(x,y):=\frac{(1-x)^{d}(1-y)^{d}}{d}\sum^{d-1}_{j=0}\sum^{d-j-1}_{k=0}(d-j-k)f^{(k)}_{d}x^{k}f^{(j)}_{d}y^{j}.$
(87)

######  Lemma 5.

In the limit $d\rightarrow\infty$, the function $I_{d}(x,y)$ becomes

$\lim_{d\to\infty}I_{d}(x,y)=\left\\{\begin{matrix}1-\frac{x}{1-x}-\frac{y}{1-y},&\frac{y}{1-y}<1-\frac{x}{1-x},\\\
0,&\frac{y}{1-y}\geq 1-\frac{x}{1-x}.\end{matrix}\right.$ (88)

###### Proof.

Let
$r_{d}^{(j)}(x):=\frac{(1-x)^{d}}{d}\sum^{d-j-1}_{k=0}(d-j-k)f^{(k)}_{d}x^{k}$,
and then,

$I_{d}(x,y)=(1-y)^{d}\sum^{d-1}_{j=0}r_{d}^{(j)}(x)f^{(j)}_{d}y^{j}.$ (89)

By definition, the function $r_{d}^{(j)}(x)$ is calculated as

$\displaystyle r_{d}^{(j)}(x)$ $\displaystyle=$
$\displaystyle\frac{d-j}{d}L^{(d-j-1)}_{d}(x)-K^{(d-j-1)}_{d}(x),$ (90)
$\displaystyle=$
$\displaystyle\big{(}\frac{d-j}{d}-\frac{x}{1-x}\big{)}L^{(d-j-1)}_{d}(x).$

In the limit of large $d$, by Lemma 4, we have

$\lim_{d\to\infty}r_{d}^{(j)}(x)=\left\\{\begin{matrix}\frac{d-j}{d}-\frac{x}{1-x},&0\leq
j<j_{x},\\\ 0,&j_{x}<j\leq d-1.\end{matrix}\right.$ (91)

where $j_{x}=\frac{(d-1)(1-2x)}{1-x}$.

Submitting the above expression for $r_{d}^{(j)}(x)$ into Eq. (89), we have
for $d\rightarrow\infty$,

$\displaystyle I_{d}(x,y)$ $\displaystyle=$
$\displaystyle(1-y)^{d}\sum^{d-1}_{j=0}\big{(}\frac{d-j}{d}-\frac{x}{1-x}\big{)}L^{(d-j-1)}_{d}(x)f^{(j)}_{d}y^{j},$
$\displaystyle=$ $\displaystyle(1-y)^{d}\sum^{\left\lfloor
j_{x}\right\rfloor}_{j=0}\big{(}\frac{d-j}{d}-\frac{x}{1-x}\big{)}f^{(j)}_{d}y^{j},$
$\displaystyle=$ $\displaystyle(1-\frac{x}{1-x})L^{(\left\lfloor
j_{x}\right\rfloor)}_{d}(y)-K^{(\left\lfloor j_{x}\right\rfloor)}_{d}(y),$
$\displaystyle=$ $\displaystyle(1-\frac{x}{1-x}-\frac{y}{1-y})L^{(\left\lfloor
j_{x}\right\rfloor)}_{d}(y).$

For $y<\frac{\left\lfloor j_{x}\right\rfloor}{\left\lfloor
j_{x}\right\rfloor+d-1}$, or equivalently, $\frac{y}{1-y}<1-\frac{x}{1-x}$,
$L^{(\left\lfloor j_{x}\right\rfloor)}_{d}(y)=1$, and otherwise,
$L^{(\left\lfloor j_{x}\right\rfloor)}_{d}(y)=0$. Therefore, we arrive at Eq.
(88). ∎

##  Appendix B QUBIT STATE TRANSFORMATIONS UNDER $\bm{\mathrm{MMTP}^{(d)}}$

###  B.1 A motivating example

Here we consider a simplest example, where the target system $S$ is a qubit
with Hamiltonian $H_{S}=E|e\rangle\langle e|$ and initially in state
$[p_{0},1-p_{0}]$, and the memory is a qubit with Hamiltonian $H_{A}=0$ and
initially in state $[\frac{1}{2},\frac{1}{2}]$. The proposal for attaining the
ground population $p^{(2)}$ as in Eq. (9) goes as follows. The initial state
of $SA$ can be expressed as a probability distribution on basis
$\\{|g0\rangle,|g1\rangle,|e0\rangle,|e1\rangle\\}$, i.e.,

$\bm{p}=\frac{1}{2}[p_{0},p_{0};1-p_{0},1-p_{0}].$ (93)

The global MTP consist of two steps. In the first step, full elementary
thermalizations are implemented subsequently between $|g0\rangle$ and each of
$|ej\rangle$ ($j=1,2$), and the state becomes

$\displaystyle\bm{p_{1}}=T^{g0,e1}T^{g0,e0}\bm{p}=\frac{1}{2}[\gamma(1-p_{0}+\gamma),p_{0};1-\gamma,(1-\gamma)(1-p_{0}+\gamma)].$
(94)

Similarly, the second step consists of $T^{g1,e0}$ and $T^{g1,e1}$ in turn,
and output state reads

$\displaystyle\bm{p_{2}}=T^{g1,e1}T^{g1,e0}\bm{p_{1}}=[a_{2}^{(1)},a_{2}^{(2)};b_{1}^{(2)},b_{2}^{(2)}],$
(95)

with

$\displaystyle a_{2}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\gamma(1-p_{0}+\gamma),$ $\displaystyle a_{2}^{(2)}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\gamma[1+(2\gamma-1)(p_{0}-\gamma)],$ $\displaystyle
b_{1}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(1-\gamma)[p_{0}+1-\gamma],$ $\displaystyle
b_{2}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(1-\gamma)[1+(2\gamma-1)(p_{0}-\gamma)].$ (96)

Then, full thermalization is applied to the auxiliary system $A$, leaving the
target state in $\rho^{(2)}$ with $p^{(2)}=a_{2}^{(1)}+a_{2}^{(2)}$, which
equals to Eq. (9).

###  B.2 Simulation of $\beta$-swap by $\bm{\mathrm{MMTP}^{(d)}}$

We propose a proposal for simulating $\beta$-swap via $\mathrm{MMTP}^{(d)}$.
Our results show that, for $d$ large enough, we can use $\mathrm{MMTP}^{(d)}$
to approximate any qubit state transformation enabled by TP. Our proposal
consists of $d^{2}$ elementary full thermalizations and a full thermalization
locally acting on $A$ as a last step. Precisely, the protocol goes as follows.

Initially, the state of $SA$ is prepared as

$\rho=\big{[}p_{0}|g\rangle\langle g|+(1-p_{0})|e\rangle\langle
e|\big{]}\otimes\frac{1}{d}\sum_{l=1}^{d}|l\rangle\langle l|.$ (97)

On basis $\\{|g1\rangle,\dots,|gd\rangle,|e1\rangle,\dots,|ed\rangle\\}$, the
probability distribution vector of the initial state is written as

$\bm{p_{0d}}=[a_{0}^{(1)},\dots,a_{0}^{(d)};b_{1}^{(0)},\dots,b_{d}^{(0)}],$
(98)

with

$a_{0}^{(k)}=\frac{p_{0}}{d},\ b_{j}^{(0)}=\frac{1-p_{0}}{d}.$ (99)

In the $k$th step ($k=1,\dots,d$), full elementary thermalizations are
implemented subsequently between $|gk\rangle$ and each of $|ej\rangle$
($j=1,\dots,d$). Then we have

$\bm{p_{kj}}=\big{\\{}\begin{array}[]{cc}T^{gk,e1}\bm{p_{k-1,d}},&j=1,\\\
T^{gk,ej}\bm{p_{k,j-1}},&2\leq j\leq d,\end{array}$ (100)

where

$\displaystyle\bm{p_{kj}}=[a_{d}^{(1)},\dots,a_{d}^{(k-1)},a_{j}^{(k)},a_{0}^{(k+1)},\cdots,a_{0}^{(d)};b_{1}^{(k)},\dots,b_{j}^{(k)},b_{j+1}^{(k-1)},\dots,b_{d}^{(k-1)}].$
(101)

Direct calculation leads to

$\displaystyle
a_{j}^{(k)}=\gamma(a_{j-1}^{(k)}+b_{j}^{(k-1)}),b_{j}^{(k)}=(1-\gamma)(a_{j-1}^{(k)}+b_{j}^{(k-1)}).$
(102)

Therefore, we have

$b_{j}^{(k)}=\frac{1-\gamma}{\gamma}a_{j}^{(k)},$ (103)

and Eq. (102) becomes

$\displaystyle a_{j}^{(k)}$ $\displaystyle=$ $\displaystyle\gamma
a_{j-1}^{(k)}+(1-\gamma)a_{j}^{(k-1)},$ $\displaystyle a_{0}^{(k)}$
$\displaystyle=$ $\displaystyle\frac{p_{0}}{d},\
a_{j}^{(0)}=\frac{\gamma}{1-\gamma}b_{j}^{(0)}=\frac{\gamma}{1-\gamma}\frac{1-p_{0}}{d}.$
(104)

Now we make the substitution

$a_{j}^{(k)}=\frac{1}{d}\gamma^{j}(1-\gamma)^{k-1}s_{j}^{(k)},$ (105)

and the above problem is equivalent to

$\displaystyle s_{j}^{(k)}$ $\displaystyle=$ $\displaystyle
s_{j-1}^{(k)}+s_{j}^{(k-1)},$ $\displaystyle s_{0}^{(k)}$ $\displaystyle=$
$\displaystyle p_{0}(1-\gamma)^{-k+1},\ s_{j}^{(0)}=(1-p_{0})\gamma^{-j+1}.$
(106)

This is in turn equivalent to

$s_{j}^{(k)}=(1-\gamma)^{-k}[(1-p_{0})\gamma^{-j+1}-(\gamma-
p_{0})\sum_{k^{\prime}=0}^{k-1}f_{j}^{(k^{\prime})}(1-\gamma)^{k^{\prime}}],$
(107)

where $f_{j}^{(k)}:=\sum_{j^{\prime}=1}^{j}f_{j^{\prime}}^{(k-1)},\
f_{j}^{(0)}=1$, or equivalently, $f_{j}^{(k)}=C_{j+k-1}^{k}$. Here
$C_{n}^{m}:=\frac{n!}{m!(n-m)!}$. The proof goes as follows. For $k=0$, Eq.
(107) reduces to $s_{j}^{(0)}=(1-p_{0})\gamma^{-j+1}$, which is equivalent to
Eq. (106). It means that Eq. (107) holds for the case with $k=0$ and
$j=1,\dots,d$. Now assume it holds for the case with $k-1$ and $j=1,\dots,d$,
and we will prove it also holds for $k$. From Eq. (106), we have
$s_{j}^{(k)}-s_{j-1}^{(k)}=s_{j}^{(k-1)}$. Taking the summation over all $j$
gives

$\sum_{j^{\prime}=1}^{d}[s_{j}^{(k)}-s_{j-1}^{(k)}]=\sum_{j^{\prime}=1}^{d}s_{j}^{(k-1)}.$
(108)

It follows that

$\displaystyle s_{j}^{(k)}$ $\displaystyle=$ $\displaystyle
s_{0}^{(k)}+(1-\gamma)^{1-k}\big{[}(1-p_{0})\sum_{j^{\prime}=1}^{j}\gamma^{1-j^{\prime}}-(\gamma-
p_{0})\sum_{k^{\prime}=0}^{k-2}\sum_{j^{\prime}=1}^{j}f_{j^{\prime}}^{(k^{\prime})}(1-\gamma)^{k^{\prime}}\big{]}$
(109) $\displaystyle=$ $\displaystyle
s_{0}^{(k)}+(1-\gamma)^{1-k}\big{[}(1-p_{0})\frac{1-\gamma^{-j}}{1-\gamma^{-1}}-(\gamma-
p_{0})\sum_{k^{\prime}=1}^{k-1}f_{j^{\prime}}^{(k^{\prime})}(1-\gamma)^{k^{\prime}-1}\big{]}$
$\displaystyle=$ $\displaystyle(1-\gamma)^{-k}[(1-p_{0})\gamma^{-j+1}-(\gamma-
p_{0})\sum_{k^{\prime}=0}^{k-1}f_{j}^{(k^{\prime})}(1-\gamma)^{k^{\prime}}].$

In the last line, we used Eq. (106).

After the $d^{2}$ elementary full thermalizations, the state becomes

$\displaystyle\bm{p_{dd}}=[a_{d}^{(1)},\dots,a_{d}^{(d)};b_{1}^{(d)},\dots,b_{d}^{(d)}].$
(110)

It follows from Eq. (107) that

$a_{d}^{(k)}=\frac{1}{d}\frac{\gamma}{1-\gamma}(1-p_{0})-\frac{\gamma-
p_{0}}{d(1-\gamma)}\gamma^{d}\sum_{k^{\prime}=0}^{k-1}f_{d}^{(k^{\prime})}(1-\gamma)^{k^{\prime}}.$
(111)

Then full thermalization is implemented locally on $A$, and the state of the
qubit becomes $[p^{(d)},1-p^{(d)}]$ with
$p^{(d)}\equiv\sum_{k=1}^{d}a_{d}^{(k)}$. Direct calculations lead to

$p^{(d)}=\frac{\gamma}{1-\gamma}(1-p_{0})-\frac{\gamma-
p_{0}}{1-\gamma}I_{d}^{(d-1)}(1-\gamma),$ (112)

where the function $I_{d}^{(d-1)}(x)$ with $x\in(0,1)$ is defined in Appendix
A. It decreases with $d$ and for finite $d$ we have
$I_{d}^{(d-1)}(1-\gamma)\in(\frac{2\gamma-1}{\gamma},\gamma]$. From Lemma 3 in
Appendix A, we have

$p^{(d)}=1-p_{0}e^{-\beta E}-(\gamma-p_{0})\delta_{d}(\gamma),$ (113)

where $\delta_{d}(\gamma)<o[(4\gamma(1-\gamma))^{d}d^{-3/2}]$. Therefore,
$\beta$-swap can be simulated by our protocol with error exponentially
decreasing with the dimension of the assisted memory.

##  Appendix C Three-dimension target system and the comparision between
$\bm{\mathrm{MMTP}^{(d)}}$ and ETP

When the dimension of target system is more than 2, ETP is a strict subset of
TP. Therefore, it is of interest to investigate the comparison between
$\mathrm{MMTP}^{(d)}$ and ETP, especially when $d$ is not large.

Consider a three-dimensional system with Hamiltonian $H_{S}=0|g\rangle\langle
g|+E(|e_{1}\rangle\langle e_{1}|+|e_{2}\rangle\langle e_{2}|)$ and initially
in the ground state $|g\rangle$. The occupation probability vector of the
initial state reads $[1,0,0]$.

The boundaries of states that can be obtained from $|g\rangle$ via TP, ETP, as
well as MTP, are plotted in Fig. 4 . Now we show that with the assistance of a
qubit memory, MTP can achieve state transformation that cannot be achieved by
ETP.

The vertices can be obtained as
Vertex
$A_{s}$($s=1,2$):$T_{M}T^{g2,e_{s}2}T^{g2,e_{s}1}T^{g1,e_{s}2}T^{g1,e_{s}1}$.
Vertex $B_{s}$
($s=1,2$):$T_{M}T^{g2,e_{\bar{s}}2}T^{g2,e_{\bar{s}}1}T^{g1,e_{\bar{s}}2}T^{g1,e_{\bar{s}}1}T^{g2,e_{s}2}T^{g2,e_{s}1}T^{g1,e_{s}2}T^{g1,e_{s}1}$.
Other states on the boundary can be obtained by elementary thermalizations
from the vertices. The states inside can be obtained by thermalization from
boundary states.

Figure 4: Comparison of the sets of state that can be reached via TP, ETP,
MTP,and $\bm{\mathrm{MMTP}^{(2)}}$ from a qutrit state. The Hamiltonian of the
qutrit reads $H_{S}=0|g\rangle\langle g|+E(|e_{1}\rangle\langle
e_{1}|+|e_{2}\rangle\langle e_{2}|)$, and the qutrit is initially in its
ground state $|g\rangle$. This example shows that via $\mathrm{MMTP}^{(2)}$
one can realize state transformations that cannot be achieved by ETP.

##  Appendix D Details in the memory-assisted protocol for work extraction

The protocol is summarized in the following scheme

$\displaystyle[0,0,1,0]\otimes[\overbrace{\frac{1}{d},...,\frac{1}{d}}^{d\
times}]$
$\displaystyle=[{0,...,0};0,...,0;\frac{1}{d},...,\frac{1}{d};0,...,0]$
$\displaystyle\stackrel{{\scriptstyle\textbf{(\@slowromancap
i@)}}}{{\longmapsto}}[0,...,0;a_{d}^{(1)},...,a_{d}^{(d)};b_{1}^{(d)},...,b_{d}^{(d)};0,...,0]$
$\displaystyle\stackrel{{\scriptstyle\textbf{(\@slowromancap
ii@)}}}{{\longmapsto}}[0,...,0;a_{d}^{(1)},...,a_{d}^{(d)};\epsilon_{1},...,\epsilon_{d};c_{1}^{(d)},...,c_{d}^{(d)}].$

Step i@ consists of $d$ subroutines. In the $j$th subroutine, full elementary
thermalizations are consequently implemented between $|e0j\rangle$ and each of
$|g1k\rangle$ with $k=1,\dots,d$. This is equivalent to simulated $\beta$-swap
between $|e0\rangle_{SW}$ and $|g1\rangle_{SW}$. By putting Eqs. (103), (105)
and (107) together, we obtain

$b_{j}^{(d)}=\frac{1}{d}[1-\gamma_{\delta}^{j}\sum^{d-1}_{k=0}f^{(k)}_{j}(1-\gamma_{\delta})^{k}].$
(114)

Further, because

$\displaystyle\sum^{j-1}_{j^{\prime}=0}f^{(k)}_{j^{\prime}+1}\gamma_{\delta}^{j^{\prime}}-f^{(k+1)}_{j}\gamma_{\delta}^{j}$
$\displaystyle=$ $\displaystyle
f^{(k+1)}_{1}+\sum^{j-1}_{j^{\prime}=1}(f^{(k+1)}_{j^{\prime}+1}-f^{(k+1)}_{j^{\prime}})\gamma_{\delta}^{j^{\prime}}-f^{(k+1)}_{j}\gamma_{\delta}^{j}=\sum^{j-1}_{j^{\prime}=0}f^{(k+1)}_{j^{\prime}+1}\gamma_{\delta}^{j^{\prime}}-\sum^{j-1}_{j^{\prime}=0}f^{(k+1)}_{j^{\prime}+1}\gamma_{\delta}^{j^{\prime}+1}$
(115) $\displaystyle=$
$\displaystyle(1-\gamma_{\delta})\sum^{j-1}_{j^{\prime}=0}f^{(k+1)}_{j^{\prime}+1}\gamma_{\delta}^{j^{\prime}},$

we arrive at

$\displaystyle b^{(d)}_{j}$ $\displaystyle=$
$\displaystyle\frac{1}{d}[(1-\gamma_{\delta}^{j})-\gamma_{\delta}^{j}\sum^{d-1}_{k=1}f^{(k)}_{j}(1-\gamma_{\delta})^{k}]=\frac{1-\gamma_{\delta}}{d}[\sum^{j-1}_{j^{\prime}=0}f^{(0)}_{j^{\prime}+1}\gamma_{\delta}^{j^{\prime}}-f^{(1)}_{j}\gamma_{\delta}^{j}-\gamma_{\delta}^{j}\sum^{d-1}_{k=2}f^{(k)}_{j}(1-\gamma_{\delta})^{k-1}]$
(116) $\displaystyle=$
$\displaystyle\frac{(1-\gamma_{\delta})^{2}}{d}[\sum^{j-1}_{j^{\prime}=0}f^{(1)}_{j^{\prime}+1}\gamma_{\delta}^{j^{\prime}}-f^{(2)}_{j}\gamma_{\delta}^{j}-\gamma_{\delta}^{j}\sum^{d-1}_{k=3}f^{(k)}_{j}(1-\gamma_{\delta})^{k-2}]=...=\frac{(1-\gamma_{\delta})^{d-1}}{d}[\sum^{j-1}_{j^{\prime}=0}f^{(d-2)}_{j^{\prime}+1}\gamma_{\delta}^{j^{\prime}}-f^{(d-1)}_{j}\gamma_{\delta}^{j}]$
$\displaystyle=$
$\displaystyle\frac{(1-\gamma_{\delta})^{d}}{d}\sum^{j-1}_{j^{\prime}=0}f^{(d-1)}_{j^{\prime}+1}\gamma_{\delta}^{j^{\prime}}.$

Notice that $b^{(d)}_{j}$ is increasing in $j$.

Step ii@ also consists of $d$ subroutines. In the following, we will label the
action of Step (ii@) as $\hat{\mathcal{T}}$, and thus,

$[b_{1}^{(d)},...,b_{d}^{(d)};0,...,0]\stackrel{{\scriptstyle\hat{\mathcal{T}}}}{{\longmapsto}}[\epsilon_{1},...,\epsilon_{d};c_{1}^{(d)},...,c_{d}^{(d)}].$
(117)

Each subroutine of $\hat{\mathcal{T}}$ is realized as
$\hat{\mathcal{T}}^{(k)}\equiv T^{e0k,e1d}\dots T^{e0k,e11}$. In order to
decide the order of $\hat{\mathcal{T}}^{(k)}$ in $\hat{\mathcal{T}}$, we first
consider the following problem.

Let $H_{eff}=0|\alpha\rangle\langle\alpha|+W\sum_{j=1}^{d}|\beta
j\rangle\langle\beta j|$ be the effective Hamiltonian on a $(d+1)$-dimensional
subspace spanned by $\\{|\alpha\rangle;|\beta 1\rangle,\dots,|\beta
d\rangle\\}$, and $\mathcal{T}_{eff}=T^{\alpha,\beta d}\dots T^{\alpha,\beta
2}T^{\alpha,\beta 1}$ be a Markovian thermal process acting on this subspace.
We define the parameters $d_{k}$ and $c_{j}^{(k)}$ ($j,k=1,\dots,d$) as

$\displaystyle[1;0,\dots,0]$
$\displaystyle\stackrel{{\scriptstyle\mathcal{T}_{eff}}}{{\longmapsto}}[d_{1};c_{1}^{(1)},\dots,c_{d}^{(1)}],$
(118) $\displaystyle[0;c_{1}^{(k-1)},\dots,c_{d}^{(k-1)}]$
$\displaystyle\stackrel{{\scriptstyle\mathcal{T}_{eff}}}{{\longmapsto}}[d_{k};c_{1}^{(k)},\dots,c_{d}^{(k)}],\
k=2,\dots,d.$ (119)

Then $d_{k}$ is calculated as follows. Because

$\displaystyle[1;\overbrace{0,...,0}^{d\ times}]\stackrel{{\scriptstyle
T^{\alpha,\beta 1}}}{{\longmapsto}}[\gamma_{{}_{\tiny W}};1-\gamma_{{}_{\tiny
W}},\overbrace{0,...,0}^{d-1\ times}]\stackrel{{\scriptstyle T^{\alpha,\beta
2}}}{{\longmapsto}}[\gamma_{{}_{\tiny W}}^{2};1-\gamma_{{}_{\tiny
W}},\gamma_{{}_{\tiny W}}(1-\gamma_{{}_{\tiny W}}),\overbrace{0,...,0}^{d-2\
times}]\ \stackrel{{\scriptstyle}}{{\longmapsto}}\dots$
$\displaystyle\stackrel{{\scriptstyle T^{\alpha,\beta
d}}}{{\longmapsto}}[\gamma_{{}_{\tiny W}}^{d};1-\gamma_{{}_{\tiny
W}},\gamma_{{}_{\tiny W}}(1-\gamma_{{}_{\tiny W}}),...,(1-\gamma_{{}_{\tiny
W}})\gamma_{{}_{\tiny W}}^{d-1}],$

from Eq. (118), we have $d_{1}=\gamma_{{}_{\tiny W}}^{d}$,
$c^{(1)}_{j}=(1-\gamma_{{}_{\tiny W}})\gamma_{{}_{\tiny W}}^{j-1}$. Similarly,
from Eq. (119), we have

$\displaystyle d_{k+1}$ $\displaystyle=$ $\displaystyle\gamma_{{}_{\tiny
W}}^{d}c_{1}^{(k)}+\gamma_{{}_{\tiny
W}}^{d-1}c_{2}^{(k)}+...+\gamma_{{}_{\tiny W}}c_{d}^{(k)},$ (120)
$\displaystyle c^{(k+1)}_{j}$ $\displaystyle=$
$\displaystyle(1-\gamma_{{}_{\tiny W}})(\gamma^{j-1}_{{}_{\tiny
W}}c^{(k)}_{1}+\gamma^{j-2}_{{}_{\tiny W}}c^{(k)}_{2}+...+c^{(k)}_{j}).$ (121)

Now let $c^{(k)}_{j}=(1-\gamma_{{}_{\tiny W}})^{k}\gamma_{{}_{\tiny
W}}^{j-1}\chi_{j}^{k}$, and we have
$\chi_{j}^{k+1}=\sum^{j}_{j^{\prime}=1}\chi_{j^{\prime}}^{k}$, and
$c^{(1)}_{j}=(1-\gamma_{{}_{\tiny W}})\gamma_{{}_{\tiny
W}}^{j-1}=(1-\gamma_{{}_{\tiny W}})\gamma_{{}_{\tiny W}}^{j-1}\chi_{j}^{1}$.
Hence, $\chi_{j}^{k}=f_{j}^{(k-1)}$,$c^{(k)}_{j}=(1-\gamma_{{}_{\tiny
W}})^{k}\gamma_{{}_{\tiny W}}^{j-1}f_{j}^{(k-1)}$. In turn, we get
$d_{k+1}=\sum^{d}_{j=1}\gamma_{{}_{\tiny
W}}^{d+1-j}c^{(k)}_{j}=\gamma_{{}_{\tiny W}}^{d}(1-\gamma_{{}_{\tiny
W}})^{k}f^{k}_{d}$. Therefore, we arrive at

$d_{k}=\gamma_{{}_{\tiny W}}^{d}(1-\gamma_{{}_{\tiny W}})^{k-1}f^{(k-1)}_{d}.$
(122)

Now consider a $2d$-subspace spanned by $\\{|\alpha 1\rangle,\dots,|\alpha
d\rangle;|\beta 1\rangle,\dots,|\beta d\rangle\\}$, and let
$\vec{x}=[x_{1},x_{2},...,x_{d};0,...,0]$ be a probability distribution in
this subspace. (Notice that we do $not$ require $\sum_{j=1}^{d}x_{j}=1$.) Here
$|\alpha j\rangle$ and $|\beta k\rangle$ are energy eigenstates with energy
$E_{\alpha}$ and $E_{\beta}$, respectively, which satisfy
$E_{\beta}-E_{\alpha}=W$. We label
$\mathcal{T}=\mathcal{T}^{(d)}\dots\mathcal{T}^{(1)}$, where
$\mathcal{T}^{(k)}=T^{\alpha k,\beta d}\dots T^{\alpha k,\beta 1}$, and
calculate the output probability distribution after the action of
$\mathcal{T}$ on $\vec{x}$.

The operator $\mathcal{T}^{(1)}$ acts on the subspace spanned by $\\{|\alpha
1\rangle;|\beta 1\rangle,\dots,|\beta d\rangle\\}$. In this subspace, the
action of $\mathcal{T}^{(1)}$ is effectively the action of
$\mathcal{T}_{eff}$, so the state transforms as

$x_{1}[1;0,...,0]\stackrel{{\scriptstyle\mathcal{T}^{(1)}}}{{\longmapsto}}x_{1}[d_{1};c_{1}^{(1)},...,c_{d}^{(1)}].$
(123)

Similarly, $\mathcal{T}^{(2)}$ acts on the subspace spanned by $\\{|\alpha
2\rangle;|\beta 1\rangle,\dots,|\beta d\rangle\\}$, and also effectively
equivalent to $\mathcal{T}_{eff}$. The input state of $\mathcal{T}^{(2)}$ in
this subspace reads

$[x_{2};x_{1}c_{1}^{(1)},...,x_{1}c_{d}^{(1)}]=x_{2}[1;0,...,0]+x_{1}[0;c_{1}^{(1)},...,c_{d}^{(1)}].$
(124)

It follows that

$[x_{2};x_{1}c_{1}^{(1)},...,x_{1}c_{d}^{(1)}]\stackrel{{\scriptstyle\mathcal{T}^{(2)}}}{{\longmapsto}}x_{2}[d_{1};c_{1}^{(1)},...,c_{d}^{(1)}]+x_{1}[d_{2};c_{1}^{(2)},...,c_{d}^{(2)}].$
(125)

Analogously, $\mathcal{T}^{(k)}$ acts on the subspace spanned by $\\{|\alpha
k\rangle;|\beta 1\rangle,\dots,|\beta d\rangle\\}$. The input state of
$\mathcal{T}^{(k)}$ in this subspace reads

$x_{k}[1;0,...,0]+x_{k-1}[0;c_{1}^{(1)},...,c_{d}^{(1)}]+\dots+x_{1}[0;c_{1}^{(k-1)},...,c_{d}^{(k-1)}].$
(126)

After the action of $\mathcal{T}^{(k)}$, the above state becomes

$x_{k}[d_{1};c_{1}^{(1)},...,c_{d}^{(1)}]+x_{k-1}[d_{1};c_{1}^{(1)},...,c_{d}^{(1)}]+\dots+x_{1}[d_{k};c_{1}^{(k)},...,c_{d}^{(k)}].$
(127)

Therefore, the population on $|\alpha k\rangle$ after the action of
$\mathcal{T}^{(k)}$ reads

$\epsilon_{k}=d_{1}x_{k}+d_{2}x_{k-1}+...+d_{k}x_{1},$ (128)

which preserves during the actions of subsequent operations
$\mathcal{T}^{(k^{\prime})}$ with $k^{\prime}>k$. Finally, the total
population on $|\alpha 1\rangle,\dots,|\alpha d\rangle$ reads

$\displaystyle\epsilon^{(d)}$ $\displaystyle=$
$\displaystyle\sum^{d}_{k=1}\epsilon_{k}=d_{1}x_{d}+...+\sum_{k=1}^{d+1-i}d_{k}x_{i}+...+\sum_{k=1}^{d}d_{k}x_{1}=\sum_{i=1}^{d}\left(\sum_{k=1}^{d+1-i}d_{k}\right)x_{i}.$
(129)

Because $d_{k}\geq 0$, $\epsilon^{(d)}$ is minimized when $x_{i}$ is
increasing in $i$. Therefore, we choose $b_{j}^{(d)}=x_{j}$ in the original
problem as in Eq. (117), and the optimal order of $\hat{\mathcal{T}}^{(k)}$ is
$\hat{\mathcal{T}}=\hat{\mathcal{T}}^{(d)}\dots\hat{\mathcal{T}}^{(2)}\hat{\mathcal{T}}^{(1)}$.

Finally, the error of work extraction under our memory-assisted protocol is
calculated as

$\displaystyle\epsilon^{(d)}$ $\displaystyle=$
$\displaystyle\sum^{d}_{j=1}\sum^{d-j+1}_{k=1}d_{k}b^{(d)}_{j}=\sum^{d}_{j=1}\sum^{d-j+1}_{k=1}\gamma_{{}_{\tiny
W}}^{d}(1-\gamma_{{}_{\tiny
W}})^{k-1}f^{(k-1)}_{d}\frac{(1-\gamma_{\delta})^{d}}{d}\sum^{j-1}_{j^{\prime}=0}f^{(d-1)}_{j^{\prime}+1}\gamma_{\delta}^{j^{\prime}}$
(130) $\displaystyle=$
$\displaystyle\frac{(\gamma_{{}_{W}}(1-\gamma_{\delta}))^{d}}{d}\sum^{d}_{j=1}\sum^{d-j}_{k=0}(1-\gamma_{{}_{\tiny
W}})^{k}f^{(k)}_{d}\sum^{j-1}_{j^{\prime}=0}f^{(d-1)}_{j^{\prime}+1}\gamma_{\delta}^{j^{\prime}}.$

where
$f^{(d-1)}_{j^{\prime}+1}=C^{d-1}_{j^{\prime}+d-1}=C^{j^{\prime}}_{j^{\prime}+d-1}=f^{(j^{\prime})}_{d}$,
and we have the following summation exchange

$\displaystyle\sum^{d}_{j=1}\sum^{d-j}_{k=0}\sum^{j-1}_{j^{\prime}=0}\to\sum^{d-1}_{j^{\prime}=0}\sum^{d}_{j=j^{\prime}+1}\sum^{d-j}_{k=0}\to\sum^{d-1}_{j^{\prime}=0}\sum^{d-j^{\prime}-1}_{k=0}\sum^{d-k}_{j=j^{\prime}+1}$
(131)

All the terms in the summation are independent of $j$, so the last summation
$\sum^{d-k}_{j=j^{\prime}+1}$ equivalent to multiplying each term by
$(d-k-j^{\prime})$. Hence,

$\displaystyle\epsilon^{(d)}$ $\displaystyle=$
$\displaystyle\frac{(1-\gamma_{\delta})^{d}}{d}\sum^{d-1}_{j^{\prime}=0}f^{(j^{\prime})}_{d}\gamma_{\delta}^{j^{\prime}}\gamma_{{}_{\tiny
W}}^{d}\sum^{d-j^{\prime}-1}_{k=0}(d-k-j^{\prime})f^{(k)}_{d}(1-\gamma_{{}_{\tiny
W}})^{k}$ (132) $\displaystyle=$
$\displaystyle\frac{(1-\gamma_{\delta})^{d}}{d}\sum^{d-1}_{j^{\prime}=0}(d-j^{\prime})f^{(j^{\prime})}_{d}(\gamma_{\delta}\gamma_{{}_{\tiny
W}})^{j^{\prime}}\frac{\gamma^{n-j^{\prime}}_{{}_{\tiny
W}}}{d-j^{\prime}}\sum^{d-j^{\prime}-1}_{k=0}(d-j^{\prime}-k)f^{(k)}_{d}(1-\gamma_{{}_{\tiny
W}})^{k}=I_{d}(\gamma_{\delta},1-\gamma_{{}_{\tiny W}}).$

Here $I_{d}(x,y)$ is the function defined in Eq. (88). This completes the
proof of Theorem 4.

## References

  * Maruyama  _et al._ [2009] K. Maruyama, F. Nori, and V. Vedral, Colloquium: The physics of maxwell’s demon and information, Rev. Mod. Phys. 81, 1 (2009).
  * Beyer  _et al._ [2019] K. Beyer, K. Luoma, and W. T. Strunz, Steering heat engines: A truly quantum maxwell demon, Phys. Rev. Lett. 123, 250606 (2019).
  * Park  _et al._ [2013] J. J. Park, K.-H. Kim, T. Sagawa, and S. W. Kim, Heat engine driven by purely quantum information, Phys. Rev. Lett. 111, 230402 (2013).
  * 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).
  * Cavaliere  _et al._ [2022] F. Cavaliere, M. Carrega, G. De Filippis, V. Cataudella, G. Benenti, and M. Sassetti, Dynamical heat engines with non-markovian reservoirs, Phys. Rev. Res. 4, 033233 (2022).
  * Manzano  _et al._ [2021] G. Manzano, D. Subero, O. Maillet, R. Fazio, J. P. Pekola, and E. Roldán, Thermodynamics of gambling demons, Phys. Rev. Lett. 126, 080603 (2021).
  * Naghiloo  _et al._ [2018] M. Naghiloo, J. J. Alonso, A. Romito, E. Lutz, and K. W. Murch, Information gain and loss for a quantum maxwell’s demon, Phys. Rev. Lett. 121, 030604 (2018).
  * Admon  _et al._ [2018] T. Admon, S. Rahav, and Y. Roichman, Experimental realization of an information machine with tunable temporal correlations, Phys. Rev. Lett. 121, 180601 (2018).
  * Hu  _et al._ [2020] C.-K. Hu, A. C. Santos, J.-M. Cui, Y.-F. Huang, D. O. Soares-Pinto, M. S. Sarandy, C.-F. Li, and G.-C. Guo, Quantum thermodynamics in adiabatic open systems and its trapped-ion experimental realization, npj Quantum Inf. 6, 73 (2020).
  * Boykin  _et al._ [2002] P. O. Boykin, T. Mor, V. Roychowdhury, F. Vatan, and R. Vrijen, Algorithmic cooling and scalable nmr quantum computers, Proceedings of the National Academy of Sciences 99, 3388 (2002).
  * Schulman  _et al._ [2005] L. J. Schulman, T. Mor, and Y. Weinstein, Physical limits of heat-bath algorithmic cooling, Phys. Rev. Lett. 94, 120501 (2005).
  * Raeisi and Mosca [2015] S. Raeisi and M. Mosca, Asymptotic bound for heat-bath algorithmic cooling, Phys. Rev. Lett. 114, 100404 (2015).
  * Clivaz  _et al._ [2019] F. Clivaz, R. Silva, G. Haack, J. B. Brask, N. Brunner, and M. Huber, Unifying paradigms of quantum refrigeration: A universal and attainable bound on cooling, Phys. Rev. Lett. 123, 170605 (2019).
  * Brandão  _et al._ [2013] F. G. S. L. Brandão, M. Horodecki, J. Oppenheim, J. M. Renes, and R. W. Spekkens, Resource theory of quantum states out of thermal equilibrium, Phys. Rev. Lett. 111, 250404 (2013).
  * Lostaglio [2019] M. Lostaglio, An introductory review of the resource theory approach to thermodynamics, Reports on Progress in Physics 82, 114001 (2019).
  * Horodecki and Oppenheim [2013] M. Horodecki and J. Oppenheim, Fundamental limitations for quantum and nanoscale thermodynamics, Nature Communications 4, 2059 (2013).
  * Francesco and Lorenza [2014] T. Francesco and V. Lorenza, Quantum resources for purification and cooling: fundamental limits and opportunities, Scientific Reports 4, 5192 (2014).
  * Allahverdyan  _et al._ [2011] A. E. Allahverdyan, K. V. Hovhannisyan, D. Janzing, and G. Mahler, Thermodynamic limits of dynamic cooling, Phys. Rev. E 84, 041109 (2011).
  * Wilming and Gallego [2017] H. Wilming and R. Gallego, Third law of thermodynamics as a single inequality, Phys. Rev. X 7, 041033 (2017).
  * Taranto  _et al._ [2023] P. Taranto, F. Bakhshinezhad, A. Bluhm, R. Silva, N. Friis, M. P. Lock, G. Vitagliano, F. C. Binder, T. Debarba, E. Schwarzhans, F. Clivaz, and M. Huber, Landauer versus nernst: What is the true cost of cooling a quantum system?, PRX Quantum 4, 010332 (2023).
  * Alhambra  _et al._ [2019] Á. M. Alhambra, M. Lostaglio, and C. Perry, Heat-Bath Algorithmic Cooling with optimal thermalization strategies, Quantum 3, 188 (2019).
  * Korzekwa and Lostaglio [2022] K. Korzekwa and M. Lostaglio, Optimizing thermalization, Phys. Rev. Lett. 129, 040602 (2022).
  * Lostaglio and Korzekwa [2022] M. Lostaglio and K. Korzekwa, Continuous thermomajorization and a complete set of laws for markovian thermal processes, Phys. Rev. A 106, 012426 (2022).
  * Ćwikliński  _et al._ [2015] P. Ćwikliński, M. Studziński, M. Horodecki, and J. Oppenheim, Limitations on the evolution of quantum coherences: Towards fully quantum second laws of thermodynamics, Phys. Rev. Lett. 115, 210403 (2015).
  * Lostaglio  _et al._ [2015] M. Lostaglio, K. Korzekwa, D. Jennings, and T. Rudolph, Quantum coherence, time-translation symmetry, and thermodynamics, Phys. Rev. X 5, 021001 (2015).
  * Lostaglio  _et al._ [2018] M. Lostaglio, Á. M. Alhambra, and C. Perry, Elementary Thermal Operations, Quantum 2, 52 (2018).
  * Korzekwa [2016] K. Korzekwa,  _Cohernce, thermodynamics and uncertainty relations_ , Ph.D. thesis, Imperial College London (2016).
  * Brandão  _et al._ [2015] F. Brandão, M. Horodecki, N. Ng, J. Oppenheim, and S. Wehner, The second laws of quantum thermodynamics, Proceedings of the National Academy of Sciences 112, 3275 (2015).
  * Müller [2018] M. P. Müller, Correlating thermal machines and the second law at the nanoscale, Phys. Rev. X 8, 041051 (2018).
  * Lipka-Bartosik and Skrzypczyk [2021] P. Lipka-Bartosik and P. Skrzypczyk, All states are universal catalysts in quantum thermodynamics, Phys. Rev. X 11, 011061 (2021).
  * Takagi and Shiraishi [2022] R. Takagi and N. Shiraishi, Correlation in catalysts enables arbitrary manipulation of quantum coherence, Phys. Rev. Lett. 128, 240501 (2022).
  * Henao and Uzdin [2021] I. Henao and R. Uzdin, Catalytic transformations with finite-size environments: applications to cooling and thermometry, Quantum 5, 547 (2021).
  * Ding  _et al._ [2021] F. Ding, X. Hu, and H. Fan, Amplifying asymmetry with correlating catalysts, Phys. Rev. A 103, 022403 (2021).
  * Patryk  _et al._ [2023] L.-B. Patryk, W. Henrik, and N. H. Y. Ng, Catalysis in quantum information theory, arXiv:2306.00798 (2023).
  * Son and Ng [2022] J. Son and N. H. Y. Ng, Catalysis in action via elementary thermal operations, arXiv:2209.15213 (2022).
  * Son and Ng [2023] J. Son and N. H. Y. Ng, A hierarchy of thermal processes collapses under catalysis, arXiv:2303.13020 (2023).
  * Czartowski  _et al._ [2023] J. Czartowski, A. d. O. Junior, and K. Korzekwa, Thermal recall: Memory-assisted markovian thermal processes, arXiv:2303.12840 (2023).