The in-plane gradient magnetic field induced vortex lattices in spin-orbit coupled Bose-Einstein condensations
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] N ov The in-plane gradient magnetic field induced vortex lattices in spin-orbit coupledBose-Einstein condensations
Xiang-Fa Zhou,
1, 2
Zheng-Wei Zhou,
1, 2
Congjun Wu, and Guang-Can Guo
1, 2 Key Laboratory of Quantum Information, University of Science and Technology of China,Chinese Academy of Sciences, Hefei, Anhui 230026, China Synergetic Innovation Center of Quantum Information and Quantum Physics,University of Science and Technology of China, Hefei, Anhui 230026, China Department of Physics, University of California, San Diego, CA 92093, USA
We consider the ground-state properties of the two-component spin-orbit coupled ultracold bosonssubject to a rotationally symmetric in-plane gradient magnetic field. In the non-interacting case,the ground state supports giant-vortices carrying large angular momenta without rotating the trap.The vorticity is highly tunable by varying the amplitudes and orientations of the magnetic field.Interactions drive the system from a giant-vortex state to various configurations of vortex latticestates along a ring. Vortices exhibit ellipse-shaped envelops with the major and minor axes deter-mined by the spin-orbit coupling and healing lengths, respectively. Phase diagrams of vortex latticeconfigurations are constructed and their stabilities are analyzed.
PACS numbers: 03.75.Mn, 03.75.Lm, 03.75.Nt, 67.85.Fg
I. INTRODUCTION
Spin-orbit (SO) coupling plays an important role incontemporary condensed matter physics, which is linkedwith many important effects ranging from atomic struc-tures, spintronics, to topological insulators [1–3]. It alsoprovides a new opportunity to search for novel states withultracold atom gases which cannot be easily realized incondensed matter systems. In usual bosonic systems,the ground state condensate wavefunctions are positive-definite known as the “no-node” theorem [4, 5]. However,the appearance of SO coupling invalidates this theorem[6]. The ground state configurations of SO coupled Bose-Einstein condensations (BEC) have been extensively in-vestigated and a rich structure exotic phases are obtainedincluding the ferromagnetic and spin spiral condensations[6–9], spin textures of the skyrmion type [6, 10–13], andquantum quasi-crystals [14], etc. On the experiment side,since the pioneering work in the NIST group [15], it hasreceived a great deal of attention, and various furtherprogresses have been achieved [16–20]. Searching fornovel quantum phases in this highly tunable system isstill an on-going work both theoretically and practically[21–28], which has been reviewed in [29–33].On the other hand, effective gradient magnetic fieldshave been studied in various neutral atomic systems re-cently. For instance, it has been shown in Ref. [34, 35]that SO coupling can be simulated by applying a se-quence of gradient magnetic field pulses without involv-ing complex atom-laser coupling. In optical lattices, the-oretic and experimental progresses show that SO cou-pling and spin Hall physics can be implemented with-out spin-flip process by employing gradient magnetic field[36, 37]. This represents the cornerstone of exploring richmany-body physics using neutral ultracold atoms. Ad-ditionally, introducing gradient magnetic fields has alsobeen employed to create various topological defects in- cluding Dirac monopoles [38] and knot solitons [39]. Itwould be very attractive to investigate the exotic physicsby combining both SO coupling and the gradient mag-netic field together in ultracold quantum gases.In this work, we consider the SO coupled BECs sub-ject to an in-plane gradient magnetic field in a 2D ge-ometry. Our calculation shows that this system supporta variety of interesting phases. The main features aresummarized as follows. First, the single-particle groundstates exhibit giant vortex states carrying large angularmomenta. It is very different from the usual fast-rotatingBEC system, in which the giant vortex state appears onlyas meta-stable states [40, 41]. Second, increasing the in-teraction strength causes the phase transition into thevortex lattice state along a ring plus a giant core. Thecorresponding distribution in momentum space changesfrom a symmetric structure at small interaction strengthsto an asymmetric one as the interaction becomes strong.Finally, the size of a single vortex is determined by twodifferent length scales, namely, the SO coupling strengthtogether with the healing length. Therefore, the vor-tex exhibits an ellipse-shaped envelope with the princi-ple axes determined by these two scales. This is differentfrom the usual vortex in rotating BECs [42–47], wherean axial symmetric density profile is always favored.The rest of this article is organized as follows. In Sect.II, the model Hamiltonian is introduced. The single par-ticle wavefunctions are described in Sect. III. The phasetransitions among different vortex lattice configurationsare investigated in Sect. IV. The possible experimentalrealizations are discussed in Sect. V. Conclusions arepresented in Sect. VI.
II. THE MODEL HAMILTONIAN
We consider a quasi-2D SO coupled BEC subject toa spatially dependent magnetic field with the followingHamiltonian as H = Z d~r ˆ ψ ( ~r ) † n ~p m + Λ r (cos θ ˆ r + sin θ ˆ ϕ ) · ~σ + 12 mω r o ˆ ψ ( ~r ) + H soc + H int , (1)where ˆ r = ~r/r with ~r = ( x, y ); ~σ = ( σ x , σ y ) are the usualPauli matrices; m is the atom mass; ω is the trappingfrequency; Λ is the strength of the magnetic field, and θ denotes the relative angle between the magnetic fieldand the radial direction ˆ r . Physically, this quasi-2D sys-tem can be implemented by imposing a highly anisotropicharmonic trap potential V H = m ( ω r + ω z z ). When ω z ≫ ω , atoms are mostly confined in the xy -plane,and the wavefunction along z axis is determined as aharmonic ground state with the characteristic length a z = p ~ / ( mω z ).For simplicity, the SO coupling employed below hasthe following symmetric form as H soc = Z d~r ˆ ψ ( ~r ) † h λm ( p x σ x + p y σ y ) i ˆ ψ ( ~r )with λ the SO coupling strength. We note that due tothis term, the magnetic fields which couples to spin canbe employed as a useful method to control the orbit de-gree of freedom of the cloud. The interaction energy iswritten as H int = g D Z d~r ˆ ψ ( ~r ) † ˆ ψ ( ~r ) † ˆ ψ ( ~r ) ˆ ψ ( ~r ) . (2)Here the contact interaction between atoms in bulk is g = 4 π ~ a s /m , where a s is the scattering length. Forthe quasi-2D geometry that we focus on, the effectiveinteraction strength is modified as g D = g D / ( √ πa z ). III. SINGLE-PARTICLE PROPERTIES
The physics of Eq. 1 can be illustrated by consider-ing the single-particle properties first. After introduc-ing the characteristic length scale of the confining trap l T = p ~ /mω , the dimensionless Hamiltonian is rewrit-ten as H ~ ω = Z d~ρ ˆ φ ( ~ρ ) † n − ~ ∇ βρ (cos θ ˆ r + sin θ ˆ ϕ ) · ~σ + α~k · ~σ + 12 ρ o ˆ φ ( ~ρ ) , (3)where α = λ/ ( mωl T ) and β = Λ l T / ( ~ ω ) are the dimen-sionless SOC and magnetic field strengths, respectively;the normalized condensates wave-function is defined as φ ( ~ρ ) = l T √ N Ψ( ~r = l T ~ρ )with N the total number of atoms; E / ( ¯ h ω ) α = 10, θ = 0 (a) β E / ( ¯ h ω ) α = 10, θ = π/ (b) β FIG. 1: The single-particle dispersion of the Hamiltonian Eq.(3) with lower energy branch as a function of the reducedmagnetic fields β for fixed α = 10 and different values of θ = 0(a), and π (b). The inset in (b) shows that the ground statescrossing for certain values of β at θ = π/ = 0, while there isno crossing in (a) at θ = 0. Since the total angular momentum ~ j z = ~ l z + ~ σ z is conserved for this typical Hamiltonian, we can use itto label the single-particle states. If the magnetic fieldalong the radial direction, i.e., θ = lπ , the Hamiltonianalso supports a generalized parity symmetry describedby iσ y P x , namely [ H , iσ y P x ] = 0 , (4)with P x the reflection operation about the y -axis satisfy-ing P x : ( x, y ) → ( − x, y ). Therefore for given eigenstates φ m = [ f ( ρ ) e imϕ , g ( ρ ) e i ( m +1) ϕ ] T with j z = ( m + 1 / { φ m , ( iσ y P x ) φ m } are degenerate for H ( θ = lπ ). Thissymmetry is broken when θ = lπ .Due to the coupling between the real space magneticfield and momentum space SO coupling, the single par-ticle ground states exhibit interesting properties at largevalues of α and β . In momentum space, the low energystate moves to a circle with the radius determined by α .The momentum space single-particle eigenstates breakinto two bands ψ ± ( ~k ) with the corresponding eigenvalues E ± ~k / ( ~ ω ) = ( | ~k | ± α | ~k | ) and eigenstates √ [1 , ± e iθ k ] T , FIG. 2: The density and phase profiles of the single-particleground states for fixed α = 6, β = 1, and different θ = π (a), π (b). From left to right: the density and phase profilesfor spin-up and spin-down components, respectively. respectively. For the lower band which we focus on, thespin orientation is h ~σ i = ( − cos θ ~k , − sin θ ~k ), which isanti-parallel to ~k . On the other hand, in the real space,for a large value of β , the potential energy in real space isminimized around the circle with the the radius r/l T = β with a spatial dependent spin polarization. Thereforearound this space circle, the local wavevector at a posi-tion ~r is aligned along the direction of the local magneticfield to minimize the energy. The projection of the localwavevector along the tangent direction of the ring givesrise to the circulation, and thus the ground state carrieslarge angular momentum m which is estimated as m ≃ πβ sin θ/ (2 π/α ) = αβ sin θ. (5)Therefore, by varying the angle θ , a series of groundstates are obtained with their angular momentum rang-ing from 0 to αβ ≫
1. This is very different from theusual method to generate giant vortex, where fast rotat-ing the trap is needed [42].For β ≫
1, the low energy wavefunctions mainly dis-tribute around the circle ρ = β . As shown in AppendixA, the approximated wavefunctions for the lowest band( n = 1) is written as φ n =1 ,j z ( ρ, ϕ ) ≃ π ρ e − ( ρ − β )22 e iρα cos θ × (cid:20) e i [ mϕ − θ ] − e i [( m +1) ϕ + θ ] (cid:21) , (6)where ϕ is the azimuthal angle. The corresponding en-ergy dispersion is approximated as E n,j z ≈ n + 1 − α − β j z − αβ sin θ ) α cos θ + β ) . (7)For given values of α and β , E n,j z is minimized at j z ≃ αβ sin θ , which is consistent with the above discussion.In the case of θ = lπ , two states with m = l and − ( l +1) are degenerated due to the symmetry defined in Eq.4. Interestingly, Eq. 7 also indicates that for integer FIG. 3: The profiles of the condensate wavefunctions of thespin-up component for α = 11, β = 6, and θ = π . The inter-action parameters are g = 15 ( a ), 35 ( b ), 75 ( c ), and 100 ( d ),respectively. We note that (c) and (d) exhibit similar profilesbut with different q . From top to bottom: the density andphase profiles in real space, and the momentum distributionswhich mainly are located around the circle | k | = α . αβ sin θ = l , an approximate degeneracy occurs for m = l and l − β for different values of θ . For θ = 0, the dispersionwith different j z never cross each other Fig. 1 (1a). Thevalues of j z for the ground state are always j z = or − due to the symmetry Eq. 4. When θ = π/ = 0,the spectra cross at certain parameter values, and theground-state can be degenerate even without additionalsymmetries as shown in Fig. 1 (1b), which is consis-tent with above discussions. For β ≫
1, the proba-bility density of the ground state single particle wave-function mainly distributes around a ring with ρ = β .Interestingly, the phase distribution exhibits the typicalArchimedean spirals with the equal-phase line satisfying ρ ∼ mϕ (or ρ ∼ ( m + 1) ϕ ) (see Fig. 2 for details). IV. PHASE TRANSITIONS INDUCED BYINTERACTION
In this section, we consider the interaction effect whichwill couple single-particle eigenstates with different val-ues of j z . It is interesting to consider the possible vortexconfigurations in various parameter regimes, which hasbeen widely considered in the case of the fast rotatingBECs.If the dimensionless interaction parameter g = g D N/ ( ~ ωl T ) is small, it is expected that the groundstate still remains in a giant-vortex state, which is similarto the non-interacting case. The envelope of the varia- FIG. 4: Ground state profiles of the condensates for α = 11, β = 6, θ = π/ g = 85(a), 105(b),125(c), and 145(d) respectively. From top to bottom: den-sity and phase profiles of the spin-up component, momentumdistribution in the lower band along the circle | k | = α . Theorientation of the ellipse-shaped vortices is determined by θ .See text for details. tional wave-function is approximated as φ j z ( ρ, ϕ ) ∼ π √ σρ e − ( ρ − β )22 σ e iρα cos θ (cid:20) e i [ mϕ − θ ] − e i [( m +1) ϕ + θ ] (cid:21) with σ the radial width of the condensates. Around athin ring inside the cloud with the radius ρ , in orderto maintain the overall phase factor e imϕ , the magni-tude of the local momentum along the azimuth directionis determined by k ϕ = m/ρ . Depending on the width σ of the cloud, the linewidth of k ϕ is proportional to δk ϕ = mσ/β . In momentum space, this leads to the ex-pansion of the distribution around the ring with | k | = α .The increasing of the kinetic energy mainly comes fromthe term ˆ E ϕ = ( j z /ρ − α sin θ ) /
2, which is estimatedas h ˆ E ϕ i j z . Details derivation of various energy contribu-tions can be found in Appendix B.Increasing the interaction strength g expands the cloudand leads to larger width σ and δk ϕ , which makes theabove variational state energetically unfavorable. In or-der to minimize the total energy, the condensates tend toinvolve additional vortices such that the local momentummainly distributes around the circle | k | = α with smaller δk ϕ . Fig. (3) and (4) show the typical ground-stateconfigurations for selected parameters. The phase accu-mulations around the inwards and outwards boundariesof the cloud are 2 πm + and 2 πm − respectively. Therefore,there are q = m + − m − vortices involved and distributedsymmetrically inside the condensates. Between two near-est vortices, the local wavefunction can be approximatelydetermined as a plane-wave state. Therefore, their cor-responding distribution in momentum space is also com-posed of q peaks located symmetrically around the circle | k | = α . FIG. 5: Enlarged density and phase profiles around singlevortex. The two figures (a) and (b) are the correspondingparts adapted from Fig. (3c) and (4d) respectively.
As further increase of the interaction strength, the con-densates break into more pieces by involving more vor-tices. The number of the vortices is qualitatively deter-mined by the competition of the azimuthal kinetic en-ergy and the kinetic energy introduced by the vortices.Specifically, if q vortices locate in the middle of the cloudaround the circle ρ ≃ β , then for the inwards part of thecondensates with ρ < ρ , the mean value of the angularmomentum can be approximated as j z, − ≈ j z − q/
2, whilefor the regime with ρ > ρ , we have j z, + ≈ j z + q/
2. Thecorresponding kinetic energy along the azimuthal direc-tion is modified as h ˆ E ϕ i = h ˆ E ϕ i j z + q β ( q − σα sin θ √ π ) . (8)This indicates that, to make the vortex-lattice state fa-vorable, we must have ( q − σα sin θ √ π ) <
0. In the limit casewith θ = 0, this condition is always violated. Therefore,the ground state remains to be an eigenstate of j z with j z = ± even for large interaction strength.We note the vortices display an ellipse-like shape withtwo main axis, as shown in Fig. 5. The phase profileis twisted, and the constant phase front exhibits a dislo-cation around vortex cores. Along the direction of localwavevector ~k , the vortex density profile is determined bythe length scale 2 πl T /α . While perpendicular to the di-rection of local ~k , the vortex profile is dominated by thehealing length ξ due to interaction. Therefore, the vor-tex density distribution is determined by two differentlength scales in mutually orthogonal directions, whichresults in ellipse-shaped vortices. Changing the interac-tion strength and SO coupling alerts the ratio of the twolength scales, thus changes the eccentricity of the ellipses.Additionally, changing the angle θ also changes the di-rection of local magnetic fields, and thus modifies theorientation of the vortices, as shown in Fig. 3 and Fig.4. On the other hand, the introduction of vortices leadto the increase of kinetic energy due to the modifica-tion of the density profile. This can be estimated as Β = g = Π Π Π Π Α Θ c (a) Β= Α= Β= Α=
100 100 200 300 400 500 600 Π Π Π g Θ c (b) FIG. 6: (a) Critical angle θ c as a function of SO couplingstrength α for fixed values of β = 4 and g = 800. (b) θ c decreases as the increase of interaction parameter g for fixed α = 10, 6, and β = 4. . q/ ( βσαξ ), where ξ = 1 / √ gn is the dimensionlesshealing length with n the bulk density of the clouds. Thetotal energy changing due to the presence of the vorticescan be written as∆ E = q β ( q − σα sin θ √ π ) + 0 . qβσαξ . (9)Several interesting features can be extracted from Eq.9. For fixed parameters g , α , and β , there always existsa critical θ c such that ∆ E = 0 is satisfied. When θ < θ c ,then ∆ E >
0, which indicates that a giant-vortex groundstate is always favored. As increasing α , θ c satisfying∆ E = 0 becomes smaller. At θ > θ c , the ground stateexhibits a lattice-type structure along the ring with agiant vortex core. The values of q is determined by mini-mizing ∆ E with respect to g , α , and β , respectively. Fig.(6) shows θ c as a function of SO strength α at which thetransition from a giant-vortex state to a vortex-latticestate occurs. When α is small, a giant-vortex state is fa-vored for all values of θ . As α increases, θ c drops quicklyinitially and decreases much slower when α becomes largeas shown in Fig. (6) (a). In Fig. (6) (b), it shows that asincreasing the interaction g , it is becomes easier to drivethe system into the vortex lattice state.Fig. (7) shows the phase diagram in the α − g planefor a fixed β = 6 for different values of θ . For a fixed α and at small values of g , the system remains to bea giant-vortex state until g reaches its critical value g c .
70 80 90 100 110 120 130 1408910111213 g Α GV IM (a)
20 40 60 80 100 120 140789101112 g Α IMGV (b)
FIG. 7: Phase diagram in the α − g plane for β = 6 with differ-ent θ = π/ π/ M means thatthe condensates support a vortex-lattice-type ground statewith M momentum peaks along the circle | k | = α . The regimewith shadow in (b) indicates the ground state shows multi-layer structure as increasing the interaction strength. Otherphases are defined as follows: GV(giant-vortex state),IM (in-termediate regime). When g > g c , the system enters into an intermediateregime in which vortices start to enter into the conden-sates from boundaries. The momentum distribution alsobreaks into several disconnected segments. More sin-gle quantum vortices are generated in the condensatesas further increasing the interaction strength. The vor-tices distribute symmetrically along the ring and sepa-rate the condensates into pieces. Between two neigh-boring vortices, the condensates are approximated by lo-cal plane-wave states. The momentum distribution com-poses of multi-peaks symmetrically located around thecircle | k | = α . Increasing g also increases the number ofthe single quantum M inside the condensates, hence in-creases the number of peaks in momentum space. For asmaller value θ = π/
3, the critical g c is increased, whichmeans that stronger interactions are needed to drive thesystem into the vortex-lattice states. Interestingly, theintermediate regime is also greatly enlarged. This is con-sistent with the limit case θ = 0, where the system re-mains to be a giant-vortex state even in the case of largeinteraction strength.More ellipse-shaped vortices are formed as furtherincreasing the interaction strength, which are self-organized into a multiple layered ring structure, as shownin Fig.(8). Around each ring, vortices distributed sym- FIG. 8: Density and momentum distributions about theground states of the condensates for θ = π/ g =400((a)), 1000((b)), and θ = π/ g = 400((c)),1000((d)). Other parameters are the same with figure (4).We note that since the two spin components share almost thesame profiles, only the densities of the spin-up component areshown for simplicity. metrically. The number of the vortices between differ-ent layers can be not equal due to their different radius.Therefore the distribution in momentum space becomesasymmetric, and exhibits complex multi-peak structuresaround the circle | k | = α . V. EXPERIMENTAL CONSIDERATION
The Hamiltonian Eq. 1 considered above can be dy-namically generated on behalf of a series of gradient mag-netic pulses [34, 35]. Starting with the typical single-particle Hamiltonian H s = p m + mω r , in the firsttime step, we employ a pair of magnetic pulses U and U † , defined as as U = e iλ ( xσ x + yσ y ) / ~ , (10)at time t = 2 nτ , (2 n + 1) τ respectively. Secondly, atypical effective gradient coupling,Λ[( x cos θ − y sin θ ) σ x + ( y cos θ + x sin θ ) σ y ] , (11)is applied during the whole time duration [(2 n +1) τ, n +1) τ ]. Combining these two time steps, an effective dy-namical evolution U = e − iH τ , which implements the de-sired dynamics. In practice, the gradient magnetic pulse in the first cycle can be simulated with quadrupole fieldsas ~B = ( x, y, − z ). When the condensates is stronglyconfined in the xy plane, the influence of the nonzero gra-dient along z -axis can be neglected. The effective gradi-ent coupling in the second cycle can be implemented withthe help of atom-laser coupling. For instance, a standardtwo-set Raman beams with blue-detuning [48] can realizean effective couplingΩ[sin( ~k · ~r ) σ x + sin( ~k · ~r ) σ y ] , (12)where the wavevectors ~k and ~k in the xy plane can bechosen as ~k = k (cos θ, − sin θ ) and ~k = k (sin θ, cos θ ).When 2 π/k is much larger than the trap length l T , therequired effective coupling is approximately obtained. Fi-nally, the phases discussed in the context can be detectedby monitoring their corresponding density and momen-tum distributions using the setup of time of flight. VI. CONCLUDING REMARKS
To summarize, we have discussed the ground statephase diagram of SO coupled BECs subject to gradientmagnetic fields. Theoretical and numerical analyses indi-cate that the system supports various interesting vortexphysics, including the single-particle giant-vortex stateswith tunable vorticity, multiple layered vortex-lattice-ring states, and the ellipse-shaped vortex profiles. There-fore, the combination of SO coupling and the gradientmagnetic fields provides a powerful method to engineervarious vortex states without rotating the trap. We hopeour work will stimulate further research of searching forvarious novel states in SO coupled bosons subject to ef-fective gradient magnetic fields.
VII. ACKNOWLEDGEMENT
X.F. Z., Z.W. Z., and G.C. G. acknowledge the supportby NSFC (Grant Nos. 11004186, 11474266,11174270),National Basic Research Program of China (GrantsNo. 2011CB921204 and No. 2011CBA00200). C. W.is supported by the NSF DMR-1410375 and AFOSRFA9550-14-1-0168, and also acknowledges the supportfrom the National Natural Science Foundation of China(11328403). [1] I. ˇZuti´c, J. Fabian, and S. Das Sarma, Rev. Mod. Phys. , 323 (2004).[2] M.Z. Hasan and C.L. Kane, Rev. Mod. Phys. , 3045(2010).[3] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. , 1057(2011).[4] R.P. Feynman, Statistical Mechanics, A Set of Lec- tures (Addison-Wesley Publishing Company, ADDRESS,1972).[5] C. Wu, Mod. Phys. Lett. B , 1 (2009).[6] C. Wu, I. Mondragon-Shem, arXiv:0809.3532; C. Wu, I.Mondragon-Shem, and X.-F. Zhou, Chin. Phys. Lett. ,097102 (2011).[7] T.D. Stanescu, B. Anderson, and V. Galitski, Phys. Rev. A , 023616 (2008).[8] C. Wang, C. Gao, C.M. Jian, and H. Zhai, Phys. Rev.Lett. , 160403 (2010).[9] T.L. Ho and S. Zhang, Phys. Rev. Lett. , 150403(2011).[10] H. Hu, B. Ramachandhran, H. Pu, and X.J. Liu, Phys.Rev. Lett. , 10402 (2012).[11] S. Sinha, R. Nath, and L. Santos, Phys. Rev. Lett. ,270401 (2011).[12] Y. Li, X. Zhou, and C. Wu, arXiv:1205.2162 (2012).[13] T. Kawakami, T. Mizushima, M. Nitta, and K. Machida,Phys. Rev. Lett. , 015301 (2012).[14] S. Gopalakrishnan, I. Martin, and E.A. Demler, Phys.Rev. Lett. , 185304 (2013).[15] Y. Lin, R.L. Compton, K. Jim´enez-Garc´ıa, J.V. Porto1,and I.B. Spielman, Nature , 628 (2009); Y.-J. Lin,R.L. Compton, A.R. Perry, W.D. Phillips, J.V. Porto,and I.B. Spielman, Phys. Rev. Lett. , 130401 (2009);Y. Lin, K. Jimenez-Garcia, and I. Spielman, Nature ,83 (2011).[16] J.-Y. Zhang, S.-C. Ji, Z. Chen, L. Zhang, Z.-D. Du, B.Yan, G.-S. Pan, B. Zhao, Y.-J. Deng, H. Zhai, S. Chen,and J.-W. Pan, Phys. Rev. Lett. , 115301 (2012).[17] P. Wang, Z.-Q. Yu, Z. Fu, J. Miao, L. Huang, S. Chai, H.Zhai, and J. Zhang, Phys. Rev. Lett. , 095301 (2012).[18] L.W. Cheuk, A.T. Sommer, Z. Hadzibabic, T. Yefsah,W.S. Bakr, and M.W. Zwierlein, Phys. Rev. Lett. ,095302 (2012)[19] C. Qu, C. Hamner, M. Gong, C. Zhang, and P. Engels,Phys. Rev. A , 021604(R) (2013).[20] A.J. Olson, S.-J. Wang, R.J. Niffenegger, C.-H. Li, C.H.Greene, Y.P. Chen, Phys. Rev. A , 013616 (2014).[21] R.M. Wilson, B.M. Anderson, and C.W. Clark, Phys.Rev. Lett. , 185303 (2013)[22] V. Achilleos, D.J. Frantzeskakis, P.G. Kevrekidis, andD.E. Pelinovsky, Phys. Rev. Lett. , 264101 (2013).[23] Y.V. Kartashov, V.V. Konotop, and F.K. Abdullaev,Phys. Rev. Lett. , 060402 (2013).[24] T. Ozawa and G. Baym, Phys. Rev. Lett. , 085304(2013).[25] Y. Deng, J. Cheng, H. Jing, C.-P. Sun, and S. Yi, Phys.Rev. Lett. , 125301 (2012).[26] Y. Zhang, L. Mao, and C. Zhang, Phys. Rev. Lett. ,035302 (2012).[27] V. E. Lobanov, Y.V. Kartashov, and V.V. Konotop,Phys. Rev. Lett. , 180403 (2014).[28] X. Luo, L. Wu, J. Chen, R. Lu, R. Wang, and L. You, arXiv:1403.0767v2; L. He, A. Ji, and W. Hof-stetter, arXiv:1404.0970v1; Zeng-Qiang Yu, 1407.0990v1;W. Han, G. Juzeli¨unas, W. Zhang, and W.-M. Liu,arXiv:1407.2972v1.[29] X. Zhou, Y. Li, Z. Cai, and C. Wu, J. Phys. B: At. Mol.Opt. Phys. , 134001 (2013).[30] J. Dalibard, F. Gerbier, G. Juzeli¨unas, and P. Ohberg,Rev. Mod. Phys. , 1523 (2011).[31] H. Zhai, Int. J. Mod. Phys. B. , 1230001 (2012); H.Zhai, arXiv:1403.8021.[32] V. Galitski and I. B. Spielman, Nature , 49 (2013).[33] N. Goldman, G. Juzeli¨unas, P. Ohberg, I. B. Spielman,arXiv: 1308.6533.[34] B.M. Anderson, I.B. Spielman, and G. Juzeli¨unas, Phys.Rev. Lett. , 125301 (2013).[35] Z.-F. Xu, L. You, and M. Ueda, Phys. Rev. A , 063634(2013).[36] C. J. Kennedy, Georgios A. Siviloglou, H. Miyake, W.C.Burton, and W. Ketterle, Phys. Rev. Lett. , 185301(2013).[38] V. Pietila, and M. Mottonen, Phys. Rev. Lett. ,030401 (2009).[39] Y. Kawaguchi, M. Nitta, and M. Ueda, Phys. Rev. Lett. , 180403 (2008).[40] V. Schweikhard, I. Coddington, P. Engels, S. Tung, andE.A. Cornell, Phys. Rev. Lett. , 210403 (2004).[41] E.J. Mueller and T.-L. Ho, Phys. Rev. Lett. , 180403(2002).[42] A.L. Fetter, Rev. Mod. Phys. , 647 (2009).[43] X. F. Zhou, J. Zhou, and C. Wu, Phys. Rev. A , 063624(2011).[44] X.-Q. Xu and J. H. Han, Phys. Rev. Lett. , 200401(2011).[45] J. Radi´c, T.A. Sedrakyan, I.B. Spielman, and V. Galitski,Phys. Rev. A , 063604 (2011).[46] A. Aftalion and P. Mason, Phys. Rev. A , 023610(2013).[47] A.L. Fetter, Phys. Rev. A , 023629 (2014).[48] X.-J. Liu, K.T. Law, and T.K. Ng, Phys. Rev. Lett. ,086401 (2014).[49] C.J. Pethick and H. Smith, BoseCEinstein Condensationin Dilute Gases, Chapter 7, Cambridge University Press,Cambridge, 2001. Appendix A: Single particle eigenstates for large β We start with the dimensionless Hamiltonian H ~ ω = Z d~ρ ˆ φ ( ~ρ ) † n − ~ ∇ βρ (cos θ ˆ r + sin θ ˆ ϕ ) · ~σ + α~k · ~σ + 12 ρ o ˆ φ ( ~ρ ) . (A1)Since the total angular momentum is conserved, the single particle eigenstates can be written as φ m =[ f ( ρ ) e imϕ , g ( ρ ) e i ( m +1) ϕ ] T with j z = ( m + 1 / ( ˆ p ρ j z ρ + ( ρσ x − β ) − α [(ˆ p ρ cos θ + j z ρ sin θ ) σ x +( j z ρ cos θ − ˆ p ρ sin θ ) σ y ] + j z σ z ρ (cid:27) (cid:20) ˜ f ( ρ )˜ g ( ρ ) (cid:21) = E (cid:20) ˜ f ( ρ )˜ g ( ρ ) (cid:21) , where ˜ f ( ρ ) = f ( ρ ) e iθ/ , ˜ g ( ρ ) = g ( ρ ) e − iθ/ , ˆ p ρ = − i ( ∂∂ρ + ρ ) is the momentum operator along the radical direction. Forlarge β ≫
1, these functions mainly distribute around the circle ρ = β in the plane, so we consider the superposition F ± ( ρ ) = [ ˜ f ( ρ ) ± ˜ g ( ρ )]), which satisfies the following approximated equations as ˆ p ρ ∓ α cos θ ˆ p ρ + j z ρ ∓ α sin θ j z ρ + ρ ∓ βρ ! F ± ( ρ ) ± iα (cid:18) ˆ p ρ sin θ − j z ρ cos θ (cid:19) F ∓ ( ρ ) = E j z F ± ( ρ ) . The above equation indicates that to minimize the kinetic energy, we need h ~p ρ i ≃ α cos θ . Around ρ = β , wehave the approximated solutions as F ± ( ρ ) ∼ H n ( ρ ± β ) e − ( ρ ± β ) / e ± iα cos θρ with H n ( r ) the usual n -th Hermitepolynomial. Therefore, F + is negligible since we always have ρ >
0. The solution now can be written as ˜ f ( ρ ) ≃ ˜ g ( ρ ) ∝ H n ( ρ − β ) e − ( ρ − β ) / e iα cos θρ . So we obtain the approximated wavefunctions for the lowest band ( n = 1) as φ n =1 ,j z ≃ π ) ρ e − ( ρ − β )22 e iρα cos θ (cid:20) e i [ mϕ − θ ] − e i [( m +1) ϕ + θ ] (cid:21) . (A2)The dispersion is estimated as E n,j z = n + 1 − α − β j z − αβ sin θ ) α cos θ + β ) , (A3)which is minimized when j z ≃ αβ sin θ , so for the kinetic term along the tangential direction ˆ E ϕ = ( j z ρ − α sin θ ) / Appendix B: Energy estimation of vortix lattice states around the ring
For weak interaction, the condensates expands along the radial direction as the parameter g is increased. When g is large enough, to lower the kinetic energy, the system tends to involve vortices located around a ring inside thecondensates, which separate the wavefunction into two parts. Inside the vortex-ring, the wavefunction for the spin-upcomponent is approximated as a giant vortex with the phase factor e i πm − φ , while outside the ring, the mean angularmomentum carried by single particle is approximated as m + ~ . The difference q = m + − m − represents the vortexnumber inside the condensates. Therefore the variational ground-state can be approximated as follows φ ( ρ, φ ) = (cid:26) φ j z − q/ ( ρ, φ ) when ρ ∈ (0 , β ) ,φ j z + q/ ( ρ, φ ) when ρ ∈ ( β, ∞ ) . (B1)We also assumes that around the circle ρ = β , vortices are involved and self-organized to compensate the phasemismatch so that the the whole wavefunction is well-defined.The increase of the kinetic energy mainly comes from the term ˆ E ϕ = ( j z ρ − α sin θ ) /
2. For ground state withoutlattices, the energy is estimated and denoted as h ˆ E ϕ i j z = Z ∞ dρ Z π dϕρφ † j z ( j z /ρ − α sin θ ) φ j z , (B2)where we use hi j z to denote the mean values over the trivial variational function φ j z Similarly, for GV states andlarge β ≫
1, the corresponding energy increasing is estimated as h ˆ E ϕ i = Z ∞ β dρ Z π dϕρφ † j z + q [( j z + q ) /ρ − α sin θ ] φ j z + q + Z β dρ Z π dϕρφ † j z − q [( j z − q ) /ρ − α sin θ ] φ j z − q ≈ h ˆ E ϕ i j z + h q ρ i j z + Z ∞ dρ Z π dϕρφ † j z q ρ ( j z ρ − α sin θ ) φ j z ≈ h ˆ E ϕ i j z + q β ( q − σα sin θ √ π ) , (B3)which gives the Eq. (8) shown in the context.Around each vortex core, both the density and phase are twisted as shown in Fig. 5. To take into account theenergy contribution of these vortices, we choose the following approximated local wavefunction as (see Fig.(8a) for θ = π/ ψ ≃ φ (cid:20) e − ik + x e −√ y/ξ + 1 + e − i ( k − x + π ) e √ y/ξ + 1 (cid:21) (cid:20) (cid:21) , where for simplicity we have set the origin to the vortex core, φ is the bulk wavefunction away from the vortex coresand is estimated as | φ | = √ n = [2 π √ σρ ] − , k ± represent the corresponding local wavevectors of the condensatesinwards and outwards the vortex ring. ξ = (2 gn ) − / is usually called as the healing length as can be seen byconsidering the above wavefunction along the line x = 01 e −√ y/ξ + 1 − e √ y/ξ + 1 = tanh[ y/ ( √ ξ )] , which is consistent with the estimation shown in [49]. We note that the approximated wavefuntion is only valid when k − x ∈ [ − π, π ] and y ∈ [ − ξ, ξ ]. We also requires k + ≃ k − ≃ α such that a vortex is formed around the origin asindicated by numerics. Since energy gain along the tangential direction has been taking into account in ˆ E ϕ , in thistypical case, the changing of kinetic energy along the direction perpendicular to the local wavevector ~k can then beapproximated by Z k − x = πk − x = − π dx Z ξ − ξ dyψ † h − ∇ y i ψ = √ √ √ πβσαξ ≈ . βσαξ . We note that the above analysis also applies for θ = 0. Taking into account all vortices inside the two-componentscondensates, we obtain the total energy introduced due to fluctuation of the vortex density profile as 0 . q/ ( βσαξβσαξ