Spin–orbit-coupled Bose–Einstein-condensed atoms confined in annular potentials

A spin–orbit-coupled Bose–Einstein-condensed cloud of atoms confined in an annular trapping potential shows a variety of phases that we investigate in the present study. Starting with the non-interacting problem, the homogeneous phase that is present in an untrapped system is replaced by a sinusoidal density variation in the limit of a very narrow annulus. In the case of an untrapped system there is another phase with a striped-like density distribution, and its counterpart is also found in the limit of a very narrow annulus. As the width of the annulus increases, this picture persists qualitatively. Depending on the relative strength between the inter- and the intra-components, interactions either favor the striped phase, or suppress it, in which case either a homogeneous, or a sinusoidal-like phase appears. Interactions also give rise to novel solutions with a nonzero circulation.


I. INTRODUCTION
The effects associated with spin-orbit coupling have long been known and studied in various physical systems, including atoms, solids, quantum dots, etc. [1][2][3][4][5][6][7], and nuclei [8].Spin-orbit coupling often plays an important role in semiconductor nanostructures, where the field of spintronics has led to many new applications [9].Concepts known from semiconductor nanoscience are nowadays often found to be transferable to confined atomic quantum gases, with an ever increasing interest in the properties of "atomtronic" devices that make use of cold atoms in a similar way like electrons in solids [10][11][12][13][14][15][16][17][18].Semiconductor "spintronic" devices are to a large extent controlled by the Rashba or Dresselhaus spin-orbit coupling in the heterostructure, giving rise to a number of intriguing spin-dependent transport phenomena [19].By means of laser-coupling techniques an analogue to spintronics has now become possible with the achievement of artificially-induced spin-orbit coupling in ultracold quantum gases of neutral atoms [20] with the exciting possibility of creating "atom-spintronic" devices that may also make use of the superfluid properties of quantum gases.Recently, such systems have been experimentally realized in Bose [21,22], as well as in Fermi [23] gases.Their unique properties and their fundamental and applicationrelated prospects have inspired intense theoretical efforts .Furthermore, in a series of other experiments, it has become possible to trap atoms in annular/toroidal traps, see, e.g., Refs.[52][53][54].Motivated by the above experiments, we investigate the lowest-energy states that appear in a spin-orbitcoupled Bose-Einstein-condensed cloud of atoms that is trapped in an annular potential, using the mean-field approximation.A lot of work has been done in homogeneous [26,27,[55][56][57][58] and in harmonically-trapped sys-tems [28,33,36,39,40,44,[59][60][61].In the presence of an annular potential, however, this problem becomes surprisingly challenging.As we show below, even in the absence of interactions the eigenvalue problem due to the spin-orbit coupling is non-trivial.
We proceed in three steps: First of all, we ignore the interactions and consider the case of a very tight annulus, developing an effectively one-dimensional theory.Within this model we integrate over the transverse degrees of freedom, which not only simplifies the problem numerically, but it also gives some additional insight into the properties of the system.Furthermore (Sec.III), the resulting equations have some limiting analytic solutions.We then solve the eigenvalue problem in an annulus of a finite width (Sec.IV).The interactions are incorporated in the last step (Sec.V).We pay special attention to the interplay between the (single-particle) effect of spin-orbit coupling and the atom-atom interactions; we stress that the parameters that are associated with both of them are realistic and controllable experimentally in atomic systems.

II. MODEL AND METHOD
We consider a Bose-Einstein-condensed cloud of atoms, confined to two dimensions.Thus, the single-particle Hamiltonian of the system is The external confining potential has the form of an annular potential, , where ρ is the radial coordinate, M is the atom mass, ω is the frequency of the potential, and R is the mean radius of the annu-lus.The term V SO associated with the spin-orbit coupling that we consider here is the one that corresponds to the experiments of the Spielman group [22].In this case the so-called Rashba-like and Dresselhaus-like spinorbit contributions have an equal strength, and thus V SO takes the one-dimensional form Here k 0 is the wavenumber of the Raman laser beams, p x is the linear momentum of the atoms in the x direction, σ j (with j = x, y, z) are the Pauli matrices, Ω is the Rabi frequency and δ is the detuning.The experimental studies of spin-orbit coupling generally choose small values of δ [21,22] and thus we ignore the corresponding term in what follows below.Furthermore, we here orient the spin axes as in Ref. [22], but most recent theoretical work [24,56] uses a cyclic spin rotation σ y → σ z , etc. Therefore, the system is described as an effective pseudospin-1/2, and the condensate wave function is written as a two-component vector [Φ ↑ , Φ ↓ ] T ("up" and "down").
The interactions are modelled via a short-ranged swave potential of the form While in general there are three different interaction strengths between the "up" and "down" components, g ↑↑ , g ↓↓ , and g ↑↓ , in the present problem we set g ↑↑ = g ↓↓ = g and tune the ratio g/g ↑↓ (both assumed to be positive).Our total Hamiltonian is thus H = H 0 + V int , which yields two coupled Gross-Pitaevskii-like equations for the two components of the order parameter, where and µ is the chemical potential, which is determined from the condition that the total occupancy of the two pseudospin components N ↑ and N ↓ is equal to the total number of atoms N .
For the specific form of the spin-orbit coupling considered and in the absence of any trapping potential and of interparticle interactions, the density distribution of the gas may be either striped, or homogeneous [27,34,[62][63][64].This is determined by the ratio between hΩ and the recoil energy E R = h2 k 2 0 /(2M ) [22].If the ratio hΩ/(4E R ) is smaller or larger than unity, the system is in the striped, or in the homogeneous phase, respectively.Furthermore, the expectation value of The term proportional to k 2 0 is constant.The term proportional to Ω is spatially-independent and is responsible for the population imbalance between the two components.The most "interesting" is the spatially-dependent term, which is proportional to k 0 .It is this term that gives rise to the striped phase.

III. NON-INTERACTING PROBLEM IN THE LIMIT OF QUASI-ONE-DIMENSIONAL MOTION
For very strong transverse confinement, the atoms reside in the ground state associated with the transverse degrees of freedom for motion along the annulus.This fact allows us to make an ansatz for Φ ↑/↓ , with a Gaussian density distribution in the transverse direction, i.e., n 0 (ρ , where a 0 is the oscillator length, a 0 = (h/M ω) 1/2 .Integrating over the transverse degrees of freedom and assuming that R ≫ a 0 we thus develop the following quasi-one-dimensional eigenvalue problem for Ψ ↑ and Ψ ↓ (i.e., two coupled differential equations), where E is the eigenenergy.
In the absence of interactions there are three energy scales, h2 /(M R 2 ), h2 k 0 /(M R), and hΩ.The above coupled equations have an analytic solution in the limit where h2 k 0 /(M R) is much smaller either than h2 /(M R 2 ), or than hΩ (we stress that for k 0 = 0 the two equations decouple, while for "small" values of k 0 R they couple perturbatively).Defining the two dimensionless parameters k 0 R ≪ 1 (i.e., the ratio between the second scale and the first), and λ ≡ hk 0 /(M ΩR) (i.e., the ratio between the second scale and the third), we find that for The eigenenergy is where ǫ ≡ 2M R 2 Ω/h is defined as the ratio between hΩ and h2 while the eigenenergy is E/(hΩ) = −1/2 − λ 2 /8.Expanding the eigenfunctions in the basis of the planewave states, we have also solved the eigenvalue problem numerically.The eigenfunctions of lowest energy have the general form Ψ ↑ = m c 2m cos 2mθ and Ψ ↓ = m d 2m+1 cos(2m + 1)θ, or vice versa, i.e., with Ψ ↑ exchanged with Ψ ↓ .From these expressions some generic features follow for Ψ ↑ and Ψ ↓ , which are also valid for the solutions found in Eqs. ( 7) and ( 8): (i) These solutions are even (the Hamiltonian is invariant under the transformation θ → −θ and thus the eigenstates are parity eigenstates).(ii) The density n ↓ has nodes at θ = π/2 and 3π/2, while at the same points n ↑ has maxima, or vice versa.(iii) The periodicity of Ψ ↑ is π and of Ψ ↓ it is 2π, or vice versa.(iv) Finally, the density of both components has a periodicity of π.
Turning to the degeneracy of these solutions, the lowest-energy eigenstates are even.For all current experiments, the laser beam has a wavelength of order 1 µm, which means that k 0 ∼ 2π (µm) −1 .For a typical value of R ∼ 10 µm, k 0 R is thus at least 10, or larger.In the typical limiting case k 0 R ≫ 1, the system is in the striped-like phase, the density variations are located around opposite ends of the ring, while the density effectively vanishes elsewhere [see, e.g., Fig. 1(e)].There is thus a corresponding antisymmetric solution which becomes nearly degenerate with the symmetric one in this limit (with a corresponding two-fold degeneracy).In addition, the transformation Ψ ↑ → Ψ ↓ , Ψ ↓ → −Ψ ↑ and Ω → −Ω leaves the eigenvalue problem unaffected.This gives rise to another two-fold (near) degeneracy for "small" Ω, which is exact for Ω = 0.
The physical picture that emerges from the above calculation is that for k 0 R ≪ 1 (and any value of Ω), i.e., for sufficiently small values of k 0 and/or of R, as well as for hk 0 /(M ΩR) ≪ 1, i.e., for sufficiently small k 0 and/or large values of Ω and/or large values of R, the density variation of both components is sinusoidal.As k 0 R increases (with all the other parameters fixed), the two components start to localize around θ = π/2 and θ = 3π/2, with a pronounced density variation around these two points, which is the analogue of the striped phase.This transition is seen in Fig. 1, where we plot the one-dimensional density of the two components for an increasing value of k 0 , k 0 R = 1.0 (a), 1.5 (b) and 1.6 (c), 2.8 (d), and 6 (e), and for a fixed value of h/(M ΩR 2 ) = 1.The population imbalance also decreases with increasing k 0 , since the term in the Hamiltonian that is proportional to Ω becomes smaller compared to the one that is proportional to k 0 .It is also remarkable that an increase of the value of k 0 R of less than a factor of two results in a rather dramatic change in the density distribution of the two components, as seen, e.g., between Figs. 1 (b) and 1 (c).As we show also below, this indicates a very rich structure, which is a generic feature of this problem.

IV. EIGENVALUE PROBLEM IN THE CASE OF AN ANNULAR POTENTIAL
We now investigate the effect of the finiteness of the width of the annulus, still in the absence of interactions.In this case, the width of the annulus (the oscillator length a 0 ) provides an extra length scale, which is a fraction of R. For a narrow annulus (R ≫ a 0 ) one goes back to the problem of quasi-one-dimensional motion discussed above.When a 0 ≈ R then the problem reduces to that of a harmonic trapping potential [28,33,36,39,40,44,59-61] with a small "hole" in the center of the trap.The most interesting case is thus the intermediate one when a 0 < ∼ R. We show in panel (b) of Fig. 2 the two-dimensional density distribution of the two components that we have obtained numerically using the method of imaginary-time propagation (for the details of this calculation see the Appendix), considering an annular potential with R/a 0 = 4 and some representative values of k 0 and Ω.In the first column, of "small" k 0 , k 0 R = 0.4, there is a large population imbalance between the two components -here hΩ/E R = 200.While the "up" component has a pronounced axial asymmetry, the dominant "down" component also lacks axial symmetry with the maximum of the density being at θ = 0 and π (which is very weak, and is hardly visible in the plot).In the second column with a "large" value of k 0 , k 0 R = 10, both components show a striped phase and have an almost equal population; here hΩ/E R = 8/25.Finally the phase shown in the third column with Ω/ω = 10 and k 0 R = 10 resembles in a sense the phase with a sinusoidal density distribution found within the quasi-one-dimensional model, with a large population imbalance (here hΩ/E R = 16/5).

V. EFFECT OF THE INTERACTIONS
Having established a rather complete picture about the single-particle eigenvalue problem, we now turn to FIG. 2: Two-dimensional density of the two pseudospin components, "up" (higher) and "down" (lower).Panel (a) shows the result of the quasi-one-dimensional model with a Gaussian transverse profile.Panel (b) shows the density in an annular potential and in the absence of interactions.Panel (c) shows the result in an annular potential in the presence of interactions for g < g ↑↓ and panel (d) for g > g ↑↓ .The color scale is not the same in all the plots, since in some of them there is a large population imbalance.
the effect of the interactions.In a homogeneous system when g ↑↓ > g the striped phase (where also the density minima/maxima of the one component coincide with the maxima/minima of the other) is energetically more favourable, since this configuration minimizes the overlap between the two components [due to the last term in Eq. ( 9) below].This may be seen if one writes the FIG.3: Two-dimensional density (black contour lines) and phase of the two pseudospin components, "up" (left) and "down" (right) for the data of the third column of Fig. 2 (d).The color scale shows the phase with blue as zero and red as 2π.
interaction energy in the following form (9) where n(r) = n ↑ (r)+n ↓ (r) is the total density.According to the numerical results that we have obtained using the method of imaginary-time propagation, when g ↑↓ > g the interactions "combine" with the (spatially dependent) term of the spin-orbit Hamiltonian that is proportional to k 0 (and favours the formation of stripes, as we argued earlier).Both terms make it energetically favourable for the system to form a striped-like phase (similar to the one that we described in the case of an annular potential with spin-orbit coupling, but without interactions).Panel (c) of Fig. 2 shows n ↑ and n ↓ , for g/g ↑↓ = 1/2, and N g/(2πRa 0 hω) = 25/(4π), or N g/(2πRa 0 E R ) = 1250/π in the first column and N g/(2πRa 0 E R ) = 2/π in the second and in the third column.While in the first two columns the density is slightly affected by the interactions [i.e., as compared to the density in Fig. 2 (b)], in the third column the interactions expand the two components more around the annulus.
In panel (d) of Fig. 2 we have chosen g/g ↑↓ = 2.The first column indicates that the density is not affected drastically by the interactions.On the other hand, in the second column we see that the density is essentially axially symmetric, in contrast to the striped-like phase, which is seen in the absence of interactions in the second column of panel (b).Thus, this phase is strongly affected by the interactions.
At this point it is interesting to make contact between the results of Fig. 2 and those of no trapping potential.As we mentioned also above, in the absence of interactions when the ratio hΩ/(4E R ) is smaller or larger than unity, the system is in the striped, or in the homogeneous phase, respectively [22].The actual value of this ratio is 50 (first column), 0.08 (second column), and 0.8 (third column).According to the above criterion, the first and the second columns of Fig. 2 (b) are in the homogeneous and in the striped phases, respectively, which is consistent with our data.For the interacting problem, the critical value for the ratio hΩ/(4E R ) is 0.19/4 = 0.0475 in this case [22,34,[62][63][64].Therefore, the first and the third columns of Fig. 2(c) and (d) are clearly in the homogeneous phase, which again is consistent with our data.
In the third column of Fig. 2 (d) we show a solution of nonzero circulation, which is only possible in the presence of interactions and for a wide enough annulus.For the same data, Fig. 3 shows the phase of the two order parameters along with the corresponding density.As seen in this plot, there are two phase singularities in the "down" component and three in the "up".In the "down" component the corresponding vortices reside in the central region of exponentially small density, while in the "up" component one of the three is in the central region and two in the region of nonzero density (for this reason the deviation of the density of the "up" component from axial symmetry is more pronounced than that of the "down" component).For the same reason, while all the other plots of Fig. 2 have a mirror symmetry with respect to the x axis, this symmetry is broken.We should also mention that the sense of circulation of the state of the third column of Fig. 2 (d) is determined by the initial condition that is given in the (iterative) method of imaginary-time propagation.Furthermore, with a "mean" density n 0 = N/(2πRa 0 ), for the healing length ξ 0 that corresponds to the coupling constant g, ξ 0 /a 0 = hω/(2n 0 g) = 2π/25 ≈ 1/2.Therefore, ξ 0 is roughly one half of the width of the annulus, since a 0 sets this scale.Clearly these length scales play a crucial role, since the vortex states we have found have to "fit" within the annulus and this becomes possible provided that the width of the annulus is wide enough, i.e., it is comparable, or larger than the coherence length ξ 0 .
The vortex states we have found in Fig. 3 involve two large parameters: Rabi coupling Ω/ω = 10 and Raman wavenumber k 0 R = 10, as well and the important effect of interactions with g > g ↑↓ .Each of these plays an important role, since the vortex states appear only in the right column of Fig. 2(d).As seen in Fig. 3, there are persistent currents in both the upper component (arising from the phase singularity in the central region) and the lower component (arising from the two phase singularities in the central region).In addition, the upper component has two vortices in the annular region with finite number density.Such states [59,60,65,66] have the potential to serve as cold-atom analogs of ring currents and definitely merit further more detailed investigation.

VI. CONCLUSIONS
Spin-orbit coupled Bose-Einstein-condensed atoms confined in an annular potential show a variety of phases as a result of the (single-particle effect of the) spin-orbit coupling and of the (many-body effect of the) coupling between the two atoms, while the non-trivial topology of the annular potential is another novel aspect of the problem that we have solved.
Figure 1 summarizes the results for the solutions of the eigenvalue problem in the case of a very narrow an-nulus, as k 0 increases (from top to bottom).The homogeneous phase that is present in an untrapped system is replaced by a sinusoidal density variation in the case of a very narrow annulus, which evolves continuously into a striped-like phase as k 0 increases.As the width of the annulus increases, this picture persists qualitatively, as seen in Fig. 2, where interactions have also been included.Depending on the relative strength of g and g ↑↓ the interactions either favor the striped phase, or suppress it, while they may also give rise to nonlinear solutions with a nonzero circulation that require nonzero interactions.The very interesting implication of these solutions is that in a spin-orbit coupled system vortex states might spontaneously be created in one or both components, sometimes in the physical region with nonzero density and sometimes in the central region with zero density (these latter vortices lead to quantized circulation around the annulus).
While our study does not in any way exhaust all the possible phases, it gives a sense of the richness of this problem.A natural future project is a more systematic study of the solutions presented in Fig. 2 (d) and Fig. 3. Finally, the dimensional reduction we have performed provides a general method, which (with the trivial inclusion of interactions) may be useful in investigating questions related with e.g., nonlinear, solitary-wave solutions in quasi-one-dimensional spin-orbit coupled systems.
The numerical method used to solve the two coupled Gross-Pitaevskii-like equations for the two pseudospin components of the order parameter, given by Eq. ( 4) of the main text, is based on a fourth-order split-step Fourier method within an imaginary-time propagation technique.The details of the method were previously discussed in Ref. [67].Therefore here we give only a brief overview of the method.
The starting point, which forms the basis of the imaginary-time propagation method, is to consider the time-dependent Schrödinger equation in imaginary time, i.e., τ = −it (in the expressions below we use atomic units for convenience) where H is the Hamiltonian, and to propagate the wave function ψ under the action of a time evolution operator exp(−τ H), If the initial state ψ(0) is expanded in the eigenstates φ i of the Hamiltonian, with eigenvalues ε i , then Therefore, the time-evolution operator leads to an exponential decay of the eigenstates, where the one with the lowest energy ε 0 decays with the slowest rate and ψ(τ ) eventually converges to the ground state φ 0 as τ → ∞.
To evaluate the time evolution given in Eq. (A.3), we use a forward fourth-order algorithm called split-step Fourier method: ψ(∆τ ) = e −(1/6)∆τ V (∆τ ) e −(1/2)∆τ T e −(2/3)∆τ Ṽ (∆τ /2) × × e −(1/2)∆τ T e −(1/6)∆τ V (0) ψ(0). (A.4) In our study V = H − K, where K is the kinetic energy, serves as an effective potential and it consists of the external trapping potential, the interaction and the spinorbit coupling terms.The mid-point effective potential in Eq. (A.4) is given by The order parameter is discretized on a square mesh of 2 N points (N is an integer) in both x and y directions with separation ∆l and its Fourier transform is computed using the discretized fast Fourier transform (FFT).
In our calculations we have used 256 2 grid points with ∆l = 0.1, where we have observed that larger grid sizes do not have a significant effect on the obtained results.The time step ∆τ had to be chosen small enough to ensure that the Hamiltonian of the system does not change dramatically between the iterations.The optimum value of ∆τ , which satisfies this condition and also provides a good convergence in our study, has been determined as 0.01.On the other hand, we have observed that the use of smaller time steps does not improve the accuracy of our results.
The time evolution procedure requires an initial condition for the wave function that is propagated in imaginary time until a steady-state solution with the (local or absolute) minimum of energy is reached after a sufficiently large number of iterations.For example, to obtain the solution with nonzero circulation given in the third column of Fig. 2 (d) and in Fig. 3, we have performed calculations using ten different initial states.We have found that all the results exhibit density distributions with a nonzero circulation and the result presented in Fig. 2 (d) and in Fig. 3 is the energetically most favorable state.