Spin current generation and relaxation in a quenched spin-orbit-coupled Bose-Einstein condensate.

Understanding the effects of spin-orbit coupling (SOC) and many-body interactions on spin transport is important in condensed matter physics and spintronics. This topic has been intensively studied for spin carriers such as electrons but barely explored for charge-neutral bosonic quasiparticles (including their condensates), which hold promises for coherent spin transport over macroscopic distances. Here, we explore the effects of synthetic SOC (induced by optical Raman coupling) and atomic interactions on the spin transport in an atomic Bose-Einstein condensate (BEC), where the spin-dipole mode (SDM, actuated by quenching the Raman coupling) of two interacting spin components constitutes an alternating spin current. We experimentally observe that SOC significantly enhances the SDM damping while reducing the thermalization (the reduction of the condensate fraction). We also observe generation of BEC collective excitations such as shape oscillations. Our theory reveals that the SOC-modified interference, immiscibility, and interaction between the spin components can play crucial roles in spin transport.

Introduction

Spin, an internal quantum degree of freedom of particles, is central to many condensed matter phenomena such as topological insulators and superconductors[1,2] and technological applications such as spintronics[3] and spin-based quantum computation[4]. Recently, neutral bosonic quasiparticles (such as exciton-polaritons and magnons) or their condensates[5-7] have attracted great interest for coherent manipulation of the spin information. For example, spin currents have been generated using exciton-polarions[8] and excitons[9] in semiconductors and magnons[10,11] in a magnetic insulator. In spin-based devices, SOC and many-body interactions are key factors for spin current manipulations. SOC can play a particularly crucial role as it may provide a mechanism (such as spin Hall effect) to control the spin, however, it can also cause spin (current) relaxation, leading to loss of spin information. Studying the effects of SOC and many-body interactions on spin relaxation is thus of great importance but also challenging due to uncontrolled disorders and the lack of experimental flexibility in solid state systems.

Cold atomic gases provide a clean and highly-controllable[12] platform for simulating and exploring many condensed matter phenomena[12-16]. For example, the generation of synthetic electric[17] and magnetic[18] fields allows neutral atoms to behave like charged particles. The synthetic magnetic and spin-dependent magnetic fields have been realized to demonstrate respectively the superfluid Hall[19] and spin Hall effects[20] in BECs. The creation of synthetic SOC in bosonic[21-25] and fermionic[26-29] atoms further paves the way to explore diverse phenomena such as topological states[30] and exotic condensates and superfluids[16,31-35]. Here, we study the effects of one-dimensional (1D) synthetic SOC on the spin relaxation in a disorder-free atomic BEC using a condensate collider, in which the SDM[36] of two BECs of different (pseudo) spin states constitute an alternating (AC) spin current. The SDM is initiated by applying a spin-dependent synthetic electric field to the BEC via quenching the Raman coupling that generates the spin-orbit-coupled (SO-coupled) band structure. Similar quantum gas collider systems (without SOC[37-41]) have been used to study physics that are difficult to access in other systems.

Charge or mass currents are typically unaffected by interactions between particles because the currents are associated with the total momentum that is unaffected by interactions. In contrast, spin currents can be intrinsically damped due to the friction resulting from the interactions between different spin components. In electronic systems, such a friction has been referred to as the spin Coulomb drag[42,43]. In atomic systems, previous studies have shown that a similar spin drag[44,45] also exists. Even in the absence of SOC, the relaxation of spin currents can be nontrivial due to, for example, interactions[36,39,46-49] and quantum statistical effects[45,50]. In one previous experiment[20], bosonic spin currents have been generated in a SO-coupled BEC using the spin Hall effect. However, how the spin currents may relax in the presence of SOC and interactions has not been explored. Here, we observe that SOC can significantly enhance the relaxation of a coherent spin current in a BEC while reducing the thermalization during our experiment. Moreover, our theory, consistent with the observations, discloses that the interference, immiscibility, and interaction between the two colliding spin components can be notably modified by SOC and play an important role in spin transport.

Results

Experimental setup

In our experiments, we create 3D 87Rb BECs in the F = 1 hyperfine state in an optical dipole trap with condensate fraction fc > 0.6 containing condensate atom number Nc ~ 1–2 × 104. As shown in Fig. 1a, counter-propagating Raman lasers with an angular frequency difference ΔωR couple bare spin and momentum states and to create synthetic 1D SOC (so called equal Rashba–Dresselhaus SOC) along [24], where the bare spin states and are Zeeman split by ħωZ ≈ ħΔωR using a bias magnetic field . Here, ħk↓ = ħ(q + kr) (ħk↑ = ħ(q − kr)) is the mechanical momentum in the y direction of the bare spin component , where ħq is the quasimomentum. The photon recoil momentum ħkr = 2πħ/λ and recoil energy are set by the Raman laser at the “magic” wavelength λ ~ 790 nm[51], where ħ is the reduced Planck constant and m is the atomic mass of 87Rb. The state can be neglected in a first-order approximation due to the quadratic Zeeman shift (see Methods). The single-particle SOC Hamiltonian, HSOC, can be written in the basis of bare spin and momentum states as[21]:where Ω is the Raman coupling (tunable by the Raman laser intensity), δR = ħ(ΔωR − ωZ) is the Raman detuning (tunable by B) and is zero in our main measurements (see Methods). A dressed state is an eigenstate of Eq. (1), labeled by q, and is a superposition of bare spin and momentum states. The q-dependent eigenvalues of (1) define the ground and excited energy bands. When Ω is below a critical Ωc, the ground band exhibits double wells, which we associate with the dressed spin up and down states. The double minima at quasimomentum ħq can be identified with the light-induced spin-dependent vector potentials (controllable by Ω), where σ labels or [20] (see Methods). The double minima merge into a single minimum as Ω increases beyond Ωc, as shown in the dashed line trajectories in Fig. 1c.

Fig. 1

Experimental setup and timing diagram used for the spin-dipole mode (SDM) experiments. a Linearly polarized Raman beams with orthogonal polarizations (indicated by the double-headed arrows along and ) counter-propagating along couple m hyperfine sublevels (bare spin states) of 87Rb atoms. The sublevels are Zeeman split by ħωZ ≈ ħΔωR = h × (3.5 MHz) using a bias magnetic field , which controls the Raman detuning δR = ħ(ΔωR − ωZ). b Experimental timing diagram: Raman coupling Ω (with an experimental uncertainty of <10%) is slowly ramped up in 80 ms to an initial value ΩI and held for 100 ms to prepare the BEC around the single minimum of the ground band at ΩI as shown in c. Then, Ω is quickly lowered to a final coupling ΩF in time tE and held for some time thold, during which we study the dynamics of the BEC in the dipole trap. Subsequently, the atoms are released for absorption imaging after a 15 ms time of flight (TOF), at the beginning of which a Stern–Gerlach process is performed for 9 ms to separate atoms of different bare spin states. c The ground band (solid lines) of synthetic SOC is calculated for a few representative Ω at δR = 0. A higher band calculated for Ω = 1.3 Er is shown as dashed lines. The colors indicate the spin compositions, with red for and blue for . The ground band minima in quasimomentum marked by dots are identified with spin-dependent vector potentials (A), which shift in opposite directions as Ω is lowered into the double minima regime during tE. This generates spin-dependent synthetic electric fields E and thus excites the SDM and an AC spin current along the SOC direction in a trapped BEC. The upper (lower) dashed circle represents the region around q = 0 in the double minima band at an exemplary ΩF = 0 (ΩF = 1.3 Er), from which the two (dressed) spin components of the BEC roll down towards the corresponding band minima in response to the application of E

We prepare a BEC around the single minimum of the ground dressed band at ΩI ( = 5.2 Er for this work) and δR = 0 by ramping on Ω slowly in 80 ms and holding it for 100 ms (Fig. 1b, c, see Methods for details). Then, we quickly lower Ω from ΩI to a final value ΩF into the “double minima” regime in time tE. The tE = 1 ms used in this work is slow enough to avoid higher band excitations but is fast compared to the trap frequencies. The dotted lines in Fig. 1c trace the opposite trajectories of and during tE. This quench process drives the system across the single minimum to double minima phase transition and generates spin-dependent synthetic electric fields  = . Consequently, atoms in different dressed spin components move off in opposite directions from the trap center (or from the region around q = 0 in the quasimomentum space as shown in Fig. 1c as dashed circles for two representative ΩF = 0,1.3 Er) and then undergo out-of-phase oscillations, thus exciting the SDM and an AC spin current. Approximately equal populations in the two dressed (or bare) spin components are maintained by keeping δR = 0 as Ω is changed from ΩI to ΩF (see Methods). After the application of E, the Raman coupling is maintained at ΩF during the hold time (thold). We then abruptly turn off both the Raman lasers and the dipole trap for time of flight (TOF) absorption imaging, measuring the bare spin and momentum composition of the atoms (Fig. 1b). Experiments are performed at various thold to map out the time evolution in the trap.

Measurements of the spin-dipole mode (SDM) and its damping

Figure 2 presents SDM measurements for a bare BEC (at ΩF = 0) and a dressed (or SO-coupled) BEC (at ΩF = 1.3 Er), with select TOF images taken after representative thold in the trap. Two TOF images labeled by thold = −1 ms are taken right before the application of E. In the bare case (Fig. 2a), the images taken at increasing thold show several cycles of relative oscillations (SDM) between the two spin components in the momentum space, accompanied by a notable reduction in the BEC fraction. We refer to the reduction of condensate fraction in this paper as thermalization. In the dressed case at ΩF = 1.3 Er (Fig. 2b), despite the fact that A are nearly the same as that for the bare case, the SDM is now strongly damped without completing one period. Besides, we observe higher BEC fraction remaining at the end of the measurement compared with the bare case. This can be seen in the narrower momentum distribution of thermal atoms with a more prominent condensate peak in Fig. 2b. From the TOF images, we fit the atomic cloud of each bare spin component (or dominant bare spin component of a dressed spin component) to a 2D bimodal distribution to extract the center-of-mass (CoM) momentum ħk↑(↓) or other (dressed) spin-dependent quantities (see Methods). The relative mechanical momentum between the two spin components in the SDM is then determined by ħkspin = ħ(k↑ − k↓).

Fig. 2

SDM of a bare or a dressed BEC. Select TOF images showing the bare spin and momentum compositions of atoms taken after applying spin-dependent synthetic electric fields E with ΩF = 0 (bare BEC) in a and ΩF = 1.3 Er (SO-coupled BEC) in b, followed by various hold times (thold) in the dipole trap. The TOF images labeled by thold = −1 ms are taken right before the application of E. The bare spin components (labeled by m, with in red and in blue) are separated along the horizontal axis. The vertical axis shows the atoms’ mechanical momentum ħk along the SOC direction . The color scale reflects the measured optical density (OD, see Methods). The total condensate atom number of the initial state at ΩI is Nc ~ (1–2) × 104 with trap frequencies ω ~ 2π × (37 ± 5) Hz and ω ~ ω ~ 2π × (205 ± 15) Hz. The TOF images (and associated analyzed quantities presented later) are typically the average of a few repetitive measurements

Figure 3a–e presents measurements of ħkspin versus thold at various ΩF. We see that the initial amplitude (2ħkr) of ħkspin is larger than the width of the atomic momentum distribution (<ħkr), and ħkspin damps to around zero at later times. The observed ħkspin as a function of thold is fitted to a damped sinusoid (see Methods) to extract the decay time constant τdamp. The SDM damping is then quantified by the inverse quality factor 1/Q = ttrap/(πτdamp), where 1/ttrap is the trap frequency along taking into account of the effective mass for the dressed case (see Methods). We observe that the damping (1/Q) is higher for larger ΩF, summarized by the purple data in Fig. 3f. Additionally, we have performed two control experiments, which suggest that SOC alone cannot cause momentum damping and thermalization if there are no collisions between the two dressed spin components. Only when there is SDM would notable thermalization be observed within the time of measurement. First, we measure the dipole oscillations[17,22] of a SO-coupled BEC with a single dressed spin component prepared in at various ΩF. This gives a spin current as well as a net mass current. We observe (e.g., Supplementary Fig. 1 in Supplementary Note 1) that these single-component cases exhibit very small damping (1/Q < 0.05, summarized by the red square data in Fig. 3f) and negligible thermalization. In another control experiment, we generate only an AC mass current without a spin current by exciting in-phase dipole oscillations of two dressed spin components of a SO-coupled BEC without relative collisions (SDM). This experiment also reveals very small damping and negligible thermalization (see Supplementary Fig. 2 in Supplementary Note 1).

Fig. 3

Momentum damping at different ΩF, for SDM and the dipole mode of a single dressed spin component. a–e Relative momentum oscillations in SDM, ħkspin, as a function of thold at various ΩF. The experimental data (scatters) are fitted to a damped sinusoidal function (line) to extract the inverse quality factor 1/Q of the oscillations. f Momentum damping (quantified by 1/Q) versus ΩF. The error bar of 1/Q is the standard error of the fit. The purple circle data correspond to the SDM (illustrated by g) and the red square data correspond to the dipole mode of a BEC with a single dressed spin component prepared in (illustrated by h). In g, h, the representative band structure is calculated at Ω = 1.0 Er

Thermalization and spin current

We now turn our attention to the thermalization, i.e., the reduction of condensate fraction due to collisions between the two spin components. To quantitatively describe the observed thermalization, the integrated optical density of the atomic cloud in each spin component is fitted to a 1D bimodal distribution to extract the total condensate fraction fc = Nc/N (see Methods) with N being the total atom number and Nc the total condensate atom number (including both spin states). The time (thold) evolution of the measured fc is plotted for the bare (ΩF = 0) and dressed (ΩF = 1.3 Er and 2.1 Er) cases in Fig. 4a. In all the cases, we observe that fc first decreases with time before it no longer changes substantially (within the experimental uncertainty) after some characteristic thermalization time (τtherm). To capture the overall behavior of the thermalization, we fit the smoothed thold-dependent data of fc to a shifted exponential decay fc(thold) = fs + (fi − fs)exp(−thold/τtherm), where τtherm represents the time constant for the saturation of the decreasing condensate fraction and fs the saturation condensate fraction (see Methods). We obtain τtherm = 3.8(4) ms, 2.4(3) ms, and 0.4(1) ms for ΩF = 0, 1.3 Er, and 2.1 Er, respectively. Besides, a notably larger condensate fraction (fs) is left for a larger ΩF, where fs ~ 0.2, 0.3, and 0.4 for ΩF = 0, 1.3 Er, and 2.1 Er, respectively. Since thermalization during our measurement time is induced by the SDM, the observation that a larger ΩF gives rise to a smaller τtherm and a larger fs (Fig. 4b) thus less thermalization is understood as due to the stronger SDM damping (smaller τdamp) at larger ΩF, stopping the relative collision between the two spin components thus the collision-induced thermalization earlier.

Fig. 4

Thermalization and spin current. a The measured condensate fraction fc = Nc/N as a function of thold for SDM in the bare case (no SOC, ΩF = 0) and the dressed cases (with SOC, ΩF = 1.3 Er and ΩF = 2.1 Er). Representative error bars show the average percentage of the standard error of the mean. The solid curves are the shifted exponential fits to the smoothed fc (see Methods). The initial condensate fraction (not shown) at ΩI (measured at thold = −1 ms) is ~0.6–0.7 for all the cases. b The saturation time constant τtherm of the decreasing fc and the saturation condensate fraction fs versus ΩF, where the vertical error bar is the standard error of the fit. c Spin current Is (normalized by v/λ = 7.4 × 103/s, where v ~ 6 mm/s is the recoil velocity) as a function of thold for ΩF = 0 and 1.3 Er. The solid curves are fits (see text)

The coherent spin current is phenomenologically defined as Is = I↑ − I↓ (see Methods), where I is given by:Here, σ labels the physical quantities associated with the spin component σ, L is the in situ BEC size along the current direction, and v = ħk/m. We exclude the contribution from the thermal atoms as only the condensate atoms participate in the coherent spin transport. In our experiments, N↑/L↑ ≈ N↓/L↓ is not observed to decrease significantly with thold, and , thus the relaxation of Is is mainly controlled by that of . Therefore, the SDM damping (reduction of v↑ − v↓) and thermalization (reduction of fc) provide the two main mechanisms for the relaxation of coherent spin current.

Figure 4c shows the normalized Is as a function of thold extracted (see Methods) for ΩF = 0 and 1.3 Er. In the bare case, the spin current oscillates around and decays to zero. In the dressed case, the spin current relaxes much faster to zero without completing one oscillation. Fitting Is versus thold to a damped sinusoidal function for ΩF = 0 or to an exponential decay for ΩF = 1.3 Er (with no observable Is oscillations) allows us to extract the spin current decay time constant τspin, which is 5.1(8) ms and 0.5(0) ms, respectively. In the dressed case Is decays much faster compared to the bare case because both τdamp and τtherm are much smaller due to stronger SDM damping. In the bare case, the thermalization plays a more important role in the relaxation of Is due to the larger reduction of condensate fraction (fi − fs) compared to the dressed case.

Observation of deformed atomic clouds and BEC shape oscillations

In addition to the SDM damping and thermalization, the atomic clouds can exhibit other rich dynamics after the application of E. We observe deformation of atomic clouds at early stages of the SDM, as shown in Fig. 5a–d. Figure 5b, d shows the observation of an elongated atomic cloud at thold = 0.5 ms in the dressed case at ΩF = 2.1 Er, in comparison with the atomic cloud at thold = 0.5 ms in the bare case shown in Fig. 5a, c. Figure 5c, d shows the integrated optical density (denoted by OD) of the atomic cloud versus the y direction, obtained by integrating the measured optical density over the horizontal direction in TOF images. The momentum distribution of the atoms at ΩF = 2.1 Er has lower OD and is more elongated without a sharp peak along the SOC direction, in comparison with the bare case that has higher OD and a more prominent peak momentum. Furthermore, we observe that the relaxation of the spin current is accompanied by BEC shape oscillations[52-54] (Fig. 5e, f), which remain even after the spin current is fully damped. These additional experimental observations are closely related to the spin current relaxation, as discussed below.

Fig. 5

Observation of deformed atomic clouds and BEC shape oscillations. a–d Observation of deformed atomic clouds at early stages of the SDM. a, b TOF images for ΩF = 0 and ΩF = 2.1 Er at thold = 0.5 ms are shown for comparison. The corresponding integrated optical density (OD) versus the momentum in the SOC direction for the spin down and up components is shown respectively in c, d. e, f Observation of BEC shape oscillations. The data showing the aspect ratio W/W (see Methods) of the condensate measured at various thold are extracted from the SDM measurements in Fig. 3, except for the additional measurements #2 and #3 in e. e For the three independent measurements in the bare case, the observed oscillations possess a complicated behavior without having a well-defined frequency given the error bars and the fluctuation in the data. Select TOF images for measurement #1 are shown in Fig. 2a. f In the dressed cases, aspect ratio oscillations with a well-defined frequency are observed in measurements at three different ΩF. The average frequency of the three aspect ratio oscillations obtained from the damped sinusoidal fit is around 58 Hz, consistent with the expected frequency for the m = 0 quadrupole mode for a cigar-shape BEC in the limit of ω/ω << 1[52]. Note that ω is not modified by Raman lasers and thus does not depend on ΩF. Select TOF images for ΩF = 1.3 Er are shown in Fig. 2b. The representative error bars in (e, f) are standard deviation of at least three measurements. The dashed lines indicate thold ~ 2τdamp after which the SDM is fully damped out

GPE simulations and interpretations

We have performed numerical simulations for the SDM based on the 3D time-dependent Gross-Pitaevskii equation (GPE), using similar parameters as in the experiments. The ΩF-dependent 1/Q extracted from the GPE-simulated SDM (Fig. 6a–c) shows qualitative agreement with the experimental measurements (Fig. 6d, e). Quantitatively, we notice that the GPE simulation generally underestimates the momentum damping compared to the experimental observation (Fig. 6e), especially at low ΩF (including the bare case). This is possibly related to the fact that our GPE simulation cannot treat thermalization (which is more prominent at low ΩF) and effects of thermal atoms. Nonetheless, the in situ (real space) spin-dependent density profiles (Fig. 6f–j) of the BECs calculated from the GPE simulations have provided important insights to understand why SOC can significantly enhance the SDM damping. Figure 6f shows that the initial BEC (just before applying E) in the trap is in an equal superposition of bare spin up and down states. Figure 6g–j shows the density profiles of the BECs at thold = 1.5 ms (after applying E) in the trap with four different ΩF (see Supplementary Movies 2, 4 and 5 in Supplementary Note 3). For the bare case, the two spin components fully separate in the real space within the trap. As ΩF becomes larger, we observe that only a smaller portion of atoms in each spin component is well separated, as marked by the white arrows. Concomitantly, a larger portion of atoms appears to get stuck around the trap center and form a prominent standing wave pattern, which we interpret as density modulations arising from the interference between the BEC wavefunctions of the two dressed spin components when and are no longer orthogonal in the presence of SOC (see Fig. 7a)[21,55-58]. Compared to the bare case, the formation of density modulations in the dressed case can lead to more deformed clouds in both the real and momentum spaces at early stages in the SDM, as revealed by the GPE simulations (Fig. 6a, b, f–j; Supplementary Movies 2, 4 and 5 in Supplementary Note 3). This is consistent with our experimental observation of a highly elongated momentum distribution of the atomic cloud along the SOC direction at early instants in the SDM of a SO-coupled BEC (Fig. 5b, d).

Fig. 6

GPE simulated SDM at various ΩF and the extracted SDM damping compared with experiment. a, b GPE simulations of the 1D momentum-space density distributions of the two bare spin components as a function of thold for the SDM at ΩF = 0 and ΩF = 1.3 Er, respectively. The 1D momentum density ρ(k) is obtained by integrating the 3D momentum density along k and k, i.e., . Then, these integrated 1D atomic momentum densities for sequential hold times (thold) are combined to show the atomic density in momentum space along the SOC direction versus thold. c GPE simulations of the SDM damping versus thold at various ΩF. The violet lines are the ħkspin (defined as the difference between the CoM momenta of the two spin components) as a function of thold for various ΩF. The CoM momentum (ħk↑,↓) of each bare spin component (at a given thold) is calculated by taking a density-weighted average of the corresponding 1D momentum density distributions such as those shown in a, b. The black lines are damped sinusoidal fits for the calculated ħkspin to extract the corresponding SDM damping (1/Q) which is shown in d along with the experimental data reproduced from Fig. 3f. e Replotting of d with 1/Q shown in logarithmic scale. f–j In situ (real space) atomic densities calculated from GPE simulations. f Initial in situ 2D density at Ω = ΩI (right before applying spin-dependent electric fields E). g–j In situ 2D density at thold = 1.5 ms (after the application of E) for ΩF = 0, 0.4 Er, 0.9 Er, and 1.3 Er, respectively. For f–j, the density is designated by brightness and the bare spin polarization by colors (red: ↓, blue: ↑, white: equal spin populations). The 2D densities ρ(x, y) in f–j are obtained by integrating the 3D atomic density along z, i.e., . In this figure, the simulations used the following parameters representative of our experiment: ΩI = 5.2 Er, δR = 0, Nc = 1.6 × 104, ω = 2π × 37 Hz, ω = ω = 2π × 205 Hz, tE = 1.0 ms

Fig. 7

Calculated nonorthogonality, effective interaction parameters, and immiscibility for two dressed spin states.

In (a–f), the calculations consider and located respectively at ħq and −ħq. a When Ω = 0, the nonorthogonality is zero because the two bare spin components are orthogonal. When Ω ≠ 0, either increasing Ω or decreasing q would increase 〈↑′|↓′〉, giving rise to stronger interference and more significant density modulations in the spatially overlapped region of the two dressed spin components. b, c Effective interspecies (g↑′↓′) and intraspecies (g↑′↑′, g↓′↓′) interaction parameters versus quasimomentum at Ω = 0.1 Er and 1.26 Er, respectively. When Ω increases or q decreases, g↑′↓′ increases while g↑′↑′ and g↓′↓′ almost remain at the bare values. As q → 0 at any finite Ω, g↑′↓′ → 2g↑′↑′ or 2g↓′↓′, which is the upper bound of g↑′↓′ (see Methods). The inset of b, c zooms out to show the maximum. d shows the immiscibility metric in Eq. (13) (see Methods) versus ħq corresponding to b. η < 0 means miscible, and η > 0 means immiscible. Over the range of plotted ħq, d can be miscible or immiscible depending on ħq. The inset of d zooms in to focus on the sign change of η. The vertical dotted line in (b–d) indicates ħq corresponding to the Ω in each case. The calculations are performed in the two-state picture described by Eq. (1) with δR = 0. e, f Immiscibility metric η versus Ω for various q. In e, as Ω becomes larger or q becomes smaller, the two dressed spin components can become more immiscible until η reaches the maximum value set by the upper bound of g↑′↓′ (see also b, c). f Zoom-in of e showing the miscible to immiscible transition (indicated by the gray dashed line at η = 0) as a function of Ω for various q. The red dot-dashed line corresponds to two dressed spin components located respectively at the band minima q, showing the well-known miscible to immiscible transition around 0.2 Er for a stationary SO-coupled BEC. In the dynamical case studied here, BECs can be located away from the band minima and approach q = 0, becoming immiscible even when Ω < 0.2 Er for small enough q

In addition to density modulations, our GPE simulation also reveals complex spatial modulation in the phase of the BEC wavefunctions (see Supplementary Fig. 9 and Supplementary Movies 3 and 6 in Supplementary Note 3). Such distortions of BEC wavefunctions in the amplitude (which determines the density) and the phase contribute to quantum pressure[59] and local current kinetic energy (see Methods) respectively, two forms of the kinetic energy that do not contribute to the global translational motion (or CoM kinetic energy) of each spin component. The sum of the CoM kinetic energy, quantum pressure, and local current kinetic energy is the total kinetic energy (see Methods). We have used GPE to calculate the time evolution of these different parts of kinetic energy for the dressed case, showing that the damping of the CoM kinetic energy (which decays to zero at later times) is accompanied by (thus likely related to) prominent increase of the quantum pressure and the local current kinetic energy (both remain at some notable finite values at later times) (see Fig. 8e–h). The increasing quantum pressure and local current kinetic energy may reflect the emergence of excitations that do not have the CoM kinetic energy. This is consistent with the experimentally observed generation of BEC shape oscillations (Fig. 5e, f), whose kinetic energy can be accounted for by the quantum pressure and the local current kinetic energy. Note that the excitation of BEC shape oscillations may also be understood by the observation of deformed clouds at early stages of the SDM (Fig. 5a–d), because the deformed shape of the BEC is no longer in equilibrium with the trap and thus initiates the shape oscillations. The observed BEC shape oscillations remain even after the SDM is completely damped in both bare and dressed cases. This indicates that the BECs are still excited even after the CoM relaxes to the single-particle band minima within the time of measurement.

Fig. 8

Time (thold) evolution of different forms of energies per particle at different ΩF as calculated by GPE. a The total energy is the sum of the total Raman energy, total potential energy, total interaction energy, and the total KE. The result in a confirms that the total energy is conserved during thold. b Total Raman energy. c Total potential energy. d Total interaction energy, sum of the bare interaction energies in i–k. e Total KE, sum of different types of kinetic energies in f–h. f Total CoM KE. g Total QP. h Total LC KE. i g↑↑ interaction energy. j g↓↓ interaction energy. k g↑↓ interaction energy

Discussion

Previous studies in stationary SO-coupled BECs (located at ground dressed band minima) have found that increasing Ω drives a miscible to immiscible phase transition at Ω ~ 0.2 Er due to the increased effective interspecies interaction (characterized by the interaction parameter g↑′↓′)[21,55-57,60]. In the miscible phase, the two dressed spin components have substantial spatial overlap, where density modulations form. It is important to note that the effective interactions, immiscibility and interference between the two dressed spin components depend on the quasimomentum (ħq) and ΩF (Fig. 7, see Methods for details). Therefore, in the dynamical case studied here, these properties vary with time and can be notably different from those in the stationary case. During the SDM, the two dressed spin components are forced to collide due to E. This can give rise to interference-induced density modulations in their spatially overlapped region even when they are immiscible. In addition, the BECs during the SDM can be located away from the band minima and approach q = 0. For the two dressed spin components with quasimomenta ±ħq, either increasing ΩF or decreasing |q| (towards 0) would increase 〈↑′|↓′〉 (Fig. 7a), giving rise to stronger interference and more significant density modulations. Such increased non-orthogonality between the two dressed spin states also notably increases the effective interspecies interaction (g↑′↓′) to become even larger than the effective intraspecies interactions (g↑′↑′ ≈ g↓′↓′) (Fig. 7b, c), enhancing further the immiscibility (Fig. 7d–f). For example, Fig. 7d shows the calculated immiscibility metric (see Methods), , versus ħq corresponding to Fig. 7b. Notice that when Ω is large enough, and can become immiscible in the whole range of quasimomentum that a BEC can access during the SDM. Figure 7e shows η versus Ω at various ħq. We see that as Ω becomes larger or q becomes smaller, the two dressed spin components can become more immiscible (i.e., η becomes more positive) until η reaches the maximum value set by the upper bound of g↑′↓′. Figure 7f zooms in the region of small Ω in Fig. 7e to focus on the sign change of η from negative to positive, which indicates the miscible to immiscible transition. Note that the red dot-dashed line (for q = q) corresponds to two dressed spin components located respectively at the band minima q, showing the well-known miscible to immiscible transition around 0.2 Er for a stationary SO-coupled BEC. In the dynamical case studied here, BECs can be located away from the band minima and approach q = 0, becoming immiscible even when Ω < 0.2 Er for small enough q.

We have performed several additional control GPE simulations, showing that the presence or the enhancement of any of these three factors can increase the damping of the relative motion between two colliding BECs: (1) interference (Supplementary Fig. 5 and Supplementary Movie 1 in Supplementary Note 3), (2) immiscibility (Supplementary Fig. 4 and Supplementary Table 1 in Supplementary Note 3), and (3) interactions (Supplementary Figs. 4, 6, 7 and 8 and Supplementary Table 1 in Supplementary Note 3), presumably by distorting the BEC wavefunctions (see Supplementary Movies 1–6 in Supplementary Note 3) irreversibly in the presence of interactions to decrease the CoM kinetic energy while increasing the quantum pressure and the local current kinetic energy. Therefore, enhanced immiscibility, interference, and interactions can all increase the damping of the SDM. For simulations in the absence of interactions, we do not observe irreversible damping within the simulation time of 100 ms (Supplementary Figs. 6–8 in Supplementary Note 3), suggesting that the interactions play an essential role for the damping mechanisms.

The physical mechanisms and processes revealed in our work may provide insights to understand spin transport in interacting SO-coupled systems. Our experiment also provides an exemplary study of the evolution of a quantum many-body system, including the generation and decay of collective excitations, following a non-adiabatic parameter change (quench). Such quench dynamics has been of great interest to study many outstanding questions in many-body quantum systems. For example, how such a system, initially prepared in the ground state but driven out of equilibrium due to a parameter quench that drives the system across a quantum phase transition, would evolve to the new ground state or thermalize has attracted great interests (see, e.g., a recent study where coherent inflationary dynamics has been observed for BECs crossing a ferromagnetic quantum critical point[61]). In our case, the sudden reduction of Ω in the Hamiltonian Eq. (1) excites the coherent spin current, whose relaxation is strongly affected by SOC and is related to the SDM damping as well as thermalization. Besides, the relaxation may be accompanied by the generation of other collective excitations such as BEC shape oscillations. Furthermore, compared to the bare case, the SOC-enhanced damping of the SDM notably reduces the collision-induced thermalization of the BEC, resulting in a higher condensate fraction left in the BEC. This condensate part exhibits a more rapid localization of its CoM motion, which may be more effectively converted to other types of excitations (associated with the SOC-enhanced distortion of the BEC wavefunctions). These features suggest that SOC opens pathways for our interacting quantum system to evolve that are absent without interactions, in our case providing new mechanisms for the spin current relaxation. Experiments on SO-coupled BECs, where many parameters can be well controlled in real time and with the potential of adding other types of synthetic gauge fields, may offer rich opportunities to study nonequilibrium quantum dynamics[62], such as Kibble–Zurek physics while quenching through quantum phase transitions[63], and superfluidity[16,33] in SO-coupled systems.

Methods

Spin-dependent vector potentials

In Eq. (1), the eigenenergies at δR = 0 are given by:For Ω < Ωc, the ground band of the energy-quasimomentum dispersion has two minima at:

The state of the atoms associated with each minimum at q can be regarded as a dressed spin state. For a double minima band structure, we thus have two dressed spin components and that constitute a pseudo spin-1/2 system (when Ω = 0, and become the bare spin and , respectively). The energy spectrum expanded around each q as E(q) = ħ2(q − q)2/(2m*) is analogous to the Hamiltonian describing a charged particle with charge Q in a magnetic vector potential A, , where m* is the effective mass of a dressed atom and m is the mass of the charged particle. Therefore, we can identify the quasimomentum (ħq) with the canonical momentum , and ħq with the light-induced spin-dependent vector potentials (A, by setting Q = 1 for our case[20]). The velocity operator corresponding to the mechanical momentum, , thus corresponds to ħ(q − q)/m*. These spin-dependent vector potentials A (represented by ħq) are tunable by Ω. For example, as seen in Fig. 1c, we can decrease Ω to separate the two ħq or increase Ω to combine them in the quasimomentum space.

Effects of the neglected state

We apply an external bias magnetic field (~5 gauss) to Zeeman split the energies E−1, E0, and E+1 of the , , and sublevels respectively (in the F = 1 hyperfine state of 87Rb atoms), where E−1 − E0 = ħωZ, E0 − E+1 = ħωZ − 2ε, ħ is the reduced Planck constant and ε = (E−1 + E+1)/2 − E0 is the quadratic Zeeman shift. The frequency difference between the two Raman lasers is ΔωR/(2π) = 3.5 MHz. The Raman detuning δR = ħ(ΔωR − ωZ) is controlled by B that controls ħωZ. In a first-order approximation, the third state can be excluded in Eq. (1) due to the quadratic Zeeman shift (2ε ~ 0.9 Er) from B but can be included in the following three-state Hamiltonian:

In our SDM experiments, we always maintain approximately equal spin populations in the and states both in the initial dressed state prepared at ΩI and in the final dressed state at ΩF (with approximately equal populations also achieved in and at ΩF). In Eq. (1) based on the two-state picture in the main text, δR = 0 can give rise to such balanced (dressed/bare) spin populations at any given Ω. However, in Eq. (5) with δR = 0, a finite Ω can lead to unbalanced (dressed/bare) spin populations. Therefore, in our experiment δR at a given Ω has to be changed to δ′(Ω,ε) to achieve the balanced spin populations (note that in the double minima regime of Eq. (5), this requirement is in a good approximation equivalent to the so-called balanced band condition where the two minima in the ground dressed band have equal energy). Such an effect is addressed in details in ref. [21]. In our case, also note that including the third state in Eq. (5) would cause the actual transition from the double minima to single minimum to occur at Ωc ~ 4.7 Er rather than at Ωc = 4.0 Er as would be predicted by Eq. (1). Additionally, Eq. (5) is used for plotting Figs. 1c and 3g, h, which more precisely means δR = δ′(Ω, ε) to achieve the balanced spin populations for the corresponding Ω. In the following, we use Eq. (5) to describe the initial state preparation process.

Initial state preparation, spin population balance, and imaging process

We create spin-polarized 87Rb BECs in in an optical dipole trap consisting of three cross laser beams (with a third beam added to the double beam dipole trap described in ref. [64]). To prepare the initial state of the BEC at the single minimum of the ground dressed band at ΩI = 5.2 Er (at δR = δ′(ΩI, ε), shown in Fig. 1c), first the Raman coupling Ω is ramped on slowly from 0 to ΩI in 80 ms (slow enough compared to the trap period and any inter-band excitation process) with δR ~ −ε in Eq. (5), such that the dominant bare spin component of the dressed BEC at any finite Ω during the ramping process remains in . Subsequently, while holding Ω at ΩI, we adjust B to change the Raman detuning from δR ~ −ε to δR = δ′(ΩI, ε) in 80 ms, and then we hold both Ω and δR for another 20 ms to let the system equilibrate. Note that adjusting δR to δ′(ΩI, ε) has to be empirically achieved by realizing the balanced spin populations, with the reasons addressed in the next paragraph. When the BEC is successfully prepared in the initial state at ΩI, equal populations in the and states can be achieved and seen in TOF images measured at thold = −1 ms.

In addition to the change in the band structure when going from the two-state picture to the three-state picture as discussed in the previous section, there are several other experimental factors that can lead to unbalanced spin populations. First, the slow drift in Ω can tilt (therefore unbalance) the band at a fixed δR. Second, a slow drift in B would give rise to a drift in δR. Third, sometimes there may still be excitations (for example, small-amplitude collective dipole oscillations of a dressed BEC) at the end of the initial state preparation[65], making the quasimomentum of the dressed BEC deviate slightly from the quasimomentum of the band minimum. As a result, the dressed BEC can have a nonzero group velocity and unbalanced spin populations at ΩI (before applying the spin-dependent electric fields E). Hence, this can lead to unbalanced spin populations after the application of E, and the spin polarization P of atoms is not maintained around zero during thold. Here, we define P = (N↑ − N↓)/(N↑ + N↓), where N↑(↓) is the total atom number of the atomic cloud (measured in the TOF images) for the bare spin component ↑(↓). Fourth, the quench process from the single minimum to double minima bands during tE (Fig. 1c) may also give rise to unbalanced spin populations, presumably because of the access to the magnetic phase in the double minima regime where the ground state is the occupation of a single dressed spin state (the two occupied dressed spin states are metastable states).

The above effects are avoided in our experiments by making sure that the balanced spin populations are empirically achieved throughout our experiment (with occasional adjustment of δR, and discarding runs with notably unbalanced spin populations). More specifically, we first make sure that balanced spin populations can be achieved at ΩI, assuring δR = δ′(ΩI, ε) after the initial preparation described above. Then, we linearly ramp δR from δ′(ΩI, ε) to δ′(ΩF, ε) as we change Ω from ΩI to ΩF in tE, and subsequently hold δR at δ′(ΩF, ε) for various thold. Here, δR = δ′(ΩF, ε) is empirically achieved by realizing balanced spin populations at Ω = ΩF for various thold. Therefore, when we state δR = 0 at a given Ω in the main text, it more precisely means that we realize balanced spin populations (as would be achieved at δR = 0 in the 2-state picture described by Eq. (1)).

The above-mentioned procedure of realizing δR = δ′(ΩF, ε) is further experimentally verified by observing balanced spin populations using the same bias magnetic fields but with tE = 15 ms and thold = 30 ms (slow enough to not to excite notable SDM). This suggests that such a choice of δR = δ′(ΩF, ε) approximates a balanced double minima band (with two equal-energy minima) at ΩF.

For the SDM measurements (e.g., Fig. 3), we make sure that the typical spin polarization is close to zero, with |P| = 0.05 ± 0.04, where 0.05 is the mean and 0.04 is the standard deviation of the data. Note that we used the total atom numbers N↑(↓) instead of condensate atom numbers to obtain P due to the less fluctuation in the fitted N↑(↓). Typically images with such small P, indicating good spin population balance for the whole atomic cloud, also do not exhibit notable spin population imbalance in their condensate parts.

After holding the atoms in the trap at ΩF for various thold, we turn off all lasers abruptly and do a 15-ms TOF, which includes a 9-ms Stern–Gerlach process in the beginning to separate the atoms of different bare spin states. Then, the absorption imaging is performed at the end of TOF to obtain the bare spin and momentum compositions of atoms. We then extract the physical quantities such as the mechanical momentum, condensate and thermal atom numbers of the atomic cloud in each spin state from such TOF images.

Analysis of momentum damping

Since the propagation direction of our imaging laser is ~27° with respect to the x-axis in the x − z plane (see Fig. 1a), the TOF images are in the y − z′ plane (where is perpendicular to in the x − z plane). The atomic cloud of each (dominant) bare spin component in the TOF images is fitted to a 2D bimodal distribution:where the first term corresponds to the condensate part according to the Thomas-Fermi approximation and the second term corresponds to the thermal part. Note that we only fit the majority bare spin cloud component when there is a distinguishable minority bare spin cloud component (which belongs to the same dressed spin state, but has a population <9% of the majority component in our experiments). This convention also applies to the analysis of the spin polarization defined above, condensate fraction, and the coherent spin current (see below). In the spin current or SOC directions , we obtain the relative mechanical momentum between the two bare spin components ħkspin = ħ(k↑ − k↓) from the difference between the center-of-mass positions of their condensate parts and the calibration of 2ħkr in TOF images (for example, 2ħkr can be calibrated from the distance between different bare spin components ↑ and ↓ that are in the same dressed spin state ↑′). To obtain the damping (1/Q) of the relative momentum oscillations in SDM (Fig. 3), ħkspin as a function of thold is fitted to a damped sinusoidal function , where τdamp is the momentum decay time constant. The data have a small offset B0 because we only use the majority bare spin component in each dressed spin component when extracting ħk↑,↓. We extract τdamp to obtain the inverse quality factor 1/Q = ttrap/(πτdamp), where is the trap period along the y direction taking into account of the effective mass meff for the dressed band around q, m is the bare atomic mass, and ω/(2π) is the trap frequency along the y direction in the absence of Raman lasers. Note that the effective masses around the two minima in the dressed ground band are nearly the same so we take their average as the meff. The standard error of the fit (95% confidence intervals) is obtained for determining the uncertainty of 1/Q shown in Fig. 3 in the main text.

For the dipole oscillations of a BEC with a single dressed spin component prepared in the state (see Supplementary Note 1), we fit ħk↓ (mechanical momentum of the dominant bare spin component ) as a function of thold to a damped sinusoidal function to extract τdamp and thus to obtain 1/Q. The minority bare spin component oscillates in phase with the dominant component with similar damping, and thus is not taken into account for determining 1/Q.

Analysis of condensate fraction

During the SDM, the atomic cloud can be significantly deformed along due to the interference between the two dressed spin components (see e.g., Fig. 5). Therefore, in order to extract the total condensate fraction (fc = Nc/N) of atoms to study the thermalization behavior as shown in Fig. 4a, b, the measured optical density (OD) of each bare spin component σ in the y − z′ plane is integrated along the y direction (the direction of SOC and the spin current as well as the direction along which the cloud can be significantly distorted) to obtain an integrated optical density versus z′ (denoted by OD). We fit OD of each bare spin component σ to a 1D bimodal distribution  + , where the first term corresponds to the condensate part according to the Thomas-Fermi approximation and the second term corresponds to the thermal part, to get the corresponding condensate and thermal atom numbers, and , respectively. The total condensate fraction is calculated as  = , shown as the scatters (unsmoothed raw data) in Fig. 4a.

To quantitatively describe the thermalization, we fit the smoothed total condensate fraction versus thold to a shifted exponential decay fc(t) = fs + (fi − fs)exp(−t/τtherm), where τtherm represents the time constant for the thermalization to stop and for the decreasing condensate fraction to saturate, with fs being the saturation condensate fraction. Because the large fluctuations in the unsmoothed data can give erroneous fitting results, each fitted curve shown as a solid line in Fig. 4a is the average of the three fits performed on the smoothed data, obtained using different levels (M = 1, 2, 3) of smoothing, where the smoothing is done by taking the average of the raw data within the nearest M time intervals.

Notice that the heating effect due to our Raman lasers (such as from spontaneous emission) is negligible within the time scale of the experiments (30 ms), because the lifetime of our BEC in the presence of the Raman lasers (with the Raman coupling considered in this work) is measured to be hundreds of ms. For example, the control experiment in Supplementary Fig. 1 shows no observable thermalization within 30 ms for dipole oscillations of a BEC with a single dressed spin component in the presence of the Raman lasers.

Coherent spin current

The I in Eq. (2) reflects the number of BEC atoms of a specific spin state passing through a cross section per unit time, and can be related to JA, where J = ncv is the current density along the SOC direction with the effective number density nc = Nc/(LA), v is the corresponding velocity, and A is an effective cross sectional area (the spin index σ is dropped in this discussion for simplicity in notations). The in situ length in the y direction, L, of each bare spin component can be estimated from the measured length of the BEC after TOF by  =  for a cigar-shape interacting BEC with ω >> ω and in the Thomas-Fermi approximation[66], where L(tTOF) is defined as 2R in Eq. (6) and L(tTOF = 0) = L. For example, for a typical L(tTOF = 15 ms) = 88 μm measured for one bare spin component of a dressed BEC prepared at ΩI, we get L = 4.5 μm for ω = 2π × 205 Hz. The two spin components have similar L when the spin populations are balanced. The in situ length L is thold-dependent during the dynamics and calculated from the thold-dependent TOF size, and is then used to obtain the thold-dependent spin current in Fig. 4b.

In the Thomas-Fermi approximation, we can also calculate L for the initial state at ΩI from the condensate atom number and trap frequencies. For example, we obtain L = 4.7 μm using Nc = 1.6 × 104 and ω = 2π × 205 Hz by , where , , a is the s-wave scattering length, and . In addition, the GPE-simulated L is 4.7 μm. These results of L are consistent with the value calculated from the TOF width.

Analysis of BEC shape oscillations

We characterize a condensate’s shape oscillations (in the y − z′ plane) of the bare spin component σ by its aspect ratio , where the condensate widths and (respectively along the y and z′ directions) are obtained from Eq. (6). We take the average of the aspect ratios of the two spin components , and plot W/W as a function of thold in Fig. 5. In Fig. 5e, caution has to be paid because the prominent thermalization in the bare case can make it challenging to fit the 2D cloud and extract the aspect ratio. The notable distortion of the cloud at the early stages of SDM can also make it difficult to perform the 2D Thomas–Fermi fit. Therefore, in Fig. 5f, we choose the thold-dependent W/W data after the corresponding dashed line (indicating thold ~ 2τdamp after which the SDM is fully damped) to fit to a damped sinusoidal function to extract the frequency of the aspect ratio oscillations.

In our experiments, there is no external modulation of the trapping potentials or shapes of the BECs to intentionally excite the shape oscillations. However, it is worth noting that shape oscillations can be induced via a non-adiabatic change in the internal energy of atomic clouds[67,68], which can take place when Ω is quickly changed or when the two spin components collide within the trap. On the other hand, we notice that in the dressed case the formation of density modulations can significantly deform the shape of a BEC (Fig. 5b, d; Supplementary Movies 2, 4 and 5 in Supplementary Note 3) and may thus also induce energetically allowed BEC shape oscillations, because the modified shape of the atomic cloud is no longer in equilibrium with the trap. Note that such a shape deformation can also change the internal energy. The m = 0 quadrupole mode excitation observed in our experiments has the lowest mode frequency among all possible quadrupole modes given our trap geometry and thus is the most energetically favorable (its mode frequency is also lower than the SDM frequency ~ω/(2π) for our trap parameters). Such nonresonant mode excitation is quite different from most previous studies, in which a collective mode of an atomic cloud is efficiently excited when it matches with the external modulation or perturbation of the trap[53,69] spatially and also spectrally (resonant with the modulation frequency). Compared to the dressed case, the bare case has less damped SDM and more significant thermalization, thus may complicate the shape oscillations due to more repeated SDM collisions and more atom loss[67,68]. We expect that the energy of the shape oscillations may eventually be converted to the energy of thermal atoms, leading to decay of the collective modes.

To further verify the excitation of the m = 0 quadrupole mode in the dressed case, we used another set of trap frequencies (see Supplementary Fig. 3 in Supplementary Note 2), and measured the condensate’s aspect ratio as a function of thold. The extracted frequency for the aspect ratio oscillations is again consistent with the predicted frequency for the m = 0 quadrupole mode.

Calculation of nonorthogonality, effective interaction parameters, and immiscibility

The interactions between atoms in bare spinor BECs are characterized by the interspecies (g↑↓, g↓↑) and intraspecies (g↑↑, g↓↓) interaction parameters, where , , c2 = −0.46a0, and c0 = 100.86a0 (a0 is the Bohr radius) for 87Rb atoms in our case. For a dressed BEC, in which is at some quasimomentum ħq (>0) and is at −ħq in the ground dressed band at Ω (in the two-state picture described by Eq. (1) in the main text with δR = 0), the effective interspecies (g↑′↓′ = g↓′↑′) and intraspecies (g↑′↑′, g↓′↓′) interaction parameters can be expressed in terms of the bare interaction g-parameters:where  = . The dressed spin states at −ħq and at ħq in the ground dressed band can be expressed asin the bare spin basis of . Using Eqs. (10) and (11), we can further obtainwhich characterizes the nonorthogonality (and thus the interference) between the two dressed spin states (where is located at ħq and is located at −ħq in the ground dressed band at Ω). Figure 7a plots such nonorthogonality versus quasimomentum for various Ω.

Note that (which is between 0 and π/2 in our case) characterizes the degree of bare spin mixing of a single dressed spin state (Eqs. (10) and (11)) as well as the nonorthogonality (due to the bare spin mixing, see Eq. (12)) between the two dressed spin states. As we can see, either decreasing Ω or increasing q would decrease (or increase ). When (or ), all the dressed spin states would approach the corresponding bare spin states, i.e., and , thus the nonorthogonality 〈↑′|↓′〉 → 0. In addition, all the effective interaction parameters would approach the corresponding bare values, i.e., g↑′↑′ → g↑↑, g↓′↓′ → g↓↓, and g↑′↓′ → g↑↓.

On the other hand, either increasing Ω or decreasing q would increase towards π/2 (or decrease ), thus enhancing the bare spin mixing, nonorthogonality and g↑′↓′. When (or ), , , and . Therefore, g↑′↓′ → 2g↑′↑′ or 2g↓′↓′, which is the upper bound of the effective interspecies interaction parameter. Figure 7b, c shows the effective interaction parameters normalized by g↑↑ versus quasimomentum ħq at Ω = 0.1 Er and Ω = 1.26 Er, respectively. When Ω increases or q decreases, g↑′↓′ increases while g↑′↑′ and g↓′↓′ almost remain at the bare values. As q → 0 at any finite Ω, g↑′↓′ approaches the upper limit 2g↑′↑′ or 2g↓′↓′.

In the case of SDM, assume that in the ground dressed band at Ω, is located at ħq and is located at −ħq at thold, we may use the immiscibility metric[70]to understand how Ω may modify the miscibility (η < 0) or immiscibility (η > 0) between and .

GPE simulations

The dynamical evolution of a BEC is simulated by the 3D time-dependent GPE[71]. To compare with the experimental data, we conduct simulations with similar parameters as those used in our experiment. The GPE of a SO-coupled BEC can be written in the following form:where is the momentum operator along , and HSOC is the (two-state) single-particle Hamiltonian Eq. (1), with q replaced by . Vtrap is the external trapping potential:where ω is the angular trap frequency along the spatial coordinate x(y, z). The wavefunction (order parameter) of a spinor BEC can be written in the formwhere ψ↓ and ψ↑ are the respective condensate wavefunctions of the two components, n↓(n↑) is the condensate density, ϕ↓(ϕ↑) is the phase of the wavefunction, r is the position, and t is time. The spatial integration of (n↓ + n↑) gives the total atom number N. The two-body interactions between atoms are described by the nonlinear interaction term Vint, which can be written in the basis of {ψ↓, ψ↑}:The interaction parameters are given byandThe spin-dependent s-wave scattering lengths for 87Rb atoms are c0 and c0 + c2, where c2 = −0.46a0 and c0 = 100.86a0 (a0 is the Bohr radius). The initial state of the SO-coupled BEC is obtained by using the imaginary time propagation method. Next we change ΩI to a final value ΩF in tE = 1.0 ms to simulate the spin-dependent synthetic electric fields. Equation (14) is used to simulate the dynamics of the BECs. The momentum space wavefunctions are calculated from the Fourier transformation of the real space wave functions. The squared amplitude of the momentum space wavefunctions is used to obtain the time-dependent momentum space density distributions shown in e.g., Fig. 6a, b.

For the GPE simulations in Fig. 6, we have checked that moderate variations in these parameters (as in our experimental data) do not affect our conclusions (while they can slightly change the 1/Q values, for example, larger 1/Q found for higher Nc). The simulations also reveal additional interesting features, such as the appearance of the opposite momentum (back-scattering) peak for each spin component in Fig. 6a, b, which are not well resolved in our experimental data.

Different forms of energies in GPE simulations

Using Eq. (16), the total energy density ε (the spatial integration of which gives the total energy of the system) can be expressed as the sum of several terms[35,59]:In the above equations, and . We will introduce ε1 to ε5 one by one in the following. The expression of ε1 in Eq. (21) is the density of the total (including two spin components) quantum pressure (QP), which is a type of kinetic energy (KE) associated with the spatial variation of the condensate density. An imaginary term appearing in the derivation of ε1 is not shown in Eq. (21) as its spatial integration (for a confined system) is zero and thus has no contribution to the energy. The expression of ε2 in Eq. (22) is the density of the sum of two types of KE, the total CoM KE (sum of the CoM KE of both bare spin components) and the total local current kinetic energy (LC KE). Both the CoM KE and LC KE are associated with the spatial variation of the phase of wavefunctions. The sum of the three types of kinetic energy (total QP, total CoM KE, and total LC KE) gives the total KE. That is, the sum of ε1 and ε2 is the density of the total KE. In the following, we derive explicit expressions for the CoM KE and LC KE. For CoM KE, it is nonzero only in the y direction because the SDM is along the y direction. Thus, the expression of CoM KE is:where is the momentum operator along for the spin down (up) component, and the last term in Eq. (27) is simply . Recall that ε2 in Eq. (22) is the density of the sum of CoM KE and LC KE. Thus, the expression of LC KE can be obtained by subtracting the expression of CoM KE in Eq. (27) from the spatial integration of ε2 (Eq. (22)):where Δ(∇ϕ↓,↑) is the standard deviation of ∇ϕ↓,↑, and note 〈∇ϕ↓,↑〉 = 0. Thus, if the wavefunction is a plane wave with a phase ϕ = qy, its LC KE is zero. For collective modes that do not have the CoM KE (for example, the quadrupole modes), the associated motional (kinetic) energy can be accounted for by LC KE and QP. The expression of ε3 in Eq. (23) is the density of the Raman energy, associated with the Raman coupling Ω. The expression of ε4 in Eq. (24) is the density of the sum of the bare intraspecies and interspecies interaction energies. The expression of ε5 in Eq. (25) is the density of the total potential energy.

To calculate the time (thold) evolution of the various forms of energies, we can in principle integrate the corresponding time-dependent energy densities over the real space. In practice, for the kinetic energy part we only perform spatial integration of ε2 (given by Eq. (22)). For convenience of computation, the total KE, total CoM KE, total LC KE, and total QP are calculated using a different approach taking advantages of the (quasi)momentum space representation of the quantum mechanical wavefunctions and operators. Specifically, the total KE is calculated by  +  in the quasimomentum space, where is the momentum operator along for the spin down (up) component, and ψ↓,↑(q, t) is the momentum-space representation of the wavefunctions (in the two directions not affected by SOC, x and z, we simply have q = p and q = p). Similarly, the total CoM KE is calculated in the quasimomentum space using Eq. (26). The total LC KE is calculated by subtracting the calculated total CoM KE from the spatial integration of ε2 (Eq. (22)). The total QP is calculated indirectly by subtracting the spatial integration of ε2 from the total KE.

The total Raman energy is calculated by the spatial integration of ε3 (Eq. (23)). The total bare intraspecies (g↑↑ and g↓↓) and interspecies (g↑↓) interaction energies are calculated by the spatial integration of the corresponding terms in ε4 (Eq. (24)). The total interaction energy is calculated as the sum of the bare intraspecies and interspecies interaction energies. The total potential energy is calculated by the spatial integration of ε5 (Eq. (25)). Lastly, the total energy of the system is calculated as the sum of the total Raman energy, total potential energy, total interaction energy, and total KE.

We note that even though our GPE simulations do not treat thermalization and thermal energies, the calculated different forms of condensate energies and their time evolution still provide valuable insights to understand the dynamical processes involved in the SDM. The GPE calculated different forms of energies shown in Fig. 8 and discussed in the associated texts below refer to the energies per particle (i.e., the calculated energies divided by the total atom number N).

In Fig. 8a, the total energy is a constant during thold, confirming the conservation of the total energy. In Fig. 8b, the total Raman energy has relatively small variations during thold. In Fig. 8c, the total potential energy in dressed cases has smaller variations during thold compared with that in the bare case. In Fig. 8d, the time evolution of the total interaction energy at different ΩF possesses a complicated behavior, mainly due to the complicated dynamics of the densities of the two spin components as well as their spatial overlap (see Supplementary Movies 2, 4 and 5 in Supplementary Note 3).

Figure 8e–h shows the time evolution of the calculated total KE, total CoM KE, total QP, and total LC KE at different ΩF, respectively. When ΩF is larger, the total CoM KE (Fig. 8f) exhibits a faster damping while QP as well as LC KE exhibit a faster increase (Fig. 8g, h, focusing on the relatively early stage of SDM) presumably due to the enhancement of the interference, immiscibility, and effective interaction between the two dressed spin components.

Figure 8i–k shows the time evolution of the calculated intraspecies and interspecies interaction energies at different ΩF. Note that the interaction energies are relatively small compared to other forms of energies, but are essential for the damping mechanisms as discussed in the main text.