Vortex lattice solutions to the Gross-Pitaevskii equation with spin-orbit coupling in optical lattices Hidetsugu Sakaguchi and Ben Li Department of Applied Science for Electronics and Materials, Interdisciplinary Graduate School of Engineering Sciences, Kyushu University, Kasuga, Fukuoka 816-8580, Japan Effectivespin-orbitcouplingcanbecreatedincoldatomsystemsusingatom-lightinteraction. We 3 study the BECs in an optical lattice using the Gross-Pitaevskii equation with spin-orbit coupling. 1 Blochstatesforthelinearequationarenumericallyobtained,andcomparedwithstationarysolutions 0 totheGross-Pitaevskiiequationwithnonlinearterms. Variousvortexlatticestatesarefoundwhen 2 thespin-orbit coupling is strong. n a PACSnumbers: 03.75.-b,03.75.Mn,05.30.Jp,67.85.Hj J 6 1 Recently, Bose-Einstein condensates (BECs) with effective spin-orbit coupling were created in cold atom systems usingatom-lightinteraction[1]. Thespin-orbit-coupledBECsareactivelystudiedtheoretically[2]. Wangetal. found ] that the mean-field ground state has two different phases: plane-wave and stripe phases depending on the nonlinear s a interactions [3]. Half vortex states were found in a spin-orbit coupled BECs confined in a harmonic potential [4, 5]. g Exotic spin textures were predicted in Bose-Hubbard models corresponding to spin-orbit coupled BECs in the Mott- - insulator phase [6, 7]. t n The GRoss-Pitaevskii (GP) equation is a mean-field approximation for the BECs with the spin-orbit coupling. a There are some studies for the GP equation with spin-orbitcoupling in optical lattices [8, 9]. In this paper, we study u vortex lattice solutions to the GP equation in a square type optical lattice. The model equation is expressed as q . mat i∂∂ψt+ = −12∇2ψ++(g|ψ+|2+γ|ψ−|2)ψ+−ǫ{cos(2πx)+cos(2πy)}ψ++λ(cid:18)∂∂ψx− −i∂∂ψy−(cid:19), d- i∂ψ− = 1 2ψ−+(g ψ− 2+γ ψ+ 2)ψ− ǫ cos(2πx)+cos(2πy) ψ−+λ ∂ψ+ i∂ψ+ , (1) n ∂t −2∇ | | | | − { } (cid:18)− ∂x − ∂y (cid:19) o c where ψ = (ψ+,ψ−) denotes the wave function of the spinor BECs, ǫ is the strength of the optical lattice, g and γ [ express the strengths of interactions respectively between the same and the different kinds of atoms, and λ denotes the strength of the Rashba spin-orbit coupling. We have assumed that the wavelength of the optical lattice is 1. 1 If g and γ are zero, Eq. (1) becomes linear equations with spatially-periodic potential. The Bloch states are v 5 stationary solutions to the linear equations, which are expressed as 6 5 ψ+(x,y,t)=φ+(x,y)exp(ikxx+ikyy iµt), ψ−(x,y,t)=φ−(x,y)exp(ikxx+ikyy iµt), (2) − − 3 1. where φ+ and φ− are periodic functions of wavelength 1. Therefore, φ+ and φ− satisfy 0 13 µφ+ = −21∇2φ++ 12(kx2+ky2)φ+−ikx∂∂φx+ −iky∂∂φy+ −ǫ{cos(2πx)+cos(2πy)}φ+ : v ∂φ− ∂φ− +λ i +ik φ +k φ , i (cid:18) ∂x − ∂y x − y −(cid:19) X 1 1 ∂φ ∂φ r µφ = 2φ + (k2+k2)φ ik − ik − ǫ cos(2πx)+cos(2πy) φ a − −2∇ − 2 x y −− x ∂x − y ∂y − { } − ∂φ+ ∂φ+ +λ i ikxφ++kyφ+ . (3) (cid:18)− ∂x − ∂y − (cid:19) The eigenvalue µ and the eigen function φ+ and φ− can be numerically obtained from the stationary solution of the 2 (b) (c) (a) -12.8(cid:13) -1.6(cid:13) -4(cid:13) -13(cid:13) -1.9(cid:13) -5(cid:13) -13.2(cid:13) -2(cid:131)˚.2(cid:13) (cid:131)˚-13.4(cid:13) (cid:131)˚-6(cid:13) -2.5(cid:13) -13.6(cid:13) -7(cid:13) -2.8(cid:13) -13.8(cid:13) -8(cid:13) 0(cid:13) 1(cid:13) 2(cid:13) 3(cid:13) 4(cid:13) 5(cid:13) 6(cid:13) 0(cid:13) 1(cid:13) 2(cid:13) 3(cid:13) 4(cid:13) 5(cid:13) 6(cid:13) 0(cid:13) 1(cid:13) 2(cid:13) 3(cid:13) 4(cid:13) 5(cid:13) 6(cid:13) k(cid:13) k(cid:13) k(cid:13) x x x FIG. 1: Eigenvalue µvs. kx for (a) λ=π/2,ky =0, (b) λ=3π/2,ky =0, and (c) λ=π,ky =kx. Dashed curvein Fig.1(a) is plotted usingEq. (8) and the dashed curvein Fig. 1(c) is obtained using Eq. (6). linear equation [10]: ∂∂φt+ = 12∇2φ+− 21(kx2+ky2)φ++ikx∂∂φx+ +iky∂∂φy+ +ǫ{cos(2πx)+cos(2πy)}φ+ ∂φ ∂φ − − λ i +ikxφ−+kyφ− +µφ+, − (cid:18) ∂x − ∂y (cid:19) ∂φ 1 1 ∂φ ∂φ − = 2φ (k2+k2)φ +ik − +ik − +ǫ cos(2πx)+cos(2πy) φ ∂t 2∇ −− 2 x y − x ∂x y ∂y { } − ∂φ+ ∂φ+ λ i ikxφ1+kyφ+ +µφ−, − (cid:18)− ∂x − ∂y − (cid:19) dµ = α(N0 N), (4) dt − where α>0 is a parameterand fixed to be 5 in our numericalsimulation. N = 01 01(|φ+|2+|φ−|2)dxdy is the total norm, and N0 is fixed to be 1 by the normalization condition. The time evoluRtioRn of the dissipative equation (4) leads to a stationary state and the total norm N approaches N0 =1. The eigenvalue µ in Eq. (3) is obtained as µ in Eq. (4) at the stationary state. In this numerical method, the ground state for fixed values of k and k is obtained x y atthe stationarystate,starting frommostinitial conditions,because the totalenergydecreasesinthe time evolution of Eq. (4). Excited states are obtained by removing the ground state by the method of orthogonalization. Figure 1(a) shows µ(k ) as a function of k for k = 0,ǫ = 5, and λ = π/2. µ(k ) is a periodic function of k with period x x y x x 2π. There are peaks near k = 0,π and 2π and minima at k π/2 and 3π/2. The peak point at k = π is a cusp x x x ∼ point, where two µ(k ) curves corresponding to the ground state and the excited state cross, although the branch x of the excited state is not shown. For λ = 0, µ(k ) increases monotonously as k increases from 0 and reaches the x x maximum at the edge of the Brillouin zone at k =π. If there is no optical lattice, i.e., ǫ=0, µ(k) takes a minimum x at k =λ where k = k2+k2 [2, 3]. The minimum point for λ=π/2 locates near k =λ, and the peak corresponds x y x q to the edge of the Brillouin zone. Figure 1(b) shows µ(k ) at λ = (3/2)π. There are a large peak at k = 0 and x x 2π and a small peak at k = π and minima at k = 2 and k = 4.3. Figure 1(c) shows µ(k) as a function of k for x x x x k = k ,λ =π,ǫ = 5. There is a large peak at k = 0 and 2π, a small peak at k =π, and minima at k 2.5 and x y x x x ∼ 3.8. The wavenumber k 2.5 is close to λ/√2 2.22 by the simplest approximation k = λ but slightly deviated. x ∼ ∼ The approximation kmin =λ for the minimum point of µ(k) becomes worse for large λ. The small peaks correspond to the edge of the Brillouin zone. Figure 2(a) shows φ+(x,y) and φ−(x,y) as a function of y in the section x=0 at kx =3π/2 for λ=3π/2. The | | | | modulus φ+ and φ− take maximum at different positions. The minimum value is almost zero, which implies the | | | | existence of vortices. Figure 2(b) shows a contour plot of φ+ for the same parameter. The locations of vortices for φ+ can be calculated from the phase distribution θ+(x,y)|=s|in−1(Imφ+(x,y)/φ+(x,y)). There exist a vortex at a | | point,ifthe integralofthe phasegrandientalongananticlockwisepathencirclingthe pointis anontrivialmultiple of 2π. We have counted the path integral by discretizing the (x,y) space with ∆x =1/64. Figure 2(c) shows positions of vortices of vorticity 1 with square and marks. The vortex cores locate near (0, 0.28) and (0, 0.48) for φ+. ± × − − In generic cases, there is a vortex of vorticity 1 or -1 at a position satisfying φ+ = 0, where a line of Reφ+ = 0 | | intersects with a line of Imφ+ =0 [11]. We do not show explicitly the positions of vortices later in Fig. 3 and Fig. 4, however, we have checked the existence of vortices of vorticity 1 or -1 at positions satisfying φ =0 by calculating ± | | 3 (a) (b) (c) (d) 0.4(cid:13) 1.5 0.06 0.2(cid:13) |(cid:131)(cid:211)| 1 y(cid:13) 0(cid:13) (cid:211)|m0.04 |(cid:131) 0.5 -0.2(cid:13) 0.02 -0.4(cid:13) 0 0 -0.25 0 0.25 -0.4(cid:13) -0.2(cid:13) 0(cid:13) 0.2(cid:13) 0.4(cid:13) 3.6 3.8 4 4.2 4.4 y x(cid:13) (cid:131)(cid:201) FIG. 2: (a) φ+ (solid curve) and φ− (dashed curve) along the line x = 0 for λ = 3π/2, and kx = 3π/2. (b) Contour plot | | | | of φ+. (c) Square shows a vortex with vorticity 1 and shows a vortex with vorticity -1. (d) Minimum values of φ+ as a | | × | | function of λ for kx =λ and ky =0. the phase distribution. Figure 2(d) shows the minimum value of φ+ as a function of λ for kx = λ,ky = 0 at ǫ = 5. | | The minimum value becomes zero and a vortex-antivortexpair appears for λ>λ 4.2. c ∼ Because φ are periodic functions with wavelength 1, φ can be expressed as the simplest approximation: ± ± φ+ = C0++C1+e2πix+C2+e−2πix+C3+e2πiy +C4+e−2πiy, φ− = C0−+C1−e2πix+C2−e−2πix+C3−e2πiy+C4−e−2πiy. (5) Substitution of this ansatz into Eq. (3) yields µC0± = (kx2+ky2)C0±/2−(ǫ/2)(C1±+C2±+C3±+C4±)+λ(±ikx+ky)C0∓, µC1± = {(kx+2π)2+ky2}C1±/2−(ǫ/2)C0±+λ{±i(kx+2π)+ky}C1∓, µC2± = {(kx−2π)2+ky2}C2±/2−(ǫ/2)C0±+λ{±i(kx+2π)+ky}C2∓, µC3± = {kx2+(ky +2π)2}C3±/2−(ǫ/2)C0±+λ{±ikx+(ky +2π)}C3∓ µC4± = {kx2+(ky −2π)2}C4±/2−(ǫ/2)C0±+λ{±ikx+(ky −2π)}C4∓. (6) For ky =0, C0− =iC0+,C1− =iC1+,C2− =iC2+ are satisfied, and then (ǫ/2)C0+ C1+ = µ (k +−2π)2/2+λ(k +2π), x x − (ǫ/2)C0+ C2+ = µ (k −2π)2/2+λ(k 2π), x x − − − C3+ = [−(ǫ/2){µµ−((kkx22++44ππ22))//22)}2+(λǫ/2(2k)λ2(+kx4−π22)πi)]C0+, { − x } − x C4+ = [−(ǫ/2){µµ−((kkx22++44ππ22))//22)}2+(λǫ/2(2k)λ2(+kx4+π22)πi)]C0+, (7) { − x } − x where µ is given by a solution of the equation k2 ǫ2/4 ǫ2/4 µ = x λk + + 2 − x µ (k +2π)2/2+λ(k +2π) µ (k 2π)2/2+λ(k 2π) x x x x − − − − ǫ2 2µ (k2+4π2) 2λk + − x − x . (8) 4 µ (k2+4π2)/2 2 λ2(k2+4π2) { − x } − x Furthermore, C3− = iC3∗+,C4− = iC4∗+ are satisfied. Here, ∗ denotes the complex conjugate. The dashed curve in Fig. 1(a)denotes µ(k ) by Eq.(8) at λ=π/2. The approximationis goodfor λ=π/2 but is not so goodfor large λ, x because the higher harmonics is necessary for the expansion in Eq. (5). We can assume that C0+,C1+ and C2+ are real numbers and C4+ =C3∗+ =ReC3+ iImC3+. Then, φ+ and φ− are expressed as − φ+ = C0++(C1++C2+)cos(2πx)+i(C1+ C2+)sin(2πx)+2ReC3+cos(2πy) 2ImC3+sin(2πy), − − φ− = iC0++i(C1++C2+)cos(2πx) (C1+ C2+)sin(2πx)+2iImC3+sin(2πy)+2iReC3+cos(2πy). (9) − − 4 (a) (b) (c) (d) 0.06 0.04 m | (cid:211) |(cid:131) 0.02 0 2 2.2 2.4 2.6 2.8 3 3.2 3.4 (cid:131)(cid:201) FIG.3: (a) Contourplot of ψ+ forλ=π,g=1,γ =0.5 andL=8. kx =ky areevaluated as3π/4. (b)Contourplot of ψ− . | | | | (c) Contourplot of φ+ tothelinear equationEq.(3) for λ=π,kx =ky =3π/4. (d)Minimum valuesof φ+ asafunction of | | | | λ for kx =ky =λ/√2. (a) (b) (c) (d) FIG. 4: (a) Contour plot of ψ+ for λ=π,g =1,γ =2 and L=8. kx =ky are evaluated as 3π/4. (b) Contour plot of ψ− . (c) Contourplot ofthesuper|pos|ition of (φ+++φ+−)/√2 tothelinear equation Eq.(3)for λ=π and kx =ky = 3π/4|. (d|) Contour plot of the superposition of (φ+|++φ+−)/√2 to|thelinear equation Eq. (3) for λ=1 and kx =ky = π/±4 | | ± Imφ+ =0 and Reφ− =0 are satisfied on the line x=0. When λ is small, the minimum values of Reφ+ and Imφ− are positive and there are no vortices. When λ is increased the minimum values decrease and reach 0, and then a vortex pair is created. A vortex core of φ+ is located at a point on the line x= 0 where Reφ+(0,y)=0 is satisfied, and similarly a vortex core of φ is located at a point on the line x=0 where Imφ (0,y)=0 is satisfied. − − Even for g and γ is not zero, the Bloch state is a good approximation for the stationary state for γ < g. We haveperformednumericalsimulationofEq.(1) by the imaginarytime evolutionmethodsimilar to Eq.(4)andfound stationary solutions. The system size is Lx×Ly = L×L and the total norm N = 0L 0L(|ψ+|2+|ψ−|2)dxdy is set to be L2 in this paper. Periodic boundary conditions are imposed. The potential iRs shRifted as U = ǫ[cos 2π(x − { − 1/2) +cos 2π(y 1/2) ] by (1/2,1/2)to confine the wave pattern in the range of [0,L] [0,L]. } { − } × Figure 3(a) and (b) show contour plots of ψ+ and ψ− at g = 1,γ = 0.5,L = 8,λ = π, and ǫ = 5. The contour | | | | plot is drawn in the region [0,2] [0,2], and the contour lines are drawn for ψ = 0.025,0.05,0.075,0.1,1and 1.5. ± × | | Vortex pairs exist in each cell of size 1 for this parameter, and a vortex lattice is constructed as a whole. Vortex lattices were experimentally found first in rotating BECs [12] and recently in BECs under synthetic magnetic fields by atom-light interaction [13]. In our model equation, vortices are spontaneously created by the spin-orbit coupling. The wavevector(kx,ky) is evaluated at (3π/4,3π/4). Positions of vortex cores for ψ+ and ψ− are mutually deviated. Figure 3(c) shows a contour plot of φ+ for the linear equation corresponding to g =0,γ =0 for kx =ky =3π/4 at | | λ = π and ǫ = 5. The eigenvalue µ takes a minimum at (k ,k ) = (3π/4,3π/4) in the finite size system of L = 8, x y where k (k ) takes a discrete value 2πn /L (2πn /L) with integer n (n ). The contour plot is almost the same as x y x y x y Fig. 3(a). It means that the Bloch wave is a good approximation for the solution to the GP equation. Figure 3(d) shows the minimum values of φ+ for the linear equation as a function of λ for kx = ky = λ/√2 at ǫ = 5. The | | minimum value becomes zero and vortices appear for λ>2.9. It is related to the existence of vortices at λ=π. Stripe wave states are expected to appear for γ >g. The superposition of Bloch waves of (k ,k ) and ( k , k ) x y x y − − is a simple approximation for γ >g. Figure 4(a) and (b) show contour plots of ψ+ and ψ− at g =1,γ =2,L=8, | | | | 5 (a) (b) (c) 8(cid:13) 8(cid:13) 8(cid:13) 6(cid:13) 6(cid:13) 6(cid:13) j(cid:13) 4(cid:13) j(cid:13) 4(cid:13) j(cid:13) 4(cid:13) 2(cid:13) 2(cid:13) 2(cid:13) 0(cid:13) 0(cid:13) 2(cid:13) 4(cid:13) 6(cid:13) 8(cid:13) 2(cid:13) 4(cid:13) 6(cid:13) 8(cid:13) 2(cid:13) 4(cid:13) 6(cid:13) 8(cid:13) i(cid:13) i(cid:13) i(cid:13) FIG.5: (a)Spinconfigurationof(sx(i,j),sy(i,j))atλ=π,g=1,γ =0.5andL=8. (b)Spinconfigurationof(sx(i,j),sy(i,j)) at λ=π,g=1,γ =2 and L=8. (c) Spin configuration of sz(i,j) at λ=π,g=1,γ =2 and L=8. and λ = π. The wavevector is evaluated as (k ,k ) = (3π/4,3π/4) in this case, too. The contour lines are drawn x y for ψ =0.025,0.05,0.075,0.1,1and 1.5. Vortex cores exist in dark pointed regions. The vortex lattice structure is ± | | rather complicated. The circular contour lines correspond to peak regions of ψ . The peak regions stand in a line ± | | in the direction of angle π/4 and the peak lines for ψ+ and ψ− alternates in the diagonal direction of angle π/4. − Figure 4(c) shows a contour plot of a linear combination (φ++ +φ+−)/√2 of two Bloch waves φ++ and φ+− with | | (k ,k ) = (3π/4,3π/4) and ( 3π/4, 3π/4) for the + component at λ = π. The superposition of the Bloch waves x y − − is a good approximation for the stationary solution to the GP equation. The superposition of two Bloch waves with opposite wavevectors generates a standing wave. For plane waves, the amplitude becomes zero at the nodal lines. The nodal lines are perturbed by the optical lattice and vortices are generated. A vortex lattice structure therefore appears even for small λ in case of γ > g. Figure 4(d) shows a vortex lattice pattern with k = k = π/4 at λ = 1 x y and ǫ=5. For large λ, a vortex pair is created in a single Bloch wave and the superposition of the two Bloch waves make the vortex lattice structure more complicated as shown in Fig. 4(c). The complicated patterns might be simplified, if a spin representation is used, which was discussed in the Bose- Hubbard model [6, 7]. The whole system is divided into cell regions of [i 1,i] [j 1,j]. The spin variables − × − s (i,j),s (i,j) and s (i,j) are defined for each cell labeled by (i,j) as x y z i j i j sx(i,j) = Zi−1Zj−1ψ†σxψdxdy =Zi−1Zj−1(ψ+∗ψ−+ψ−∗ψ+)dxdy, i j i j sy(i,j) = Zi−1Zj−1ψ†σyψdxdy =Zi−1Zj−1(−iψ+∗ψ−+iψ−∗ψ+)dxdy, i j i j sz(i,j) = ψ†σzψdxdy = (ψ+ 2 ψ− 2)dxdy, (10) Zi−1Zj−1 Zi−1Zj−1 | | −| | where σ ,σ and σ are the Paulimatrix, and † denotes the complex conjugate transpose. Figure 5(a) shows (s ,s ) x y z x y correspondingto the pattern in Figs. 3(a) and (b) for g =1,γ =0.5,λ=π and ǫ=5. The vector (s (i,j),s (i,j)) is x y expressed as an arrow on each lattice point at (i 1/2,j 1/2). The pattern is interpreted as a ferromagnetic state − − in the (x,y) plane in this spin representation. The spin s is zero for this pattern. Figures 5(b) and (c) show spin z configurationsrespectivelyfor(s ,s )ands forthepatternatg =1andγ =2showninFigs.4(a)and(b). Thespin x y z configurationis also rather complicated. The wavelengthof the spin configurationis 4 both in the i and j directions. An anti-ferromagnetic order is seen in the diagonal direction of angle π/4 and a ferromagnetic order appears in its orthogonaldirectionofangle π/4both for the (s ,s ) ands patterns. The (s ,s ) componentappearsat the sites x y z x y − where the s component vanishes, and the s component appears at the sites where the (s ,s ) component vanishes. z z x y Tosummarize,wehavestudiedtheGross-Pitaevskiiequationwithspin-orbitcouplinginanopticallattice. Wehave found that a vortex lattice structure appears for large λ in case of γ <g. A vortex lattice structure appears even for smallλincaseofγ >g,becausethenodallinesinthestripewavepatternareperturbedbytheopticallattice. Wehave found a complicated spin configurationin a case of γ >g. The complicated patterns can be qualitatively understood by the correspondingBlochwaves. The Blochwavesarefurther approximatedby a Fourierseries expansionwith five modes to understand the formation of the vortices. We have obtained various spin configurations by changing the 6 parameter λ. The detailed phase diagrams by changing various parameters are under study. [1] Y.-J. Lin,K. Jim´enez-Garc´ıa, and I. B. Spielman, Nature471, 83 (2011). [2] H.Zhai, Int.J. Mod. Phys. 26, 1230001 (2012). [3] C. Wang, C. Gao, C. -M. Jian, and H.Zhai, Phys.Rev.Lett. 105, 160403 (2010). [4] B. Ramachandhran,B. Opanchuk,X.-J. Liu, H.Pu, P. D.Drummond,and H.Hu,Phys. Rev.A 85, 023606 (2012). [5] Y.Zhang, L. Mao, and C. Zhang, Phys.Rev. Lett.108, 035302 (2012). [6] J. Radic, A. DiCiolo, K. Sun,and V.Galitski, Phys.Rev.Lett. 109, 085303 (2012). [7] W. S.Cole, S. Zhang, A.Paramekanti, and N. Trivedi, Phys.Rev.Lett. 109. 085302 (2012). [8] J. Larson and E. Sjo¨qvist, Phys.Rev. A.79, 043627 (2009). [9] Y.Zhang and C. Zhang, arXiv:1203.2389 (2012). [10] H.Sakaguchi and H.Takeshita, J. Phys.Soc. Jpn. 77, 054003 (2008). [11] A.Ohta, R.Kashiwa, and H. Sakaguchi, Phys.Rev.A 82, 055602 (2010). [12] J. E. Williams and M. J. Holland, Nature401, 568 (1999). [13] Y.-L.Lin, R.L. Compton, K.Jim´enez-Garc´ıa, J. V.Porto,and I.B. Spielman, Nature462, 628 (2009)