# Thermal recall: Memory-assisted Markovian thermal processes

Jakub Czartowski∗ Doctoral School of Exact and Natural Sciences, Jagiellonian
University, 30-348 Kraków, Poland  Faculty of Physics, Astronomy and Applied
Computer Science, Jagiellonian University, 30-348 Kraków, Poland.  A. de
Oliveira Junior∗ Kamil Korzekwa  Faculty of Physics, Astronomy and Applied
Computer Science, Jagiellonian University, 30-348 Kraków, Poland.

###### Abstract

We develop a resource-theoretic framework that allows one to bridge the gap
between two approaches to quantum thermodynamics based on Markovian thermal
processes (which model memoryless dynamics) and thermal operations (which
model arbitrarily non-Markovian dynamics). Our approach is built on the notion
of memory-assisted Markovian thermal processes, where memoryless thermodynamic
processes are promoted to non-Markovianity by explicitly modelling ancillary
memory systems initialised in thermal equilibrium states. Within this setting,
we propose a family of protocols composed of sequences of elementary two-level
thermalisations that approximate all transitions between energy-incoherent
states accessible via thermal operations. We prove that, as the size of the
memory increases, these approximations become arbitrarily good for all
transitions in the infinite temperature limit, and for a subset of transitions
in the finite temperature regime. Furthermore, we present solid numerical
evidence for the convergence of our protocol to any transition at finite
temperatures. We also explain how our framework can be used to quantify the
role played by memory effects in thermodynamic protocols such as work
extraction. Finally, our results show that elementary control over two energy
levels at a given time is sufficient to generate all energy-incoherent
transitions accessible via thermal operations if one allows for ancillary
thermal systems.

##  I Introduction

00footnotetext: ∗ These authors contributed equally to this work.

Information has become ubiquitous in thermodynamics. It all started with
Maxwell’s seminal inquiry [1]: _what would happen if we had knowledge of a
system’s state?_ The ramifications of this hypothesis led to potential
violations of the second law of thermodynamics and a century-long puzzle [2,
3]. Ultimately, it was found that thermodynamics imposes physical restrictions
on information processing [4, 5], resulting in the development of frameworks
devoted to incorporating information into thermodynamics [6, 7, 8, 9, 10, 11].
A crucial concept at the intersection between these two fields is _memory_ , a
thermodynamic resource for storing, processing, and erasing information. In
particular, memory effects can bring numerous advantages, including enhanced
cooling [12], generation of entanglement [13, 14] or improved performance of
heat engines and refrigerators [15, 16, 17]. However, realistic quantum
mechanical systems are open and governed by non-unitary time evolution, which
encompasses the irreversible phenomena such as energy dissipation, relaxation
to thermal equilibrium or stationary non-equilibrium states, and the decay of
correlations [18, 19]. Hence, assumptions like weak coupling, large bath size,
and fast decaying correlations are commonly made in modelling such systems,
thus neglecting memory effects. This raises the question of how memoryless
processes get modified when system-bath memory effects become non-negligible,
i.e., how to assess and quantify the role of memory in thermodynamic processes
[20].

The resource theory of thermodynamics [21, 22, 23, 24] is a relatively recent
framework allowing one to address foundational questions in thermodynamics. By
relying on the notion of thermal operations [21, 22], a set of transformations
that can be carried out without an external source of work or coherence, it
offers a complete set of laws for characterising general state transformations
under thermodynamic constraints. The downsides of this formalism are twofold.
Firstly, it focuses only on snapshots of the evolution, making it hard to
discuss how the processes are realised in time. Secondly, it may require
precise control over the system and the bath. The first problem was addressed
by developing a hybrid framework that reconciles resource theory and master
equation approaches [25, 26], where the concept of a Markovian thermal process
was introduced. This new set of operations refines the notion of thermal
operations by encoding all relevant constraints of a Markovian evolution. The
second problem was partially addressed in Ref. [27] by introducing the concept
of elementary thermal operations, i.e., a subset of transformations that can
be decomposed into a series of thermal operations, each acting only on two
energy levels of the system. Such decompositions offer a method to bypass the
need for a complete control over interactions between the system and the
environment, and the approach was recently generalised to also include
catalytic transformations [28]. While elementary operations require only a
limited control, they still rely on non-Markovian effects, and so the question
of quantifying memory effects in the resource theory of thermodynamics remains
open.

Figure 1: Memory-assisted Markovian thermal processes. Schematic
representation of the general setting. (a) Initially, the main system (large
blue circle) is coupled to a heat bath at inverse temperature $\beta$ (small
red circles) and their interaction is Markovian, so that the bath is in
thermal equilibrium at each moment in time. (b)-(e) Then, the control is
extended to parts of the environment (small blue circles with blue background)
that do not instantaneously thermalise to equilibrium after interactions with
the system, and can thus lead to non-Markovian dynamics of the main system.

In this work, we make a step towards bridging the gap between thermal
operations and Markovian thermal processes for energy-incoherent states by
introducing and investigating _memory-assisted Markovian thermal processes_
(MeMTPs). These are defined by extending the Markovian thermal processes
framework [25, 26] with ancillary memory systems, allowing one to interpolate
between memoryless dynamics and the one with full control. More specifically,
we demonstrate that energy-incoherent states achievable from a given initial
state via thermal operations can be approached arbitrarily well by repeatedly
interacting the main system with a memory that is initialised in a thermal
equilibrium state (and therefore is thermodynamically resourceless) via an
algorithmic procedure composed of Markovian thermal processes. Physically,
this can be seen as a partial control over the bath degrees of freedom, where
the bath can be thought of as a large, discrete, collection of smaller thermal
units, and one can control the interactions of the main system with a small
number of thermal subsystems (see Fig. 1).

Following this idea, we introduce a family of memory-extended Markovian
thermodynamic protocols that require minimal control and are valid for any
temperature regime. In the infinite temperature limit, we prove that our
protocol can arbitrarily well simulate any state transition that can be
achieved via thermal operations. More precisely, our first main result states
that when memory grows and the number of interactions goes to infinity, the
full set of states achievable by thermal operations can be reached by MeMTPs.
We also provide analytic expressions for the convergence rates, which scale
either polynomially with the number of degrees of freedom or exponentially
with the number of memory subsystems. Moreover, based on strong numerical
evidence, we propose a conjecture for a more accurate approximation of
arbitrary state transformation (i.e., converging faster with the growing size
of the memory) through sequences of truncated versions of our protocol. Our
second main result extends these considerations to the finite-temperature
regime. Here, we first show analytic convergence of our MeMTP protocol to a
particular subset of state transformations that can be achieved by thermal
operations. These include all the so-called $\beta$-swaps, as well as
$\beta$-cycles, which form a thermodynamic equivalent of cyclic permutations.
Then, based on numerical simulations, we conjecture that actually an arbitrary
state reachable via thermal operations can be obtained using a proper sequence
of truncated protocols.

With these results at hand, we then proceed to discussing their applicability.
First, we explain how to assess the role played by memory in the performance
of thermodynamic protocols by investigating work extraction in the
intermediate regime of limited memory. We thus interpolate between the two
extremes of no memory and complete control, and quantify environmental memory
effects with the amount of extractable work from a given non-equilibrium
state. Second, we consider the task of cooling a two-level system using a two-
dimensional memory characterised by a non-trivial Hamiltonian. This example
represents a minimal model requiring the manipulation of two two-level
systems. Various experimental proposals are available across distinct
platforms suitable for realising this specific model. Such platforms encompass
quantum dots [29, 30], superconducting circuits [31, 32], and atom-cavity
systems [33]. We then clarify that all transitions achievable via thermal
operations can be performed using a subset of thermal operations that only
affect two energy levels at the same time. This may seem to contradict the
results of Refs. [27, 34], where it was proven that thermal operations
constrained to just two energy levels of the system are not able to generate
all thermodynamically allowed transitions. We resolve this apparent
contradiction by noting that in our case we require the control over two
levels of the joint system-memory state, and not just the system state.
Finally, we discuss the behaviour of free energy during a non-Markovian
evolution, explaining the role of the memory as a free energy storage.

The paper is structured as follows. First, in Sec. II we recall the frameworks
of thermal operations and Markovian thermal processes. Next, in Sec. III, we
introduce the central notion of this paper, the memory-assisted Markovian
thermal processes, and then describe the protocol that employs thermal memory
states to approximate non-Markovian thermodynamic state transitions with
Markovian thermal processes. We then explain how this approximation
convergences to the full set of transitions achievable via thermal operations
as the size of the memory grows. Sec. IV contains discussion and application
of our results. Finally, in Sec. V, we conclude and provide outlook for future
research.

##  II Framework

In this work our aim is to investigate how memory effects affect
thermodynamics of a finite-dimensional quantum system coupled to a heat bath
at inverse temperature $\beta=1/k_{B}T$ with $k_{B}$ denoting the Boltzmann
constant. The investigated system of dimension $d$ is described by a
Hamiltonian $H=\sum_{i}E_{i}\outerproduct{E_{i}}{E_{i}}$ and is prepared in an
initial state $\rho$. The thermal environment, with a Hamiltonian $H_{E}$, is
assumed to be in a thermal equilibrium state,

$\gamma_{E}=\frac{e^{-\beta H_{E}}}{\tr(e^{-\beta H_{E}})},$ (1)

and the system together with the heat bath start initially in an uncorrelated
state $\rho\otimes\gamma_{E}$. We will limit our considerations to _energy-
incoherent_ states $\rho$, i.e., the ones that commute with the Hamiltonian,
$[\rho,H]=0$. In that case, the state of the system can be equivalently
described by a probability vector $\bm{p}$ of eigenvalues of $\rho$,
corresponding to populations in the energy eigenbasis. The thermal equilibrium
state of the system, given by Eq. (1) with $H_{E}$ replaced by $H$, is then
represented by a vector of thermal populations $\bm{\gamma}$. We will denote
the populations of generic input and output states by bold lowercase letters,
$\bm{p}$ and $\bm{q}$, respectively. When referring to the components of these
vectors, we use non-bold symbols with a lower index numbering the components,
e.g., the $i$-th component of $\bm{p}$ is denoted as $p_{i}$. The crucial
point now is how the dynamics of the system interacting with the thermal bath
is modelled. In what follows, we will review two frameworks, which can be seen
as extreme cases with or without any memory effects involved.

_Thermal operations_ (TOs) framework [21, 22, 24] uses minimal assumptions on
the joint system-bath dynamics by only assuming that the joint system is
closed and thus evolves unitarily, and that this unitary evolution is energy-
preserving. Formally, a set of thermal operations consists of completely
positive trace-preserving maps that transform the state $\rho$ in the
following way:

${\cal
E}(\rho)=\Tr_{E}\left[U\left(\rho\otimes\gamma_{E}\right)U^{\dagger}\right],$
(2)

where $U$ is a joint unitary that commutes with the total Hamiltonian of the
system and the bath

$[U,H\otimes\mathbbm{1}_{E}+\mathbbm{1}\otimes H_{E}]=0,$ (3)

and the environmental Hamiltonian $H_{E}$ is arbitrary. Since there are no
further constraints on $U$, arbitrarily strong correlations can build up
between the system and the bath, and one can expect non-Markovian memory
effects to come into play. At the same time, from the perspective of control
theory, generating an arbitrary TO may require very complex and fine-tuned
control over system-bath interactions [27].

_Markovian thermal processes_ (MTPs) framework [25, 26], on the other hand,
uses typical assumptions of the theory of open quantum systems (weak coupling,
large bath size, quickly decaying correlations, etc.) [18], to argue that the
system undergoes an open dynamics described by a Lindblad master equation [35,
36, 37],

$\frac{d\rho(t)}{dt}=-i\quantity[H,\rho(t)]+\mathcal{L}_{t}\quantity(\rho(t)).$
(4)

In the above, $[\cdot,\cdot]$ denotes the commutator and $\mathcal{L}_{t}$ is
the Lindbladian with the following general form:

$\mathcal{L}_{t}(\rho)=\sum_{i}r_{i}(t)\left[L_{i}(t)\rho
L_{i}(t)^{\dagger}-\frac{1}{2}\Bigl{\\{}L_{i}(t)^{\dagger}L_{i}(t),\rho\Bigl{\\}}\,\\!\right]\\!,\\!$
(5)

with $\\{\cdot,\cdot\\}$ denoting the anticommutator, $L_{i}(t)$ being time-
dependent jump operators, and $r_{i}(t)\geq 0$ being time-dependent non-
negative jump rates. Moreover, the thermal state of the system is a stationary
solution of the dynamics, $\mathcal{L}_{t}(\gamma)=0$, and the Lindbladian
$\mathcal{L}_{t}$ commutes with the generator of the Hamiltonian dynamics
$-i[H,\cdot]$ for all times $t$. Formally, an MTP is then any quantum channel
${\cal E}$ that results from integrating Eq. (4) between 0 and $\tau\geq 0$.
Since the dynamics generated by an MTP arises explicitly from a Markovian
model, there are no memory effects. Also, as shown in Ref. [25], the universal
set of controls that allows one to generate any incoherent state
transformation achievable via MTPs consists only of two-level partial
thermalisations. These transform the populations of two energy levels, $i$ and
$j$, in the following way

$\displaystyle p_{i}$
$\displaystyle\rightarrow(1-\lambda)p_{i}+\lambda\frac{p_{i}+p_{j}}{\gamma_{i}+\gamma_{j}}\gamma_{i},$
(6a) $\displaystyle p_{j}$
$\displaystyle\rightarrow(1-\lambda)p_{j}+\lambda\frac{p_{i}+p_{j}}{\gamma_{i}+\gamma_{j}}\gamma_{j},$
(6b)

where $\lambda\in[0,1]$.

Figure 2: Thermal operations vs Markovian thermal processes. Sets of states
that a three-level system with an equidistant energy spectrum $\bm{E}=(0,1,2)$
and prepared in an energy-incoherent state $\bm{p}=(0.7,0.2,0.1)$ (depicted by
a black dot $\bullet$) can be transformed to [green region $C_{+}(\bm{p})$] or
transformed from [blue region $C_{-}(\bm{p})$] by (a) thermal operations and
(b) Markovian thermal processes with respect to inverse temperature
$\beta=0.3$. In (c) we show the overlap of sets of achievable states via TOs
(dark green) and MTPs (light green). The thermal state of the system is
depicted by a black star at the intersection of the dashed lines.

The set of states $C^{\textrm{TO}}_{+}(\bm{p})$ achievable via thermal
operations from a given incoherent initial state $\bm{p}$ can be fully
characterised using the notion of thermomajorisation [22] (see Appendix A-2
for more details). The so-called future thermal cone
$C^{\textrm{TO}}_{+}(\bm{p})$ [38] is a convex set that consists of at most
$d!$ extreme points, the construction of which was given in Lemma 12 of Ref.
[27] (see Fig. 2a for an example with a three-level system). On the other
hand, the set of states $C^{\textrm{MTP}}_{+}(\bm{p})$ achievable via
Markovian thermal processes from a state $\bm{p}$ was recently characterised
using the notion of continuous thermomajorisation [25] (see Appendix A-3 for
an extended discussion). The future Markovian thermal cone
$C^{\textrm{MTP}}_{+}(\bm{p})$ is not convex (as illustrated in Fig. 2b for a
three-level system), but Theorem 4 of Ref. [25] provides a construction of its
extreme points using sequences of two-level full thermalisations (i.e.,
transformations from Eqs. (6a)-(6b) with $\lambda=1$). As can be seen in Fig.
2c, $C^{\textrm{MTP}}_{+}(\bm{p})\subset C^{\textrm{TO}}_{+}(\bm{p})$ and the
difference between these two sets of thermodynamically accessible states
arises purely from memory effects.

A recap of _majorisation_ , _thermomajorisation_ , and _continuous
thermomajorisation_ , three notions needed to fully understand the
convertibility of states under TOs and MTPs, is presented in Appendix A. We
will refrain from restating it here. Instead, we provide a summary in Fig. 3,
illustrating the existence of a given thermodynamic transformation and its
conditions as expressed by these partial order relations.

##  III Bridging the gap with memory

We begin this section by explaining the main building block of this work,
namely the notion of memory-assisted Markovian thermal processes. Then, we
demonstrate how energy-incoherent states achievable from a given initial state
$\bm{p}$ via thermal operations [i.e., any $\bm{q}\in
C_{+}^{\mathrm{TO}}(\bm{p})$] can be approached arbitrarily well using memory-
assisted Markovian thermal processes with large enough memory. We will start
by simplifying the problem and showing that it is sufficient to only consider
the achievability of the extreme points of $C_{+}^{\mathrm{TO}}(\bm{p})$.
Next, we will introduce MeMTP protocols that will serve us to approach extreme
points of $C^{\textrm{TO}}_{+}(\bm{p})$ using MTPs acting on the system and
memory state, $\bm{p}\otimes\bm{\gamma}_{M}$. Finally, we will analyse the
performance of these protocols, i.e., we will show how well they approximate
the desired transformations as the size of the memory $N$ grows. Due to
structural differences, we will do this separately for the infinite
temperature limit and the case of finite temperatures.

###  III.1 Memory-assisted Markovian thermal processes

In this work we want to interpolate between the two extreme regimes of
arbitrarily strong and no memory effects described by TO and MTP frameworks.
We will achieve this by focusing on the more restrictive MTP framework and
extending it by explicitly modelling memory effects by bringing ancillary
systems in thermal states, that will be discarded at the end. More precisely,
we consider MTPs acting on a composite system consisting of the main
$d$-dimensional system in a state $\rho$ and an $N$-dimensional memory system
prepared in its thermal state $\gamma_{M}$ (i.e., given by Eq. (1) with
$H_{E}$ replaced by the Hamiltonian $H_{M}$ of the memory system, which can be
arbitrary). The thermality of the ancillary system $M$ is crucial, as this way
we ensure that no extra thermodynamic resources are brought in unaccounted,
and the only role played by $M$ is to bring extra dimensions that can act as a
memory. As already explained in the introduction, this can also be viewed as
having control over the small $N$-dimensional part of the bath. Formally, we
define the following set of quantum channels.

######  Definition 1 (Memory-assisted MTPs).

A quantum channel ${\cal E}$ is called a memory-assisted Markovian thermal
process (MeMTP) with memory of size $N$, if it can be written as

${\cal E}(\rho)=\Tr_{M}[{\cal E}_{\mathrm{MTP}}(\rho\otimes\gamma_{M})],$ (7)

where ${\cal E}_{\mathrm{MTP}}$ is a Markovian thermal process acting on the
original system extended by an $N$-dimensional ancillary system $M$ prepared
in a thermal state $\gamma_{M}$.

Figure 3: Interpolating between extreme regimes. Arrows between distributions
represent the existence of specific thermodynamic transformations, whose
existence is determined by partial-order relations: thermomajorisation
$\succ_{\beta}$ for thermal operations and continuous thermomajorisation β for
Markovian thermal processes. Our results demonstrate that the gap between
these two frameworks can be bridged with the use memory-assisted MTPs
employing ancillary memory systems of growing dimension $N$ and prepared in
thermal states.

As already mentioned, in this work we will focus on transformations between
energy-incoherent states. Our aim is to show that the sets of states
achievable from a given $\bm{p}$ via memory-assisted MTPs interpolate between
$C^{\textrm{MTP}}_{+}(\bm{p})$ (for $N=1$) and $C^{\textrm{TO}}_{+}(\bm{p})$
(for $N\to\infty$). The illustration of the growing strength of memory effects
captured by our framework is presented in Fig. 3. As a final note, observe
that this framework can be formally related to a particular kind of catalytic
transformations [39, 23]. This is because the ancillary memory system can
always be thermalised at the end of the process and this way be brought to the
initial state.

###  III.2 Simplification to extreme points

We start by recalling the following notion that is crucial for our analysis.

######  Definition 2 ($\beta$-ordering).

Let $\bm{p}$ be an arbitrary energy-incoherent state of a $d$-dimensional
system, and $\bm{\gamma}$ denote the corresponding thermal Gibbs state. The
$\beta$-ordering of $\bm{p}$ is defined as the permutation $\pi_{\bm{p}}$ that
arranges the vector $\quantity(p_{1}/\gamma_{1},...,p_{d}/\gamma_{d})$ in a
non-increasing order, i.e.,

$\bm{p}^{\,\beta}=\quantity(p_{\bm{\pi}^{-1}_{\bm{p}}(1)},...,p_{\bm{\pi}^{-1}_{\bm{p}}(d)}).$
(8)

The $d$-dimensional matrix representation of $\pi_{\bm{p}}$ will be denoted by
$\Pi_{\bm{p}}$, i.e., $\Pi_{\bm{p}}\bm{p}=\bm{p}^{\beta}$.

The future thermal cone $C_{+}^{\mathrm{TO}}(\bm{p})$ is a polytope with at
most $d!$ extreme points, one for each possible $\beta$-order $\pi$ [27]. We
will denote them by $\bm{p}^{\pi}$ (in particular it means that
$\bm{p}^{\pi_{\bm{p}}}=\bm{p}$). Now, we will use two crucial observations.
First, in Ref. [27] it was shown that

$\bm{q}\in C_{+}^{\mathrm{TO}}(\bm{p})\Rightarrow\bm{q}\in
C_{+}^{\mathrm{TO}}\quantity(\bm{p}^{\pi_{\bm{q}}}),$ (9)

meaning that all states with a $\beta$-order $\pi$ that can be achieved from
$\bm{p}$ via thermal operations can also be achieved starting from
$\bm{p}^{\pi}$. And second, it was shown in Ref. [25] that

$\left[\bm{q}\in
C_{+}^{\mathrm{TO}}(\bm{p})~{}~{}\mathrm{and}~{}~{}\pi_{\bm{q}}=\pi_{\bm{p}}\right]\Rightarrow\bm{q}\in
C_{+}^{\mathrm{MTP}}(\bm{p}),$ (10)

meaning that within the same $\beta$-order as the initial state, the subsets
of states achievable via TOs and via MTPs do coincide. As a result, if one can
construct memory-assisted MTPs that reach all the extreme points of
$C_{+}^{\mathrm{TO}}(\bm{p})$, then one can also get to every state in
$C_{+}^{\mathrm{TO}}(\bm{p})$ via MeMTPs. This is done by simply first
transforming $\bm{p}$ to a given extreme point $\bm{p}^{\pi}$ of
$C_{+}^{\mathrm{TO}}(\bm{p})$, and then using MTPs to get from $\bm{p}^{\pi}$
to every state with a $\beta$-order $\pi$ in $C_{+}^{\mathrm{TO}}(\bm{p})$.

Note that the upper index on probability distributions serves a dual purpose.
First, it is used to denote the $\beta$-ordered version of the state. For
example, $\bm{p}^{\beta}$ represents the $\beta$-ordered counterpart of
$\bm{p}$. Second, an upper index $\pi$ is used to denote particular extreme
points of the state’s future. More precisely, $\bm{p}^{\pi}$ represents the
extreme point of the future of $\bm{p}$, with $\beta$-order given by $\pi$.
Notably, $\pi_{\bm{p}}$ represents the $\beta$-order of the state $\bm{p}$.

In order to quantify how well a given state in $C_{+}^{\mathrm{TO}}(\bm{p})$
can be approximated, we will use the total variation distance defined by

$\delta\quantity(\bm{p},\bm{q}):=\frac{1}{2}\sum_{i=1}^{d}\absolutevalue{p_{i}-q_{i}}.$
(11)

From the discussion above, it should be clear that if we can construct MeMTP
protocols approximating every extreme point with an error at most $\epsilon$,
then

$\forall\bm{q}\in C_{+}^{\mathrm{TO}}\quantity(\bm{p}):\quad\min_{{\cal
P}\in\mathrm{MeMTP}}\delta\quantity({\cal P}(\bm{p}),\bm{q})\leq\epsilon.$
(12)

A particular subset of extreme points of $C_{+}^{\mathrm{TO}}(\bm{p})$ that we
will investigate in more detail is given by those extreme states that can be
achieved via sequences of $\beta$-swaps. A $\beta$-swap $\Pi_{ij}^{\beta}$ can
be seen as a thermodynamic analogue of a population swap between levels $i$
and $j$ [27]:

$\Pi^{\,\beta}_{ij}:=\begin{bmatrix}1-e^{-\beta{(E_{j}-E_{i})}}&1\\\
e^{-\beta{(E_{j}-E_{i})}}&0\end{bmatrix}\oplus\mathbbm{1}_{{\backslash}(ij)},$
(13)

with $E_{i}\leq E_{j}$ and $\mathbbm{1}_{{\backslash}(ij)}$ denoting the
$(d-2)\times(d-2)$ identity matrix on the subspace of all energy levels except
$i,j$. Note that in the infinite temperature limit ($\beta=0$), the above
recovers a transposition on levels $i$ and $j$, which we will simply denote by
$\Pi_{ij}$. In this limiting case, all extreme points of
$C_{+}^{\mathrm{TO}}(\bm{p})$ can be obtained by sequences of transpositions
(that is because an extreme point in that case is of the form $\Pi\bm{p}$ for
a permutation matrix $\Pi$, and every $\Pi$ can be constructed from
transpositions).

For finite temperatures, a $\beta$-swap transforms $\bm{p}$ into the extreme
point $\bm{p}^{\pi}$ if the $\beta$-orders $\pi_{\bm{p}}$ and $\pi$ differ
only by a transposition of adjacent elements [27, 38]. In other words, it
happens for $\Pi_{ij}^{\beta}$ when $\pi_{\bm{p}}(i)=\pi_{\bm{p}}(j)\pm 1$.
More generally, a sequence of $\beta$-swaps with non-overlapping supports will
also produce an extreme point of $C^{\mathrm{TO}}_{+}(\bm{p})$, and so a total
number of extreme points that can be achieved by sequences of $\beta$-swaps
for dimension $d$ (including the starting point) is given by $F(d+1)$, where
$F(k)$ is the $k$-th Fibonacci number [40]. Finally, we will also make use of
the notion of $\beta$-cycles that we now define. For a state $\bm{p}$,
consider a $k$-dimensional subset of energy levels
$\quantity{i_{1},\ldots,i_{k}}$ neighbouring in the $\beta$-order, i.e.,
$\pi_{\bm{p}}(i_{j+1})=\pi_{\bm{p}}(i_{j})+1$. Denote by $\pi$ a cyclic
permutation on this subset, i.e., either $\pi(i_{j})=i_{j+1\text{ mod }k}$, or
$\pi(i_{j})=i_{j-1\text{ mod }k}$. Then, a thermal operation mapping $\bm{p}$
to its extreme point $\bm{p}^{\pi^{\prime}}$ is called a $\beta$-cycle, if
$\Pi^{\prime}=\Pi\Pi_{\bm{p}}$ (here $\Pi,\Pi^{\prime}$ and $\Pi_{\bm{p}}$
denote matrix representations of permutations $\pi,\pi^{\prime}$ and
$\pi_{\bm{p}}$). To emphasise that a given $\beta$-cycle acts on $k$ levels,
we will sometimes refer to it as a $\beta$-$k$-cycle.

To summarise the notion used, we denote a specific permutation as lowercase
letter $\pi$, which acts on integers. Its matrix representation is represented
by capital $\Pi$, and, in particular, the transposition of elements $i$ and
$j$ has its matrix counterpart $\Pi_{ij}$. We use the upper index $\beta$,
e.g., $\Pi^{\beta}$, when referring to an extreme operation that permutes the
$\beta$-order of the state by a permutation $\pi$.

###  III.3 Memory-assisted protocols

The basic building blocks of all our protocols are given by two-level
elementary thermalisations that are formally defined as follows.

######  Definition 3 (Two-level thermalisations).

Consider a system in a state $\bm{p}$ with the corresponding thermal state
$\bm{\gamma}$. Then, an MTP transformation

$\\{p_{i},p_{j}\\}\to\left\\{\frac{p_{i}+p_{j}}{\gamma_{i}+\gamma_{j}}\gamma_{i},\frac{p_{i}+p_{j}}{\gamma_{i}+\gamma_{j}}\gamma_{j}\right\\}$
(14)

is called a two-level thermalisation between levels $i$ and $j$, and the
corresponding matrix acting on probability vectors will be denoted by
$T_{ij}$. Moreover, if $\pi_{\bm{p}}(i)=\pi_{\bm{p}}(j)\pm 1$, then $T_{ij}$
is called a neighbour thermalisation.

Let us note that the importance of neighbour thermalisations and the reason we
employ them in our protocols stems from the fact that their sequences produce
the extreme points of the Markovian thermal cone $C_{+}^{\mathrm{MTP}}$ [25].
Intuitively, one can expect that in order to approximate extreme states of
$C_{+}^{\mathrm{TO}}(\bm{p})$ using MeMTPs, one should get to the extreme
points of $C_{+}^{\mathrm{MTP}}(\bm{p}\otimes\bm{\gamma}_{M})$, and these can
be achieved by neighbour thermalisations of the composite system-memory state.

Before delving into the full details of our protocol for approximating
$\beta$-swaps, let us start with a high-level description to provide some
insight into our investigation. We begin with the simplest case of a two-level
system and a two-dimensional memory, drawing an analogy between continuous
(thermo)majorisation and connected vessels (see Fig. 4 for a schematic
representation). Considering two vessels—one filled with liquid and one
empty—the most one can do when they are connected is to equalise the levels of
the liquid between them. However, by adhering to the simple schematic provided
in Fig. 4, it is possible to exceed this intuitively unbeatable limit.

Figure 4: Simplest example using connected vessels analogy. Continuous
(thermo)majorisation on $d$-level probability vectors is equivalent to a task
of redistributing the content of $d$-ordered vessels that are connected
pairwise. Adding a memory in a Gibbs state is then akin to multiplying glasses
– one empty and one full glass become $N$ pairs of full and empty glasses. The
process involves five steps, read from top to bottom and left to right, which
represent the simplest protocol that allows shifting more than half of the
liquid from full to empty glasses. The final distribution after the protocol
is applied is $(3/8,5/8)$.

The sequence in Fig. 4 should be read from panels (a) to (f), in a top-down
and left-to-right manner. We begin with four vessels: two are half-filled and
two are empty. We can connect these vessels in pairs, thus equalising the
fluid levels. The first two steps involve connecting the half-filled vessels
sequentially to the first empty vessel. Likewise, the next two steps connect
both initially half-filled vessels to the second empty one. The final step,
which involves equalising the fluid levels in pairs, is analogous to
thermalising the memory. An astute reader can confirm that $5/8$ of the total
fluid ends up in the vessels that were initially empty, thereby surpassing the
$1/2$ limit.

We now describe our proposition for a MeMTP protocol approximating the
$\beta$-swap $\Pi^{\,\beta}_{ij}$ between the $i$-th and $j$-th energy levels
of the main system. It involves a sequence of $N^{2}$ two-level
thermalisations of the state of the composite system (see Fig. 5), which
includes the main system and a memory starting at thermal equilibrium, i.e., a
state $\bm{p}\otimes\bm{\gamma}_{M}$. In particular, we focus on the
populations
$[\bm{p}\,\otimes\,\bm{\gamma}_{M}]_{N(i-1)+1},\ldots,[\bm{p}\,\otimes\,\bm{\gamma}_{M}]_{Ni}$
corresponding to the $i$-th level of the main system and similarly for the
$j$-th level. The protocol can be split into a sequence of $N$ rounds
$\mathcal{R}_{k}^{(ij)}$ with $k=1,\ldots,N$ consisting of $N$ steps each
(shaded area in Fig. 5). In the $k$-th round, we select the entry
$[\bm{p}\otimes\bm{\gamma}_{M}]_{N(i-1)+k}$ and thermalise it sequentially
with all the levels corresponding to the level $j$ of the main system:

$\mathcal{R}^{(ij)}_{k}(\bm{p}\otimes\bm{\gamma}_{M}):=\quantity(\prod_{l=1}^{N}T_{(i-1)N+k,\,(j-1)N+l})(\bm{p}\otimes\bm{\gamma}_{M}).$
(15)

Note that if $\pi_{\bm{p}}(i)=\pi_{\bm{p}}(j)\pm 1$ (i.e., the $\beta$-orders
of $\bm{p}$ and $\Pi^{\beta}_{ij}\bm{p}$ differ by a transposition of adjacent
elements), then all thermalisations performed are neighbour thermalisations.
Using the above, we can now define the action of the truncated protocol
$\widetilde{\mathcal{P}}^{(ij)}$:

$\widetilde{\mathcal{P}}^{(ij)}(\bm{p}\otimes\bm{\gamma}_{M}):=\mathcal{R}^{(ij)}_{N}\circ\ldots\circ\mathcal{R}^{(ij)}_{1}(\bm{p}\otimes\bm{\gamma}_{M}).$
(16)

Figure 5: $\beta$-swap protocol. A two-level subsystem of a generically
$d$-dimensional system, represented by blue and red squares, is extended by an
$N$-dimensional memory represented by grey squares. The composite
$2N$-dimensional system undergoes $N$ rounds of processing, where the $k$-th
round involves $N$ sequential two-level thermalisations of the first $N$
entries with the $(N+k)$-th entry (represented by the shaded colour around the
squares). After the final thermalisation step, the memory can be discarded.

The final step is to decouple the main system from the memory using a full
thermalisation ${\cal T}$ of the memory system $M$, which acts on a general
joint state $\bm{Q}$ as:

$\mathcal{T}(\bm{Q})=\bm{q}\otimes\bm{\gamma}_{M},\qquad
q_{i}=\sum_{j=1}^{N}Q_{N(i-1)+j}.$ (17)

Thus, the full protocol approximating a $\beta$-swap $\Pi_{ij}^{\beta}$ is
given by

$\mathcal{P}^{(ij)}(\bm{p}\otimes\bm{\gamma}_{M})=\mathcal{T}\circ\widetilde{\mathcal{P}}^{(ij)}(\bm{p}\otimes\bm{\gamma}_{M}).$
(18)

Putting the whole protocol in a simple words – we take all “full” levels,
connect them to the first “empty” level sequentially, and then repeat the
process for all the empty levels.

We will also employ more general protocols that aim at approximating the
transformation of the initial state $\bm{p}$ into an extreme state of
$C_{+}^{\mathrm{TO}}(\bm{p})$ given by $\bm{p}^{\pi^{\prime}}$. Denote matrix
representations of $\beta$-orders of these states by $\Pi_{\bm{p}}$ and
$\Pi^{\prime}=\Pi\Pi_{\bm{p}}$ for some permutation matrix $\Pi$. Moreover,
let us decompose $\Pi$ into neighbour transpositions with respect to $\bm{p}$,
i.e., we write $\Pi=\Pi_{i_{m}j_{m}}\dots\Pi_{i_{1}j_{1}}$ with every
consecutive transposition $\Pi_{i_{k}j_{k}}$ changing the $\beta$-order of the
state $\Pi_{i_{k-1}j_{k-1}}\dots\Pi_{i_{1}j_{1}}\bm{p}$ only by a
transposition of adjacent elements. Then, we define the following two
protocols to approximate $\bm{p}^{\pi^{\prime}}$:

$\displaystyle\mathcal{P}^{{\Pi}}$
$\displaystyle:=\mathcal{P}^{(i_{m}j_{m})}\circ\dots\circ\mathcal{P}^{(i_{1}j_{1})},$
(19a) $\displaystyle\widetilde{\mathcal{P}}^{{\Pi}}$ $\displaystyle:={\cal
T}\circ\widetilde{\mathcal{P}}^{(i_{m}j_{m})}\circ\dots\circ\widetilde{\mathcal{P}}^{(i_{1}j_{1})}.$
(19b)

Note that, by construction, all two-level thermalisations performed in the
above protocols are neighbour thermalisations.

Let us summarise the notation that we introduced for the protocols. Rounds of
the protocol for approximating a swap between levels $i$ and $j$ of the main
system are denoted by $\mathcal{R}^{(ij)}$. Similarly, a full protocol with
and without final thermalization is denoted by $\mathcal{P}^{(ij)}$ and
$\tilde{\mathcal{P}}^{(ij)}$, respectively. Final thermalisation is denoted by
$\mathcal{T}$.

As a last remark, it should be noted that the protocol discussed in this
section is not unique. Alternative protocols achieving equivalent $\beta$-swap
approximations can be found in Appendix F. Additionally, we provide a
comparison of the convergence rates of different algorithms toward a specified
target state.

###  III.4 Achieving extreme points of $C_{+}^{\mathrm{TO}}$ for $\beta=0$

We are now ready to state our main results concerning the power of memory-
assisted Markovian thermal processes in the infinite temperature limit. Let us
recall that we focus on a $d$-level system and an $N$-dimensional memory in
energy-incoherent states represented by probability distributions $\bm{p}$ and
$\bm{\gamma}_{M}$. Since $\beta=0$, the thermal state of the memory is
described by a uniform distribution $\bm{\eta}_{M}$ with every entry equal to
$1/N$. We start with the following lemma, the proof of which can be found in
Appendix C (with the necessary background on the mathematical tools used
presented in Appendix B).

######  Lemma 1 (Memory-assisted transposition).

In the infinite temperature limit, $\beta=0$, and for an $N$-dimensional
memory, the MeMTP protocol $\mathcal{P}^{(ij)}$ acts as

$\mathcal{P}^{(ij)}(\bm{p}\otimes\bm{\eta}_{M})=\bm{q}\otimes\bm{\eta}_{M},$
(20)

with

$\bm{q}=\left(\Pi_{ij}+\epsilon\quantity(\mathbbm{1}-\Pi_{ij})\right)\bm{p},$
(21)

and $\epsilon$ given by

$\epsilon=(\pi
N)^{-1/2}+o\quantity(N^{-1/2})\overset{N\rightarrow\infty}{\longrightarrow}0.$
(22)

It is well known that any permutation of $d$ elements can be decomposed into a
product of at most $d$ transpositions or $\binom{d}{2}$ neighbour
transpositions. Therefore, by employing Lemma 1, we can demonstrate that an
arbitrary permutation can be achieved using a composition of our approximate
protocols.

######  Theorem 2 (Memory-assisted permutation).

In the infinite temperature limit, $\beta=0$, and for an $N$-dimensional
memory, $\Pi$ can be approximated by the MeMTP protocol $\mathcal{P}^{{\Pi}}$
as follows:

$\mathcal{P}^{{\Pi}}(\bm{p}\otimes\bm{\eta}_{M})=\bm{q}\otimes\bm{\eta}_{M},$
(23)

where

$\bm{q}=\left(\Pi+\epsilon\bm{\Delta}+o\quantity(N^{-1/2})\right)\bm{p}\overset{N\rightarrow\infty}{\rightarrow}\Pi\bm{p},$
(24)

with $\epsilon=(\pi N)^{-1/2}$ and the operator $\bm{\Delta}$ defined in terms
of transpositions appearing in the definition of ${\cal P}^{\Pi}$ in Eq.
(19a):

$\bm{\Delta}=\sum_{l=1}^{m}\quantity(\prod_{k=l+1}^{m}\Pi_{i_{k}j_{k}})\quantity(\mathbbm{1}-\Pi_{i_{l}j_{l}})\quantity(\prod_{k=1}^{l-1}\Pi_{i_{k}j_{k}}).$
(25)

###### Proof.

From the definition of $\mathcal{P}^{\Pi}$ and Lemma 1 we get

$\bm{q}=\quantity[\Pi_{i_{m}j_{m}}+\epsilon\quantity(\mathbbm{1}-\Pi_{i_{m}j_{m}})]\dots\quantity[\Pi_{i_{1}j_{1}}+\epsilon\quantity(\mathbbm{1}-\Pi_{i_{1}j_{1}})]\bm{p}.$
(26)

Clearly, the leading term is given by $\Pi\bm{p}$, whereas the next leading
term, proportional to $\epsilon$, is given by $\epsilon\Delta\bm{p}$. All
higher order terms scale at least as $\epsilon^{2}$, so are of the order
$o(N^{-1/2})$. ∎

The above theorem can then be directly used to obtain the bound on how close
one can get from a given $\bm{p}$ to any state $\bm{q}\in
C_{+}^{\mathrm{TO}}(\bm{p})$ using MeMTPs with $N$-dimensional memory. We
explain how to derive such a bound for a given $\bm{p}$ in Appendix D, whereas
below we present a weaker, but much simpler, bound that is independent of
$\bm{p}$.

######  Corollary 3.

Consider states $\bm{p}$ and $\bm{q}\in C_{+}^{TO}(\bm{p})$. Then, in the
infinite temperature limit, $\beta=0$, and for an $N$-dimensional memory,
there exists a MeMTP protocol $\mathcal{P}$ such that

${\cal P}(\bm{p}\otimes\bm{\eta}_{M})=\bm{q}^{\prime}\otimes\bm{\eta}_{M},$
(27)

with

$\delta(\bm{q}^{\prime},\bm{q})\leq\frac{d(d-1)}{2\sqrt{\pi
N}}+o\quantity(N^{-1/2}).$ (28)

###### Proof.

First, define $\Pi$ as a permutation that changes the $\beta$-order of
$\bm{p}$ to that of $\bm{q}$. In other words, the $\beta$-order of $\Pi\bm{p}$
is $\pi_{\bm{q}}$. Then, using Theorem 2, we have that

${\cal P}^{\Pi}(\bm{p}\otimes\bm{\eta}_{M})=\bm{r}\otimes\bm{\eta}_{M}$ (29)

with

$\delta(\bm{r},\Pi\bm{p})\simeq\frac{1}{2\sqrt{\pi
N}}\sum_{i=1}^{d}|(\Delta\bm{p})_{i}|\lesssim\frac{m}{\sqrt{\pi
N}}\lesssim\frac{d(d-1)}{2\sqrt{\pi N}},$ (30)

where $\simeq$ and $\lesssim$ denote the equalities and inequalities up to
$o(N^{-1/2})$. In the above, we have used the triangle inequality and the fact
that one can always decompose $\Pi$ into at most $d(d-1)/2$ neighbour
transpositions. Next, from Eq. (10), we know that there exists an MTP protocol
${\cal P}^{\prime}$ mapping $\Pi\bm{p}$ to $\bm{q}$. Using the contractiveness
of the total variation distance under stochastic processing, we then have

$\displaystyle\delta({\cal P}^{\prime}(\bm{r}),\bm{q})$
$\displaystyle=\delta({\cal P}^{\prime}(\bm{r}),{\cal
P}^{\prime}(\Pi\bm{p}))\leq\delta(\bm{r},\Pi\bm{p})\lesssim\frac{d(d-1)}{2\sqrt{\pi
N}}.$ (31)

We thus conclude that by choosing ${\cal P}=({\cal
P}^{\prime}\otimes\mathcal{I}_{M})\circ{\cal P}^{\Pi}$, Eqs. (27)-(28) are
satisfied. ∎

Furthermore, we present the following conjecture for a better approximation of
arbitrary permutations.

######  Conjecture 1 (Improved convergence).

In the infinite temperature limit, $\beta=0$, and for an $N$-dimensional
memory, $\widetilde{\mathcal{P}}^{{\Pi}}$ gives a better approximation of a
permutation $\Pi$ than $\mathcal{P}^{\Pi}$:

$\delta\left(\Pi\bm{p},\tilde{\bm{q}}\right)\leq\delta\left(\Pi\bm{p},\bm{q}\right),$
(32)

where $\tilde{\bm{q}}$ and $\bm{q}$ are defined via

$\widetilde{\mathcal{P}}^{{\Pi}}(\bm{p}\otimes\bm{\eta}_{M})=\tilde{\bm{q}}\otimes\bm{\eta}_{M},\qquad\mathcal{P}^{{\Pi}}(\bm{p}\otimes\bm{\eta}_{M})=\bm{q}\otimes\bm{\eta}_{M}.$
(33)

The conjecture is solidified by strong numerical evidence (see Fig. 6 for an
example considering $d=6$). We note that the convergence is better, but the
overall character of $O(N^{-1/2})$ is still preserved. More specifically, we
observe that for permutations given by $\beta$-$k$-cycles with $k\leq d$,
there is no advantage to removing the intermediate thermalisations (i.e., no
advantage of $\widetilde{{\cal P}}^{\Pi}$ over ${{\cal P}}^{\Pi}$). The
advantage already appears for a composition of $\beta$-$d$-cycle with
$\beta$-$(d-1)$-cycle, leading to the $\beta$-order $(d,d-1,1,\ldots,d-2)$
(here, without loss of generality, we assumed that the initial $\beta$-order
is given by $(1,2,\dots,d)$). In general, the advantage grows with the number
of composed $\beta$-cycles (see Fig. 6, where different colours and markers
correspond to different length compositions of $\beta$-cycles). In particular,
we verified that for a permutation $(16)(25)(34)$, which is composed of
$\beta$-cycles of length $6$ through $2$ (or $15=\binom{6}{2}$ neighbour
transpositions), both $\mathcal{P}^{\Pi}$ and $\widetilde{\mathcal{P}}^{\Pi}$
converge to the actual extreme point $\Pi\bm{p}$. Surprisingly, we find that
all the other possible permutations fall within the convergence advantage
class of one of the aforementioned $\beta$-cycle compositions. This includes,
in particular, the cases when the last $\beta$-cycle in the sequence is
incomplete, i.e., it is shortened by the final subsequence of $\beta$-swaps of
any length.

Figure 6: Convergence rates at infinite temperature. Log-log plot of the total
variation distance $\delta$ between the extreme points $\bm{p}^{\pi}\in
C_{+}^{\mathrm{(TO)}}(\bm{p})$ and the states obtained from $\bm{p}$ via the
algorithm $\mathcal{P}^{\Pi}$ (left panel) and $\widetilde{\mathcal{P}}^{\Pi}$
(right panel), as a function of the memory size $N$. Here,
$\bm{p}=\quantity(0.37,0.24,0.16,0.11,0.07,0.05)$, $\beta=0$, and different
colours correspond to families of extreme points $\bm{p}^{\pi}$ with different
convergence rates (from bottom to top $\bm{p}^{\pi}$ is obtained from $\bm{p}$
via a $\beta$-6-cycle, a composition of a $\beta$-6-cycle with a
$\beta$-5-cycle, and so on). All convergences behave as $O(N^{-1/2})$ as given
in Eq. (28), which can be seen by the comparison with the function
$0.4/\sqrt{N}$ (dashed black line), with multiplicative advantage for
$\widetilde{\mathcal{P}}^{\Pi}$ over $\mathcal{P}^{\Pi}$.

###  III.5 Achieving extreme points of $C_{+}^{\mathrm{TO}}$ for $\beta\neq 0$

Our second main result concerns the power of memory-assisted Markovian thermal
processes at finite temperatures. We start with the following generalisation
of Lemma 1, the proof of which can be found in Appendix C.

######  Theorem 4 (Memory-assisted $\beta$-swap).

For a finite temperature, $\beta\neq 0$, and for an $N$-dimensional memory
described by a trivial Hamiltonian (so that its thermal state is
$\bm{\eta}_{M}$), $\Pi_{ij}^{\beta}$ can be approximated by the MeMTP protocol
${\cal P}^{(ij)}$ as follows:

${\cal P}^{(ij)}(\bm{p}\otimes\bm{\eta}_{M})=\bm{q}\otimes\bm{\eta}_{M},$ (34)

with

$\delta(\bm{q},\Pi_{ij}^{\beta}\bm{p})=\frac{\quantity(4\Gamma_{i}\Gamma_{j})^{N}}{\quantity(\Gamma_{i}-\Gamma_{j})^{2}}\quantity[\frac{\absolutevalue{p_{i}\Gamma_{j}-p_{j}\Gamma_{i}}}{(N+1)\sqrt{\pi
N}}+o\quantity(N^{-3/2})]$ (35)

where we have used $\Gamma_{i}=\gamma_{i}/(\gamma_{i}+\gamma_{j})$ and
likewise for $\Gamma_{j}$.

By using the above theorem, one can approximate with arbitrary precision a
total of $F(d+1)$ extreme points achievable by a composition of non-
overlapping $\beta$-swaps (recall that $F(k)$ is the $k$-th Fibonacci number).
However, since $F(d+1)\leq d!$ for $d\geq 3$, not all extreme points of
$C_{+}^{\mathrm{TO}}(\bm{p})$ can be obtained this way. Nevertheless, we
conjecture that using MeMTP protocols $\widetilde{{\cal P}}^{\Pi}$ that are
composed of blocks imitating $\beta$-swaps, just without intermediate
thermalisations, one can reach all the extreme points of
$C_{+}^{\mathrm{TO}}(\bm{p})$.

######  Conjecture 2 (Extreme points of $C_{+}^{\mathrm{TO}}$).

Consider a state $\bm{p}$ and the extreme point of
$C_{+}^{\mathrm{TO}}(\bm{p})$ given by $\bm{p}^{\pi^{\prime}}$, with matrix
representations of $\beta$-orders of these states satisfying
$\Pi^{\prime}=\Pi\Pi_{\bm{p}}$ for some permutation $\Pi$. Then, for a finite
temperature, $\beta\neq 0$, and for an $N$-dimensional memory described by a
trivial Hamiltonian (so that its thermal state is $\bm{\eta}_{M}$), the MeMTP
protocol $\widetilde{{\cal P}}^{\Pi}$ acts as

$\widetilde{{\cal
P}}^{\Pi}(\bm{p}\otimes\bm{\eta}_{M})=\bm{q}\otimes\bm{\eta}_{M},$ (36)

with

$\delta(\bm{q},\bm{p}^{\pi^{\prime}})\overset{N\rightarrow\infty}{\longrightarrow}0.$
(37)

The conjecture is solidified by the following two pieces of evidence, which
utilize the truncated protocols $\tilde{\mathcal{P}}^{\Pi}$. First, we provide
an analytical proof for convergence of the truncated protocols to a subset of
extreme points beyond the ones achievable with simple $\beta$-swaps.

######  Theorem 5 (Memory-assisted $\beta$-3-cycle).

Consider a state $\bm{p}$ with entries $i_{1},i_{2},i_{3}$ being neighbours in
the $\beta$-order (i.e.,
$\pi_{\bm{p}}(i_{1})=\pi_{\bm{p}}(i_{2})+1=\pi_{\bm{p}}(i_{3})+2$), and the
extreme point of $C_{+}^{\mathrm{TO}}(\bm{p})$ given by
$\bm{p}^{\pi^{\prime}}$, with matrix representations of $\beta$-orders of
these states satisfying $\Pi^{\prime}=\Pi\Pi_{\bm{p}}$ for
$\Pi=\Pi_{i_{1}i_{3}}\Pi_{i_{2}i_{3}}$. Then, for a finite temperature,
$\beta\neq 0$, and for an $N$-dimensional memory described by a trivial
Hamiltonian (so that its thermal state is $\bm{\eta}_{M}$), the MeMTP protocol
$\widetilde{{\cal P}}^{\Pi}$ acts as

$\widetilde{{\cal
P}}^{\Pi}(\bm{p}\otimes\bm{\eta}_{M})=\bm{q}\otimes\bm{\eta}_{M}$ (38)

with

$\delta(\bm{q},\bm{p}^{\pi^{\prime}})\overset{N\rightarrow\infty}{\rightarrow}0.$
(39)

The proof of the above theorem can be found in Appendix E, and potentially the
same proving techniques can be applied to higher-order cycles. This would then
provide a general method for simulating $\beta$-cycles, as well as any
combinations of non-overlapping $\beta$-cycles, with arbitrary precision
through MeMTPs.

Figure 7: Convergence rates at finite temperature. Log-log plot of the total
variation distance $\delta$ between the states obtained from $\bm{p}$ via the
algorithm $\widetilde{\mathcal{P}}^{\widetilde{\Pi}_{ij}}$ and the
corresponding extreme points $\bm{p}^{\pi}\in C_{+}^{\mathrm{(TO)}}(\bm{p})$,
as a function of the memory size $N$. Here,
$\bm{p}=\quantity(0.37,0.24,0.16,0.11,0.07,0.05)$, $\beta=0.1$, and different
colours correspond to extreme points $\bm{p}^{\pi}$ with matrix representation
of the $\beta$-order $\pi$ given by $\Pi=\widetilde{\Pi}_{ij}\Pi_{\bm{p}}$.
For all curves, the convergence is better than $O(N^{-1/2})$ in agreement with
Eq. (41), as can be seen by the comparison with the limiting line $1/\sqrt{\pi
N}$ for $\beta=0$ (dot-dashed black line).

The second piece of evidence is based on extensive numerical simulations
demonstrating convergence to arbitrary extreme points beyond both
$\beta$-swaps and aforementioned $\beta$-cycles. Let us consider any state
$\bm{p}$ and a set of permutations defined by a recurrence formula

$\widetilde{\Pi}_{ij}:=\quantity(\prod_{k=1}^{j}\Pi_{\pi_{\bm{p}}(k),\,\pi_{\bm{p}}(k+1)})\widetilde{\Pi}_{i-1,\,d+1-i}$
(40)

with $1\leq j\leq d-i$ and assuming the starting condition
$\widetilde{\Pi}_{0d}=\mathbbm{1}$. Note that $\widetilde{\Pi}_{1,d-1}$
represents a full $\beta$-$d$-cycle, $\widetilde{\Pi}_{i,d+1-i}$ a composition
of $i$ $\beta$-cycles of length from $d$ to $d+1-i$, and finally
$\widetilde{\Pi}_{d,1}$ is a permutation which fully reverses the
$\beta$-order of $\bm{p}$. For each such permutation, we have considered the
action of the protocol
$\widetilde{\mathcal{P}}^{\widetilde{\Pi}_{ij}}\quantity(\bm{p})$ and its
convergence to the respective extreme point $\bm{p}^{\pi}$ with
$\Pi=\widetilde{\Pi}_{ij}\Pi_{\bm{p}}$ (recall that $\Pi$ is a matrix
representation of $\pi$). In each case, we have observed the convergence of
the form from Eq. (37) that is better than $N^{-1/2}$. Results for an
exemplary state in dimension $d=6$ are presented in Fig. 7, where a total of
$15$ different curves are shown to lie below the $N^{-1/2}$ limit and
diverging from it.

Finally, based on Theorem 4, the proof of Theorem 5, and numerical evidence,
one can reasonably strengthen Conjecture 2 to make the following statement on
the convergence:

$\delta(\bm{q},\bm{p}^{\pi^{\prime}})=O\quantity(\frac{e^{-A({\Pi})N}}{N^{3/2}}),$
(41)

where $A(\Pi)=O(1)$ is a permutation-dependent exponent.

##  IV Discussion and applications

In Section III, we demonstrated a method of achieving an arbitrary state from
the future cone of TO using MeMTPs through MTP operations acting upon the
system extended by memory, initiated in the thermal state $\gamma$. In the
following sections, we will apply our protocol to study information-based
quantum thermodynamic processes, such as work extraction and cooling. Next, we
revisit the question of the sufficiency of two-level control for TOs. Finally,
we provide a brief discussion of the behaviour of the free energy and
correlations with the progression of our protocol. This sheds light on how
non-Markovian effects arise in the memory-assisted protocol.

###  IV.1 Work extraction

Generally, definitions of work rely either on the control and manipulation of
external parameters that determine the system’s Hamiltonian [41, 42] or on the
assumption that work is a random variable and a controlled Hamiltonian
evolution is employed to determine the work statistics [43, 44, 45, 46]. In
contrast, the resource-theoretic approach avoids the presence of any external
agent and does not involve changes to the Hamiltonian. Furthermore, this
framework differs from traditional approaches by shifting the focus from
average and higher moments of the work distribution to the single-shot regime.
In this regime, the question is posed as to what is the maximum amount of work
that can be extracted while allowing for a probability of failure $\epsilon$.

Figure 8: $\epsilon$-deterministic work extraction with MeMTPS. Transformation
error $\epsilon$ as a function of the work $W$ extracted from a two-level
system with energy splitting $\Delta$ prepared in a thermal state at
temperature $1/\beta_{S}$ smaller than the environmental temperature $1/\beta$
with parameters $\beta_{S}\Delta=2$ and $\beta\Delta=1$. System-environment
interactions are modelled by TOs (dashed black curve), MTPs (dashed red curve)
and memory-assisted Markovian thermal process with a memory of size $2$, $4$ ,
$8$, $16$, $32$, $64$ and $128$, respectively.

The so-called $\epsilon$-deterministic work extraction, which typically
involves an out-of-equilibrium system $S$, a thermal bath at inverse
temperature $\beta$, and a battery $B$ initially in an energy eigenstate
$E_{0}$ [47, 22, 48]. The aim is to increase the energy of $B$ by an amount
$W$ by exciting it from $E_{0}$ to $E_{1}=E_{0}+W$ with a success probability
$1-\epsilon$. The optimal error $\epsilon$ for a given $W$ can be obtained via
thermomajorisation condition for transformations given by thermal operations
[22] and through continuous thermomajorisation relations when transformations
are given by Markovian thermal processes [26]. Our framework allows one to
interpolate between the two extremes by including a memory system with varying
dimension $N$.

Consider a two-level system $S$ and a two-level battery $B$ with energy levels
$(0,\Delta)$ and $(0,W)$, respectively. Assume that the initial state of the
joint system is given by $\bm{p}_{SB}=\bm{p}\otimes(1,0)$. One can then select
the extreme point
$\bm{p}^{\pi^{\prime}}_{SB}\in{C}^{\mathrm{TO}}_{+}\quantity(\bm{p}_{SB})$
from the future thermal cone of the composite system for which the following
relation is satisfied with the minimum value of $\epsilon_{\mathrm{TO}}$:

$\bm{\gamma}\otimes(\epsilon_{\mathrm{TO}},1-\epsilon_{\mathrm{TO}})\in{C}^{\mathrm{TO}}_{+}\quantity(\bm{p}^{\pi^{\prime}}_{SB}).$
(42)

In other words, $\bm{p}^{\pi^{\prime}}_{SB}$ is an intermediate state from
which one can achieve minimal error for extracting $W$ work from $\bm{p}$ via
any thermal operation. We can now define $\Pi$ as a permutation that maps the
matrix representation of the initial $\beta$-order of $\bm{p}_{SB}$ to the
final $\beta$-order of $\bm{p}_{SB}^{\pi^{\prime}}$. Then, by using the
algorithm $\widetilde{\mathcal{P}}^{\Pi}$, we can transform $\bm{p}_{SB}$ into
a state $\bm{q}_{SB}$ that approximates $\bm{p}^{\pi^{\prime}}_{SB}$. Finally,
due to Eq. (10), we can use standard thermomajorisation to find the minimal
value of $\epsilon_{N}$ for which the state $\bm{q}_{SB}$ can be transformed
to $\bm{\gamma}\otimes(\epsilon_{N},1-\epsilon_{N})$ via MTPs. Note that
$\epsilon_{N}$ then corresponds to the probability of failure of extracting
work $W$ from $\bm{p}$ using a memory of size $N$. Numerical simulations of
this procedure (see Fig. 8) show that as $N$ grows, $\epsilon_{N}$ decreases,
allowing us to conjecture that
$\lim_{N\rightarrow\infty}\epsilon_{N}=\epsilon_{\mathrm{TO}}$. However, note
that the convergence is not uniform: it is the slowest around $W=0$ and the
kink at $W=1/\beta\log(1+e^{-\beta\Delta})$. Nevertheless, Fig. 8 clearly
shows that even a small size memory can significantly improve the quality of
the extracted work.

###  IV.2 Cooling a two-level system using a two-dimensional memory with
nontrivial Hamiltonian

As a second application of our findings, we consider the task of cooling a
two-level system with the aid of a two-dimensional memory. The setup involves
a two-level system with energy gap $E_{S}$, extended by a memory system with
energy gap $E_{M}$. The joint system’s energy level structure is depicted in
Fig. 9. We assume that the difference between energy gaps is such that it
allows one to selectively couple with the bath, i.e., $E_{S}-E_{M}\neq E_{M}$.
This enables us to separately address transitions
$\ket{01}\leftrightarrow\ket{10}$, $\ket{00}\leftrightarrow\ket{11}$ together
with two coupled pairs of the form $\ket{0i}\leftrightarrow\ket{1i}$ and
$\ket{i0}\leftrightarrow\ket{i1}$. We will refer to these operations
thermalise these levels as operation 1, 2, 3 and 4, respectively.

Figure 9: Energy level structure. Schematic diagram of a two-level system,
consisting of a main system and a memory with energy gaps $E_{S}$ and $E_{M}$,
respectively. The energy gap of the composite system is such that it can be
selectively coupled to the thermal bath.

Let us now assume, for simplicity, that the system starts in an excited state
extended by Gibbs memory, $\bm{p}\otimes\bm{\gamma}_{M}$ with
$p_{i}=\delta_{i1}$. If we consider the main system alone with access only to
MTPs, one can cool it down only to the ambient temperature. In this case, the
system will reach thermal equilibrium, and the resulting distribution is given
by

$\bm{\gamma}_{S}=\quantity(\frac{1}{1+e^{-\beta E_{S}}},\frac{e^{-\beta
E_{S}}}{1+e^{-\beta E_{S}}}).$ (43)

However, by implementing our protocol, which can be realised as a sequence of
operations, $1\rightarrow 2\rightarrow 3\rightarrow 4$ and discarding
(thermalising) the memory, we arrive at
$\mathcal{P}(\bm{p}\otimes\bm{\gamma}_{M})=\bm{q}\otimes\bm{\gamma}_{M}$ with

$\bm{q}=\matrixquantity[\frac{e^{\beta
E_{M}}+e^{\beta\left(E_{M}+E_{S}\right)}+e^{\beta\left(2E_{M}+E_{S}\right)}+e^{\beta\left(E_{M}+2E_{S}\right)}+e^{\beta
E_{S}}}{\left(e^{\beta
E_{S}}+1\right)\left(e^{\beta\left(E_{M}-E_{S}\right)}+1\right)\left(e^{\beta\left(E_{M}+E_{S}\right)}+1\right)}\\\
\frac{e^{\beta E_{M}}+e^{\beta\left(2E_{M}+E_{S}\right)}+e^{\beta
E_{S}}}{\left(e^{\beta E_{S}}+1\right)\left(e^{\beta E_{M}}+e^{\beta
E_{S}}\right)\left(e^{\beta\left(E_{M}+E_{S}\right)}+1\right)}].$ (44)

The distance of this state from the Gibbs state $\bm{\gamma}_{S}$ at ambient
temperature in terms of the 1-norm is given by

$\norm{\bm{q}-\bm{\gamma}_{S}}_{1}=\frac{1}{\left(e^{-\beta
E_{S}}+1\right)\left[\cosh\left(\beta E_{M}\right)+\cosh\left(\beta
E_{S}\right)\right]},$ (45)

which is positive for every non-zero value of $E_{M}$ and $E_{S}$. This means
that despite non-triviality of the memory’s spectrum, our simple memory-
extended protocol achieves a cooling advantage over Markovian processes.

###  IV.3 Two-level control is sufficient for thermal operations

Figure 10: Convergence to states inaccessible via elementary thermal
operations. Log-log plot of the total variation distance $\delta$ between the
state $\bm{q}$ from Eq. (46) and the state obtained from
$\bm{p}=(1,0,\dots,0)$ via the algorithm
$\widetilde{\mathcal{P}}^{\Pi_{\bm{q}}}$, as a function of the memory size
$N$. Left: systems with energy spectra $E_{i}=i$ with $d=3,\ldots,10$ (bottom
to top) and for $\beta=1.1\log(2)>\beta_{\mathrm{crit}}$. For all presented
dimensions the convergence is better than $1/(2\sqrt{N})$ (black dashed line).
Right: systems with energy spectra $(0,E_{1},E_{2},1)$ taken from a grid with
interval $\Delta=1/64$ (translucent blue lines) for
$\beta=1.1\cdot\beta_{\text{crit}}$. The red lines represent the extreme cases
of convergence, while the thick blue line is the average convergence. Note
that the top red line almost agrees with $1/(2\sqrt{N})$ (black dashed line),
which well approximates the expected convergence for $\beta=0$ and agrees with
the fact that it is obtained for almost completely degenerate levels.

In Ref. [27] it has been proved that there exist thermodynamic state
transformations that cannot be decomposed into the so-called _elementary
thermal operations_ , i.e., thermal operations acting only on two levels of
the system at the same time. Then, in Ref. [34], for any dimension $d$, an
explicit final state $\bm{q}\in C_{+}^{\mathrm{TO}}(\bm{p})$ was given such
that it cannot be achieved (even approximately) starting from the ground state
$\bm{p}=(1,0,\dots,0)$ using convex combinations of sequences of elementary
thermal operations. More precisely, given the energy spectrum of the system
with $E_{i+1}\geq E_{i}$, this final state is given by

$\displaystyle\bm{q}$ $\displaystyle=\left(1-\sum_{i=2}^{d}e^{-\beta
E_{i}},e^{-\beta E_{2}},\dots,e^{-\beta E_{d}}\right)$ (46)

with $\beta\geq\beta_{\text{crit}}$ such that
$1-\sum_{i=2}^{d}e^{-\beta_{\text{crit}}E_{i}}=0$. It was then proven by the
authors of Ref. [34] that there exists $\epsilon>0$ such that any
$\bm{q}^{\prime}$ achievable from $\bm{p}$ satisfies
$\delta(\bm{q},\bm{q}^{\prime})\geq\epsilon$.

Given the above, one might conclude that being able to selectively couple to
the bath just two energy levels at once is highly restrictive and does not
allow one to induce all the transitions possible via general thermal
operations. This conclusion, however, would be incorrect, as the restriction
only arises when one is limited to coupling only two levels of the _system_ at
a given time. As we have seen in this paper, when one is allowed to bring an
auxiliary $N$-level system in a thermal equilibrium state $\gamma_{M}$, then
the ability to selectively couple to the bath just two energy levels of the
joint system allows one to induce all transitions of the main system possible
via thermal operations as $N\to\infty$. Crucially, the operation

${\cal E}(\rho)=\rho\otimes\gamma_{M}$ (47)

is a thermal operation for every $N$. Thus, ${\cal E}$ followed by a sequence
of elementary thermal operations on the joint system, followed by discarding
the system $M$ at the end, can induce any energy-incoherent state transition
of the system possible via general thermal operations. In other words,
elementary control over two energy levels at a given time is sufficient to
generate all thermodynamically possible transitions if we allow ancillary
thermal systems.

Figure 11: Evolution of non-equilibrium free energies and correlations. Non-
equilibrium free energies of the main system [(a) and (e)], the memory [(b)
and (f)], and the joint system [(c) and (g)], as well as the mutual
information between the system and memory [(d) and (h)], as a function of the
step number $n$ of the protocol ${\cal P}^{(ij)}$. Here, the composite system
consists of a three-level system initialised in a state $\bm{p}=(0.7,0.2,0.1)$
and a 16-dimensional degenerate memory prepared in a thermal (maximally mixed)
state, and the plots are presented for two inverse temperatures, $\beta=0$ and
$\beta=0.5$.

We illustrate the above with the following numerical examples, showing that
our MeMTP protocol $\widetilde{\mathcal{P}}^{\Pi_{\bm{q}}}$ (which consists of
only two-level operations) is able to transform $\bm{p}$ into
$\bm{q}^{\prime}$ that approximates $\bm{q}$ arbitrarily well (i.e.,
$\delta(\bm{q}^{\prime},\bm{q})\to 0$ as $N\to\infty$). In order to focus
attention, we chose a constant $\beta=1.1\cdot\beta_{\text{crit}}$. First, we
considered systems of varying dimension $d$, up to $d_{\text{max}}=80$, and
fixed the energy structure to $E_{i}=i$, corresponding to quantum harmonic
oscillator. We observed that the convergence for all these dimensions scales
according to the predictions from Conjecture 2, which can be seen in the left
panel of Fig. 10 for $d=3,\ldots,10$ and memory sizes up to $N=2^{8}$.
Moreover, in order to ascertain that the convergence does not depend on the
energy structure of the system, we fixed $d=4$ and considered energy levels
$(0,E_{1},E_{2},1)$ with $E_{1}<E_{2}$ taken from a grid with spacing $\Delta
E=2^{-6}$, resulting in $1953$ uniformly distributed points. For each of these
points, we have considered the protocol with memory size up to $N=2^{8}$. As
demonstrated in the right panel of Fig. 10, it turns out again that the
convergence, independently from the energy structure, is better than
$1/\sqrt{N}$, in accordance with Conjecture 2.

###  IV.4 Non-equilibrium free energy evolution

To understand how non-Markovian effects arise in the memory-assisted protocol,
we will now examine the evolution of the system and memory during the protocol
$\mathcal{P}^{(ij)}$. More precisely, let us denote the joint state of the
system and memory after the $n$-th two-level thermalisation step of the
protocol by $\bm{p}_{SM}^{(n)}$. Similarly, let $\bm{p}^{(n)}$ and
$\bm{p}_{M}^{(n)}$ denote the reduced states of the system and memory after
the $n$-th step. Then, in the spirit of the analysis performed for elementary
thermal operations in Ref. [28], we will examine the behaviour of the
following entropic quantities. First, we will look at the relative entropy
between $\bm{p}^{(n)}$ and the thermal state of the system $\bm{\gamma}$,

$D\left(\bm{p}^{(n)}\|\bm{\gamma}\right)=\sum^{d}_{i=1}p^{(n)}_{i}\log\frac{p^{(n)}_{i}}{\gamma_{i}},$
(48)

which is a thermodynamic monotone, as it decreases under (Markovian) thermal
operations, and is directly related to the non-equilibrium free energy [49].
We will also look into the behaviour of the analogous quantities for the joint
system and the memory system. Moreover, to track the correlations that build
up between the system and memory, we will investigate the mutual information
between them, which is given by
$D(\bm{p}^{(n)}_{SM}\|\bm{p}^{(n)}\otimes\bm{p}^{(n)}_{M})$.

We use a three-level system and a $16$-dimensional memory as an illustrative
example. We consider the joint system undergoing a $\beta$-swap protocol
$\mathcal{P}^{(ij)}$ for $\beta=0$ and $\beta=0.5$. As shown in Figs. 11a and
11e, the non-equilibrium free energy of the main system initially decreases to
a minimum, and then increases until it reaches a level that closely
approximates the target state (a swap/$\beta$-swap). It is important to note
that this observed increase is only possible because of the presence of the
memory system. In contrast, note that the global non-equilibrium free energy
decreases after each step, as depicted in Figs. 11c and 11g. However, during
the process, a fraction of the main system’s non-equilibrium free energy is
transferred to the memory, which acts as a free energy storage. As such, it
later enables the system to increase its local free energy again, hence
allowing it to achieve the final state. More interestingly, the free energy of
the memory, presented in Figs. 11b and 11f, exhibits a comb-like structure
consisting of $d-1$ teeth with $(d+1)$ steps each. Specifically, within the
$k$-th tooth, the first $d-k$ steps increase the free energy, while the
remaining $k$ steps decrease it. Note that as long as the memory is not
thermalised, its non-equilibrium free energy does not go to zero. However, for
$\beta=0$, it approaches a value very close to zero, but there are still
correlations between the memory and the system, which are illustrated in Figs.
11d and 11e.

##  V Conclusions and outlook

In this paper, we proposed a novel approach to investigate memory effects in
thermodynamics by introducing the concept of memory-assisted Markovian thermal
processes. These were defined by extending the framework of Markovian thermal
processes with ancillary memory systems brought in thermal equilibrium states.
Our construction allowed us to interpolate between the regime of memoryless
dynamics and the one with full control over all degrees of freedom of the
system and the bath. Using a family of protocols composed of Markovian thermal
processes, we demonstrated that energy-incoherent states achievable from a
given initial state via thermal operations could be approximated arbitrarily
well via our algorithmic procedure employing memory. Furthermore, we analysed
the convergence of our protocols in the infinite memory limit, finding
polynomial and exponential convergence rates for infinite and finite
temperatures, respectively. In the infinite temperature limit, we provided
analytic convergence to the entire set of states accessible via thermal
operations. For finite temperatures, we proved the convergence to a subset of
accessible states and, based on extensive numerical evidence, we conjectured
that a modified version of our protocol can realise arbitrary transitions
achievable via thermal operations with an exponential convergence rate that
grows with memory size. Our model-independent approach can be seen as a
significant step forward in understanding ultimate limits of the Markovian
evolution in general, which should be contrasted with the model-specific
approaches to the so-called Markovian embedding [50, 51, 52]. On the other
hand, it may be seen as far less general than the approach taken in Ref. [53],
where our work would correspond to a step towards simulating arbitrary
evolution with Markov-Stinespring curves.

Our framework should be contrasted with previous investigations into the
embeddability of Markov chains into continuous Markov processes via hidden
states [54, 55, 56]. These approaches consider a system with a set of ”visible
states,” which are occupied and subjected to operations, along with ”hidden
states” that enable an embedding of operations that would otherwise be
unfeasible. These hidden states could correspond to unpopulated levels within
the main system—an option that remains viable for implementing our protocols.
Alternatively, we could choose to extend the system by introducing a memory
initialized in a specific energy state $\bm{r}$ with populations
${r_{i}=\delta_{ij}}$ for some $j$. Nevertheless, this approach carries two
caveats. The first involves the cost associated with bringing the non-
equilibrium state $\bm{r}$ into play, as opposed to a Gibbs state. The second
pertains to the catalytic nature of the operation—whether operations of the
form $\bm{p}\otimes\bm{r}\rightarrow\bm{q}\otimes\bm{r}$, where $\bm{r}$ is
retrievable, can indeed expand the space of accessible states $\bm{q}$.

We also explained how our results can be employed to quantitatively assess the
role of memory for the performance of thermodynamic protocols. In this
context, we discussed the dependence on the memory size of the amount and
quality of work extracted from a given non-equilibrium state. In addition, we
introduced a minimal model designed to cool a two-level system below ambient
temperature using a two-dimensional memory. However, the method can be used as
well to investigate other thermodynamic protocols, such as information erasure
or thermodynamically free encoding of information [57]. Furthermore, we
revealed that all transitions accessible via thermal operations can be
accomplished using a restricted set of thermal operations that exclusively
affect only two energy levels (of the system extended by a memory) at any
given time. These findings carry important implications, not only for the
development of efficient thermodynamic protocols, such as optimal cooling and
Landauer erasure, but also for the exploration of novel avenues of research
focused on characterising memory effects in thermodynamics. Finally, we also
commented on the role played by the memory system as a free energy storage
that enables non-Markovian effects.

Our results offer many possibilities for generalisation and further research.
First, one can try proving that the future thermal cone for memory-assisted
Markovian thermal processes agrees with that of thermal operations in the
limit of infinite memory,
$\lim_{N\rightarrow\infty}C^{+}_{\text{MeMTP}}=C^{+}_{\text{TO}}$, as
suggested by Conjecture 2. This can be built upon the proofs for $\beta$-swaps
(Theorem 4) and $\beta$-$3$-cycles (Theorem 5) presented in this work. Second,
one may also attempt to show that the convergence of the proposed protocols
$\mathcal{P}^{\Pi}$ and $\widetilde{\mathcal{P}}^{\Pi}$ is optimal with
respect to the memory size. In other words, one could investigate the upper-
bound on the power of memory-assisted Markovian thermal processes with a given
size of memory $N$. Third, from a more practical point of view, it may be
worthwhile to explore MeMTPs involving finite and infinite memory with non-
trivial energy level structure. The practical relevance of this direction can
be understood by considering the introduction of non-degenerate splitting of
the levels for the full system, which would allow the level pairs to be
addressed independently.

In addition to the above, there are also less clear-cut goals for future
efforts, such as expanding the studies beyond energy-incoherent states into
the full range of quantum states. Furthermore, while our work focused on a
single main system, an interesting avenue for future work could be to
investigate many non-interacting subsystems. This extension could shed light
on the combined consequences of finite-size and memory effects, providing
valuable insights into the behaviour of larger, more complex systems.
Specifically, characterising such effects could help to identify strategies
for improving the efficiency of thermodynamic protocols in practical
applications. Finally, one can also consider memory composed of many
equivalent systems (such as a multi-qubit memory), and analyse the potential
challenges arising from energy-level degeneration in such a setting.

Finally, the feasibility of the introduced algorithm can be studied from a
control perspective, following the approaches outlined in [58, 59]. The first
approach introduces the notion of a space-time trade-off, which refers to the
minimal amount of memory and time steps required to classically implement a
given process. The second approach deals with control complexity, defining it
as the number of levels a given operation non-trivially acts versus the time
steps needed to implement that process. Our algorithm has specific time and
memory requirements, namely $N$-dimensional memory and $N^{2}$ time steps.
Furthermore, it is limited to the simplest two-level processes at any given
time, meaning its control complexity is as low as possible. Nonetheless,
future work might explore variations of our protocol (or any of the variants
presented in Appendix F). Such explorations could focus on enabling
parallelisation of certain steps by expanding available memory, thereby
illustrating the space-time trade-off.

###### Acknowledgements.

KK would like to thank Paul Skrzypczyk for asking a stimulating question at
QIP2022 that started this project. The authors acknowledge financial support
from the Foundation for Polish Science through the TEAM-NET project (contract
no. POIR.04.04.00-00-17C1/18-00). During the process of preparing this
manuscript, we were made aware of the overlap with Ref. [60], and we thank
Jeongrak Son and Nelly H. Y. Ng for insightful discussions on this topic.

##  Appendix A Thermodynamic evolution of energy-incoherent states

This work focuses on energy-incoherent states that can be represented by
$d$-dimensional probability vectors of their eigenvalues, and their evolution
is described by stochastic matrices. As a result, this specificity allows one
to replace density operators and quantum channels with probability vectors and
stochastic matrices, respectively. Consequently, given two states, $\rho$ and
$\sigma$, with eigenvalues $\bm{p}$ and $\bm{q}$, the existence of a thermal
operation between $\rho$ and $\sigma$ is equivalent to the existence of a
Gibbs-preserving stochastic matrix between $\bm{p}$ and $\bm{q}$ [21, 22].
Surprisingly, determining whether a given state $\bm{p}$ can be transformed
into a target state $\bm{q}$ reduces to checking a finite list of conditions
expressed by a partial-order relation between the initial and target state. In
the infinite-temperature limit, these rules are encoded by the majorisation
relation [61], and in the finite temperature case by thermomajorisation [62,
22]. If one further constrains the evolution to be Markovian, then the
aforementioned rules are expressed by the notion of continuous
thermomajorization [24]. In this appendix, we review and summarise well-known
results concerning these partial order relations.

###  A.1 Majorisation

To formulate the solution underlying the thermodynamic evolution of energy-
incoherent states under thermal operation, we first need to recall the concept
of majorisation [61] (see also Ref. [24] for a detailed discussion).

######  Definition 4 (Majorisation).

Given two $d$-dimensional probability distributions $\bm{p}$ and $\bm{q}$, we
say that $\bm{p}$ _majorises_ $\bm{q}$, and denote it by $\bm{p}\succ\bm{q}$,
if and only if the following condition holds:

$\sum_{i=1}^{k}p_{i}^{\downarrow}\geq\sum_{i=1}^{k}q_{i}^{\downarrow}\quad\text{for
all}\quad k\in\\{1\dots d\\},$ (49)

where $\bm{p}^{\downarrow}$ denotes the vector $\bm{p}$ rearranged in a non-
increasing order.

Equivalently, the majorisation relation between distributions can be expressed
in terms of the _majorisation curve_. Given a distribution $\bm{p}$, we define
a piecewise linear curve $f_{\bm{p}}(x)$ in $\mathbb{R}^{2}$. This curve is
obtained by joining the origin $(0,0)$ and the points
$\left(\sum_{i=1}^{k}\eta_{i},\sum_{i=1}^{k}p^{\downarrow}_{i}\right)$, for
each $k$ in the set $\\{1,...,d\\}$. Then, $\bm{p}$ majorises $\bm{q}$ if, and
only if, the majorisation curve $f_{\bm{p}}(x)$ of $\bm{p}$ is always above
that of $\bm{q}$:

$\bm{p}\succ\bm{q}\iff\forall x\in\left[0,1\right]:~{}f_{\bm{p}}(x)\geq
f_{\bm{q}}(x)\,.$ (50)

Majorisation can be interpreted as a formalisation of the notion of disorder
with respect to the uniform distribution $\bm{\eta}$. Note that the uniform
state $\bm{\eta}$ is majorised by any other distribution, while every
distribution is majorised by the sharp state $\bm{s}=(1,0,...,0)$.
Furthermore, if $\bm{p}\succ\bm{q}$, then all Rényi entropies associated with
$\bm{p}$ are smaller than those of $\bm{q}$ [63].

Figure 12:  _Thermomajorisation curve_. For three different states $\bm{p}$,
$\bm{q}$ and $\bm{c}$, with a thermal Gibbs state
$\bm{\gamma}\propto(1,e^{-\beta},e^{-2\beta})$, and $\beta>0$, we plot their
thermomajorisation curves $f^{\,\beta}_{\bm{p}}(x)$ [black curve],
$f^{\,\beta}_{\bm{q}}(x)$ [dashed red curve] and $f^{\,\beta}_{\bm{c}}(x)$
[dot green curve], respectively. While $\bm{p}$ thermomajorises $\bm{q}$
[since $f^{\,\beta}_{\bm{p}}(x)$ is never below $f^{\,\beta}_{\bm{q}}(x)$],
both states are incomparable with $\bm{c}$, as their thermomajorisation curves
cross with $f^{\,\beta}_{\bm{c}}(x)$.

The most general transformation between two probability distributions,
$\bm{p}$ and $\bm{q}$, is accomplished by a stochastic matrix, where $\Lambda$
fulfills $\Lambda_{ij}\geq 0$ and $\sum_{i}\Lambda_{ij}=1$. In the infinite-
temperature limit, thermodynamic transformations are accomplished by
stochastic matrices that preserve the uniform distribution $\bm{\eta}$.
Mathematically, this means that $\Lambda$ is bistochastic, i.e., it
additionally satisfies $\sum_{j}\Lambda_{ij}=1$. As every bistochastic matrix
can be written as a convex combination of permutation matrices, this implies
that the set of $d\times d$ bistochastic matrices is a convex polytope with
$d!$ vertices, one for each permutation in $\mathcal{S}_{d}$. Finally, since
the existence of a bistochastic matrix connecting $\bm{p}$ to $\bm{q}$ is
equivalent to $\bm{p}\succ\bm{q}$, we conclude that the set of states
$C^{\textrm{TO}}_{+}(\bm{p})$ achievable via thermal operations from a given
incoherent initial state $\bm{p}$ is given by

$C^{\textrm{TO}}_{+}(\bm{p})=\operatorname{conv}\left[\left\\{\Pi\bm{p}\,,\mathcal{S}_{d}\ni\pi\mapsto\Pi\right\\}\right],$
(51)

where $\Pi$ denotes a permutation matrix corresponding to the permutation
$\bm{\pi}$ with $d$ elements, and $\operatorname{conv[S]}$ the convex hull of
the set $S$. A detailed discussion about the construction of the future
$C^{\textrm{TO}}_{+}$ and its properties can be found in Sec. III of Ref.
[38].

###  A.2 Thermomajorisation

Thermomajorisation is a measure of disorder relative to the thermal
distribution $\bm{\gamma}$. Mathematically, given $\bm{\gamma}$, we say that
$\bm{p}$ _thermomajorise_ $\bm{q}$, and denote by $\bm{p}\succ_{\beta}\bm{q}$,
if there exists a stochastic matrix $\Lambda^{\bm{\beta}}$, which leaves the
vector $\bm{\gamma}$ invariant and maps $\bm{p}$ onto $\bm{q}$.

The starting point of defining the thermodynamic equivalent of the
majorisation is done by presenting the concept of $\beta$-ordering (see
Definition 2), and then introducing the notion of thermomajorisation curves
[22]. Recall that the $\beta$-ordering of $\bm{p}$ is defined by a permutation
$\bm{\pi}_{\bm{p}}$ that sorts $p_{i}/\gamma_{i}$ in a non-increasing order.
The $\beta$-ordered version of $\bm{p}$, denoted as $\bm{p}^{\beta}$, is
obtained by arranging the elements of $\bm{p}$ according to the permutation
$\bm{\pi}_{\bm{p}}$

$\bm{p}^{\,\beta}=\left(p_{\bm{\pi}_{\bm{p}}^{-1}(1)},\dots,p_{\bm{\pi}_{\bm{p}}^{-1}(d)}\right).$
(52)

A thermomajorisation curve is defined as a piecewise linear curve composed of
linear segments connecting the point $(0,0)$ and the points defined by
consecutive subsums of the $\beta$-ordered form of the probability
$\bm{p}^{\beta}$ and the Gibbs state $\bm{\gamma}^{\beta}$,

$\left(\sum_{i=1}^{k}\gamma^{\,\beta}_{i},~{}\sum_{i=1}^{k}p^{\,\beta}_{i}\right):=\left(\sum_{i=1}^{k}\gamma_{\bm{\pi}^{-1}_{\bm{p}}(i)},~{}\sum_{i=1}^{k}p_{\bm{\pi}^{-1}_{\bm{p}}(i)}\right),$
(53)

for $k\in\\{1,\dots,d\\}$. Thus, given two $d$-dimensional probability
distributions $\bm{p}$ and $\bm{q}$, and a fixed inverse temperature $\beta$,
we say that $\bm{p}$ _thermomajorises_ $\bm{q}$ and denote it as
$\bm{p}\succ_{\beta}\bm{q}$, if the thermomajorisation curve
$f^{\,\beta}_{\bm{p}}$ is above $f^{\,\beta}_{\bm{q}}$ everywhere, i.e.,

$\bm{p}\succ_{\beta}\bm{q}\iff\forall x\in[0,1]:~{}f^{\,\beta}_{\bm{p}}(x)\geq
f^{\,\beta}_{\bm{q}}(x)\,.$ (54)

See Fig. 12 for an example considering a three-level system.

For finite temperatures, general thermodynamic transformations between two
probability distributions, $\bm{p}$ and $\bm{q}$, are accomplished by
stochastic matrices preserving the Gibbs state, i.e., matrices $\Lambda$ such
that $\Lambda\bm{\gamma}=\bm{\gamma}$, which are commonly referred to as
Gibbs-preserving (GP) matrices in the literature. As the rules governing state
transformations are no longer described by a majorisation relation, but
instead by thermomajorisation, the characterisation of the future thermal cone
is no longer given by Eq. (51). Nevertheless, the set of Gibbs-preserving
matrices is still a convex set, and the extreme points of the future thermal
cone can be obtained by the following lemma.

######  Lemma 6 (Lemma 12 of Ref. [27]).

Given $\bm{p}$, consider the following distributions $\bm{p}^{\,\bm{\pi}}\in
C^{\textrm{TO}}_{+}(\bm{p})$ constructed for each permutation
$\bm{\pi}\in\S_{d}$. For $i\in\left\\{1,\dots,d\right\\}$:

  1. 1.

Let $x_{i}^{\bm{\pi}}=\sum_{j=0}^{i}e^{-\beta
E_{\bm{\pi}^{-1}\left(j\right)}}$ and
$y_{i}^{\bm{\pi}}=f^{\,\beta}_{\bm{p}}\left(x_{i}^{\bm{\pi}}\right)$.

  2. 2.

Define
$p^{\bm{\pi}}_{i}:=y^{\bm{\pi}}_{\bm{\pi}(i)}-y^{\bm{\pi}}_{\bm{\pi}(i)-1}$,
with $y_{0}:=0$.

Then, all extreme points of $C^{\textrm{TO}}_{+}(\bm{p})$ have the form
$\bm{p}^{\bm{\pi}}$ for some ${\bm{\pi}}$. In particular, this implies that
$C^{\textrm{TO}}_{+}(\bm{p})$ has at most $d!$ extremal points.

The above lemma allows one to characterise the future thermal cone of $\bm{p}$
by constructing states $\bm{p}^{\bm{\pi}}$ for each
$\bm{\pi}\in\mathcal{S}_{d}$, and taking their convex hull.

###  A.3 Continuous thermomajorisation

If we enforce that the transformation is Markovian, the conditions are
captured by the notion of continuous thermomajorisation [25]. Formally, one
asks if there exists a continuous path within the probability simplex that
connects these two distributions, such that the preceding distribution is
thermomajorised by the succeeding one at any two points along this path. Such
a notion is defined as follows.

######  Definition 5 (Continuous thermomajorisation).

A distribution $\bm{p}$ _continuously thermomajorises_ $\bm{q}$, denoted
$\bm{p}{}_{\beta}\bm{q}$, if there exists a continuous path of probability
distributions $\bm{r}(t)$ for $t\in[0,t_{f})$ such that

  1. 1.

$\bm{r}(0)=\bm{p}$,

  2. 2.

$\forall~{}t_{1},t_{2}\in[0,t_{f}):\quad t_{1}\leq
t_{2}\Rightarrow\bm{r}(t_{1})\succ_{\beta}\bm{r}(t_{2})$,

  3. 3.

$\bm{r}(t_{f})=\bm{q}$.

The path $\bm{r}(t)$ is called _thermomajorising trajectory_ from $\bm{p}$ to
$\bm{q}$.

Let us make a few comments about the above definition. Firstly, it is worth
noting that when $\beta=0$ and the thermal state is replaced by the uniform
fixed point, $\bm{\gamma}=\bm{\eta}$, the above definition corresponds to a
continuous version of standard majorisation, denoted by . Secondly,
determining whether a given initial state continuously thermomajorises a
target state is a difficult problem, and unlike the other variants of
majorisation presented so far, there is no continuous thermomajorisation curve
for this type of majorisation that would facilitate a quick check.
Nevertheless, the necessary and sufficient conditions are known [24]. These
comprise a complete set of entropy production inequalities that can be reduced
to a finitely verifiable set of constraints. Lastly, the notion of continuous
thermomajorisation encapsulates all constrains of memoryless thermal processes
on population dynamics.

##  Appendix B Regularised incomplete beta function

The content of this appendix is based on Refs. [64, 65].

###  B.1 Definition and properties

The finite-size corrections to a $\beta$-swap and its compositions are
determined by the cumulative distribution function (CDF) known as the
regularised incomplete beta function. This function is closely related to the
well-known beta function $B(a,b)$ and is widely used in deriving our results.
To provide the necessary background and present its key properties, we first
recall the definition and properties of the beta function:

$B(a,b)=\int_{0}^{1}t^{a-1}(1-t)^{b-1}dt,$ (55)

with $a,b\in\mathbbm{C}$. The beta function relates to the gamma function in
the following way

$B(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}.$ (56)

The incomplete beta function $B_{x}(a,b)$ is defined by changing the upper
limit of integration in Eq. (55) to an arbitrary variable, i.e.,

$B_{x}(a,b)=\int_{0}^{x}t^{a-1}(1-t)^{b-1}dt.$ (57)

Finally, we define the regularised incomplete beta function $I_{x}(a,b)$
(regularised beta function for short) by normalising the incomplete beta
function,

$I_{x}(a,b)=\frac{B_{x}(a,b)}{B(a,b)}.$ (58)

We present plots of the regularised beta function for a few selected values of
$x$ in Fig. (13).

Throughout this work, we assume that $a,b>0$ and $0\leq x\leq 1$. It is easily
noted that $I_{0}(a,b)=0$, $I_{1}(a,b)=1$, and $I_{0}(a,b)\leq I_{x}(a,b)\leq
I_{1}(a,b)$, thus making it a proper CDF. Furthermore, for
$a,\,b\in\mathbb{Z}$, $I_{x}(a,b)$ can be written in terms of a binomial
function

$I_{x}(a,b)=(1-x)^{b}\sum_{j=a}^{\infty}\binom{b+j-1}{j}x^{j}.$ (59)

From this equation, by using the geometric series and its derivatives, one can
conclude that

$I_{x}(0,b)=1.$ (60)

Figure 13: Regularised beta function. Plots of the regularised incomplete beta
function as a function of (a) $a$ and (b) $b$.

Moreover,

$(1-x)^{b}\sum_{j=a}^{n}\binom{b+j-1}{j}x^{j}=I_{x}(a,b)-I_{x}(n+1,b).$ (61)

If $a=1$, then $I_{x}(1,b)=1-(1-x)^{b}$ and Eq. (59) simplifies to

$(1-x)^{b}\sum_{j=1}^{n}\binom{b+j-1}{j}x^{j}=1-(1-x)^{b}-I_{x}(n+1,b).$ (62)

The next two useful properties of $I_{x}(a,b)$ are the symmetry relation

$I_{x}(a,b)=1-I_{1-x}(b,a),$ (63)

and the relation for an equal argument,

$I_{x}(a,a)=\frac{1}{2}I_{4x(1-x)}\quantity(a,\frac{1}{2}),$ (64)

when $0\leq x\leq 1/2$. Finally, there are two recurrence relations which
allow one to shift either of the arguments of the function by one,

$\displaystyle I_{x}(a,b)$
$\displaystyle=I_{x}(a+1,b)+\frac{x^{a}(1-x)^{b}}{aB(a,b)},$ (65a)
$\displaystyle I_{x}(a,b)$
$\displaystyle=I_{x}(a,b+1)-\frac{x^{a}(1-x)^{b}}{aB(a,b)},$ (65b)
$\displaystyle I_{x}(a,b)$
$\displaystyle=I_{x}(a+1,b-1)+\frac{x^{a}(1-x)^{b-1}}{aB(a,b)},$ (65c)
$\displaystyle I_{x}(a,b)$
$\displaystyle=I_{x}(a-1,b+1)-\frac{x^{a-1}(1-x)^{b}}{aB(a,b)}.$ (65d)

###  B.2 General relations involving the gamma function

Next, let us present general relations and properties of the gamma function
that are extensively used in our proofs. First, recall that for every positive
integer $n$

$\Gamma(n)=(n-1)!.$ (66)

Using the above, one can derive a simple formula for the following expression
that appears when dealing with the regularised beta function:

$\frac{1}{nB(n,n)}=\frac{\Gamma(2n)}{n\Gamma(n)^{2}}=\frac{\Gamma(2n)}{n\Gamma(n)^{2}}\frac{2n}{2n}=\frac{1}{2}\binom{2n}{n}.$
(67)

Other important functional equation for the gamma function is the Legendre
duplication formula:

$\Gamma(n)\Gamma\quantity(n+\frac{1}{2})=2^{1-2n}\sqrt{\pi}\Gamma(2n).$ (68)

The factorial terms can be approximated using Stirling’s approximation

$n!=\sqrt{2\pi n}\quantity(\frac{n}{e})^{n}\left(1+O(n^{-1})\right).$ (69)

###  B.3 Asymptotic analysis

We will be interested in the asymptotic behaviour of the regularised beta
function. So, in this section, we introduce important relations and identities
that will be useful for proving our main theorems. Let us begin by considering
the regularised beta function $I_{x}(a,b)$, where $x$ and $b$ fixed. For
$a\to\infty$, we have the following asymptotic expansion:

$\displaystyle I_{x}(a,b)=$ $\displaystyle\>\>\Gamma(a+b)x^{a}(1-x)^{b-1}$
$\displaystyle\>\>\times\Biggr{[}\sum_{k=0}^{n-1}\frac{1}{\Gamma(a+k+1)\Gamma(b-k)}\quantity(\frac{x}{1-x})^{k}$
$\displaystyle\>+O\quantity(\frac{1}{\Gamma(a+n+1)})\Biggr{]}.$ (70)

Note that $O$-term vanishes in the limit only if $n\geq b$. Furthermore, for
each $n=0,1,2,...$ If $b=1,2,3,...$, and $n>b$, the $O$-term can be omitted,
as the result is exact. In this work, Eq. (B.3) will be expanded up to the
second order. Specifically, for the values of $a=N,b=1/2$, we have the
following equation

$\displaystyle\\!I_{x}\quantity(N,\frac{1}{2})$
$\displaystyle\simeq\\!\frac{x^{N}}{\sqrt{1-x}}\quantity[\frac{\Gamma\quantity(N+\frac{1}{2})}{\Gamma\quantity(N+1)\Gamma\quantity(\frac{1}{2})}+\frac{\Gamma\quantity(N+\frac{1}{2})}{\Gamma\quantity(N+2)\Gamma\quantity(-\frac{1}{2})}\quantity(\frac{x}{1-x})]$
$\displaystyle=\\!\frac{x^{N}}{\sqrt{1-x}}\frac{(2N)!}{4^{N}(N!)^{2}}\quantity[1-\frac{x}{2(N+1)(1-x)}],$
(71)

where $\simeq$ symbol hides the terms of the order $O(x^{N}/N^{2})$.

Next, we will consider sums of the regularised beta functions over the second
argument and their limit as $a\to\infty$. We start with

$\displaystyle\sum_{i=1}^{a}I_{x}(a,i+1)$
$\displaystyle=\sum_{i=1}^{a}\frac{B_{x}(a,i+1)\Gamma(a+1+i)}{\Gamma(a)\Gamma(i+1)}$
$\displaystyle=\sum_{i=1}^{a}\frac{(a+i)!}{(a-1)!i!}\int_{0}^{x}dt\>t^{a-1}(1-t)^{i}$
$\displaystyle=\int_{0}^{x}dt\sum_{i=1}^{a}\binom{a+i}{i}t^{a-1}(1-t)^{i},$
(72)

where we have used definitions of the beta function and the gamma function for
integer arguments. Using Eq. (61), the above expression can be recast as

$\displaystyle\sum_{i=1}^{a}I_{x}(a,i+1)$
$\displaystyle=\int_{0}^{x}\frac{dt}{t^{2}}\>[I_{x}(1,a-1)-I_{1-x}(a-1,a-1)]$
$\displaystyle=\int_{0}^{x}\frac{dt}{t^{2}}\>\quantity[\frac{1}{2}I_{4x(1-x)}\quantity(a-1,\frac{1}{2})-I_{x}(a-1,1)]$
(73)

for $0\leq x\leq 1/2$.

Finally, using the asymptotic expansion from Eq. (B.3), we get

$\displaystyle\int_{0}^{x}dt\>\Biggl{\\{}-t^{a-3}$
$\displaystyle+\frac{1}{2x^{2}}\Biggl{[}\frac{\Gamma(a-1/2)}{\Gamma(a+1)}\frac{[4x(1-x)]^{a-1}}{\sqrt{1-4x(1-x)}}+$
$\displaystyle
O\quantity(\frac{(4x(1-x))^{a-1}\Gamma(a-1/2)}{\Gamma(a+1)})\Biggl{]}\Biggl{\\}}.$
(74)

Thus, we see that this term vanishes in the limit of $a\to\infty$, and we find
that the following sum vanishes in the limit:

$\lim_{a\to\infty}\sum_{i=1}^{a}I_{x}(a,i+1)=0.$ (75)

Moreover, we will need the following sum:

$\displaystyle\sum_{i=1}^{N}$ $\displaystyle
I_{x}(N+1-i,N)=\sum_{i=1}^{N}\frac{B_{x}(N+1-i,N)}{B(N+1-i,N)}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\frac{N\Gamma(2N+1-i)}{N\Gamma(N+1-i)\Gamma(N)}\int_{0}^{x}dt\,t^{N-i}(1-t)^{N-1}\times$
$\displaystyle\times\int_{0}^{x}dtN(1-t)^{N-1}\sum_{i=1}^{N}\binom{2N-i}{N-i}t^{N-i}$
$\displaystyle=$
$\displaystyle\int_{0}^{x}dtN\frac{(1-t)^{N-1}}{(1-t)^{N+1}}(1-t)^{N+1}\sum_{j=0}^{N-1}\binom{(N+1)+j-1}{j}t^{j}$
$\displaystyle=$
$\displaystyle\int_{0}^{x}dt\frac{N}{(1-t)^{2}}\Bigg{(}I_{1-x}(0,N+1)-I_{1-x}(N,N+1)\Bigg{)}$
$\displaystyle=$
$\displaystyle\int_{0}^{x}dt\frac{N}{(1-t)^{2}}\Bigg{(}1-I_{1-x}(N,N)-\frac{t^{N}(1-t)^{N}}{NB(N,N)}\Bigg{)}$
$\displaystyle=$ $\displaystyle
N\int_{0}^{x}dt\frac{1}{(1-t)^{2}}=N\frac{x}{1-x}.$ (76)

##  Appendix C Proofs of Lemma 1 and Theorem 4

To prove Lemma 1 and Theorem 4, we will consider a composite system consisting
of the main $d$-dimensional system and an $N$-dimensional memory system.
Without loss of generality, we can assume that the main system is a two-level
system with $i=1$ and $j=2$, whose state is described by unnormalised
probability vector. We begin by deriving an expression that describes how the
composite system evolves under the memory-assisted protocol
$\mathcal{P}^{(12)}$. Next, to gain insight into the behaviour of the joint
system as the memory size grows, we will prove the asymptotic result. Finally,
in the last subsection, we will show how this result implies the desired
convergence rates.

###  C.1 Dynamics induced by two-level thermalisations

Consider a two-level system, described by a Hamiltonian
$H=E_{1}\outerproduct{E_{1}}{E_{1}}+E_{2}\outerproduct{E_{2}}{E_{2}}$ and
initially prepared in an energy-incoherent state $\bm{p}=(b,c)$, together with
a memory system, described by a trivial Hamiltonian $H_{M}=0$ and prepared in
a maximally mixed state $\bm{\eta}_{M}=\quantity(1,\ldots,1)/N$. The joint
state $\bm{r}:=\bm{p}\otimes\bm{\eta}_{M}$ of the composite system is then
given by

$\bm{r}^{(0)}\equiv\bm{r}=\frac{1}{N}\Big{(}\underbrace{b,\ldots,b}_{N\,\text{times}}|\underbrace{c,\ldots,c}_{N\,\text{times}}\Big{)},$
(77)

whereas the joint thermal distribution
$\bm{\Gamma}:=\bm{\gamma}\otimes\bm{\eta}_{M}$ is given by

$\displaystyle\bm{\Gamma}$ $\displaystyle=\frac{1}{N(e^{-\beta
E_{1}}+e^{-\beta E_{2}})}\quantity[e^{-\beta E_{1}},\cdots,e^{-\beta
E_{1}}|e^{-\beta E_{2}},\cdots,e^{-\beta E_{2}}]$
$\displaystyle=\frac{1}{N}\quantity(\gamma_{1},...,\gamma_{1}|\gamma_{2},...,\gamma_{2}).$
(78)

For the sake of brevity, we introduce the rescaled Gibbs factors, which will
be used extensively in the proofs:

$\Gamma_{ij}=\frac{\gamma_{i}}{\gamma_{i}+\gamma_{j}}.$ (79)

Under a series of two-level thermalisations, the joint state of the composite
system at the $k$-th round is given by

$\\!\\!\\!\mathcal{R}_{k}\quantity(\bm{r}^{(k-1)})\equiv\bm{r}^{(k)}\\!=\\!\quantity(b_{1}^{(k)},\ldots,b_{N}^{(k)}\,|\,c_{1}^{(N)},\ldots,c^{(N)}_{k},c,\ldots,c),$
(80)

where $b^{(k)}_{j}$ and $c^{(N)}_{j}$ satisfy the following recurrence
relations

$\displaystyle b^{(k)}_{j}$
$\displaystyle=\Gamma^{j}_{21}\,\quantity(\frac{\Gamma_{12}}{\Gamma_{21}}c+\Gamma_{12}\sum_{i=1}^{j}b^{(k-1)}_{i}\Gamma^{\,-i}_{21}),$
(81) $\displaystyle c^{(N)}_{j}$
$\displaystyle=\Gamma_{21}^{\,N}c+\Gamma_{21}^{\,N+1}\sum_{i=1}^{N}b^{(j-1)}_{i}\Gamma^{\,-i}_{21},$
(82)

with $b^{(0)}_{j}=b$. Equations (81) and (82) can be understood by noting that
during the $k$-th round of two-level thermalisations, an additional $c$ is
added to the previous $(k-1)$-th entry and the resulting state is again
thermalised. By iterating this process for the first $k$ rounds, we can derive
a closed-form expression for the entries $b^{(k)}_{j}$ and $c^{(N)}_{j}$:

$\displaystyle\\!\\!\\!b_{j}^{(k)}$
$\displaystyle=\Gamma_{12}\Gamma_{21}^{j-1}c\sum_{i=0}^{k-1}\binom{j\\!+\\!i\\!-\\!1}{i}\Gamma_{12}^{i}+b\Gamma_{12}^{k}\sum_{i=1}^{j}\binom{j\\!+\\!k\\!-\\!1\\!-\\!i}{k\\!-\\!1}\Gamma_{21}^{j-i},$
(83) $\displaystyle\\!\\!\\!c^{(N)}_{j}$
$\displaystyle=c\Gamma^{N}_{21}\sum_{i=0}^{j-1}\binom{N\\!+\\!i\\!-\\!1}{i}\Gamma_{12}^{i}+b\Gamma_{21}\Gamma^{j}_{12}\sum_{i=0}^{N-1}\binom{i\\!+\\!j}{j}\Gamma_{21}^{i}.$
(84)

After $N^{2}$ rounds of two-level thermalisations, the composite final state
$\bm{r}^{(N)}$ is given by

$\displaystyle\widetilde{\mathcal{P}}^{(12)}(\bm{r})$
$\displaystyle\equiv\bm{r}^{(N)}=\bm{q}\otimes\bm{\eta}_{M}$
$\displaystyle=\frac{1}{N}{\small\quantity[b^{(N)}_{1},\ldots,b^{(N)}_{N}\,\bigg{|}\,c^{(N)}_{1},\ldots,c^{(N)}_{N}]}.$
(85)

###  C.2 Infinite memory limit

Up till now, we have obtained a closed-form expression describing the action
of the truncated protocol $\widetilde{\mathcal{P}}^{(12)}$. In order to prove
Lemma 1 and Theorem 4, the next step consists of decoupling the main system
from memory using a full thermalisation $\mathcal{T}$:

$\mathcal{P}^{(12)}(\bm{p}\otimes\bm{\gamma}_{M})=\mathcal{T}\circ\widetilde{\mathcal{P}}^{(12)}(\bm{p}\otimes\bm{\gamma}_{M}).$
(86)

Since we are initially focused on demonstrating the asymptotic results, our
aim is to show that as $N$ approaches infinity, we achieve a $\beta$-swap:

$\displaystyle\lim_{N\to\infty}\frac{\sum_{i=1}^{N}b^{(N)}_{i}}{N}$
$\displaystyle=c+b(1-e^{-\beta(E_{2}-E_{1})}),$ (87a)
$\displaystyle\lim_{N\to\infty}\,\frac{\sum_{i=1}^{N}c_{i}^{(N)}}{N}$
$\displaystyle=be^{-\beta(E_{2}-E_{1})}.$ (87b)

To prove the limits in Eqs. (87a) and (87b), we will begin by using the
conservation of probability,

$\lim_{N\to\infty}\sum_{i=1}^{N}\frac{b^{(N)}_{i}}{N}+\lim_{N\to\infty}\sum_{i=1}^{N}\frac{c_{i}^{(N)}}{N}=b+c.$
(88)

This allows us to express the first limit as a function of the second:

$\lim_{N\to\infty}\sum_{i=1}^{N}\frac{b^{(N)}_{i}}{N}=b+c-\lim_{N\to\infty}\sum_{i=1}^{N}\frac{c^{(N)}_{i}}{N}.$
(89)

We will now examine the non-trivial term on the right hand side of Eq. (89).
We start by using Eq. (61) to re-write the first term of $c^{(N)}_{j}$
appearing in Eq. (84) in terms of the regularised beta function

$\displaystyle\evaluated{c^{(N)}_{j}}_{\begin{subarray}{c}b=0\\\
c=1\end{subarray}}=\Gamma^{N}_{21}\sum_{i=0}^{j-1}\binom{N+i-1}{i}\Gamma_{12}^{i}$
$\displaystyle=[I_{\Gamma_{12}}(0,N)-I_{\Gamma_{12}}(j,N)]$
$\displaystyle=I_{\Gamma_{21}}(N,j),$ (90)

where we also used Eq. (63) to invert the arguments of the regularised beta
function. Thus, using (75), we conclude that

$\lim_{N\to\infty}\frac{1}{N}\sum_{i=1}^{N}\evaluated{c^{(N)}_{j}}_{\begin{subarray}{c}b=0\\\
c=1\end{subarray}}=\lim_{N\to\infty}\frac{1}{N}\sum_{i=1}^{N}I_{\Gamma_{21}}(N,i)=0.$
(91)

Taking into account the second term of $c^{(N)}_{j}$ appearing in Eq. (84),
one can immediately evaluate the sum

$\displaystyle\lim_{N\rightarrow\infty}\evaluated{c^{(N)}_{j}}_{\begin{subarray}{c}b=1\\\
c=0\end{subarray}}=\Gamma_{21}\Gamma^{j}_{12}\sum_{i=0}^{\infty}\binom{i+k}{k}\Gamma_{21}^{i}=\frac{\Gamma_{21}}{\Gamma_{12}},$
(92)

and therefore obtain that

$\lim_{N\to\infty}\,\frac{\sum_{i=1}^{N}c^{(N)}_{i}}{N}=be^{-\beta(E_{2}-E_{1})}.$
(93)

Substituting this result to Eq. (89), we prove that

$\lim_{N\to\infty}\frac{\sum_{i=1}^{N}b^{(N)}_{i}}{N}=c+b[1-e^{-\beta(E_{2}-E_{1})}].$
(94)

As a result, in the limit of $N\to\infty$, the protocol $\mathcal{P}^{(12)}$
achieves a $\beta$-swap.

###  C.3 Finite memory convergence rates

To complete the proofs of Lemma 1 and Theorem 4, we will now analyse what the
approximation error is for a finite size of the memory $N$. This will tell us
how quickly the initial state convergences to a $\beta$-swap as a function of
$N$. First, we define two functions governing the convergence

$\displaystyle\\!\\!\\!\mathbb{E}:=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\evaluated{c^{(N)}_{i}}_{c=1,b=0}\\!\\!=\frac{\Gamma^{N}_{21}}{N}\sum_{i=1}^{N}\sum_{j=0}^{i-1}\binom{N\\!+\\!j\\!-\\!1}{j}\Gamma^{j}_{12},$
(95a) $\displaystyle\\!\\!\\!\mathbb{F}:=$
$\displaystyle\frac{1}{N}\sum_{j=1}^{N}\evaluated{b_{j}^{(N)}}_{c=0,b=1}\\!\\!=\frac{\Gamma^{N}_{12}}{N}\sum_{j=1}^{N}\sum_{i=1}^{j}\binom{j\\!+\\!N\\!-\\!1\\!-\\!i}{N\\!-\\!1}\Gamma^{j-i}_{21},\\!$
(95b)

which allow us to write the final state of the system as

$\bm{q}=b\matrixquantity(\mathbb{F}\\\
1-\mathbb{F})+c\matrixquantity(1-\mathbb{E}\\\ \mathbb{E}).$ (96)

Next, we will asymptotically expand Eq. (95a). The starting point is to reduce
the double sum into one, then convert the binomial sums into regularised beta
functions and use its properties given by Eqs. (65c), (64) and (67):

$\displaystyle\mathbb{E}$
$\displaystyle=\frac{\Gamma^{N}_{21}}{N}\sum_{i=1}^{N}(N-i+1)\binom{N+i-2}{i-1}\Gamma_{12}^{i-1}$
$\displaystyle=\Gamma^{N}_{21}\sum_{i=0}^{N-1}\binom{N+i-1}{i}\Gamma_{12}^{i}-\frac{\Gamma^{N}}{N}\sum_{i=1}^{N-1}i\>\binom{N+i-1}{i}\Gamma_{12}^{i}$
$\displaystyle=I_{\Gamma_{21}}(N,N)-\frac{\Gamma_{12}}{\Gamma_{21}}I_{\Gamma_{21}}(N+1,N-1)$
$\displaystyle=\quantity(1-\frac{\Gamma_{12}}{\Gamma_{21}})I_{\Gamma_{21}}(N,N)+\frac{\Gamma_{12}}{\Gamma_{21}}\frac{\Gamma_{21}^{N}\Gamma_{12}^{N-1}}{NB(N,N)}$
$\displaystyle=\frac{1}{\Gamma_{21}}\quantity[(\Gamma_{21}-\Gamma_{12})\frac{1}{2}I_{4\Gamma_{21}\Gamma_{12}}\quantity(N,\frac{1}{2})+\frac{(\Gamma_{21}\Gamma_{12})^{N}}{N}\frac{\Gamma(2N)}{\Gamma(N)^{2}}].$
(97)

Now we expand Eq. (97) up to second order using Eq. (B.3). Recall that such an
expansion is an approximation up to terms of the order $O(x^{N}/N^{2})$ (where
$x$ will be given by $4\Gamma_{21}\Gamma_{12}$), which will be dropped since
our final approximation will be up to the order $o(x^{N}/N^{3/2})$ or
$o(1/N^{1/2})$ for finite and infinite temperatures, respectively. Then, the
gamma functions appearing in the expansions are simplified using Eq. (68) and
Eq. (67). Finally, using Stirling’s approximation [Eq. (69)], we arrive at

$\displaystyle\mathbbm{E}=$
$\displaystyle\frac{(4\Gamma_{21}\Gamma_{12})^{N}}{2\Gamma_{21}\sqrt{\pi
N}}\Bigg{[}\frac{(\Gamma_{21}-\Gamma_{12})}{\sqrt{1-4\Gamma_{21}\Gamma_{12}}}-\frac{(\Gamma_{21}-\Gamma_{12})(4\Gamma_{21}\Gamma_{12})}{2(1-4\Gamma_{21}\Gamma_{12})^{3/2}(N+1)}$
$\displaystyle\quad\qquad\qquad+1+o\left(N^{-1}\right)\Bigg{]}[1+O(N^{-1})].$
(98)

As the last step, we simplify the above expression for $\Gamma_{12}>1/2$
(finite temperature case) by using the fact that $\Gamma_{12}=1-\Gamma_{21}$,
to arrive at

$\\!\\!\\!\\!\\!\\!\evaluated{\mathbb{E}}_{\Gamma_{12}>\frac{1}{2}}\\!\\!\\!\\!\\!\\!\\!=(4\Gamma_{12}\Gamma_{21})^{N}\left[\frac{\Gamma_{12}}{(\Gamma_{12}\\!-\\!\Gamma_{21})^{2}\sqrt{\pi
N}(N\\!+\\!1)}+o\left(N^{-3/2}\right)\right].$ (99)

The infinite-temperature limit, $\beta=0$, is similarly analysed by using Eq.
(97) and plugging $\Gamma_{12}=\Gamma_{21}=1/2$. In this case, the scaling is
slightly different and is given by

$\evaluated{\mathbb{E}}_{\Gamma_{12}=\frac{1}{2}}=\frac{1}{\sqrt{\pi
N}}+o\left(N^{-1/2}\right).$ (100)

The other half of estimating the convergence rate stems from considering Eq.
(95b). Proceeding in the same manner as before, we re-write it as

$\displaystyle\mathbb{F}$
$\displaystyle=\frac{\Gamma^{N}_{12}}{N}\sum_{k=1}^{N}(N-k+1)\binom{N+k-2}{N-1}\Gamma_{21}^{k-1}$
$\displaystyle=\Gamma^{N}_{12}\sum_{k=0}^{N-1}\binom{N+k-1}{k}\Gamma_{21}^{k}-\frac{\Gamma^{N}_{12}}{N}\sum_{k=1}^{N-1}k\>\binom{N+k-1}{k}\Gamma_{21}^{k}$
$\displaystyle=[1-I_{\Gamma_{21}}(N,N)]-\frac{\Gamma_{21}}{\Gamma_{12}}[1-I_{\Gamma_{21}}(N-1,N+1)]$
$\displaystyle=[1-I_{\Gamma_{21}}(N,N)]-\frac{\Gamma_{21}}{\Gamma_{12}}\quantity[1-I_{\Gamma_{21}}(N,N)-\frac{(\Gamma_{21}\Gamma_{12})^{N}}{\Gamma_{21}NB(N,N)}]$
$\displaystyle=\frac{(\Gamma_{12}-\Gamma_{21})}{\Gamma_{12}}\quantity[1-I_{4\Gamma_{12}\Gamma_{12}}\quantity(N,\frac{1}{2})]+\frac{(4\Gamma_{21}\Gamma_{12})^{N}}{2\Gamma_{12}}\frac{1}{\sqrt{\pi
N}}.$ (101)

Again, we expand Eq. (101) up to second order using Eq. (B.3). Then, the gamma
functions appearing in the expansions are simplified using Eq. (68), the
remaining ones are simplified by using Eq. (67). Finally, using Stirling’s
approximation [Eq. (69)], we arrive at

$\displaystyle\mathbbm{F}=\>$
$\displaystyle\quantity(\frac{\Gamma_{12}-\Gamma_{21}}{2\Gamma_{12}})\Biggl{[}2-\frac{(4\Gamma_{21}\Gamma_{12})^{N}}{\sqrt{1-4\Gamma_{21}\Gamma_{12}}}\frac{1}{\sqrt{\pi
N}}$
$\displaystyle\\!\\!\times\quantity(1-\frac{4\Gamma_{21}\Gamma_{12}}{2(N+1)(1-4\Gamma_{21}\Gamma_{12})})\Biggr{]}$
$\displaystyle+\frac{(4\Gamma_{21}\Gamma_{12})^{N}}{2\Gamma_{12}}\left[\frac{1}{\sqrt{\pi
N}}+o(N^{-3/2})\right].$ (102)

For $\Gamma_{12}>1/2$ (finite temperature case), we use the fact that
$\Gamma_{12}=1-\Gamma_{21}$, to simplify the above as

$\mathbb{F}\simeq\underbrace{\frac{\Gamma_{12}-\Gamma_{21}}{\Gamma_{12}}}_{=1-e^{-\beta(E_{2}-E_{1})}}+\underbrace{\Gamma_{21}\frac{(4\Gamma_{12}\Gamma_{21})^{N}}{(N+1)\sqrt{\pi
N}\quantity(\Gamma_{12}-\Gamma_{21})^{2}}}_{\equiv\mathbb{G}},$ (103)

where $\simeq$ hides the $o$-terms. The infinite-temperature limit, $\beta=0$,
is obtained by using Eq. (101) and plugging $\Gamma_{12}=\Gamma_{21}=1/2$.
This yields the following convergence

$\evaluated{\mathbb{F}}_{\Gamma_{12}=\frac{1}{2}}=\frac{1}{\sqrt{\pi
N}}+o\left(N^{-1/2}\right)$ (104)

with the expression for $\mathbb{G}$ modified to $(\pi N)^{-1/2}$.

As a final step to prove Lemma 1 and Theorem 4, we calculate explicitly the
state of the primary system after the protocol $\mathcal{P}^{(12)}$. This is
done by substituting the results for $\mathbb{E}$ and $\mathbb{F}$ to Eq.
(96), yielding

(105)

Thus, the distance between $\bm{q}$ and the target $\Pi^{\beta}_{12}\bm{p}$ is
given by

$\displaystyle\delta\quantity(\Pi^{\beta}_{12}\bm{p},\,\bm{q})$
$\displaystyle=\absolutevalue{b\mathbb{G}-c\mathbb{E}}.$ (106)

For $\beta=0$ case, the above gives

$\displaystyle\delta\quantity(\Pi^{\beta}_{12}\bm{p},\,\bm{q})$
$\displaystyle=\frac{\absolutevalue{b-c}}{\sqrt{\pi
N}}+o\left(N^{-1/2}\right),$ (107)

whereas $\beta\neq 0$ case, it gives

$\displaystyle\delta\quantity(\Pi^{\beta}_{12}\bm{p},\,\bm{q})$
$\displaystyle=(4\Gamma_{12}\Gamma_{21})^{N}\Bigg{[}\frac{\absolutevalue{b\Gamma_{21}-c\Gamma_{12}}}{(\Gamma_{12}-\Gamma_{21})^{2}\sqrt{\pi
N}(N+1)}$ $\displaystyle\qquad\qquad\qquad+o(N^{-3/2})\Bigg{]}.$ (108)

These prove Lemma 1 and Theorem 4 after going to the notation used therein,
i.e., $b\rightarrow p_{1},\,c\rightarrow
p_{2},\,\Gamma_{12}\rightarrow\Gamma_{1}$ and
$\Gamma_{21}\rightarrow\Gamma_{2}$. ∎

##  Appendix D Strengthening Corollary 3

Corollary 3 deals with a very general approach to bounding the distance
between the target state $\bm{p}^{\pi}$ and its approximation obtained from
$\bm{p}$ via the MeMTP protocol $\mathcal{P}^{\Pi}$. However, it can be
improved by taking into account the set of indices on which the permutation
acts.

######  Corollary 7.

Consider states $\bm{p}$ and $\bm{q}\in C_{+}^{TO}(\bm{p})$. Then, in the
infinite temperature limit, $\beta=0$, and for an $N$-dimensional memory,
there exists a MeMTP protocol $\mathcal{P}$ such that

${\cal P}(\bm{p}\otimes\bm{\eta}_{M})=\bm{q}^{\prime}\otimes\bm{\eta}_{M},$
(109)

with

$\delta(\bm{q}^{\prime},\bm{q})\leq\frac{1}{2\sqrt{\pi
N}}\sum_{\begin{subarray}{c}k,l\\\ k\neq
l\end{subarray}}\absolutevalue{p_{i_{k}}-p_{i_{l}}}+o\quantity(N^{-1/2})=:\epsilon,$
(110)

where $\quantity{i_{1},\ldots,i_{d^{\prime}}}\subset\quantity{1,\ldots,d}$ is
a subset of indices neighbouring in the $\beta$-order,
$\pi_{\bm{p}}(i_{j})+1=\pi_{\bm{p}}(i_{j+1})$, such that
$\Pi_{\bm{q}}=\prod_{i=1}\Pi_{j_{i}k_{i}}$ with
$j_{i},k_{i}\in\quantity{i_{1},\ldots,i_{d^{\prime}}}$.

###### Proof.

First, we consider a target state to be an extreme point $\bm{p}^{\pi^{*}}$
such that the $\beta$-order $\Pi_{\bm{p}^{*}}\equiv\Pi^{*}$ can be decomposed
into the maximal number of $d^{\prime}(d^{\prime}-1)/2$ neighbour swaps on the
levels $i_{1}$ through $i_{d^{\prime}}$. Taking explicitly Eq. (25) from
Theorem 2, one finds that

$\delta\left(\bm{p}^{\pi^{*}},\bm{r}^{*}\right)=\frac{1}{2\sqrt{\pi
N}}\sum_{\begin{subarray}{c}k,l\\\ k\neq
l\end{subarray}}\absolutevalue{p_{i_{k}}-p_{i_{l}}}+o\quantity(N^{-1/2}),$
(111)

where for convenience we used
$\mathcal{P}^{\Pi^{*}}(\bm{p}\otimes\eta_{M})=\bm{r}^{*}\otimes\eta_{M}$. We
note that the above expression in fact provides a general upper bound for any
permutation $\Pi$ on the aforementioned subset of $d^{\prime}$ levels –
defining $\mathcal{P}^{\Pi}(\bm{p}\otimes\eta_{M})=\bm{r}\otimes\eta_{M}$ we
find that

$\delta\left(\Pi\bm{p},\bm{r}\right)\leq\epsilon.$ (112)

Now, there are two cases to be considered. First, take a state $\bm{q}$ which
is in the future of the approximation point $\bm{r}$, $\bm{q}\in
C_{+}^{MTP}\quantity(\bm{r})$, from which it follows that

$\exists\mathcal{O}\in\text{MTP}:\delta\left(\bm{q},\mathcal{O}(\bm{r})\right)=0.$
(113)

Otherwise, $\bm{q}$ is not in the future cone of $\bm{r}$. In this case, we
first note that there exists a ball $B(\bm{r},\epsilon^{\prime})\ni\Pi\bm{p}$
with radius $\epsilon^{\prime}\leq\epsilon$ with respect to
$\delta(\cdot,\cdot)$. Thanks to the planarity of the boundaries $\partial
C_{+}^{MTP}\quantity(\bm{r})$ and $\partial C_{+}^{MTP}\quantity(\Pi\bm{p})$
when restricted to a fixed $\beta$-order, we can consider the extreme case

$\bm{q}\in\partial
C_{+}^{MTP}\quantity(\Pi\bm{p})\\!\Rightarrow\\!\exists\bm{r}^{\prime}\\!\in\\!\partial
C_{+}^{MTP}\quantity(\bm{r}):\delta\quantity(\bm{q},\bm{r}^{\prime})\\!\leq\\!\epsilon^{\prime}\\!\leq\\!\epsilon$
(114)

and the same argument applies for any $\bm{q}\in
C_{+}^{MTP}\quantity(\Pi\bm{p})\backslash C_{+}^{MTP}\quantity(\bm{r})$, thus
concluding the proof. ∎

The bound presented in Corollary 7 can be further improved by taking into
account the possibility of dividing the set
$\quantity{i_{1},\ldots,i_{d^{\prime}}}$ into subsets that are not mixed at
any step when considering the decomposition of $\Pi_{\bm{q}}$ into neighbour
transpositions. Finally, we point out that, in agreement with Eq. (113), there
will exist such states $\bm{q}$ that are attainable exactly, and moreover,
their volume will increase together with the size of memory $N$.

##  Appendix E Proof of Theorem 5

To prove Theorem 5, we will consider a composite system consisting of the main
$d$-dimensional system and an $N$-dimensional memory system. Without loss of
generality, we can simply assume that the main system has three levels with
$i_{1}=1,i_{2}=2$ and $i_{3}=3$, and its state is described by an unnormalised
probability vector $\bm{p}=(a,b,c)$. The Hamiltonian is then given by
$H=\sum_{i=1}^{3}E_{i}\outerproduct{E_{i}}{E_{i}}$, while the memory system is
described by a trivial Hamiltonian $H_{M}=0$ and prepared in a maximally mixed
state $\bm{\eta}_{M}=(1/N,...,1/N)$. The joint state of the composite system,
$\bm{r}:=\bm{p}\otimes\bm{\eta}_{M}$, is then given by

$\bm{r}^{(0)}\equiv\bm{r}=\frac{1}{N}\Big{(}\underbrace{a,...,a}_{N\text{
times}}|\underbrace{b,...,b}_{N\text{ times}}|\underbrace{c,...,c}_{N\text{
times}}\Big{)},$ (115)

and the joint thermal state is given by

$\displaystyle\bm{\Gamma}$ $\displaystyle=\frac{1}{ZN}[e^{-\beta
E_{1}},...,e^{-\beta E_{1}}|e^{-\beta E_{2}},...,e^{-\beta E_{2}}|e^{-\beta
E_{3}},...,e^{-\beta E_{3}}]$
$\displaystyle=\frac{1}{N}\quantity(\gamma_{1},...,\gamma_{1}|\gamma_{2},...,\gamma_{2}|\gamma_{3},...,\gamma_{3}),$
(116)

where $Z=\sum_{i=1}^{3}e^{-\beta E_{i}}$.

As before, the starting point consists of understanding how the joint state of
the composite system changes under the action of the composite protocol
$\widetilde{\mathcal{P}}^{(13)}_{N}\circ\widetilde{\mathcal{P}}^{(23)}_{N}$,
whose action is summarised in two steps:

  1. 1.

Two-level thermalisation between second and third energy levels.

  2. 2.

Two-level thermalisation between first and third energy levels.

The final state $\bm{r}^{(N)}$ is then given by

$\displaystyle\\!\\!\\!\widetilde{{\cal P}}^{\Pi}(\bm{r})$
$\displaystyle\equiv\bm{r}^{(N)}=\bm{q}\otimes\bm{\eta}_{M}$
$\displaystyle=\frac{1}{N}{\small\quantity[a^{(N)}_{1},\ldots,a^{(N)}_{N}\,\bigg{|}\,b^{(N)}_{1},\ldots,b^{(N)}_{N}\,\bigg{|}\,c^{(N)}_{1},\ldots,c^{(N)}_{N}]},$
(117)

where $\Pi=\Pi_{\Gamma_{23}}\Pi_{\Gamma_{23}}$. Note that due to probability
conservation, characterising the second and third entries of Eq. (117) is
sufficient.

After the first protocol $\widetilde{\mathcal{P}}^{23}$, the second energy
level remains “untouched” and, as a result, its entries are given by Eq. (83)
(with $\Gamma_{12}$ and $\Gamma_{21}$ replaced by $\Gamma_{23}$ and
$\Gamma_{32}$, respectively, as defined in (79)). The other two entries are
obtained in a similar way as Eqs. (81)-(82). The action of the protocol
generates a recurrence formula that allows us to write the last entry
$c^{(N)}_{k}$ as

$\displaystyle\\!\\!\\!\\!c^{(N)}_{k}=a\Gamma^{k-1}_{13}\sum_{i=1}^{N}$
$\displaystyle\binom{N+k-1-i}{k-1}\Gamma_{31}^{N+1-i}$
$\displaystyle\quad\quad\quad\quad+\Gamma^{N}_{31}\sum_{l=0}^{N-1}\Gamma_{13}^{l}\binom{N+l-1}{l}c_{k-l},$
(118)

where $c_{k}$ is given by

$c_{k}=c\Gamma^{N}_{32}\sum_{i=0}^{k-1}\binom{N+i-1}{i}\Gamma^{i}_{23}+b\Gamma_{32}\Gamma^{k}_{23}\sum_{i=0}^{N-1}\binom{i+k}{k}\Gamma_{32}^{i}.$
(119)

Since, without loss of generality, we assumed that $\bm{p}$ has
$\beta$-ordering $(123)$, the proof boils down to demonstrating that
$\widetilde{{\cal P}}^{\Pi}(\bm{r})$ sends $\bm{p}$ to the following extreme
point

$\bm{p}^{(321)}=\quantity[a+\frac{\Gamma_{32}}{\Gamma_{23}}b-a\frac{\Gamma_{31}}{\Gamma_{13}},c+b\left(1-\frac{\Gamma_{32}}{\Gamma_{23}}\right),\frac{\Gamma_{31}}{\Gamma_{13}}a].$
(120)

Therefore, we need to prove the following limits

$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{i=1}^{N}b^{(N)}_{i}$
$\displaystyle=c+b\,\quantity(1-\frac{\Gamma_{32}}{\Gamma_{23}}),$ (121a)
$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{i=1}^{N}c^{(N)}_{i}$
$\displaystyle=a\frac{\Gamma_{31}}{\Gamma_{13}}.$ (121b)

###  E.1 Proof of limit (121a)

We start by recalling that $b^{(N)}_{j}$ is given by:

$b_{j}^{(N)}\\!=c\\!\frac{\Gamma_{23}}{\Gamma_{32}}\Gamma_{32}^{j}\sum_{i=0}^{N-1}\binom{j\\!+\\!i\\!-\\!1}{i}\Gamma_{23}^{i}+b\Gamma_{23}^{N}\sum_{i=1}^{j}\binom{j\\!+\\!N\\!-\\!1\\!-\\!i}{N\\!-\\!1}\Gamma_{32}^{j-i}.$
(122)

Comparing Eqs. (122) and the right-hand side of (121a), we see that in order
to prove Eq. (121a), we need to prove the following two limits:

$\displaystyle\lim_{N\to\infty}\sum_{j=1}^{N}\evaluated{\frac{b^{(N)}_{j}}{N}}_{\begin{subarray}{c}b=0\\\
c=1\end{subarray}}$
$\displaystyle=\lim_{N\rightarrow\infty}\frac{\Gamma_{23}}{\Gamma_{32}}\frac{1}{N}\sum_{j=1}^{N}\Gamma_{32}^{j}\sum_{i=0}^{N-1}\binom{j+i-1}{i}\Gamma_{23}^{i}$
$\displaystyle=1$ (123)

and

$\displaystyle\lim_{N\to\infty}\sum_{j=1}^{N}\evaluated{\frac{b^{(N)}_{j}}{N}}_{\begin{subarray}{c}b=1\\\
c=0\end{subarray}}$
$\displaystyle=\lim_{N\rightarrow\infty}\Gamma_{23}^{N}\frac{1}{N}\sum_{j=1}^{N}\Gamma_{32}^{j}\sum_{i=1}^{j}\binom{j+N-1-i}{N-1}\Gamma_{32}^{-i}$
$\displaystyle=\quantity(1-\frac{\Gamma_{32}}{\Gamma_{23}}).$ (124)

We begin by proving Eq. (E.1). First, we rewrite this expression as:

$\displaystyle\sum_{j=1}^{N}\evaluated{\frac{b^{(N)}_{j}}{N}}_{\begin{subarray}{c}b=1\\\
c=0\end{subarray}}$
$\displaystyle=\frac{1}{N}\frac{\Gamma_{23}}{\Gamma_{32}}\sum_{i=0}^{N-1}\Gamma_{23}^{i}\sum_{j=1}^{N}\Gamma_{32}^{j}\binom{j+i-1}{i}$
$\displaystyle=\frac{1}{N}\Gamma_{23}\sum_{i=0}^{N-1}\Gamma_{23}^{i}\sum_{j=0}^{N-1}\Gamma_{32}^{j}\binom{j+i}{i}.$
(125)

We can evaluate the second sum in Eq. (125) as follows:

$\displaystyle\sum_{j=0}^{N-1}\Gamma_{32}^{j}\binom{j+i}{i}$
$\displaystyle=(\Gamma_{23})^{-i-1}\Bigg{(}1-\frac{B_{\Gamma_{32}}(N,i+1)}{B(N,i+1)}\Bigg{)}$
$\displaystyle=(\Gamma_{23})^{-i-1}[1-I_{\Gamma_{23}}(N,i+1)].$ (126)

Thus, substituting Eq. (E.1) into Eq. (125), we obtain

$\displaystyle\sum_{j=1}^{N}\evaluated{\frac{b^{(N)}_{j}}{N}}_{\begin{subarray}{c}b=1\\\
c=0\end{subarray}}=1-\frac{1}{N}\sum_{i=0}^{N-1}I_{\Gamma_{32}}(N,i+1).$ (127)

Using Eq. (75), we conclude that the second term in Eq. (127) vanishes in the
limit of $N\rightarrow\infty$, and therefore

$\lim_{n\to\infty}\sum_{j=1}^{N}\evaluated{b^{(N)}_{j}}_{\begin{subarray}{c}b=0\\\
c=1\end{subarray}}=1,$ (128)

so that we have proved Eq. (E.1).

To prove Eq. (E.1), we begin by manipulating it so that we can express it in a
simpler form:

$\displaystyle\sum_{j=1}^{N}\evaluated{\frac{b^{(N)}_{j}}{N}}_{\begin{subarray}{c}b=1\\\
c=0\end{subarray}}\\!\\!\\!\\!$
$\displaystyle=\frac{\Gamma_{23}^{N}}{N}\sum_{j=1}^{N}\Gamma_{32}^{j}\sum_{i=1}^{j}\binom{N+j-1-i}{N-1}\Gamma_{32}^{-i}$
$\displaystyle=\frac{\Gamma_{23}^{N}}{N}\sum_{j=0}^{N-1}(N-j)\binom{N+j-1}{j}\Gamma_{32}^{j}$
$\displaystyle=\Gamma_{23}^{N}\sum_{j=0}^{N-1}\binom{N\\!+\\!j\\!-\\!1}{j}\Gamma_{32}^{j}-\frac{\Gamma_{23}^{N}}{N}\sum_{j=0}^{N-1}j\,\binom{N\\!+\\!j\\!-\\!1}{j}\Gamma_{32}^{j}.$
(129)

Applying Eq. (61) to transform the first term of Eq. (E.1) into a difference
of regularised beta functions, and then using its asymptotic expansion, we
obtain

$\Gamma_{23}^{N}\sum_{j=0}^{N-1}\binom{N+j-1}{j}\Gamma_{32}^{j}=I_{\Gamma_{32}}(0,N)-I_{\Gamma_{32}}(N,N)\simeq
1.$ (130)

Next, we consider the second term in Eq. (E.1), which can be directly
evaluated as

$\displaystyle-\frac{\Gamma_{23}^{N}}{N}$
$\displaystyle\sum_{j=0}^{N-1}j\binom{N+j-1}{j}\Gamma_{32}^{j}$
$\displaystyle=-\frac{\Gamma_{32}}{\Gamma_{23}}\Gamma_{23}^{N+1}\sum_{j=-1}^{N-2}\binom{(N+1)+j-1}{j}\Gamma_{32}^{j}$
$\displaystyle=-\frac{\Gamma_{32}}{\Gamma_{23}}\quantity[(I_{\Gamma_{32}}(0,N+1)-I_{\Gamma_{32}}(N-1,N+1)]$
$\displaystyle\simeq-\frac{\Gamma_{32}}{\Gamma_{23}},$ (131)

where in the last line we used the asymptotic expansion of $I_{x}(a,b)$ to
approximate the difference between regularised beta functions. Collecting all
the terms, we conclude that the limit is given by

$\lim_{N\to\infty}\sum_{j=1}^{N}\evaluated{\frac{b^{(N)}_{j}}{N}}_{\begin{subarray}{c}b=1\\\
c=0\end{subarray}}=1-\frac{\Gamma_{32}}{\Gamma_{23}}.$ (132)

Therefore, combining the above with Eq. (128), we get the desired limit:

$\lim_{N\to\infty}\frac{1}{N}\sum_{j=1}^{N}\frac{b^{(N)}_{j}}{N}=c+b\,\quantity(1-\frac{\Gamma_{32}}{\Gamma_{23}}).$
(133)

###  E.2 Proof of limit (121b)

As before, in order to prove Eq. (121b), we will also need to prove two other
limits. Recall that $c^{(N)}_{j}$ is given by

$\displaystyle\\!\\!\\!\\!c^{(N)}_{j}=a\Gamma^{j-1}_{13}\sum_{i=1}^{N}$
$\displaystyle\binom{N+j-1-i}{j-1}\Gamma_{31}^{N+1-i}$
$\displaystyle\quad\quad\quad\quad+\Gamma^{N}_{31}\sum_{l=0}^{N-1}\Gamma_{13}^{l}\binom{N+l-1}{l}c_{j-l},$
(134)

with

$c_{j}=c\Gamma^{N}_{32}\sum_{i=0}^{j-1}\binom{N+i-1}{i}\Gamma^{i}_{23}+b\Gamma_{32}\Gamma^{j}_{23}\sum_{i=0}^{N-1}\binom{i+j}{j}\Gamma_{32}^{i}.$
(135)

Since $c_{j-l}$ does not depend on $a$, the problem reduces to showing that

$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{j=1}^{N}\evaluated{c^{(N)}_{j}}_{\begin{subarray}{c}a=1\\\
b,c=0\end{subarray}}$
$\displaystyle=\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{j=1}^{N}\Gamma_{13}^{j}\sum_{i=1}^{N}\binom{N+j-i}{j}\Gamma_{31}^{N+1-i}$
$\displaystyle=\frac{\Gamma_{31}}{\Gamma_{13}},$ (136)

and

$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{j=1}^{N}\evaluated{c^{(N)}_{j}}_{\begin{subarray}{c}a=0\end{subarray}}$
$\displaystyle=\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{j=1}^{N}\Gamma_{31}^{N}\sum_{l=0}^{N-1}\Gamma_{13}^{l}\binom{N+l-1}{l}c_{j-l}$
$\displaystyle=0.$ (137)

Let us start by proving Eq. (136). First, we manipulate $c^{(N)}_{j}$ and
rewrite it in terms of the incomplete beta function as follows:

$\displaystyle\evaluated{c^{(N)}_{j}}_{\begin{subarray}{c}a=1\\\
b,c=0\end{subarray}}$
$\displaystyle=\Gamma_{13}^{j-1}\sum_{i=1}^{N}\binom{N+j-1-i}{j-1}\Gamma_{31}^{N+1-i}$
$\displaystyle=\sum_{i=1}^{N}\quantity[(I_{\Gamma_{13}}(0,N+1-i)-I_{\Gamma_{13}}(N,N+1-i)]$
$\displaystyle=\sum_{i=1}^{N}I_{\Gamma_{31}}(N+1-i,N).$ (138)

Using Eq. (76), we obtain

$\displaystyle\sum_{i=1}^{N}I_{\Gamma_{31}}(N+1-i,N)=N\frac{\Gamma_{31}}{\Gamma_{13}}.$
(139)

Thus, collecting all the terms, we get the desired limit

$\lim_{N\to\infty}\frac{1}{N}\sum_{j=1}^{N}\evaluated{c^{(N)}_{j}}_{\begin{subarray}{c}a=1\\\
b,c=0\end{subarray}}=\frac{\Gamma_{31}}{\Gamma_{13}}.$ (140)

The final step is to show that the remaining limit from Eq. (137) is zero,
namely

$\lim_{N\rightarrow\infty}\frac{1}{N}\Gamma_{31}^{N}\sum_{k=1}^{N}\sum_{l=0}^{N-1}\Gamma_{13}^{l}c_{k-l}=0.$
(141)

Since $c_{k-l}$ has two contributions, one needs to show that both limits go
to zero. Treating each separately, we first write the first term of Eq. (135)
in terms of the incomplete beta function,

$\evaluated{c_{k}}_{\begin{subarray}{c}c=1,b=0\end{subarray}}=\Gamma^{N}_{32}\sum_{i=0}^{k-1}\binom{N+i-1}{i}\Gamma^{i}_{23}=I_{\Gamma_{32}}(N,k)\leq
1,$ (142)

where the upper bound comes from the incomplete beta function being a CDF.
Thus,

$\displaystyle\sum_{k=1}^{N}\frac{\Gamma^{N}_{31}}{N}\sum_{l=0}^{N-1}\Gamma_{13}^{l}$
$\displaystyle\binom{N+l-1}{l}\evaluated{c_{k-l}}_{\begin{subarray}{c}c=1,b=0\end{subarray}}$
$\displaystyle\leq\sum_{k=1}^{N}\frac{\Gamma^{N}_{31}}{N}\sum_{l=0}^{N-1}\Gamma_{13}^{l}\binom{N+l-1}{l}$
$\displaystyle=I_{\Gamma_{31}}(N,N),$ (143)

and this term goes to zero for $\Gamma_{31}\leq 1/2$. This can be seen from
the asymptotic expansion of $I_{\Gamma_{31}}(N,N)$.

Finally, we need to show that the second term of Eq. (135) is zero. To do so,
we first re-write the second term as

$\displaystyle\\!\\!\\!\\!\sum_{k=1}^{N}\evaluated{c^{(N)}_{k}}_{\begin{subarray}{c}b=1\\\
a,c=0\end{subarray}}\\!\\!$
$\displaystyle=\frac{1}{N}\Gamma_{31}^{N}\sum_{l=0}^{N-1}\Gamma_{13}^{l}\binom{N+l-1}{l}\sum_{k=1}^{N}\evaluated{c_{k-l}}_{\begin{subarray}{c}c=0\\\
b=1\end{subarray}}$
$\displaystyle=\frac{1}{N}\Gamma_{31}^{N}\sum_{l=0}^{N-1}\Gamma_{13}^{l}\binom{N+l-1}{l}\sum_{k=l+1}^{N}\evaluated{c_{k-l}}_{\begin{subarray}{c}c=0\\\
b=1\end{subarray}}$
$\displaystyle=\frac{1}{N}\Gamma_{31}^{N}\sum_{l=0}^{N-1}\Gamma_{13}^{l}\binom{N+l-1}{l}\sum_{k=1}^{N-l}\evaluated{c_{k}}_{\begin{subarray}{c}c=0\\\
b=1\end{subarray}}.$ (144)

Notice that the above expression can be bounded by

$\displaystyle\sum_{k=1}^{N}\evaluated{c^{(N)}_{k}}_{\begin{subarray}{c}b=1\\\
a,c=0\end{subarray}}$
$\displaystyle\leq\frac{1}{N}\quantity(\Gamma_{31}^{N}\sum_{l=0}^{N-1}\Gamma_{13}^{l}\binom{N+l-1}{l})\quantity(\sum_{k=1}^{N}\evaluated{c_{k}}_{\begin{subarray}{c}c=0\\\
b=1\end{subarray}}),$

and the right-hand side of equation Eq. (E.2) can be expressed in terms of the
regularised beta function as follows:

$\displaystyle\frac{1}{N}\Gamma_{31}^{N}$
$\displaystyle\sum_{l=0}^{N-1}\Gamma_{13}^{l}\binom{N+l-1}{l}\quantity(\sum_{k=1}^{N}\evaluated{c_{k}}_{\begin{subarray}{c}c=0\\\
b=1\end{subarray}})$
$\displaystyle=\Gamma_{32}\frac{1}{N}I_{\Gamma_{31}}(N,N)\sum_{k=1}^{N}\Gamma^{k}_{23}\sum_{i=0}^{N-1}\binom{i+k}{k}\Gamma_{32}^{i}$
$\displaystyle=\Gamma_{32}\frac{1}{N}I_{\Gamma_{31}}(N,N)\sum_{k=1}^{N}\Gamma^{k}_{23}\Gamma_{23}^{-1-k}\quantity[1-I_{\Gamma_{32}}(N,k+1)]$
$\displaystyle=\frac{\Gamma_{32}}{\Gamma_{23}}\frac{1}{N}I_{\Gamma_{31}}(N,N)\quantity(\frac{1}{N}\sum_{k=1}^{N}I_{\Gamma_{23}}(k+1,N)).$
(145)

Now, note that the first term goes to zero when $\Gamma_{31}<1/2$, whereas the
second term is also bounded by one as

$\frac{1}{N}\sum_{k=1}^{N}I_{\Gamma_{23}}(k+1,N)\leq\frac{1}{N}\sum_{k=1}^{N}1=1.$
(146)

Since

$0\leq\lim_{N\to\infty}\sum_{k=1}^{N}\frac{1}{N}\evaluated{c^{(N)}_{k}}_{\begin{subarray}{c}b=1\\\
a,c=0\end{subarray}}\leq 0,$ (147)

we conclude that the resulting limit is zero. Therefore, combining (143) and
(147),

$\lim_{N\rightarrow\infty}\frac{1}{N}\Gamma_{31}^{N}\sum_{k=1}^{N}\sum_{l=0}^{N-1}\Gamma_{13}^{l}c_{k-l}=0.$
(148)

Collecting together Eqs. (133), (140) and (148) we conclude that Theorem 5 is
proved.∎

##  Appendix F Alternative protocols realising equivalent $\beta$-swap
approximation

Let us revisit the truncated protocol $\tilde{\mathcal{P}}^{(ij)}$ with
$d$-dimensional memory as introduced in Section III, where it was defined in
Eq. (16). This protocol is composed of $d^{2}$ two-level elementary
thermalisations, denoted as $T_{kl}$. Each thermalisation $T_{kl}$ can be
represented as a point $(k,l)$ on a plane, and an algorithm can be represented
as an arrow pointing from the previous thermalisation to the next one. For
instance, we present below the diagram of $\tilde{\mathcal{P}}$ for $d=7$:

Visually, it is obvious that we iterate through an entire column before
shifting to the next one. In other words, all the ’filled’ levels are used
sequentially to fill up the first ’empty’ level, and the same process is
repeated for all the subsequent ’empty’ levels. After investigating the
following two algorithms

The generalisation of the Blue algorithm can be summarised in the following
way, starting with $i=1$:

  1. 1.

Iterate through $i$-th column, starting from the first unvisited point.

  2. 2.

Iterate through $i$-th row, starting from the first unvisited point.

  3. 3.

Set $i\rightarrow i+1$ and go back to step 1 if any unvisited point remains.

The Red algorithm can be most easily understood as the reverse of the Blue
algorithm. Instead of decreasing the length of vertical and horizontal
stretches, they are gradually increased in the Red algorithm. This allows the
Red algorithm to be recursively implemented, taking into account gradually
more and more levels of memory, as indicated by dashed lines. The protocol for
$d$-dimensional memory is implemented by extending the $d-1$-dimensional
version with an additional row and column.

We furthermore investigated the Cyan family, which mimics the Blue algorithm
and can be defined in the following manner:

  1. 1.

Iterate through row or column, starting from the first unvisited point.

  2. 2.

Repeat step 1. until no unvisited points remain.

Moreover, we considered Orange family related to the Cyan family with an
analogous reversal as between Blue and Red.

Figure 14: Convergence to the target state $\bm{q}=(0,1)$ of different
algorithms acting on the state $\bm{p}=(1,0)$ extended by memory with
dimension $d=30$. Note that all algorithms. Note that all algorithms finish at
the same value of $\norm{\bm{p}-\bm{q}}_{1}$.

For these algorithms, we find the following properties we have observed from
explicit implementation for a range of dimensions and inverse temperatures
$\beta$, but we have not been able to prove them analytically:

  * •

All of the aforementioned algorithms acting on an initial state
$\bm{p}\otimes\gamma_{d}$ result in the same state as
$\tilde{\mathcal{P}}(\bm{p}\otimes\bm{\gamma}_{d})$.

  * •

Blue and Red algorithms are slowest and fastest algorithms, respective,
according to the convergence to the $\beta$-swap with respect to the 1-norm.

  * •

Each algorithm in the Cyan and Orange family provide slower and faster
convergence than the original algorithm $\tilde{\mathcal{P}}$, respectively.

The statement reinforces the claims mentioned earlier and explains that the
1-norm between the intermediate states of the system and the target state (the
$\beta$-swapped counterpart of $\bm{p}$ with $\beta=0$) is plotted for
simplicity (see Fig. 14). The memory dimension is specified as $d=30$.

## References

  * Maxwell [1872] J. Maxwell,  _Theory of Heat_, Text-books of science (Longmans, Green, and Company, 1872).
  * Szilard [1929] L. Szilard, Zeitschrift für Physik 53, 840 (1929).
  * Brillouin [1951] L. Brillouin, J. Appl. Phys. 22, 334 (1951).
  * Landauer [1961] R. Landauer, IBM J. Res. Dev. 5, 183 (1961).
  * Bennett [1982] C. H. Bennett, Int. J. Theor. Phys. 21, 905 (1982).
  * Maruyama  _et al._ [2009] K. Maruyama, F. Nori, and V. Vedral, Rev. Mod. Phys. 81, 1 (2009).
  * Seifert [2012] U. Seifert, Rep. Prog. Phys. 75, 126001 (2012).
  * Sagawa [2012] T. Sagawa, Prog. Theor. Exp. Phys. 127, 1 (2012).
  * Goold  _et al._ [2016] J. Goold, M. Huber, A. Riera, L. del Rio, and P. Skrzypczyk, J. Phys. A Math. Theor. 49, 143001 (2016).
  * Binder  _et al._ [2018] F. Binder, L. A. Correa, C. Gogolin, J. Anders, and G. Adesso, Fundam. Theor. Phys. 195, 1 (2018).
  * Deffner and Campbell [2019] S. Deffner and S. Campbell,  _Quantum Thermodynamics. An introduction to the thermodynamics of quantum information_ (Morgan Claypool Publishers, 2019).
  * Taranto  _et al._ [2020] P. Taranto, F. Bakhshinezhad, P. Schüttelkopf, F. Clivaz, and M. Huber, Phys. Rev. Appl. 14, 054005 (2020).