---
abstract: 'Heat-bath algorithmic cooling (HBAC) provides algorithmic ways to improve the purity of quantum states.  These techniques are complex iterative processes that change from each iteration to the next and this poses a significant challenge to implementing these algorithms.  Here we introduce a new technique that on a fundamental level, shows that it is possible to do algorithmic cooling and even reach the cooling limit without any knowledge of the state and using only a single fixed operation, and on a practical level, presents a more feasible and robust alternative for implementing HBAC.  We also show that our new technique converges to the asymptotic state of HBAC, however, this could require exponentially many iterations and while the cooling algorithm can be efficiently implemented, reaching the cooling limit remains impractical.  This brings HBAC to the realm of feasibility and makes it a viable option for realistic application in quantum technologies.' 
author:
- Sadegh Raeisi
- Mária Kieferová
- Michele Mosca
bibliography:
- 'Cooling\_Bound\_AC.bib'
title: Novel Technique for Robust Optimal Algorithmic Cooling
---

Many quantum effects and quantum technologies rely on fragile quantum fluctuations that can easily be overwhelmed by thermal fluctuations.  This is why often these quantum technologies rely on techniques for suppressing thermal fluctuations, e.g. cooling the quantum system in a cryostat refrigerator or a dilution fridge.  There are however algorithmic techniques for cooling that more surgically extract energy from the spins of interest and enable temperatures to be lowered beyond what would be feasible with conventional cooling of the entire system. 

Heat-bath algorithmic cooling (HBAC) are techniques that operate on an ensemble of qubits and effectively cools down and purifies a target subset of qubits.  HBAC drives the system out of equilibrium by transferring the entropy from the target qubits to the rest of the ensemble.  This is often referred to as “compression” since it uses information theoretical techniques to compress the entropy to the non-target elements of the ensemble and effectively cools down the target qubits.  The target and the refrigerant qubits are often referred to as the “computation” and the “reset” qubits respectively. 

This idea was first introduced in [@schulman_scalable_1998] for a closed system, for which, the cooling is limited by the Shannon bound for information compression.  This process heats up the reset qubits beyond their initial temperature.  It was later proposed to use a heat-bath to recycle the reset qubits and enhance the cooling beyond the Shannon bound [@boykin_algorithmic_2002].  In this setting, the reset qubits, through the interaction with a heat-bath, are cooled down to the heat-bath temperature again.  This is known as the “reset step”. 

Even with the help of a heat-bath, it is not possible to extract all the entropy from the target qubits [@schulman_physical_2005].  The optimal technique was introduced by Schulman et al.  in [@schulman_physical_2005] and is known as the partner pairing algorithm (PPA). 

The existence of the limit was proved by Schulman et al.  in [@schulman_physical_2005].  Later Raeisi and Mosca [@raeisi2015asymptotic] established the asymptotic limit of PPA, with the corresponding asymptotic state, and proved that the process asymptotically approaches the cooling limit.  We refer to the optimal asymptotic cooling state as OAS. 

One of the main challenges with HBAC techniques, especially the ones that converge to OAS, is that they are highly complex iterative processes.  The operations change from each iteration to the next and in many cases, there is no recipe for implementing the operations in each step. 

For instance, PPA sorts the diagonal elements of the density matrix in each iteration.  But this means that not only one needs to know the state in each iteration, but also the operation for implementing the sort would change in each iteration as the state changes. 

An obvious question is whether it would be possible to reach the OAS and the cooling limit with a fixed state-independent operation in each iteration.  Note that the state is constantly changing through the cooling process, and naturally, the compression should change too, as is the case with PPA.  Reaching the OAS with a fixed operation seems even more non-trivial.  In other words, for a HBAC technique with fixed operation, the compression should be tuned such that without knowing the state, not only it can extract entropy from the state, but through the repetition of the process, it would push the state into the OAS.  Here we answer this question and show that this is in fact possible.  We introduce an compression operation that can push the state to the OAS and reach the cooling limit of HBAC, without any knowledge of the state.  Further, we show that this fixed operation can be implemented efficiently and give a recipe for building the quantum circuit. 

Besides the fundamental significance, this result could have a critical impact on the feasibility of HBAC techniques.  First, in contrast to techniques such as PPA, our algorithm can be efficiently implemented.  We however show that reaching the OAS would require exponentially many iterations.  Second, the state-independence of operations, makes our algorithm simple and more robust and turns HBAC to a viable option for generating large scale supplies of high-purity quantum states.  We start by introducing our algorithm and then present the complexity analysis.  Next, we analyse the complexity of PPA and argue that the classical and quantum gate complexity are at least exponential in the number of qubits.  We also investigate the robustness of PPA and TSAC and prove that our technique always increases the purity, even in the presence of noise, whereas PPA, with noise, could end up decreasing the polarization. 

We use the standard framework of HBAC techniques where we have an ensemble of $n+1$ qubits.  We assume that the last qubit is the reset qubit and the rest are computation qubits and use the subscript $R$ and $C$ to refer to the reset and the computation qubits.  Also we assume that the Hilbert space is ordered as $\mathcal{H}_{C}\otimes\mathcal{H}_{R}$, i.e. the first part are the computation qubits and the last part are the reset qubits. 

In our technique, instead of sorting the diagonal elements, we apply the following unitary for compression in each iteration

$$U_{TS}=\left(\begin{array}{ccccc}
1\\
 & X\\
 &  & \ddots\\
 &  &  & X\\
 &  &  &  & 1
\end{array}\right),\label{eq:Bi-sort_U}$$

where $X$ is the Pauli X operator and the first and the last elements of the matrix are the number one.  The matrix is $2^{n+1}\times2^{n+1}$ and acts on both the computation and the reset qubits.  We refer to the unitary $U_{TS}$ as the two-sort operator and to our technique as two-sort algorithmic cooling (TSAC).  The unitary $U_{TS}$ swaps every two neighboring elements on the diagonal of the density matrix, except for the first and the last elements.  Intuitively, this is a partial sort that acts locally on the density matrix.  This is the golden operation that makes it possible to cool the state to OAS without knowing the state. 

After compression, the reset qubit is reset which is equivalent to $\mathcal{R}\left[\rho\right]=\text{Tr}_{R}\left(\rho\right)\otimes\rho_{R}$, where $\text{Tr}_{R}$ is the partial trace over the reset qubit and $\rho_{R}$ is the “reset state” $$\rho_{R}=\frac{1}{z}\left(\begin{array}{cc}
e^{\epsilon} & 0\\
0 & e^{-\epsilon}
\end{array}\right),\label{eq:reset-state}$$ with $z=\left(e^{\epsilon}+e^{-\epsilon}\right)$. The parameter $\epsilon$ is called the polarization. 

Mathematically, each iteration applies the following channel on the full density matrix $$\label{eq:TSAC_Iteration}
\mathcal{C}_{TS}\left[\rho\right]=U_{TS}^{\dagger}\left(\text{Tr}_{R}\left(\rho\right)\otimes\rho_{R}\right)U_{TS}.$$

This process is clearly independent of the iteration or the state and it can be described by a time-homogeneous Markov process.  We find the transfer matrix and use its spectrum to prove that the process converges to the OAS and to give an upper-bound for the required number of iterations. 

The sequences of the elements on the diagonal of the density matrix form a Markov chain.  We use a vector with $2^{n+1}$ elements $\left\{ \lambda^t\right\}$ to represent the state after the $t$th iteration.  We use a similar notation for the density matrix of the computation qubits (without the reset qubit) and use $\left\{ p^t\right\} $ instead. 

Figure (\[fig:Pictorial\_Iterations\]) gives a pictorial description of the process in each iteration.  It starts with the sequence $\{ \lambda^{t} \}$, the diagonal elements of the density matrix of the $n$ computation and one reset qubit in the $t^{\text{th}}$ iteration.  First, there is the reset step which takes the reset qubit to the state in Eq. ( \[eq:reset-state\]).  This takes every two neighboring elements $\lambda^{t}_{2k+1}$ and $\lambda^{t}_{2k+2}$ to $p^t_k= \lambda^{t}_{2k+1}+\lambda^{t}_{2k+2}$ and then splits them into $\zeta^{t}_{2k+1}=p^t_k e^{ \epsilon}/z$ and $\zeta^{t}_{2k+2}=p^t_k e^{- \epsilon}/z$. Now the two-sort unitary is applied and rearranges the array to $\{ \lambda^{t+1} \}$ such that $\lambda^{t+1}_{2k}=\zeta^{t}_{2k+1}$ and $\lambda^{t+1}_{2k+1}=\zeta^{t}_{2k}$.

![ \[fig:Pictorial\_Iterations\] The pictorial description of an iteration.  The input is the list of the diagonal elements of the full density matrix, $\{ \lambda^{t} \}$. The reset step first merges every two neighbouring elements (partial trace) and then replaces the reset qubit with $\rho_R$ which splits each element into two elements again.  Next, the $U_{TS}$ swaps all the neighbouring elements except for the first and the last one and gives the final state. ]( Figure2.pdf){width="\columnwidth"}

For simplicity, we focus on the computation qubits and trace out the reset qubit.  This gives the following update rule for the diagonal elements of the computation qubits $$p_{i}^{t+1}=p_{i-1}^{t}\frac{e^{-\epsilon}}{z}+p_{i+1}^{t}\frac{e^{\epsilon}}{z}.\label{eq:Update_Rule}$$ for $1<i<2^n$.

Similarly for the first and the last elements, the update rules are $p_{1}^{t+1}=(p_{1}^{t}+p_{2}^{t})\frac{e^{\epsilon}}{z}$ and $p_{2^{n}}^{t+1}=(p_{2^{n}-1}^{t}+p_{2^{n}}^{t})\frac{e^{-\epsilon}}{z}$.

These update rules give the following transition matrix for the Markov process $$T=\frac{1}{z}\left(\begin{array}{ccccc}
e^{\epsilon} & e^{\epsilon} & 0 & \cdots & 0\\
e^{-\epsilon} & 0 & e^{\epsilon} & \cdots & 0\\
0 & e^{-\epsilon} & 0 & \cdots & 0\\
0 & 0 & \cdots & \ddots & \vdots\\
0 & 0 & \cdots & e^{-\epsilon} & e^{-\epsilon}
\end{array}\right).\label{eq:MarkovMatrix}$$

It is easy to verify that $\vec{\{p^{t+1}\}} =T \vec{\{p^t \}}$ and gives the update rules above.  The matrix $T$ has a unique eigenvalue $1$ and the remaining eigenvalues are ${\Upsilon}_{k}=\frac{2\cos{\frac{l\pi}{2^n}}}{z}$ for $k=1,2, \cdots, 2^n-1$. The eigenstate corresponding to eigenvalue one is $$\rho =p_0 \lbrace1,e^{-2\epsilon}, e^{-4\epsilon}, \cdots\rbrace,
\label{OAS}$$ which is the OAS and $p_0$ is the normalization factor [@raeisi2015asymptotic].  For the detailed calculation of the eigensystem, see the Supplemental Material(SM).  Since all the other eigenvalues lie in the interval $(1, -1)$, the Markov chain asymptotically converges to $\rho$. This proves that our technique asymptotically achieves the cooling limit of HBAC. 

Figure (\[fig:Two-Sort\_Unitary\_Circuit\]) gives a $O\left(n^{2}\right)$ circuit for the implementation of the two-sort unitary.  We first shift the basis of the two-sort unitary by one, which transforms it to $\text{Toff}(n+1) \sigma_x^{(n+1)}$. Then a Pauli $\sigma_x$ on the reset qubit turns the matrix to a multiple-control Toffoli gate. ( See SM for details).  This shows that our technique can be efficiently implemented.  However, to reach the OAS, we need to investigate how many iteration would be required.  The mixing time of a Markov chain gives the number of iterations required to get within $\xi$ distance of the asymptotic state.  We use the spectral gap, i.e. the difference between 1 and the second largest eigenvalue, to upper-bound the mixing time.  The mixing time is upper-bounded by $$\label{eq:mixingtime_upperbound}
t_{mix}\left(\xi\right)\leq\log\left(\frac{1}{\xi\, l}\right)\frac{1}{\Delta},$$ where $\Delta$ is the spectral gap and $l=p_0 e^{-(2^n-1)\epsilon}$ is the smallest element of the array in Eq. ( \[OAS\]) [@levin2009markov]. 

The spectral gap is $\Delta = 1-\left(2\cos{\frac{\pi}{2^n}}\right)/\left(e^{\epsilon}+e^{-\epsilon}\right)\geq \frac{z-2}{z}$. This reduces the bound to $$t_{mix} (\xi) \leq \left( \log\left(\frac{1}{\xi \, p_0 e^{-(2^n-1)\epsilon}}\right) \left(\frac{z}{z-2}\right) \right).$$ It is easy to check that $t_{mix} \in O(2^{n})$. This gives $O(n^2 2^n)$ for the overall complexity of TSAC. 

![ \[fig:Two-Sort\_Unitary\_Circuit\]The circuit for implementing the two-sort operation. ]( Figure3.pdf){width="0.5\columnwidth"}

Now, we compare our technique with PPA.  The key ingredient of PPA is sorting the diagonal elements of the density matrix in each iteration.  This transfers as much entropy as possible from the computation elements to the reset qubit [@schulman_physical_2005].  Then the reset qubit is reset back to its equilibrium state. 

Typically for the cost analysis of PPA, only the number of iterations is taken into account and the implementation cost of the sort operations is not addressed.  We will present a more careful complexity analysis of PPA.  There are three different cost that needs to be considered, first the classical cost of finding the circuits for the sort operations, second, the gate complexity of the sort operations, and last, how well the process performs in the presence of noise. 

Each iteration of PPA takes a density matrix as input and sort the diagonal elements.  This requires sorting elements of a $2^{n+1}$ array.  Even assuming that that finding the sort operation is easy we need to find quantum circuits to implement them.  This would scale as $O(e^n)$ at best.  This is only the classical cost of the algorithm.  Without this, it would not be even possible to start implementing PPA.  For systems as large as 20-40 qubits, such as the one in [@pande2017strong], it is not even possible to get the process started.  This is in contrast to our technique, where each iteration is already known and there is a specific circuit for implementing it. 

Next is the gate complexity of PPA.  Due to the complexity of the sort operations, it is difficult to bound the number of gates required for PPA.  Naively, there are $(2^n)!$ permutation matrices of size $2^n$. Assuming a finite number of one- and two-qubit gates, there are only $\left(\text{poly}(n)\right)^d$ circuits with at most $d$ gates on $n$ qubits.  Taking $d$ to be $\text{poly}(n)$, one can see that only a small fraction of permutations can be implemented efficiently. 

The permutation operations can be decomposed into separate cycles that form disjoint blocks in the permutation matrix.  These cycles could have different sizes and the size of each cycle determines the number of states it permutes cyclically.  These are known as $k$-cycles, where $k$ is the size of the block, or equivalently the number of states that are permuted cyclically by the corresponding block. 

Assume that for all the $k$-cycles in the permutation, we can find an efficient circuit.  Also assume that, given a certain state ${\left| j\right\rangle}$, it is possible to efficiently determine which block the state belongs to.  Furthermore, we assume that cycles of equal size can be implemented in parallel efficiently.  These assumptions may not be true, but any lower bound established with these assumptions still hold when any of these assumptions are weakened or dropped.  Under these assumptions, the cost would depend on the number of $k$-cycles with distinct $k$ values.  This is the number of blocks in the permutation matrix that have different size.  We refer to this quantity as NBDS. 

Implementation of each sort operation requires the implementation of all the blocks.  Blocks of different size cannot be fully parallelised and for switching between each two blocks of unequal size, some quantum operation would be required.  This sets the number of blocks of different size, NBDS, as a lower bound for the complexity of any sort operation. 

Figure (\[fig:NBDS\]) shows the simulation results of PPA for different number of computation qubits, $n$ and indicates that NBDS grows exponentially with $n$ and implies that our lower-bound for the gate complexity of PPA scales exponentially with $n$. Here for any value of $n$, we get a sequence of permutation matrices and pick the permutation that has the largest NBDS. 

Last, there is the fragility to practical imperfections.  Ideally, all the operations in PPA can be pre-computed.  However, with even a small deviation from the assumed input state for an iteration, the pre-computed sort operation could end up heating the computation qubits, instead of cooling them.  To account for the imperfections, techniques like quantum state tomography are required to monitor the state which are too expensive and challenging and even then, the success of the process cannot be guaranteed. 

For instance, imagine that state tomography is used to determine the state.  Experimentally the tomography is not perfect and involves some estimation errors.  Figure (\[fig:Noisy\_Tomography\]) shows how these estimation errors can affect the cooling process and how sensitive the process is to these imperfections.  These simulations are for HBAC with $n=2$ and one reset qubit and the reset polarization of $\epsilon = .02$ and include a Gaussian noise where each element of the diagonal of the density matrix is determined with a Gaussian error.  The parameter $\sigma$ is the width of the Gaussian distribution.  As $\sigma$ increases, the process becomes random and would not approach the cooling limit any more.  Figure (\[fig:Noisy\_Tomography\]) also shows the result for TSAC which regardless of the noise, would always converge to the OAS.  This is not a generic noise model, but is relevant for techniques like PPA and shows that with noise, PPA can heat the state.  Detailed investigation of more realistic noise models for HBAC techniques is required, but is out of the scope of this work.  Next we show that TSAC monotonically pushes the state towards the OAS. 

\[thm:Robustness\] Given some state $\rho$ and a reset state $\rho_R$ with polarization $\epsilon$, if the polarization of the first computation qubit is less than the HBAC limit, each iteration of TSAC, as in equation (\[eq:TSAC\_Iteration\]) would increase the polarization of the first qubit. 

The polarization of the first qubit is determined by the first half of the diagonal elements of $rho$ and we need to show that under TSAC, it would increase.  The assumption that the state is hotter than the HBAC limit, translates to $$p_{2^{n-1}}e^{-\epsilon} < p_{2^{n-1}+1}e^{\epsilon}.$$ After the iteration, these two are swapped by the $U_{TS}$ which increases the sum of the first half of the diagonal elements and as a result would increase the polarization of the first qubit. 

Note that this can be extended to other computation qubits.  It is just easier to show for the first qubit. 

Since the process does not depend on the state, even when the state is perturbed from the ideal one, the process continues to cool it down.  Note that this does not imply that our algorithm is robust to all imperfections.  Specifically, with faulty operations, no algorithm can guarantee the convergence to OAS.  But this is true for any algorithm without a fault-tolerant implementation. 

Table (\[tab:Comparison\]) gives a comparison between TSAC and PPA.  Here the “noise sensitivity” is the sensitivity to deviations from the expected state in the process.  The table demonstrates that TSAC outperforms PPA in almost every aspect and presents a more realistic option in practice. 

In conclusion, our work presents a novel viable technique for optimal HBAC which shows that optimal HBAC is possible without any knowledge of the state and without changing the operation through the process.  It also brings using these techniques for realistic applications to the realm of possibility.  The new technique, in contrast to PPA, uses a fixed operation in every iteration, which addresses the fragility issues in previous works.  More precisely, the new technique is robust against imperfections and noise in the state.  It is also possible to combine our work with other HBAC techniques such as $3BC$ [@kaye2007cooling] or $PAC$[@fernandez2004algorithmic] to further reduce the cost.  This could speed up the convergence and reduce the gate complexity.  However, it remains open to see how far the complexity may be reduced. 

We thank Alex Parent for helpful discussions.  This work was supported by ERC Starting Grant OPTOMECH, Canada’s NSERC, MPrime, CIFAR, CFI and the research grant system of Sharif University of Technology (G960219).  IQC and Perimeter Institute are supported in part by the Government of Canada and the Province of Ontario. 

Supplementary Material
======================

Spectrum of the transfer matrix {#spectrum-of-the-transfer-matrix .unnumbered}
-------------------------------

We solve the eigenvalue equation, $T \Phi^{(p)}={\Upsilon}\Phi^{(p)}$, indexing the eigenvectors by $p$. Using the sparsity and the structure of $T$, we can rewrite the eigenvalue equations as

$$\begin{aligned}
{\Upsilon}^{(p)}\Phi_{1}^{(p)} &=(\Phi_{1}^{(p)}+\Phi_{2}^{(p)})\frac{e^{\epsilon}}{z} \label{eq:FirstT},\\
{\Upsilon}^{(p)}\Phi_{k}^{(p)} &= \Phi_{k-1}^{(p)}\frac{e^{-\epsilon}}{z}+\Phi_{k+1}^{(p)}\frac{e^{\epsilon}}{z}\label{eq:Update_RuleT},\\
{\Upsilon}^{(p)}\Phi_{2^{n}}^{(p)} &=(\Phi_{2^{n}-1}^{(p)}+\Phi_{2^{n}}^{(p)})\frac{e^{-\epsilon}}{z}\label{eq:LastT},\end{aligned}$$

for $1<k<2^n$. We use the ansatz $$\Phi_k^{(p)} = e^{(ip-\epsilon)k} + \alpha e^{(-ip-\epsilon)k}$$ with arbitrary complex parameters $\alpha$ and $p$. This ansatz automatically satisfies \[eq:Update\_RuleT\] with eigenvalue $${\Upsilon}^{(p)}=\frac{2\cos{p}}{e^{\epsilon}+e^{-\epsilon}}.$$

We set the value of $\alpha$ by solving \[eq:FirstT\] and obtain $\alpha=\frac{e^{ip}-e^{-\Delta}}{e^{-\Delta}-e^{-ip}}$. Note that this result forbids $p=0$ because it gives $\Phi_k^{(0)}=0$.

At last, we satisfy \[eq:LastT\] and obtain allowed values of $p$. The solution $ip=\pm\epsilon$ gives eigenvalue $1$ and corresponds to the eigenvector $\Phi_{k}^{(p)} = e^{-2\Delta k}$. This is the asymptotic state of PPA as was proved in [@raeisi2015asymptotic]. 

The remaining eigenvalues are of form $\frac{2\cos{\frac{j\pi}{2^n}}}{e^{\epsilon}+e^{-\epsilon}}$ for $1\leq j < 2^n$. All these eigenvalues lie in the range $(-1,1)$. In other words, the Markov chain has a unique eigenvalue one and all the other eigenvalues are smaller than one.  Therefore the Markov chain defined by the transition matrix $T$ converges to the +1 eigenvector, which is OAS. 

The convergence rate is determined by the difference between $1$ and the second largest eigenvalue, ${\Upsilon}^{(2)}=\frac{2\cos{\frac{\pi}{2^n}}}{e^{\epsilon}+e^{-\epsilon}}$. We can bound the gap as $\Delta = 1-\left(2\cos{\frac{\pi}{2^n}}\right)/\left(e^{\epsilon}+e^{-\epsilon}\right)\geq \frac{z-2}{z}$. The mixing time is then upper-bounded by $$\label{eq:SM_mixingtime_upperbound}
t_{mix}\left(\xi\right)\leq\log\left(\frac{1}{\xi\, l}\right)\frac{1}{\Delta} \leq \left( \log\left(\frac{1}{\xi\, l}\right) \left(\frac{z}{z-2}\right) \right) \leq c_1 \log\left(\frac{1}{ l}\right)+c2,$$ where $c_1=\left(\frac{z}{z-2}\right)$ and $c_2= \left(\frac{z}{z-2}\right) \log\left(\frac{1}{\xi}\right)$ are both constant with respect to $n$. To find the scaling of the upper-bound, we need to calculate the

$$\log\left(\frac{1}{ l}\right)=\log\left(\frac{1}{ p_0 e^{-(2^n-1)\epsilon}}\right)=\log\left(\frac{(1-(e^{-2\epsilon})^{2^n})}{ (1-e^{-2\epsilon}) e^{-(2^n-1)\epsilon}}\right).$$ For the values of $l$ and $p_0$ see [@raeisi2015asymptotic]. 

To understand the scaling, we take $n\gg1$ which simplifies the bound to $$\log\left(\frac{1}{ e^{-(2^n-1)\epsilon}}\right)-c_3,$$ where $c_3=\log\left(1-e^{-2\epsilon}\right)$. So the scaling of the upper-bound is $O(2^n)$.

The Circuit for $U_{TS}$ {#the-circuit-for-u_ts .unnumbered}
------------------------

The unitary that sorts the density matrix lexicographically is $$\begin{aligned}
U _{TS}=
 \begin{pmatrix}
  1 & 0 & 0 & 0 & 0 &\cdots & 0 & 0 & 0 \\
  0 & 0 & 1 & 0 & 0 &\cdots & 0 & 0 & 0 \\
  0 & 1 & 0 & 0 & 0 &\cdots & 0 & 0 & 0 \\
  0 & 0 & 0 & 0 & 1 &\cdots & 0 & 0 & 0 \\
  0 & 0 & 0 & 1 & 0 &\cdots & 0 & 0 & 0 \\
  \vdots  & \vdots & \vdots   & \vdots & \vdots & \ddots &   \vdots  & \vdots & \vdots  \\
  0 & 0 & 0 & 0 & 0 &\cdots & 0 & 1 & 0 \\
  0 & 0 & 0 & 0 & 0 &\cdots & 1 & 0 & 0 \\
  0 & 0 & 0 & 0 & 0 &\cdots & 0 & 0 & 1
\end{pmatrix}.  \label{SM:sorting_unitary}\end{aligned}$$ $U$ is a block diagonal $2^n\times 2^n$ matrix with $1\times1$ blocks in the upper left and lower right corners and the other blocks are $X=\begin{pmatrix} 0 & 1 \\ 1 & 0\end{pmatrix}$.

The first and the last operations are SHIFT$_m$ operators which are defined as $$\text{SHIFT}_m{\left| x_1x_2\dots x_n\right\rangle}={\left| \left(x_1x_2\dots x_n + m\right) mod 2^n\right\rangle}. 
\label{SM:Shift}$$ This notation is mostly symbolic but should be understood as binary addition between strings labeling states on $n$ qubits.  For example, SHIFT$_1{\left| 101\right\rangle}={\left| 110\right\rangle}$. The other operations are multiple-control-Toffoli, Toff$_n$, which is a controlled NOT with $n-1$ controls and NOT (Pauli X) on the last qubit. 

NOTs are easy to implement and for SHIFT$_1$ and Toff$_n$ we use the construction in [@saeedi2013linear]. 

After Toff$_n$ we apply NOT on the last qubit.  Multiple-control-Toffoli and NOT give together the unitary $$\begin{aligned}
U'=
 \begin{pmatrix}
  0 & 1 & 0 & 0 &\cdots & 0 & 0 & 0 & 0\\
  1 & 0 & 0 & 0 &\cdots & 0 & 0 & 0  & 0\\
  0 & 0 & 0 & 1 &\cdots & 0 & 0 & 0  & 0\\
  0 & 0 & 1 & 0 &\cdots & 0 & 0 & 0  & 0\\
  \vdots  & \vdots   & \vdots & \vdots & \ddots &   \vdots & \vdots  & \vdots & \vdots  \\
  0 & 0 & 0 & 0 &\cdots & 0 & 1 & 0 & 0\\
  0 & 0 & 0 & 0 &\cdots & 1 & 0 & 0 & 0\\
  0 & 0 & 0 & 0 &\cdots & 0 & 0 & 1 & 0\\
  0 & 0 & 0 & 0 &\cdots & 0 & 0 & 0 & 1\\
\end{pmatrix}.  \label{middle_unitary}\end{aligned}$$ Therefore, to apply U, we must shift all the rows and columns of $U'$ cyclically.  This is what SHIFT$_{+1}$ and its conjugate transpose SHIFT$_{-1}$ do.  We can implement SHIFT$_{+1}$ with Quantum Fourier Transform and rotations [@beth2001quantum].  In Fig (\[shift\]), we show how to implement a more general operation, SHIFT$_{+m}$.

In this circuit, first, QFT$^{-1}$ transforms the bit-strings from registers to phases $${\left| x_{1}\right\rangle}{\left| x_{2}\right\rangle}\dots{\left| x_{n}\right\rangle}\rightarrow \frac{1}{\sqrt{2^{n}}}
\left({\left| 0\right\rangle}+e^{2\pi i\frac{x_{1}x_{2}\dots x_{n}}{2^{n}}}{\left| 1\right\rangle}\right)
\left({\left| 0\right\rangle}+e^{2\pi i\frac{x_{2}\dots x_{n}}{2^{n-1}}}{\left| 1\right\rangle}\right) \cdots
\left({\left| 0\right\rangle}+e^{2\pi i\frac{x_{n}}{2}}{\left| 1\right\rangle}\right).$$

In the next step, we apply rotations around $z$ on each qubit.  This yields $$\frac{1}{\sqrt{2^n}}\left({\left| 0\right\rangle} + e^{2\pi i\frac{ x_1x_2\dots x_n + m }{2^n}}{\left| 1\right\rangle} \right)\left({\left| 0\right\rangle} + e^{2\pi i\frac{ x_2\dots x_n + m}{2^{n-1}}}{\left| 1\right\rangle} \right)\cdots \left({\left| 0\right\rangle} + e^{2\pi i\frac{x_n + m}{2}}{\left| 1\right\rangle} \right).$$ After applying the Quantum Fourier Transform, we get the desired state ${\left| x_1x_2\dots x_n +m\right\rangle}$.

![ The circuit for the SHIFT operator in Eq. ( \[SM:Shift\]).[]{data-label="shift"}](Figure4_SM.pdf){width="0.4\columnwidth"}

Complexity of the sort operations in PPA {#complexity-of-the-sort-operations-in-ppa .unnumbered}
----------------------------------------

First, note that the implementation cost of a sort operation could in general be exponential in $n$. Here we give a lower bound for the implementation cost that takes into account the cost of the sort operations.  We use numerical evidence to show that the lower-bound scales exponentially with $n$.

A permutation matrix basically permutes the basis and we need to find the circuit with the minimum number of gates that gives the same basis transformation.  Clearly, mapping every element of the basis would be exponential in $n$ since there are $2^n$ basis to map.  So we need to parallelize as many mappings as possible. 

As explained in the text, we assume that all the cycles of the permutation have efficient implementation and for cycles of the same size, it is easy to implement them in parallel.  Without parallelizing the implementation of blocks that have different sizes, the complexity would be proportional to the number of blocks with different size (NBDS). 

The numerical evidence presented in figure (5) in the text, shows that NBDS scales exponentially with $n$.

Non-diagonal density matrices {#non-diagonal-density-matrices .unnumbered}
-----------------------------

It is possible that the density matrix becomes non-diagonal during the cooling process due to the noise in the system.  Here we show that our technique can handle non-diagonal density matrices as well and would take them to the OAS. 

We show that the presence of any non-diagonal elements will not affect the action of $U_{TS}$ on the diagonal of the density matrix.  In other words, the same permutation is applied on the diagonal elements even if the density matrix is non-diagonal in computational basis. 

We can write $U_{TS}$ as $$U_{TS} = {\left| 0\right\rangle}{\left\langle 0\right|} + {\left| 2^{n}-1\right\rangle}{\left\langle 2^{n}-1\right|} + \sum_{i\in odd}^{2^{n}-2}{({\left| i+1\right\rangle}{\left\langle i\right|}+{\left| i\right\rangle}{\left\langle i+1\right|})}.$$ Now, if we apply the evolution to some arbitrary density matrix $\rho = \sum_{i,j}{\rho_{ij} {\left| i\right\rangle}{\left\langle j\right|} }$ , for the diagonal elements we get $${\left\langle i\right|} U_{TS}\rho U^{\dagger}_{TS}{\left| i\right\rangle} = \rho_{0,0}{\left| 0\right\rangle}{\left\langle 0\right|} + 
\rho_{2^{n}-1,2^{n}-1}{\left| 2^{n}-1\right\rangle}{\left\langle 2^{n}-1\right|} + 
\sum_{i\in odd}{(\rho_{i,i}{\left| i+1\right\rangle}{\left\langle i+1\right|}+\rho_{i+1,i+1}{\left| i\right\rangle}{\left\langle i\right|})}.$$ The first two terms show that the first and last elements of the diagonal elements of the density matrix are preserved, while the last term indicates that for all odd values of $i$, we get $\rho_{i,i} \leftrightarrow \rho_{i+1,i+1}$. This gives the same permutation for the diagonal element as for a diagonal density matrix.  This means that the diagonal element would again follow a time-independent Markov chain given by the same transfer matrix, T (Eq. ( 5) of the main text). 

Note that polarization only depends on the diagonal elements of the density matrix, so it suffices to focus on the dynamics of these diagonal elements.  As was shown here, this dynamics would be the same regardless of off-diagonal elements of the density matrix.