# Heat-Bath Algorithmic Cooling with optimal thermalization strategies

Álvaro M. Alhambra  Perimeter Institute for Theoretical Physics, Waterloo, ON
N2L 2Y5, Canada  Matteo Lostaglio  ICFO-Institut de Ciencies Fotoniques, The
Barcelona Institute of Science and Technology, Castelldefels (Barcelona),
08860, Spain  Christopher Perry  QMATH, Department of Mathematical Sciences,
University of Copenhagen, Universitetsparken 5, 2100 Copenhagen, Denmark

###### Abstract

Heat-Bath Algorithmic Cooling is a set of techniques for producing highly pure
quantum systems by utilizing a surrounding heat-bath and unitary interactions.
These techniques originally used the thermal environment only to fully
thermalize ancillas at the environment temperature. Here we extend HBAC
protocols by optimizing over the thermalization strategy. We find, for any
$d$-dimensional system in an arbitrary initial state, provably optimal cooling
protocols with surprisingly simple structure and exponential convergence to
the ground state. Compared to the standard ones, these schemes can use fewer
or no ancillas and exploit memory effects to enhance cooling. We verify that
the optimal protocols are robusts to various deviations from the ideal
scenario. For a single target qubit, the optimal protocol can be well
approximated with a Jaynes-Cummings interaction between the system and a
single thermal bosonic mode for a wide range of environmental temperatures.
This admits an experimental implementation close to the setup of a micromaser,
with a performance competitive with leading proposals in the literature. The
proposed protocol provides an experimental setup that illustrates how non-
Markovianity can be harnessed to improve cooling. On the technical side we 1.
introduce a new class of states called _maximally active states_ and discuss
their thermodynamic significance in terms of optimal unitary control, 2.
introduce a new set of thermodynamic processes, called _$\beta$
-permutations_, whose access is sufficient to simulate a generic
thermalization process, 3. show how to use abstract toolbox developed within
the resource theory approach to thermodynamics to perform challenging
optimizations, while combining it with open quantum system dynamics tools to
approximate optimal solutions within physically realistic setups.

Cooling is a central problem in quantum physics and in realizing technologies
for quantum information processing. The ability to produce a set of highly
pure, ‘cold’, quantum states is vital for the construction of a quantum
computer [DiVincenzo(2000)]. More generally, the observation of quantum
effects often requires cooling and, as such, many techniques have been
developed to cool systems efficiently in platforms ranging from cavity
optomechanics [Liu et al.(2013)Liu, Xiao, Luan, and Wong] to NMR [Schulman and
Vazirani(1999), Boykin et al.(2002)Boykin, Mor, Roychowdhury, Vatan, and
Vrijen], ion traps [Eschner et al.(2003)Eschner, Morigi, Schmidt-Kaler, and
Blatt] and superconducting qubits [Valenzuela et al.(2006)Valenzuela, Oliver,
Berns, Berggren, Levitov, and Orlando].

Here we use powerful techniques, developed within the resource theory approach
to thermodynamics [Lostaglio(2018)], to greatly extend an important class of
cooling algorithms known as Heat-Bath Algorithmic Cooling (HBAC) [Schulman and
Vazirani(1999), Boykin et al.(2002)Boykin, Mor, Roychowdhury, Vatan, and
Vrijen]. The goal of these is to maximize the purity of a target system $S$ in
a given number of cooling rounds. Each round of the algorithm starts with a
unitary applied to the target together with several auxiliary systems $A$
initialized in a thermal state, with the aim of pumping entropy away from the
target. Next, the auxiliary systems are re-thermalized through coupling with a
heat-bath, before the entire process is repeated in the next round. The
asymptotically optimal protocol of this form (in terms of the purity reached
in infinitely many rounds) is the Partner Pairing Algorithm (PPA), introduced
in [Schulman et al.(2005)Schulman, Mor, and Weinstein], whose asymptotic
performance has been recently derived for a single target qubit starting in a
maximally mixed state [Rodríguez-Briones and Laflamme(2016), Raeisi and
Mosca(2015)].

Here, we generalize HBAC in the following sense: instead of using the thermal
environment only to refresh the auxiliary systems through a _complete_
thermalization, we allow strategies involving an _incomplete_ thermalization
of system and ancillas. In particular, we optimize the protocol over _every_
‘dephasing thermalization’, that is any quantum map on $SA$ which 1. leaves
the thermal state on $SA$ fixed and 2. dephases the input state in the energy
eigenbasis. This generalizes the ‘rethermalization’ used in previous
protocols, and thus extends HBAC to a larger set that will be called Extended
Heat-Bath Algorithmic Cooling (xHBAC). For example, the recent protocol
introduced in Ref. [Rodriguez-Briones et al.(2017)Rodriguez-Briones, Li, Peng,
Mor, Weinstein, and Laflamme], which proposes the use of the heat bath to
implement a non-local thermalization ‘state-reset’ (SR) process related to the
Nuclear Overhauser Effect [Overhauser(1953)], can already be seen as a (non-
optimal) protocol within our extended family of xHBAC.

For _every_ finite dimensional target state in an _arbitrary_ initial state,
we optimize the cooling performance over _every_ xHBAC for _any_ given number
of rounds. We give an analytical form for the optimal cooling operations,
uncovering their elegant structure, and show that the ground state population
goes to $1$ exponentially fast in the number of rounds. This opens up a new
avenue in the experimental realization of algorithmic cooling schemes, one
requiring control over fewer ancillas, but better control of the interaction
with the environment.

Our results suggest that xHBAC schemes provide new cooling protocols in
practically relevant settings. We show that, for a single qubit target and no
ancillas, the optimal protocol in xHBAC can be approximated by coupling the
qubit by a Jaynes-Cummings (JC) interaction to a single bosonic thermal mode,
itself weakly coupled to the external bath. The memory effects present in the
JC interaction are crucial in approximating the theoretical optimal cooling.
The performance of the ideal cooling protocol is robust to certain kinds of
noise and imperfections, and it outperforms HBAC schemes with a small number
of auxiliary qubits. This makes it, in our opinion, the most promising
proposal for an experimental demonstration of cooling through xHBAC, as well
as an experimental demonstration that non-Markovianity can be harnessed to
improve cooling.

##  1 Results and discussion

###  1.1 A general cooling theorem

A general xHBAC will consist of a number of rounds and manipulate two types of
systems, the target system $S$ to be purified and the auxiliary systems $A$.
Furthermore, we will assume we can access a thermal environment at inverse
temperature $\beta=(kT)^{-1}$, with $k$ Boltzmann’s constant. We denote by
$H_{S}=\sum_{i=0}^{d-1}E_{i}|i\rangle\langle i|$, $E_{0}\leq\dots\leq
E_{d-1}$, the Hamiltonian of $S$ and by $\rho^{\left(k\right)}_{S}$ the state
after round $k$. Also, we denote by $H_{A}$ the Hamiltonian of the auxiliary
systems, initially in state $\rho_{A}$, and by $\tau_{X}=e^{-\beta
H_{X}}/\textnormal{Tr}{[e^{-\beta H_{X}}]}$ the thermal state on $X=S,A$. A
protocol is made of $k$ rounds, each allowing for (Fig. 1):

  1. 1.

Unitary. Any unitary $U^{(k)}$ applied to $SA$.

  2. 2.

Dephasing thermalization. Any $\Lambda^{(k)}$ applied to $SA$, where
$\Lambda^{(k)}$ is any quantum map such that _i)_
$\Lambda^{(k)}(\tau_{S}\otimes\tau_{A})=\tau_{S}\otimes\tau_{A}$ (thermal
fixed point); _ii)_ If $|E\rangle$, $|E^{\prime}\rangle$ are states of
distinct energy on $SA$, $\langle
E|\Lambda^{(k)}(\rho_{SA})|E^{\prime}\rangle=0$ (dephasing).

At the end of the round, a refreshing stage returns the auxiliary systems $A$
to their original state $\rho_{A}$. We denote the set of protocols whose
rounds have this general form by $\mathcal{P}_{\rho_{A}}$. Typically
$\rho_{A}=\tau_{A}$, and this operation is simply a specific kind of dephasing
thermalization. However, more generally, $\rho_{A}$ may also be the output of
some previous cooling algorithm, subsequently used as an auxiliary system.
Note that $A$ is distinguished by the rest of the environment in that it is
under complete unitary control.

Some comments about xHBAC protocols are in order. First, if we skipped the
dephasing thermalizations and $\rho_{A}=\tau_{A}$, we would get back to a
standard HBAC scheme. xHBAC protocols include thermalization of subsets of
energy levels, as in the mentioned SR protocol of Ref. [Rodriguez-Briones et
al.(2017)Rodriguez-Briones, Li, Peng, Mor, Weinstein, and Laflamme]. But they
are by no means limited to these strategies. Second, for every $k$ we wish to
find a protocol maximizing the ground state population
$p^{\left(k\right)}_{0}=\langle 0|\rho^{\left(k\right)}_{S}|0\rangle$ over all
sequences of $k$-round operations. Note that often the literature on cooling
is restricted to finding asymptotically optimal protocols (e.g., the PPA
[Schulman et al.(2005)Schulman, Mor, and Weinstein]), but here we find optimal
$k$ rounds protocols for every $k$. Finally, we note in passing that within
this work we will not make use of auxiliary ‘scratch qubits’ [Rodríguez-
Briones and Laflamme(2016), Raeisi and Mosca(2015)], i.e. $S$ itself is the
target system to be cooled.

To give the analytical form of the optimal protocols we need to introduce two
notions. First, that of _maximally active_ states. Given a state $\rho$ with
Hamiltonian $H$, the maximally active state $\hat{\rho}$ is formed by
diagonalizing $\rho$ in the energy eigenbasis and ordering the eigenvalues in
increasing order with respect to energy. That is, $\hat{\rho}$ is the most
energetic state in the unitary orbit of $\rho$.

Figure 1: A generic xHBAC protocol.

Second, we define the _thermal polytope_ in the following way. Given an
initial state $\rho$ with populations in the energy eigenbasis p, it is the
set of populations $\boldsymbol{p}^{\prime}$ of $\Lambda(\rho)$, with
$\Lambda$ an arbitrary dephasing thermalization. Without loss of generality,
we assume that every energy eigenspace has been diagonalized by energy
preserving unitaries. The thermal polytope is convex and has a finite number
of extremal points (Lemma 12, [Lostaglio et al.(2018)Lostaglio, Alhambra, and
Perry]). It is particularly useful in thermodynamics to classify out of
equilibrium distributions according to their $\beta$-order [Horodecki and
Oppenheim(2013)]. Given a vector $\boldsymbol{p}$ and a correspondent
Hamiltonian with energy levels $E_{i}$, the $\beta-$order of $\boldsymbol{p}$
is the permutation of $\pi$ of $0,1,2,\dots$ such that the Gibbs-rescaled
population is sorted in non-increasing order: $p_{\pi(0)}e^{\beta
E_{\pi(0)}}\geq p_{\pi(1)}e^{\beta E_{\pi(1)}}\geq\dots$.

Crucially, for any dimension $d$, we provide an explicit construction of a set
of dephasing thermalizations, denoted $\Lambda^{(\pi,\alpha)}$ and called
$\beta$-_permutations_ , able to map any given $\boldsymbol{p}$ to every
extremal point of its thermal polytope (Methods, Sec. 2.1 and Fig. 2). $\pi$
and $\alpha$ are two vectors of integers, each ranging over all permutations
of $\\{0,...,(d-1)(r-1)\\}$, if $r$ is the dimension of $A$. Their
significance is the following. If we are given an initial state
$\boldsymbol{p}$ with $\beta$-order $\pi$, the state
$\boldsymbol{q}^{\alpha}:=\Lambda^{(\pi,\alpha)}(\boldsymbol{p})$ is the
‘optimal’ state among all those with $\beta$-order $\alpha$. More precisely,
there is no other state in the thermal polytope of $\boldsymbol{p}$ which has
$\beta$-order $\alpha$ and can be transformed into $\boldsymbol{q}^{\alpha}$
by dephasing thermalizations. This is in analogy with the role of permutations
under doubly stochastic maps. In fact, in the limit $\beta\rightarrow+\infty$,
$\beta$-permutations are permutations and the $\beta$-ordering is the standard
sorting.

A particularly important $\beta$-permutation, denoted by $\beta^{\rm opt}$, is
the one that maximizes the ground state population of $S$ among all dephasing
thermalizations and, furthermore, achieves the largest partial sums
$\sum_{i=0}^{l}p^{(k)}_{i}$, $l=0,...,d-1$, $p^{(k)}_{i}=\langle
i|\rho^{(k)}_{S}|i\rangle$. $\beta^{\rm opt}$ is constructed as follows:
$\beta^{\rm opt}=\Lambda^{(\pi,\alpha)}$ with $\pi$ an ordering of energy
levels of $SA$ such that, if $q^{(k)}_{m}$ are populations of $SA$ at round
$k$ and $E^{SA}_{m}$ are the energies of $H_{S}+H_{A}$,
$q^{(k)}_{\pi(m)}e^{\beta E^{SA}_{\pi(m)}}$ are sorted in non-increasing order
in $m$; and

$\displaystyle\alpha=$
$\displaystyle\big{\\{}\left(0,r-1\right),\left(0,r-2\right),\dots,\left(0,0\right),$
$\displaystyle\left(1,r-1\right),\left(1,r-2\right),\dots,\left(1,0\right),\dots,\left(d-1,0\right)\big{\\}}.$
(1)

This identifies an optimal dephasing thermalization. But what is the optimal
unitary control that we need to apply? The following lemma answers this
question:

######  Lemma 1.

Given a state $\rho$ with population $\boldsymbol{p}$, the thermal polytope of
the maximally active state $\hat{\rho}$ contains the thermal polytope of
$U\rho U^{\dagger}$ for every unitary $U$.

The proof can be found in the Methods, Sec. 2.3. With these concepts in place,
we can give the optimal cooling protocol for any given set of auxiliary
systems in $A$.

######  Theorem 1.

Assume $S$ is a $d$-level system with $E_{d-1}>E_{0}$ and $\beta>0$. Without
loss generality (by an initial diagonalizing unitary) we can take the initial
state of $S$ to be diagonal in the energy eigenbasis with
$p^{\left(0\right)}_{0}\geq\dots\geq p^{\left(0\right)}_{d-1}$. Then, for a
given auxiliary state $\rho_{A}$, the optimal cooling protocol in
$\mathcal{P}_{\rho_{A}}$ is such that in each round $k$:

  1. 1.

The unitary mapping $\rho^{\left(k-1\right)}_{S}\otimes\rho_{A}$ to the
corresponding maximally active state is applied to $SA$.

  2. 2.

$\beta^{\rm opt}$ is applied to $SA$.

The optimal protocol achieves $p_{0}^{(k)}\rightarrow 1$ at least
exponentially fast in $k$, even with no ancilla.

The intuition behind this protocol is simple (see Fig. 2). In a single round,
the unitary maximizes the amount of energy in $SA$, in accordance to Lemma 1.
Next, the optimal $\beta$-permutation is applied. For the proof, see Methods
(Sec. 2.4).

Figure 2: An optimal cooling round on a $d=3$ system. Map the initial state
after round $k-1$ to its maximally active state by a unitary $U$ (continuous
black arrow). The thermal polytope of the maximally active state then
describes all possible final states achievable by arbitrary thermalizations.
$\beta$-permutations map to the extremal points of these polytope (dotted red
arrows), and $\beta^{\rm opt}$ achieves the one closest to the ground state
(dash-dot red arrow). The distance to the ground state decreases exponentially
fast in the number $k$ of rounds.

In analyzing the performance of the protocols $\mathcal{P}_{\rho_{A}}$ it is
important to keep in mind the cost of preparing and controlling the auxiliary
systems $A$, especially if $\rho_{A}\neq\tau_{A}$. Remarkably, however, the
optimal xHBAC protocol that uses no auxiliary systems $A$ still has
$p^{(k)}_{0}\rightarrow 1$ exponentially in $k$. Furthermore, such protocol
has the further advantage that the cooling operations do not change with the
round $k$. Also note that the initial state-dependent unitary can always be
replaced by a complete thermalization while maintaining the same cooling
scaling, so that no prior knowledge of $\rho^{(0)}_{S}$ is needed to unlock a
strong cooling performance. Finally, in the Methods section 2.5, we construct
for any $d$ an explicit protocol that not only uses no ancillas, and can be
realized by a sequence of two level interactions: first, a unitary swaps the
populations of the ground and most excited states; then, a sequence of two-
level dephasing thermalizations between energies $(d-1,d-2)$, $(d-2,d-3)$,
$\dots$, $(1,0)$ is performed, each maximising the net population transfer
$i\mapsto i-1$ (these are $\beta$-swaps, discussed in more detail in the next
section). Every $d-1$ repetitions, the performance takes an elegant form:

$p^{\left(k\left(d-1\right)\right)}_{0}=1-e^{-k\beta\left(E_{d-1}-E_{0}\right)}\left(1-p^{\left(0\right)}_{0}\right).$
(2)

###  1.2 The qubit case

To show how the general results can be used in a specific case, let us analyze
in detail the single qubit case with no auxiliary systems $A$ (denoted
$\mathcal{P}_{\emptyset}$) and energy gap $E$. The only nontrivial
$\beta$-permutation is $\Lambda_{\beta}=\Lambda^{(\pi,\alpha)}$ with
$\pi=\\{1,0\\}$ and $\alpha=\\{0,1\\}$, which induces transition probabilities
$\langle 0|\Lambda_{\beta}(|1\rangle\langle 1|)|0\rangle=1$, $\langle
1|\Lambda_{\beta}(|0\rangle\langle 0|)|1\rangle=e^{-\beta E}$. This is the
$\beta$-swap introduced in [Lostaglio et al.(2018)Lostaglio, Alhambra, and
Perry], which can be realized by the dephasing thermalization

$\displaystyle\Lambda_{\beta}(\rho_{S})=$
$\displaystyle\sigma_{-}\rho_{S}\sigma_{+}+e^{-\beta
E}\sigma_{+}\rho_{S}\sigma_{-}$ (3) $\displaystyle+(1-e^{-\beta
E})\sigma_{-}\sigma_{+}\rho_{S}\sigma_{-}\sigma_{+}$

where $\sigma_{+}=|1\rangle\langle 0|$, $\sigma_{-}=\sigma_{+}^{\dagger}$.
Since the unitary mapping the state to its correspondent maximally active
state is the Pauli $X$ unitary, Theorem 1 reads as follows:

######  Corollary 1.

Assume $E>0$, $\beta>0$. Without loss of generality (by making use of an
initial diagonalizing unitary), we can take the initial state of the system to
be diagonal in the energy basis with $p_{0}^{(0)}\geq p_{1}^{(0)}$. The
optimal cooling protocol in $\mathcal{P}_{\emptyset}$ is such that in each
round $k$:

  1. 1.

The Pauli $X$ unitary is applied to $S$.

  2. 2.

The $\beta$-swap $\Lambda_{\beta}$ is applied to $S$.

The population of the ground state after round $k$ is:

$p^{\left(k\right)}_{0}=1-e^{-k\beta E}\left(1-p^{\left(0\right)}_{0}\right),$
(4)

and $p_{0}^{\left(k\right)}\rightarrow 1$ as $k\rightarrow\infty$.

The performance can be obtained by direct computation, or by setting $d=2$ in
Eq. (2). Note the simple structure of the optimal protocol, in particular the
fact that the same operation is applied iteratively at each round.
Furthermore, note that $\beta$-swaps cannot be realized by Markovian
interactions with the thermal environment. In fact, if we optimized over
Markovian interactions only, the optimal protocol with no auxiliary systems
would be the trivial thermalization at the environment temperature
$\rho_{S}\mapsto\tau_{S}$. Then, at best we would achieve a ground state
population $1/(1+e^{-\beta E})$. The proof is given in the Methods (Sec. 2.6).
This shows that memory effects can be used to great advantage in cooling
scenarios. In particular, differently from standard HBAC protocols, our
protocol achieves exponential convergence to the ground state. The ‘price’ we
pay is the need for greater control over the thermalization steps; however, we
need no auxiliary systems in the unitary stage and the protocol iterates the
same transformation at every round, which simplifies the implementation.

We presented an optimal protocol for qubits, but how is it realized by an
explicit interaction and environment? Here we answer this question. It is
known that an infinite dimensional environment is necessary to increase the
purity of the target system to 1 in the absence of initial system-environment
correlations [Ticozzi and Viola(2014), Silva et al.(2016)Silva, Manzano,
Skrzypczyk, and Brunner]; however, we do not need a complex environment: it
was shown in [Lostaglio et al.(2018)Lostaglio, Alhambra, and Perry] that a
single bosonic mode in a thermal state $\tau_{B}=\left(1-e^{-\beta
E}\right)\sum_{n=0}^{\infty}e^{-n\beta E}|n\rangle\langle n|$ suffices to
implement the required $\beta$-swap. Specifically, the unitary needed is

$U^{\beta}_{SB}=|0,0\rangle\langle
0,0|+\sum_{n=1}^{\infty}\left(|0,n\rangle\langle 1,n-1|+|1,n-1\rangle\langle
0,n|\right).$ (5)

What is perhaps more remarkable is that, as we prove, we do not need to
rethermalize (or otherwise refresh) the thermal mode at every round of the
protocol. The _same_ mode can be reused in each round, in spite of the
correlation build-up and the back-reaction, with no re-thermalisation needed.
This provides a simplified and explicit version of Corollary 1 (proof in
Methods, Sec. 2.7):

######  Theorem 2.

Under the assumptions of Corollary 1, the optimal protocol in
$\mathcal{P}_{\emptyset}$ has the initial state
$\rho^{\left(0\right)}_{S}\otimes\tau_{B}$ and is such that in each round $k$:

  1. 1.

The Pauli $X\otimes\mathbb{I}$ unitary is applied to
$\rho^{\left(k-1\right)}_{SB}$, the state of system-bath after round $k-1$.

  2. 2.

The unitary $U^{\beta}_{SB}$ is applied to
$X\otimes\mathbb{I}(\rho^{\left(k-1\right)}_{SB})X\otimes\mathbb{I}$.

###  1.3 Robustness to imperfections

Theorem 2 holds even beyond some of the (standard) idealizations we made:

  1. 1.

If the bosonic mode gap does not perfectly match the system gap,
$U^{\beta}_{SB}$ still realises the cooling of Eq. (4), with the caveat that
some work flows at each $\beta$-swap.

  2. 2.

We can consider a typical imperfection, i.e. a small anharmonicity in the
ladder of $B$. In the Methods (Sec. 2.8) we show that taking
$E_{n+1}-E_{n}=E(1-(n+1)\tau^{2})+o(\tau^{2})$, the protocol of Theorem 2
performs very closely to the ideal scenario of Eq. (4) even for moderate
anharmonicity (within $0.005\%$ for $\beta E$=1, $\tau=0.05$).

Furthermore, up to now we considered an idealized scenario in which one can
perform the required $\beta$-swap exactly. In a real experiment, however, one
will only realize an approximation of it (we will see an explicit model in the
next section). Consider a noise model in which, instead of the de-excitation
probability $1$ required by the $\beta$-swap, one can only realize dephasing
thermalizations $\Lambda^{(k)}$ with an induced de-excitation probability
$\langle 0|\Lambda^{(k)}(|1\rangle\langle 1|)|0\rangle\leq 1-\epsilon$. Denote
this set by $\mathcal{P}^{\epsilon}_{\emptyset}$. Define an $\epsilon$-_noisy
$\beta$-swap_ as the dephasing thermalization $\Lambda^{\epsilon}_{\beta}$
with $\langle 0|\Lambda^{\epsilon}_{\beta}(|1\rangle\langle
1|)|0\rangle=1-\epsilon$ (when $\epsilon=0$ this is the $\beta$-swap). Then
one can derive the following noise-robust version of Theorem 1:

######  Theorem 3.

Under the assumptions of Corollary 1 and given
$\epsilon\leq\frac{1}{1+e^{\beta E}+e^{2\beta E}}$, the optimal nontrivial
cooling protocol in $\mathcal{P}^{\epsilon}_{\emptyset}$ is such that in each
round $k$:

  1. 1.

The Pauli $X$ unitary is applied to $S$.

  2. 2.

The $\epsilon$-noisy $\beta$-swap $\Lambda^{\epsilon}_{\beta}$ is applied to
$S$.

The population of the ground state after round $k$ is:

$\displaystyle p_{\rm 0}^{(k)}=$ $\displaystyle
1-\frac{\epsilon}{2-(1-\epsilon)Z}$ (6)
$\displaystyle\,\,-\left((1-\epsilon)Z-1\right)^{k}\left(1-\frac{\epsilon}{2-(1-\epsilon)Z}-p_{0}^{(0)}\right),$

where $Z=1+e^{-\beta E}$ and $p^{\left(k\right)}_{0}\rightarrow
1-\frac{\epsilon}{2-(1-\epsilon)Z}$ as $k\rightarrow\infty$.

In words: even in the presence of (moderate) noise, the optimal strategy is to
perform the best approximation to the ideal protocol; hence, the simple
structure of the optimal protocols is robust to imperfections. Theorem 3 is
proved in Section 2.9 by a tedious but straightforward optimization. There, we
also prove another form of ‘robustness’: if we optimize over the _larger_ set
of thermalization models known as _thermal operations_ [Brandão et
al.(2013)Brandão, Horodecki, Oppenheim, Renes, and Spekkens], which include
protocols where some amount of superposition in the energy basis survives the
thermalization step, Theorem 3 holds unchanged. Surprisingly, the extra
coherent control in the thermalization phase is not necessary for optimal
cooling.

###  1.4 An experimental proposal

We can now provide an explicit experimental proposal. More details about the
calculations and simulations involved can be found in Sec. 3.1 -3.2.

The unitary $U^{\beta}_{SB}$ of Eq. (5) can be realized exactly with an
intensity-dependent Jaynes-Cummings model [Naderi et al.(2005)Naderi,
Soltanolkotabi, and Roknizadeh, Åberg(2014)],
$\tilde{H}_{\textrm{JC}}=g(\sigma_{+}\otimes(aa^{\dagger})^{-1/2}a+\sigma_{-}\otimes(aa^{\dagger})^{-1/2}a^{\dagger})$.
However, a perhaps more promising avenue is to approximate the $\beta$-swap
steps with a resonant Jaynes-Cummings (JC) coupling with a thermal bosonic
mode, $H_{\textrm{JC}}=g(\sigma_{+}\otimes a+\sigma_{-}\otimes a^{\dagger})$.
Assuming good control of the interaction time $s$, we can numerically optimize
$s$ to realize an $\epsilon$-noisy $\beta$-swap with the small $\epsilon$. If
the thermal mode is reset at each round, we can compute the performance of
this implementation using Theorem 3.

In Fig. 3 we compare this protocol with leading proposals in the literature,
for an initially thermal target. The JC protocol outperforms the PPA with 2
ancillas [Schulman et al.(2005)Schulman, Mor, and Weinstein] (even with some
limitations in the timing accuracy and the total available waiting time), with
the exclusion of the very high temperature regime; when $\beta E$ is not too
small, it performs comparably (if potentially slightly worse) to the ideal
version of the non-local thermalization scheme with 1 ancilla proposed in
[Rodriguez-Briones et al.(2017)Rodriguez-Briones, Li, Peng, Mor, Weinstein,
and Laflamme]. The optimal $\beta$-swap protocol outperforms all these
protocols, but requires an intensity-dependent JC model whose implementation
we leave as an open question.

Figure 3: Ground state occupation $p^{(k)}_{0}$ at step $k$ for different
environment temperatures for: ideal optimal cooling protocol from Theorem 1
(blue), upper and lower bounds on a JC realization of this protocol (red, see
Sec. 3.1 for details on how to calculate them), SR$\Gamma$ protocol from
[Rodriguez-Briones et al.(2017)Rodriguez-Briones, Li, Peng, Mor, Weinstein,
and Laflamme] run with 1 ancilla qubit (green), PPA protocol from [Schulman et
al.(2005)Schulman, Mor, and Weinstein] run with 2 ancilla qubits (purple). The
initial state for all protocols is thermal.

These considerations assume that at every round the bosonic mode is reset to
the thermal state. However we show, using a standard master equation (Eq. (51)
in Sec. 3.2), that the reset can be substituted by a more realistic slow
rethermalization of the single mode with an external environment. Since
reasonably high cooling is achieved after 2 rounds, this suggests the
following implementation: a stream of slowly fired atoms passes through two
identical cavities resonant with the qubits we are trying to cool, each
supporting a single mode initially thermal; at the entrance of each cavity we
perform a Pauli $X$ operation, and let each qubit interact with the cavity
mode for a chosen time (see Fig. 4 for a schematic description). Numerics show
that if re-thermalization of the cavity mode is sufficiently quick compared to
the firing rate, the cooling performance settles to a constant as many atoms
are cooled (see Fig. 5).

This is, to our knowledge, the most appealing setting to experimentally
implement the protocol and is, in fact, highly reminiscent of a micromaser
[Filipowicz et al.(1986)Filipowicz, Javanainen, and Meystre, Walther et
al.(2006)Walther, Varcoe, Englert, and Becker]. This device consists of a
cavity with a harmonic oscillator in an initially thermal state, and it works
provided we are able to keep this oscillator out of thermal equilibrium. The
firing of an excited atom through a cavity in the micromaser can be seen as a
single instance of our proposed imperfect implementation of the optimal
cooling protocol. However, the figures of merit in each case are different: in
the micromaser, we need very pure atoms in order to excite the cavity
efficiently, while in the cooling protocol the aim is to obtain these very
pure atoms in the first place. Nevertheless, this suggests that experimental
settings where the micromaser has been shown to be possible might be good
platforms in which to test our algorithm. To our knowledge, this currently
includes both cavity QED [Filipowicz et al.(1986)Filipowicz, Javanainen, and
Meystre, Walther et al.(2006)Walther, Varcoe, Englert, and Becker] and solid-
state settings [Rodrigues et al.(2007)Rodrigues, Imbers, and Armour].

In conclusion, often HBAC techniques made the implicit assumption that the
best way of exploiting the external environment to pump entropy away from the
system is to thermalize the auxiliary ancillas. Here we show how optimizing
over the thermalization strategy provides new cooling protocols breaking
previously established limits, in particular allowing exponential convergence
to the ground state in the ideal scenario.

Figure 4: The optimal protocol for cooling can be approximated by one in which
a population inversion (Pauli $X$) is applied to the qubit before it enters
the cavity, where it interacts with a resonant mode initially at temperature
$\beta$. If the interaction parameters are appropriately chosen, such that a
small $\epsilon$ is achieved, the outgoing qubit has a much lower temperature.
Figure 5: Cooling achieved on each atom by two rounds of the Pauli/Jaynes-
Cummings protocols, as a function of the number of atoms already fired. Each
curve represents a different ratio between the re-thermalization rate $A$ and
the rate $r$ at which the atoms are fired. The atoms are prepared in an
initially thermal state with $\beta E=1$, the coupling strength is $g=1$ and
time of interaction $t=98.92$.

##  2 Methods

###  2.1 $\beta$-permutations and the thermal polytope

Here we discuss how to construct $\beta$-permutations and how they
characterize the thermal polytope. The first key definition is a
generalization of the concept of majorization. Recall that, given two
$d-1$-dimensional probabilities $\boldsymbol{p}$ and
$\boldsymbol{p}^{\prime}$, $\boldsymbol{p}$ majorizes
$\boldsymbol{p}^{\prime}$, denoted
$\boldsymbol{p}\succeq\boldsymbol{p}^{\prime}$, if
$\sum_{i=0}^{k}p_{i}^{\downarrow}\geq\sum_{i=0}^{k}{p^{\prime}}_{i}^{\downarrow}$
for $k=0,...,d-2$, where $\boldsymbol{x}^{\downarrow}$ denotes the vector
$\boldsymbol{x}$ arranged in descending order. $\rho\succeq\rho^{\prime}$ is
defined as the same relation among the corresponding eigenvalues. Then define

######  Definition 1 (Thermo-majorization [Horodecki and Oppenheim(2013)]).

Given a state $\rho$ with Hamiltonian $H=\sum_{i=0}^{d-1}E_{i}|i\rangle\langle
i|$, let $\boldsymbol{p}=\left(p_{0},\dots,p_{d-1}\right)$ where
$p_{i}=\langle i|\rho|i\rangle$ and
$\boldsymbol{E}=\left(E_{1},\dots,E_{d}\right)$. The thermo-majorization curve
of $\rho$ is formed by:

  1. 1.

Applying a permutation $\pi$ of $\\{0,...,d-1\\}$ to both $\boldsymbol{p}$ and
$\boldsymbol{E}$ such that $p_{\pi(i)}e^{\beta E_{\pi(i)}}$ is in non-
increasing order. We refer to the vector $\pi$ with elements $\pi(i)$ as the
_$\beta$ -order_ of $\boldsymbol{p}$.

  2. 2.

Plotting the ordered ‘elbow’ points $(0,0)$,
$\left\\{\left(\sum_{i=0}^{k}e^{-\beta
E_{\pi(i)}},\sum_{i=0}^{k}p_{\pi(i)}\right)\right\\}_{k=0}^{d-1}$ and
connecting them piecewise linearly to form a concave curve - the thermo-
majorization curve of $\boldsymbol{p}$.

Given two probability distributions $\boldsymbol{p}$ and
$\boldsymbol{p}^{\prime}$ associated with the same energy levels, we say that
$\boldsymbol{p}$ thermo-majorizes $\boldsymbol{p}^{\prime}$ if the thermo-
majorization curve of $\boldsymbol{p}$ is never below that of
$\boldsymbol{p}^{\prime}$. We denote this by
$\boldsymbol{p}\succeq_{\textrm{th}}\boldsymbol{p}^{\prime}$. Also, if $\rho$
and $\rho^{\prime}$ are states with population vectors $\boldsymbol{p}$ and
$\boldsymbol{p}^{\prime}$, the notation
$\rho\succeq_{\textrm{th}}\rho^{\prime}$ denotes
$\boldsymbol{p}\succeq_{\textrm{th}}\boldsymbol{p}^{\prime}$.

Note that $\succeq_{th}$ becomes $\succeq$ in the infinite temperature limit
$\beta\rightarrow 0$. One can show that
$\boldsymbol{p}\succeq_{\textrm{th}}\boldsymbol{p}^{\prime}$ is equivalent to
the existence of a Gibbs stochastic matrix (a stochastic matrix $G$ with
$G\boldsymbol{g}=\boldsymbol{g}$, if $g_{i}\propto e^{-\beta E_{i}}$) such
that $G\boldsymbol{p}=\boldsymbol{p}^{\prime}$ [Ruch et al.(1980)Ruch,
Schranner, and Seligman]. Furthermore, a dephasing thermalization $\Lambda$
acts on the population vector as a Gibbs-stochastic matrix $G_{j|i}=\langle
j|\Lambda(|i\rangle\langle i|)|j\rangle$; conversely, every Gibbs-stochastic
matrix can be realized by dephasing thermalizations. As such, the thermal
polytope of $\boldsymbol{p}$ coincides with the set of probability
distributions $\boldsymbol{p}^{\prime}$ such that
$\boldsymbol{p}\succeq_{\textrm{th}}\boldsymbol{p}^{\prime}$. That this is a
convex polytope follows from the fact that the set of Gibbs-stochastic
matrices also is. To characterize the thermal polytope of $\boldsymbol{p}$
explicitly, we construct a set of maps, called $\beta$-permutations, that take
$\boldsymbol{p}$ to each one of its extremal points.

_Algorithm to construct $\Lambda^{(\pi,\alpha)}$._ By the discussion above
each $\beta$-permutation $\Lambda^{(\pi,\alpha)}$ is fully characterized by a
matrix of transition probabilities $P^{(\pi,\alpha)}_{j|i}=\langle
j|\Lambda^{(\pi,\alpha)}(|i\rangle\langle i|)|j\rangle$. Let $\pi$ and
$\alpha$ each be one of the $d!-1$ possible $\beta$-orders. Order the rows of
$P^{(\pi,\alpha)}$ according to $\alpha$ and the columns according to $\pi$:

$G=\begin{pmatrix}P^{(\pi,\alpha)}_{\alpha(0)\pi(0)}&\dots&P^{(\pi,\alpha)}_{\alpha(0)\pi(d-1)}\\\
\vdots&&\vdots\\\
P^{(\pi,\alpha)}_{\alpha(d-1)\pi(0)}&\dots&P^{(\pi,\alpha)}_{\alpha(d-1)\pi(d-1)}\end{pmatrix}.$

Then $G$ is constructed as follows.

_Row 0_ : If $e^{-\beta E_{\alpha\left(0\right)}}<e^{-\beta
E_{\pi\left(0\right)}}$, set:

$\displaystyle k_{0}=1,\;\;G_{00}=\frac{e^{-\beta
E_{\alpha\left(0\right)}}}{e^{-\beta
E_{\pi\left(0\right)}}},\;\;G_{0j}=0,\quad j\in[1,d-1],$

else, let $k_{0}$ be the smallest integer such that:

$\sum_{j=0}^{k_{0}}e^{-\beta E_{\pi\left(j\right)}}\geq e^{-\beta
E_{\alpha\left(0\right)}}.$

and set:

$\displaystyle G_{0j}$ $\displaystyle=1,\quad 0\leq j\leq k_{0}-1$
$\displaystyle G_{0k_{0}}$ $\displaystyle=\frac{e^{-\beta
E_{\alpha\left(0\right)}}-\sum_{j=0}^{k_{0}-1}e^{-\beta
E_{\pi\left(j\right)}}}{e^{-\beta E_{\pi\left(k_{0}\right)}}}$ $\displaystyle
G_{0j}$ $\displaystyle=0,\quad j\in[k_{0}+1,d-1]$

_Row $m\geq 1$_: If $\sum_{i=0}^{m}e^{-\beta
E_{\alpha\left(i\right)}}<\sum_{j=0}^{k_{m-1}}e^{-\beta
E_{\pi\left(j\right)}}$, set:

$\displaystyle k_{m}$ $\displaystyle=k_{m-1},\;\;G_{mj}=0,\quad 0\leq j\leq
k_{m-1}-1$ $\displaystyle G_{mk_{m-1}}$ $\displaystyle=\frac{e^{-\beta
E_{\alpha\left(m\right)}}}{e^{-\beta
E_{\pi\left(k_{m-1}\right)}}},\;G_{mj}=0,\;j\in[k_{m-1}+1,d-1]$

else, let $k_{m}$ be the smallest integer such that:

$\sum_{j=0}^{k_{m}}e^{-\beta E_{\pi\left(j\right)}}\geq\sum_{i=0}^{m}e^{-\beta
E_{\alpha\left(i\right)}}.$ (7)

and set:

$\displaystyle G_{mj}=0,\;j$
$\displaystyle\in[0,k_{m-1}-1],\;G_{mk_{m-1}}=1-\sum_{i=0}^{m-1}G_{ik_{m-1}}$
$\displaystyle G_{mj}$ $\displaystyle=1,\quad j\in[k_{m-1}+1,k_{m}-1]$
$\displaystyle G_{mk_{m}}$ $\displaystyle=\frac{\sum_{i=0}^{m}e^{-\beta
E_{\alpha\left(i\right)}}-\sum_{j=0}^{k_{m}-1}e^{-\beta
E_{\pi\left(j\right)}}}{e^{-\beta E_{\pi\left(k_{m}\right)}}}$ $\displaystyle
G_{mj}$ $\displaystyle=0,\quad j\in[k_{m}+1,d-1].$

In the next section, we prove the important fact that, if $\pi$ is the
$\beta$-order of $\boldsymbol{p}$, the set
$\\{P^{(\pi,\alpha)}\boldsymbol{p}\\}_{\alpha}$, for $\alpha$ varying over all
possible $\beta$-orders, includes all extremal points of the thermal polytope
of $\boldsymbol{p}$ (Lemma 2).

###  2.2 Proof of properties of $\beta$-permutations

In the following we show that the $\beta$-permutations, defined Section 1.1,
convert $\boldsymbol{p}$ into any of the extremal points of the thermal
polytope $\boldsymbol{p}^{\alpha}$, while being independent of the specific
form of $\boldsymbol{p}$ (the set of $\beta$-permutation that need to be
applied only depends on the $\beta$-order of $\boldsymbol{p}$). This is in
complete analogy with permutations and their role within the theory of doubly-
stochastic matrices as determined by Birkhoff’s theorem [Birkhoff(1946)] (in
fact, $\beta$-permutations become permutations in the infinite temperature
limit).

Specifically, we will prove the following facts:

  1. 1.

$P^{(\pi,\alpha)}$ is a Gibbs-stochastic map.

  2. 2.

Given $\boldsymbol{p}$ with $\beta$-order $\pi$,
$\boldsymbol{p}^{\alpha}:=P^{(\pi,\alpha)}\boldsymbol{p}$ has $\beta$-order
$\alpha$ and is such that there is no $\boldsymbol{q}$ with $\beta$-order
$\alpha$ such that
$\boldsymbol{p}\succeq_{th}\boldsymbol{q}\succeq_{th}\boldsymbol{p}^{\alpha}$
(Lemma 2).

  3. 3.

The set of $\boldsymbol{p}^{\alpha}$ for varying $\alpha$ includes all the
extremal points of the themal polytope of $\boldsymbol{p}$ (the latter is
equivalently defined as all $\boldsymbol{q}$ such that
$\boldsymbol{p}\succeq_{th}\boldsymbol{q}$).

####  2.2.1 Proof of claim 1

Recall that we reordered the rows of $P^{(\pi,\alpha)}$ according to $\alpha$
and the columns according to $\pi$:

$G=\begin{pmatrix}P^{(\pi,\alpha)}_{\alpha(0)\pi(0)}&\dots&P^{(\pi,\alpha)}_{\alpha(0)\pi(d-1)}\\\
\vdots&&\vdots\\\
P^{(\pi,\alpha)}_{\alpha(d-1)\pi(0)}&\dots&P^{(\pi,\alpha)}_{\alpha(d-1)\pi(d-1)}\end{pmatrix}.$

Hence, the condition of Gibbs-stochasticity can be rewritten as

$\displaystyle G_{ij}$ $\displaystyle\geq 0$ (8)
$\displaystyle\sum_{i=0}^{d-1}G_{ij}$ $\displaystyle=1$ (9)
$\displaystyle\sum_{j=0}^{d-1}G_{ij}e^{-\beta E_{\pi\left(j\right)}}$
$\displaystyle=e^{-\beta E_{\alpha\left(i\right)}}.$ (10)

For simplicity, we repeat here the algorithm for constructing
$P^{(\pi,\alpha)}$. Row 0 is populated as follows:

if $e^{-\beta E_{\alpha\left(0\right)}}<e^{-\beta E_{\pi\left(0\right)}}$ then

Set:  $\displaystyle k_{0}$ $\displaystyle=1$ (11) $\displaystyle G_{00}$
$\displaystyle=\frac{e^{-\beta E_{\alpha\left(0\right)}}}{e^{-\beta
E_{\pi\left(0\right)}}}$ (12) $\displaystyle G_{0j}$ $\displaystyle=0,\quad
1\leq j\leq d-1$ (13)

else

Let $k_{0}$ be the smallest integer such that:  $\sum_{j=0}^{k_{0}}e^{-\beta
E_{\pi\left(j\right)}}\geq e^{-\beta E_{\alpha\left(0\right)}}.$ (14)

Set:  $\displaystyle G_{0j}$ $\displaystyle=1,\quad 0\leq j\leq k_{0}-1$ (15)
$\displaystyle G_{0k_{0}}$ $\displaystyle=\frac{e^{-\beta
E_{\alpha\left(0\right)}}-\sum_{j=0}^{k_{0}-1}e^{-\beta
E_{\pi\left(j\right)}}}{e^{-\beta E_{\pi\left(k_{0}\right)}}}$ (16)
$\displaystyle G_{0j}$ $\displaystyle=0,\quad k_{0}+1\leq j\leq d-1$ (17)

end if

Figure 6: Here we illustrate the notion of a _maximal $\beta$-order state_.
The state $\boldsymbol{p}$, shown in blue, has $\beta$-order
$\pi=\left(0,1,2\right)$. The state $\boldsymbol{p}^{\alpha}$ for
$\alpha=\left(2,0,1\right)$ is shown in red and has a simple geometrical
interpretation. The choice of $\alpha$ (Condition 1 in Sec. 2.2.2) fixes the
$x$-axis points to be $x_{0}=e^{-\beta E_{2}}$, $x_{1}=e^{-\beta
E_{2}}+e^{-\beta E_{0}}$, $x_{2}=e^{-\beta E_{2}}+e^{-\beta E_{0}}+e^{-\beta
E_{1}}$. Then Condition 2 is equivalent to the request that the curve in red
touches the blue curve at these points.

Let us first check that Eqs. (8) and (10) are fulfilled in each clause of row
$m=0$. In the first clause, it is clear that $G_{00}\geq 0$. Gibbs
preservation follows as:

$\sum_{j=0}^{d-1}G_{0j}e^{-\beta E_{\pi\left(j\right)}}=G_{00}e^{-\beta
E_{\pi\left(0\right)}}=e^{-\beta E_{\alpha\left(0\right)}}.$

For the second clause, it is again true that we have $0\leq G_{0j}\leq 1$, for
all $j$ (if $G_{0k_{0}}$ were greater than $1$, then this would contradict the
definition of $k_{0}$ in Eq. (14)). We then note that:

$\displaystyle\sum_{j=0}^{d-1}G_{0j}e^{-\beta
E_{\pi\left(j\right)}}=\sum_{j=0}^{k_{0}-1}e^{-\beta
E_{\pi\left(j\right)}}+e^{-\beta
E_{\alpha\left(0\right)}}-\sum_{j=0}^{k_{0}-1}e^{-\beta
E_{\pi\left(j\right)}}=e^{-\beta E_{\alpha\left(0\right)}}$

so Gibbs preservation is satisfied.

With row zero in place, row $m$ ($m\in\left\\{1,\dots d-1\right\\}$) is
populated as follows:

if $\sum_{i=0}^{m}e^{-\beta
E_{\alpha\left(i\right)}}<\sum_{j=0}^{k_{m-1}}e^{-\beta
E_{\pi\left(j\right)}}$ then

Set:  $\displaystyle k_{m}$ $\displaystyle=k_{m-1}$ (18) $\displaystyle
G_{mj}$ $\displaystyle=0,\quad 0\leq j\leq k_{m-1}-1$ (19) $\displaystyle
G_{mk_{m-1}}$ $\displaystyle=\frac{e^{-\beta
E_{\alpha\left(m\right)}}}{e^{-\beta E_{\pi\left(k_{m-1}\right)}}}$ (20)
$\displaystyle G_{mj}$ $\displaystyle=0,\quad k_{m-1}+1\leq j\leq d-1$ (21)

else

Let $k_{m}$ be the smallest integer such that:  $\sum_{j=0}^{k_{m}}e^{-\beta
E_{\pi\left(j\right)}}\geq\sum_{i=0}^{m}e^{-\beta E_{\alpha\left(i\right)}}$
(22)

Set:  $\displaystyle G_{mj}$ $\displaystyle=0,\quad 0\leq j\leq k_{m-1}-1$
(23) $\displaystyle G_{mk_{m-1}}$
$\displaystyle=1-\sum_{i=0}^{m-1}G_{ik_{m-1}}$ (24) $\displaystyle G_{mj}$
$\displaystyle=1,\quad k_{m-1}+1\leq j\leq k_{m}-1$ (25) $\displaystyle
G_{mk_{m}}$ $\displaystyle=\frac{\sum_{i=0}^{m}e^{-\beta
E_{\alpha\left(i\right)}}-\sum_{j=0}^{k_{m}-1}e^{-\beta
E_{\pi\left(j\right)}}}{e^{-\beta E_{\pi\left(k_{m}\right)}}}$ (26)
$\displaystyle G_{mj}$ $\displaystyle=0,\quad k_{m}+1\leq j\leq d-1$ (27)

end if

We now check that Eq. (8) and Eq. (10) are fulfilled in each of the clauses
for each row in the matrix. For the first clause $0\leq G_{mk_{m-1}}\leq 1$
follows from:

$\displaystyle e^{-\beta E_{\alpha(m)}}$
$\displaystyle\leq\sum_{i=0}^{k_{m-1}}e^{-\beta
E_{\pi(i)}}-\sum_{i=0}^{m-1}e^{-\beta E_{\alpha(i)}}$ $\displaystyle=e^{-\beta
E_{\pi(k_{m-1})}}+\sum_{i=0}^{k_{m-1}-1}e^{-\beta
E_{\pi(i)}}-\sum_{i=0}^{m-1}e^{-\beta E_{\alpha(i)}}$ $\displaystyle\leq
e^{-\beta E_{\pi(k_{m-1})}},$

where in the last inequality we used the definition of $k_{m-1}$. Finally,
note that:

$\sum_{j=0}^{d-1}G_{mj}e^{-\beta E_{\pi\left(j\right)}}=G_{mk_{m-1}}e^{-\beta
E_{\pi\left(k_{m-1}\right)}}=e^{-\beta E_{\alpha\left(m\right)}}$

and hence Eq. (10) holds.

For the second clause, the definitions of $k_{m-1}$ and $k_{m}$ ensure that
both $0\leq G_{mk_{m-1}}\leq 1$ and $0\leq G_{mk_{m}}\leq 1$ hold. We now show
Gibbs preservation. First note that:

$\displaystyle\sum_{j=0}^{d-1}G_{mj}e^{-\beta E_{\pi\left(j\right)}}$
$\displaystyle=G_{mk_{m-1}}e^{-\beta
E_{\pi\left(k_{m-1}\right)}}+\sum_{j=k_{m-1}+1}^{k_{m}-1}e^{-\beta
E_{\pi\left(j\right)}}+G_{mk_{m}}e^{-\beta E_{\pi\left(k_{m}\right)}}$
$\displaystyle=\left(1-\sum_{i=0}^{m-1}G_{ik_{m-1}}\right)e^{-\beta
E_{\pi\left(k_{m-1}\right)}}+\sum_{j=k_{m-1}+1}^{k_{m}-1}e^{-\beta
E_{\pi\left(j\right)}}+\sum_{i=0}^{m}e^{-\beta
E_{\alpha\left(i\right)}}-\sum_{j=0}^{k_{m}-1}e^{-\beta
E_{\pi\left(j\right)}}$ (28)

where we have used Eq. (24) and Eq. (26) in the last line. Consider now
$\sum_{i=0}^{m-1}G_{ik_{m-1}}$. Note that there exists an integer $r$ such
that $0\leq r\leq m-1$ and:

$\displaystyle G_{ik_{m-1}}$ $\displaystyle=0,\quad i<r$ $\displaystyle
G_{rk_{m-1}}$ $\displaystyle=\frac{\sum_{i=0}^{r}e^{-\beta
E_{\alpha\left(i\right)}}-\sum_{j=0}^{k_{r}-1}e^{-\beta
E_{\pi\left(j\right)}}}{e^{-\beta E_{\pi\left(k_{r}\right)}}}$ $\displaystyle
G_{ik_{m-1}}$ $\displaystyle=\frac{e^{-\beta
E_{\alpha\left(i\right)}}}{e^{-\beta E_{\pi\left(k_{i}\right)}}},\quad r+1\leq
i\leq m-1$

with $k_{i}=k_{r}$ for $r\leq i\leq m-1$. Using these expressions in Eq. (28),
we see that:

$\displaystyle\sum_{j=0}^{d-1}G_{mj}e^{-\beta E_{\pi\left(j\right)}}$
$\displaystyle=\left(1-\frac{\sum_{i=0}^{m-1}e^{-\beta
E_{\alpha\left(i\right)}}-\sum_{j=0}^{k_{m-1}-1}e^{-\beta
E_{\pi\left(j\right)}}}{e^{-\beta
E_{\pi\left(k_{m-1}\right)}}}\right)e^{-\beta E_{\pi\left(k_{m-1}\right)}}$
$\displaystyle\quad\quad+\sum_{j=k_{m-1}+1}^{k_{m}-1}e^{-\beta
E_{\pi\left(j\right)}}+\sum_{i=0}^{m}e^{-\beta
E_{\alpha\left(i\right)}}-\sum_{j=0}^{k_{m}-1}e^{-\beta
E_{\pi\left(j\right)}}$ $\displaystyle=\sum_{j=0}^{k_{m}-1}e^{-\beta
E_{\pi\left(j\right)}}+e^{-\beta
E_{\alpha\left(m\right)}}-\sum_{j=0}^{k_{m}-1}e^{-\beta
E_{\pi\left(j\right)}}$ $\displaystyle=e^{-\beta E_{\alpha\left(m\right)}}$

and hence the Gibbs distribution is preserved.

Finally, let us show that the matrix constructed is stochastic. Using the
expressions for $G_{ik_{m-1}}$ given above, the fact that $k_{i}=k_{r}$ for
all $r\leq i\leq m-1$ and the definition of $k_{m-1}$,

$\sum_{i=0}^{m-1}G_{ik_{m-1}}=\frac{\sum_{i=0}^{m-1}e^{-\beta
E_{\alpha(i)}}-\sum_{i=0}^{k_{m-1}-1}e^{-\beta E_{\pi(i)}}}{e^{-\beta
E_{\pi(k_{m-1})}}}<1.$

Together with Eq. (24) and the already proved fact that $G_{ij}\in[0,1]$ for
all $i,j\in\left\\{0,\dots,d-1\right\\}$, we get that $G$ is stochastic.

####  2.2.2 Proof of claim 2

This can be stated as the following Lemma:

######  Lemma 2.

If $\boldsymbol{p}$ has $\beta$-order $\pi$ and
$\boldsymbol{p}^{\alpha}:=P^{(\pi,\alpha)}\boldsymbol{p}$, then
$\boldsymbol{p}^{\alpha}$ has $\beta$-order $\alpha$ and is ‘maximal’ in the
sense that there is no $\boldsymbol{q}$ with $\beta$-order $\alpha$ such that
$\boldsymbol{p}\succeq_{\textrm{th}}\boldsymbol{q}\succeq_{\textrm{th}}\boldsymbol{p}^{\alpha}.$

We now construct the proof. $\beta$-permutations have a simple geometrical
description in terms of thermo-majorization curves. In row zero, we want to
maximize the population of $E_{\alpha\left(0\right)}$ subject to the thermo-
majorization constraints. To do this, we compare $e^{-\beta
E_{\alpha\left(0\right)}}$ and $e^{-\beta E_{\pi\left(0\right)}}$. If
$e^{-\beta E_{\alpha\left(0\right)}}<e^{-\beta E_{\pi\left(0\right)}}$, then
we cannot move all of the population in $E_{\pi\left(0\right)}$ to
$E_{\alpha\left(0\right)}$ without violating thermo-majorization and must
instead move only a fraction of it, as given by Eq. (12). This case is
illustrated in Fig. 7a. On the other hand, if $e^{-\beta
E_{\alpha\left(0\right)}}\geq e^{-\beta E_{\pi\left(0\right)}}$, then we can
move all of the population in $E_{\pi\left(0\right)}$ to
$E_{\alpha\left(0\right)}$ and set $G_{00}=1$. We then try to move population
from $E_{\pi\left(1\right)}$ to $E_{\alpha\left(0\right)}$ and repeat this
process until we reach an energy level whose population we cannot move into
$E_{\alpha\left(0\right)}$ without violating thermo-majorization. This energy
level is defined through Eq. (14) and given by $E_{\pi\left(k_{0}\right)}$.
From $E_{\pi\left(k_{0}\right)}$ we can only move a fraction of the population
into $E_{\alpha\left(0\right)}$, given by Eq. (16). When we reach this point,
we have moved as much population as possible into $E_{\alpha\left(0\right)}$.
This case is illustrated in Fig. 7b.

In the $m$th row we are determining the population of
$E_{\alpha\left(m\right)}$ after the transformation. At this point in the
construction, all of the population from energy level $E_{\pi\left(j\right)}$
for $0\leq j\leq k_{m-1}-1$ has been transferred already into the set of
energies $\left\\{E_{\alpha\left(i\right)}\right\\}_{i=0}^{m-1}$. We thus have
$G_{mj}=0$ for $0\leq j\leq k_{m-1}-1$. We now check to see how much of the
remaining population in $E_{\pi\left(k_{m-1}\right)}$ we can move to
$E_{\alpha\left(m\right)}$ subject to the thermo-majorization constraints. If
we can only move a fraction of it, we follow the ‘if’ clause in the above
algorithm and determine $G_{mk_{m-1}}$ to be given by Eq. (20). This is
illustrated in Fig. 8a. Alternatively, if we can move all of it, we follow the
‘else’ clause and $G_{mk_{m-1}}$ is given by Eq. (24). The rest of the
construction follows a similar line of argument to the first row. We move all
of the population from the energy levels after $E_{\pi\left(k_{m-1}\right)}$
to $E_{\alpha\left(m\right)}$ until we reach an energy level
$E_{\pi\left(k_{m}\right)}$, where this is not possible due to thermo-
majorization. This is illustrated in Fig. 8b.

Let $c_{\boldsymbol{p}}:\left[0,Z_{S}\right]\rightarrow\left[0,1\right]$ be
the function such that $c_{\boldsymbol{p}}(x)$ is the height of the thermo-
majorisation curve of $\boldsymbol{p}$ at $x$. The action of the
$\beta$-permutation matrix described above is such that the thermo-
majorization curve of $\boldsymbol{p}^{\alpha}=P^{(\pi,\alpha)}\boldsymbol{p}$
is constructed as follows. For $i\in\left\\{0,\dots,d-1\right\\}$, denoting by
$(x^{\alpha}_{i},y^{\alpha}_{i})$ the points that are piecewise linearly
connected to give the thermo-majorization curve of $\boldsymbol{p}^{\alpha}$,
one has:

  1. 1.

Let $x_{i}^{\alpha}=\sum_{j=0}^{i}e^{-\beta E_{\alpha^{-1}\left(j\right)}}$
and $y_{i}^{\alpha}=c_{\boldsymbol{p}}\left(x_{i}^{\alpha}\right)$.

  2. 2.

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

By construction, there is no $\boldsymbol{q}$ with $\beta$-order $\alpha$ such
that
$\boldsymbol{p}\succeq_{th}\boldsymbol{q}\succeq_{th}\boldsymbol{p}^{\alpha}$.

####  2.2.3 Proof of claim 3

That all extremal points of the thermal polytope associated to
$\boldsymbol{p}$ have the above form is stated and proved in Lemma 12 of Ref.
[Lostaglio et al.(2018)Lostaglio, Alhambra, and Perry].

(a) Row 0, case 1

(b) Row 0, case 2

Figure 7: Thermo-majorization diagrams illustrating the construction of row
zero of $P^{(\pi,\alpha)}$. In case 1, $e^{-\beta
E_{\alpha\left(0\right)}}<e^{-\beta E_{\pi\left(0\right)}}$ holds, while in
case 2 it does not.

(a) Row $m$, case 1

(b) Row $m$, case 2

Figure 8: Thermo-majorization diagrams illustrating the construction of the
row $m$ of $P^{(\pi,\alpha)}$. In case 1, $\sum_{i=0}^{m}e^{-\beta
E_{\alpha\left(i\right)}}<\sum_{j=0}^{k_{m-1}}e^{-\beta
E_{\pi\left(j\right)}}$ holds, while in case 2 it does not.

###  2.3 Proof of Lemma 1: optimal unitary control

First we recall the notion of _maximally active state_ associated to a given
$\rho$, defined as:

######  Definition 2 (Maximally active state).

Let $\rho$ be the state of a system with associated Hamiltonian
$H=\sum_{i=0}^{d-1}E_{i}|i\rangle\langle i|$, where $E_{i}\leq E_{i+1}$. The
maximally active state associated to $\rho$ is

$\hat{\rho}=\sum_{i=0}^{d-1}\lambda^{\uparrow}_{i}|i\rangle\langle i|,$

where $\left\\{\lambda^{\uparrow}_{i}\right\\}_{i=0}^{d-1}$ are the
eigenvalues of $\rho$ arranged in ascending order.

This corresponds to the state with the highest energy along the whole unitary
orbit of $\rho$, i.e. $\hat{\rho}=\arg\max_{U}\textnormal{Tr}[U\rho
U^{\dagger}H]$. We can now prove Lemma 1 (note that we use facts from Section
2.1 in this proof):

###### Proof.

The claim is equivalent to proving $\hat{\rho}\succeq_{\textrm{th}}U\rho
U^{\dagger}$ for every unitary $U$. We will first show that for any
permutation $\pi$ of $\\{1,...,d\\}$:

$\boldsymbol{\lambda}^{\uparrow}\succeq_{\textrm{th}}\pi\boldsymbol{\lambda}^{\uparrow},$
(29)

where $\boldsymbol{\lambda}^{\uparrow}$ denotes the vector of eigenvalues of
$\hat{\rho}$ arranged in ascending order. To show it, we will construct a
sequence of $d-1$ Gibbs-stochastic matrices converting
$\boldsymbol{\lambda}^{\uparrow}$ into $\pi\boldsymbol{\lambda}^{\uparrow}$,
passing through the intermediate states
$\boldsymbol{q}^{\left(1\right)},\dots,\boldsymbol{q}^{\left(d-1\right)}$,
with $\boldsymbol{q}^{\left(d-1\right)}=\pi\boldsymbol{\lambda}^{\uparrow}$
and $\boldsymbol{q}^{\left(0\right)}=\boldsymbol{\lambda}^{\uparrow}$.

We start with $\boldsymbol{q}^{\left(0\right)}$. We wish to move the
probability $\lambda^{\uparrow}_{0}$ associated with energy level $E_{0}$ to
energy level $E_{\pi\left(0\right)}$. To do this, we transpose in sequence the
populations of energy levels $0\leftrightarrow 1$, $1\leftrightarrow 2$,
$\dots$, $\pi(0)-1\leftrightarrow\pi(0)$. Each of these is possible under
Gibbs-stochastic matrices, because
$\lambda^{\uparrow}_{0}\leq\lambda^{\uparrow}_{j+1}$, $\forall j$ and
$E_{j}\leq E_{j+1}$. In fact, they are achieved with the matrices:

$G_{\left(j,j+1\right)}=\begin{pmatrix}1-\mu
e^{-\beta\left(E_{j+1}-E_{j}\right)}&\mu\\\ \mu
e^{-\beta\left(E_{j+1}-E_{j}\right)}&1-\mu\end{pmatrix}\oplus\mathbb{I}_{\backslash\left(j,j+1\right)},$

where $\mathbb{I}_{\backslash\left(j,j+1\right)}$ denotes the identity on all
levels different from $(j,j+1)$ and

$\mu=\frac{\lambda^{\uparrow}_{j+1}-\lambda^{\uparrow}_{0}}{\lambda^{\uparrow}_{j+1}-\lambda^{\uparrow}_{0}e^{-\beta\left(E_{j+1}-E_{j}\right)}}\geq
0.$

Hence it is always possible to move population $\lambda^{\uparrow}_{0}$ to
energy level $E_{\pi\left(0\right)}$ by setting

$\boldsymbol{q}^{\left(1\right)}=G_{(\pi(0)-1,\pi(0))}\dots
G_{(1,2)}G_{(0,1)}\boldsymbol{q}^{(0)}.$

$\boldsymbol{q}^{\left(1\right)}$ coincides with
$\pi\boldsymbol{\lambda}^{\uparrow}$ on element $\pi(0)$. Now, if we truncate
element $\pi(0)$ from $\boldsymbol{q}^{\left(1\right)}$, we obtain a vector
satisfying

$\displaystyle q^{\left(1\right)}_{0}=\lambda^{\uparrow}_{1}\leq
q^{\left(1\right)}_{1}=\lambda^{\uparrow}_{2}\leq\dots$ $\displaystyle\leq
q^{\left(1\right)}_{\pi(0)-1}\leq q^{\left(1\right)}_{\pi(0)+1}\leq\dots\leq
q^{\left(1\right)}_{d-1}=\lambda^{\uparrow}_{d-1}.$

Reasoning as before, we find a second sequence of Gibbs-stochastic matrices
acting on the truncated vector, that transposes the populations in adjacent
positions and moves $\lambda^{\uparrow}_{1}$ to the energy level
$E_{\pi\left(1\right)}$. Applied to $\boldsymbol{q}^{(1)}$, this sequence give
a state $\boldsymbol{q}^{\left(2\right)}$ which coincides with
$\pi\boldsymbol{\lambda}^{\uparrow}$ in elements $\pi(0),\pi(1)$. It should be
clear that this construction can be repeated on every intermediate
$\boldsymbol{q}^{(m)}$, each time applying it to the distribution in which we
ignore the energy levels $E_{\pi(0)},...,E_{\pi(m-1)}$, since these have the
correct occupation probability. This provides a sequence of states culminating
in $\boldsymbol{q}^{\left(d-1\right)}=\pi\boldsymbol{\lambda}^{\uparrow}$ as
required.

Having shown that Eq. (29) holds, we now argue for
$\hat{\rho}\succeq_{\textrm{th}}U\rho U^{\dagger}$ for all unitaries $U$.
Given $U$, let $\tilde{\boldsymbol{p}}=\textrm{diag}\left[U\rho
U^{\dagger}\right]$ (remembering to first diagonalize the degenerate energy
subspaces if necessary) and we want to prove
$\boldsymbol{\lambda}^{\uparrow}\succeq_{th}\tilde{\boldsymbol{p}}$. Let
$\tilde{\boldsymbol{p}}^{\uparrow}$ be the state formed by arranging the
elements of $\tilde{\boldsymbol{p}}$ in ascending order. By Eq. (29), we have
that
$\tilde{\boldsymbol{p}}^{\uparrow}\succeq_{\textrm{th}}\tilde{\boldsymbol{p}}$.
In addition, by the Schur-Horn theorem, we have that
$\boldsymbol{\lambda}^{\uparrow}\succeq\tilde{\boldsymbol{p}}^{\uparrow}$.
Combining this with the fact that both $\boldsymbol{\lambda}^{\uparrow}$ and
$\tilde{\boldsymbol{p}}^{\uparrow}$ are ordered in terms of increasing
occupation probability and have the same $\beta$-order, it follows that:

  1. 1.

If the thermo-majorization curve of $\boldsymbol{\lambda}^{\uparrow}$ is
constructed from points $\left\\{(x_{m},y_{m})\right\\}_{m=0}^{d-1}$ and that
of $\tilde{\boldsymbol{p}}^{\uparrow}$ from points
$\left\\{(x^{\prime}_{m},y^{\prime}_{m})\right\\}_{m=0}^{d-1}$, we have
$x_{m}=x^{\prime}_{m}$ for each $m$.

  2. 2.

$y_{m}\geq y^{\prime}_{m}$ for each $m$, since
$\sum_{i=0}^{m}\lambda^{\uparrow}_{i}\geq\sum_{i=0}^{m}\tilde{p}^{\uparrow}_{i}$
for all $m=1,\dots,d-1$.

Hence
$\boldsymbol{\lambda}^{\uparrow}\succeq_{th}\tilde{\boldsymbol{p}}^{\uparrow}\succeq_{th}\tilde{\boldsymbol{p}}$,
which implies the claim. ∎

###  2.4 Proof of Theorem 1

First, let us focus on optimality in a single round. That the initial unitary
can be taken to be the one mapping
$\rho^{\left(k-1\right)}_{S}\otimes\rho_{A}$ to the corresponding maximally
active state follows immediately from Lemma 1. Next, we perform the dephasing
thermalization that maximizes the population of the ground state of $S$ within
the thermal polytope. If $\pi_{k}$ is the $\beta$-order of the state after the
unitary, this can be chosen to be the $\beta$-permutation
$\Lambda^{(\pi_{k},\alpha)}$ with $\alpha$ given in Eq. (1). That this is an
ordering that maximizes the ground state population of $S$ follows from the
concavity of thermo-majorization curves. Furthermore, it maximizes
$\sum_{i=0}^{l}p^{\downarrow\left(k\right)}_{i}$,
$l\in\left\\{0,\dots,d-1\right\\}$ where
$\left\\{p^{\downarrow\left(k\right)}_{i}\right\\}_{i=0}^{d-1}$ are the
populations of $\rho^{\left(k\right)}_{S}$ arranged in descending order. This
follows from Lemma 2 and the fact that the $\beta$-permutation of Eq. (1) is
the one that maximises the $x$-axis coordinates of the elbow points of the
output thermo-majorization curve.

We now formally show that the concatenation of such rounds forms an optimal
protocol in $\mathcal{P}_{\rho_{A}}$. Suppose that at the beginning of round
$k$ we have one of two diagonal states $\rho^{\left(k-1\right)}_{S}$ and
$\tilde{\rho}^{\left(k-1\right)}_{S}$ such that
$\rho^{\left(k-1\right)}_{S}\succeq\tilde{\rho}^{\left(k-1\right)}_{S}$.
$\rho^{(k)}_{S}$ will represent the trajectory followed by the state when we
apply the claimed optimal protocol, whereas
$\tilde{\rho}^{{\left(k\right)}}_{S}$ will be the trajectory followed by a
generic protocol (hence, $\rho^{(0)}_{S}=\tilde{\rho}^{(0)}_{S}$). Our goal is
to show that:

$\displaystyle\rho^{\left(k\right)}_{S}$
$\displaystyle=\textnormal{Tr}_{A}\left[\Lambda^{(\pi_{k},\alpha)}\circ\mathcal{U}_{\textrm{m.a.}}^{\left(k\right)}\left(\rho^{\left(k-1\right)}_{S}\otimes\rho_{A}\right)\right]\succeq\tilde{\rho}^{\left(k\right)}_{S}=\textnormal{Tr}_{A}\left[\Lambda^{\left(k\right)}\circ\mathcal{V}^{\left(k\right)}\left(\tilde{\rho}^{\left(k-1\right)}_{S}\otimes\rho_{A}\right)\right]$
(30)

for all choices of $\rho_{A}$, $\mathcal{V}^{(k)}$ and $\Lambda^{(k)}$. Here
$\mathcal{U}_{\textrm{m.a.}}^{\left(k\right)}$ denotes the unitary creating
the maximally active state on $SA$, $\mathcal{V}^{\left(k\right)}$ is an
arbitrary unitary and $\Lambda^{\left(k\right)}$ a dephasing thermalization.
This implies not only that our protocol maximizes the ground state population
in a given round but also that deviating from it can only have an adverse
effect on the achievable population in subsequent rounds.

To see that Eq. (30) holds, first note that as
$\rho^{\left(k-1\right)}_{S}\succeq\tilde{\rho}^{\left(k-1\right)}_{S}$, then:

$\rho^{\left(k-1\right)}_{S}\otimes\rho_{A}\succeq\tilde{\rho}^{\left(k-1\right)}_{S}\otimes\rho_{A},\quad\forall\rho_{A}.$
(31)

This gives:

$\displaystyle\mathcal{U}_{\textrm{m.a.}}^{\left(k\right)}\left(\rho^{\left(k-1\right)}_{S}\otimes\sigma_{A}\right)\succeq_{\textrm{th}}\mathcal{U}_{\textrm{m.a.}}^{\left(k\right)}\left(\tilde{\rho}^{\left(k-1\right)}_{S}\otimes\sigma_{A}\right)\succeq_{\textrm{th}}\mathcal{V}^{\left(k\right)}\left(\tilde{\rho}^{\left(k-1\right)}_{S}\otimes\sigma_{A}\right),\quad\forall\mathcal{V}^{\left(k\right)}.$
(32)

The first line follows from Eq. (31) by comparison of the thermomajorization
curves, since maximally active states have the same $\beta$-orders. The second
line follows from Lemma 1. Now, by definition of $\Lambda^{(\pi_{k},\alpha)}$
and Eq. (32),

$\displaystyle\textnormal{Tr}_{A}$
$\displaystyle\left[\Lambda^{(\pi_{k},\alpha)}\circ\mathcal{U}_{\textrm{m.a.}}^{\left(k\right)}\left(\rho^{\left(k-1\right)}_{S}\otimes\sigma_{A}\right)\right]\succeq\textnormal{Tr}_{A}\left[\Lambda^{\left(k\right)}\circ\mathcal{V}^{\left(k\right)}\left(\tilde{\rho}^{\left(k-1\right)}_{S}\otimes\sigma_{A}\right)\right].$

###  2.5 A qudit protocol with $\beta$-swaps and no ancillas

The $\beta$-permutations affecting only two energy levels $(i,j)$ at once are
called _$\beta$ -swaps_ [Lostaglio et al.(2018)Lostaglio, Alhambra, and Perry]
and have the form

$\beta_{i,j}=\begin{pmatrix}1-e^{-\beta\left(E_{j}-E_{i}\right)}&1\\\
e^{-\beta\left(E_{j}-E_{i}\right)}&0\end{pmatrix}\oplus\mathbb{I}_{\backslash\left(i,j\right)}.$
(33)

where $\mathbb{I}_{\backslash(i,j)}$ is the identity on every level $l\neq
i,j$.

We now define the following Gibbs-stochastic matrix, which is a many-level
generalization of the matrix $A_{6}$ from Ref. [Mazurek and Horodecki(2018)]:

$\displaystyle\tilde{G}=\begin{pmatrix}1-e^{\beta\Delta_{01}}&1-e^{\beta\Delta_{12}}&...&1-e^{\beta\Delta_{d-1\,d-2}}&1\\\
e^{\beta\Delta_{01}}&0&...&0&0\\\ 0&e^{\beta\Delta_{12}}&\ddots&0&0\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\ 0&0&...&0&0\\\
0&0&...&e^{\beta\Delta_{d-1\,d-2}}&0\end{pmatrix}.$

where $\Delta_{ij}=E_{i}-E_{j}$. $\tilde{G}$ can be generated as a sequence of
$\beta$-swaps: $\tilde{G}=\prod_{i=0}^{d-2}\beta_{i+1,i}$. Now consider the
protocol in which, at each round, we perform the following steps:

  1. 1.

A unitary transformation $U$ is performed, that swaps the population of the
ground state and the most excited state.

  2. 2.

$\tilde{G}$ is performed.

Note that $U$ induces the permutation $S_{(0,d-1)}$ that flips the states
$0\leftrightarrow d-1$. Next, define the resulting cooling stochastic matrix
$C=\tilde{G}S_{(0,d-1)}$. Also denote by
$\Omega=\sum^{d-2}_{i=0}(E_{i+1}-E_{i})=E_{d-1}-E_{0}$. Then one can verify:

$\displaystyle
C^{d-1}=\begin{pmatrix}1&1-e^{-\beta\Omega}&1-e^{-\beta\Omega}&\dots&1-e^{-\beta\Omega}&\\\
0&e^{-\beta\Omega}&0&\dots&0\\\ 0&0&e^{-\beta\Omega}&\ddots&0\\\
\vdots&\vdots&0&\ddots&\vdots\\\
0&0&\dots&\dots&e^{-\beta\Omega}\end{pmatrix}.$

Let $p^{(k)}_{i}$ be occupations after $k$ steps of the protocol, and
$R^{(k)}=1-p^{(k)}_{0}$. One has

$p^{(k\left(d-1\right))}_{0}=1-R^{(k\left(d-1\right))},\quad
R^{(k\left(d-1\right))}=e^{-\beta\Omega k}S^{(1)}.$

By assumption $\Omega>0$ and $\beta>0$, and Eq. (2) follows.

###  2.6 Proof of trivial Markovian cooling for $d=2$

Markovian dephasing thermalizations are defined as those dephasing
thermalizations that can be written as the solution of a master equation with
a (possibly time-dependent) generator in Lindblad form. Here we show that the
lowest achievable temperature by Markovian dephasing thermalizations without
ancillary systems is limited to that of the environment.

As discussed in Sec. 2.1 of the Methods, dephasing thermalizations act on the
population vector $\boldsymbol{p}$ as Gibbs-stochastic matrices, in this case
$2\times 2$. Using the results of Ref. [Wolf et al.(2008)Wolf, Eisert, Cubitt,
and Cirac], a direct computation shows that the action on $\boldsymbol{p}$ of
Markovian dephasing thermalizations can be written as

$G_{\rm Mark}=(1-\lambda)\mathbb{I}+\lambda\beta_{01},\quad 0\leq\lambda\leq
1/(1+e^{-\beta E})$ (34)

with $\beta_{01}$ given by Eq. (33).

Let $(p,1-p)$, $(q,1-q)$ $(s,1-s)$ be the energy distribution at the start,
after the unitary, and after the Markovian dephasing thermalization,
respectively. A direct calculation shows that the thermalization relates
$\boldsymbol{s}$ and $\boldsymbol{q}$ by

$s=\left(1-\lambda e^{-\beta E}\right)q+\lambda(1-q).$ (35)

Since $\boldsymbol{p}$ and $\boldsymbol{q}$ are related by a unitary, we have
that $p(1-p)\leq q(1-q)$. Assume that the system is initially hotter than the
bath, $p<1/(1+e^{-\beta E})$. By direct inspection the maximum is found at
$\lambda=1/(1+e^{-\beta E})$, where one achieves $s=1/(1+e^{-\beta E})$, the
thermal ground state population at temperature $\beta$. The same would be true
if one optimized over all Markovian thermal operations (using again the
results of Ref. [Wolf et al.(2008)Wolf, Eisert, Cubitt, and Cirac] and the
proof in Sec. 2.9).

###  2.7 Proof of Theorem 2

We prove this theorem via a direct calculation of the ground state probability
after a repeated application of the unitaries. Note that the state
$\rho^{(k)}_{S}$ at the beginning of round $k+1$ of the optimal protocol is
incoherent for every $k\geq 1$. As such, the system-bath state after round $k$
will have the general form:

$\rho^{\left(k\right)}_{SB}=\sum_{n=0}^{\infty}\left(p^{\left(k\right)}_{0,n}|0,n\rangle\langle
0,n|+p^{\left(k\right)}_{1,n}|1,n\rangle\langle 1,n|\right),$

for some occupation probabilities $p^{(k)}_{i,n}$ with
$i\in\left\\{0,1\right\\}$ and $n\in\left\\{0,1,\dots\right\\}$. Applying $X$
to $S$, followed by $U^{\beta}_{SB}$ gives

$\displaystyle\rho^{\left(k\right)}_{SB}\mapsto\rho^{\left(k+1\right)}_{SB}=p^{\left(k\right)}_{1,0}|0,0\rangle\langle
0,0|+p^{\left(k\right)}_{0,0}|0,1\rangle\langle 0,1|$
$\displaystyle+\sum_{n=1}^{\infty}\left(p^{\left(k\right)}_{1,n}|1,n-1\rangle\langle
1,n-1|+p^{\left(k\right)}_{0,n}|0,n+1\rangle\langle 0,n+1|\right).$

The relation between the occupations at step $k$ and those at step $k+1$ is
(see Fig. 9)

$\displaystyle p^{\left(k+1\right)}_{0,0}$
$\displaystyle=p^{\left(k\right)}_{1,0},$ $\displaystyle
p^{\left(k+1\right)}_{0,n}$
$\displaystyle=p^{\left(k\right)}_{0,n-1},\quad\quad n\geq 1,$ $\displaystyle
p^{\left(k+1\right)}_{1,n}$
$\displaystyle=p^{\left(k\right)}_{1,n+1},\quad\quad n\geq 0.$

Figure 9: The circulation of populations induced by each round (Pauli $X$
followed by $U^{\beta}_{SB}$) of the optimal protocol. Arrows indicate
complete transfer of population. The picture gives an intuitive understanding
of the cooling mechanism of the optimal protocol.

Solving these equations recursively, we find the populations at step $k$ as a
function of the initial occupation probabilities (at $k=0$):

$\displaystyle p^{\left(k\right)}_{1,n}$
$\displaystyle=p^{\left(0\right)}_{1,n+k}$ $\displaystyle
p^{\left(k\right)}_{0,n}$
$\displaystyle=\begin{cases}p^{\left(0\right)}_{0,n-k},\quad&n\geq k,\\\
p^{\left(0\right)}_{1,k-n+1},\quad&n<k.\end{cases}$

As before, we denote by $p^{(k)}_{0}$ the ground state population of the
system after $k$ steps of the protocol. We have

$\displaystyle
p^{(k)}_{0}=p_{0,0}^{(k)}+\sum_{n=1}^{\infty}p^{(k)}_{0,n}=p^{(k-1)}_{1,0}+\sum_{n=1}^{\infty}p^{(k-1)}_{0,n-1}=p^{(k-1)}_{1,0}+p^{(k-1)}_{0}.$

We now use the relation

$p^{(k-1)}_{1,0}=p^{(k-2)}_{1,1}=\dots=p^{(0)}_{1,k-1}=p^{(0)}t^{(0)}_{k},$

where $t^{(0)}_{k}=\left(1-e^{-\beta E}\right)e^{-k\beta E}$ denotes the
occupation of the $k$th energy level of the bosonic mode at the beginning of
the protocol. This finally leads to

$p^{\left(k\right)}_{0}=p^{\left(0\right)}_{0}+\left(1-p^{\left(0\right)}\right)\sum_{n=0}^{k-1}t^{\left(0\right)}_{n}.$
(36)

Direct substitution shows that this expression is identical to that derived in
Theorem 1 for the optimal protocol.

###  2.8 Robustness of Theorem 2 to anharmonicity

We model an anharmonicity in the oscillator as a correction to the energy gap
between the energy levels $n$ and $n+1$ (which we now label as
$E^{\text{an}}_{n}$)

$E^{\text{an}}_{n+1}-E^{\text{an}}_{n}=E(1-(n+1)\tau^{2})+o(\tau^{2}),$ (37)

where $\tau\ll 1$. To analize its effect, notice that the ground state from
Eq. (36) is

$p^{\left(k\right)}_{0}=p^{\left(0\right)}_{0}+\left(1-p^{\left(0\right)}\right)\sum_{n=0}^{k-1}t^{\left(0\right)}_{n},$
(38)

and that this does not rely on any particular form for the populations
$t_{k}^{(0)}$. Thus, for any $\tau$, we have to compare the two following
sums,

$\sum_{n=0}^{k-1}t^{\left(0\right),\text{an}}_{n}=\frac{\sum_{n=0}^{k-1}e^{-\beta
E^{\text{an}}_{n}}}{\sum_{n=0}^{\infty}e^{-\beta
E^{\text{an}}_{n}}}\quad\sum_{n=0}^{k-1}t^{\left(0\right)}_{n}=\frac{\sum_{n=0}^{k-1}e^{-\beta
E_{n}}}{\sum_{n=0}^{\infty}e^{-\beta E_{n}}}.$ (39)

A significantly large value $\tau=0.05$ gives a distortion $E^{\rm
an}_{n}/E_{n}$ of $5\%$ at $n=8$. Also,
$\frac{\sum_{n=0}^{k-1}t^{\left(0\right),\text{an}}_{n}}{\sum_{n=0}^{k-1}t^{\left(0\right)}_{n}}$
approaches $1$ as $k$ grows and the largest deviation from
$\frac{\sum_{n=0}^{k-1}e^{-\beta E_{n}}}{\sum_{n=0}^{\infty}e^{-\beta E_{n}}}$
peaks at less than $0.005\%$. This shows that these protocols are robust to
realistic anharmonicities.

###  2.9 Proof of Theorem 3 and its extension to thermal operations

We give a proof of a more general result than Theorem 3, optimizing the
thermalization step over the set of thermal operation (TO) $\mathcal{T}$ on a
qubit, which strictly include all dephasing thermalizations.

####  2.9.1 Preliminaries

First we review the definition of thermal operations [Janzing et
al.(2000)Janzing, Wocjan, Zeier, Geiss, and Beth, Brandão et al.(2013)Brandão,
Horodecki, Oppenheim, Renes, and Spekkens], which are channels acting on a
quantum system $S$ in state $\rho_{S}$ and with Hamiltonian $H_{S}$ as:

$\mathcal{T}\left(\rho_{S}\right)=\textnormal{Tr}_{B}\left[U\left(\rho_{S}\otimes\frac{e^{-\beta
H_{B}}}{\textnormal{Tr}\left[e^{-\beta
H_{B}}\right]}\right)U^{\dagger}\right],$ (40)

where $\beta=1/(kT)$ is the inverse temperature of the surrounding heat-bath,
$H_{B}$ is an arbitrary Hamiltonian and $U$ is an energy preserving unitary
that satisfies $\left[U,H_{S}+H_{B}\right]=0$.

Let $S$ be a qubit with Hamiltonian $H_{S}=E|1\rangle\langle 1|$. Given
$\rho_{S}$, in addition to the occupation probabilities
$\boldsymbol{p}_{S}=\left(p,1-p\right)^{T}$, let $\rho_{01}=\langle
0|\rho_{S}|1\rangle$ denote the off-diagonal term. Similarly, for a second
state $\sigma_{S}$, let $\boldsymbol{q}_{S}=\left(\langle
0|\sigma_{S}|0\rangle,\langle 1|\sigma_{S}|1\rangle\right)^{T}$ and
$\sigma_{01}=\langle 0|\sigma_{S}|1\rangle$. Then any thermal operation on a
qubit can be described by:

  1. 1.

The action on the population (given by a Gibbs-stochastic matrix $G$).

  2. 2.

The action on the off-diagonal term.

Taken together this gives:

$\mathcal{T}(\rho_{S})=\sigma_{S}\,\Leftrightarrow\begin{cases}\boldsymbol{q}_{S}=G\boldsymbol{p}_{S}=\begin{pmatrix}1-\lambda
e^{-\beta E}&\lambda\\\ e^{-\beta
E}\lambda&1-\lambda\end{pmatrix}\boldsymbol{p}_{S},\quad\lambda\in[0,1]\\\
|\sigma_{01}|=c|\rho_{01}|,\quad 0\leq c\leq\sqrt{(1-e^{-\beta
E}\lambda)(1-\lambda)}\end{cases},$ (41)

where the expression for $c$ follows from the constraint
$|c|\leq\sqrt{G_{00}G_{11}}$ [Ćwikliński et al.(2015)Ćwikliński, Studziński,
Horodecki, and Oppenheim, Lostaglio et al.(2015)Lostaglio, Korzekwa, Jennings,
and Rudolph] which is a consequence of the complete positivity of
$\mathcal{T}$. One can also apply a phase $e^{i\theta}$ to $\rho_{01}$ by a
unitary, so that in general $c\not\in\mathbb{R}$, but since this will be
irrelevant in the following considerations, without loss of generality we set
$\theta=0$ and take $c\geq 0$ (one can always reversibly transform to any
$\theta\neq 0$ by an energy-preserving unitary, which is a thermal operation).
Hence, for our purposes a qubit thermal operation $\mathcal{T}$ is defined by
the two parameters:

$\boldsymbol{v}\left[\mathcal{T}\right]=\left(\lambda,c\right).$ (42)

With this notation, the $\beta$_-swap_ introduced in the main text is a
thermal operation $\mathcal{T}_{\beta}$ that removes all quantum coherences
and has $G=G^{\beta\textrm{-swap}}$ where

$G^{\beta\textrm{-swap}}=\begin{pmatrix}1-e^{-\beta E}&1\\\ e^{-\beta
E}&0\end{pmatrix}$ (43)

corresponds to
$\boldsymbol{v}\left[\mathcal{T}_{\beta}\right]=\left(1,0\right)$. Dephasing
thermalizations are the subset of thermal operations with $c=0$.

####  2.9.2 Proof of an extended version of Theorem 3

Define the set of $\epsilon$-noisy thermal operations as the set of thermal
Operation $\mathcal{T}_{\epsilon}$ such that

$\boldsymbol{v}[\mathcal{T}_{\epsilon}]=(\lambda,c),\quad\lambda\leq
1-\epsilon.$ (44)

Let $\mathcal{T}^{\epsilon}_{\emptyset}$ denote the set of cooling protocols
using no ancilla in which at each round $k$

  1. 1.

A unitary $U^{(k)}$ is applied to $S$,

  2. 2.

An $\epsilon$-noisy thermal operation $\mathcal{T}^{(k)}_{\epsilon}$ is
applied to $S$.

  3. 3.

A unitary $V^{(k)}$ is applied to $S$

Then,

######  Theorem 4.

Under the assumptions of Corollary 1 and given
$\epsilon\leq\frac{1}{1+e^{\beta E}+e^{2\beta E}}$, the optimal nontrivial
cooling protocol in $\mathcal{T}^{\epsilon}_{\emptyset}$ is such that in each
round $k$:

  1. 1.

The Pauli $X$ unitary is applied to $S$.

  2. 2.

Any $\epsilon$-noisy thermal operation $\mathcal{T}_{\epsilon}$ with
$\boldsymbol{v}(\mathcal{T}_{\epsilon})=(1-\epsilon,c)$ is applied to $S$.

The population of the ground state after round $k$ is:

$\displaystyle p_{\rm
0}^{(k)}=1-\frac{\epsilon}{2-(1-\epsilon)Z}\,\,-\left((1-\epsilon)Z-1\right)^{k}\left(1-\frac{\epsilon}{2-(1-\epsilon)Z}-p_{0}^{(0)}\right),$

where $Z=1+e^{-\beta E}$ and $p^{\left(k\right)}_{0}\rightarrow
1-\frac{\epsilon}{2-(1-\epsilon)Z}$ as $k\rightarrow\infty$.

Before we prove this theorem, let us discuss the consequences. Note that a
particular choice for the optimal cooling protocol is $\mathcal{T}_{\epsilon}$
with $\boldsymbol{v}(\mathcal{T}_{\epsilon})=(1-\epsilon,0)$. This is just the
$\epsilon$-noisy $\beta$-swap defined in the main text. Since
$\mathcal{T}^{\epsilon}_{\emptyset}\subset\mathcal{P}^{\epsilon}_{\emptyset}$,
this means that the above theorem implies Theorem 3 of the main text as an
immediate corollary. Also note, that setting $\epsilon=0$, we obtain as a
simple consequence a direct proof of Corollary 1 that does not goes through
Theorem 1.

Furthermore, note it does not matter what is the choice for $c$ in the
$\epsilon$-noisy thermal operation: it could even vary from round to round.
The Jaynes-Cummings implementation described in the main text is an
$\epsilon$-noisy thermal operation. On the population it acts as an
$\epsilon$-noisy $\beta$-swap, but it has $c\neq 0$. However, the above
theorem shows that this does not make any difference, since having control on
the coherent part of the evolution does not provide any advantage and the
performance is independent by the choice of $c$ at each step. In other words,
as claimed in the main text, we can safely ignore the coherent evolution and
focus on the induced population dynamics, at least for $\epsilon$ small
enough.

###### Proof.

The most general protocol can be written as

$\rho^{(k)}_{S}=\mathcal{V}^{(k)}\mathcal{T}^{(k)}\mathcal{U}^{(k)}\dots\mathcal{V}^{(1)}\mathcal{T}^{(1)}\mathcal{U}^{(1)}\mathcal{U}_{\textrm{di}}(\rho_{S}),$
(45)

where $\mathcal{U}^{(i)}(\cdot)=U^{(i)}(\cdot)U^{(i)\,\dagger}$ and
$\mathcal{V}^{(i)}(\cdot)=V^{(i)}(\cdot)V^{(i)\,\dagger}$ denote general
unitaries applied to $S$ in round $i$ and $\mathcal{T}^{(i)}$ is a thermal
operation on $S$ applied at round $i$. Here
$\mathcal{U}_{\textrm{di}}=U_{\textrm{di}}(\cdot)U^{\dagger}_{\textrm{di}}$ is
the unitary that transforms $\rho_{S}$ into a diagonal form with
$p_{0}^{(0)}\geq p_{1}^{(0)}$. This can be done without loss of generality,
since any protocol in $\mathcal{P}$ starts with an arbitrary unitary
$\mathcal{U}^{(1)}$.

We want to maximize the ground state population over all choices of
$\mathcal{U}^{(i)}$, $\mathcal{T}^{(i)}$ and $\mathcal{V}^{(i)}$. As we
perform an arbitrary unitary at the beginning of every round, we can assume
without loss of generality that the state of $S$ at the start of round $k+1$
is diagonal in the energy eigenbasis:

$\rho^{\left(k\right)}_{S}=\begin{pmatrix}p^{\left(k\right)}&0\\\
0&1-p^{\left(k\right)}\end{pmatrix}$

and that $p^{\left(k\right)}\geq\frac{1}{2}$ (here we drop the subscript 0 to
simplify the notation). In round $k+1$ of the protocol, we have:

$\displaystyle\rho^{\left(k\right)}_{S}=\begin{pmatrix}p^{\left(k\right)}&0\\\
0&1-p^{\left(k\right)}\end{pmatrix}\stackrel{{\scriptstyle\mathcal{U}^{\left(k\right)}}}{{\longrightarrow}}$
$\displaystyle\rho^{\prime\left(k\right)}_{S}=\begin{pmatrix}q^{\left(k\right)}&a^{(k)}\\\
a^{*(k)}&1-q^{\left(k\right)}\end{pmatrix}$
$\displaystyle\stackrel{{\scriptstyle\mathcal{T}^{\left(k\right)}}}{{\longrightarrow}}$
$\displaystyle\rho^{\prime\prime\left(k\right)}_{S}=\begin{pmatrix}s^{\left(k\right)}&b^{(k)}\\\
b^{*(k)}&1-s^{\left(k\right)}\end{pmatrix}$
$\displaystyle\stackrel{{\scriptstyle\mathcal{V}^{\left(k\right)}}}{{\longrightarrow}}$
$\displaystyle\rho^{\left(k+1\right)}_{S}=\begin{pmatrix}p^{\left(k+1\right)}&0\\\
0&1-p^{\left(k+1\right)}\end{pmatrix}$

and our goal is to maximize $p^{\left(k+1\right)}$.

As $\rho^{\left(k+1\right)}_{S}$ and $\rho^{\prime\prime\left(k\right)}_{S}$
are related by a unitary, maximizing $p^{\left(k+1\right)}$ corresponds to
minimizing the determinant of $\rho^{\prime\prime\left(k\right)}_{S}$. Among
all thermal operations associated to a fixed Gibbs-stochastic matrix $G$ with
$G\boldsymbol{q}^{(k)}:=\boldsymbol{s}^{(k)}$ (with
$\boldsymbol{q}^{\left(k\right)}$ and $\boldsymbol{s}^{\left(k\right)}$
defined through the above matrices), optimality of the protocol imposes that
we choose $\mathcal{T}^{(k)}$ to be a thermal operation that maximizes the
absolute value of $b^{(k)}$ \- i.e. it preserves the maximum possible amount
of coherence. This is achieved by the thermal operation that maximizes the
parameter $c$ for given $\lambda$, i.e., from Eq. (41), $c=(1-\lambda
e^{-\beta E})(1-\lambda)$ or, with the parametrization of Eq. (42),
$\boldsymbol{v}[{\mathcal{T}^{(k)}}]=(\lambda,(1-\lambda e^{-\beta
E})(1-\lambda))$.

What remains to be done is to perform an optimization over all possible Gibbs-
stochastic matrices $G$, parametrized by $\lambda\in[0,\lambda_{\max}]$ as in
Eq. (41), where $\lambda_{\max}\geq 1-\frac{1}{1+e^{\beta E}+e^{2\beta E}}$
(which corresponds to $\epsilon\leq\frac{1}{1+e^{\beta E}+e^{2\beta E}}$). For
each $G$, the relation between $s^{\left(k\right)}$, $b^{(k)}$ and
$q^{\left(k\right)}$ and $a^{(k)}$ is

$\displaystyle s^{\left(k\right)}$ $\displaystyle=\left(1-\lambda e^{-\beta
E}\right)q^{\left(k\right)}+\lambda\left(1-q^{\left(k\right)}\right)$
$\displaystyle|b^{\left(k\right)}|^{2}$
$\displaystyle=|a^{\left(k\right)}|^{2}\left(1-\lambda e^{-\beta
E}\right)\left(1-\lambda\right).$

Finally, using the unitarity of $U^{\left(k\right)}$, we can relate
$q^{\left(k\right)}$ and $a^{(k)}$ to $p^{\left(k\right)}$ via:

$p^{\left(k\right)}\left(1-p^{\left(k\right)}\right)=q^{\left(k\right)}\left(1-q^{\left(k\right)}\right)-\left|a^{(k)}\right|^{2}.$

The determinant of $\rho^{\prime\prime\left(k\right)}_{S}$ is thus given by

$\displaystyle
f_{p^{\left(k\right)},E}\left(q^{\left(k\right)},\lambda\right)=$
$\displaystyle\left(\lambda-1\right)\left(e^{-\beta
E}\lambda-1\right)\left(-\left(p^{\left(k\right)}-q^{\left(k\right)}\right)\right)\left(p^{\left(k\right)}+q^{\left(k\right)}-1\right)$
$\displaystyle\quad-\left(q^{\left(k\right)}\left(e^{-\beta
E}\lambda+\lambda-1\right)-\lambda\right)\left(q^{\left(k\right)}\left(e^{-\beta
E}\lambda+\lambda-1\right)-\lambda+1\right),$ (46)

for fixed $E$ and $p^{\left(k\right)}$. Note that the equation is quadratic in
both $\lambda$ and $q^{\left(k\right)}$. The coefficient of $\lambda^{2}$ is

$-\left[\left(1-\left(q^{\left(k\right)}\right)^{2}\right)-q^{\left(k\right)}e^{-\beta
E}\right]^{2}-e^{-\beta
E}\left[q^{\left(k\right)}\left(1-q^{\left(k\right)}\right)-p^{\left(k\right)}\left(1-p^{\left(k\right)}\right)\right],$

which is clearly negative as
$q^{\left(k\right)}\left(1-q^{\left(k\right)}\right)=p^{\left(k\right)}\left(1-p^{\left(k\right)}\right)+|a^{\left(k\right)}|^{2}$,
and thus the minimum values will be obtained at either $\lambda=0$ or
$\lambda=\lambda_{\max}$. The case $\lambda=0$ corresponds to not implementing
the TO and leads to $p^{\left(k+1\right)}=p^{\left(k\right)}$. Taking
$\lambda=\lambda_{{\rm max}}$, we still need to show that the Pauli $X$
unitary at each step is optimal.

Let us rewrite the function
$f_{p^{\left(k\right)},E}\left(q^{\left(k\right)},\lambda_{\rm max}\right)$ as
a polynomial in $q^{(k)}$,

$f_{p^{\left(k\right)},E}\left(q^{\left(k\right)},\lambda_{\rm
max}\right)\equiv f^{(2)}_{p^{\left(k\right)},E}\left(\lambda_{\rm
max}\right)(q^{(k)})^{2}+f^{(1)}_{p^{\left(k\right)},E}\left(\lambda_{\rm
max}\right)q^{(k)}+f^{(0)}_{p^{\left(k\right)},E}\left(\lambda_{\rm
max}\right).$

Given that it is a quadratic equation, the location of its minimum depends on
the sign of the coefficient $f^{(2)}_{p^{\left(k\right)},E}\left(\lambda_{\rm
max}\right)=\lambda_{\rm max}(1+e^{-\beta E}-\lambda_{\rm max}(1+e^{-\beta
E}+e^{-2\beta E}))$. To see when the Pauli $X$ is optimal, we want to find the
cases in which the solution is either $q^{(k)}=p^{(k)}$ or $q^{(k)}=1-p^{(k)}$
(that is, on the boundary of the range, given unitarity), which occurs when
$f^{(2)}_{p^{\left(k\right)},E}\left(\lambda_{\rm max}\right)<0$. This is
equivalent to

$\lambda_{{\rm max}}>\frac{1+e^{-\beta E}}{1+e^{-\beta E}+e^{-2\beta
E}}=1-\frac{1}{1+e^{\beta E}+e^{2\beta E}},$ (47)

which is true by assumption. On top of this, we find that

$f_{p^{\left(k\right)},E}\left(p^{\left(k\right)},\lambda_{\rm
max}\right)-f_{p^{\left(k\right)},E}\left(1-p^{\left(k\right)},\lambda_{\rm
max}\right)=\lambda_{{\rm max}}(1-e^{-\beta E})(\lambda_{{\rm
max}}(1+e^{-\beta E})-1)(2p^{\left(k\right)}-1),$

so the minimum is at
$q^{(k)}=\min\left\\{p^{\left(k\right)},1-p^{\left(k\right)}\right\\}=1-p^{(k)}$
as long as $\lambda_{{\rm max}}>\frac{1}{1+e^{-\beta E}}$ (which is implied by
Eq. (47)).

Thus, the optimal protocol in $\mathcal{T}^{\epsilon}_{\emptyset}$ is one in
which, for every $k$, $\mathcal{T}^{(k)}$ is a thermal operation
$(\lambda_{{\rm max}},c_{k})$ with $\lambda_{\text{max}}=1-\epsilon$ (the
“best approximation” of the $\beta$-swap) and
$\mathcal{U}^{(k)}(\cdot)=\mathcal{X}(\cdot):=X(\cdot)X^{\dagger}$. Note that
we do not specify how each $\mathcal{T}^{(k)}$ acts on the off-diagonal
element of the quantum state (i.e., the parameter $c_{k}$) simply because the
input state contains no coherence. As such, the protocol is also optimal for
the set of dephasing thermalizations, since at step $k$ one can perform any TO
with $G_{01}=\lambda_{{\rm max}}$, without any control required on $c_{k}$.
The optimal ground state population achieved by the above protocol satisfies

$p^{\left(k+1\right)}=\left(1-\lambda_{{\rm max}}e^{-\beta
E}\right)\left(1-p^{\left(k\right)}\right)+\lambda_{{\rm
max}}p^{\left(k\right)},$

as one can verify by a direct computation. Solving this recursion relation
gives Eq. (4). One recovers the scaling of Theorem 1 when $\lambda_{{\rm
max}}=1$. Furthermore, since $\lambda_{{\rm max}}Z-1\leq e^{-\beta E}$ one has
exponential convergence to
$1-\frac{1-\lambda_{\textrm{max}}}{1-2\lambda_{\textrm{max}}Z}$. Note that the
trivial protocol in which we do not do a thermal operation is optimal only
when $p^{(0)}\geq
1-\frac{1-\lambda_{\textrm{max}}}{1-2\lambda_{\textrm{max}}Z}$ (that is, when
the initial ground state population is higher than the optimal asymptotic
value).

Finally, let us show that the concatenation of optimal rounds is optimal
overall. To this end, let $p^{(k)}$ and $\tilde{p}^{(k)}$ be two ground state
populations with $p^{(k)}\geq\tilde{p}^{(k)}$. One can compute

$p^{(k+1)}-\tilde{p}^{(k+1)}\geq(\lambda_{\textrm{max}}Z-1)(p^{(k)}-\tilde{p}^{(k)})\geq
0,$

since $\lambda_{\textrm{max}}>1/Z$. Hence the optimal protocol is a
concatenation of the optimal single round protocol, and is also independent of
the initial state of the system (with exclusion of the $k=0$ unitary). ∎

##  3 Experimental proposal

###  3.1 Bounds on the achievable cooling in the Jaynes-Cummings model

We compare the performance of a protocol that uses a Jaynes-Cummings (JC)
interaction to implement the $\beta$-swap against the ideal unitary
$U^{\beta}_{SB}$ of Eq. (5). Consider a resonant JC Hamiltonian in rotating
wave approximation

$H_{\textrm{JC}}=g(\sigma_{+}\otimes a+\sigma_{-}\otimes a^{\dagger}),$ (48)

coupling $S$ with a single-mode bosonic bath prepared in a thermal state. Here
$a^{\dagger}$ and $a$ are creation and annihilation operators on $B$ and
$\sigma_{+}=|1\rangle\langle 0|$ and $\sigma_{-}=|0\rangle\langle 1|$. The de-
excitation probability is

$G_{0|1}(s)=\frac{1}{Z_{B}}\sum_{n=1}^{\infty}\sin^{2}(s\sqrt{n})e^{-\beta
E(n-1)},$ (49)

where $Z_{B}=(1-e^{-\beta E})^{-1}$ and $s$ is the normalized interaction time
($s=gt$, if $U_{JC}=e^{-iH_{JC}t}$). To realize a dephasing thermalization,
strictly speaking one should dephase in the energy basis. However, since in
the JC interaction the evolutions of populations and coherences are decoupled
and the subsequent Pauli $X$ simply inverts the populations, we can skip this
step (furthermore, the results of Sec. 2.9 show that these coherences cannot
be exploited by changing the unitary step). Hence at each round $k$:

  1. 1.

A Pauli $X$ operation is performed on $S$.

  2. 2.

The interaction Hamiltonian $H_{\rm JC}$ operates for a time $s$ that
maximizes Eq. (49).

  3. 3.

The bosonic mode is reset to a thermal state.

In Ref. [Lostaglio et al.(2018)Lostaglio, Alhambra, and Perry], temperature-
dependent upper bounds for the maximum achievable transition probability
$G^{{\rm max}}(s)$ were derived as a function of $\bar{\beta}=\beta E$:

$\displaystyle G^{{\rm
max}}(\bar{\beta})\leq\left\\{\begin{array}[]{ll}\frac{1}{16}\left(8e^{-\bar{\beta}}-e^{2\bar{\beta}}+e^{3\bar{\beta}}+8\right),\;\textrm{for
}\bar{\beta}\in[0,\frac{\log(4)}{3}],\\\
e^{-4\bar{\beta}}-e^{-3\bar{\beta}}+1,\quad\quad\quad\quad\quad\textrm{for
}\bar{\beta}\geq\frac{\log(4)}{3}.\end{array}\right.$

Substituting these upper bounds in Theorem 3, we obtain upper bounds on the
ground state population achieved in the JC models after $k$ rounds, presented
in Fig. 3. In particular for the asymptotic population

$p^{(\infty)}_{0}\leq\left\\{\begin{array}[]{ll}\frac{1}{e^{\bar{\beta}}+\frac{16e^{2\bar{\beta}}}{-16e^{\bar{\beta}}+e^{3\beta}-8}+1},\;\textrm{for
}\bar{\beta}\in[0,\frac{\log(4)}{3}],\\\
\frac{1}{\frac{e^{\bar{\beta}}}{e^{4\bar{\beta}}+1}+1},\quad\quad\quad\quad\quad\textrm{for
}\bar{\beta}\geq\frac{\log(4)}{3}.\end{array}\right.$ (50)

To obtain an explicit protocol whose performance lower bounds the optimal, one
can numerically optimize $G_{0|1}(s)$ over $s$ within a finite domain. We take
$s\in[0,5\times 10^{3}]$, and obtain the curves in Fig. 3. As an illustrative
example, for $\bar{\beta}=1$ one has that the optimal JC asymptotic cooling is
$p^{(\infty)}_{0}\in[0.9401,0.9534]$. As we know from Theorem 3, the
convergence is exponential and can be computed explicitly from Eq. (4) with
$\epsilon=1-G_{0|1}(\tilde{s})$ for the chosen interaction time $\tilde{s}$.

While a detailed analysis is beyond the scope of the current work and would
require fixing specific experimental parameters, it is worth noticing that,
excluding the high temperature regime, the JC protocol appears to be superior
to the ideal PPA protocol with 2 ancillas even taking into account two kinds
of time limitation:

  1. 1.

Limited waiting time in the cavity, i.e. a bound on the maximum available $s$;

  2. 2.

Limited timing accuracy, i.e. the achieved $s$ fluctuates around a target
value.

As a case study, set $\bar{\beta}=1$ and limit $s$ to $s=gt\leq 10$. Then the
best available approximation to the $\beta$-swap is realized for $s\approx
7.87$, with the performance monotonically decreasing in a neighbourhood of
this value. We then consider an accurate ‘time-limited JC model’ where we set
$s=7.87$; and on top of this we allow various degrees of inaccuracy: errors on
$s$ of $\pm 0.1$, $\pm 0.2$ and $\pm 0.3$. The result are summarized in Fig.
10. The time-limited JC performs better than PPA with 2 ancillas. Adding
timing errors on $s$, we obtain a worst-case cooling performance above PPA
till around $\pm 0.2$.

Figure 10: In the time-limited JC model (blue curve on top) we optimise the
interaction time $s$ over the limited period $s\in[0,10]$. To this we add
increasing errors in the accuracy with which $s$ is achieved: $\pm 0.1$, $\pm
0.2$ and $\pm 0.3$ (dashed curves in red, orange and green). We see that the
JC models still outperform PPA till an accuracy on $s$ of around $\pm 0.2$.

###  3.2 Jaynes-Cummings model without refreshing the thermal mode

In the previous result the mode is always reset back to the thermal state
after every step. Theorem 2 shows that when using $U^{\beta}_{SB}$ the same
mode can be used repeatedly. Does this still hold (at least approximately)
when the interaction is via the Jaynes-Cummings model?

We explore this question via a numerical simulation of a modified algorithm
where now, at each round $k$, instead of rethermalizing the mode completely,
the bosonic mode is partially re-thermalized via a dissipation process. This
is done with a standard master equation, which models the evolution of state
$\rho_{B}$ of the cavity mode due to the interaction with an external thermal
field [Scala et al.(2007)Scala, Militello, Messina, Piilo, and Maniscalco]:

$\displaystyle\frac{\text{d}\rho_{B}}{\text{d}t}=-iE\left[a^{\dagger}a,\rho_{B}\right]-\frac{1}{2}An\left[aa^{\dagger}\rho_{B}-2a^{\dagger}\rho_{B}a+\rho_{B}aa^{\dagger}\right]-\frac{1}{2}A\left(n+1\right)\left[a^{\dagger}a\rho_{B}-2a\rho_{B}a^{\dagger}+\rho_{B}a^{\dagger}a\right].$
(51)

Here $A$ is the rate of loss of cavity photons (controlling the strength of
the re-thermalization) and $n=\frac{1}{e^{\beta E}-1}$ is the average number
of reservoir quanta with energy $E$.

We compare the cooling achieved after $k$ rounds in the case of full reset at
each round (Eq. (4) with $\epsilon=1-G_{0|1}(\tilde{s})$), with the cooling
achieved for various finite re-thermalization times. In each case, we fix
$\beta E=1$ and a particular interaction time between the qubit and the cavity
mode, $\tilde{s}=g\tilde{t}=98.92$. Note that the same procedure can be
applied for any $s$, or even taking $s$ to be a random variable to simulate
imperfections in the timing. The results are shown in Fig. 11. We find that,
unlike in the case of implementing $U^{\beta}_{SB}$, one has to reset the mode
back to the thermal state as much as possible in order for the algorithm to
work efficiently with a fixed interaction time.

Figure 11: Cooling achieved by the Pauli/Jaynes-Cummings protocols as a
function of the number of rounds, when the mode is re-thermalized for various
times $t_{\text{th}}$. The infinite time corresponds to the mode being
completely re-thermalized at every step (light blue curve at the top), while
$t_{\text{th}}=0$ corresponds to no re-thermalization at all (red curve with
wide oscillations). We see that re-thermalization is needed to achieve the
greatest cooling. For this figure, the parameters are $\beta E=1$, $A=1$,
$g=1$ and the Jaynes-Cummings interaction is turned on for a period of time
$\tilde{t}=98.92$.

Reasonably high cooling can be achieved by interrupting the protocol after $2$
rounds. The above findings suggest the experimental setup discussed in the
main text, where atoms are slowly fired inside two identical cavities resonant
with the qubit transitions we are cooling. The protocol on each atom then
consists of:

  1. 1.

First Pauli $X$ applied.

  2. 2.

First JC interaction applied, for some time $\tilde{s}$.

  3. 3.

Second Pauli $X$ applied.

  4. 4.

Second JC interaction applied, for the same time $\tilde{s}$.

  5. 5.

The cavity modes undergo re-thermalization for a finite time according to Eq.
(51).

While in the first step the JC interaction achieves a de-excitation
probability of $G_{0|1}(\tilde{s})$ given by Eq. (49), every subsequent atom
interacts with only partially re-thermalized cavities, for which Eq. (49) does
not hold. Thus, the protocol may not achieve the same cooling on every atom.
Nevertheless, one may expect that, by firing the atoms slowly enough, the re-
thermalization of the cavities due to losses will be sufficient to make the
cooling performance almost constant. This intuition is confirmed in Fig. 5. We
take the atoms to be initially in a thermal state with $\beta E=1$. We then
plot the final ground state population achieved by each atom passing through
the two cavities, as a function of the number of atoms already cooled. The
various curves represent different choices for the ratio $A/r$ between the
strength of the re-thermalization and the rate at which the atoms are fired
(with the caveat that we assume $r$ small enough so that two atoms are never
present at the same time in a cavity). When $A/r=\infty$, the single mode has
time to re-thermalize perfectly, the performance is the same for each atom and
given by Eq. (4) with $k=2$ and $\epsilon=1-G_{0|1}(98.92)$. In realistic
scenarios, however, we see that the incomplete thermalization negatively
impacts upon the performance. However, we also see that the performance
stabilizes to a constant after a small number of atoms are fired, so for a
high enough ratio $A/r$ cooling of any number of atoms is possible. This may
be understood as the creation of a steady state in the cavity field.

The authors would like to thank Philippe Faist, Amikam Levy, Mohammad
Mehboudi, Nayeli A. Rodríguez-Briones, Raam Uzdin and Marcus Huber for useful
discussions and feedback. Research at Perimeter Institute is supported by the
Government of Canada through the Department of Innovation, Science and
Economic Development and by the Province of Ontario through the Ministry of
Research, Innovation and Science. ML acknowledges financial support from the
the European Union’s Marie Sklodowska-Curie individual Fellowships
(H2020-MSCA-IF-2017, GA794842), Spanish MINECO (Severo Ochoa SEV-2015-0522 and
project QIBEQI FIS2016-80773-P), Fundacio Cellex and Generalitat de Catalunya
(CERCA Programme and SGR 875). CP acknowledges financial support from the
European Research Council (ERC Grant Agreement no. 337603) and VILLUM FONDEN
via the QMATH Centre of Excellence (Grant no. 10059).

## References

  * [DiVincenzo(2000)] David P DiVincenzo.  The physical implementation of quantum computation.  _Fortschritte der Physik: Progress of Physics_ , 48(9-11):771–783, 2000.  doi: 10.1002/1521-3978(200009)48:9/11<771::AID-PROP771>3.0.CO;2-E.
  * [Liu et al.(2013)Liu, Xiao, Luan, and Wong] Yong-Chun Liu, Yun-Feng Xiao, Xingsheng Luan, and Chee Wei Wong.  Dynamic dissipative cooling of a mechanical resonator in strong coupling optomechanics.  _Phys. Rev. Lett._ , 110:153606, Apr 2013.  doi: 10.1103/PhysRevLett.110.153606.
  * [Schulman and Vazirani(1999)] Leonard J Schulman and Umesh V Vazirani.  Molecular scale heat engines and scalable quantum computation.  In _Proceedings of the thirty-first annual ACM symposium on Theory of computing_ , pages 322–329. ACM, 1999.  doi: 10.1145/301250.301332.
  * [Boykin et al.(2002)Boykin, Mor, Roychowdhury, Vatan, and Vrijen] P Oscar Boykin, Tal Mor, Vwani Roychowdhury, Farrokh Vatan, and Rutger Vrijen.  Algorithmic cooling and scalable nmr quantum computers.  _Proceedings of the National Academy of Sciences_ , 99(6):3388–3393, 2002.  doi: 10.1073/pnas.241641898.
  * [Eschner et al.(2003)Eschner, Morigi, Schmidt-Kaler, and Blatt] Jürgen Eschner, Giovanna Morigi, Ferdinand Schmidt-Kaler, and Rainer Blatt.  Laser cooling of trapped ions.  _JOSA B_ , 20(5):1003–1015, 2003.  doi: 10.1364/JOSAB.20.001003.
  * [Valenzuela et al.(2006)Valenzuela, Oliver, Berns, Berggren, Levitov, and Orlando] Sergio O Valenzuela, William D Oliver, David M Berns, Karl K Berggren, Leonid S Levitov, and Terry P Orlando.  Microwave-induced cooling of a superconducting qubit.  _Science_ , 314(5805):1589–1592, 2006.  doi: https://doi.org/10.1126/science.1134008.
  * [Lostaglio(2018)] Matteo Lostaglio.  Thermodynamic laws for populations and quantum coherence: A self-contained introduction to the resource theory approach to thermodynamics.  _arXiv preprint arXiv:1807.11549_ , 2018.
  * [Schulman et al.(2005)Schulman, Mor, and Weinstein] Leonard J Schulman, Tal Mor, and Yossi Weinstein.  Physical limits of heat-bath algorithmic cooling.  _Phys. Rev. Lett._ , 94(12):120501, 2005.  doi: 10.1103/PhysRevLett.94.120501.
  * [Rodríguez-Briones and Laflamme(2016)] Nayeli Azucena Rodríguez-Briones and Raymond Laflamme.  Achievable polarization for heat-bath algorithmic cooling.  _Phys. Rev. Lett._ , 116(17):170501, 2016.  doi: 10.1103/PhysRevLett.116.170501.
  * [Raeisi and Mosca(2015)] Sadegh Raeisi and Michele Mosca.  Asymptotic bound for heat-bath algorithmic cooling.  _Phys. Rev. Lett._ , 114(10):100404, 2015.  doi: 10.1103/PhysRevLett.114.100404.
  * [Rodriguez-Briones et al.(2017)Rodriguez-Briones, Li, Peng, Mor, Weinstein, and Laflamme] Nayeli A Rodriguez-Briones, Jun Li, Xinhua Peng, Tal Mor, Yossi Weinstein, and Raymond Laflamme.  Heat-bath algorithmic cooling with correlated qubit-environment interactions.  _New Journal of Physics_ , 19(11):113047, 2017\.  doi: 10.1088/1367-2630/aa8fe0.
  * [Overhauser(1953)] Albert W Overhauser.  Aw overhauser, phys. rev. 89, 689 (1953).  _Phys. Rev._ , 89:689, 1953.  doi: 10.1103/PhysRev.89.689.
  * [Lostaglio et al.(2018)Lostaglio, Alhambra, and Perry] Matteo Lostaglio, Álvaro M. Alhambra, and Christopher Perry.  Elementary Thermal Operations.  _Quantum_ , 2:52, February 2018.  ISSN 2521-327X.  doi: 10.22331/q-2018-02-08-52.
  * [Horodecki and Oppenheim(2013)] M. Horodecki and J. Oppenheim.  Fundamental limitations for quantum and nanoscale thermodynamics.  _Nat. Commun._ , 4:2059, June 2013.  doi: 10.1038/ncomms3059.
  * [Ticozzi and Viola(2014)] Francesco Ticozzi and Lorenza Viola.  Quantum resources for purification and cooling: fundamental limits and opportunities.  _Scientific Reports_ , 4:5192, 2014.  doi: 10.1038/srep05192.
  * [Silva et al.(2016)Silva, Manzano, Skrzypczyk, and Brunner] Ralph Silva, Gonzalo Manzano, Paul Skrzypczyk, and Nicolas Brunner.  Performance of autonomous quantum thermal machines: Hilbert space dimension as a thermodynamical resource.  _Phys. Rev. E_ , 94(3):032120, 2016.  doi: 10.1103/PhysRevE.94.032120.
  * [Brandão et al.(2013)Brandão, Horodecki, Oppenheim, Renes, and Spekkens] Fernando G. S. L. Brandão, Michał Horodecki, Jonathan Oppenheim, Joseph M. Renes, and Robert W. Spekkens.  Resource theory of quantum states out of thermal equilibrium.  _Phys. Rev. Lett._ , 111:250404, Dec 2013.  doi: 10.1103/PhysRevLett.111.250404.
  * [Naderi et al.(2005)Naderi, Soltanolkotabi, and Roknizadeh] MH Naderi, M Soltanolkotabi, and R Roknizadeh.  A theoretical scheme for generation of nonlinear coherent states in a micromaser under intensity-dependent jaynes-cummings model.  _The European Physical Journal D-Atomic, Molecular, Optical and Plasma Physics_ , 32(3):397, 2005.  doi: 10.1140/epjd/e2004-00197-8.
  * [Åberg(2014)] Johan Åberg.  Catalytic coherence.  _Phys. Rev. Lett._ , 113:150402, Oct 2014.  doi: 10.1103/PhysRevLett.113.150402.
  * [Filipowicz et al.(1986)Filipowicz, Javanainen, and Meystre] P Filipowicz, J Javanainen, and P Meystre.  Theory of a microscopic maser.  _Phys. Rev. A_ , 34(4):3077, 1986.  doi: 10.1103/PhysRevA.34.3077.
  * [Walther et al.(2006)Walther, Varcoe, Englert, and Becker] Herbert Walther, Benjamin TH Varcoe, Berthold-Georg Englert, and Thomas Becker.  Cavity quantum electrodynamics.  _Reports on Progress in Physics_ , 69(5):1325, 2006.  doi: 10.1088/0034-4885/69/5/R02.
  * [Rodrigues et al.(2007)Rodrigues, Imbers, and Armour] DA Rodrigues, J Imbers, and AD Armour.  Quantum dynamics of a resonator driven by a superconducting single-electron transistor: A solid-state analogue of the micromaser.  _Phys. Rev. Lett._ , 98(6):067204, 2007.  doi: 10.1103/PhysRevLett.98.067204.
  * [Ruch et al.(1980)Ruch, Schranner, and Seligman] Ernst Ruch, Rudolf Schranner, and Thomas H Seligman.  Generalization of a theorem by hardy, littlewood, and polya.  _Journal of Mathematical Analysis and Applications_ , 76(1):222–229, 1980.  doi: 10.1016/0022-247X(80)90075-X.
  * [Birkhoff(1946)] Garrett Birkhoff.  Tres observaciones sobre el algebra lineal.  _Univ. Nac. Tucumán Rev. Ser. A_ , 5:147–151, 1946\.
  * [Mazurek and Horodecki(2018)] Paweł Mazurek and Michał Horodecki.  Decomposability and convex structure of thermal processes.  _New Journal of Physics_ , 20(5):053040, 2018\.  doi: 10.1088/1367-2630/aac057.
  * [Wolf et al.(2008)Wolf, Eisert, Cubitt, and Cirac] Michael Marc Wolf, J Eisert, TS Cubitt, and J Ignacio Cirac.  Assessing non-markovian quantum dynamics.  _Physical review letters_ , 101(15):150402, 2008\.  doi: 10.1103/PhysRevLett.101.150402.
  * [Janzing et al.(2000)Janzing, Wocjan, Zeier, Geiss, and Beth] D. Janzing, P. Wocjan, R. Zeier, R. Geiss, and Th. Beth.  Thermodynamic cost of reliability and low temperatures: Tightening Landauer’s principle and the second law.  _Int. J. Theor. Phys._ , 39(12):2717–2753, 2000\.  doi: 10.1023/A:1026422630734.
  * [Ćwikliński et al.(2015)Ćwikliński, Studziński, Horodecki, and Oppenheim] Piotr Ćwikliński, Michał Studziński, Michał Horodecki, and Jonathan Oppenheim.  Limitations on the evolution of quantum coherences: Towards fully quantum second laws of thermodynamics.  _Phys. Rev. Lett._ , 115:210403, Nov 2015.  doi: 10.1103/PhysRevLett.115.210403.
  * [Lostaglio et al.(2015)Lostaglio, Korzekwa, Jennings, and Rudolph] Matteo Lostaglio, Kamil Korzekwa, David Jennings, and Terry Rudolph.  Quantum coherence, time-translation symmetry, and thermodynamics.  _Phys. Rev. X_ , 5(2):021001, 2015.  doi: 10.1103/PhysRevX.5.021001.
  * [Scala et al.(2007)Scala, Militello, Messina, Piilo, and Maniscalco] M Scala, B Militello, A Messina, J Piilo, and S Maniscalco.  Microscopic derivation of the jaynes-cummings model with cavity losses.  _Phys. Rev. A_ , 75(1):013811, 2007.  doi: 10.1103/PhysRevA.75.013811.