Algorithmic Cooling of Nuclear Spins using Long-Lived Singlet Order

Algorithmic cooling methods manipulate an open quantum system in order to lower its temperature below that of the environment. We achieve signiﬁcant cooling of an ensemble of nuclear spin-pair systems by exploiting the long-lived nuclear singlet state, which is an antisymmetric quantum superposition of the ”up” and ”down” Zeeman states. The eﬀect is demonstrated by nuclear magnetic resonance (NMR) experiments on a molecular system containing a coupled pair of near-equivalent 13 C nuclei. The populations of the system are subjected to a repeating sequence of cyclic permutations separated by relaxation intervals. The long-lived nuclear singlet order is pumped well beyond the unitary limit. The pumped singlet order is converted into nuclear magnetization which is enhanced by 21% relative to its thermal equilibrium value.


I. INTRODUCTION
A quantum system which is in contact with a thermal reservoir eventually comes into equilibrium, so that the populations of the quantum states are described by the Boltzmann distribution at the reservoir temperature. In the case of homonuclear spin systems at ordinary temperatures in a high magnetic field B 0 along the z-axis of the laboratory axis system, the high-temperature (or weak order) approximation applies [1][2][3] , so that the thermal equilibrium state of the spin ensemble may be described by a density operator of the form where N H is the dimension of Hilbert space (given by the number of independent quantum states of each spin system), ω 0 = −γB 0 is the nuclear Larmor frequency, γ is the magnetogyric ratio, and I z is the spin angular momentum operator along the field. For an ensemble of N -spin-1/2 systems, the dimension of Hilbert space is N H = 2 N . The second term in equation 1 represents the nuclear spin order generated by thermal equilibration in a strong magnetic field. In the vast majority of nuclear magnetic resonance (NMR) experiments, this thermal spin order is manipulated by resonant radiofrequency pulses, and is eventually converted into transverse magnetization, which precesses in the magnetic field, and induces an electrical signal by Faraday induction. It is a prominent paradigm in NMR spectroscopy that the initial spin order established by thermal equilibration cannot be increased by the application of resonant radiofrequency pulses, but only transformed into different spin order modes, with the total quantity of spin order being either conserved, or decreased by irreversible decoherence effects.
There are, however, some situations in which resonant radiofrequency irradiation does increase the amount of nuclear spin order, at least locally. A seminal example is the steady-state Overhauser effect, in which radiofrequency irradiation which is resonant with one spin species may lead to an increase in the magnetization of a second spin species 4,5 . However, the Overhauser effect does not lead to decrease in entropy for the spin system as a whole, and also relies on the dominance of the dipole-dipole relaxation mechanism. More recently, a more general set of techniques called algorithmic cooling have been demonstrated. These methods systematically build up the order of an open quantum system through a deterministic sequence of events [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] . Algorithmic cooling methods have been applied with success to nuclear spin systems [15][16][17][18][19][20][21][22] . These methods utilise the ongoing contact of the nuclear spin system with the thermal environment during the experimental procedure. The nuclear spin system is an open quantum system which constantly exchanges energy and entropy with the thermal environment, also during a radiofrequency pulse sequence. The algorithmic cooling of nuclear spins emerged from the field of NMR quantum computation 20,23,24 , and treats the individual nuclear spins as independent qubits which may be addressed individually by frequency-selective radiofrequency fields.
One example of algorithmic cooling has been demonstrated on a spin system containing two different nuclear spin species, one of which relaxes rapidly (and is therefore strongly coupled to the thermal environment) while the second species relaxes slowly (and therefore has a much weaker environmental coupling) [15][16][17][18][19]21 . These spin species have been called "refrigerant" and "computation" qubits 9 . The magnetization of the slowly relaxing species is used to temporarily store the initial spin order, while the rapidly relaxing "refrigerant" species acquires a new batch of spin order through contact with the thermal environment. The total entropy of the nuclear spin system is thereby decreased below its initial value, leading to an increase in the total nuclear spin magnetization. However, this particular scheme is restricted to spin systems with two species having very different relaxation properties, and with sufficiently well-resolved spectra that each species may be addressed individually by controlling the frequencies of the radiofrequency pulses. These are very restrictive conditions.
In this paper we demonstrate the algorithmic cooling of an ensemble of spin-1/2 pair systems which exhibit nearmagnetic-equivalence [25][26][27][28] , meaning that a small difference in chemical environments induces a difference in Larmor frequencies which is much smaller than the scalar spin-spin coupling. This appears to be an unpromising substrate for algorithmic cooling for the following reasons: (i) the individual spins-1/2 (qubits) are tightly coupled, do not display resolved spectral resonances, and cannot be addressed selectively using resonant irradiation; (ii) two of the four energy eigenstates are entangled superpositions of the "up-down" and "down-up" qubit states; (iii) the dominant relaxation mechanisms, which couple the quantum system to the environment, are the dipole-dipole (DD) coupling and chemical shift anisotropy (CSA) mechanisms, which affect both spins equally. Unlike previous algorithmic cooling demonstrations, no distinction may be made between fast-relaxing "refrigerant" and slowly relaxing "computation" qubits 9 . Despite these obstacles, we succeeded in a demonstration of algorithmic cooling on a near-equivalent spin-pair system, by exploiting the extended lifetime of the singlet spin order.

II. THEORY
Consider an ensemble of spin systems, each consisting of two coupled spins-1/2 in a molecular environment such that they exhibit near magnetic equivalence 25-28 , meaning that a small difference in chemical environments induces a difference in Larmor frequencies which is smaller than the scalar spin-spin coupling. Under these conditions, the four Hamiltonian eigenstates are approximately equal to the singlet and triplet states of the spin-1/2 pair, numbered as follows: In many cases, spin-1/2 pairs of this type relax predominantly through the motional modulation of the dipoledipole coupling of the two spins (the DD relaxation mechanism). Since the DD coupling is mutual, the DD relaxation mechanism affects both spins equally, investing both spins with an identical relaxation time. Hence there seems to be no possibility of achieving algorithmic cooling. Furthermore, in the case of near magnetic equivalence, the individual spins do not have resolved spectra and cannot be regarded as independent qubits which may be manipulated selectively by radiofrequency fields.
Algorithmic cooling therefore appears to be doubly impossible for this system. Nevertheless, we demonstrate in this article an algorithmic cooling procedure which does work for strongly coupled homonuclear spin-1/2 pairs. The algorithmic cooling procedure exploits the singlet state |1 , which is an antisymmetric entangled state of the two spins-1/2. The singlet order of the system, defined as the mean population difference between the singlet state |1 and triplet states {|2 , |3 , |4 }, is protected against DD relaxation and may have a much longer lifetime than population differences within the triplet manifold 25,[29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48] . Techniques exist for coherent interconversion of singlet and triplet spin order, even in the near-equivalence limit 25,35 . Spin order established by the initial contact with the thermal molecular environment is stored temporarily as singlet order. Thermal contact of the spin system with the environment establishes a new batch of nuclear spin order during a subsequent interval, so that the total entropy of the system is lower than after the initial thermal equilibration step. Careful manipulations allow this enhanced spin order to be converted into enhanced nuclear magnetization. The result is that by the end of the pulse sequence, the nuclear magnetization exceeds that established initially by thermal equilibration, breaking the standard paradigm of most NMR experiments. This method also breaks the paradigm of most algorithmic cooling methods, which rely on the existence of individually accessible qubits.

A. A thought experiment
The basic principle is illustrated by the "thought experiment" in figure 1. Assume that the spin-lattice relaxation time constant T 1 is much shorter than the singlettriplet relaxation time constant T S , so that the populations of the triplet manifold {|2 , |3 , |4 } come rapidly into thermal equilibrium with the molecular environment, while the population of the singlet state |1 remains isolated. Figure 1(a) depicts the thermal equilibrium position of the spin ensemble, adopted when the spin system is left in contact with the lattice for a sufficiently long time. The filled balls in figure 1(a) represent positive deviations from a uniform population distribution, while the empty circles indicate negative population deviations. The excess population in the lower triplet state |2 , and depleted population in the upper triplet state |4 , indicate a finite magnetization along the field direction. An arbitrary population difference of +6 units is shown between the two extreme triplet states. This number is chosen for arithmetic convenience (see below).
The algorithmic cooling procedure starts by an exchange of populations between states |1 and |2 , corresponding to the operations: where φ and φ are arbitrary phase factors. In practice the population swap operation π 12 may be implemented with high fidelity by a sequence of applied radiofrequency fields, as described later. The consequence of swapping the populations of |1 and |2 is shown in figure 1 Note that the effect of the population swap is to deplete slightly the overall population of the triplet manifold {|2 , |3 , |4 }. The mean population of the triplet manifold in figure 1(b) is −1 unit, (relative to a uniform population distribution), as opposed to figure 1(a), where the mean deviation in triplet populations is zero. At the same, the relatively large population deviation of +3 units is "stored" in the long-lived singlet state |1 . The next step is to leave the spin system undisturbed in the magnetic field for a time long compared to T 1 , but short compared to T S . This allows the triplet manifold to come back into thermal equilibrium with the environment, while the singlet population is left approximately undisturbed. The result is the population distribution shown in figure 1(c). This is the same as the initial population distribution shown in figure 1(a), but superposed on the mean population depletion of −1 for the triplet manifold. The population difference between the lower and upper triplet states is restored to +6 units.
If the population swap π 12 is again executed, the resulting population distribution is as shown in figure 1(d).
Notably, the population difference between the lower and upper triplet states is now +7 units, i.e. a factor of 7/6 larger than its thermal equilibrium value. This is the same Zeeman polarization as would be exhibited by a system with a spin temperature equal to 6/7 of the environmental temperature. Although the effect is not large, this simple thought experiment establishes the principle of algorithmic cooling using long-lived spin states.
A more formal analysis of the thought experiment in figure 1 is provided by considering the transformations of a vector of state populations, defined as follows: where r p r = 1. The thermal equilibrium state at the beginning of the sequence may be written, in the hightemperature limit, as showing that the lower triplet state is slightly enhanced in population relative to the upper triplet state. The Zeeman and singlet orders of the spin ensemble are given in terms of the populations as follows: These quantities may be derived from the population vector p at any time through Z = Z · p, and S = S · p, The Zeeman order Z is proportional to the population difference between the lower and upper triplet states, while the singlet order S is given by the difference between the singlet population and the mean triplet population.
In thermal equilibrium, the Zeeman order is given by using the thermal equilibrium populations in equation 5.
Both sets of spins have the same spin temperature. There is no nuclear singlet order in thermal equilibrium: The first step in figure 1 may be represented by an application of the following permutation matrix to the population vector: The second step in figure 1 represents an idealised relaxation process, which brings the three triplet populations back into thermal equilibrium with the environment, while leaving the singlet population unchanged. We call this manipulation a triplet thermal reset, and represent it by the following matrix: In practice, the triplet thermal reset may be achieved approximately by leaving an appropriate relaxation interval, in the case that T S T 1 . Matrix-vector multiplications lead to the end result of the operations in figure 1: p = π 12 · Θ reset T · π 12 · p eq 1 12 to first order in . This corresponds to the final Zeeman order: Hence, the protocol in figure 1 enhances the Zeeman order by a factor of 7/6, relative to thermal equilibrium.
B. Cooling cycles 1. Pumping of singlet order A greater cooling effect is established by applying a repeated sequence of cooling cycles to the spin system. Several possible cooling cycles exist: the most effective one we have found may be written as follows (in left-toright chronological order): The symbol τ denotes a relaxation delay, while π ijk denotes a unitary transformation with the following properties: where {ϕ, ϕ , ϕ } are arbitrary phase factors. The cyclic permutations of three spin state populations may be represented by the following matrices: The cooling cycle in equation 17 involves alternating cyclic permutation operations on the states {|1 , |2 , |4 } in the opposite sense, separated by relaxation delays τ , which in the idealized case implement triplet thermal reset operations, as described by the matrix in equation 14.
The effect of the cycle C on the spin-state populations, in the case of idealized triplet thermal reset operations, may be expressed by the matrix to first order in .
If the cycle C is applied repeatedly to the system, a steady-state population vector is eventually established. This steady-state population vector p ss is defined by the condition C · p ss = p ss (21) This implies that p ss lies in the nullspace of C − 1, where 1 is the unit matrix. Linear algebra analysis leads to the following vector of steady-state populations: to first order in . The singlet order in the steady state is given by while the Zeeman order in the steady state is given by Hence the cooling cycles do not immediately enhance the Zeeman order.

Entropy
The cooling cycles reduce the overall entropy of the spin system. A useful quantity is the von Neumann entropy 49 , defined as follows: For small values of spin polarization, the population vector in equation 4 corresponds to the following von Neumann entropy: The first term is the entropy of a completely disordered four-level system. The last term represents the reduction in entropy achieved by thermal contact with the finitetemperature environment.
In the steady state established by repeated pumping cycles, the von Neumann entropy is given by applying equation 25 to the population vector in equation 22. This leads to the following expression: In ideal circumstances, the reduction in entropy achieved by the cooling cycles is therefore a factor of 5/2 more than that achieved by thermal equilibration.

Zeeman order enhancement
The entropy reduction of the spin system may be cashed in as an enhancement of the Zeeman magnetization by applying an additional π 12 permutation operation. This leads to the following population vector: which exhibits the following Zeeman order: Hence, repeated cooling cycles C followed by a π 12 permutation operation lead to an enhancement of the Zeeman order with respect to thermal equilibrium. A comparison with equation 11 indicates that the enhancement factor is 3/2 in the ideal case. In the language of spin thermodynamics 50 an equivalent amount of Zeeman order would be generated if the spin temperature were reduced by a factor of 2/3. However, it should be noted that the population distribution in equation 28 does not conform to a Boltzmann distribution under a Zeeman Hamiltonian and hence cannot strictly be described by a spin temperature.
Note that the enhancement in Zeeman order applies to both sets of spins at the same time, unlike the nuclear Overhauser effect, which enhances the Zeeman order of one set of spins while destroying the order of a different set.
It is shown in appendix D that the cycle C implements a "greedy" algorithmic cooling protocol, as defined in reference 8 .

C. Liouville space analysis
An alternative insight into the algorithmic cooling protocol is achieved by using a geometric interpretation in the Liouville space of orthogonal spin operators 51,52 . In the following discussion, superoperators are denoted by hats (ˆ), while operators go bare-headed and are recognized by context.
Express the density operator as a sum of orthogonal operator terms: where the dimension of Liouville space is N L = N 2 H , ρ q are coefficients and the operators {Q q } are orthonormal: The Liouville bracket is defined as follows 52 Bounds on unitary evolution. The projection of the density operator onto the subspace spanned by the singlet order operator (horizontal axis) and Zeeman order operator (vertical axis) is represented as a point on the plane with coordinates (ρS, ρZ ), defined by equations 34 and 35. The thermal equilibrium density operator is represented by the point (0, ρ eq Z ) (blue dot). The entropy bound on unitary evolution corresponds to a circle of radius ρ eq Z = /(2 √ 2) (dashed line). The Sørensen bound is represented by a non-equilateral hexagon (yellow shape), derived in appendix A. Unitary evolution, starting from thermal equilibrium, leads to a trajectory which is bounded by the hexagon (blue line).
Define the first two operators in the orthogonal set as follows: The projection of a spin density operator ρ onto the {Q S , Q Z } Liouville subspace may be represented as a point with coordinates (ρ S , ρ Z ), defined as follows: As indicated above, the coefficients ρ S and ρ Z are proportional to the singlet order S and Zeeman order Z of the spin ensemble, respectively. For example, the thermal equilibrium density operator corresponds to a point (0, ρ eq Z ), where the equilibrium Zeeman coordinate is given by The thermal equilibrium state is indicated by the blue dot in figure 2.
Trajectory of the spin density operator in the subspace spanned by the singlet order (horizontal axis) and Zeeman order (vertical axis), starting from thermal equilibrium (blue dot). Three repetitions of the cycle C (equation 17) lead to the steady-state density operator indicated by the red dot, generating singlet order well beyond the unitary bound (right edge of the yellow polygon). Enhanced magnetization is generated by applying a π12 permutation to the steadystate density operator. This leads to a final density operator indicated by the green dot, clearly exceeding the thermal equilibrium magnetization. The displayed trajectory is simulated for the following parameters: T1 = 7.3 s, TS = 240 s, and assumes ideal permutation operations. All trajectories were calculated using SpinDynamica 53 .

Unitary Bounds
Coherent evolution of the spin ensemble leads to a unitary transformation preserving the Frobenius norm ρ = r,s | r| ρ |s | 2 . Hence coherent evolution starting from the thermal equilibrium position (0, ρ eq Z ) leads to a point (ρ S , ρ Z ) within a circle of radius |ρ eq Z |, called the entropy bound 54 (dashed circle in figure 2).
Sørensen 54 showed that the bounds on unitary evolution are in general more restrictive than the entropy bound. The Sørensen bounds on unitary evolution, starting from thermal equilibrium, may be written where Q is a normalized Hermitian operator, Q Q = 1, and Λ(Q) is a vector whose elements are the eigenvalues of Q, arranged in ascending order. The regions of Liouville space which are excluded by equation 37 may be derived by a graphical construction 55 . As shown in Appendix A, the bounds on unitary evolution restrict the point (ρ S , ρ Z ) to the interior of a non-equilateral hexagon, shown in yellow in figure 2. Unitary propagation of the thermal equilibrium state follows a trajectory which lies entirely within the yellow hexagon (or touches its boundary).
For any unitary transformation, the two-way conversion amplitude is bounded by a U MSM ≤ a U,max MSM , where the theoretical maximum is given by and the subscript may be read as "magnetization-tosinglet-order-to-magnetization". Hence only ∼ 66.7% of the nuclear magnetization may be passed through singlet order and back again. This theoretical limit applies to any coherent manipulation of the spin system, and has been approached closely in practice by using robust adiabatic transformations 42,43,45,47 . Figure 3 shows numerical simulations of the density operator trajectory under repeated cooling cycles C (equation 17, starting from thermal equilibrium (blue dot). The blue arrows show the trajectory of the density operator under three repetitions of C. After an initial phase of rather irregular motion, the evolution settles into a repeated "bow tie" loop. The upwards-pointing arrows represent the relaxation during the τ intervals, while the cyclic permutations π 124 and π 142 provide the diagonal lines which form the "wings" of the "bow tie".

Pumping of Singlet Order
The red circle in figure 3 indicates the steady-state position of the spin density operator, which returns exactly to the same position after application of a further C cycle. The steady-state point has coordinates (ρ ss S , ρ ss Z ). Figure 3 shows that this point is essentially reached after three repetitions of C, starting from thermal equilibrium. The application of further cooling cycles serve no purpose.
The steady-state point indicated by the red circle lies not only outside the yellow hexagon, but even beyond the dashed circle in figure 3. This indicates that three cooling cycles decrease the entropy of the system beyond the Frobenius bound.
The same steady-state loop is attained by the application of repeated cycles C to any initial state.
Note that the steady state loop never corresponds to an enhancement of Zeeman order beyond the thermal equilibrium value. The application of repeated cycles C pumps the singlet order beyond the unitary bound, but not the Zeeman order. The behaviour of the steady state singlet order may be investigated in detail with the help of SpinDynamica software 53 . Define the relaxation rate constants as R 1 = T −1 1 and R S = T −1 S , with T 1 and T S being the relaxation time constants for Zeeman magnetization and singlet order respectively. In the cases of interest here, singlet order decays much more slowly than longitudinal magnetization, R S R 1 , and T S T 1 . An analytical expression for the steady-state point (ρ ss S , ρ ss Z ) under the cooling cycle C may be obtained by making suitable approximations and assumptions (see Appendix B). This leads to the following expression for the steady-state singlet order: where e 1 = e R1τ /2 , e S = e R S τ /2 , and The dependence of the steady-state singlet order on the ratio of relaxation times T S /T 1 and the cooling cycle intervals τ is shown by the contour plot in figure 4. The cooling effect is most effective for large values of T S /T 1 and intervals τ which are long compared to T 1 but short relative to T S . The limit of equation 41 for ideal condi- In the limiting case, singlet order is enhanced by a factor of 3/2 relative to that achievable by unitary transformation of a thermal equilibrium state (equation 38).

Enhanced Zeeman Order
Although the steady state loop contains points of enhanced singlet order, the Zeeman order is not enhanced. Zeeman order enhancement may be achieved by applying a further π 12 permutation to the steady-state system. The trajectory of the density operator under the entire sequence is shown in figure 3, with the transformation under the terminal π 12 permutation shown by the green arrow. The final state is shown by the green dot. Under ideal conditions, the thermal equilibrium magnetization is enhanced by a factor of 3/2, with respect to thermal equilibrium.
In the language of spin thermodynamics 50 , this quantity of Zeeman order corresponds to the cooling of the nuclear spin system to a temperature equal to 2/3 of the sample temperature.

III. MATERIALS AND METHODS
Experimental demonstration of this algorithmic cooling procedure requires (1) a suitable molecular system displaying near magnetic equivalence and a large relaxation time ratio T S /T 1 , and (2) a robust and high-fidelity experimental procedure for implementing the population permutations.

A. Sample
A suitable molecular system is given by the 13 C 2labelled naphthalene derivative shown in Figure 5 and referred to hereafter as 13 C 2 -I. This system has been shown to support long-lived 13 C 2 singlet order with a decay time constant T S exceeding 1 hour in low magnetic field 38 . The experiments described below were performed at a high magnetic field of 16.4 T on a sample of 13 C 2 -I dissolved in acetone-d 6 . Under these conditions, the relaxation time constants for the 13 C pair were determined to be T 1 = 7.36 ± 0.02 s and T S = 214 ± 5 s respectively, corresponding to a relaxation time ratio of T S /T 1 = 29.1 ± 0.7.
The 13 C NMR spectrum of 13 C 2 -I in acetone-d 6 at a field of 16.4 T is shown in Figure 5. The spin system is well-described by a 13 C-13 C scalar coupling J = 54.141 ± 0.001 Hz and a small chemical shift difference ∆δ = 0.057 ± 0.001 ppm, due to the deliberate chemical asymmetry of 13 C 2 -I. The difference in the chemically shifted Larmor frequencies for the two 13 C sites in the 16.4 T magnetic field is ω ∆ = ω 0 ∆δ = 2π × 10.01 Hz. Since |2πJ| |ω ∆ |, the 13 C pair comprises a nearequivalent AB spin system giving rise to a four-peak pattern with two strong inner peaks and two weak outer peaks, as shown in Figure 5. The splitting between the strong inner peaks is given by ∼ ω 2 ∆ /(8π 2 J) 0.9 Hz. The outer peaks are very weak but may be enhanced by specialised pulse techniques 56 .
The work described here used a solution of 13 C 2 -I in acetone-d 6 sealed into a micro-cell of length 15 mm, inserted into a standard 5 mm NMR tube and positioned about 1 cm above the bottom of the tube. This procedure avoids convection 41 and provides high homogeneity for both the static field and the radiofrequency field. The experimental demonstrations use APSOC (Adiabatic Passage Spin Order Conversion) pulses [42][43][44][45] . These are amplitude-modulated radiofrequency pulses of carefully optimized shape, applied at a frequency slightly displaced from the centre of the spectral pattern. The shape of the APSOC pulse is optimized for the given spin system parameters using previously described algorithms [42][43][44][45] , and is specified in appendix C. The APSOC pulses are applied with phase φ = 0 but off-resonance by ±35 Hz, relative to the centre of the NMR spectrum, taking into account the sign of the Larmor frequency 57 , as depicted in Figure 6(a,b).
As explained in refs. [42][43][44][45] , an APSOC pulse induces a high-fidelity population exchange between the singlet state |1 and one of the outer triplet states |2 or |4 . However, care is needed, since these states are not defined in the laboratory frame, but in a frame rotated by ±π/2 about the rotating-frame y-axis. The propagator for the two varieties of APSOC pulse may be written as follows: Here I rs y is a single-transition operator 62,63 and the operators Φ (±) induce phase shifts of the individual spin states: where φ (±) r are uncontrolled phase factors, which are unimportant in the experiments described here, since only the spin state populations are involved.

Cyclic Permutation Sequences
The pumping cycle C (equation 17) requires cyclic permutation operations π 124 and π 142 for the spin state populations, defined by the transformations in equation 18. Each cyclic permutation is realised by applying an AP-SOC pulse immediately followed by a π/2 pulse of appropriate phase.
For example, consider an APSOC(-) pulse, immediately followed by a strong π/2 pulse with phase π/2. This combination has a propagator given through equation 45 by: The APSOC-pulse combination transforms the four spin states as follows: This shows that the combination of an APSOC(-) pulse with an appended (π/2) y pulse fulfils the properties of the cyclic permutation operator π 124 defined in equation 18. Similarly, an APSOC(+) pulse with an appended (π/2) −y pulse has the functionality of the cyclic permutation operator π 142 .
In practice, high fidelity for the π/2 rotations is ensured by using composite pulses 24,58-61 of the form 180 ±30 90 ±150 , where the notation β φ indicates the pulse flip angle β and the phase φ, both in units of degrees (see figure 6(a,b)). These sequences are derived by taking the second halves of composite π-pulses of the form 90 ±150 360 ±30 90 ±150 , and are very well-compensated for errors in the radio-frequency field amplitude 59 . The phases of the composite pulses take into the sign of the Larmor frequency and the radio-frequency mixing scheme 57,64 . The APSOC(∓)-180 ±30 90 ±150 sequences provide high-fidelity cyclic permutations of the spin-state populations, well-compensated for radiofrequency field errors.

A. Pumping of Singlet Order
The build-up of singlet order under repeated pumping cycles C was assessed by the protocols shown in figure 7. A fine-grained study is achieved by varying the number of permutations n P , instead of the number of cycles n C . Since each cycle C contains two cyclic permutations, even values of n P are accessed by completing n P /2 cycles (figure 7(a)). Odd values of n P are accessed by appending an additional delay and permutation element to (n P − 1)/2 cycles ( figure 7(b)).
The intervals τ are chosen such that T 1 τ T S . Under these circumstances, each relaxation interval τ provides an acceptable approximation to the triplet thermal reset Θ reset T (equation 14), allowing the triplet populations to approach thermal equilibrium with the environment, while the singlet population is left approximately unchanged. In practice the value τ = 28.0 s was used.
The pumped singlet order was allowed to evolve for an optional interval τ ev , followed by a T 00 signal filtration element, as specified in figure 6(c). The T 00 filter suppresses all signal components that do not pass through singlet order 37 . A final APSOC pulse efficiently converts the singlet order to transverse magnetization, which induces an NMR signal in the subsequent interval of free precession. The integrated amplitude of the NMR signal is therefore a reliable measure of the singlet order generated by the repeated cycles C. Fig. 8 shows the singlet-filtered signal amplitude, normalized with respect to the thermal equilibrium signal induced by a π/2 pulse, as a function of the number of permutations n P , using a relaxation interval τ = 28.0 s, and zero evolution interval, τ ev = 0. The singlet-filtered signal amplitude builds up rapidly, reaching a steadystate of around 0.85 times the thermal equilibrium signal after approximately 6 permutations. Further permutations do not increase the observed singlet order.
The unitary bound on passing thermal equilibrium magnetization through singlet order is given by the factor 2/3 (equation 40). This two-way magnetization-tosinglet-order-to-magnetization bound is shown by the red dashed line in Fig. 8. The results show that the algorithmic cooling cycles successfully generated singlet spin order with an amplitude 1.27 times the unitary limit. The loss with respect to the ideal enhancement factor of 3/2 (equation 43) may be attributed to non-ideality of the triplet thermal reset operations, due to the finite ratio of T S to T 1 , and experimental infidelities in the population permutations.
The singlet pumping efficiency depends on the relaxation delays τ , as shown in Figure 9(a). If the delays are too short, there is not enough time to establish fresh polarization in the triplet manifold before the next permutation; if the delays are too long, singlet-triplet transitions allow the loss of singlet order during the delays. In the case considered, the generation of singlet order is optimised for delays of duration τ = 28.0 s.
The generation of pumped singlet order was verified by allowing an evolution time for the singlet order to decay, before conversion into transverse magnetization for detection. As shown in Figure 9(b), the singlet order decays with the characteristically long time constant of T S = 209 ± 6 s, which greatly exceeds the spin-lattice relaxation time constant of T 1 = 7.36 ± 0.02 s. This behaviour verifies that the state generated by the pumping protocol does indeed correspond to long-lived nuclear singlet order.

B. Enhancement of Magnetization
The low-entropy state achieved by pumping the singlet order may be used to enhance the nuclear magnetization. An enhanced NMR signal may be generated from the steady state by an additional relaxation delay τ followed by an APSOC pulse ( figure 7(c)). The optimal value of τ was found empirically to be 18.0 s, which is smaller than the optimal pumping interval τ = 28.0 s. Figure 10 shows experimental results demonstrating the enhancement of the nuclear spin magnetization with respect to its thermal equilibrium value. Figure 10(b) shows a 13 C NMR signal of 13 C 2 -I, enhanced by 21% with respect to the thermal equilibrium spectrum in figure 10(a). The same enhancement in nuclear magnetization would be accomplished by a reduction in spin temperature by a factor of 0.83. This result proves that it is possible to enhance the nuclear magnetization by algorithmic cooling, even on a strongly coupled spin system which lacks any distinction between fast-relaxing and slowly relaxing spins.  Fig. 6. Singlet order is pumped by nC = nP /2 repetitions of the cycle C (equation 17), followed by an optional evolution interval τev, a T00 filter (figure 6(c)) and an APSOC(+) pulse for generation of the NMR signal. (b) Protocol for an odd number of permutations nP . The number of cycles is nC = (nP − 1)/2, and an additional delay and permutation is appended; (c) Protocol for detecting enhanced nuclear magnetization for an even number of permutations nP . After pumping the singlet order by nP /2 cycles, the system is allowed to relax for a delay τ , followed by an APSOC(+) pulse to generate observable transverse magnetization. In all cases, a singlet-order-destruction (SOD) sequence was applied after signal detection, in order to allow faster repetition of the experiments 65 (not shown).

V. CONCLUSIONS
Since the achieved enhancement of nuclear magnetization is modest, we do not claim that this method is in serious competition with true hyperpolarization techniques such as dynamic nuclear polarization 66 . Nevertheless the phenomenon is of interest, since it transcends several paradigms which have played a prominent role in the theory and implementation of algorithmic cooling. The algorithmic cooling of nuclear spin systems has generally been framed in the language of NMR quantum computation 23 , in which a prevailing assumption is that the individual nuclear spins may be treated, to a good approximation, as independent qubits which are indepen- dently addressable by suitable radiofrequency fields [15][16][17][18][19][20][21][22] . The current work has shown that (i) algorithmic cooling may be achieved on systems which lack independently addressable qubits; (ii) algorithmic cooling may exploit entangled superpositions of the qubit states (the singlet state in the current case); (iii) algorithmic cooling does not necessarily require qubits with distinct relaxation behaviour; (iv) algorithmic cooling may be used to enhance modes of spin order which do not correspond to spin polarization. The last point is salient: The primary target of the cooling algorithm in the current work is not the establishment of a low effective spin temperature in the spin ensemble but the enhancement of non-magnetic singlet order. This is an anticorrelated state of the spin ensemble which does not correspond to macroscopic spin polarization and which cannot be described by a spin temperature. Nevertheless, as shown here, singlet order may be pumped by an algorithmic procedure to a level which significantly exceeds that which is accessible by a unitary transformation of the thermal equilibrium state. We also show that pumped singlet order may be converted in a subsequent step to nuclear magnetization that is substantially enhanced with respect to its thermal equilibrium value, and which does therefore correspond to a reduced spin temperature.
This work exploits a very simple system: an ensemble of strongly coupled spin-1/2 pairs. It is not currently known how the phenomenon described here translates to more complex quantum systems, including those which use substrates other than nuclear spins. The fundamental limits on algorithmic cooling of the type described here, in which the relevant system eigenstates are not discrete qubit states but entangled superposition states, are not yet known. It is also unknown whether the theoretical bounds on algorithmic cooling in discrete qubit systems 9-13 apply here too. We expect the current report to stimulate interest, and further development, in the NMR and quantum information communities, and elsewhere.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.
which may by specified by the notation (r, θ), where the polar coordinates are r = {x 2 + y 2 } 1/2 and θ = arctan(x/y). The thermal equilibrium density operator corresponds to a point with polar coordinates (ρ eq Z , 0) (filled blue disk in figure 11). Unitary evolution of the density operator, starting from the thermal equilibrium point (ρ eq Z , 0), leads to a new point in the Liouville subspace, with polar coordinates (r, θ). We now derive the constraints on the physically accessible points (r, θ), imposed by the Sørensen bound (equation 37), by adapting the reasoning in reference 55 .
Define the operator Q θ as a superposition of the normalized Zeeman order operator Q Z and singlet order operator Q S with mixing angle θ: This corresponds to a point with polar coordinates (1, θ). The eigenvalues of Q θ are given by As the angle θ is varied, the eigenvalues in equation A3 cross. Eigenvalue degeneracies occur at the crossing points, which are found at angles θ given by θ ∈ {0, Θ, π/2, π − Θ, π, π + Θ, 3π/2, 2π − Θ} (A4) where Θ = arctan From the Sørensen bound (equation 37), the projection of ρ onto Q θ is subject to the following constraint: Here Λ(Q) represents a vector containing the eigenvalues of the operator Q, arranged in ascending order. The set of points (r B (θ), θ), ∀ 0 ≤ θ < 2π, defines a boundary line B, shown by the black curve in figure 11. The curve B displays discontinuities at the eigenvalue crossings, which occur for θ values given in equation A4, with the exception of θ ∈ {0, π}.
All points (r, θ) which are outside the curve B are physically inaccessible by unitary transformation of the thermal equilibrium state. However, not all points within the curve B are physically accessible. Consider, for example, the boundary point (r B (Θ), Θ), which is shown by the filled orange disk in figure 11. The orange dashed line passes through this point and the origin. The point (r, θ) is only physically accessible if its projection onto the orange dashed line is less than, or equal to, r B (Θ). This excludes all points (r, θ) in the half-plane shown by the orange-shaded region in figure 11. Part of this excluded region lies within the boundary line B.
Since the curve B has re-entrant projections at the six θ angles {Θ, π/2, π−Θ, π+Θ, 3π/2, 2π−Θ}, the constraints imposed by these six values of θ delineate the boundaries of the allowed region of Liouville space. Excluding all forbidden regions leaves the central hexagonal region in . This region is physically inaccessible, starting from thermal equilibrium. Note that there is a region of Liouville subspace which lies within B but which is nevertheless physically inaccessible. natural choice; however, pure DD relaxation would lead to an infinite value of T S . A realistic relaxation model would need to be relatively complex, including DD, chemical shift anisotropy, and external random field mechanisms 28,34 . For simplicity, a partially-correlated randomfield mechanism was used in the current work.
The analysis of algorithmic cooling requires the use of relaxation superoperators which take into account the finite temperature of the environment. Several methods have been proposed for the construction of such superoperators 3,52,67,68 . In this work, we use a Lindbladian formulation, which is rigorously consistent with open quantum systems theory 49 , avoiding the pitfalls of the other methods 3 .
The thermalized relaxation superoperator is given in the extreme-narrowing limit by 3 1µ ] (B1) Here the two spins-1/2 are denoted I j and I k , ω ran is the root-mean-square amplitude of the fluctuating random fields, expressed in frequency units, and τ c is the correlation time of the random fields. The amplitude of the random fields is assumed to be the same on both spins, and their correlation coefficient is −1 ≤ κ ≤ 1. The symbol T (j) λ,µ represents a component of the rankλ irreducible spherical tensor operator for spin I j . The Lindbladian dissipatorD[A, B] is defined as follows 3 : where the "bullet" symbol indicates the position of the operand, for example: The relaxation superoperator in equation B1 leads to the following expressions for the characteristic relaxation rate constants: The matrix elements given above are conveniently evaluated using SpinDynamica software 53 .
In the high-temperature and high-field approximations, the superoperator in equation B1 leads to the following matrix for the thermalized transition probabilities per unit time, between the four spin eigenstates: The evolution of spin populations over the cycle C in equation 17 is given by the propagation matrix C = π 142 exp{W θ τ }π 124 exp{W θ τ } (B6) and the permutation matrices π 124 and π 142 are given in equation 19. The steady-state population vector p ss of C is found by determining the null space of the matrix C − 1, followed by normalisation so that the sum of all elements is equal to 1. The full expression for the steady-state vector is rather lengthy, but leads after some manipulation to the steady-state singlet order given by equation 41.

Appendix C: APSOC pulses
The APSOC shape used in this work is wellapproximated by a polynomial of the form The maximum value of the rf field amplitude, expressed as a nutation frequency, is ω max APSOC = 2π × 181 Hz, and the duration of the APSOC pulse is given by T APSOC = 0.36 s. The shape of the APSOC pulse is specified by the following polynomial coefficients: In this appendix we show that the cooling cycle C implements a "greedy" optimisation protocol that aims to maximise the total order of the system at each iteration step 8 . The total order of the system is quantified in terms of its purity, which is closely related to the entropy of the system 8,49 . The von Neumann entropy (equation 25) may be approximated to lowest order as follows: The quantity Tr(ρ 2 ) represents the purity (γ) of the system 49 γ = Tr(ρ 2 ), and has been utilised previously in the study of algorithmic cooling protocols 8,11,13 . In contrast to the entropy S, the purity γ is maximised for pure states and minimised for fully mixed states. A pure state is characterised by a purity of γ = 1, whereas a fully mixed state is characterised by a purity of γ = N −1 H . The iterative application of the cooling cycle C leads to the following recurrence relation for the density matrix of the system ρ n = C(ρ n−1 ), ρ 0 = ρ eq , with the corresponding set of purities γ n = Tr(ρ 2 n ), γ 0 = Tr(ρ 2 eq ) = γ eq .
The equilibrium purity γ eq for a weakly polarised spin-1/2 pair at moderate magnetic fields is given by The connection between a greedy protocol and the cooling strategy C may be established by considering an arbitrary two-step cooling cycleC C = (π a − τ − π b − τ ). (D6) The population vector p 1 thus fulfils the same assumptions that have been used to determine the optimal strategy for the 0'th step. Since there is nothing special about the first step, it follows that for a general step j one has δ j 4 > 0 (D20) and (δ j This implies that at each step of the cooling process the application of the cyclẽ results in a local optimisation of the purity. The cycle given by equation D22 is slightly different to the cycle C given in equation 17. However, in the limit of infinitely long singlet and infinitely short triplet relaxation times the two cycles coincide with each other. Inspection of the triplet thermal reset Θ reset T reveals that it is right invariant under elements of S 3 with π j ∈ S 3 acting on the indices {2, 3, 4}. It is therefore permissible to replace π 14 by π 24 π 14 = π 124 and π 12 by π 24 π 12 = π 142 . With this choice the two cooling cycles are identical C =C and the cycle C does indeed follow a "greedy" protocol. The situation for finite relaxation times τ may be analysed in a similar fashion. Here, the permutation decompositions π 24 π 14 = π 124 and π 24 π 12 = π 142 may be given a physical interpretation. Application of the permutations π 14 and π 12 , respectively, increases the purity of the system. The permutation π 24 , on the other hand, exchanges triplet populations in order to facilitate quicker thermal equilibration within the triplet manifold. For the idealised triplet thermal reset this permutation is superfluous as the reset operation is considered to restore triplet thermal equilibrium "instantaneously".