Superfluid is a special kind of fluid which does not suffer dissipation when flowing through a tube. It loses its superfluidity only when its speed exceeds a certain critical value. The superfluidity of liquid helium 4 was first explained by Landau.^[10] He considered a superfluid moving inside a stationary tube with velocity v. Since the system is invariant under the Galilean transformation, this scenario is equivalent to a stationary fluid inside a moving tube. If the elementary excitation in a stationary superfluid with momentum q has energy ε ₀(q), then the energy of the same excitation in the background of a moving fluid with v is ε _v(q) = ε ₀(q) + v· q. A fluid experiences friction only through emitting elementary excitations, and it is a superfluid if these elementary excitations are energetically unfavorable. In other words, a superfluid satisfies the constraint ε _v(q) > 0. It readily leads to the well-known Landau’ s criterion for superfluid,

Here v_c is the critical velocity of the superfluid, which is determined by the smallest slope of the excitation spectrum of a stationary superfluid.

Another way of deriving the formula of critical velocity is from the point view of Cerenkov radiation. Consider a macroscopic impurity moving in the superfluid generates an excitation. According to the conservations of both momentum and energy, we should have

where m₀ is the mass of the impurity, v_i and v_f are the initial and final velocities of the impurity, respectively. The above two conservations (2) and (3) cannot be satisfied simultaneously when

The critical velocity v_c here has the same expression as that obtained from Landau’ s criterion. If the excitations are phonons, i.e., ε ₀(q) = c| q| , then v_c < c. This means that the impurity could not generate phonons in the superfluid and would not experience any viscosity when its speed was smaller than the sound speed. This is in fact nothing but the Cerenkov radiation, ^{[12– 14]} where a charged particle radiates only when its speed exceeds the speed of light in the medium. There is also interesting work discussing the possible deviation from this relation.^[16]

These two different ways of derivation are equivalent when the system has the Galilean invariance. By transforming to another reference frame illustrated in Fig.  1(b), the superfluid can be viewed as being dragged by a moving tube. We then replace the moving tube with a macroscopic impurity moving inside the superfluid as shown in Fig.  1(c). Consequently, for systems with Galilean invariance, the critical velocity of a superfluid moving inside a tube without experiencing friction is just the critical velocity of an impurity moving in a superfluid without generating excitations. For systems without the Galilean invariance, such an equivalence is lost as we will see later with the spin– orbit coupled BEC.

Fig.  1.

	Figure Option ViewDownloadNew Window
	Fig.  1. (a) A superfluid moves inside a stationary tube. (b) The superfluid is dragged by a tube moving at the speed of v. (c) An impurity moves at v in the superfluid. The two-way arrow indicates the equivalence between different scenarios.

For ultracold bosonic gases, the low-energy excitation is phonon, linear with respect to momentum q, thus the critical speed of this superfluid is just the sound speed. Taking into account the non-uniformity of the trapped gas, the critical velocity measured in experiments agrees well with the value predicted by Landau’ s theory.^[17] For liquid helium 4, as there is another kind of elementary excitations called rotons, the critical speed of superfluid helium 4 is largely determined by the roton excitation, much smaller than its sound speed. Nevertheless, the critical velocity measured in experiments is still one order of magnitude smaller than the value predicted by the theory.^[18] Therefore, it is remarkable that Landau’ s prediction of critical velocity was experimentally confirmed with BEC almost six decades after its invention.

One remark is warranted on Landau’ s theory of superfluidity. Landau’ s criterion  (1) of critical speed does not apply for many superfluids. However, Landau’ s energetic argument for superfluidity is very general and can be applied to all the cases considered in this review. We shall use this argument to determine the critical speeds of various superfluids.

As liquid helium 4 is a strongly-interacting system, the theoretical calculation of its excitation spectrum is challenging. However, for a dilute ultracold bosonic gases, convenient yet precise approximations can be made to determine its excitation in theory. Due to the dilute nature and short-range interaction, the interaction between atoms can be approximated by a contact interaction. Furthermore, for low energy scattering of bosonic atoms, the s-wave scattering channel dominates. Hence, as a good approximation, the actual complex interaction between atoms is replaced by an effective s-wave contact interaction, i.e., U(r₁, r₂) = 4π ħ ²aδ (r₁ − r₂)/M, where a is the s-wave scattering length. At zero temperature, assume all the particles condense into the same orbit ψ (r), then its evolution is governed by the mean field Gross– Pitaevskii (GP) equation, ^[19]

where V(r) is the trapping potential or other external potential, and c = 4π ħ ²a/M is the interaction parameter or coupling constant. The GP energy functional reads

Note that here the wave function ψ is normalized to the total particle number N. The approximation we make here is the mean-field approximation, and the interaction between atoms is replaced by an effective mean-field potential. The validity of the approximation usually depends on the condensed fraction. For weakly interacting dilute bosonic gases near zero temperature, the condensed fraction can be more than 90 percent. In these situations, the GP equation works well. It has been used to calculate the collective excitations of trapped BEC, as well as vortex dynamics. Pretty good agreement is achieved between theory and experiment. With a stochastic term describing the effect of thermal atoms, the modified stochastic GP equation can also simulate BEC systems at finite temperature.^{[20, 21]} In addition, for atoms with spin, one can also derive a multi-component GP equation to describe a spinor superfluid.^[22]

The stationary state the GP equation describes is a zeroth-order approximation in some sense. It only takes into account the interaction between the condensed particles. The first-order approximation is to take into account the interaction between the condensed particles and un-condensed particles. This is settled by the Bogoliubov theory of elementary excitations.^[23] The standard Bogoliubov theory is to diagonalize the mean-field Hamiltonian by the Bogoliubov transformation. Here we adopt an equivalent yet more convenient way to deal with this problem in inhomogeneous systems. For a stationary state ψ (r, t) that satisfies the GP equation, we consider some small time-dependent perturbations δ ψ (r, t) added to this stationary state. The perturbed state Ψ = ψ + δ ψ also satisfies the GP equation,

Expanding this equation to first order in the perturbation δ ψ , one arrives at the equation of motion of the perturbation δ ψ (r, t),

and its complex conjugate partner. Assume the state before perturbation is a stationary state with the wave function , and write the perturbations as , then we arrive at a pair of equations that u(r) and v(r) satisfy,

where the matrix 𝓜 is given by

with H₀ = − ħ ²∇ ²/2M + V(r) − μ + 2cn(r). By diagonalizing 𝓜 , one obtains two sets of solutions, but only the solution satisfying the constraint represents physical excitations. Note that since the characteristic matrix 𝓜 is not hermitian, its eigenvalues are not necessarily real. If its eigenvalue ε = ħ ω has nonzero imaginary part, the perturbation will grow exponentially with time, signaling that the state before perturbation suffers dynamical instability. If its eigenvalue is real but negative, elementary excitations associated with the perturbations will be energetically favorable and superfluidity is lost, which is called Landau instability. For the excitations of repulsive Bose gases without external potentials, the low-energy excitations must be non-negative and gapless from the Hugenholtz– Pines theorem.^[24]

As a specific example, we apply the above formalism to the uniform Bose gas, namely, putting V(r) = 0. The wave function before perturbation is just , and the perturbation has this form , and . Plugging these wave functions into the above equations, one immediately obtains the excitation spectrum for a uniform Bose gas,

where is the single particle spectrum. At long-wavelength or small momentum q, the excitation has the asymptotic form of This is nothing but phonon excitation, and is just the speed of sound. From Landau’ s theory of critical velocity, we conclude that the critical velocity of a uniform superfluid Bose gas is just the sound speed.

The method of mean-field approximation and Bogoliubov transformation is very general, and applies in other more complicated situations. In the following discussion, we use this method to study three types of superfluids, superfluid in a periodic potential, superfluid with spin– orbit coupling and superfluid of pure spin current.

Now we study the superfluidity of a BEC in a periodic potential, ^{[29, 30]} which is provided by the optical lattice in cold atom experiments. For simplicity, we consider a quasi-one-dimensional BEC, confined in a cigar-shaped trap. We treat the system with the mean-field theory and obtain the grand-canonical GP Hamiltonian as

where all the variables are scaled to be dimensionless with respect to a set of characteristic parameters of the system, the atomic mass M, the wave number k_L of the laser light generating the optical lattice, and the average density n₀ of the BEC. The chemical potential μ and the strength of the periodic potential v are in units of , the wave function ψ is in units of x is in units of k_L/2, and t is in units of . The interaction constant is given by , where a_s > 0 is the s-wave scattering length.

For non-interacting case (c = 0), diagonalizing the Hamiltonian will give the standard Bloch waves and energy bands. When the mean-field interaction is turned on (c ≠ 0), in principle the Hamiltonian allows for other types of solutions which have no counterpart in the non-interacting case.^{[31, 32]} Here we focus on the solutions which still have the form of Bloch waves, i.e., ψ _k(x) = e^ikxϕ _k(x), where ϕ _k(x) has the same period with the optical lattice. ϕ _k(x) can be found by extremizing the Hamiltonian above.^[30] The solution found in this way should satisfy the stationary GP equation with periodic potential,

To determine the superfluidity of these Bloch states, we must consider elementary excitations around these Bloch states, and check whether the excitation energy is always positive. Positive excitation energy indicates that the Bloch state is a local energy minimum, and it is stable against small perturbations. Due to the periodicity of the Bloch wave, the perturbations can be decomposed into different decoupled modes labeled by q,

where q ranges between − 1/2 and 1/2, and the perturbation functions u_k and v_k are of periodicity of 2π .

Following the similar method in Section  2.2, we linearize the GP equation above to obtain the Bogoliubov equation that u_k and v_k satisfy

where

with

This eigenvalue equation has two sets of solutions, with one corresponding to physical excitations, which is mostly phonon excitation, and the other called anti-phonon that is not physical. If the physical excitation ε _k(q) is positive, the Bloch wave ψ _k is a local minimum, and the system will have superfluidity. Otherwise, the system suffers Landau instability or dynamical instability, depending on whether ε _k(q) is real negative or complex.

In the case of free space v = 0, the Bloch state ψ _k becomes a plane wave e^ikx. Then the operator 𝓜 _k(q) becomes

and we recover the excitations in the uniform case

We immediately see that the excitation energy is always real, which means that the BEC flows in free space are always dynamically stable.

When there is a periodic potential, the situation is dramatically different, where the excitation energy can have imaginary part, signaling the dynamical instability of the system. By numerically solving the Bogoliubov equation above, we show the stability phase diagrams for BEC Bloch waves in the panels of Fig.  2, where different values of v and c are considered. The results have reflection symmetry in k and q, so we only show the parameter region, 0 ≤ k ≤ 1/2 and 0 ≤ q ≤ 1/2. In the shaded area (light or dark) of each panel of Fig.  2, the excitation energy is negative, and the corresponding Bloch states ψ _k are saddle points. For those values of k outside the shaded area, the Bloch states are local energy minima and represent superfluids. The superfluid region expands with increasing atomic interaction c, and occupies the entire Brillouin zone for sufficiently large c. On the other hand, the lattice potential strength v does not affect the superfluid region very much as we see in each row. The phase boundaries for v ≪ 1 are well reproduced from the analytical expression for v = 0, which is plotted as triangles in the first column.

Fig.  2.

Figure Option ViewDownloadNew Window

Fig.  2. Stability phase diagrams of BEC Bloch states in optical lattices. k is the wave number of BEC Bloch waves; q denotes the wave number of perturbation modes. In the shaded (light or dark) area, the perturbation mode has negative excitation energy; in the dark shaded area, the mode grows or decays exponentially in time. The triangles in (a1– a4) represent the boundary, q2/4 + c = k2, of saddle point regions at v = 0. The solid dots in the first column are from the analytical results of Eq.  (20). The circles in (b1) and (c1) are based on the analytical expression  (21). The dashed lines indicate the most unstable modes for each Bloch wave ψ k.[29]

If ε _k(q) is complex, the system suffers dynamical instability, which is shown by the dark-shaded areas in Fig.  2. The dynamical instability is the result of the resonance coupling between a phonon mode and an anti-phonon mode by first-order Bragg scattering. The matrix 𝓜 _k(q) is real in the momentum representation, meaning that its complex eigenvalue can appear only in conjugate pairs and they must come from a pair of real eigenvalues that are degenerate prior to the coupling. Degeneracies or resonances within the phonon spectrum or within the anti-phonon spectrum do not give rise to dynamical instability; they only generate gaps in the spectra. Based on this general conclusion, we consider two special cases, which allow for simple explanations of the onset of dynamical instability.

One case is the weak periodic potential limit v ≪ 1, where we can approximate the boundary with the free space case. This case corresponds to the first column of Fig.  2. In this limit, we can approximate the phonon spectrum and the anti-phonon spectrum with the ones given by Eq.  (19). By equating them, ε ₊ (q − 1) = ε ₋ (q). For the degeneracy, we find that the dynamical instability should occur on the following curves:

These curves are plotted as solid dots in Fig.  2, and they fall right in the middle of the thin dark-shaded areas. To some extent, one can regard these thin dark-shaded areas as broadening of the curves  (20). It is noted in Ref.  [33] that the relation  (20) is also the result of ε ₊ (q − 1) + ε ₊ (− q) = 0, which involves only the physical phonons. Therefore, the physical meaning of Eq.  (20) is that one can excite a pair of phonons with zero total energy and with total momentum equal to a reciprocal wave number of the lattice.

The other case, c ≪ v, is shown in the first row of Fig.  2. The open circles along the left edges of these dark-shaded areas are given by

where E₁(k) is the lowest Bloch band of

In this linear periodic system, the excitation spectrum (phonon or anti-phonon) just corresponds to transitions from the Bloch states of energy E₁(k) to other Bloch states of energy E_n(k + q), or vise versa. The above equation is just the resonance condition between such excitations in the lowest band (n = 1). Alternatively, we can write the resonance condition as

Thus, this condition may be viewed as the energy and momentum conservation for two particles interacting and decaying into two different Bloch states E₁(k + q) and E₁(k − q). This is the same physical picture behind Eq.  (20).

One common feature of all the diagrams in Fig.  2 is that there are two critical Bloch wave numbers, k_t and k_d. Beyond k_t the Bloch waves ψ _k suffer the Landau instability; beyond k_d the Bloch waves ψ _k are dynamically unstable. The onset of instability at k_d always corresponds to q = 1/2. In other words, if we drive the Bloch state ψ _k from k = 0 to k = 1/2 the first unstable mode is always q = ± 1/2, which represents period doubling. Only for k > k_d can longer wavelength instabilities occur. The growth of these unstable modes drives the system far away from the Bloch state and spontaneously breaks the translational symmetry of the system.

Besides inducing the dynamical instability, the presence of the optical lattice has also non-trivial consequences on the concept of critical velocity. In contrast to the homogeneous superfluid which has only one critical velocity, there are two distinct critical velocities for a periodic superfluid.^[34] The first one, which we call inside critical velocity, is for an impurity to move frictionlessly in the periodic superfluid system (Fig.  3(a)); the second, which is called trawler critical velocity, is the largest velocity of the lattice for the superfluidity to maintain (Fig.  3(b)). We illustrate these two critical velocities with a BEC in a one-dimensional optical lattice.

Fig.  3.

	Figure Option ViewDownloadNew Window
	Fig.  3. (a) A macroscopic obstacle moves with a velocity of v inside a superfluid residing in a periodic potential. The superfluid and periodic potential are “ locked” together and there is no relative motion between them. (b) The lattice where a superfluid resides is slowly accelerated to a velocity of v.

The presence of the optical lattice plays a decisive role in the appearance of the two critical velocities: two very different situations can arise. The first situation is described in Fig.  3(a), where one macroscopic impurity moves inside the superfluid. The key feature in this situation is that there is no relative motion between the superfluid and the lattice. The other situation is illustrated in Fig.  3(b), where the lattice is slowly accelerated to a given velocity and there is a relative motion between the superfluid and the lattice. For these two different situations, two critical velocities naturally arise.

In the first scenario, we consider a moving impurity that generates an excitation with momentum q and energy ε ₀(q) in the BEC. According to the conservations of both momentum and energy, we should have

where m₀ is the mass of the impurity, v_i and v_f are the initial and final velocities of the impurity, respectively, G is the reciprocal vector, and ε ₀(q) is the excitation of the BEC at the lowest Bloch state k = 0. Note that in contrast to the conservation of momentum in free space, here the total momentum of the impurity and the excitation is not exactly conserved due to the presence of an optical lattice, and the momenta differing by integer multiples of reciprocal lattice vector are equivalent. For phonon excitations, i.e., ε ₀(q) = c| q| , these two conservations can always be satisfied simultaneously no matter how small the velocity of the impurity is. In other words, the critical velocity for this scenario is exactly zero.

In the other scenario, there is relative motion between the superfluid and the optical lattice, and the superfluid no longer resides in the k = 0 point of the Brillouin zone. We should examine the stability of Bloch waves with nonzero k, which is discussed in detail in the previous subsection. The critical wave number k_t mentioned above corresponds precisely to the trawler critical velocity v_t here. As k_t is not zero, v_t is not zero.

Both critical velocities can be measured with BECs in optical lattices. The inside critical velocity v_i can be measured with the same experimental setting as that in Ref.  [35], where the superfluidity of a BEC was studied by moving a blue-detuned laser inside the BEC. For the trawler critical velocity v_t, one can repeat the experiment in Ref.  [28], where a BEC is loaded in a moving optical lattice. One only needs to shift his attention from dynamical instability to superfluidity. The potential difficulty lies in that the Landau instability occurs over a much larger time scale, which may be beyond the life time of a BEC.^[28]

We use the method introduced in Section  2.2 to study the superfluidity of a BEC with SOC by computing its elementary excitations.^[57] In experiments, only the equal combination of Rashba and Dresselhaus coupling is realized. Here we use the Rashba coupling as an example as the main conclusion does not rely on the details of the SOC type.

We calculate how the elementary excitations change with the flow speed and manage to derive from these excitations, the critical speeds for the two different scenarios shown in Fig.  4(a) and 4(b). We find that there are two branches of elementary excitations for a BEC with SOC: the lower branch is phonon-like at long wavelengths and the upper branch is generally gapped. Careful analysis of these excitations indicates that the critical flowing velocity for a BEC with SOC (Fig.  4(a)) is non-zero, while the critical dragging speed is zero (Fig.  4(b)). This shows that the critical velocity depends on the reference frame for a BEC with SOC and, probably, for any superfluid that has no Galilean invariance.

Specifically, we consider a BEC with pesudospin 1/2 and the Rashba SOC. The system can be described by the Hamiltonian^{[44, 47, 58, 59]}

where γ is the SOC constant, C₁ and C₂ are interaction strengths between the same and different pesudospin states, respectively. For simplicity and easy comparison with previous theory, we focus on the homogeneous case V(r) = 0 despite that the BEC usually resides in a harmonic trap in experiments. Besides, we limit ourselves to the case C₁ > C₂, namely, in the plane wave phase. In the following discussion, for simplicity, we set ħ = M = 1 and ignore the non-essential z direction, treating the system as two-dimensional. We also assume the BEC moves in the y direction, and the critical velocity is found to be not influenced by the excitation in the z direction.

The GP equation obtained from the Hamiltonian  (32) has plane wave solutions

where tanθ _k = k_x/k_y and μ (k) = | k| ²/2 − γ | k| + (C₁ + C₂)/2. The solution ϕ _k is the ground state of the system when | k| = γ . There are another set of plane wave solutions, which have higher energies and are not relevant to our discussion.

We study first the scenario depicted in Fig.  4(a), where the BEC flows with a given velocity. Since the system is not invariant under the Galilean transformation, we cannot use Laudau’ s argument to find the excitations for the flowing BEC from the excitation of a stationary BEC. We have to compute the excitations directly. This can be done by computing the elementary excitations of the state ϕ _k with the Bogoliubov equation for different values of k.

Without loss of generality, we choose k = kŷ with k > 0. Following the standard procedure of linearizing the GP equation, ^{[29, 30]} we have the following Bogoliubov equation:

where

with

and

As usual, there are two groups of eigenvalues and only the ones whose corresponding eigenvectors satisfy | u_i| ² − | v_i| ² = 1 (i = 1, 2) are physical.

In general there are no simple analytical results. We have numerically diagonalized 𝓜 to obtain the elementary excitations. We find that part of the excitations are imaginary for BEC flows with | k| < γ . This means that all the flows with | k| < γ are dynamically unstable and therefore do not have superfluidity. For other flows with k ≥ γ , the excitations are always real and they are plotted in Fig.  5. One immediately notices that the excitations have two branches, which contact each other at a single point. Closer examination shows that the upper branch is gapped in most of the cases, while the lower branch has phonon-like spectrum at large wavelengths. These features are more apparent in Fig.  6, where only the excitations along the x axis and y axis are plotted.

Fig.  5.

	Figure Option ViewDownloadNew Window
	Fig.  5. Elementary excitations of a BEC flow with the SOC. (a) k = 1; (b) k = 3; (c) k = 4. C1 = 10, C2 = 4, and γ = 1.

Fig.  6.

	Figure Option ViewDownloadNew Window
	Fig.  6. Excitations along the x axis (the first column) and y axis (the second column) at different values of k. (a1, a2) k = 1; (b1, b2) k = 3; (c1, c2) k = 4. C1 = 10, C2 = 4, and γ = 1.

In Figs.  5(c) and 6(c2), we note that part of the excitations in the upper branch are negative, indicating that the underlying BEC flow is thermodynamically unstable and has no superfluidity. In fact, our numerical computation shows that there exists a critical value k_c: when k > k_c either part of the upper branch of excitations or part of the lower branch or both become negative. This means that the flows described by the plane wave solution ϕ _{k, −} with | k| > k_c suffer Landau instability and have no superfluidity. Together with the fact that the flows with | k| < γ are dynamically unstable, we can conclude that only the flows with γ ≤ | k| ≤ k_c have superfluidity. The corresponding critical speed is v_c = k_c − γ . We have plotted how the critical flowing velocity v_c varies with the SOC parameter γ (Fig.  7).

Fig.  7.

	Figure Option ViewDownloadNew Window
	Fig.  7. Critical flowing velocity vc of a BEC as a function of the SOC parameter γ . (a) C1 = 11, C2 = 3; (b) C1 = 14, C2 = 6.

We turn to another reference frame illustrated in Fig.  4(b), where the BEC can be viewed as being dragged by a moving tube. To simplify the discussion, we replace the moving tube with a macroscopic impurity moving inside the BEC as shown in Fig.  4(c). Correspondingly, the question “ whether the BEC will be dragged along by the moving tube?” is replaced by an equivalent question “ whether the impurity will experience any viscosity?” . Suppose that the moving impurity generates an excitation in the BEC. According to the conservations of both momentum and energy, we should have

as those for a BEC without SOC, where ε ₀(q) is the excitation of the BEC at k = γ . The critical dragging velocity derived from Eqs.  (36) and (37) is given by

If the excitations were pure phonons, i.e., ε ₀(q) = c| q| , these two conservations would not be satisfied simultaneously when v ≈ | v_i| ≈ | v_f| < c. This means that the impurity could not generate phonons in the superfluid and would not experience any viscosity when its speed was smaller than the sound speed. Unfortunately, for our BEC system, the elementary excitations ε ₀(q) are not pure phonons, as will be shown below.

When γ ≠ 0, the excitations ε ₀(q) also share two branches. Along the x axis, these two branches are

where s₁ = 2γ ⁴ + γ ²(C₁ − C₂), , , t₂ = 2s₁s₂, and t₃ = 2s₁ + (γ ² + C₂/2)². Along the y axis, the excitations of the ground state are

When γ > 0, the upper branch is always parabolic at small q_x with a gap . When expanded to the second order of q_x, the lower branch has the following form:

This shows that is parabolic at long wavelengths instead of linear as usually expected for a boson system. This agrees with the results in Ref.  [49]. This parabolic excitation has a far-reaching consequence: according to Eq.  (38), the critical dragging velocity v_c is zero, very different from the critical flowing velocity for a BEC moving in a tube. This shows that the critical velocity for a BEC with SOC is not independent of the reference frame, in stark contrast with a homogeneous superfluid without SOC. This result of course has the root in the fact that the BEC described by the SOC Hamiltonian  (32) is not invariant under the Galilean transformation.^[56]

We have also investigated the superfluidity with the general form of SOC, which is a mixture of Rashba and Dresselhaus coupling. Mathematically, this SOC term has the form α σ _xp_y − β σ _yp_x. The essential physics is the same: the critical flowing speed is different from the critical dragging speed, and therefore the critical velocity depends on the choice of the reference frame. However, the details do differ when α ≠ β . The critical dragging speed is no longer zero. Without loss of generality, we let α > β . The slope of the excitation spectrum for the ground state along the y axis is

and the slope along the x axis,

The critical dragging velocity is the smaller one of the above two slopes, both of which are nonzero.

Spin– orbit coupled BECs have been realized by many different groups^{[39– 42]} through coupling ultracold ⁸⁷Rb atoms with laser fields. The strength of the SOC in the experiments can be tuned by changing the directions of the lasers^{[39– 41]} or through the fast modulation of the laser intensities.^[60] The interaction between atoms can be adjusted by varying the confinement potential, the atom number or through the Feshbach resonance.^[6] For the scenario in Fig.  4(b), one can use a blue-detuned laser to mimic the impurity for the measurement of the critical dragging speed similar to the experiment in Ref.  [61]. For the scenario in Fig.  4(a), there are two possible experimental setups for measuring the critical flowing speed. In the first one, one generates a dipole oscillation similar to the experiment in Ref.  [41] but with a blue-detuned laser inserted in the middle of the trap. The second one is more complicated. At first, one can generate a moving BEC with a gravitomagnetic trap.^[62] One then uses Bragg spectroscopy^{[63, 64]} to measure the excitations of the moving BEC, from which the superfluidity can be inferred. For the typical atomic density of 10¹⁴– 10¹⁵  cm^{− 3} achievable in current experiments, ^[61] and the experimental setup in Ref.  [39], the critical flowing velocity is 0.2– 0.6  mm/s, while the critical dragging velocity is still very small, about 10^{− 3}– 10^{− 2}  mm/s. To further magnify the difference between the two critical velocities, one may use the Feshbach resonance to tune the s-wave scattering length.

We consider a spin-1 BEC in free space. The mean-field wave function of such a spin-1 BEC satisfies the following GP equation:^[82]

where ψ _m (m = 1, 0, − 1) are the components of the macroscopic wave function, is the total density, is the spin density vector, and S = (S_x, S_y, S_z) is the spin operator vector with S_i (i = x, y, z) being the three Pauli matrices in the spin-1 representation. The collisional interactions include a spin-independent part c₀ = 4π ħ ²(a₀ + 2a₂)/3M and a spin-dependent part c₂ = 4π ħ ²(a₂ − a₀)/3M, with a_f (f = 0, 2) being the s-wave scattering length for spin-1 atoms in the symmetric channel of total spin f.

We consider a spin current state of the above GP equation with the form

where n is the density of the uniform BEC. The wave function above describes a state with component m = 1 moving at speed ħ k₁/M, component m = − 1 moving at speed ħ k₂/M, and component m = 0 stationary. The requirement of equal chemical potential leads to | k₁| = | k₂| . In the case where k₁ = − k₂, this state carries a pure spin current: the total mass current is zero as it has equal mass counterflow while the spin current is nonzero.

To determine whether the state  (46) represents a superfluid, we need to compute its Bogoliubov excitation spectrum, also using the method introduced in Section  2.2. It is instructive to first consider the special case where there is no counterflow, i.e., k₁ = k₂ = 0. The excitation spectra are found to be , and , respectively, with ε _q = ħ ²q²/2M. Hence for antiferromagnetic interaction (c₀ > 0, c₂ > 0), all branches of the spectra are real and there is a double degeneracy in one branch of the spectra. The phonon excitations give two sound velocities, (i = 0, 2), corresponding to the speeds of density wave and spin wave, respectively. However, the existence of phonon excitation does not mean that the pure spin current (k₁ = − k₂ ≠ 0) is a superflow as we cannot obtain a current with k₁ = − k₂ ≠ 0 from a state with k₁ = k₂ = 0 by a Galilean transformation.

For the counterflow state with k₁ = − k₂ ≠ 0, the stability has been studied in Ref.  [65]. It is found that, for the antiferromagnetic interaction case (c₀ > 0, c₂ > 0), the excitation spectrum of the m = 0 component always has nonzero imaginary part in the long-wavelength limit as long as there is counterflow between the two components, and the imaginary excitations in the m = 1, − 1 components only appear for a large enough relative velocity . For the ferromagnetic interaction case (c₀ > 0, c₂ < 0), both excitation spectra of the m = 0 and m = 1, − 1 components have nonzero imaginary parts for any relative velocity. This means that the pure spin current cannot be stable in any cases.

For the general non-collinear case (k = (k₁ + k₂)/2 ≠ 0) and antiferromagnetic interaction, the excitation spectrum for the m = 0 component is found to be

We see here that as long as the momenta of the two components are not exactly parallel, i.e., k₁ is not exactly equal to k₂, then | k| < | k₁| , and there is always dynamical instability for the long-wavelength excitations.

Therefore, the spin current in Eq.  (46) is generally unstable and not a superflow. This instability originates from the interaction process described by in the second quantized Hamiltonian. This energetically favored process converts two particles in the m = 1, − 1 components, respectively, into two stationary particles in the m = 0 component. To suppress such a process and achieve a stable pure spin current, one can utilize the quadratic Zeeman effect. With the quadratic Zeeman effect of a negative coefficient, the Hamiltonian adopts an additional term λ m²(λ < 0 and m = 1, 0, − 1). This term does not change the energy of the m = 0 component, but lowers the energy of the other two components m = 1, − 1. As a result, there arises a barrier for two atoms in the m = 1, − 1 components scattering to the m = 0 component, and the scattering process is thus suppressed.

The above intuitive argument can be made more rigorous and quantitative. Consider the case k₁ = − k₂. With the quadratic Zeeman term, the excitation spectrum for the m = 0 component changes to

Therefore, as long as − λ − ħ ²| k₁| ²/2M > 0, long-wavelength excitations will be stable for the m = 0 component. From the excitation of the m = 0 component, one can obtain a critical relative velocity of the spin current, . There is another nonzero critical velocity determined by the excitations of the m = 1, − 1 components. The overall critical velocity of the system is the smaller one of v₀ and v₁. Therefore, below the critical relative velocity v_c = min{v₀, v₁}, the pure spin current is stable and a superflow. The experimental scheme to realize such a Zeeman effect is discussed in Subsection  5.3.

In the cylindrical geometry, we consider a pure spin current formed by two vortices with opposite circulation in the m = 1, − 1 components. From similar arguments, one can expect that interaction will make such a current unstable. Inspired by the quadratic Zeeman effect method above, we propose to use SOC to stabilize it. The SOC can be viewed as a momentum-dependent effective magnetic field that is exerted only on the m = 1, − 1 components. Therefore, it is possible that SOC lowers the energy of m = 1, − 1 components and, consequently, suppresses the interaction process leading to the instability.

The model of a spin-1 BEC subjecting to the Rashba SOC can be described by the following energy functional:

where ρ is the density, V(r) = Mω ²(x² + y²)/2 is the trapping potential, and γ is the strength of SOC. 〈 · · · 〉 is the expectation value taken with respect to the three component wave function ψ = (ψ ₁, ψ ₀, ψ _{− 1})^T. Here we assume the confinement in the z direction is very tight with trapping frequency ħ ω _z ≫ μ , k_BT, μ being the chemical potential and T the temperature, and have integrated the model with respect to the z direction to obtain an effective two-dimensional system. The two-dimensional coupling constants are related to the three-dimensional scattering lengths through the relations, , and , with .^[83] The SOC strength γ defines a characteristic length a_soc = ħ /Mγ , and can be rescaled to be dimensionless with respect to the harmonic oscillator length . Then we characterize the strength of SOC with the dimensionless quantity . The SOC of Rashba type here can be generated in various ways, which are discussed in the next subsection.

The above model can describe a spin-1 BEC of ²³Na confined in a pancake harmonic trap. Assume the atom number is about 10⁵. Using the estimate of scattering lengths a₀ = 50a_B, a₂ = 55a_B, ^[84] with a_B being the Bohr radius, the ground state of spin-1 ²³Na should be antiferromagnetic because c̃ ₀ > 0 and c̃ ₂ > 0.^[45] Previous studies of spin-1 BEC with Rashba SOC mostly focus on the strong SOC regime, where the ground state is found to be the plane-wave phase or the stripe phase, for ferromagnetic interaction and antiferromagnetic interaction, respectively.^[44] Here we are interested in the antiferromagnetic interaction case and the weak SOC regime (κ ≪ 1), and calculate the ground-state wave function of the energy functional with the method of imaginary time evolution.

We find that when the SOC is weak (κ ≪ 1), the ground-state wave function has the form

with χ ₁(r) = − χ _{− 1}(r) and all χ _i real. The ground state is shown in Fig.  8. Such a ground state consists of an anti-vortex in the m = 1 component and a vortex in the m = − 1 component. The m = 0 component does not carry angular momentum. Since | ψ ₁| = | ψ _{− 1}| , the net mass current vanishes.

Fig.  8.

Figure Option ViewDownloadNew Window

Fig.  8. Amplitudes (a1, b1, c1) and phase angles (a2, b2, c2) of the three component wave function ψ = (ψ 1, ψ 0, ψ − 1)T at the ground state of Hamiltonian  (49) for a BEC of 23Na confined in a pancake trap. The particle number is 105, the radial frequency of the trap is 2π × 55  Hz, and the dimensionless SOC strength is κ = 0.04. The confinement in the z direction is very tight with ω z = 2π × 4200  Hz, so the system is effectively two-dimensional. Units of the X and Y axes are ah.

The wave function  (50) can be understood at the single-particle level. In terms of the ladder operators of spin and angular momentum, the SOC term reads

where S_± is the ladder operator of spin, and is the creation operator of the left (right) circular quanta.^[85] When the SOC is very weak (κ ≪ 1), its effect can be accounted for in a perturbative way. From the ground state Ψ ⁽⁰⁾ = | 0, 0〉 , the first-order correction to the wave function for small γ is given by

where | m_s, m_o〉 denotes a state with spin quantum number m_s and orbital magnetic quantum number m_o. One immediately sees that ψ ₁ has angular momentum − ħ and ψ _{− 1} has angular momentum ħ . In addition, the amplitudes of both ψ ₁ and ψ _{− 1} are proportional to κ .

There exits a continuity equation for the spin density and spin current, which is

The spin current density tensor denotes the spin component) is defined as^{[86, 87]}

where

and c.c. means the complex conjugate. The second part in v_nl is induced by the SOC.

By the definition in Eq.  (54), the spin current density carried by the ground state  (50) is

From both analytical and numerical results of the wave function, | ψ ₁| ≪ | ψ ₀| , so roughly points in the y direction, while almost points in the − x direction. represents a flow whose amplitude has rotational symmetry. From the numerical results shown in Fig.  9, we see that is a counterclockwise flow. The amplitudes of and are of the same order, both proportional to κ , while that of , proportional to κ ², is much lower. It is evident that the state in Eq.  (50) carries no mass current and only pure spin current. Since the spin current is in the ground state, it must be stable. In this way, we have realized a superfluid of pure spin current, or a pure spin supercurrent.

Fig.  9.

	Figure Option ViewDownloadNew Window
	Fig.  9. Distribution of the spin current densities (a), (b), and (c) of the ground state shown in Fig.  8. The length of the arrows represents the strength of the spin current. The arrow length of different subfigures is not to scale. κ = 0.04. Units of the X and Y axes are ah.

In this subsection, we propose the experimental schemes to generate and detect the pure spin currents discussed in Subsections  5.1 and 5.2.

The planar pure spin current can be easily generated. By applying a magnetic field gradient, the two components m = 1, − 1 will be accelerated in opposite directions and a pure spin current is generated as done in Refs.  [70] and [71]. To stabilize this spin current, one needs to generate the quadratic Zeeman effect. We apply an oscillating magnetic field Bsinω t with the frequency ω being much larger than the characteristic frequency of the condensate, e.g., the chemical potential μ . The time averaging removes the linear Zeeman effect; only the quadratic Zeeman effect remains. The coefficient of the quadratic Zeeman effect from the second-order perturbation theory is given by λ = (gμ _BB)²/Δ E_hf, where g is the Landé g factor of the atom, μ _B is the Bohr magneton, and Δ E_hf is the hyperfine energy splitting.^[88] For the F = 2 manifold of ⁸⁷Rb, Δ E_hf < 0, so the coefficient of the quadratic Zeeman effect is negative.

The circular flow in Subsection  5.2 may find prospective realizations in two different systems: cold atoms and exciton BEC. In cold atoms, we consider a system consisting of a BEC of ²³Na confined in a pancake trap, where the confinement in the z direction is so tight that one can treat the system effectively as two-dimensional. The SOC can be induced by two different methods. One is by the exertion of a strong external electric field E in the z direction. Due to the relativistic effect, the magnetic moment of the atom will experience a weak SOC, where the strength γ = gμ _B| E| /Mc². Here M is the atomic mass and c is the speed of light. For weak SOC (small γ ), the fraction of atoms in the m = 1, − 1 components is proportional to γ ². For an experimentally observable fraction of atoms, e.g., 0.1% of 10⁵ atoms, using the typical parameters of ²³Na BEC, the estimated electric field is of the same order of magnitude as the vacuum breakdown field. For atoms with smaller mass or larger magnetic moment, the required electric field can be lowered. Another method of realizing SOC is to exploit the atom laser interaction, where strong SOC can be created in principle.^[89] In exciton BEC systems, as the effective mass of exciton is much smaller than that of atom, the required electric field is four to five orders of magnitude smaller, which is quite feasible in experiments.^{[90– 93]}

The vortex and anti-vortex in the m = 1, − 1 components can be detected by the method of time of flight. First one can split the three spin components with the Stern– Gerlach effect. The appearance of vortex or anti-vortex in the m = 1, − 1 components is signaled by a ring structure in the time-of-flight image. After a sufficiently long expansion time, the ring structure should be clearly visible.