Magnetic phases of spin-1 spin–orbit-coupled Bose gases

Phases of matter are characterized by order parameters describing the type and degree of order in a system. Here we experimentally explore the magnetic phases present in a near-zero temperature spin-1 spin–orbit-coupled atomic Bose gas and the quantum phase transitions between these phases. We observe ferromagnetic and unpolarized phases, which are stabilized by spin–orbit coupling's explicit locking between spin and motion. These phases are separated by a critical curve containing both first- and second-order transitions joined at a tricritical point. The first-order transition, with observed width as small as h × 4 Hz, gives rise to long-lived metastable states. These measurements are all in agreement with theory.

M ost magnetic systems are composed of localized particles such as electrons 1 , atomic nuclei 2 and ultracold atoms in optical lattices [3][4][5][6] , each with a magnetic moment m. By contrast, itinerant magnetism 7,8 describes systems where the magnetic particles, here ultracold neutral atoms, can themselves move freely, and for which magnetism is generally weak. Our spin-orbit-coupled Bose-Einstein condensates [9][10][11][12] (BECs) constitute a magnetically ordered itinerant system in which-unlike more conventional spinor BECs 13 -the atoms' kinetic energy explicitly drives a phase transition between two different ordered phases 10 . While the coupling between spin and momentum afforded by spin-orbit coupling (SOC) is insufficient to stabilize ferromagnetism in itinerant fermionic systems 14 , in bosonic systems it leads to magnetic phases that are not present in spinor BECs without SOC 11,12 .
Phase transitions can generally be described in terms of a free energy G(M z )-including the total internal energy along with thermodynamic contributions that are negligible for our nearly zero temperature T ¼ 0 system-that is minimized for an equilibrium system. Here the magnetization M z ¼Ŝ z =' is an order parameter, associated with the spinŜ, which changes abruptly as our system undergoes a phase transition. Figure 1c shows typical T ¼ 0 free energies: a first-order phase transition (top panel) occurs when the number of local minima in G(M z ) stays fixed, but the identity of the global minima changes; and a second-order phase transition (bottom panel) occurs when degenerate global minima merge or separate. These defining features are true both for T40 thermal and T ¼ 0 quantum phase transitions.
For spin-1/2 systems (that is, total angular momentum, f ¼ 1/2) like electrons, ferromagnetic order can be represented in terms of a magnetization vector M¼Ŝ =' . This is rooted in the fact that the three components of the spin operatorŜ transform vectorially under rotation. More specifically, any Hamiltonian describing a two level system may be expressed as H¼' O 0 þ X 1 ÁŜ, the sum of a scalar (rank-0 tensor) and a vector (rank-1 tensor) contribution. The former, described by O 0 gives an overall energy shift, and the latter takes the form of the linear Zeeman effect from an effective magnetic field proportional to X 1 . Going beyond this, fully representing a spin-1 (total angular momentum f ¼ 1 with three m F sublevels: | À 1i, |0i, and | þ 1i) Hamiltonian with angular momentumF requires an additional five-component rank-2 tensor operator-the quadrupole tensorand therefore there exist 'magnetization' order parameters that are not simply associated with any spatial direction 13,15,16 .
Pioneering studies in GaAs quantum wells 17,18 showed that material systems with equal contributions of Rashba and Dresselhaus SOC described by the term ak xFz , subject to a transverse magnetic field with Zeeman term O 1Fx , can equivalently be described as a spatially periodic effective magnetic field. Our experiments with spin-1 atomic systems use 'Raman' lasers with wavelength l to induce SOC of this form 9,19-24 with strength a ¼ 2:k R /m, where the single-photon recoil energy and momentum are E R ¼' 2 k 2 R =2m and :k R ¼ 2p:/l. This atomic system can therefore be described by the magnetic Hamiltonian describing atoms with mass m and momentum :k interacting with an effective Zeeman magnetic field x)e y helically precessing in the e x À e y ; and an additional Zeeman-like tensor coupling with strength O 2 . Here,F 2 ð Þ zz =' ¼F 2 z =' 2 À 2=3 is an element of the quadrupole tensor operator. The competing contributions between kinetic and magnetic ordering energies (interactions select between different nearly degenerate ground states, but only weakly perturb the location of the phase transitions, see Methods) make ours an archetype system for studying exotic magnetic order and understanding the associated quantum phase transitions, of which both first and second order are present in our experiments (Fig. 1b).
We can easily understand the first-order transition in the limit of infinitesimal O 1 where the tensor field favours either: a polar BEC for O 2 40 (m F ¼ 0: unmagnetized, M z ¼ 0), or a ferromagnetic BEC for O 2 o0 (m F ¼ þ 1 or À 1: magnetized, |M z | ¼ 1). As with spinor BECs 25 , these phases are separated by a first-order phase transition at O 2 ¼ 0. The ferromagnetic phase spontaneously breaks the Z 2 symmetry associated with the Hamiltonian's invariance under the exchange | À 1i2| þ 1i. The second-order transition can be intuitively understood by considering the large O 1 limit where the system forms a spin-helix BEC (with local magnetization antiparallel to X 1 :  unmagnetized, M z ¼ 0). This order increases the system's kinetic energy, leading to the second-order phase transition into the ferromagnetic phase shown in Fig. 1b as O 1 decreases. Analogues to this second-order phase transition are present in other systems with effective spin degrees of freedom such as double-leg ladders 26 or engineered optical lattices 27,28 . These two extreme limits continuously connect at the point , the green star in Fig. 1b, where the small-O 1 first-order phase transition gives way to the large-O 1 second-order transition, and together these regions constitute a curve of critical points Here we realized spin-1 spinorbit-coupled BECs and varied the magnetic coupling fields using externally applied fields. By directly measuring the system's magnetization, we studied the associated quantum phase transitions present in the phase diagram, all in quantitative agreement with theory.

Results
Experimental set-up. As shown in Fig. 1a, we realized this magnetic system by illuminating 87 Rb BECs in the f ¼ 1 ground state manifold with a pair of counter propagating and orthogonally polarized Raman lasers that coherently coupled the manifold's m F states. Physically, the spatial interference of the orthogonally polarized laser beams give rise to the helical effective magnetic field (see Methods) with period l/2. As we first showed 9 using effective f ¼ 1/2 systems, this introduces both a spin-orbit and a Zeeman term into the BEC's Hamiltonian, equivalent to equation (1). Here the quadratic Zeeman shift from a large bias magnetic field B 0 e z split the low-field degeneracy of the | À 1i2|0i and |0i2| þ 1i transitions, and we independently Raman coupled these state pairs with equal strength O 1 . We dynamically tuned the quadrupole tensor field strength O 2 by simultaneously adjusting the Raman frequency differences; as shown in Fig. 1 we selected frequencies differences where the detuning from the | þ 1i to |0i and | À 1i to |0i were both equal to O 2 (see Methods). Without this technique, only the upper half plane of the phase diagram (Fig. 1b) would be accessible: containing only an unmagnetized phase, therefore lacking any phase transitions.
In each experiment, we first prepared BECs at a desired point in the phase diagram, possibly having crossed the phase transition during preparation. A combination of trap dynamics 29,30 , collisions and evaporation 31 kept the system in or near (local) thermal equilibrium. We then made magnetization measurements directly from the Bose-condensed atoms measured in the spin resolved momentum distribution obtained using the time-of-flight (TOF) techniques described in ref. 30.
Critical line of phase transitions. Our experiment first focused on thermodynamic phase transitions. We made vertical (horizontal) scans through the phase diagram by initializing the system in the unmagnetized phase at a desired value of O 1 (O 2 ) with O 2 \0 (O 1 t10E R ), and then ramping O 2 (O 1 ) through the transition region. (As discussed in the methods our nominally horizontal scans of O 1 followed slightly curved trajectories through the phase diagram, such as the red dashed curve in Fig. 2c). Following such ramps, domains with both ±M z can rapidly form, and we therefore focus on the tensor magnetization M zz ¼F , which is sensitive to this local magnetic order.
Using horizontal scans, we crossed through the second-order phase transition O 2 oO Ã 2 À Á where the free energy evolves continuously from having one minimum (with M zz ¼ 0, for large O 1 ) to having two degenerate minima (with M zz 40, for smaller O 1 ). As shown in Fig. 2a, M zz continuously increases with decreasing O 1 , reaching its saturation value as O 1 -0. Repeating these processes for O Ã 2 oO 2 o0, we found a sharp firstorder transition. In each case, data are plotted along with theory with no adjustable parameters. Using data of this type for a range of O 2 and fitting to numeric solutions of equation (1), we obtained the critical points plotted in Fig. 2c. Because horizontal cuts through the phase diagram are nearly tangent to the transition curve for small O 2 , this produced large uncertainties in O C 1 for the first-order phase transition.
We studied the first-order phase transition with greater precision by ramping O 2 through the transition at fixed O 1 (Fig. 2b) and found near perfect agreement with theory. For all the experimentally measured critical points, see Fig. 2c top, separating the unmagnetized and ferromagnetic phase, we also measured the corresponding transition width defined as the required interval for the curve to fall from 50 to 20% of its full range. This width D decreases sharply at O Ã 1 , marking the crossover between secondand first-order phase transitions (see Fig. 2c bottom). In these data, the width of the first-order transition becomes astonishingly narrow: as small as 0.0011 This narrowness results from the energetic penalty associated with condensation into multiple modes for repulsively interacting bosons. In addition we find that ramps through the first-order transition are hysteretic, and very slow ramps (see Methods) for the system to equilibrate.
During the long equilibration times required to study this transition, spin-domains formed in our system, shown in Fig. 2e. For Bose-condensed atoms, our TOF procedure expanded the initial spin-distribution allowing us to reconstruct any in situ spin structure. Figure 2e shows that domains form as the systems enters into the magnetized phase; these domains mark the role of interactions in spontaneously breaking the Hamiltonian's Z 2 symmetry (see Methods for a comparison with sodium where the sign of the spin-dependent interactions is reversed, and the Z 2 symmetry remains unbroken for a wide range of parameters). Figure 2b shows that additional tripartite spin structure is present very near the first-order phase transition, which was not anticipated in our initial single-particle description. This tripartite mixture, predicted in ref. 32 is an in-plane ferromagnetic phase with no analogue in spinor BECs or effective spin-1/2 SOC BECs.
Metastable states. We observed that scans crossing the second-order transition typically required 50 ms to equilibrate, while for scans crossing the first-order transitions we allowed as long as 1.5 s for equilibration. Systems taken through a first-order phase transition can remain in long-lived metastable states. Here a metastable state with M z ¼ 0 persists in the ferromagnetic phase, and a pair of metastable states with M z a0 persists in the unmagnetized phase. We began our study of this metastability by quenching through the first-order transition at O 1 ¼ 0.74(8)ER with differing rates from À 0.5 to À 0.2 E R s À 1 , as shown in Fig. 3. We observed the transition width continuously decreases with decreasing absolute value ramp rate (inset to Fig. 3), consistent with slow relaxation from a metastable initial state.
We explored the full regime of metastability by initializing BECs in each of the m F states, at fixed O 2 , then rapidly ramping O 1 (t) from zero to its final value fast enough that the system did not adiabatically follow into the true ground state, yet slow enough that the quasi-equilibrium metastable state was left near its local equilibrium. We found that the rate t200E R s À 1 was a good compromise between these two requirements. For points near the first-order phase transition three metastable states exist (Fig. 4); near the second-order transition this count decreases, giving two local minima which merge to a single minimum beyond the second-order transition.
We experimentally identified the number of metastable states by using M z and its higher moments, having started in each of the three m F initial states. A small variance in M z , o0.25, indicates the final states are clustered together-associated with a single global minimum in the free energy G(M z )-and it increases when metastable or degenerate ground states are present. We distinguished systems with two degenerate magnetization states (M z E±1) from those with three states by the same method, since when M z E ± 1, the variance of |M z | is smaller than 0.25, and it distinguishably increases beyond 0.25 as a third metastable state appears with M z ¼ 0. In this way we fully mapped the system's metastable states in agreement with theory, as shown in Fig. 4 Discussion We accurately measured the two-parameter phase diagram of a spin-1 BEC, containing a ferromagnetic phase and an unmagnetized phase, continuously connecting a polar spinor BEC to a spinhelix BEC. The ferromagnetic phase in this itinerant system is stabilized by SOC, and vanishes as the SOC strength :k R goes to zero. Our observation of controlled quench dynamics through a first-order phase transition opens the door for realizing Kibble-Zurek physics 33,34 in this system, where the relevant parameters can be controlled at the individual Hz level. The quadrupole tensor field 0 /F 2 ð Þ zz studied here is the q ¼ 0 component of the rank-2 spherical tensor operatorF 2 ð Þ q , with qA{ ± 2, ± 1, 0}. The physics of this system would be further enriched by the addition of the remaining four tensor fields. The q ¼ 0 term we included is the most simple of the tensor fields to deploy, as it only required control over frequencies. The q ¼ ± 1 components are relatively simple to incorporate by radio frequency (RF)-coupling the |m F ¼ À 1i to |m F ¼ 0i and |m F ¼ þ 1i to |m F ¼ 0i transitions with different phases. The q ¼ ±2 components require direct coupling between |m F ¼ þ 1i and |m F ¼ À 1i, which is straightforward using two-photon microwave transitions, but is challenging to include with significant strength.

Methods
System preparation. We created NE4 Â 10 5 atoms 87 Rb BECs in the ground electronic state |f ¼ 1i manifold 35 , confined in the locally harmonic trapping potential formed at the intersection of two 1,064-nm laser beams propagating along e x and e y giving trap frequencies of (o x , o y , o z )/2p ¼ (33(2), 33(2), 145(5)) Hz. The quadratic contribution to the B 0 ¼ 35.468(1) G magnetic field's Eh Â 25 MHz Zeeman shift lifts the degeneracy between the |f ¼ 1, ð Þ kHz. We denote the energy differences between these states as :do À 1,0 , :do 0, þ 1 and :do À 1, þ 1 . independently Raman coupled the |m F ¼ À 1, 0i and |m F ¼ 0, þ 1i state-pairs, respectively. Furthermore we selected 2o À À o þ À 1 þ o þ þ 1 À Á ¼do À 1; þ 1 such that, after making the rotating wave approximation (RWA) the |m F ¼ ±1i states were energetically shifted by the same O  The physical basis for SOC using Raman lasers. As explained in ref. 36, the spin-dependent (vector) part of the light-matter interaction can be written in terms of an effective Zeeman field with overall strength given by u v as defined in ref. 36, where E is the total optical electric field from the Raman lasers. This field then enters into the Hamiltonian as an effective magnetic field The electric field for the laser geometry depicted in Fig. 1a is , giving rise to the time-dependent effective Zeeman coupling term where we defined do AE 1 ¼o À À o þ AE 1 . This gives the Hamiltonian term In our experiment the o Z /2p E25 MHz linear Zeeman shift is large compared with all other energy scales, so we make the RWA to arrive at the low-frequency Hamiltonian We selected our frequencies to be in four-photon resonance with the | À 1i to | þ 1i transition, giving (do , in which case the Hamiltonian is time-periodic. Equation (1) in the manuscript is obtained by making independent RWAs on the | À 1i-|0i and |0i-| þ 1i transitions separately, giving the helically precessing coupling described in the manuscript.
Floquet and polynomial shift. The frequency differences do À 1;0 % o À À o þ À 1 and do 0; þ 1 % o À À o þ þ 1 nominally Raman dress both | À 1, 0i and |0, þ 1i state pairs independently. In practice, the cross coupling may be substantial and the adjusted eigenenergies may be computed exactly from Floquet theory. For our Measurement details. In the horizontal scans in Fig. 1 of the manuscript, we ramped O 1 at a rate of E À 40E R /s, allowing the system to adiabatically track the ground state, and allowed 50 ms for equilibration before the measurement process. In contrast, for our vertical scans we found the system required between 0.2 ms and 2 s to equilibrate. We studied the metastable states by performing three separate experiments for each raw data point: one each for a system initially prepared in m F ¼ 0, ±1 state at the desired value of O 2 . We then adiabatically ramped O 1 from zero to its final value. For each resulting (O 1 , O 2 ) pair, we obtained the magnetization M z for each initial state. In all of the described procedures O 1 was turned off immediately, at the start of our 28 ms time-of-flight imaging process, which included a Stern-Gerlach gradient to separate the spin components.  Figure 3 | Quenching dynamics. The system was prepared in the unmagnetized phase with O 1 ¼ 0.74(8)E R and O 2 was ramped through the phase transition at ramp-rates dO 2 /dt ¼ À 0.2, À 0.3, À 0.4, and À 0.5E R s À 1 (blue, black, red, and green symbols, respectively). The curves are guides to the eye. The inset shows the decreasing width, defined as the required interval for the curve to fall from 50 to 20% of its full range, of the first-order transition as the ramp-rate decreases. Vertical error bars correspond to one standard deviation of up to 10 measurements, and horizontal error bars correspond to a systematic error of 0.005E R from fluctuations in the bias field. Error bars on the inset panel correspond to the 95% confidence interval on the fitting function of the quenching data. Top, Measured magnetization plotted along with theory. The system was prepared at the desired O 2 ¼ À 2E R ; O 1 (t) was then increased to its displayed final value; during this ramp O 2 also changed, and the system followed the curved trajectory in the bottom panel. Each displayed data point is an average of up to 10 measurements, and the coloured region reflects the uncertainty in theory resulting from our E5% systematic uncertainty in O 1 . Circles/crosses/stars represent data starting in m F ¼ þ 1, 0, and À 1 respectively. Bottom, state diagram: theory and experiment. Blue: two states; black: three states; white: one state. Coloured areas denote calculated regions where the colour-coded number of stable/metastable states are expected. Symbols are the outcome of experiment. Each displayed data point is an average of up to 20 measurements.
Operators. The total angular momentum f ¼ 1 spin operators in equation (1) of the main manuscript take the explicit form in the basis of the magnetic sublevels | À 1i, |0i, and | þ 1i; together these comprise the vector operatorF¼F x e x þF y e y þF z e z . Likewise the quadrupole tensor operator is expressed asF In terms of these operators it is clear than any Hamiltonian involving only1,F x ,F y andF zz is invariant under the transformation that swaps | þ 1i and | À 1i: a discrete Z 2 symmetry.
Wavefunctions. The wavefunctions in the polarized and unpolarized regimes are qualitatively different. In the unpolarized regime the wavefunction takes the general form c The value of A depends in detail on O 1 and O 2 , but two limits are clear. First, when O 1 ¼ 0 and O 2 40 the system forms a spinor BEC in the polar phase, with In the polarized regime, the wavefunction has the general form (the definition of magnetization); and second, M z ¼ À k 0 /2k R (ensuring zero center-of-mass motion, see ref. 30). In our experiment m F ¼ ± 1 are coupled at second order in O 1 , and for M z j jt1, these states are E16E R detuned. Thus, for M z t þ 1, we have A À 1 E0 with corrections at order O 2 1 , giving the wavefunction in terms of the magnetization.
Free energy and phase diagram. We obtained the free energy G(M z ) as a function of the magnetization M z by first numerically solving the system's Hamiltonian given by equation (1), obtaining the eigenenergies E s (k) and state c s k ð Þ, each identified by a momentum :k and a 'band' index sA{ À 1, 0, þ 1}. We then computed M z for each of these states (dependent on k x , but independent of k y and k z ), thereby obtaining the internal energy E(M z ) in the lowest band (s ¼ À 1).
As our BEC is very near the ground state the free energy G(M z ) ¼ E(M z ) À TS, where T is the temperature and S is the entropy is well approximated by G(M z )EE(M z ), and it is this free energy, which is plotted in Fig. 1c. We then obtained the phase diagram in Fig. 1b by numerically computing the free energy for each pair O 1 , O 2 and identifying its equilibrium magnetization.
For non-interacting systems, we found that the curves defining the phase transitions and also those bounding the region containing metastable states could all be computed in closed form. First, the critical point at which the first-and second-order phase transitions meet is at The curve defining the first-order phase transition (for 0 ! O 2 ! O Ã 2 ) is given by and the curve defining the second-order phase transition (for O Ã 2 ! O 2 ) is given by The upper boundary of the metastable regime (in the unmagnetized phase) is given by and the lower boundary of the metastable regime (in the ferromagnetic phase) is given by This is the same equation defining the second-order phase transition with an added ±, the full curve defining the boundary of the metastable regime crosses over between the þ and À solutions at Magnetic fields. Because the free energy G(M z ) is sensitive to unwanted detuning d from the four-photon resonance near the phase transitions, which contributes an added symmetry breaking field dF z to the Hamiltonian, controlling the bias magnetic field and nulling its gradients is critical. A pair of flux-gate sensors measuring the ambient magnetic field along e z , allowed us to compensate for longterm field drifts. We compensated any field gradients using four pairs of anti-Helmholtz coils in a clover leaf configuration 37 , and a conventional anti-Helmholtz pair, all aligned along e z .  Interactions. We studied the impact of interactions at the level of mean field theory using a variational approach, assuming an infinite homogenous system. For each point in the phase diagram labelled by (O 1 , O 2 ), we first located the local minima in the single-particle free energy described by the HamiltonianĤ of equation (1) in the manuscript. The free energy had from one to three minima, with energies E j and eigenstates c j E . We then considered an infinite system and minimized the mean field energy density for an arbitrary linear combination of these single particle states with amplitudes a j , where n s (x) is the local density in a given spin state s; n T (x) is the total local density; and c 0 and c 2 are the spin-independent and spin-dependent interactions, respectively. For 87 Rb87 these have the ratio c 2 /c 0 E À 0.005 and for 23 Na they are c 2 /c 0 E þ 0.05. In our minimization, we modelled our systems with a typical mean-field energy of (c 0 þ c 2 )n T E1 kHz per particle. Figure 5 shows the result of this calculation both for rubidium and sodium. In both cases the overall phase diagram (Fig. 5a,d) is shaped by the single-particle Hamiltonian; at this coarse level the rubidium phase diagram is hardly different from that predicted from single particle physics, but in the case of sodium a large swath of the expected ferromagnetic phase remains symmetry unbroken. This phase continuously connects to an equal superposition of | À 1i and | þ 1i as O 1 -0.
The situation becomes more complex as we focus on M z near the curve defining the first-order phase transition (Fig. 5b,e), for the case of rubidium a new mixed phase appears at low O 1 analogous to the striped phase in spin À 1/2 systems, but nothing new is apparent for sodium.
Last, we consider the same region, but looking at the fraction of the variational wavefunction in the m F ¼ 0 spin state. For rubidium, this allows us to identify a new state which is a three-way mixture of all three components considered in the variational calculation (with no analogue in the O 1 ¼ 0 spinor limit), and we can see the abrupt transition in sodium from a state connecting to the polar phase (O 2 40) and to the uniaxial nematic phase (O 2 o0). Each of these phases are as described in ref. 32.