observations following from embedding and projection operations. First, note
that majorisation in embedded space implies majorisation between projections,
thus
$\boldsymbol{t}^{\mathfrak{M}}\succ_{\mathfrak{M}}\boldsymbol{q}\Rightarrow\mathfrak{P}_{\boldsymbol{\pi}}(\boldsymbol{t}^{\mathfrak{M}})\succ_{\beta}\mathfrak{P}_{\boldsymbol{\pi}}(\boldsymbol{q})$
for every order $\boldsymbol{\pi}$. Second, observe that embedding preserves
the majorisation relations between the vectors, therefore
$\boldsymbol{q}\in\mathcal{T}_{\emptyset}(\boldsymbol{p})\Leftrightarrow\mathfrak{M}(\boldsymbol{q})\in\mathcal{T}_{\emptyset}(\mathfrak{M}(\boldsymbol{p}))$.
These two statements show that it is enough to consider vectors
$\boldsymbol{t}^{(n,\boldsymbol{\pi})}=\mathfrak{P}_{\boldsymbol{\pi}}(\boldsymbol{t}^{\mathfrak{M}(n)})$
and convex combinations of their future thermal cones, thus proving Lemma 2.5.
∎

####  Subsection 3 Construction of probabilistic majorisation cones

In this section, we derive the results presented in Section 2. Specifically,
we explore how to extend the notion of majorisation cones to probabilistic
ones. For the convenience of the reader, we restate the theorem concerning
probabilistic transformations:

See 1.4

To establish a direct connection to majorization, we reformulate the above
theorem as follows:

$\forall_{1\leq k\leq
d}:\mathcal{P}(\boldsymbol{p},\boldsymbol{q})\leq\frac{\sum_{j=k}^{d}p^{\downarrow}_{j}}{\sum_{j=k}^{d}q^{\downarrow}_{j}}=\frac{1-\sum_{j=1}^{k-1}p^{\downarrow}_{j}}{1-\sum_{j=1}^{k-1}p^{\downarrow}_{j}}=\frac{1-P_{k}}{1-Q_{k}}.$
(52)

By setting $\mathcal{P}(\boldsymbol{p},\boldsymbol{q})=1$, we recover the
standard majorisation condition on deterministic convertibility,

$\forall_{1\leq k\leq
d}:1\leq\frac{\sum_{j=k}^{d}p^{\downarrow}_{j}}{\sum_{j=k}^{d}q^{\downarrow}_{j}}\Leftrightarrow\boldsymbol{p}\prec\boldsymbol{q}.$
(53)

To determine the probabilistic past cone
$\mathcal{T}_{-}(\boldsymbol{p},\mathcal{P})$ at probability $\mathcal{P}$, we
consider

$\displaystyle\forall_{1\leq k\leq d}:\mathcal{P}\leq\frac{1-Q_{k}}{1-P_{k}}$
$\displaystyle\Rightarrow\mathcal{P}-\mathcal{P}P_{k}\leq 1-Q_{k}\Rightarrow
Q_{k}\leq\mathcal{P}P_{k}+(1-\mathcal{P})$ (54)
$\displaystyle\Rightarrow\boldsymbol{q}\prec\tilde{\boldsymbol{p}},$

with an auxiliary distribution

$\tilde{p}_{i}=\begin{cases}\mathcal{P}p^{\downarrow}_{1}+(1-\mathcal{P})&\text{for
}i=1,\\\ \mathcal{P}p^{\downarrow}_{i}&\text{otherwise},\end{cases}$ (55)

which is always a proper probability distribution ordered non-increasingly,
$\tilde{\boldsymbol{p}}=\tilde{\boldsymbol{p}}^{\downarrow}$, therefore
providing a proper Lorenz curve.

Following a similar procedure for the future cone
$\mathcal{T}_{+}(\boldsymbol{p},\mathcal{P})$ leads to

$\displaystyle\mathcal{P}\leq\frac{1-P_{k}}{1-Q_{k}}$
$\displaystyle\Rightarrow\mathcal{P}^{-1}-\mathcal{P}^{-1}P_{k}\geq
1-Q_{k}\Rightarrow Q_{k}\geq\mathcal{P}^{-1}P_{k}+(1-\mathcal{P}^{-1})$ (56)
$\displaystyle\Rightarrow\boldsymbol{q}\succ\hat{\boldsymbol{p}},$

with the second auxiliary distribution

$\hat{p}_{i}=\begin{cases}\mathcal{P}^{-1}p^{\downarrow}_{1}+(1-\mathcal{P}^{-1})&\text{for
}i=1,\\\ \mathcal{P}^{-1}p^{\downarrow}_{i}&\text{otherwise}.\end{cases}$ (57)

In contrast to the case of the past cone, the distribution
$\hat{\boldsymbol{p}}$ in this formulation is not ordered in a non-increasing
manner beyond a certain value of $\mathcal{P}$. At first glance, one might
think that reordering should solve the problem; however, it would be
equivalent to a decrease in the probabilistic future with decreasing
$\mathcal{P}$, which creates a contradiction. The solution is provided by
noting that Vidal’s criterion deals with rescaled entries of the Lorenz curve
rather than the probabilities themselves. Therefore, the Lorenz curve for
$\hat{\boldsymbol{p}}$ should remain convex for all values of $\mathcal{P}$
without the need for reordering.

Consider the following critical values of $\mathcal{P}$, namely,

$\mathcal{P}_{n}=(n-1)p^{\downarrow}_{n}-\sum_{i=1}^{n-1}p^{\downarrow}_{i}+1,$
(58)

for which the first $n$ entries of the distribution $\hat{\boldsymbol{p}}$
will not be ordered non-increasingly, resulting in an improper Lorenz curve.
The resulting non-convexity is controlled by replacing

$\left\\{\hat{p}_{1},\ldots,\hat{p}_{n}\right\\}\rightarrow\frac{1}{n}\sum_{i=1}^{n}\hat{p}_{i}\left\\{1,\ldots
1\right\\},$ (59)

which ensures that $\hat{\boldsymbol{p}}=\hat{\boldsymbol{p}}^{\downarrow}$.

This way, the auxiliary ordered distributions $\tilde{\boldsymbol{p}}$ and
$\hat{\boldsymbol{p}}$, together with the construction for the deterministic
majorisation cones provide the full construction of the probabilistic cones.
However, it is important to note that when considering entanglement and
coherence theories, the roles of the future and past majorisation cones are
reversed, and this should be taken into account.

####  Subsection 4 Volumes of entanglement majorisation cones

Finally, here we present the methods used to obtain the isovolumetric sets for
entanglement cones, as discussed in Section 3. Additionally, we plot these
curves for systems of unequal dimensions to compare with the case of equal
dimension.

Consider a uniform Haar distribution of pure states in a composed space
$\ket{\psi}\in\mathcal{H}^{N}\otimes\mathcal{H}^{M}$ with $N\leq M$. The
partial tracing induces a measure in the space of reduced states,
$\rho=\operatorname{Tr}_{2}\ket{\psi}\bra{\psi}$, characterised by the
distribution of eigenvalues
$\Lambda=\left\\{\lambda_{1},\ldots,\lambda_{N}\right\\}$ of the reduced state
[178]:

$P_{N,M}(\Lambda)=C_{N,M}\delta\left(1-\sum_{i}\lambda_{i}\right)\prod_{i}\lambda_{i}^{M-N}\theta(\lambda_{i})\prod_{i<j}(\lambda_{i}-\lambda_{j})^{2},$
(60)

where $\delta$ and $\theta$ are the Dirac delta and Heavyside step functions,
respectively. The normalisation constant is given by

$C_{N,M}=\frac{\Gamma(NM)}{\prod_{j=0}^{N-1}\Gamma(M-j)\Gamma(N-j+1)},$ (61)

where $\Gamma$ is the Gamma function. Before we continue with our discussion,
let us first understand the role played by all factors in Eq. (60). As
previously mentioned, $C_{N,M}$ is the normalisation. The delta function
ensures that the spectrum sums to one (is normalised), whereas the step
function will guarantee that it is positive. Now, notice that the first
product (capital pi notation) essentially does not influence the behaviour of
the distribution for $N=M$ and introduces the repelling of the faces of the
probability simplex otherwise, as it goes to zero whenever $\lambda_{i}=0$ for
any $i$. The second product is responsible for the repelling from
distributions with any two entries equal, since it goes to zero whenever
$\lambda_{i}=\lambda_{j}$ for any $i\neq j$.

Figure 10:  _Isovolumetric sets for entanglement $3\times M$ bipartite
systems_. The density $P_{3,M}(\Lambda)$ of Schmidt coefficients of pure
states for qutrit-quMit systems depends heavily on the dimension $M$ of the
second system. Panels (a-d) and (e-h) present the isovolumetric lines for
past, incomparable and future regions and the density of the states for $M=6$
and $M=30$, respectively. Note that for larger $M$ the density $P_{3,M}$ is
more and more concentrated around the regions close to the centre [compare (d
and (h)]. This affects the subset of states with large future volume, making
it smaller [(c) and (g)] as well as the set of states with large past volume,
enlarging it [(a) with (e)].

Sampling from the $P_{N,M}$ distribution, ordinarily done by generating state
vectors $\ket{\psi}\in\mathcal{H}_{NM}$ which would be computationally
prohibitive for large $M$, can be achieved using only $O(N)$ random numbers
for any dimension of the secondary system. It has been demonstrated in Ref.
[172] that the distribution $P_{N,M}$ is precisely the Laguerre unitary
ensemble generated by Wishart matrices of size $N$ and parameter $M$ [181]
and, in turn, generated using a tridiagonal method containing only $O(N)$
random real numbers [182], which indeed allows one to study the $P_{N,M}$
distributions for arbitrary high-dimensional ancillary systems.

The procedure for generating the isovolumetric lines for the majorisation
cones with a given distribution, $P_{N,M}$, proceeds as follows:

  1. 1.

Generate a sample of $n$ sets of eigenvalues
$\\{\Lambda_{1},\ldots,\Lambda_{n}\\}$ taken from the distribution $P_{N,M}$
using the tridiagonal method.

  2. 2.

Consider regularly spaced grid of points $S$ in a single chamber of the full
probability simplex $\Delta_{N}$ (e.g., $p_{1}\geq p_{2},...,p_{d})$ in order
to avoid repeated counting (achieving $N!$ decrease in operations)

  3. 3.

For each $\boldsymbol{p}\in S$ consider its majorisation cones
$\mathcal{T}_{i}(\boldsymbol{p})$ and divide $S$ into
$S_{i}\equiv\left\\{\Lambda_{i}\in\mathcal{T}_{i}(\boldsymbol{p})\right\\}$.

  4. 4.

This way we arrive at the approximations of the volumes of the three regions,

$\mathcal{V}_{i}\approx\frac{|S_{i}|}{n},$

where $|X|$ denotes the number of elements in a set $X$.

We applied this method for $N=3$ with $M=3$, displayed in the main text of
Fig. 7, and additionally with $M=6$ and $30$, as shown in Fig. 10. These two
cases show the significant dependence of the isovolumetric lines on the size
of the environment.

###  Section 6 Concluding remarks

In this chapter, we investigated the structure of the thermodynamic arrow of
time by analysing thermal cones and their behavior under different conditions.
By dividing the space of energy-incoherent states into the past, the future,
and the incomparable region, in analogy with the future, past, and spacelike
regions of Minkowski spacetime, we identified thermal cones as regions of the
probability simplex that encode the achievability of state transformations
under thermal operations. We specifically focused on energy-incoherent states
in the presence of a thermal bath at a finite temperature and in the limit of
temperature going to infinity, fully characterising and carefully analysing
the incomparable and past thermal cones in both regimes. Additionally, we
identified the volumes of the thermal cones as thermodynamic monotones and
performed a detailed analysis of their behavior. Our results can be applied
directly to the study of entanglement, as the order defined on the set of
bipartite pure entangled states by local operations and classical
communication is the opposite of the thermodynamic order in the limit of
infinite temperature. In this context, the future thermal cone becomes the
past for entanglement, and the past becomes the future. Furthermore, a similar
extension can be drawn to coherence resource theory.

There are several potential research directions for generalising and extending
the results presented in this chapter. One possible avenue is to expand the
analysis beyond energy-incoherent states to encompass the full space of
quantum states. While available tools for this purpose are comparatively
scarce [153, 179, 183], various techniques can be employed to construct
coherent thermal cones under both Gibbs-preserving and thermal operations for
qubit systems, as shown in Fig. 8 and explained in Section 4. Additionally, it
would be interesting to explore an equivalent construction of the past and
incomparable regions for continuous thermomajorisation. Extending our analysis
to this setting could provide further insight into the structure of the
thermodynamic arrow of time for memoreless processes. Futhermore, an extension
that includes many non-interacting subsystems (possibly independent and
identically distributed) could be done by defining an appropriate function to
analyse the behavior of the thermal cones. While our investigation focused on
single subsystems, such an extension could be useful in exploring the behavior
of thermodynamic quantities in a broader range of physical systems.

Recently, the framework outlined in this section was applied in Ref. [184] for
exploring the thermodynamic constraints on the pivotal task of generating
entanglement using non-equilibrium states. In that work, the authors present a
detailed construction of the future thermal cone of entanglement – the set of
entangled states that a given separable state can thermodynamically evolve to.
In a similar spirit to the discussion presented here, the properties of the
future thermal cone of entanglement were studied and related to the ability to
generate entanglement.

Finally, the fundamental limits and advantages of using a catalyst to aid
thermodynamic transformations within the framework of thermal cones were
addressed in Ref. [185]. More precisely, the set of states to which a given
initial state can thermodynamically evolve (the catalyzable future) or from
which it can evolve (the catalyzable past) with the help of a strict catalyst
was characterized. Secondly, lower bounds on the dimensionality required for
the existence of catalysts under thermal processes, along with bounds on the
catalyst’s state preparation, were found.

##  Chapter 5 Memory-assisted Markovian thermal processes

Information has become ubiquitous in thermodynamics. It all started with
Maxwell’s seminal inquiry [186]: _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 [187,
188]. Ultimately, it was found that thermodynamics imposes physical
restrictions on information processing [189, 190], resulting in the
development of frameworks devoted to incorporating information into
thermodynamics [191, 192, 193, 35, 36, 37]. 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 [194],
generation of entanglement [195, 196] or improved performance of heat engines
and refrigerators [197, 198, 199]. 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 [53, 54]. 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
[200].

As we have seen in Chapter 3, the resource theory of thermodynamics is a
relatively recent framework allowing one to address foundational questions in
thermodynamics. By relying on the notion of thermal operations 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 [60, 145], 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. [154] 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 [201]. 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.

In this chapter, 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 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).

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.

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 [202, 203], superconducting circuits [204, 205], and atom-cavity
systems [206]. 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. [154, 155], 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 chapter is structured as follows. First, in Sec. 1 we compare the
frameworks of thermal operations and Markovian thermal processes, and then
introduce the central notion of this chapter, the memory-assisted Markovian
thermal processes. Next, in Sec. 2, we 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. Section 3 contains
discussion and application of our results. Finally, in Sec. 6, we conclude and
provide outlook for future research.

###  Section 1 Thermal operations vs Markovian thermal processes

_Thermal operations_ (TOs) framework as discussed in Chapter 3 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. 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
[154].

_Markovian thermal processes_ (MTPs) framework [60, 145], on the other hand,
uses typical assumptions of the theory of open quantum systems (weak coupling,
large bath size, quickly decaying correlations, etc.) [53], to argue that the
system undergoes an open dynamics described by a Lindblad master equation [69,
142, 70]. Since the dynamics generated by an MTP arises explicitly from a
Markovian model, there are no memory effects. Also, as shown in Ref. [60], 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},$
(1a) $\displaystyle p_{j}$
$\displaystyle\rightarrow(1-\lambda)p_{j}+\lambda\frac{p_{i}+p_{j}}{\gamma_{i}+\gamma_{j}}\gamma_{j},$
(1b)

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
$\boldsymbol{E}=(0,1,2)$ and prepared in an energy-incoherent state
$\boldsymbol{p}=(0.7,0.2,0.1)$ (depicted by a black dot $\bullet$) can be
transformed to [green region $\mathcal{T}_{+}(\boldsymbol{p})$] or transformed
from [blue region $C_{-}(\boldsymbol{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 $\mathcal{T}^{\textrm{TO}}_{+}(\boldsymbol{p})$ achievable
via thermal operations from a given incoherent initial state $\boldsymbol{p}$
can be fully characterised using the notion of thermomajorisation. The so-
called future thermal cone $\mathcal{T}^{\textrm{TO}}_{+}(\boldsymbol{p})$
(see Chapter 4) is a convex set that consists of at most $d!$ extreme points,
the construction of which was given in Corollary 1.1 (see Fig. 2a for an
example with a three-level system). On the other hand, the set of states
$\mathcal{T}^{\textrm{MTP}}_{+}(\boldsymbol{p})$ achievable via Markovian
thermal processes from a state $\boldsymbol{p}$ was recently characterised
using the notion of continuous thermomajorisation. The future Markovian
thermal cone $\mathcal{T}^{\textrm{MTP}}_{+}(\boldsymbol{p})$ is not convex
(as illustrated in Fig. 2b for a three-level system), but Theorem 4.2 provides
a construction of its extreme points using sequences of two-level full
thermalisations (i.e., transformations from Eqs. (1a)-(1b) with $\lambda=1$).
As can be seen in Fig. 2c,
$\mathcal{T}^{\textrm{MTP}}_{+}(\boldsymbol{p})\subset\mathcal{T}^{\textrm{TO}}_{+}(\boldsymbol{p})$
and the difference between these two sets of thermodynamically accessible
states arises purely from memory effects.

###  Section 2 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
$\boldsymbol{p}$ via thermal operations [i.e., any $\boldsymbol{q}\in
C_{+}^{\mathrm{TO}}(\boldsymbol{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}}(\boldsymbol{p})$. Next, we will introduce MeMTP protocols
that will serve us to approach extreme points of
$C^{\textrm{TO}}_{+}(\boldsymbol{p})$ using MTPs acting on the system and
memory state, $\boldsymbol{p}\otimes\boldsymbol{\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.

####  Subsection 1 Memory-assisted Markovian thermal processes

In this chapter, our objective is to interpolate between the two extreme
regimes characterised by arbitrarily strong and no memory effects, as
described within the TO and MTP frameworks. To accomplish this, we will
concentrate on the more restrictive MTP framework and extending it by
explicitly modeling 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. (3) 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 2.1 (Memory-assisted MTPs).

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

$\mathcal{E}(\rho)=\textrm{Tr}_{M}[\mathcal{E}_{\mathrm{MTP}}(\rho\otimes\gamma_{M})],$
(2)

where $\mathcal{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}$.

As already mentioned, in this chapter we will focus on transformations between
energy-incoherent states. Our aim is to show that the sets of states
achievable from a given $\boldsymbol{p}$ via memory-assisted MTPs interpolate
between $\mathcal{T}^{\textrm{MTP}}_{+}(\boldsymbol{p})$ (for $N=1$) and
$\mathcal{T}^{\textrm{TO}}_{+}(\boldsymbol{p})$ (for $N\to\infty$). As a final
note, observe that this framework can be formally related to a particular kind
of catalytic transformations [89, 57]. 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.

[frametitle= 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.

$\begin{matrix}\text{Memoryless}&\boldsymbol{p}\xrightarrow[]{\text{MTP}}\boldsymbol{q}&\Longleftrightarrow&\boldsymbol{p}{}_{\beta}\boldsymbol{q}\\\
\biggl{|}&&\vdots\\\ \text{Finite
memory}&\boldsymbol{p}\xrightarrow[]{\text{MeMTP}}\boldsymbol{q}&\Longleftrightarrow&\boldsymbol{p}\otimes\boldsymbol{\gamma}_{M}{}_{\beta}\boldsymbol{q}\otimes\boldsymbol{\gamma}_{M}\\\
\bigg{\downarrow}&&\>\>\begin{matrix}\vdots\\\ {}_{N\to\infty}\\\
\downarrow\end{matrix}\\\ \text{Full
control}&\boldsymbol{p}\xrightarrow[]{\text{TO}}\boldsymbol{q}&\Longleftrightarrow&\boldsymbol{p}\succ_{\beta}\boldsymbol{q}\end{matrix}$

####  Subsection 2 Simplification to extreme points

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

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

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

$\boldsymbol{q}\in\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{p})\Rightarrow\boldsymbol{q}\in\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{p}^{\pi_{\boldsymbol{q}}}),$
(3)

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

$\left[\boldsymbol{q}\in\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{p})\leavevmode\nobreak\
\leavevmode\nobreak\ \mathrm{and}\leavevmode\nobreak\ \leavevmode\nobreak\
\pi_{\boldsymbol{q}}=\pi_{\boldsymbol{p}}\right]\Rightarrow\boldsymbol{q}\in\mathcal{T}_{+}^{\mathrm{MTP}}(\boldsymbol{p}),$
(4)

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
$\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{p})$, then one can also get to
every state in $\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{p})$ via MeMTPs.
This is done by simply first transforming $\boldsymbol{p}$ to a given extreme
point $\boldsymbol{p}^{\pi}$ of
$\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{p})$, and then using MTPs to get
from $\boldsymbol{p}^{\pi}$ to every state with a $\beta$-order $\pi$ in
$\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{p})$.

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

$\delta(\boldsymbol{p},\boldsymbol{q}):=\frac{1}{2}\sum_{i=1}^{d}|p_{i}-q_{i}|.$
(5)

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\boldsymbol{q}\in\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{p}):\quad\min_{{\cal
P}\in\mathrm{MeMTP}}\delta({\cal
P}(\boldsymbol{p}),\boldsymbol{q})\leq\epsilon.$ (6)

A particular subset of extreme points of
$\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{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$ [154]:

$\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)},$
(7)

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
$\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{p})$ can be obtained by sequences
of transpositions (that is because an extreme point in that case is of the
form $\Pi\boldsymbol{p}$ for a permutation matrix $\Pi$, and every $\Pi$ can
be constructed from transpositions).

For finite temperatures, a $\beta$-swap transforms $\boldsymbol{p}$ into an
extreme point of $\boldsymbol{p}^{\pi}$ if the $\beta$-orders
$\pi_{\boldsymbol{p}}$ and $\pi$ differ only by a transposition of adjacent
elements [154, 62]. In other words, it happens for $\Pi_{ij}^{\beta}$ when
$\pi_{\boldsymbol{p}}(i)=\pi_{\boldsymbol{p}}(j)\pm 1$. More generally, a
sequence of such non-overlapping $\beta$-swaps will also produce an extreme
point of $\mathcal{T}^{\mathrm{TO}}_{+}(\boldsymbol{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 [207]. Finally, we will also make use of
the notion of $\beta$-cycles that we now define. For a state $\boldsymbol{p}$,
consider a $k$-dimensional subset of energy levels $\\{i_{1},\ldots,i_{k}\\}$
neighbouring in the $\beta$-order, i.e.,
$\pi_{\boldsymbol{p}}(i_{j+1})=\pi_{\boldsymbol{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
$\boldsymbol{p}$ to its extreme point $\boldsymbol{p}^{\pi^{\prime}}$ is
called a $\beta$-cycle, if $\Pi^{\prime}=\Pi\Pi_{\boldsymbol{p}}$ (here
$\Pi,\Pi^{\prime}$ and $\Pi_{\boldsymbol{p}}$ denote matrix representations of
permutations $\pi,\pi^{\prime}$ and $\pi_{\boldsymbol{p}}$). To emphasise that
a given $\beta$-cycle acts on $k$ levels, we will sometimes refer to it as a
$\beta$-$k$-cycle.

####  Subsection 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 2.2 (Two-level thermalisations).

Consider a system in a state $\boldsymbol{p}$ with the corresponding thermal
state $\boldsymbol{\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\\}$
(8)

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_{\boldsymbol{p}}(i)=\pi_{\boldsymbol{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
$\mathcal{T}_{+}^{\mathrm{MTP}}$ [60]. Intuitively, one can expect that in
order to approximate extreme states of
$\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{p})$ using MeMTPs, one should get
to the extreme points of
$\mathcal{T}_{+}^{\mathrm{MTP}}(\boldsymbol{p}\otimes\boldsymbol{\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. 2.2 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. 2.2, it is possible to exceed this intuitively unbeatable limit.
[0.34cm]  _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. 2.2 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. 3), which
includes the main system and a memory starting at thermal equilibrium, i.e., a
state $\boldsymbol{p}\otimes\boldsymbol{\gamma}_{M}$. In particular, we focus
on the populations
$[\boldsymbol{p}\,\otimes\,\boldsymbol{\gamma}_{M}]_{N(i-1)+1},\ldots,[\boldsymbol{p}\,\otimes\,\boldsymbol{\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. 3). In the $k$-th round, we select the entry
$[\boldsymbol{p}\otimes\boldsymbol{\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}(\boldsymbol{p}\otimes\boldsymbol{\gamma}_{M}):=\quantity(\prod_{l=1}^{N}T_{(i-1)N+k,\,(j-1)N+l})(\boldsymbol{p}\otimes\boldsymbol{\gamma}_{M}).$
(9)

Note that if $\pi_{\boldsymbol{p}}(i)=\pi_{\boldsymbol{p}}(j)\pm 1$ (i.e., the
$\beta$-orders of $\boldsymbol{p}$ and $\Pi^{\beta}_{ij}\boldsymbol{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)}(\boldsymbol{p}\otimes\boldsymbol{\gamma}_{M}):=\mathcal{R}^{(ij)}_{N}\circ\ldots\circ\mathcal{R}^{(ij)}_{1}(\boldsymbol{p}\otimes\boldsymbol{\gamma}_{M}).$
(10)

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

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

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

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

Figure 3: _$\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.

We will also employ more general protocols that aim at approximating the
transformation of the initial state $\boldsymbol{p}$ into an extreme state of
$\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{p})$ given by
$\boldsymbol{p}^{\pi^{\prime}}$. Denote matrix representations of
$\beta$-orders of these states by $\Pi_{\boldsymbol{p}}$ and
$\Pi^{\prime}=\Pi\Pi_{\boldsymbol{p}}$ for some permutation matrix $\Pi$.
Moreover, let us decompose $\Pi$ into neighbour transpositions with respect to
$\boldsymbol{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}}\boldsymbol{p}$ only by a
transposition of adjacent elements. Then, we define the following two
protocols to approximate $\boldsymbol{p}^{\pi^{\prime}}$:

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

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

####  Subsection 4 Achieving extreme points of $\mathcal{T}_{+}^{\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$-level memory in energy-
incoherent states represented by probability distributions $\boldsymbol{p}$
and $\boldsymbol{\gamma}_{M}$. Since $\beta=0$, the thermal state of the
memory is described by a uniform distribution $\boldsymbol{\eta}_{M}$ with
every entry equal to $1/N$. We start with the following lemma, the proof of
which can be found in Section 2 (with the necessary background on the
mathematical tools used presented in Section 1).

######  Lemma 2.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)}(\boldsymbol{p}\otimes\boldsymbol{\eta}_{M})=\boldsymbol{q}\otimes\boldsymbol{\eta}_{M},$
(14)

with

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

and $\epsilon$ given by

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

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 2.1, we can demonstrate that an
arbitrary permutation can be achieved using a composition of our approximate
protocols.

######  Theorem 2.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}}(\boldsymbol{p}\otimes\boldsymbol{\eta}_{M})=\boldsymbol{q}\otimes\boldsymbol{\eta}_{M},$
(17)

where

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

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

$\boldsymbol{\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}}).$
(19)

###### Proof.

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

$\boldsymbol{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}})]\boldsymbol{p}.$
(20)

Clearly, the leading term is given by $\Pi\boldsymbol{p}$, whereas the next
leading term, proportional to $\epsilon$, is given by
$\epsilon\Delta\boldsymbol{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 $\boldsymbol{p}$ to any state
$\boldsymbol{q}\in\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{p})$ using MeMTPs
with $N$-dimensional memory. We explain how to derive such a bound for a given
$\boldsymbol{p}$ in Section 2, whereas below we present a weaker, but much
simpler, bound that is independent of $\boldsymbol{p}$.

######  Corollary 2.3.

Consider states $\boldsymbol{p}$ and $\boldsymbol{q}\in
C_{+}^{TO}(\boldsymbol{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}(\boldsymbol{p}\otimes\boldsymbol{\eta}_{M})=\boldsymbol{q}^{\prime}\otimes\boldsymbol{\eta}_{M},$
(21)

with

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

###### Proof.

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

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

with

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

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. (4), we know that there exists an MTP protocol
${\cal P}^{\prime}$ mapping $\Pi\boldsymbol{p}$ to $\boldsymbol{q}$. Using the
contractiveness of the total variation distance under stochastic processing,
we then have

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

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

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

######  Conjecture 2.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\boldsymbol{p},\tilde{\boldsymbol{q}}\right)\leq\delta\left(\Pi\boldsymbol{p},\boldsymbol{q}\right),$
(26)

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

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

The conjecture is solidified by strong numerical evidence (see Fig. 4 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. 4, 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\boldsymbol{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 4: _Convergence rates at infinite temperature_. Log-log plot of the
total variation distance $\delta$ between the extreme points
$\boldsymbol{p}^{\pi}\in\mathcal{T}_{+}^{\mathrm{(TO)}}(\boldsymbol{p})$ and
the states obtained from $\boldsymbol{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,
$\boldsymbol{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
$\boldsymbol{p}^{\pi}$ with different convergence rates (from bottom to top
$\boldsymbol{p}^{\pi}$ is obtained from $\boldsymbol{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})$, 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}$.

####  Subsection 5 Achieving extreme points of $\mathcal{T}_{+}^{\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 2.1, the proof of which can be found in Section 2.

######  Theorem 2.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
$\boldsymbol{\eta}_{M}$), $\Pi_{ij}^{\beta}$ can be approximated by the MeMTP
protocol ${\cal P}^{(ij)}$ as follows:

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

with

$\delta(\boldsymbol{q},\Pi_{ij}^{\beta}\boldsymbol{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})]$ (29)

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
$\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{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
$\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{p})$.

######  Conjecture 2.2 (Extreme points of $\mathcal{T}_{+}^{\mathrm{TO}}$).

Consider a state $\boldsymbol{p}$ and the extreme point of
$\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{p})$ given by
$\boldsymbol{p}^{\pi^{\prime}}$, with matrix representations of $\beta$-orders
of these states satisfying $\Pi^{\prime}=\Pi\Pi_{\boldsymbol{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 $\boldsymbol{\eta}_{M}$), the MeMTP protocol $\widetilde{{\cal
P}}^{\Pi}$ acts as

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

with

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

The conjecture is solidified by the following two evidences. First, we can
actually prove it in the following special case.

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

Consider a state $\boldsymbol{p}$ with entries $i_{1},i_{2},i_{3}$ being
neighbours in the $\beta$-order (i.e.,
$\pi_{\boldsymbol{p}}(i_{1})=\pi_{\boldsymbol{p}}(i_{2})+1=\pi_{\boldsymbol{p}}(i_{3})+2$),
and the extreme point of $\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{p})$ given
by $\boldsymbol{p}^{\pi^{\prime}}$, with matrix representations of
$\beta$-orders of these states satisfying
$\Pi^{\prime}=\Pi\Pi_{\boldsymbol{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 $\boldsymbol{\eta}_{M}$), the MeMTP
protocol $\widetilde{{\cal P}}^{\Pi}$ acts as

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

with

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

The proof of the above theorem can be found in Section 2, 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.

The second evidence is based on extensive numerical simulations. Let us
consider any state $\boldsymbol{p}$ and a set of permutations defined by a
recurrence formula

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

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 $\boldsymbol{p}$. For each such permutation, we have
considered the action of the protocol
$\widetilde{\mathcal{P}}^{\widetilde{\Pi}_{ij}}\quantity(\boldsymbol{p})$ and
its convergence to the respective extreme point $\boldsymbol{p}^{\pi}$ with
$\Pi=\widetilde{\Pi}_{ij}\Pi_{\boldsymbol{p}}$ (recall that $\Pi$ is a matrix
representation of $\pi$). In each case, we have observed the convergence of
the form from Eq. (31) that is better than $N^{-1/2}$. Results for an
exemplary state in dimension $d=6$ are presented in Fig. 5, 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 2.4, the proof of Theorem 2.5, and numerical
evidence, one can reasonably strengthen Conjecture 2.2 to make the following
statement on the convergence:

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

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

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

###  Section 3 Discussion and applications

In Section 2, 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.

####  Subsection 1 Work extraction

Figure 6: _$\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.

We start our discussion with $\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}$ [208, 56]. 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
[56] and through continuous thermomajorisation relations when transformations
are given by Markovian thermal processes [145]. 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 $\boldsymbol{p}_{SB}=\boldsymbol{p}\otimes(1,0)$. One
can then select the extreme point
$\boldsymbol{p}^{\pi^{\prime}}_{SB}\in{\mathcal{T}}^{\mathrm{TO}}_{+}\quantity(\boldsymbol{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}}$:

$\boldsymbol{\gamma}\otimes(\epsilon_{\mathrm{TO}},1-\epsilon_{\mathrm{TO}})\in{\mathcal{T}}^{\mathrm{TO}}_{+}\quantity(\boldsymbol{p}^{\pi^{\prime}}_{SB}).$
(36)

In other words, $\boldsymbol{p}^{\pi^{\prime}}_{SB}$ is an intermediate state
from which one can achieve minimal error for extracting $W$ work from
$\boldsymbol{p}$ via any thermal operation. We can now define $\Pi$ as a
permutation that maps the matrix representation of the initial $\beta$-order
of $\boldsymbol{p}_{SB}$ to the final $\beta$-order of
$\boldsymbol{p}_{SB}^{\pi^{\prime}}$. Then, by using the algorithm
$\widetilde{\mathcal{P}}^{\Pi}$, we can transform $\boldsymbol{p}_{SB}$ into a
state $\boldsymbol{q}_{SB}$ that approximates
$\boldsymbol{p}^{\pi^{\prime}}_{SB}$. Finally, due to Eq. (4), we can use
standard thermomajorisation to find the minimal value of $\epsilon_{N}$ for
which the state $\boldsymbol{q}_{SB}$ can be transformed to
$\boldsymbol{\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 $\boldsymbol{p}$ using a memory of size $N$. Numerical
simulations of this procedure (see Fig. 6) 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. 6 clearly
shows that even a small size memory can significantly improve the quality of
the extracted work.

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

Figure 7: _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.

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. 7. 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.

Let us now assume, for simplicity, that the system starts in an excited state
extended by Gibbs memory, $\boldsymbol{p}\otimes\boldsymbol{\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

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

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}(\boldsymbol{p}\otimes\boldsymbol{\gamma}_{M})=\boldsymbol{q}\otimes\boldsymbol{\gamma}_{M}$
with

$\boldsymbol{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)}].$ (38)

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

$\norm{\boldsymbol{q}-\boldsymbol{\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]},$ (39)

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.

####  Subsection 3 Two-level control is sufficient for thermal operations

In Ref. [154] 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. [155], for any dimension $d$, an
explicit final state
$\boldsymbol{q}\in\mathcal{T}_{+}^{\mathrm{TO}}(\boldsymbol{p})$ was given
such that it cannot be achieved (even approximately) starting from the ground
state $\boldsymbol{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\boldsymbol{q}$ $\displaystyle=\left(1-\sum_{i=2}^{d}e^{-\beta
E_{i}},e^{-\beta E_{2}},\dots,e^{-\beta E_{d}}\right)$ (40)

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. [155] that there exists $\epsilon>0$ such that any
$\boldsymbol{q}^{\prime}$ achievable from $\boldsymbol{p}$ satisfies
$\delta(\boldsymbol{q},\boldsymbol{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. 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

$\mathcal{E}(\rho)=\rho\otimes\gamma_{M}$ (41)

is a thermal operation for every $N$. Thus, $\mathcal{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 8: _Convergence to states inaccessible via elementary thermal
operations._ Log-log plot of the total variation distance $\delta$ between the
state $\boldsymbol{q}$ from Eq. (40) and the state obtained from
$\boldsymbol{p}=(1,0,\dots,0)$ via the algorithm
$\widetilde{\mathcal{P}}^{\Pi_{\boldsymbol{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.

We illustrate the above with the following numerical examples, showing that
our MeMTP protocol $\widetilde{\mathcal{P}}^{\Pi_{\boldsymbol{q}}}$ (which
consists of only two-level operations) is able to transform $\boldsymbol{p}$
into $\boldsymbol{q}^{\prime}$ that approximates $\boldsymbol{q}$ arbitrarily
well (i.e., $\delta(\boldsymbol{q}^{\prime},\boldsymbol{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.2, which can be seen in the left panel of Fig. 8 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. 8, it turns out again that the convergence,
independently from the energy structure, is better than $1/\sqrt{N}$, in
accordance with Conjecture 2.2.

Figure 9: _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
$\boldsymbol{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$.

####  Subsection 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 $\boldsymbol{p}_{SM}^{(n)}$. Similarly, let $\boldsymbol{p}^{(n)}$
and $\boldsymbol{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. [201], we will examine the behaviour
of the following entropic quantities. First, we will look at the relative
entropy between $\boldsymbol{p}^{(n)}$ and the thermal state of the system
$\boldsymbol{\gamma}$,

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

which is a thermodynamic monotone, as it decreases under (Markovian) thermal
operations, and is directly related to the non-equilibrium free energy [90].
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(\boldsymbol{p}^{(n)}_{SM}\|\boldsymbol{p}^{(n)}\otimes\boldsymbol{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. 9a and
9e, 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. 9c and 9g. 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. 9b and 9f, 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.
9d and 9e.

###  Section 4 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 2, where it was defined in Eq.
(10). 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 in Fig. 10a the diagram of $\tilde{\mathcal{P}}$ for
$d=7$.

Figure 10: _Alternative protocols_. Three different protocols are used to
achieve equivalent $\beta$-swap approximations. Each thermalisation is
represented by a point on a plane, with an arrow indicating the transition
from one thermalisation to the next.

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. Now, let us look at the two
algorithms depicted in Fig. 10b and Fig. 10c. The action 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
one. 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:

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
$\boldsymbol{p}\otimes\boldsymbol{\gamma}_{d}$ result in the same state as
$\tilde{\mathcal{P}}(\boldsymbol{p}\otimes\boldsymbol{\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 $\boldsymbol{p}$ with $\beta=0$) is plotted in
Fig. 11.

Figure 11: _Algorithm convergence analysis_. Convergence to the target state
$\boldsymbol{q}=(0,1)$ of different algorithms acting on the state
$\boldsymbol{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{\boldsymbol{p}-\boldsymbol{q}}_{1}$.

###  Section 5 Derivation of the results

####  Subsection 1 Regularised incomplete beta function

The content of this section is based on Refs. [209, 210].

_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,$ (43)

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)}.$ (44)

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

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

[-5.0cm]  _Regularised beta function_. Plots of the regularised incomplete
beta function as a function of (a) $a$ and (b) $b$.  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)}.$ (46)

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

Throughout this chapter, 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}.$ (47)

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

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

Moreover,

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

If $a=1$, then $I_{x}(1,b)=1-(1-x)^{b}$ and Eq. (47) 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).$ (50)

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

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

and the relation for an equal argument,

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

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)},$ (53a)
$\displaystyle I_{x}(a,b)$
$\displaystyle=I_{x}(a,b+1)-\frac{x^{a}(1-x)^{b}}{aB(a,b)},$ (53b)
$\displaystyle I_{x}(a,b)$
$\displaystyle=I_{x}(a+1,b-1)+\frac{x^{a}(1-x)^{b-1}}{aB(a,b)},$ (53c)
$\displaystyle I_{x}(a,b)$
$\displaystyle=I_{x}(a-1,b+1)-\frac{x^{a-1}(1-x)^{b}}{aB(a,b)}.$ (53d)

_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)!.$ (54)

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}.$
(55)

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).$ (56)

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).$ (57)

_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)=\Gamma(a+b)x^{a}(1-x)^{b-1}$
$\displaystyle\Biggr{[}\sum_{k=0}^{n-1}\frac{1}{\Gamma(a+k+1)\Gamma(b-k)}\quantity(\frac{x}{1-x})^{k}$
$\displaystyle\hskip 34.14322pt+O\quantity(\frac{1}{\Gamma(a+n+1)})\Biggr{]}.$
(58)

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. (1) 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)}],$
(59)

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},$
(60)

where we have used definitions of the beta function and the gamma function for
integer arguments. Using Eq. (49), 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)],$
(61)

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

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

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

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.$ (63)

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=\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\hskip
48.36958pt\times\int_{0}^{x}dtN(1-t)^{N-1}\sum_{i=1}^{N}\binom{2N-i}{N-i}t^{N-i}$
$\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=\int_{0}^{x}dt\frac{N}{(1-t)^{2}}\Bigg{(}I_{1-x}(0,N+1)-I_{1-x}(N,N+1)\Bigg{)}$
$\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=N\int_{0}^{x}dt\frac{1}{(1-t)^{2}}=N\frac{x}{1-x}.$ (64)

####  Subsection 2 Proofs of Lemma 2.1 and Theorem 2.4

To prove Lemma 2.1 and Theorem 2.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.

_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 $\boldsymbol{p}=(b,c)$,
together with a memory system, described by a trivial Hamiltonian $H_{M}=0$
and prepared in a maximally mixed state
$\boldsymbol{\eta}_{M}=\quantity(1,\ldots,1)/N$. The joint state
$\boldsymbol{r}:=\boldsymbol{p}\otimes\boldsymbol{\eta}_{M}$ of the composite
system is then given by

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

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

$\displaystyle\boldsymbol{\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}).$
(66)

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}}.$ (67)

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(\boldsymbol{r}^{(k-1)})\equiv\boldsymbol{r}^{(k)}\\!=\\!\quantity(b_{1}^{(k)},\ldots,b_{N}^{(k)}\,|\,c_{1}^{(N)},\ldots,c^{(N)}_{k},c,\ldots,c),$
(68)

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}),$
(69) $\displaystyle c^{(N)}_{j}$
$\displaystyle=\Gamma_{21}^{\,N}c+\Gamma_{21}^{\,N+1}\sum_{i=1}^{N}b^{(j-1)}_{i}\Gamma^{\,-i}_{21},$
(70)

with $b^{(0)}_{j}=b$. Equations (69) and (70) 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},$
(71) $\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}.$
(72)

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

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

_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 2.1 and Theorem 2.4, the next step consists of decoupling the main
system from memory using a full thermalisation $\mathcal{T}$:

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

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})}),$ (75a)
$\displaystyle\lim_{N\to\infty}\,\frac{\sum_{i=1}^{N}c_{i}^{(N)}}{N}$
$\displaystyle=be^{-\beta(E_{2}-E_{1})}.$ (75b)

To prove the limits in Eqs. (75a) and (75b), 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.$
(76)

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}.$
(77)

We will now examine the non-trivial term on the right hand side of Eq. (77).
We start by using Eq. (49) to re-write the first term of $c^{(N)}_{j}$
appearing in Eq. (72) 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),$ (78)

where we also used Eq. (51) to invert the arguments of the regularised beta
function. Thus, using (63), 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.$
(79)

Taking into account the second term of $c^{(N)}_{j}$ appearing in Eq. (72),
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}},$
(80)

and therefore obtain that

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

Substituting this result to Eq. (77), 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})}].$
(82)

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

_Finite memory convergence rates_

To complete the proofs of Lemma 2.1 and Theorem 2.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},$
(83a) $\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},\\!$
(83b)

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

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

Next, we will asymptotically expand Eq. (83a). 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. (53c), (52) and (55):

$\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}}].$
(85)

Now we expand Eq. (85) up to second order using Eq. (1). 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. (56) and
Eq. (55). Finally, using Stirling’s approximation [Eq. (57)], 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\hskip 99.58464pt+1+o\left(N^{-1}\right)\Bigg{]}[1+O(N^{-1})].$
(86)

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].$ (87)

The infinite-temperature limit, $\beta=0$, is similarly analysed by using Eq.
(85) 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).$ (88)

The other half of estimating the convergence rate stems from considering Eq.
(83b). 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}}.$ (89)

Again, we expand Eq. (89) up to second order using Eq. (1). Then, the gamma
functions appearing in the expansions are simplified using Eq. (56), the
remaining ones are simplified by using Eq. (55). Finally, using Stirling’s
approximation [Eq. (57)], 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}}\quantity(1-\frac{4\Gamma_{21}\Gamma_{12}}{2(N+1)(1-4\Gamma_{21}\Gamma_{12})})\Biggr{]}$
$\displaystyle\hskip
128.0374pt+\frac{(4\Gamma_{21}\Gamma_{12})^{N}}{2\Gamma_{12}}\left[\frac{1}{\sqrt{\pi
N}}+o(N^{-3/2})\right].$ (90)

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}},$ (91)

where $\simeq$ hides the $o$-terms. The infinite-temperature limit, $\beta=0$,
is obtained by using Eq. (89) 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)$ (92)

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

As a final step to prove Lemma 2.1 and Theorem 2.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.
(84), yielding

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

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

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

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

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

$\displaystyle\delta\quantity(\Pi^{\beta}_{12}\boldsymbol{p},\,\boldsymbol{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)}+o(N^{-3/2})\Bigg{]}.$ (95)

These prove Lemma 2.1 and Theorem 2.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}$. ∎

_Strengthening Corollary 2.3_

Corollary 2.3 deals with a very general approach to bounding the distance
between the target state $\boldsymbol{p}^{\pi}$ and its approximation obtained
from $\boldsymbol{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 5.1.

Consider states $\boldsymbol{p}$ and $\boldsymbol{q}\in
C_{+}^{TO}(\boldsymbol{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}(\boldsymbol{p}\otimes\boldsymbol{\eta}_{M})=\boldsymbol{q}^{\prime}\otimes\boldsymbol{\eta}_{M},$
(96)

with

$\delta(\boldsymbol{q}^{\prime},\boldsymbol{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,$
(97)

where $\quantity{i_{1},\ldots,i_{d^{\prime}}}\subset\quantity{1,\ldots,d}$ is
a subset of indices neighbouring in the $\beta$-order,
$\pi_{\boldsymbol{p}}(i_{j})+1=\pi_{\boldsymbol{p}}(i_{j+1})$, such that
$\Pi_{\boldsymbol{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
$\boldsymbol{p}^{\pi^{*}}$ such that the $\beta$-order
$\Pi_{\boldsymbol{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. (19) from Theorem 2.2, one
finds that

$\delta\left(\boldsymbol{p}^{\pi^{*}},\boldsymbol{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}),$
(98)

where for convenience we used
$\mathcal{P}^{\Pi^{*}}(\boldsymbol{p}\otimes\eta_{M})=\boldsymbol{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}(\boldsymbol{p}\otimes\eta_{M})=\boldsymbol{r}\otimes\eta_{M}$
we find that

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

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

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

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

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

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

The bound presented in Corollary 5.1 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_{\boldsymbol{q}}$ into
neighbour transpositions. Finally, we point out that, in agreement with Eq.
(100), there will exist such states $\boldsymbol{q}$ that are attainable
exactly, and moreover, their volume will increase together with the size of
memory $N$.

_Proof of Theorem 2.5_

To prove Theorem 2.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 $\boldsymbol{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 $\boldsymbol{\eta}_{M}=(1/N,...,1/N)$. The joint state of the composite
system, $\boldsymbol{r}:=\boldsymbol{p}\otimes\boldsymbol{\eta}_{M}$, is then
given by

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

and the joint thermal state is given by

$\boldsymbol{\Gamma}=\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}}],$ (103)

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 $\boldsymbol{r}^{(N)}$ is then given by

$\displaystyle\\!\\!\\!\widetilde{{\cal P}}^{\Pi}(\boldsymbol{r})$
$\displaystyle\equiv\boldsymbol{r}^{(N)}=\boldsymbol{q}\otimes\boldsymbol{\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}]},$
(104)

where $\Pi=\Pi_{13}\Pi_{23}$. Note that due to probability conservation,
characterising the second and third entries of Eq. (104) 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. (71)
(with $\Gamma_{12}$ and $\Gamma_{21}$ replaced by $\Gamma_{23}$ and
$\Gamma_{32}$, respectively). The other two entries are obtained in a similar
way as Eqs. (69)-(70). 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}+\Gamma^{N}_{31}\sum_{l=0}^{N-1}\Gamma_{13}^{l}\binom{N+l-1}{l}c_{k-l},$
(105)

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}.$
(106)

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

$\boldsymbol{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].$
(107)

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}}),$ (108a)
$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{i=1}^{N}c^{(N)}_{i}$
$\displaystyle=a\frac{\Gamma_{31}}{\Gamma_{13}}.$ (108b)

_Proof of limit ( 108a) _

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}.$
(109)

Comparing Eqs. (109) and the right-hand side of (108a), we see that in order
to prove Eq. (108a), 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$ (110)

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}}).$ (111)

We begin by proving Eq. (2). 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}.$
(112)

We can evaluate the second sum in Eq. (112) 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)].$ (113)

Thus, substituting Eq. (2) into Eq. (112), 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).$ (114)

Using Eq. (63), we conclude that the second term in Eq. (114) 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,$ (115)

so that we have proved Eq. (2).

To prove Eq. (2), 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}.$
(116)

Applying Eq. (49) to transform the first term of Eq. (2) 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.$ (117)

Next, we consider the second term in Eq. (2), 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}},$ (118)

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}}.$ (119)

Therefore, combining the above with Eq. (115), 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}}).$
(120)

_Proof of limit ( 108b)_

As before, in order to prove Eq. (108b), 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}\binom{N+j-1-i}{j-1}\Gamma_{31}^{N+1-i}+\Gamma^{N}_{31}\sum_{l=0}^{N-1}\Gamma_{13}^{l}\binom{N+l-1}{l}c_{j-l},$
(121)

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}.$
(122)

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}},$ (123)

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.$ (124)

Let us start by proving Eq. (123). 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).$ (125)

Using Eq. (64), we obtain

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

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}}.$ (127)

The final step is to show that the remaining limit from Eq. (124) 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.$
(128)

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. (122)
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_{32}(N,k)\leq
1,$ (129)

where we could bound by one because the incomplete beta function is 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}}\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),$ (130)

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. (122) 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}}.$ (131)

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. (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_{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_{31}(N,N)\sum_{k=1}^{N}\Gamma^{k}_{23}\Gamma_{23}^{-1-k}\quantity[1-I_{32}(N,k+1)]$
$\displaystyle=\frac{\Gamma_{32}}{\Gamma_{23}}\frac{1}{N}I_{31}(N,N)\quantity(\frac{1}{N}\sum_{k=1}^{N}I_{23}(k+1,N)).$
(132)

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_{23}(k+1,N)\leq\frac{1}{N}\sum_{k=1}^{N}1=1.$
(133)

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,$ (134)

we conclude that the resulting limit is zero. Therefore,

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

###  Section 6 Concluding remarks

In this chapter, 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 [211, 212, 213]. On the other
hand, it may be seen as far less general than the approach taken in Ref.
[214], where our work would correspond to a step towards simulating arbitrary
evolution with Markov-Stinespring curves.

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. However, the
method can be used as well to investigate other thermodynamic protocols, such
as information erasure or thermodynamically free encoding of information
[141]. 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.2. This can be built upon the proofs for
$\beta$-swaps (Theorem 2.4) and $\beta$-$3$-cycles (Theorem 2.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 [215, 216]. 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 Section 4). Such explorations could focus
on enabling parallelisation of certain steps by expanding available memory,
thereby illustrating the space-time trade-off.

##  Chapter 6 Fluctuation-dissipation relations for thermodynamic distillation
processes

Almost two centuries ago, Robert Brown observed that pollen seeds immersed in
water move randomly in erratic motion [217]. It was not until the 1905 papers
by Einstein and Smoluchowski [218, 219] that scientists understood that this
“Brownian” motion is induced by the bombardment of pollen particles by water
molecules. Crucially, by noting that these collisions would also create
friction for the particle being pulled through the fluid, Einstein realised
that the two processes, fluctuations of particle’s position and dissipation of
its energy, have the same origin and thus must be related. Over the years,
physicists generalised and formalised this observation into fundamental
fluctuation-dissipation relations describing the behavior of systems driven
out of equilibrium [13, 220].

Now, it is well known that near-equilibrium, linear response theory provides a
general proof of the fluctuation-dissipation theorem, which states that the
response of a given system when subject to an external perturbation is
expressed in terms of the fluctuation properties of the system in thermal
equilibrium [13]. The theoretical description underlying the fluctuation-
dissipation relations is usually expressed in terms of the stochastic
character of thermodynamic variables. This approach is strongly motivated
since it is experimentally viable [221, 222].

As introduced in the Chapter 3, the resource theory of thermodynamics aims to
go beyond the thermodynamic limit and the assumption of equilibrium. It is
often presented as an extension of statistical mechanics to scenarios with
large fluctuations, referred to as single-shot statistical mechanics. [223,
224]. A natural question is then whether fluctuation-dissipation relations are
present in such a resource-theoretic description. Although important insights
have been obtained in trying to connect the information-theoretic and
fluctuation theorem approaches [225, 226], they have, so far, not been
explicitly related to dissipation. The tools required for the analysis of free
energy dissipation in a resource-theoretic framework were developed in Refs.
[227, 103, 130, 228], where the authors investigated irreversibility of
thermodynamic processes due to finite-size effects. However, the relation
between fluctuations and actual dissipation was not derived and, moreover,
these results were obtained for quasi-classical case of energy-incoherent
states. Thus, they were not able to account for quantum effects that come into
play when dealing with even smaller systems, when fluctuations around
thermodynamic averages are no longer just thermal in their origin.

This chapter pushes towards a genuinely quantum framework characterising
optimal thermodynamic state transformations and links fluctuations with free
energy dissipation. We investigate a special case of state interconversion
processes known as _thermodynamic distillations_. These are thermodynamic
processes in which a given initial quantum system is transformed, with some
transformation error, to a pure energy eigenstate of the final system. In
particular, we focus on the initial system consisting of asymptotically many
non-interacting subsystems that are either energy-incoherent and non-identical
(in different states and with different Hamiltonians), or pure and identical.
Within this setting, our main results are given by two theorems. The first one
yields the optimal transformation error as a function of the free energy
difference between the initial and target states, and the free energy
fluctuations in the initial state. This can be seen as an extension of
previously derived results on optimal thermodynamic state transformations
[227, 103, 130, 228]. The second theorem provides a precise relation between
the free energy fluctuations of the initial state and the minimal amount of
free energy dissipated in the optimal thermodynamic distillation process. It
is conceptually novel and does not form an extension of previously known
results, and as such it constitutes our main contribution. Note that the
second theorem employs the first one as one of the building blocks.

Our results allow us for a rigorous study of important thermodynamic
protocols. First of all, we extend the analysis of work extraction to the
regime of not necessarily identical incoherent states, as well as to pure
states. By directly applying our main results, we obtain a second-order
asymptotic expression for the optimal transformation error while extracting a
given amount of work from the initial system. Moreover, we also verify the
accuracy of the obtained expression by comparing it with the numerically
optimised work extraction process. As a second application, we analyse the
optimal energetic cost of erasing $N$ independent bits prepared in arbitrary
states. In this case, we obtained the optimal transformation error for the
erasure process as a function of invested work. The last application we
consider is the optimal thermodynamically-free communication scheme, i.e., the
optimal encoding of information into a quantum system without using any extra
thermodynamic resources. Applying our theorems gives us the optimal number of
messages that can be encoded into a quantum system in a thermodynamically
free-way, which we show to be directly related to the non-equilibrium free
energy of the system. This result can be interpreted as the inverse of the
Szilard engine, as in this process we use the ability to perform work to
encode information. Furthermore, our results connect the fluctuations of free