Thermally-robust spin correlations between two atoms

Bottom-up assembly of few- and many-body systems from individual atoms could deliver robust entanglement, which is a key resource for quantum technologies\cite{Isenhower2010, Regal2015, Endres2016}. Furthermore, it allows quantum processes to be observed at the individual event level, revealing information concealed by ensemble-averaged measurements\cite{Toschek1986, Andersen2010} and providing insights into molecular processes\cite{Andersen2013, Liu2018}, thermalization\cite{Greiner2016}, and quantum thermodynamics\cite{Rosnagel:2016}. Here, we study the spin dynamics of a two-body system consisting of individually-assembled pairs of \boldsymbol{$^{85}$}Rb atoms, whose collisional properties prevent investigation in the many-body regime. The thermal spin-2 atoms show perfect pair correlation between magnetic sublevels on timescales exceeding one second, with measured relative number fluctuations \boldsymbol{$11.9\pm0.3$} dB below quantum shot noise (QSN), limited only by detection efficiency. Both microscopic simulations and experiments display relaxation dynamics, contrary to the coherent spin waves witnessed in finite-temperature many-body experiments\cite{Lett2013} and zero-temperature two-body experiments\cite{Bloch2005}. The relative relaxation rates are consistent with theoretical predictions of the \boldsymbol{$^{85}$}Rb spin-dependent interaction strengths\cite{Greene2001}. Our experiment is a versatile platform for studying two-body quantum dynamics and may provide a route to thermally-robust entanglement generation.


Load two atoms
wherer = (r x ,r y ,r z ) andp are relative position and momentum operators, respectively, µ the reduced mass, ω j the atomic oscillation frequency in the j th dimension, andĤ Z,i the Zeeman shift for the i th atom. Our experiments use thermal atoms with k B T much larger than ω j , Zeeman energies, and atomic interaction energies.
By measuring magnetic sublevels of the atomic pair for different collision times, we confirm that the spin dynamics is governed by the simple model of spin-changing collisions depicted in Fig. 1b, which yields strong correlations between the m-states in a given pair. This requires the three measurement series summarised in Fig. 2 Probability that zero, one, or two atoms remain in the optical tweezer after a given collision time. a, When atoms in |0 are expelled (immediately after a given collision time), the probability that both atoms were in |0 (and therefore ejected) decreases, while probability that both atoms remain correspondingly increases. b, Expelling atoms solely from |−1 gives only single-atom loss events, which is opposite to the result in c. c, Expelling atoms from both |−1 and |1 gives only pair loss. In all cases and throughout the collision time, the bias magnetic field was 8.5 Gauss.

. A particular |m is
detected by ejecting atoms in this state. In Fig. 2a we expel atoms in |0 after a given collision time. The probability that both atoms are in |0 (i.e. no remaining atoms) decays with increasing collision time, while the probability that both atoms remain grows correspondingly. The probability of observing one remaining atom is always negligible, implying that collisions cause both atoms to leave |0 simultaneously. In Fig. 2b we start with both atoms in |0 but eject atoms in | − 1 . The probability that one atom is in |−1 grows with collision time, but both are never |−1 , since the probability that both atoms are ejected is effectively zero. In Fig. 2c we eject atoms in both |−1 and |1 . This ejects both atoms, or none. Combining this with Fig. 2b, we conclude that when one atom is in |−1 , the other is in |1 . The populations of |−1 and |1 are therefore almost perfectly correlated. Similar data for |±2 shows these populations are also strongly correlated. The lasting pair correlation on timescales exceeding one second is facilitated by having individual atomic pairs. In contrast, in many-body experiments with spin-2 atoms, subsequent spin-changing collisions would likely deteriorate such perfect pair correlations.
We quantify the pair correlation with the relative number squeezing, ζ 2 (see Methods). Without correcting for finite detection efficiency, it is 11.9 ± 0.3 dB below QSN for the |±1 populations. Since our atomic-pair ensemble is thermal, this large pair correlation is thermally robust. ζ 2 is limited solely by our detection efficiency (see Methods); improved detection efficiency could reduce ζ 2 by a further order of magnitude. For many-body systems, the highest reported relative number squeezing via spin-changing collisions is 11.4 dB below QSN (12.4 dB after correcting for detection inefficiency) [30].
The bias magnetic field affects the spin dynamics through iĤ Z,i . Since our model conserves total magnetization, the first-order Zeeman contributions cancel for the accessible two-body states, so iĤ Z,i only contributes via secondorder terms. We investigate how iĤ Z,i affects the spin dynamics by measuring the |0, 0 population after a 40 ms collision time for different bias fields (Fig. 3). At low biases, the dynamics are highly magnetic-field dependent, whereas for higher biases the dynamics are effectively magnetic-field independent. Here typical thermal energies are much larger than second-order Zeeman energies for all biases investigated. The atom pairs therefore have sufficient thermal energy to overcome the Zeeman shift when undergoing spin-changing collisions, so, in contrast to ultracold samples, the dynamics should not necessarily quench at high biases.
To understand the spin evolution, we simulated the dynamics governed by Eq.
(1) with a simplified interaction H s = V (r) × m1,m2,m 1 ,m 2 g m 1 ,m 2 m1,m2 |m 1 , m 2 m 1 , m 2 |, where g m 1 ,m 2 m1,m2 are determined from predicted spin-dependent swave scattering lengths [12] and V (r) is a Gaussian with width chosen to reproduce the total free-space s-wave collision cross section (see Methods). A Gaussian pseudopotential moderates problems that afflict zero-length interaction potentials in tight traps [28], while still avoiding the complexity of a more completeĤ s .
The simulation was conducted by averaging over a thermal ensemble of initial states evolved using Eq. (1). The initial states were relative-motion eigenstates ofp 2 /(2µ) + j 1 2 µω 2 jr 2 j with two-particle spin state |0, 0 . Due to the prohibitively-large Hilbert space required at the experimental temperature, simulations were restricted to a lower temperature of 8.8 µK.
The simulation qualitatively captures the spin dynamics (Fig. 3). We observe a crossover from fast dynamics at low magnetic-field strengths to slow dynamics at high fields.Ĥ couples the three allowed spin modes, |0, 0 ,Ŝ |1, −1 , andŜ |2, −2 (inset, Fig. 3). When the pair is in a particular spin mode, it behaves as an effective single particle within a harmonic trap with the interaction potential placed at the trap centre. At low magnetic fields, iĤ Z,i is negligible, so any relative-motion eigenstate with a particular spin mode (e.g. |0, 0 ) is approximately degenerate to relative-motion eigenstates in other spin modes (e.g.Ŝ |1, −1 , and/orŜ |2, −2 ); the degeneracy is only lifted by the atom-atom interaction's spin-state dependence. The resulting resonant coupling efficiently transfers population between spin modes at low magnetic fields. In contrast, at high fields this degeneracy is lifted, the majority of initially-occupied states have no near-resonant coupling to other spin modes, leaving only off-resonant coupling, and the dynamics largely cease.
Our experiments and simulations display spin relaxation dynamics for any parameter investigated, contrary to finite-temperature many-body experiments [10] which exhibit high-contrast coherent oscillations between spin modes. This is expected: since coupling between spin modes depends on the initial relative-motion eigenstate, incoherent thermal averaging over initial states leads to relaxation dynamics similar to Fig. 2. Figure 3 shows a quantitative difference between simulation and experiment. In the high bias, magnetic-fieldindependent regime, the simulation gives |0, 0 population at t = 40 ms close to the t = 0 population, while in the experiment it is lower. Figure 4 demonstrates the cause of this difference. The experiment shows slow relaxation to equal populations of the three spin modes, while the simulation dynamics are quenched (no spin-changing collisions). Here equal population is not complete thermalization within states that conserve total magnetization; since atoms with different internal states can be considered distinguishable, the thermalized populations with m = ±1 and m = ±2 would be twice that of |0, 0 .
Generally, a priori calculations of thermal decoherence in colliding atomic ensembles pose a challenge for theory, often necessitating phenomenological rate-equation approaches to account for dissipation [31]. In our system, several effects that are not included in the simulations might explain the dynamics in Fig. 4. Magnetic noise might affect the dynamics or slight polarization pollution of the optical tweezer light could give a slightly m-dependent trap, the latter invalidating our separation of the pair's centre-of-mass and relative coordinates. A more realistic atom-atom interactionĤ s may also introduce new collisional timescales not captured by our simulations' simplified interaction. Finally, the five-fold temperature difference between our simulations' practical limit and the experimental temperature could play a role. However, this appears an unlikely explanation, as the simulation does not reveal long-time dynamics for any of the temperatures we investigated. on the cross section for the process, which is proportional to the squared magnitude of the coupling matrix elements. These are determined from theoretically-predicted 85 Rb spin-dependent interaction strengths [12]. Based on this, the ratio of the rates between |0, 0 Ŝ |1, −1 andŜ |1, −1 Ŝ |2, −2 is 2.34, while the rate between |0, 0 Ŝ |2, −2 is negligible. Fitting using a single overall rate as the fitting parameter matches the data very well (Fig. 4), indicating that collisional dynamics at high magnetic fields can be understood in terms of the cross sections. Figure 3 therefore displays a crossover from a quantum regime at low magnetic fields to a classical regime at high fields, where the collision dynamics do not dependent upon coherences between different spin states.
To conclude: a single pair of 85 Rb atoms in an optical tweezer displays spin dynamics that yield strong correlation between magnetic substates of the two atoms. Unlike both finite temperature many-body experiments and zerotemperature two-body experiments, our finite-temperature two-body experiments show relaxation dynamics rather than coherent spin waves. Our experiments could provide new insights into quantum relaxation processes and quantum thermodynamics. The bottom-up assembling of pairs from individual atoms allows us to study 85 Rb, whose effective attractive interactions are unfavourable for ultracold-ensemble collision experiments [29]. The record-high pair correlation measured is only limited by detection inefficiency. Improving upon this technical limitation might allow studies of unexplored effects, such as violations of total magnetization conservation due to spin-orbit coupling [29].
Our experiments indicate that spin-changing collisions could provide a useful finite-temperature entanglement resource that is robust to thermal noise. The solid curves are a fit of the measured data with spin-changing rate equations, while the ratio of the rates between |0, 0 Ŝ |1, −1 andŜ |1, −1 Ŝ |2, −2 is determined from the theoretically-predicted spin-dependent interaction strengths. The inset illustrates that the simplified theoretical model used for our simulations fails to capture the long-time relaxation dynamics in the high magnetic-field regime.

Experimental procedure.
We initially cool and trap a cloud of 85 Rb atoms using magneto-optical trapping. We then load a small number of atoms from the cloud into two optical tweezers separated by ∼ 4 µm, each with a trap width of ∼1.05 µm and depth of 58 MHz. The two optical tweezers are formed by focusing two steerable linearly polarized laser beams (λ = 1064 nm) with a high-numerical-aperture lens (NA = 0.55). We use blue-detuned light-assisted collisions to reduce the occupancy of each trap to a single atom and confirm the presence of the two isolated atoms via fluorescence imaging [5,26,27]. The probability that there are two atoms, one in each tweezer, after the loading procedure is ∼0.64.
After the loading process, the atoms are prepared in the desired f = 2, m = 0 groundstate in two steps. First, we optically pump atoms to the f = 3, m = 0 state by applying linearly-polarized optical pumping light with two frequencies corresponding to the 85 Rb D 1 f = 2 to f = 3 and the f = 3 to f = 3 transitions. During this, the bias magnetic field of 8.5 gauss defines the quantization axis for the atoms in the groundstate. This gives an atomic population of 0.99 occupying the f = 3, m = 0 state. Last, we apply a π-pulse (1.57 µs) of co-propagating Raman beams (∼36 GHz red detuned from the D 2 line) to coherently transfer the atoms from the f = 3, m = 0 state to the f = 2, m = 0 state.
Using a 20 ms frequency sweep of an acousto-optical modulator, we adiabatically bring the two tweezers closer until they are merged (the distance between the centres of the two laser beams is ∼900 nm). We then adiabatically ramp off one of the tweezers in ∼17 ms while the other is simultaneously ramped to the desired trap depth and the bias magnetic field is set to the chosen value. The procedure leaves the atoms in the same optical tweezer where the collisional interactions generate the |m population dynamics.
To observe the results shown in Fig. 2 and Fig. 4, we use the following experimental parameters: a trap depth of h × 58 MHz, oscillation frequencies 2π × 136 kHz and 2π × 22 kHz for the radial and axial dimensions, respectively, an atomic temperature of 107 µK, and a bias magnetic field of 8.5 Gauss. For Fig. 3, we use a trap depth of h × 10 MHz, oscillation frequencies 2π × 56 kHz and 2π × 9 kHz for the radial and axial dimensions, respectively, and an atomic temperature of 44 µK.
The detection of atoms in a particular |m of the f = 2 manifold is done by ejecting the atoms out of the trap. In the presence of the magnetic field, we use a Raman process to transfer only the population in the specific |m to the f = 3 manifold. We then deplete the f = 3 population using the push out technique [33] and then measure the number of remaining atoms in the trap using fluorescence detection [34]. This procedure yields that the lost atoms were in the detected |m while the remaining atoms were in the other states. In our push out technique, the detection efficiencies are 0.944 ± 0.004 and 0.997 ± 0.003 for the f = 2 and f = 3 states, respectively.

Relative number squeezing.
The correlations between the |±1 of the two atoms (shown in Fig. 2) can be characterized by the population imbalance J z = (N +1 − N −1 ) /2, and the relative number squeezing [20] ζ 2 = (∆Jz) 2 N/4 . ∆J z is the standard deviation of J z , N ±1 is number of atoms in |±1 , and N is the total number of atoms. We deduce the number squeezing from the data in Fig. 2c at the collision times of 150, 250, 350, and 500 ms. From the figure, the result of zero atoms remaining in the optical tweezer after the detection indicates that one atom is in |1 and the other is in |−1 , so J z = 0 in this case. For the result of one atom remaining in the tweezer (indicates that one atom is in |±1 while another one is in |0 , |2 , or |−2 ), J z = ±0.5. Then we can calculate ζ 2 = 0.5 2 P1 P1+P0 / N 4 , where P 1 and P 0 are the probabilities of one and zero atoms remaining in the tweezer, respectively. Here, we assume that the probability of having both atoms in |1 or |−1 is zero.
Our measurement of the correlation can be influenced by the detection efficiency since the detection error in both f = 2 and f = 3 states will contribute to the measured value of P 1 . The directly measured variance (∆J z ) 2 is 0.032 ± 0.002, while the detection error gives a variance of 0.034 ± 0.002 under the assumption that the actual (∆J z ) 2 = 0. This shows the measured degree of relative number squeezing can be entirely attributed to the detection efficiency.
Coupling strengths and rate equations.
We deduce the transition rates from the spin-dependent interaction strengths. We assume that the transition rates betweenŜ|m, −m andŜ|m , −m are incoherent and have strengths proportional to | m , −m |ŜĤ sŜ |m, −m | 2 . For low collisional energy, the interaction Hamiltonian of two atoms is approximated by [11] H s = V (r) m1,m2,m 1 ,m 2 g m 1 ,m 2 m1,m2 |m 1 , m 2 m 1 , m 2 | , wherer is the relative position. The coupling coefficient between the initial |m 1 , m 2 and final |m 1 , m 2 of the atom pair is where g F = 4π 2 a F /m with a F the s-wave scattering length for two atoms colliding in a channel with total spin F . As shown in the Supplemental, provided both F = 2 atoms are initially prepared in the m = 0 Zeeman state, there are only six unique coupling coefficients in the above sum: For 85 Rb, the theoretically-predicted s-wave scattering lengths are a 0 = −740 ± 60 a.u., a 2 = −570 ± 50 a.u., a 4 = −390 ± 20 a.u. [12]. By assuming the transition rate γ mm betweenŜ |m, −m andŜ |m , −m is proportional to −0.09 respectively. We therefore set γ 02 to zero in the following rate equations. Ignoring γ 02 , we use the following rate equation to model the experimental results in Fig. 4: where PŜ |m,−m is theŜ |m, −m population. Using the above ratio of rates, we set γ 01 = 2.34 × γ 12 and fit the entire experimental dataset in Fig. 4 using the single fitting parameter γ 12 .
Theoretical model of collisional spin dynamics.
In the low-energy regime where s-wave collisions dominate, it is customary to take V (r) = δ(r) [11]. However, in this case spin-changing dynamics only occur for eigenstates where n x , n y , n z are all even (see Supplemental). In contrast, states where (say) n x is even and n y and n z are odd never evolve. These latter kinds of states represent roughly 70% of the ensemble at 44 µK, implying that this model predicts that the population of |0, 0 never drops below 0.7, at odds with what we experimentally observe.
We wish to use a simplified atom-atom interaction model that allows for numerical calculations involving a high number of modes, while at the same time avoids the problem with the delta-function interaction model [32]. In particular, there is some evidence that the zero-range δ-function pseudopotential fails to replicate the scattering properties of the underlying physical potential in trapped systems when the magnitude of the s-wave scattering length is on the order or greater than the harmonic oscillator lengthscale [28]. Further, there is a greater discrepancy for negative scattering lengths. In our experiment a 0 /d = −0.44, a 2 /d = −0.34, and a 4 /d = −0.23, where d = /(mω) andω = (ω x ω y ω z ) 1/3 . We use a Gaussian pseudopotential V (r) = exp[−r 2 /(2w 2 )]/(2πw 2 ) 3/2 with w = (a 4 0 + a 4 2 + a 4 4 )/(a 2 0 + a 2 2 + a 2 4 ) ∼ 650 a.u., since (1) it is finite range and couples all even-parity eigenstates, (2) it gives the same total scattering cross section as the δ-function pseudopotential (see Supplemental), (3) it smoothly recovers the (regularised) δ-function in the w → 0 limit, and (4) the form of the spin-changing coupling matrix is sufficiently simple that a numeric calculation is tractable.
We numerically solve for the spin-changing dynamics by expanding ψ i (r) on a finite basis of even-parity eigenstates of H rel (r): ψ i (r, t) = εn x,ny ,nz ≤Ecut c i nx,ny,nz (t)ϕ nx (x)ϕ ny (y)ϕ nz (z), where the sum is over all eigenstates with energy nx,ny,nz less than some energy cutoff E cut . It is necessary to choose E cut sufficiently large that εn x,ny ,nz ≤Ecut P(n x , n y , n z ) ≈ 1 and coupling to the highest-energy, sparsely-occupied modes is negligible. For the computational resources at our disposal, these conditions limit our calculations to temperatures no greater than 8.8 µK -roughly one fifth the temperature of the experiment.
In this basis the state is represented by c = [c 0 , c 1 , c 2 ] , where c i is the vector of coefficients c i nx,ny,nz for modes satisfying nx,ny,nz ≤ E cut . Equations (6) Here is a diagonal matrix with energies ε nx,ny,nz along the diagonal and the coupling matrix T is defined via T mx,my,mz nx,ny,nz = I nx,mx I ny,my I nz,mz /(2πw 2 ) 3/2 , where the integrals I ni,mi = dx i ϕ ni (x i ) exp[−x 2 i /(2w 2 )]ϕ i mi (x i ) have an analytic solution in terms of Gauss hypergeometric functions (see Supplemental). Diagonalising H = UDU † gives the solution c(t) = U exp[− i Dt]U † c(0). Thus, for a given initial condition ψ 0 (r, 0) = ϕ mx (x)ϕ my (y)ϕ mz (z) we can compute the population of the jth two-boson spin state N j mx,my,mz (t) = εn x,ny ,nz ≤Ecut |c j nx,ny,nz (t)| 2 . The total population of the jth two-boson spin state assuming a thermal initial state is given by an incoherent sum over N j mx,my,mz (t) weighted by the Boltzmann probability P(m x , m y , m z ): P |j,−j (t) = εm x ,my ,mz ≤Ecut P(m x , m y , m z )N j mx,my,mz (t).
This procedure was used to generate the simulation data plotted in Fig. 3 and Fig. 4.

AUTHOR CONTRIBUTIONS
The experiments were carried out by P.S. and E.S.. The data analysis was performed by P.S.. Theoretical and numerical analysis was performed by S.S.S. and supervised by A.S.B.. All work was supervised by M.F.A.. All authors discussed the results and contributed to the manuscript.