Spin distillation cooling of ultracold Bose gases

We study the spin distillation of spinor gases of bosonic atoms and find two different mechanisms in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^{52}$$\end{document}52Cr and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{23}$$\end{document}23Na atoms, both of which can cool effectively. The first mechanism involves dipolar scattering into initially unoccupied spin states and cools only above a threshold magnetic field. The second proceeds via equilibrium relaxation of the thermal cloud into empty spin states, reducing its proportion in the initial component. It cools only below a threshold magnetic field. The technique was initially demonstrated experimentally for a chromium dipolar gas (Naylor et al. in Phys Rev Lett 115:243002, 2015), whereas here we develop the concept further and provide an in-depth understanding of the required physics and limitations involved. Through numerical simulations, we reveal the mechanisms involved and demonstrate that the spin distillation cycle can be repeated several times, each time resulting in a significant additional reduction of the thermal atom fraction. Threshold values of magnetic field and predictions for the achievable temperature are also identified.

www.nature.com/scientificreports/ some detail the physical underpinnings and conditions under which the spin distillation cooling can take place. We find that the physics of the cooling differs between chromium and sodium gases, as do the fundamental limitations on the processes in question. The paper is organized as follows: We first introduce the model in the "Model" section, including specific aspects of the dipolar 52 Cr and contact interacting 23 Na systems. Results are then presented in the "Chromium: condensate assisted dipolar scattering to higher spins" and "Sodium: cooling thanks to redistribution of the thermal cloud" sections for chromium and sodium, respectively. Each demonstrates one of the two mechanisms we have identified. General remarks and a summary are given in the "Conclusions" section.

Model
We consider two systems: 52 Cr in the s = 3 state and 23 Na in the F = 1 hyperfine state. The gases are at temperatures that are a sizable fraction of the critical temperature T c , so they are modeled by a complex classical field, which is the most tractable and accurate model under these conditions. Several complementary techniques have been developed [24][25][26][27][28][29][30][31][32] , reviewed in 33-35 . Initial thermal states. The cooling simulations start with thermal clouds of atoms in a 3d harmonic trap polarized in a single spin state, which is m s = −3 for for chromium and m F = 0 for sodium. Therefore a single-component stochastic Gross-Pitaevskii equation (SGPE) 28,36 provides a convenient route to directly obtain thermal ensembles at a chosen temperature and particle number N 35,37 . One evolves the complex field �(r) in technical time τ according to the equation until a stable fluctuation of observables in time around a mean is obtained-i.e. the stationary ensemble. Eqn.
(1) is the classical field equation for a field in thermal and diffusive contact with a reservoir at temperature T and chemical potential µ . Here m is the atomic mass, ω j the trap frequencies in directions j = {x, y, z} (with coordinates r j ), and γ is a dimensionless coupling strength to the reservoir. For our system the trap frequencies were ω j /2π = (250, 300, 215)Hz, in all cases. The η is a complex white noise field with mean zero and variances �η * (r, τ )η(r ′ , τ ′ )� = δ (3) (r − r ′ )δ(τ − τ ′ ) , �η(r, τ )η(r ′ , τ ′ )� = 0 . Interaction energy is captured by the functional H I (�) , which is with contact interactions in the 1st term, and dipolar interactions (only in chromium) in the second. For chromium, H d (�) = m s [H d (�(r))] 11 with the latter given in the "Methods" section and g = g 6 as described in the "Chromium gas in the s = 3 state" section; For sodium, H d = 0 and g = c 0 as described in the "Sodium gas in the F = 1 hyperfine state" section. The magnetic field and chromium dipoles are oriented along the z axis.
The standard SGPE approach models the quantum degenerate gas using a classical field that describes the low energy/high density part of the system on a basis set spanning modes with occupations O (1) , and a reservoir that describes all remaining high energy modes with occupations O (1) . A common practice is to restrict the low energy subspace for to all plane wave modes below a certain momentum cutoff k cut . We also follow this practice, setting the numerical lattice spacing �x = π/k cut using k cut = 0.78 √ 2πmk B T according to the optimal multi-observable choice in 3d from 38 . The box size was around 3.2R TF = (3.2/ω j ) √ 2µ/m , and was chosen to provide the experimentally measured condensate fraction. Thus, the mode space is chosen according to a balance of various observables. The simulations start from the neutral initial conditions of the vacuum state �(r, 0) = 0 , and used γ = 0.5 , an efficient value for finding equilibrium ensembles (they do not depend on γ 36 ). µ is chosen to obtain the desired ensemble average particle number. Dependence on the two control parameters T and µ strongly suggests that this approach should be equivalent to others based on the grand canonical statistical ensemble (as, for example, in 31 ). Indeed, it can be demonstrated that the stochastic field fulfills the fluctuation-dissipation relation 39 which is an indication of required thermal properties. Differences between grand canonical and canonical ensembles are negligible in the interacting gas in the regime of interest [40][41][42] .
Samples of the stabilized field �(r, τ ) in the stationary ensemble are used as the initial states for further evolution. This subsequent evolution, during which cooling cycles are carried out is best performed without the high energy bath, because of its inherently non-equilibrium nature and the long time scales involved. It proceeds then via a plain Gross-Pitaevskii equation (GPE) 33,34 . Chromium gas in the s = 3 state. The first system we focus on is a partially condensed gas of 52 Cr atoms in the s = 3 state, as per the initial experiment 18 . In the classical-field model the cloud is described by the sevencomponent spinor wave function ψ(r) = (ψ 3 (r), ψ 2 (r), ψ 1 (r), ψ 0 (r), ψ −1 (r), ψ −2 (r), ψ −3 (r)) T with one component for each value of the spin projection m s ∈ {−3, −2, . . . , 2, 3} along the direction of the applied magnetic field B . It evolves according to the multicomponent Gross-Pitaevskii equation 43 : www.nature.com/scientificreports/ where the energy functional consists of three parts. The first one represents the single-particle term which includes the kinetic, potential V trap = 1 2 m j ω j r j , as well as Zeeman energies, in which µ = g L µ B s/ . Since the magnetic field direction is constant in the z direction, −µ · B = −g L µ B m s |B| ≈ 2µ B m s B.
The second term originates from the contact interactions. H c is a 7 × 7 matrix in spinor components. The matrix elements are evaluated in a two-atom basis consisting of |s, m 1 �|s, m 2 � states, where m 1 , m 2 ∈ {−3, 2, . . . , 2, 3} , and |s, m s � is the state of a single atom with spin s and projection m s . The matrix elements are given by Here, the symbol �s, m 1 ; s, m 2 |SM� is a Clebsch-Gordan coefficient and |SM� is the state of a pair of atoms with a total spin S and spin projection M. Only S = 0, 2, 4, 6 channels are allowed, M = m 1 + m 2 = m ′ 1 + m ′ 2 and the interaction strengths characterizing the colliding atoms with a total spin S are given as g S = 4π 2 a S /m . We take the scattering lengths to be: a 0 = 91 , a 2 = −7 , a 4 = 63 , and a 6 = 102 in units of the Bohr radius 44,45 . The last term in Eq. (1) describes the dipolar interactions. Since the spin projection of a pair of atoms colliding via dipolar forces can change at most by 2 and the spin projection of a single atom changes maximally by 1, the matrix H d becomes tridiagonal 43  Evolution starts with the equilibrated stochastic field �(r) obtained with Eq. (1) placed in the m s = −3 state as ψ −3 (r) , and vacuum ψ m s (r) = 0 in the other spin states m s = −2, . . . , 3 . Subsequent evolution follows Eq. (4).
Sodium gas in the F = 1 hyperfine state. The second system we focus on is a thermal gas of 23 Na atoms. Due to their low magnetic moment, their dipolar interactions will be not considered. This is a spin-1 system described by a spinor wave function with three components ψ = (ψ 1 , ψ 0 , ψ −1 ) T of spin projection m F = 0, ±1 along the externally applied magnetic field.
The model Hamiltonian of the system is simpler than in the case of 52 Cr and contains no dipolar term. Therefore in the evolution equation (4), H d = 0 . The contact part for the jth component, coming from (6) with S = 0, 2 , can be written in which n(r) = m F n m F (r) = m F |ψ m F (r)| 2 is the local total atom density and F = (ψ † F x ψ, ψ † F y ψ, ψ † F z ψ) is the local spin density. The F x,y,z are the spin-1 matrices. The spin-independent and spin-dependent interaction coefficients are c 0 = 4π 2 (2a 2 + a 0 )/3m and c 2 = 4π 2 (a 2 − a 0 )/3m , where a S is the s−wave scattering length for colliding atoms with total spin S 46,47 . We take a 0 = 50 and a 2 = 55 in units of the Bohr radius a B , consistent with most determinations [48][49][50] .
Since the Hamiltonian of an F = 1 spinor gas is invariant with respect to a rotation of the spin vector around the direction of the magnetic field, we consider only the quadratic Zeeman effect 51 . The magnetic part of the Hamiltonian is H QZE = −q dr n 0 (r), where a constant offset has been dropped. n 0 is the atom density in the m F = 0 component, and q = α q |B| 2 is the Zeeman energy. For sufficiently low magnetic fields, α q /2π ≈ 277 Hz/G 252 . Then Similarly to the chromium case, the state with the lowest effective magnetic energy (here m F = 0 ) is initially populated by the stochastic field �(r) generated by (1). The other components start with a tiny seed of ψ ±1 (r) = 10 −5 �(r) . The dynamics of the system is subsequently described by:

Chromium: condensate assisted dipolar scattering to higher spins
We proceed as in the initial experiment 18 and begin with the chromium gas confined in the harmonic trap, with all atoms in the lowest Zeeman state of the system, m s = −3 . Dipolar collisions allow for a change of the magnetization of the system 45,[53][54][55][56] . In fact, transfer of atoms to m s ≥ −2 components is possible only due to dipolar interactions 43 . Energies available in the cloud can surmount the Zeeman energy difference |g L �m s |µ B B ∼ 2µ B B when the magnetic field magnitude B = |B| is sufficiently low, making possible the population of higher Zeeman states. What one also wants for cooling is for thermal atoms to change spin state, but condensate atoms to remain in m s = −3 . The conditions needed for this are twofold. First that thermal atoms can cross the Zeeman threshold, requiring Second that condensate atoms cannot, requiring that the ground state is ferromagnetic. The threshold magnetic field above which the ground state of the system remains polarized in m s = −3 57 is where n is the atomic density. When the condensate peak density is used for n, condition (11) prevents all condensate atoms from changing spin state. We can safely omit dipolar contributions to condition (11) because their input to condensate energy is minor 58,59 . The ratio of dipolar length ( a dd = µ 2 m/3 2 ) to contact scattering length in the m s = −3 state is ǫ dd = a dd /a 6 = 0.15 ≪ 1.
In each cooling cycle the initial cloud is allowed to evolve in-trap for 155ms. At the end of the cycle, all atoms in the higher m s ≥ −2 spin states are removed. Then a new cycle is begun starting with the remaining atoms in m s = −3 . Condensate fractions of the clouds are estimated by the spatial averaging technique 34,60 . In this case, the cloud is averaged over the z coordinate to construct an approximate one-body density matrix ρ j gives an estimate of condensate occupation in spin state j.
In Fig. 1 three successive cycles of a spin distillation procedure of partially condensed 52 Cr are presented. Initially the cloud contains N ≈ 20000 atoms, with a condensate fraction of 0.55, corresponding to a temperature of about T 0 = 235nK. These conditions are similar to the experiments 18 . In this simulation the magnetic field was kept at the same value B = k B T 0 /2µ B through all cycles, chosen as per the temperature condition in (10) with the initial temperature. After three cooling steps the condensate fraction increases very strongly up to 0.9. Notably, this simulation demonstrates that successive cooling cycles can give continued improvement. The condition (10) eventually breaks down due to falling T. In fact no further cooling is seen after the third cycle, and with this protocol the final condensate fraction is 91%, with N = 13000 atoms in total. www.nature.com/scientificreports/ It has been conjectured that much lower temperatures can be reached when the magnetic field is adapted to a running value of T 18 . Figure 2 shows the progress when the magnetic field is adapted at the beginning of each cooling cycle to the actual temperature value T(t) according to the rule is the the ideal gas critical temperature. Here, the final condensate fraction is about 0.97 with N = 12850 atoms remaining.
The dynamics is investigated in more detail in Figs. 3 and 4 which plot the evolution of several quantities over a single cooling step. The first of these shows that the number of condensed atoms in the m s = −3 state does not significantly change in time, while thermal cloud atoms migrate to higher m s components (primarily m s = −2 ), which are purely thermal. This is as planned, and in agreement with the conditions (10)- (11). Taking the Thomas-Fermi estimate of the peak density n peak −3 in a fully condensed cloud 61 as the maximum value of n, we estimate that B th = 0.29 mG for our case.
An understanding of the cooling mechanism seen can be summarized as follows: where the c and th superscripts denote condensate and thermal fractions, respectively. Higher spin states are reached through At ultracold temperatures Bose enhancement favors processes that involve the condensate, but all-condensate transitions like those seen in 18,43 are ruled out by the Zeeman energy threshold.   rates, due to contact interactions, are determined by the density n of atomic cloud, the typical velocity v, and the cross section σ through γ = nσ v . For the m s = −2 Zeeman state, at the peak atomic density and at the critical temperature γ ≈ 40s −1 ( 64 ). Hence, a timescale of hundreds of milliseconds is needed to establish equilibrium within the m s = −2 state. On the other hand, atoms in different spin states are not seen to equilibrate with each other, the kinetic energy per thermal atom does not equipartition between the Zeeman states: Fig. 4 right frame. One reason is the low rate of spin changing collisions due to the small spatial overlap of thermal clouds appearing in different spin states. This occurs because higher spin states have low population and the resulting clouds exhibit speckle-like density patterns due to thermal fluctuations. Moreover, the equipartition is not guaranteed to hold in equilibrium for energy gaps of the order of k B T or more. 4. Kinetic energy per thermal atom is higher in more energetic spin states. This can be due to the following properties of dipolar interactions. For one, atoms transferred dipolarly to a higher spin state must carry nonzero angular momentum 62 . The higher the spin state, the more angular momentum needs to be carried. Therefore, the threshold energy consists of both Zeeman and rotational contributions. The energy associated with the rotation is next dissipated into the cloud increasing the kinetic energy present. The growth of the dipolar collision rate with velocity at low magnetic fields 44,63 may also contribute to the highest energy atoms moving further up the spin ladder. Since spin-changing collisions are very low rate processes 64 , the net result is an increase of kinetic energy per thermal atom in successive spin states.
Together, the conditions (10) and (11) set a limit for this cooling approach: k B T 2µ B B th . Substituting the peak density estimate n = µ/g 6 , we obtain In particular, depending on the scattering lengths, this can be very far below µ , which is the usual limiting value for standard evaporative cooling. For our present parameters, condition (12) gives a very low value of 39nK.

Sodium: cooling thanks to redistribution of the thermal cloud
The second case of a sodium gas in the F = 1 hyperfine state, was considered briefly by 18 via a basic uniform model. In the detailed model here, we find that the cooling is also effective, and notably that the process differs significantly from the chromium case. The initial thermal state is prepared in the m F = 0 component, which has lowest quadratic Zeeman energy. Transfer of thermal atoms to the m F = ±1 components is possible via the spin-asymmetric contact interactions (the last terms proportional to c 2 in Eqs. (9)).
Zero-magnetization case. We start our analysis with a case of zero magnetization, i.e. when the populations of m F = ±1 states are equal. Additionally, we choose initial occupations of these states as almost zeroonly a tiny seed is left to allow a transfer to m F = ±1 states. Other starting conditions are possible as well and are mentioned in the next subsection. The results of a simulation are presented in Fig. 5 for Zeeman energy q = 0.011 ω x , about N = 20000 atoms and an initial temperature of about 235nK. Significant cooling is seen in each cooling cycle, via an increase of condensate fraction in m F = 0 , though notably the timescale is longer than for chromium. Cooling slows after about 1 s, which is a good time to remove atoms in the m F = ±1 states and begin a new cycle. After three cycles the m F = 0 population is N ∼ 4500 − 10000 (large Rabi oscilations) with 90% condensate.   Fig. 1 (left) and as in Fig. 2  www.nature.com/scientificreports/ A visible element of the evolution is a rapid oscillation (period ∼ 30ms). This stems from a Rabi oscillation that coherently transfers both condensate and noncondensate atoms between the m F = 0 and m F = ±1 hyperfine states as can be deduced from the top row of Fig. 6. Its period ( ∼ 30ms) corresponds well to predictions based on the energy of the spin changing term c 2 n 0 . The presence of the oscillation indicates coherent exchange between spin states, both in the condensate and the thermal cloud (at least initially). Therefore a different process is dominant than the one seen for chromium.   The obvious condition to transfer thermal atoms from m F = 0 to ±1 via the 2n 0 ⇋ n ±1 mechanism is that the thermal energy is sufficient to overcome the quadratic Zeeman energy, q, i.e. q k B T . However, the magnetic field used here is already practically negligible compared to the thermal energy ( k B T ∼ 1800q in Fig. 5). Instead, we find a limiting factor at much lower magnetic field. Figure 7 shows the first step of cooling for q/ ω x = 0.15 and 0.17. Surprisingly, no cooling is observed for the higher magnetic field, even though the thermal energy there is still enormously greater than the Zeeman energy: k B T/q ≈ 100 . Figure 8 shows the abrupt breakdown of cooling with growing Zeeman energy q in detail. There must be another limitation involved.
The spin mixing terms responsible for transfer from the initial m F = 0 to m F = ±1 have energy of the order of c 2 n 0 . Therefore, the amplitude of any spin mixing process will decay rapidly once the energy difference exceeds this level. In order to take two atoms from m F = 0 and place them in the m F = ±1 states, an energy of q per atom is required, leading to the additional condition for cooling onset. Thresholds of the same order ( q O(c 2 n) ) are commonly seen in studies of the stationary states of low temperature F = 1 anti-ferromagnetic condensates 51,65,66 . In our case, taking the peak density of a 100% condensate, condition (13) is q 0.16 ω x , which matches quite well with the threshold seen in Fig. 8. Further corroborating evidence for this interpretation of the cooling threshold includes the fact that no Rabi oscillations are seen above threshold at q = 0.2 ω x in Fig. 6 (bottom right). Moreover, if one takes a small condensed fraction of 6%, no cooling or Rabi oscillations are observed (Fig. 9, top left), nor transfer to m F = ±1 , despite the same Zeeman energy q = 0.09 ω x that cooled the larger condensate in Fig. 7. This is also consistent and 15 ω x . At the lowest q 0.025 ω x , the 1 s experimentally realistic cycle is too short for stabilisation and the randomness reflects snapshots of different stages of the system's Rabi cycle at 1 s. An example is seen in Fig. 6 www.nature.com/scientificreports/ expected since a small condensate has highly reduced density, and the condition (13) now becomes q 0.03 ω x -which is not met. What is the actual mechanism for cooling, though? It cannot be a threshold Zeeman energy that sorts atoms into high and low energy in different spin states like we saw in chromium, since the atoms can be seen flowing freely to all spin states and back. On the other hand, the population transfer process is not reversible in the long term and the populations and condensate fractions settle to a more or less stationary value. There is another collision process that does not change spin populations, but can rearrange atoms between modes. The spinpreserving and spin-exchanging collisions between m F = 0 and m F = ±1 atoms include contributions such as which can exchange thermal and condensate fractions within a spin component, and degrade the Rabi oscillations, leading to a stationary equilibrated state. Notably, the cooling is arrested when the iψ 0 ∼ (n +1 + n −1 )ψ 0 and iψ ±1 ∼ n 0 ψ ±1 terms responsible for the ψ 0 &ψ ±1 ⇋ ψ 0 &ψ ±1 collisions are removed ( Fig. 9 top right), proving the crucial role of this process.
The process can be followed in Fig. 9 (bottom). Two long timescales are evident: Most of the flow of thermal atoms from m F = 0 to m F = ±1 is completed by about 1 s, and it is also at this time that the Rabi oscillations of the thermal components die out (also Fig. 6, bottom left). The Rabi oscillations between condensates last longer, and there is a second timescale of about 3 s over which the last thermal atoms accumulate in m F = ±1 and the kinetic energy of thermal atoms in all three hyperfine states equilibrates.
The cooling procedure can in fact be understood as a re-equilibration punctuated by periodic removal of the atoms in the unwanted spin states. Once condition (13) is met, the number of accessible thermal degrees of freedom increases threefold. The thermal atoms redistribute, leaving only 1 3 of the original number in m F = 0 . Figures 6 and 9 (bottom left) confirm that the number of condensate atoms in m F = 0 reduces only slightly by about 2000, while the initial 9000 thermal atoms in m F = 0 redistribute approximately equally among the three hyperfine states, leaving only 3000 in m F = 0 . Hence, after m F = ±1 populations are removed, significant cooling has occurred.
Whether re-distribution can occur depends on whether condition (13) is fulfilled. For homogeneous stationary states it can be shown 51 that for zero magnetization miscible solutions exist only provided q < 2c 2 n , which resembles (13). Only then can the spin-changing processes ψ 0 &ψ 0 ⇋ ψ 1 &ψ −1 for condensed atoms start, so (13) works as a kind of ignition condition. After the appearance of condensed atoms in ±1 components, the real www.nature.com/scientificreports/ re-equilibration can ensue. Within this stage, the thermal atoms are redistributed to ±1 states via spin-preserving and spin-exchanging collisions ψ 0 &ψ ±1 ⇋ ψ 0 &ψ ±1 , and become ready for removal. The present study is complementary to the work on sodium in 18 (supplement) which used a Bogoliubov approximation to investigate the high condensate fraction regime. The present work considers low initial condensate fraction, and we find that the Bogoliubov regime of condensate fraction 0.9 can be reached after just several cooling cycles. The two works mesh then at intermediate temperatures, covering the whole range of cooling from near T c to very low T. Cooling cycle efficiency becomes very high in the phonon regime when T ≪ µ/k B 18 , due to a softer phonon spectrum for the m F = ±1 states ( For a density n. Then, below a given energy E max , which corresponds to wavenumbers k 0,± max = E max / √ c 0,2 n , respectively, the number of states in volume V is N 0,± ∼ (4π/3)(V /(2π) 3 )(k 0,± max ) 3 . As a result, the number of accessible states at the beginning of each cooling cycle goes from N 0 to N 0 + 2N ± . This is an increase by a factor of about 1 + 2(c 0 /c 2 ) 3/2 , rather that merely by 3 as seen here at higher temperatures. It has been argued that this kind of cooling may be ultimately only limited by quantum depletion effects at extremely low entropies and thermal fractions.
Finally, we have investigated whether cooling is seen in 23 Na via the process established in the "Chromium: condensate assisted dipolar scattering to higher spins" section. The condition (10) suggests values of q ≈ 20 ω x to be appropriate as an energy barrier that can only be overcome by thermal atoms. Figure 10 (left panel) shows the results of such a simulation. Negligible cooling is observed over this timescale. In fact, this is not unexpected, because there are no dipolar interactions. Therefore, the coherent rapid process in which a thermal atom in m F = 0 is scattered by the condensate to m F = 1 or −1 cannot take place. The incoherent scattering rate can be estimated via Ŵ = nvσ with cross-section 68 σ = 4π((a 2 − a 0 )/3) 2 , and thermal relative velocity v = 4 √ k B T/πm 69 . Integratng over space as in 69 yields a time constant for the m F = ±1 populations of τ = 2 3/2 / Ŵ ∼ 15 s, well beyond typical experimental and simulation times. We expect, though, that a nonzero initial magnetization can help because it opens up other scattering channels. The right panel of Fig. 10 shows that this is indeed the case.
Non-standard initial conditions. To consider values of initial magnetization other than zero, or initial nonzero populations in the m F = ±1 components would significantly increase the complexity of the studies, since in such a case the problem becomes dependent actually on three parameters: the quadratic Zeeman energy and the numbers of atoms in m F = +1 and m F = −1 states. On the other hand, as already discussed at the end of the previous section and shown in Fig. 10  Leaving detailed analysis of cooling processes in full three-parameter space for a future work, here we just present some intriguing results for the "non-standard initial conditions" case. We have found that if a nonzero magnetization is introduced, or if there is minority occupation of m F = ±1 states initially, then the range of q that allows cooling increases. Two examples are shown in Fig. 11 for q = 0.2 ω x , a value that is too high to lead to cooling when all atoms start in m F = 0 . Here the initial populations in m F = ±1 are obtained by coherently transferring part of the initial cloud from the m F = 0 hyperfine state at t = 0 , e.g. by applying an appropriate short Bragg pulse to the cloud 67 . Why cooling appears at higher q here is not obvious and remains an open question. It may, however provide an additional useful pathway. For example, the case with nonzero seed in m F = 1 can still be cycled by emptying only the m F = −1 spin state, or cooling at a given q could be started by producing small seeds at m F = ±1 with a Bragg pulse.

Conclusions
The spin distillation cooling mechanism as originally proposed and partially experimentally tested in 18 , was studied using advanced numerical models that allow for a realistic treatment of temperature, Bose enhancement, and nonlinear dynamics over experimental timescales. These calculations confirm that successive cycles of cooling can be applied to a spinor gas, leading to a very high final condensate fraction (97% in our case), while simultaneously www.nature.com/scientificreports/ losing only a minority of atoms. Spin distillation was observed both for systems with long-range ( 52 Cr, 7 spin states) and short-range ( 23 Na, 3 spin states) forces. We identify two different mechanisms for a cooling cycle: {1} With a high Zeeman energy barrier, as seen in the chromium example, the thermal atoms dipolarly scattering with condensed ones are able to surmount the energy barrier. Such a mechanism moves thermal energy out of the initial cloud to relatively few atoms in higher spin states, from which they and the entropy they carry with them can be cyclically removed. This process works for magnetic fields above a threshold (11), and is eventually limited to temperatures of (12). Notably, these can be very low temperatures far below the k B T ∼ µ that is possible with standard evaporative cooling. Also, since the scattered atoms are in a different mode space from the source cloud, removal does not eat into the condensate like in standard evaporation.
{2} With a very low Zeeman energy barrier, as seen in the sodium example, an initially spin polarised thermal cloud relaxes into m empty hyperfine states, leaving only 1/(m + 1) of it in the original spin state. Here the cooling can be ignited only below a maximum threshold for the magnetic field (13), and is not highly sensitive to the value of the temperature. The process can work at very small magnetic fields, far below the bounds set by (11) and appears quite universal in its simplicity.
There are similarities in the mechanisms we have identified. In both of them the existence of a condensate is crucial. For chromium, thanks to Bose enhancement, thermal atoms are scattering off condensed ones to higher spin states, state by state. For sodium, coherent transfer of already condensed atoms is fast and triggers the process of redistribution of thermal atoms. Also, in both mechanisms, the role of contact interactions in thermalization of newly populated states is similar. In the chromium case, they allow thermal atoms to reach equilibrium within, but not between, spin components. For sodium they arrest the Rabi oscillations of condensed and noncondensed atoms between components, first making oscillations of thermal atoms die out and finally decohering condensed and noncondensed atoms.
Our cooling protocol makes some idealisations. For instance, contrarily to what happens in real experiments, we assume perfect and instantaneous removal of unwanted atoms. This idealization can be partially relaxed, simply by trying simulations where a small remnant fraction of atoms is left in m s = −2, −1, . . . , 3 for chromium and m F = ±1 components for sodium after each cooling cycle. We have checked that the cooling efficiency remains unchanged to within the noise level in the cooling plots even for a fraction of retained atoms equal to about 10%.
The efficiency of evaporative cooling is commonly characterised by the parameter γ = −d(ln D fin /D ini )/d(ln N fin /N ini ) 70 , where D ini ( D fin ) is the initial (final) phase-space density of a gas and N ini and N fin are the corresponding numbers of atoms. The cooling efficiency in the first cooling cycle can be calculated assuming D fin and N fin are taken at the end of the cycle and concern chromium and sodium atoms in the m s = −3 and m F = 0 components, respectively. In our case, to estimate the phase-space density we take the product of the spatial extents of the atomic cloud in the x, y, and z directions with the respective extents in momentum space available via Fourier transform of the classical field. The classical field is averaged over a short period of time to diminish irrelevant fluctuations before determining the size of the cloud (as the full size at half maximum). As a result, we obtain a value of parameter γ ≈ 1 both for chromium and sodium atoms, which is In summary, spin distillation cooling emerges as a versatile technique to reach ultra low temperatures in spinor gases, below what is achievable via standard evaporative cooling. It appears to be applicable in a wide range of regimes and for many atomic species, acting through at least two distinct physical mechanisms. The simulations presented above verify that it can be effective under realistic conditions.

Matrix elements for the dipolar Hamiltonian
The dipolar matrix element [H d ] 11 is given in CGS units by 43 where γ Cr = g L µ B / is the gyromagnetic ratio for 52 Cr. The orthonormal set of spatial coordinates r = [x, y, z] are arranged such that z lies along the direction of the applied magnetic field B . The off-diagonal dipolar matrix element [H d ] 10 is Received: 2 December 2020; Accepted: 23 February 2021 (14) [H d (r)] 11 = ( γ Cr ) 2 dr ′ 1