Heat-Bath Algorithmic Cooling with optimal ther-
malization strategies
Álvaro M. Alhambra
1
, Matteo Lostaglio
2
, and Christopher Perry
3
1
Perimeter Institute for Theoretical Physics, Waterloo, ON N2L 2Y5, Canada
2
ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute of Science and Technology, Castelldefels (Barcelona), 08860, Spain
3
QMATH, Department of Mathematical Sciences, University of Copenhagen, Universitetsparken 5, 2100 Copenhagen, Denmark
September 19, 2019
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 orig-
inally used the thermal environment only to
fully thermalize ancillas at the environment tem-
perature. Here we extend HBAC protocols by
optimizing over the thermalization strategy. We
find, for any d-dimensional system in an arbi-
trary initial state, provably optimal cooling pro-
tocols 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 devia-
tions from the ideal scenario. For a single target
qubit, the optimal protocol can be well approxi-
mated with a Jaynes-Cummings interaction be-
tween the system and a single thermal bosonic
mode for a wide range of environmental temper-
atures. This admits an experimental implemen-
tation close to the setup of a micromaser, with a
performance competitive with leading proposals
in the literature. The proposed protocol pro-
vides 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 β-
permutations, whose access is sufficient to sim-
ulate a generic thermalization process, 3. show
how to use abstract toolbox developed within
the resource theory approach to thermodynam-
ics to perform challenging optimizations, while
combining it with open quantum system dynam-
ics tools to approximate optimal solutions within
physically realistic setups.
Cooling is a central problem in quantum physics and
in realizing technologies for quantum information pro-
cessing. The ability to produce a set of highly pure,
‘cold’, quantum states is vital for the construction of a
quantum computer [1]. 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 optomechan-
ics [2] to NMR [3, 4], ion traps [5] and superconducting
qubits [6].
Here we use powerful techniques, developed within
the resource theory approach to thermodynamics [7],
to greatly extend an important class of cooling al-
gorithms known as Heat-Bath Algorithmic Cooling
(HBAC) [3, 4]. The goal of these is to maximize the
purity of a target system S in a given number of cool-
ing rounds. Each round of the algorithm starts with a
unitary applied to the target together with several aux-
iliary 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 [8], whose asymptotic performance has
been recently derived for a single target qubit starting
in a maximally mixed state [9].
Here, we generalize HBAC in the following sense: in-
stead of using the thermal environment only to refresh
the auxiliary systems through a complete thermaliza-
tion, we allow strategies involving an incomplete ther-
malization of system and ancillas. In particular, we
optimize the protocol over every ‘dephasing thermaliza-
tion’, that is any quantum map on SA which 1. leaves
the thermal state on SA fixed and 2. dephases the in-
put state in the energy eigenbasis. This generalizes the
‘rethermalization’ used in previous protocols, and thus
Accepted in Quantum 2019-09-16, click title to verify 1
arXiv:1807.07974v3 [quant-ph] 18 Sep 2019
extends HBAC to a larger set that will be called Ex-
tended Heat-Bath Algorithmic Cooling (xHBAC). For
example, the recent protocol introduced in Ref. [10],
which proposes the use of the heat bath to implement
a non-local thermalization ‘state-reset’ (SR) process re-
lated to the Nuclear Overhauser Effect [11], can already
be seen as a (non-optimal) protocol within our extended
family of xHBAC.
For every finite dimensional target state in an arbi-
trary 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 opera-
tions, 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) interac-
tion to a single bosonic thermal mode, itself weakly cou-
pled 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 ex-
perimental 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. Further-
more, we will assume we can access a thermal environ-
ment at inverse temperature β = (kT )
1
, with k Boltz-
mann’s constant. We denote by H
S
=
P
d1
i=0
E
i
|iihi|,
E
0
··· E
d1
, the Hamiltonian of S and by ρ
(k)
S
the state after round k. Also, we denote by H
A
the
Hamiltonian of the auxiliary systems, initially in state
ρ
A
, and by τ
X
= e
βH
X
/Tr[e
βH
X
] the thermal state
on X = S, A. A protocol is made of k rounds, each
allowing for (Fig. 1):
1. Unitary. Any unitary U
(k)
applied to SA.
Figure 1: A generic xHBAC protocol.
2. Dephasing thermalization. Any Λ
(k)
applied
to SA, where Λ
(k)
is any quantum map such that
i) Λ
(k)
(τ
S
τ
A
) = τ
S
τ
A
(thermal fixed point);
ii) If |Ei, |E
0
i are states of distinct energy on SA,
hE|Λ
(k)
(ρ
SA
)|E
0
i = 0 (dephasing).
At the end of the round, a refreshing stage returns the
auxiliary systems A to their original state ρ
A
. We de-
note the set of protocols whose rounds have this general
form by P
ρ
A
. Typically ρ
A
= τ
A
, and this operation
is simply a specific kind of dephasing thermalization.
However, more generally, ρ
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 or-
der. First, if we skipped the dephasing thermaliza-
tions and ρ
A
= τ
A
, we would get back to a standard
HBAC scheme. xHBAC protocols include thermaliza-
tion of subsets of energy levels, as in the mentioned SR
protocol of Ref. [10]. But they are by no means lim-
ited to these strategies. Second, for every k we wish
to find a protocol maximizing the ground state popu-
lation p
(k)
0
= h0|ρ
(k)
S
|0i over all sequences of k-round
operations. Note that often the literature on cooling
is restricted to finding asymptotically optimal proto-
cols (e.g., the PPA [8]), 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’ [9], 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 maxi-
mally active states. Given a state ρ with Hamiltonian
H, the maximally active state ˆρ is formed by diagonal-
izing ρ in the energy eigenbasis and ordering the eigen-
values in increasing order with respect to energy. That
is, ˆρ is the most energetic state in the unitary orbit of
ρ.
Second, we define the thermal polytope in the follow-
ing way. Given an initial state ρ with populations in
the energy eigenbasis p, it is the set of populations
p
0
of Λ(ρ), with Λ an arbitrary dephasing thermaliza-
tion. Without loss of generality, we assume that ev-
ery energy eigenspace has been diagonalized by energy
preserving unitaries. The thermal polytope is convex
Accepted in Quantum 2019-09-16, click title to verify 2
and has a finite number of extremal points (Lemma 12,
[12]). It is particularly useful in thermodynamics to
classify out of equilibrium distributions according to
their β-order [13]. Given a vector p and a correspondent
Hamiltonian with energy levels E
i
, the βorder of p is
the permutation of π of 0, 1, 2, . . . such that the Gibbs-
rescaled population is sorted in non-increasing order:
p
π(0)
e
βE
π(0)
p
π(1)
e
βE
π(1)
. . . .
Crucially, for any dimension d, we provide an explicit
construction of a set of dephasing thermalizations, de-
noted Λ
(π)
and called β-permutations, able to map
any given p to every extremal point of its thermal poly-
tope (Methods, Sec. 2.1 and Fig. 2). π and α are two
vectors of integers, each ranging over all permutations
of {0, ..., (d1)(r 1)}, if r is the dimension of A. Their
significance is the following. If we are given an initial
state p with β-order π, the state q
α
:= Λ
(π)
(p) is the
‘optimal’ state among all those with β-order α. More
precisely, there is no other state in the thermal polytope
of p which has β-order α and can be transformed into
q
α
by dephasing thermalizations. This is in analogy
with the role of permutations under doubly stochastic
maps. In fact, in the limit β +, β-permutations
are permutations and the β-ordering is the standard
sorting.
A particularly important β-permutation, denoted by
β
opt
, is the one that maximizes the ground state popula-
tion of S among all dephasing thermalizations and, fur-
thermore, achieves the largest partial sums
P
l
i=0
p
(k)
i
,
l = 0, ..., d 1, p
(k)
i
= hi|ρ
(k)
S
|ii. β
opt
is constructed
as follows: β
opt
= Λ
(π)
with π an ordering of en-
ergy 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)
π(m)
e
βE
SA
π(m)
are sorted in non-increasing order in m;
and
α =
(0, r 1) , (0, r 2) , . . . , (0, 0) ,
(1, r 1) , (1, r 2) , . . . , (1, 0) , . . . , (d 1, 0)
.
(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 ρ with population p, the ther-
mal polytope of the maximally active state ˆρ contains the
thermal polytope of U ρU
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
d1
> E
0
and β > 0. Without loss generality (by
Figure 2: An optimal cooling round on a d = 3 system. Map
the initial state after round k1 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. β-permutations map
to the extremal points of these polytope (dotted red arrows),
and β
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.
an initial diagonalizing unitary) we can take the initial
state of S to be diagonal in the energy eigenbasis with
p
(0)
0
··· p
(0)
d1
. Then, for a given auxiliary state ρ
A
,
the optimal cooling protocol in P
ρ
A
is such that in each
round k:
1. The unitary mapping ρ
(k1)
S
ρ
A
to the correspond-
ing maximally active state is applied to SA.
2. β
opt
is applied to SA.
The optimal protocol achieves p
(k)
0
1 at least expo-
nentially 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 β-permutation is applied. For the
proof, see Methods (Sec. 2.4).
In analyzing the performance of the protocols P
ρ
A
it is important to keep in mind the cost of prepar-
ing and controlling the auxiliary systems A, especially
if ρ
A
6= τ
A
. Remarkably, however, the optimal xH-
BAC protocol that uses no auxiliary systems A still has
p
(k)
0
1 exponentially in k. Furthermore, such proto-
col has the further advantage that the cooling opera-
tions do not change with the round k. Also note that
the initial state-dependent unitary can always be re-
placed by a complete thermalization while maintaining
the same cooling scaling, so that no prior knowledge of
Accepted in Quantum 2019-09-16, click title to verify 3
ρ
(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 de-
phasing thermalizations between energies (d 1, d 2),
(d 2, d 3), . . . , (1, 0) is performed, each maximising
the net population transfer i 7→ i1 (these are β-swaps,
discussed in more detail in the next section). Every d1
repetitions, the performance takes an elegant form:
p
(k(d1))
0
= 1 e
(E
d1
E
0
)
1 p
(0)
0
. (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 P
) and
energy gap E. The only nontrivial β-permutation is
Λ
β
= Λ
(π)
with π = {1, 0} and α = {0, 1}, which
induces transition probabilities h0|Λ
β
(|1ih1|)|0i = 1,
h1|Λ
β
(|0ih0|)|1i = e
βE
. This is the β-swap introduced
in [12], which can be realized by the dephasing thermal-
ization
Λ
β
(ρ
S
) =σ
ρ
S
σ
+
+ e
βE
σ
+
ρ
S
σ
(3)
+ (1 e
βE
)σ
σ
+
ρ
S
σ
σ
+
where σ
+
= |1ih0|, σ
= σ
+
. Since the unitary map-
ping the state to its correspondent maximally active
state is the Pauli X unitary, Theorem 1 reads as fol-
lows:
Corollary 1. Assume E > 0, β > 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
p
(0)
1
. The
optimal cooling protocol in P
is such that in each round
k:
1. The Pauli X unitary is applied to S.
2. The β-swap Λ
β
is applied to S.
The population of the ground state after round k is:
p
(k)
0
= 1 e
E
1 p
(0)
0
, (4)
and p
(k)
0
1 as k .
The performance can be obtained by direct compu-
tation, 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 β-swaps cannot
be realized by Markovian interactions with the thermal
environment. In fact, if we optimized over Marko-
vian interactions only, the optimal protocol with no
auxiliary systems would be the trivial thermalization
at the environment temperature ρ
S
7→ τ
S
. Then,
at best we would achieve a ground state population
1/(1 + e
β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 in-
finite dimensional environment is necessary to increase
the purity of the target system to 1 in the absence of
initial system-environment correlations [14, 15]; how-
ever, we do not need a complex environment: it was
shown in [12] that a single bosonic mode in a thermal
state τ
B
=
1 e
βE
P
n=0
e
E
|nihn| suffices to im-
plement the required β-swap. Specifically, the unitary
needed is
U
β
SB
= |0, 0ih0, 0|+
X
n=1
(|0, nih1, n 1| + |1, n 1ih0, n|) .
(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 P
has the initial state ρ
(0)
S
τ
B
and is such that in each round k:
1. The Pauli X I unitary is applied to ρ
(k1)
SB
, the
state of system-bath after round k 1.
2. The unitary U
β
SB
is applied to X I(ρ
(k1)
SB
)X I.
1.3 Robustness to imperfections
Theorem 2 holds even beyond some of the (standard)
idealizations we made:
1. If the bosonic mode gap does not perfectly match
the system gap, U
β
SB
still realises the cooling of
Accepted in Quantum 2019-09-16, click title to verify 4
Eq. (4), with the caveat that some work flows at
each β-swap.
2. We can consider a typical imperfection, i.e. a small
anharmonicity in the ladder of B. In the Meth-
ods (Sec. 2.8) we show that taking E
n+1
E
n
=
E(1 (n + 1)τ
2
) + o(τ
2
), the protocol of Theo-
rem 2 performs very closely to the ideal scenario of
Eq. (4) even for moderate anharmonicity (within
0.005% for βE=1, τ = 0.05).
Furthermore, up to now we considered an idealized
scenario in which one can perform the required β-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 re-
quired by the β-swap, one can only realize dephas-
ing thermalizations Λ
(k)
with an induced de-excitation
probability h0|Λ
(k)
(|1ih1|)|0i 1 . Denote this set
by P
. Define an -noisy β-swap as the dephasing ther-
malization Λ
β
with h0|Λ
β
(|1ih1|)|0i = 1 (when = 0
this is the β-swap). Then one can derive the following
noise-robust version of Theorem 1:
Theorem 3. Under the assumptions of Corollary 1 and
given
1
1+e
βE
+e
2βE
, the optimal nontrivial cooling
protocol in P
is such that in each round k:
1. The Pauli X unitary is applied to S.
2. The -noisy β-swap Λ
β
is applied to S.
The population of the ground state after round k is:
p
(k)
0
=1
2 (1 )Z
(6)
((1 )Z 1)
k
1
2 (1 )Z
p
(0)
0
,
where Z = 1+e
βE
and p
(k)
0
1
2(1)Z
as k .
In words: even in the presence of (moderate) noise,
the optimal strategy is to perform the best approxi-
mation to the ideal protocol; hence, the simple struc-
ture 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 an-
other form of ‘robustness’: if we optimize over the larger
set of thermalization models known as thermal opera-
tions [16], which include protocols where some amount
of superposition in the energy basis survives the ther-
malization step, Theorem 3 holds unchanged. Surpris-
ingly, the extra coherent control in the thermalization
phase is not necessary for optimal cooling.
𝛽𝐸 = 0.5
Out[45]=
0
1
2
3
4
5
6
0.88
0.90
0.92
0.94
0.96
0.98
1.00
Number of rounds k
p
k
SWAP JC model lower
JC model upper SR
2
PPA
𝛽𝐸 = 0.2
𝛽𝐸 = 1
𝛽𝐸 = 2
Out[45]=
0
1
2
3
4
5
6
0.88
0.90
0.92
0.94
0.96
0.98
1.00
Number of rounds k
p
k
SWAP JC model lower
JC model upper SR
2
PPA
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 re-
alization of this protocol (red, see Sec. 3.1 for details on how
to calculate them), SRΓ protocol from [10] run with 1 ancilla
qubit (green), PPA protocol from [8] run with 2 ancilla qubits
(purple). The initial state for all protocols is thermal.
1.4 An experimental proposal
We can now provide an explicit experimental proposal.
More details about the calculations and simulations in-
volved can be found in Sec. 3.1 -3.2.
The unitary U
β
SB
of Eq. (5) can be realized exactly
with an intensity-dependent Jaynes-Cummings model
[17, 18],
˜
H
JC
= g(σ
+
(aa
)
1/2
a + σ
(aa
)
1/2
a
).
However, a perhaps more promising avenue is to ap-
proximate the β-swap steps with a resonant Jaynes-
Cummings (JC) coupling with a thermal bosonic mode,
H
JC
= g(σ
+
a + σ
a
). Assuming good control of
the interaction time s, we can numerically optimize s to
realize an -noisy β-swap with the small . If the ther-
mal 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 pro-
posals in the literature, for an initially thermal target.
The JC protocol outperforms the PPA with 2 ancillas
[8] (even with some limitations in the timing accuracy
and the total available waiting time), with the exclusion
of the very high temperature regime; when βE is not
too small, it performs comparably (if potentially slightly
worse) to the ideal version of the non-local thermaliza-
tion scheme with 1 ancilla proposed in [10]. The op-
timal β-swap protocol outperforms all these protocols,
but requires an intensity-dependent JC model whose
implementation we leave as an open question.
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
Accepted in Quantum 2019-09-16, click title to verify 5
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 β. If the interaction parameters
are appropriately chosen, such that a small is achieved, the outgoing qubit has a much lower temperature.
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 ini-
tially 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 fir-
ing rate, the cooling performance settles to a constant
as many atoms are cooled (see Fig. 5).
This is, to our knowledge, the most appealing set-
ting to experimentally implement the protocol and is,
in fact, highly reminiscent of a micromaser [19, 20]. 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 pro-
posed 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 cool-
ing protocol the aim is to obtain these very pure atoms
in the first place. Nevertheless, this suggests that exper-
imental 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 in-
cludes both cavity QED [19, 20] and solid-state settings
[21].
In conclusion, often HBAC techniques made the im-
plicit 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 es-
tablished limits, in particular allowing exponential con-
vergence to the ground state in the ideal scenario.
2 Methods
2.1 β-permutations and the thermal polytope
Here we discuss how to construct β-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
p and p
0
, p majorizes p
0
, denoted p p
0
, if
P
k
i=0
p
i
P
k
i=0
p
0
i
for k = 0, ..., d 2, where x
denotes the vector
x arranged in descending order. ρ ρ
0
is defined as the same relation among the corresponding eigenvalues. Then
define
Definition 1 (Thermo-majorization [13]). Given a state ρ with Hamiltonian H =
P
d1
i=0
E
i
|iihi|, let p =
(p
0
, . . . , p
d1
) where p
i
= hi|ρ|ii and E = (E
1
, . . . , E
d
). The thermo-majorization curve of ρ is formed by:
1. Applying a permutation π of {0, ..., d 1} to both p and E such that p
π(i)
e
βE
π(i)
is in non-increasing order.
We refer to the vector π with elements π(i) as the β-order of p.
2. Plotting the ordered ‘elbow’ points (0, 0),
n
P
k
i=0
e
βE
π(i)
,
P
k
i=0
p
π(i)
o
d1
k=0
and connecting them piecewise
Accepted in Quantum 2019-09-16, click title to verify 6
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 βE = 1, the coupling strength is g = 1 and time of
interaction t = 98.92.
linearly to form a concave curve - the thermo-majorization curve of p.
Given two probability distributions p and p
0
associated with the same energy levels, we say that p thermo-majorizes
p
0
if the thermo-majorization curve of p is never below that of p
0
. We denote this by p
th
p
0
. Also, if ρ and ρ
0
are states with population vectors p and p
0
, the notation ρ
th
ρ
0
denotes p
th
p
0
.
Note that
th
becomes in the infinite temperature limit β 0. One can show that p
th
p
0
is equivalent
to the existence of a Gibbs stochastic matrix (a stochastic matrix G with Gg = g, if g
i
e
βE
i
) such that
Gp = p
0
[22]. Furthermore, a dephasing thermalization Λ acts on the population vector as a Gibbs-stochastic matrix
G
j|i
= hj|Λ(|iihi|)|ji; conversely, every Gibbs-stochastic matrix can be realized by dephasing thermalizations. As
such, the thermal polytope of p coincides with the set of probability distributions p
0
such that p
th
p
0
. 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 p explicitly, we construct a set of maps, called β-permutations, that take p to each one of its
extremal points.
Algorithm to construct Λ
(π)
. By the discussion above each β-permutation Λ
(π)
is fully characterized by a
matrix of transition probabilities P
(π)
j|i
= hj|Λ
(π)
(|iihi|)|ji. Let π and α each be one of the d! 1 possible
β-orders. Order the rows of P
(π)
according to α and the columns according to π:
G =
P
(π)
α(0)π(0)
. . . P
(π)
α(0)π(d1)
.
.
.
.
.
.
P
(π)
α(d1)π(0)
. . . P
(π)
α(d1)π(d1)
.
Then G is constructed as follows.
Row 0 : If e
βE
α(0)
< e
βE
π(0)
, set:
k
0
= 1, G
00
=
e
βE
α(0)
e
βE
π(0)
, G
0j
= 0, j [1, d 1],
Accepted in Quantum 2019-09-16, click title to verify 7
else, let k
0
be the smallest integer such that:
k
0
X
j=0
e
βE
π(j)
e
βE
α(0)
.
and set:
G
0j
= 1, 0 j k
0
1
G
0k
0
=
e
βE
α(0)
P
k
0
1
j=0
e
βE
π(j)
e
βE
π(k
0
)
G
0j
= 0, j [k
0
+ 1, d 1]
Row m 1: If
P
m
i=0
e
βE
α(i)
<
P
k
m1
j=0
e
βE
π(j)
, set:
k
m
= k
m1
, G
mj
= 0, 0 j k
m1
1
G
mk
m1
=
e
βE
α(m)
e
βE
π
(
k
m1
)
, G
mj
= 0, j [k
m1
+ 1, d 1]
else, let k
m
be the smallest integer such that:
k
m
X
j=0
e
βE
π(j)
m
X
i=0
e
βE
α(i)
. (7)
and set:
G
mj
= 0, j [0, k
m1
1], G
mk
m1
= 1
m1
X
i=0
G
ik
m1
G
mj
= 1, j [k
m1
+ 1, k
m
1]
G
mk
m
=
P
m
i=0
e
βE
α(i)
P
k
m
1
j=0
e
βE
π(j)
e
βE
π(k
m
)
G
mj
= 0, j [k
m
+ 1, d 1].
In the next section, we prove the important fact that, if π is the β-order of p, the set {P
(π)
p}
α
, for α varying
over all possible β-orders, includes all extremal points of the thermal polytope of p (Lemma 2).
2.2 Proof of properties of β-permutations
In the following we show that the β-permutations, defined Section 1.1, convert p into any of the extremal points of
the thermal polytope p
α
, while being independent of the specific form of p (the set of β-permutation that need to
be applied only depends on the β-order of p). This is in complete analogy with permutations and their role within
the theory of doubly-stochastic matrices as determined by Birkhoff’s theorem [23] (in fact, β-permutations become
permutations in the infinite temperature limit).
Specifically, we will prove the following facts:
1. P
(π)
is a Gibbs-stochastic map.
2. Given p with β-order π, p
α
:= P
(π)
p has β-order α and is such that there is no q with β-order α such that
p
th
q
th
p
α
(Lemma 2).
3. The set of p
α
for varying α includes all the extremal points of the themal polytope of p (the latter is equivalently
defined as all q such that p
th
q).
Accepted in Quantum 2019-09-16, click title to verify 8
2.2.1 Proof of claim 1
Recall that we reordered the rows of P
(π)
according to α and the columns according to π:
G =
P
(π)
α(0)π(0)
. . . P
(π)
α(0)π(d1)
.
.
.
.
.
.
P
(π)
α(d1)π(0)
. . . P
(π)
α(d1)π(d1)
.
Hence, the condition of Gibbs-stochasticity can be rewritten as
G
ij
0 (8)
d1
X
i=0
G
ij
= 1 (9)
d1
X
j=0
G
ij
e
βE
π(j)
= e
βE
α(i)
. (10)
For simplicity, we repeat here the algorithm for constructing P
(π)
. Row 0 is populated as follows:
if e
βE
α(0)
< e
βE
π(0)
then
Set:
k
0
= 1 (11)
G
00
=
e
βE
α(0)
e
βE
π(0)
(12)
G
0j
= 0, 1 j d 1 (13)
else
Let k
0
be the smallest integer such that:
k
0
X
j=0
e
βE
π(j)
e
βE
α(0)
. (14)
Set:
G
0j
= 1, 0 j k
0
1 (15)
G
0k
0
=
e
βE
α(0)
P
k
0
1
j=0
e
βE
π(j)
e
βE
π(k
0
)
(16)
G
0j
= 0, k
0
+ 1 j d 1 (17)
end if
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
0. Gibbs preservation follows as:
d1
X
j=0
G
0j
e
βE
π(j)
= G
00
e
βE
π(0)
= e
βE
α(0)
.
For the second clause, it is again true that we have 0 G
0j
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:
d1
X
j=0
G
0j
e
βE
π(j)
=
k
0
1
X
j=0
e
βE
π(j)
+ e
βE
α(0)
k
0
1
X
j=0
e
βE
π(j)
= e
βE
α(0)
so Gibbs preservation is satisfied.
With row zero in place, row m (m {1, . . . d 1}) is populated as follows:
Accepted in Quantum 2019-09-16, click title to verify 9
𝐸
2
𝐸
2
𝐸
0
𝐸
0
𝐸
1
𝐸
1
Figure 6: Here we illustrate the notion of a maximal β-order state. The state p, shown in blue, has β-order π = (0, 1, 2). The
state p
α
for α = (2, 0, 1) is shown in red and has a simple geometrical interpretation. The choice of α (Condition 1 in Sec. 2.2.2)
fixes the x-axis points to be x
0
= e
βE
2
, x
1
= e
βE
2
+ e
βE
0
, x
2
= e
βE
2
+ e
βE
0
+ e
βE
1
. Then Condition 2 is equivalent to
the request that the curve in red touches the blue curve at these points.
if
P
m
i=0
e
βE
α(i)
<
P
k
m1
j=0
e
βE
π(j)
then
Set:
k
m
= k
m1
(18)
G
mj
= 0, 0 j k
m1
1 (19)
G
mk
m1
=
e
βE
α(m)
e
βE
π
(
k
m1
)
(20)
G
mj
= 0, k
m1
+ 1 j d 1 (21)
else
Let k
m
be the smallest integer such that:
k
m
X
j=0
e
βE
π(j)
m
X
i=0
e
βE
α(i)
(22)
Set:
G
mj
= 0, 0 j k
m1
1 (23)
G
mk
m1
= 1
m1
X
i=0
G
ik
m1
(24)
G
mj
= 1, k
m1
+ 1 j k
m
1 (25)
G
mk
m
=
P
m
i=0
e
βE
α(i)
P
k
m
1
j=0
e
βE
π(j)
e
βE
π(k
m
)
(26)
G
mj
= 0, k
m
+ 1 j 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
Accepted in Quantum 2019-09-16, click title to verify 10
first clause 0 G
mk
m1
1 follows from:
e
βE
α(m)
k
m1
X
i=0
e
βE
π(i)
m1
X
i=0
e
βE
α(i)
= e
βE
π(k
m1
)
+
k
m1
1
X
i=0
e
βE
π(i)
m1
X
i=0
e
βE
α(i)
e
βE
π(k
m1
)
,
where in the last inequality we used the definition of k
m1
. Finally, note that:
d1
X
j=0
G
mj
e
βE
π(j)
= G
mk
m1
e
βE
π
(
k
m1
)
= e
βE
α(m)
and hence Eq. (10) holds.
For the second clause, the definitions of k
m1
and k
m
ensure that both 0 G
mk
m1
1 and 0 G
mk
m
1
hold. We now show Gibbs preservation. First note that:
d1
X
j=0
G
mj
e
βE
π(j)
= G
mk
m1
e
βE
π
(
k
m1
)
+
k
m
1
X
j=k
m1
+1
e
βE
π(j)
+ G
mk
m
e
βE
π(k
m
)
=
1
m1
X
i=0
G
ik
m1
!
e
βE
π
(
k
m1
)
+
k
m
1
X
j=k
m1
+1
e
βE
π(j)
+
m
X
i=0
e
βE
α(i)
k
m
1
X
j=0
e
βE
π(j)
(28)
where we have used Eq. (24) and Eq. (26) in the last line. Consider now
P
m1
i=0
G
ik
m1
. Note that there exists an
integer r such that 0 r m 1 and:
G
ik
m1
= 0, i < r
G
rk
m1
=
P
r
i=0
e
βE
α(i)
P
k
r
1
j=0
e
βE
π(j)
e
βE
π(k
r
)
G
ik
m1
=
e
βE
α(i)
e
βE
π
(
k
i
)
, r + 1 i m 1
with k
i
= k
r
for r i m 1. Using these expressions in Eq. (28), we see that:
d1
X
j=0
G
mj
e
βE
π(j)
=
1
P
m1
i=0
e
βE
α(i)
P
k
m1
1
j=0
e
βE
π(j)
e
βE
π
(
k
m1
)
!
e
βE
π
(
k
m1
)
+
k
m
1
X
j=k
m1
+1
e
βE
π(j)
+
m
X
i=0
e
βE
α(i)
k
m
1
X
j=0
e
βE
π(j)
=
k
m
1
X
j=0
e
βE
π(j)
+ e
βE
α(m)
k
m
1
X
j=0
e
βE
π(j)
= e
βE
α(m)
and hence the Gibbs distribution is preserved.
Finally, let us show that the matrix constructed is stochastic. Using the expressions for G
ik
m1
given above, the
fact that k
i
= k
r
for all r i m 1 and the definition of k
m1
,
m1
X
i=0
G
ik
m1
=
P
m1
i=0
e
βE
α(i)
P
k
m1
1
i=0
e
βE
π(i)
e
βE
π(k
m1
)
< 1.
Together with Eq. (24) and the already proved fact that G
ij
[0, 1] for all i, j {0, . . . , d 1}, we get that G is
stochastic.
Accepted in Quantum 2019-09-16, click title to verify 11
2.2.2 Proof of claim 2
This can be stated as the following Lemma:
Lemma 2. If p has β-order π and p
α
:= P
(π)
p, then p
α
has β-order α and is ‘maximal’ in the sense that there
is no q with β-order α such that p
th
q
th
p
α
.
We now construct the proof. β-permutations have a simple geometrical description in terms of thermo-
majorization curves. In row zero, we want to maximize the population of E
α(0)
subject to the thermo-majorization
constraints. To do this, we compare e
βE
α(0)
and e
βE
π(0)
. If e
βE
α(0)
< e
βE
π(0)
, then we cannot move all of the
population in E
π(0)
to E
α(0)
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
βE
α(0)
e
βE
π(0)
, then we can move
all of the population in E
π(0)
to E
α(0)
and set G
00
= 1. We then try to move population from E
π(1)
to E
α(0)
and
repeat this process until we reach an energy level whose population we cannot move into E
α(0)
without violating
thermo-majorization. This energy level is defined through Eq. (14) and given by E
π(k
0
)
. From E
π(k
0
)
we can only
move a fraction of the population into E
α(0)
, given by Eq. (16). When we reach this point, we have moved as much
population as possible into E
α(0)
. This case is illustrated in Fig. 7b.
In the mth row we are determining the population of E
α(m)
after the transformation. At this point in the
construction, all of the population from energy level E
π(j)
for 0 j k
m1
1 has been transferred already into
the set of energies
E
α(i)
m1
i=0
. We thus have G
mj
= 0 for 0 j k
m1
1. We now check to see how much of
the remaining population in E
π(k
m1
)
we can move to E
α(m)
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
m1
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
m1
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
π(k
m1
)
to E
α(m)
until we reach an energy level E
π(k
m
)
,
where this is not possible due to thermo-majorization. This is illustrated in Fig. 8b.
Let c
p
: [0, Z
S
] [0, 1] be the function such that c
p
(x) is the height of the thermo-majorisation curve of p at x.
The action of the β-permutation matrix described above is such that the thermo-majorization curve of p
α
= P
(π)
p
is constructed as follows. For i {0, . . . , d 1}, denoting by (x
α
i
, y
α
i
) the points that are piecewise linearly connected
to give the thermo-majorization curve of p
α
, one has:
1. Let x
α
i
=
P
i
j=0
e
βE
α
1
(j)
and y
α
i
= c
p
(x
α
i
).
2. Define p
α
i
:= y
α
α(i)
y
α
α(i1)
, with y
α(1)
:= 0.
By construction, there is no q with β-order α such that p
th
q
th
p
α
.
2.2.3 Proof of claim 3
That all extremal points of the thermal polytope associated to p have the above form is stated and proved in
Lemma 12 of Ref. [12].
2.3 Proof of Lemma 1: optimal unitary control
First we recall the notion of maximally active state associated to a given ρ, defined as:
Definition 2 (Maximally active state). Let ρ be the state of a system with associated Hamiltonian
H =
P
d1
i=0
E
i
|iihi|, where E
i
E
i+1
. The maximally active state associated to ρ is
ˆρ =
d1
X
i=0
λ
i
|iihi|,
where
n
λ
i
o
d1
i=0
are the eigenvalues of ρ arranged in ascending order.
This corresponds to the state with the highest energy along the whole unitary orbit of ρ, i.e. ˆρ =
arg max
U
Tr[UρU
H]. We can now prove Lemma 1 (note that we use facts from Section 2.1 in this proof):
Accepted in Quantum 2019-09-16, click title to verify 12
(a) Row 0, case 1
(b) Row 0, case 2
Figure 7: Thermo-majorization diagrams illustrating the construction of row zero of P
(π,α)
. In case 1, e
βE
α(0)
< e
βE
π(0)
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
(π,α)
. In case 1,
P
m
i=0
e
βE
α(i)
<
P
k
m1
j=0
e
βE
π(j)
holds, while in case 2 it does not.
Proof. The claim is equivalent to proving ˆρ
th
UρU
for every unitary U. We will first show that for any
permutation π of {1, ..., d}:
λ
th
πλ
, (29)
where λ
denotes the vector of eigenvalues of ˆρ arranged in ascending order. To show it, we will construct a sequence
of d 1 Gibbs-stochastic matrices converting λ
into πλ
, passing through the intermediate states q
(1)
, . . . , q
(d1)
,
with q
(d1)
= πλ
and q
(0)
= λ
.
We start with q
(0)
. We wish to move the probability λ
0
associated with energy level E
0
to energy level E
π(0)
.
To do this, we transpose in sequence the populations of energy levels 0 1, 1 2, . . . , π(0) 1 π(0). Each of
these is possible under Gibbs-stochastic matrices, because λ
0
λ
j+1
, j and E
j
E
j+1
. In fact, they are achieved
with the matrices:
G
(j,j+1)
=
1 µe
β(E
j+1
E
j
)
µ
µe
β(E
j+1
E
j
)
1 µ
I
\(j,j+1)
,
where I
\(j,j+1)
denotes the identity on all levels different from (j, j + 1) and
µ =
λ
j+1
λ
0
λ
j+1
λ
0
e
β(E
j+1
E
j
)
0.
Accepted in Quantum 2019-09-16, click title to verify 13
Hence it is always possible to move population λ
0
to energy level E
π(0)
by setting
q
(1)
= G
(π(0)1(0))
. . . G
(1,2)
G
(0,1)
q
(0)
.
q
(1)
coincides with πλ
on element π(0). Now, if we truncate element π(0) from q
(1)
, we obtain a vector satisfying
q
(1)
0
= λ
1
q
(1)
1
= λ
2
. . .
q
(1)
π(0)1
q
(1)
π(0)+1
··· q
(1)
d1
= λ
d1
.
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 λ
1
to the energy level E
π(1)
. Applied to q
(1)
, this
sequence give a state q
(2)
which coincides with πλ
in elements π(0), π(1). It should be clear that this construction
can be repeated on every intermediate q
(m)
, each time applying it to the distribution in which we ignore the energy
levels E
π(0)
, ..., E
π(m1)
, since these have the correct occupation probability. This provides a sequence of states
culminating in q
(d1)
= πλ
as required.
Having shown that Eq. (29) holds, we now argue for ˆρ
th
UρU
for all unitaries U . Given U , let
˜
p = diag
UρU
(remembering to first diagonalize the degenerate energy subspaces if necessary) and we want to prove λ
th
˜
p.
Let
˜
p
be the state formed by arranging the elements of
˜
p in ascending order. By Eq. (29), we have that
˜
p
th
˜
p.
In addition, by the Schur-Horn theorem, we have that λ
˜
p
. Combining this with the fact that both λ
and
˜
p
are ordered in terms of increasing occupation probability and have the same β-order, it follows that:
1. If the thermo-majorization curve of λ
is constructed from points {(x
m
, y
m
)}
d1
m=0
and that of
˜
p
from points
{(x
0
m
, y
0
m
)}
d1
m=0
, we have x
m
= x
0
m
for each m.
2. y
m
y
0
m
for each m, since
P
m
i=0
λ
i
P
m
i=0
˜p
i
for all m = 1, . . . , d 1.
Hence λ
th
˜
p
th
˜
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
ρ
(k1)
S
ρ
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
π
k
is the β-order of the state after the unitary, this can be chosen to be the β-permutation Λ
(π
k
)
with α 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
P
l
i=0
p
(k)
i
, l {0, . . . , d 1} where
n
p
(k)
i
o
d1
i=0
are the
populations of ρ
(k)
S
arranged in descending order. This follows from Lemma 2 and the fact that the β-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 P
ρ
A
. Suppose that
at the beginning of round k we have one of two diagonal states ρ
(k1)
S
and ˜ρ
(k1)
S
such that ρ
(k1)
S
˜ρ
(k1)
S
. ρ
(k)
S
will represent the trajectory followed by the state when we apply the claimed optimal protocol, whereas ˜ρ
(k)
S
will
be the trajectory followed by a generic protocol (hence, ρ
(0)
S
= ˜ρ
(0)
S
). Our goal is to show that:
ρ
(k)
S
= Tr
A
h
Λ
(π
k
)
U
(k)
m.a.
ρ
(k1)
S
ρ
A
i
˜ρ
(k)
S
= Tr
A
h
Λ
(k)
V
(k)
˜ρ
(k1)
S
ρ
A
i
(30)
for all choices of ρ
A
, V
(k)
and Λ
(k)
. Here U
(k)
m.a.
denotes the unitary creating the maximally active state on SA, V
(k)
is an arbitrary unitary and Λ
(k)
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 ρ
(k1)
S
˜ρ
(k1)
S
, then:
ρ
(k1)
S
ρ
A
˜ρ
(k1)
S
ρ
A
, ρ
A
. (31)
Accepted in Quantum 2019-09-16, click title to verify 14
This gives:
U
(k)
m.a.
ρ
(k1)
S
σ
A
th
U
(k)
m.a.
˜ρ
(k1)
S
σ
A
th
V
(k)
˜ρ
(k1)
S
σ
A
, ∀V
(k)
. (32)
The first line follows from Eq. (31) by comparison of the thermomajorization curves, since maximally active states
have the same β-orders. The second line follows from Lemma 1. Now, by definition of Λ
(π
k
)
and Eq. (32),
Tr
A
h
Λ
(π
k
)
U
(k)
m.a.
ρ
(k1)
S
σ
A
i
Tr
A
h
Λ
(k)
V
(k)
˜ρ
(k1)
S
σ
A
i
.
2.5 A qudit protocol with β-swaps and no ancillas
The β-permutations affecting only two energy levels (i, j) at once are called β-swaps [12] and have the form
β
i,j
=
1 e
β(E
j
E
i
)
1
e
β(E
j
E
i
)
0
I
\(i,j)
. (33)
where I
\(i,j)
is the identity on every level l 6= i, j.
We now define the following Gibbs-stochastic matrix, which is a many-level generalization of the matrix A
6
from
Ref. [24]:
˜
G =
1 e
β
01
1 e
β
12
... 1 e
β
d1 d2
1
e
β
01
0 ... 0 0
0 e
β
12
.
.
.
0 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 0 ... 0 0
0 0 ... e
β
d1 d2
0
.
where
ij
= E
i
E
j
.
˜
G can be generated as a sequence of β-swaps:
˜
G =
Q
d2
i=0
β
i+1,i
. Now consider the protocol
in which, at each round, we perform the following steps:
1. A unitary transformation U is performed, that swaps the population of the ground state and the most excited
state.
2.
˜
G is performed.
Note that U induces the permutation S
(0,d1)
that flips the states 0 d 1. Next, define the resulting cooling
stochastic matrix C =
˜
GS
(0,d1)
. Also denote by Ω =
P
d2
i=0
(E
i+1
E
i
) = E
d1
E
0
. Then one can verify:
C
d1
=
1 1 e
β
1 e
β
. . . 1 e
β
0 e
β
0 . . . 0
0 0 e
β
.
.
.
0
.
.
.
.
.
. 0
.
.
.
.
.
.
0 0 . . . . . . e
β
.
Let p
(k)
i
be occupations after k steps of the protocol, and R
(k)
= 1 p
(k)
0
. One has
p
(k(d1))
0
= 1 R
(k(d1))
, R
(k(d1))
= e
βk
S
(1)
.
By assumption > 0 and β > 0, and Eq. (2) follows.
Accepted in Quantum 2019-09-16, click title to verify 15
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 p as Gibbs-
stochastic matrices, in this case 2 ×2. Using the results of Ref. [25], a direct computation shows that the action on
p of Markovian dephasing thermalizations can be written as
G
Mark
= (1 λ)I + λβ
01
, 0 λ 1/(1 + e
βE
) (34)
with β
01
given by Eq. (33).
Let (p, 1p), (q, 1 q) (s, 1s) 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 s and q by
s =
1 λe
βE
q + λ(1 q). (35)
Since p and q are related by a unitary, we have that p(1 p) q(1 q). Assume that the system is initially hotter
than the bath, p < 1/(1 + e
βE
). By direct inspection the maximum is found at λ = 1/(1 + e
βE
), where one
achieves s = 1/(1 + e
βE
), the thermal ground state population at temperature β. The same would be true if one
optimized over all Markovian thermal operations (using again the results of Ref. [25] 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 ρ
(k)
S
at the beginning of round k + 1 of the optimal protocol is incoherent for every
k 1. As such, the system-bath state after round k will have the general form:
ρ
(k)
SB
=
X
n=0
p
(k)
0,n
|0, nih0, n| + p
(k)
1,n
|1, nih1, n|
,
for some occupation probabilities p
(k)
i,n
with i {0, 1} and n {0, 1, . . . }. Applying X to S, followed by U
β
SB
gives
ρ
(k)
SB
7→ ρ
(k+1)
SB
= p
(k)
1,0
|0, 0ih0, 0| + p
(k)
0,0
|0, 1ih0, 1|
+
X
n=1
p
(k)
1,n
|1, n 1ih1, n 1| + p
(k)
0,n
|0, n + 1ih0, n + 1|
.
The relation between the occupations at step k and those at step k + 1 is (see Fig. 9)
p
(k+1)
0,0
= p
(k)
1,0
,
p
(k+1)
0,n
= p
(k)
0,n1
, n 1,
p
(k+1)
1,n
= p
(k)
1,n+1
, n 0.
Solving these equations recursively, we find the populations at step k as a function of the initial occupation
probabilities (at k = 0):
p
(k)
1,n
= p
(0)
1,n+k
p
(k)
0,n
=
(
p
(0)
0,nk
, n k,
p
(0)
1,kn+1
, n < k.
As before, we denote by p
(k)
0
the ground state population of the system after k steps of the protocol. We have
p
(k)
0
= p
(k)
0,0
+
X
n=1
p
(k)
0,n
= p
(k1)
1,0
+
X
n=1
p
(k1)
0,n1
= p
(k1)
1,0
+ p
(k1)
0
.
Accepted in Quantum 2019-09-16, click title to verify 16
Figure 9: The circulation of populations induced by each round (Pauli X followed by U
β
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.
We now use the relation
p
(k1)
1,0
= p
(k2)
1,1
= ··· = p
(0)
1,k1
= p
(0)
t
(0)
k
,
where t
(0)
k
=
1 e
βE
e
E
denotes the occupation of the kth energy level of the bosonic mode at the beginning
of the protocol. This finally leads to
p
(k)
0
= p
(0)
0
+
1 p
(0)
k1
X
n=0
t
(0)
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
an
n
)
E
an
n+1
E
an
n
= E(1 (n + 1)τ
2
) + o(τ
2
), (37)
where τ 1. To analize its effect, notice that the ground state from Eq. (36) is
p
(k)
0
= p
(0)
0
+
1 p
(0)
k1
X
n=0
t
(0)
n
, (38)
and that this does not rely on any particular form for the populations t
(0)
k
. Thus, for any τ, we have to compare
the two following sums,
k1
X
n=0
t
(0),an
n
=
P
k1
n=0
e
βE
an
n
P
n=0
e
βE
an
n
k1
X
n=0
t
(0)
n
=
P
k1
n=0
e
βE
n
P
n=0
e
βE
n
. (39)
A significantly large value τ = 0.05 gives a distortion E
an
n
/E
n
of 5% at n = 8. Also,
P
k1
n=0
t
(0),an
n
P
k1
n=0
t
(0)
n
approaches 1 as
k grows and the largest deviation from
P
k1
n=0
e
βE
n
P
n=0
e
β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) T on a qubit, which strictly include all dephasing thermalizations.
2.9.1 Preliminaries
First we review the definition of thermal operations [16, 26], which are channels acting on a quantum system S in
state ρ
S
and with Hamiltonian H
S
as:
T (ρ
S
) = Tr
B
U
ρ
S
e
βH
B
Tr [e
βH
B
]
U
, (40)
Accepted in Quantum 2019-09-16, click title to verify 17
where β = 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 [U, H
S
+ H
B
] = 0.
Let S be a qubit with Hamiltonian H
S
= E|1ih1|. Given ρ
S
, in addition to the occupation probabili-
ties p
S
= (p, 1 p)
T
, let ρ
01
= h0|ρ
S
|1i denote the off-diagonal term. Similarly, for a second state σ
S
, let
q
S
= (h0|σ
S
|0i, h1|σ
S
|1i)
T
and σ
01
= h0|σ
S
|1i. Then any thermal operation on a qubit can be described by:
1. The action on the population (given by a Gibbs-stochastic matrix G).
2. The action on the off-diagonal term.
Taken together this gives:
T (ρ
S
) = σ
S
q
S
= Gp
S
=
1 λe
βE
λ
e
βE
λ 1 λ
!
p
S
, λ [0, 1]
|σ
01
| = c|ρ
01
|, 0 c
p
(1 e
βE
λ)(1 λ)
, (41)
where the expression for c follows from the constraint |c|
G
00
G
11
[27, 28] which is a consequence of the complete
positivity of T . One can also apply a phase e
to ρ
01
by a unitary, so that in general c 6∈ R, but since this will be
irrelevant in the following considerations, without loss of generality we set θ = 0 and take c 0 (one can always
reversibly transform to any θ 6= 0 by an energy-preserving unitary, which is a thermal operation). Hence, for our
purposes a qubit thermal operation T is defined by the two parameters:
v [T ] = (λ, c) . (42)
With this notation, the β-swap introduced in the main text is a thermal operation T
β
that removes all quantum
coherences and has G = G
β-swap
where
G
β-swap
=
1 e
βE
1
e
βE
0
(43)
corresponds to v [T
β
] = (1, 0). 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 -noisy thermal operations as the set of thermal Operation T
such that
v[T
] = (λ, c), λ 1 . (44)
Let T
denote the set of cooling protocols using no ancilla in which at each round k
1. A unitary U
(k)
is applied to S,
2. An -noisy thermal operation T
(k)
is applied to S.
3. A unitary V
(k)
is applied to S
Then,
Theorem 4. Under the assumptions of Corollary 1 and given
1
1+e
βE
+e
2βE
, the optimal nontrivial cooling
protocol in T
is such that in each round k:
1. The Pauli X unitary is applied to S.
2. Any -noisy thermal operation T
with v(T
) = (1 , c) is applied to S.
The population of the ground state after round k is:
p
(k)
0
= 1
2 (1 )Z
((1 )Z 1)
k
1
2 (1 )Z
p
(0)
0
,
where Z = 1 + e
βE
and p
(k)
0
1
2(1)Z
as k .
Accepted in Quantum 2019-09-16, click title to verify 18
Before we prove this theorem, let us discuss the consequences. Note that a particular choice for the optimal
cooling protocol is T
with v(T
) = (1 , 0). This is just the -noisy β-swap defined in the main text. Since
T
P
, this means that the above theorem implies Theorem 3 of the main text as an immediate corollary. Also
note, that setting = 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 -noisy thermal operation: it could even
vary from round to round. The Jaynes-Cummings implementation described in the main text is an -noisy thermal
operation. On the population it acts as an -noisy β-swap, but it has c 6= 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 small enough.
Proof. The most general protocol can be written as
ρ
(k)
S
= V
(k)
T
(k)
U
(k)
. . . V
(1)
T
(1)
U
(1)
U
di
(ρ
S
), (45)
where U
(i)
(·) = U
(i)
(·)U
(i)
and V
(i)
(·) = V
(i)
(·)V
(i)
denote general unitaries applied to S in round i and T
(i)
is a
thermal operation on S applied at round i. Here U
di
= U
di
(·)U
di
is the unitary that transforms ρ
S
into a diagonal
form with p
(0)
0
p
(0)
1
. This can be done without loss of generality, since any protocol in P starts with an arbitrary
unitary U
(1)
.
We want to maximize the ground state population over all choices of U
(i)
, T
(i)
and 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:
ρ
(k)
S
=
p
(k)
0
0 1 p
(k)
and that p
(k)
1
2
(here we drop the subscript 0 to simplify the notation). In round k + 1 of the protocol, we have:
ρ
(k)
S
=
p
(k)
0
0 1 p
(k)
U
(k)
ρ
0(k)
S
=
q
(k)
a
(k)
a
(k)
1 q
(k)
T
(k)
ρ
00(k)
S
=
s
(k)
b
(k)
b
(k)
1 s
(k)
V
(k)
ρ
(k+1)
S
=
p
(k+1)
0
0 1 p
(k+1)
and our goal is to maximize p
(k+1)
.
As ρ
(k+1)
S
and ρ
00(k)
S
are related by a unitary, maximizing p
(k+1)
corresponds to minimizing the determinant of
ρ
00(k)
S
. Among all thermal operations associated to a fixed Gibbs-stochastic matrix G with Gq
(k)
:= s
(k)
(with
q
(k)
and s
(k)
defined through the above matrices), optimality of the protocol imposes that we choose 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 λ, i.e., from
Eq. (41), c = (1 λe
βE
)(1 λ) or, with the parametrization of Eq. (42), v[T
(k)
] = (λ, (1 λe
βE
)(1 λ)).
What remains to be done is to perform an optimization over all possible Gibbs-stochastic matrices G, parametrized
by λ [0, λ
max
] as in Eq. (41), where λ
max
1
1
1+e
βE
+e
2βE
(which corresponds to
1
1+e
βE
+e
2βE
). For each G,
the relation between s
(k)
, b
(k)
and q
(k)
and a
(k)
is
s
(k)
=
1 λe
βE
q
(k)
+ λ
1 q
(k)
|b
(k)
|
2
= |a
(k)
|
2
1 λe
βE
(1 λ) .
Finally, using the unitarity of U
(k)
, we can relate q
(k)
and a
(k)
to p
(k)
via:
p
(k)
1 p
(k)
= q
(k)
1 q
(k)
a
(k)
2
.
Accepted in Quantum 2019-09-16, click title to verify 19
The determinant of ρ
00(k)
S
is thus given by
f
p
(k)
,E
q
(k)
, λ
= (λ 1)
e
βE
λ 1
p
(k)
q
(k)

p
(k)
+ q
(k)
1
q
(k)
e
βE
λ + λ 1
λ
q
(k)
e
βE
λ + λ 1
λ + 1
, (46)
for fixed E and p
(k)
. Note that the equation is quadratic in both λ and q
(k)
. The coefficient of λ
2
is

1
q
(k)
2
q
(k)
e
βE
2
e
βE
h
q
(k)
1 q
(k)
p
(k)
1 p
(k)
i
,
which is clearly negative as q
(k)
1 q
(k)
= p
(k)
1 p
(k)
+ |a
(k)
|
2
, and thus the minimum values will be obtained
at either λ = 0 or λ = λ
max
. The case λ = 0 corresponds to not implementing the TO and leads to p
(k+1)
= p
(k)
.
Taking λ = λ
max
, we still need to show that the Pauli X unitary at each step is optimal.
Let us rewrite the function f
p
(k)
,E
q
(k)
, λ
max
as a polynomial in q
(k)
,
f
p
(k)
,E
q
(k)
, λ
max
f
(2)
p
(k)
,E
(λ
max
) (q
(k)
)
2
+ f
(1)
p
(k)
,E
(λ
max
) q
(k)
+ f
(0)
p
(k)
,E
(λ
max
) .
Given that it is a quadratic equation, the location of its minimum depends on the sign of the coefficient
f
(2)
p
(k)
,E
(λ
max
) = λ
max
(1 + e
βE
λ
max
(1 + e
βE
+ e
2β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
(k)
,E
(λ
max
) < 0. This is equivalent to
λ
max
>
1 + e
βE
1 + e
βE
+ e
2βE
= 1
1
1 + e
βE
+ e
2βE
, (47)
which is true by assumption. On top of this, we find that
f
p
(k)
,E
p
(k)
, λ
max
f
p
(k)
,E
1 p
(k)
, λ
max
= λ
max
(1 e
βE
)(λ
max
(1 + e
βE
) 1)(2p
(k)
1),
so the minimum is at q
(k)
= min
p
(k)
, 1 p
(k)
= 1 p
(k)
as long as λ
max
>
1
1+e
βE
(which is implied by Eq.
(47)).
Thus, the optimal protocol in T
is one in which, for every k, T
(k)
is a thermal operation (λ
max
, c
k
) with
λ
max
= 1 (the “best approximation” of the β-swap) and U
(k)
(·) = X(·) := X(·)X
. Note that we do not specify
how each 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
= λ
max
, without any control required on c
k
. The optimal ground
state population achieved by the above protocol satisfies
p
(k+1)
=
1 λ
max
e
βE
1 p
(k)
+ λ
max
p
(k)
,
as one can verify by a direct computation. Solving this recursion relation gives Eq. (6). One recovers the scaling of
Theorem 1 when λ
max
= 1. Furthermore, since λ
max
Z 1 e
βE
one has exponential convergence to 1
1λ
max
12λ
max
Z
.
Note that the trivial protocol in which we do not do a thermal operation is optimal only when p
(0)
1
1λ
max
12λ
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 ˜p
(k)
be
two ground state populations with p
(k)
˜p
(k)
. One can compute
p
(k+1)
˜p
(k+1)
(λ
max
Z 1)(p
(k)
˜p
(k)
) 0,
since λ
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).
Accepted in Quantum 2019-09-16, click title to verify 20
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 β-swap
against the ideal unitary U
β
SB
of Eq. (5). Consider a resonant JC Hamiltonian in rotating wave approximation
H
JC
= g(σ
+
a + σ
a
), (48)
coupling S with a single-mode bosonic bath prepared in a thermal state. Here a
and a are creation and annihilation
operators on B and σ
+
= |1ih0| and σ
= |0ih1|. The de-excitation probability is
G
0|1
(s) =
1
Z
B
X
n=1
sin
2
(s
n)e
βE(n1)
, (49)
where Z
B
= (1 e
β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. A Pauli X operation is performed on S.
2. The interaction Hamiltonian H
JC
operates for a time s that maximizes Eq. (49).
3. The bosonic mode is reset to a thermal state.
In Ref. [12], temperature-dependent upper bounds for the maximum achievable transition probability G
max
(s)
were derived as a function of
¯
β = βE:
G
max
(
¯
β)
(
1
16
8e
¯
β
e
2
¯
β
+ e
3
¯
β
+ 8
, for
¯
β [0,
log(4)
3
],
e
4
¯
β
e
3
¯
β
+ 1, for
¯
β
log(4)
3
.
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
()
0
1
e
¯
β
+
16e
2
¯
β
16e
¯
β
+e
3β
8
+1
, for
¯
β [0,
log(4)
3
],
1
e
¯
β
e
4
¯
β
+1
+1
, for
¯
β
log(4)
3
.
(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 [0, 5 × 10
3
], and obtain the curves in Fig. 3. As an illustrative example,
for
¯
β = 1 one has that the optimal JC asymptotic cooling is p
()
0
[0.9401, 0.9534]. As we know from Theorem 3,
the convergence is exponential and can be computed explicitly from Eq. (6) with = 1 G
0|1
(˜s) for the chosen
interaction time ˜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. Limited waiting time in the cavity, i.e. a bound on the maximum available s;
2. Limited timing accuracy, i.e. the achieved s fluctuates around a target value.
As a case study, set
¯
β = 1 and limit s to s = gt 10. Then the best available approximation to the β-swap is
realized for s 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 ±0.1, ±0.2 and ±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 ±0.2.
Accepted in Quantum 2019-09-16, click title to verify 21
Figure 10: In the time-limited JC model (blue curve on top) we optimise the interaction time s over the limited period s [0, 10].
To this we add increasing errors in the accuracy with which s is achieved: ±0.1, ±0.2 and ±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 ±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
β
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 ρ
B
of the cavity mode due to the
interaction with an external thermal field [29]:
dρ
B
dt
= iE
a
a, ρ
B
1
2
An
aa
ρ
B
2a
ρ
B
a + ρ
B
aa
1
2
A (n + 1)
a
B
2
B
a
+ ρ
B
a
a
. (51)
Here A is the rate of loss of cavity photons (controlling the strength of the re-thermalization) and n =
1
e
β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. (6) with = 1G
0|1
(˜s)),
with the cooling achieved for various finite re-thermalization times. In each case, we fix βE = 1 and a particular
interaction time between the qubit and the cavity mode, ˜s = g
˜
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
β
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.
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. First Pauli X applied.
2. First JC interaction applied, for some time ˜s.
3. Second Pauli X applied.
Accepted in Quantum 2019-09-16, click title to verify 22
0 2 4 6 8 10 12
Number of rounds
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Ground state population
t
th
=0
t
th
=0.1
t
th
=0.5
t
th
=1
t
th
=
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
th
. The infinite time corresponds to the mode being completely re-thermalized at every step (light
blue curve at the top), while t
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 βE = 1, A = 1, g = 1 and the
Jaynes-Cummings interaction is turned on for a period of time
˜
t = 98.92.
4. Second JC interaction applied, for the same time ˜s.
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
(˜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
β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 = , the
single mode has time to re-thermalize perfectly, the performance is the same for each atom and given by Eq. (6)
with k = 2 and = 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.
Acknowledgments
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
Accepted in Quantum 2019-09-16, click title to verify 23
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
[1] 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.
[2] 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.
[3] 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.
[4] 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.
[5] 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.
[6] 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.
[7] 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.
[8] 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.
[9] 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.
[10] 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.
[11] Albert W Overhauser. Aw overhauser, phys. rev. 89, 689 (1953). Phys. Rev., 89:689, 1953. DOI: 10.1103/Phys-
Rev.89.689.
[12] 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.
[13] M. Horodecki and J. Oppenheim. Fundamental limitations for quantum and nanoscale thermodynamics. Nat.
Commun., 4:2059, June 2013. DOI: 10.1038/ncomms3059.
[14] Francesco Ticozzi and Lorenza Viola. Quantum resources for purification and cooling: fundamental limits and
opportunities. Scientific Reports, 4:5192, 2014. DOI: 10.1038/srep05192.
[15] 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.
[16] 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.
[17] 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.
[18] Johan Åberg. Catalytic coherence. Phys. Rev. Lett., 113:150402, Oct 2014. DOI: 10.1103/Phys-
RevLett.113.150402.
Accepted in Quantum 2019-09-16, click title to verify 24
[19] 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.
[20] Herbert Walther, Benjamin TH Varcoe, Berthold-Georg Englert, and Thomas Becker. Cavity quantum elec-
trodynamics. Reports on Progress in Physics, 69(5):1325, 2006. DOI: 10.1088/0034-4885/69/5/R02.
[21] 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.
[22] 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.
[23] Garrett Birkhoff. Tres observaciones sobre el algebra lineal. Univ. Nac. Tucumán Rev. Ser. A, 5:147–151, 1946.
[24] Paw 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.
[25] 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.
[26] D. Janzing, P. Wocjan, R. Zeier, R. Geiss, and Th. Beth. Thermodynamic cost of reliability and low tem-
peratures: Tightening Landauer’s principle and the second law. Int. J. Theor. Phys., 39(12):2717–2753, 2000.
DOI: 10.1023/A:1026422630734.
[27] 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.
[28] 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.
[29] 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.
Accepted in Quantum 2019-09-16, click title to verify 25