Microscopic Wave Functions of Spin Singlet and Nematic Mott States of Spin-One Bosons in High Dimensional Bipartite Lattices Michiel Snoek and Fei Zhou ∗ ITP, Utrecht University, Minnaert building, Leuvenlaan 4, 3584 CE Utrecht, The Netherlands (Dated: February 2, 2008) We present microscopic wave functions of spin singlet Mott insulating states and nematic Mott insulating states. We also investigate quantum phase transitions between the spin sin- glet Mott phase and the nematic Mott phase in both large-N limit and small-N limit (N being the number of particles per site) in high dimensional bipartite lattices. In the mean field approximationemployedinthisarticlewefindthatphasetransitionsaregenerallyweaklyfirstorder. PACS numbers: 03.75.Mn, 75.10.Jm, 75.45.+j 4 0 0 I. INTRODUCTION dimensional Mott states for both even and odd numbers 2 of particles per site are presented in Ref. 6. n The recent observation of correlated states of bosonic Theorganizationisasfollows. InsectionII,wepresent a atoms in optical lattices has generated much interest.1,2 the general setting for the study of spin order-disorder J quantum phase transitions. In section III, we present As known for a while, when bosons in lattices interact 5 mean field results on the quantum phase transitions in with each other repulsively, they can be localized and form a Mott insulating state instead of a condensate.3,4 both small N and large N limits. In section IV, we dis- 2 cuss issues which are to be understood in the future. v Thisphenomenonhasbeenobservedintheopticallattice 8 experiment. By varying laser intensities of optical lat- 9 tices, Greiner et al. have successfully investigated Mott 1 states of spinless bosons by probing spin polarized cold II. ALGEBRA AND SETTING 6 atoms in optical lattices with a large potential depth.1,2 0 A. The microscopic Hamiltonian in the dilute limit We are interested in spin correlated Mott insulating 3 statesofspin-onebosons,especiallyspin-onebosonswith 0 / antiferromagneticinteractions. Someaspectsofspincor- The microscopic lattice Hamiltonian we employ to t a relatedMottinsulatingstates wereinvestigatedrecently. study spin correlated states of spin-one bosons is: m For an even number of particles per site, both spin sin- - glet Mott insulators and nematic Mott insulators were Hmicroscopic = −t (ψk†,mψl,m+h.c.) nd found in certain parameter regimes, while for high di- Xhkli o msiteensoinonlyalnleamttaicteics winistuhlaatningodsdtantuesmbweerreofprpoaprtoisceldes.5pIenr + ψk†,mψk,mUρ(k,l)ψl†,m′ψl,m′ (1) c k,l : one-dimensionallattices,itwasdemonstratedthatforan X iv odd number of particles per site, Mott states should be + ψk†,mSmγnψknUS(k,l)ψl†,m′Smγ′n′ψl,n′ X dimerized valence-bond-crystals,which support interest- k,l X ingfractionalizedquasi-excitations.6Effectsofspincorre- r a lationsonMottinsulator-superfluidtransitionshavebeen Here ψk†,m is the creation operator of a spin-one particle studied and remain to be fully understood.7 at site k with spin-index m = 0, 1. kl indicates that ± h i Inthisarticle,weanalyzethemicroscopicstructuresof the sum should be taken over nearest neighbors and Sγ spin singlet Mott insulating states (SSMI) and nematic (γ =x,y,z) are spin-one matrix operators given as: Mott insulating states (NMI). We study, quantitatively, 0 1 0 0 i 0 1 0 0 quantum phase transitions between these two phases in Sx= 1 1 0 1 Sy= 1 i −0 i Sz= 0 0 0 . high-dimensionalbipartitelattices. Inthemeanfieldap- √2 √2 − 0 1 0 0 i 0 0 0 1 proximation, we demonstrate that for an even number − of particles per site, the transitions are weakly first or- Uρ(k,l) and Us(k,l) are, respectively, spin-independent der. We should emphasize that results obtained in this and spin-dependent interaction parameters between two paper are only valid in high dimensions. In one dimen- bosons at site k and l. sional lattices, nematic order does not survive long wave In the dilute limit, which is defined as a limit where length quantum fluctuations; detailed discussions onlow ρ¯a3 1(aisthescatteringlengthandρ¯theaverageden- ≪ sity),atomsscatterins-wavechannels. Fortwospin-one atoms,thescatteringtakesplaceinthetotalspinS =0,2 ∗Permanentaddress: DepartmentofPhysicsandAstronomy, channels, with scattering lengths a0,2. Interactions be- tweenatomscanbeapproximatedasspin-dependentcon- UniversityofBritishColumbia,6224AgricultureRoad,Vancouver, B.C.V6T1Z1,Canada tact interactions.8 In the lattice model introduced here, 2 calculations yield Of particular interest is the singlet creation operator Uρ(k,l)=E δ and US(k,l)=E δ . (2) 1 1 c kl s kl √6ψk†,αψk†,α = √6(ψk†,0ψk†,0−2ψk†,1ψk†,−1). (10) The parameters E and E are given by: c s We find the following properties for this operator: 4π~2ρ¯(2a +a ) 4π~2ρ¯(a a ) 2 0 2 0 Ec = 3MN c˜, Es = 3MN− c˜, (3) [Sˆαk,ψl†,αψl†,α] = [Sˆαk,ψl,αψl,α]=0; (11a) where N is the average number of atoms per site, M is [ψk,αψk,α,ψl†,β] = 2δklψk,β; (11b) the mass of atoms and c˜is a constant. [ψk,αψk,α,ψl†,βψl†,β] = δkl(4ρˆk+6). (11c) B. Algebras C. The on-site dynamics For the study of spin correlatedstates in lattices, it is The total spin operator can be expressed as: rather convenient to introduce the following operators: 1 Sˆ2k =ρˆk(ρˆk+1)−ψk†,αψk†,αψk,βψk,β. (12) ψk†,x = √2(ψk†,−1−ψk†,1) (4a) So, eigenstates of the total spin operator have to i be eigenstates of the ”singlet counting operator” ψk†,y = √2(ψk†,−1+ψk†,1) (4b) ψk†,αψk†,αψk,βψk,β. Defining the state Ψn such that: ψk†,z = ψk†,0 (4c) k,0 wherek againlabelsalatticesite. Inthisrepresentation: ρˆkΨnk,0 =nΨnk,0, ψk†,αψk†,αψk,βψk,βΨnk,0 =0, (13) we find that wave functions of these eigenstates are: ψk†,mSmαnψk,n =Sˆαk ≡−iǫαβγψk†,βψk,γ (5) Ψnk,m =C(ψk†,αψk†,α)mΨnk,−02m (14) α,β,γ x,y,z . Thedensityoperatorcanbeexpressed ∈{ } where C is a normalization constant. From Eq. 11c it in a usual way: follows that: ρˆk ≡ψk†,γψk,γ. (6) ψk†,αψk†,αψk,βψk,βΨnk,m =(4m(n−m)+2m)Ψnk,m (15) Consequently the Hamiltonian is given as: Using that ρˆ Ψn =nΨn we derive: k k,m k,m Hlatt=−t˜Xhkli(ψk†,αψl,α+h.c.)+Xk Ecρˆ2k+EsSˆ2k−ρˆkµ(7) Sˆ2kΨnk,m =(n−2m)(n−2m+1)Ψnk,m. (16) So S = n 2m. Now if n is even, S is also even and k k where we have introduced the chemical potential µ; Sˆ2 whenniso−dd,S isoddtoo. Foranevennumberofpar- k k is the total spin operator Sˆ Sˆ . ticles per site N the states labeled by S =0,2,4,...,N k,α k,α k ψ (α = x,y,z) are bosonic operators obeying the are present, whereas for an odd number of particles per k,α following commutation relations: site S =1,3,5,...,N are allowed. This reflects the ba- k sic property of the many body wavefunction of spin-one [ψk,α,ψl,β]=[ψk†,α,ψl†,β]=0, [ψk,α,ψl†,β]=δklδαβ. (8) bosons,whichhastobesymmetricundertheinterchange of two particles. Taking into account Eqs. 5,6,8, one can verify the fol- Solutions for spin correlated condensates with finite lowing algebras: numbers of particles were previously obtained9; in the thermodynamical limit, these states evolve into polar [Sˆαk,ψl,β] = δkliǫαβγψk,γ (9a) condensates.8,10,11 Also there,two-bodyscatteringswere [Sˆαk,ψl†,β] = δkliǫαβγψk†,γ (9b) showntoleadtoeither”antiferromagnetic”or”ferromag- netic” spin correlations in condensates. Spin correlated [Sˆαk,Sˆβl] = δkliǫαβγSˆγk (9c) condensates have been investigated in experiments.12,13 [ρˆ ,ψ ] = δ ψ (9d) k l,α kl k,α − D. The effective Hamiltonian for Mott states [ρˆk,ψl†,α] = δklψk†,α (9e) In the limit when t˜ E , atoms are localized and c [Sˆα,ρˆ] = 0 (9f) onlyvirtualexchangepr≪ocessesareallowed. Aneffective k l 3 Hamiltonianinthislimitcanbederivedinasecondorder III. PHASE TRANSITIONS BETWEEN SSMI’S perturbative calculation of the Hamiltonian in Eq.7: AND NMI’S HMott = S2ˆI2k −J˜ex ψk†,αψl,αψl†,βψk,β +h.c. . A. Two particles per site Xk Xhkli(cid:16) (cid:17)(17) Inthecaseoftwoparticlespersite,theon-siteHilbert Here J˜ = t˜2 . In deriving Eq. 17, we have taken into space is six-dimensional, including a spin singlet state ex 2Ec account that E E . s c ≪ Tofacilitatediscussions,weintroducethefollowingop- ψ ψ erator: S =0,S =0 = η† η† 0 , (23) z | i √6 | i 1 Qˆk,αβ =ψk†,αψk,β − 3δαβψk†,γψk,γ, (18) and five spin S =2-states whose expectation value √3 Qˆ |Qηξi= 2 Qηξψη†ψξ†|0i (24) Q˜ = h αβi (19) Qˆ where Q is a symmetric and traceless tensor with five αβ ref ηξ h i independentelements. AllstatesintheHilbertspaceare is the nematic order parameter. The reference state symmetric under the interchangeofbosons;as expected, ψ = (nαψk†,α)N 0 is a maximally ordered state. the states Qηξ are orthogonal to S = 0,Sz = 0 . It is ref k √N! | i convenient|tochioosethe followingr|epresentationoifQ : Choosing n=e , we obtain: ηξ Q z 1 0 0 hQˆαβiref =N−03 −31 0 . (20) Qηξ(n)=nηnξ− 13δηξ, (25) 0 0 2 3 with the director n as a unit vector living on S2. States Q˜ varies in a range of [ 1,1]. defined by the director n form an over-complete set in −2 IntermsoftheoperatorQˆ ,theeffectiveMottHamil- the subspace spanned by five S =2 states. αβ tonian can be rewritten as (up to an energy shift): Whenthe hopping iszero,onenoticesthatthe Hamil- tonianinEq.17commuteswithSˆ2;thegroundstatewave k Heff =Es Tr[QˆkQˆk−QˆkQˆ†k]−J˜ex Tr[QˆkQˆl]. (21) function is Xk Xhkli Ψ = S =0,S =0 . (26) z k | i | i Finally we define k Y η˜= zJ˜ex (22) nOinanthcoemomthuetrehsawnidth, wTrh[eQˆn EsQˆgoes]taonzdertoh,etghreouHnadmsitlatote- E k,αβ l,βα s wave function can be confirmed as: as a dimensionless parameter, which can be varied con- tinuously; z is the coordination number of lattice. 2 1 Ψ = Q(n) + S =0,S =0 . (27) k z k | i 3| i √3| i k r Y E. The range of the physical parameters for any choice of the director n. To study spin nematic or spin singlet Mott states at From Eq. 3 it is clear that E and E depend on the an arbitrary η˜, we introduce a trial wave function which s c density, number of atoms, the mass of atoms and scat- is a linear superposition of singlet states and symmetry tering lengths. However,their ratio depends only on the breaking states: scattering lengths. According to current estimates14,15, for sodium atoms this ratio is given as Es 9 10 2. In Ec ≈ · − thiTshpeappaerra,mweetaerret˜icnatnerbeestveadriiendtihnedelpimenitdEenstl≪ybEycc.hang- |Ψiθ = cosθ|S =0,Sz =0ik+sinθ|Qηξ(n)ik. (28) k Y ingthedepthoftheopticallattice. Awiderangeisexper- imentallyaccessible;onecanvaryfromtheregimewhere Here θ is a variable to be determined by the variational t˜ E to a regime where t˜ E . We limit ourselves to method. c s M≫ottstates(t˜ E ), where≪allbosonsarelocalized,but Astraightforwardcalculationleadstothefollowingre- c ≪ the ratio η˜can have arbitrary values. sults: 4 E/Es Q˜ 0.015 0.8 0.01 0.005 0.6 Q˜ 0.4 0.2 0.4 0.6 0.8 -0.005 0.2 -0.01 η˜ 0.5 1 1.5 2 -0.015 -0.2 FIG.1: Energy(measuredinunitsofEs)versusQ˜forvarious η˜ for N = 2. Curves from top to bottom are for η˜ = 0.97, 0.99, 1.0, 1.02. FIG. 2: (Color online) The nematic order parameter as a function of η˜for N =2. The phase transition takes place at η˜ = 1. Data along the black lines represent ground states; the red (light) lines are for metastable states. Spheres with double-headedarrowsareintroducedtorepresentorderingin director n defined in Eq. 25 in different Mott states. In spin singlet states, the director n is uncorrelated; in rod-like ne- E(θ) = Ψ Ψ (29) h |H| iθ maticstates,thedirector nisordered andindisk-likestates, 2 = 6E sin2θ zJ˜ ((2√2cosθsinθ+sin2θ)2 theaxis of theeasy plane of the director nis ordered. s ex − 3 sin2θ Q˜ = (√2cosθsinθ+ ) (30) 2 It is worth emphasizing that a positive Q˜ corresponds In terms of Q˜, the energy can be expressed as: to a rod-like nematic state; for Q˜ =1 the state is micro- scopically given by: 4 2 4 8 E =6E + Q˜ 2Q˜2+Q˜+1 zJ˜ Q˜2, s ex 9 9 − 9 − − 3 (cid:18) q (cid:19) (nαψα†)2 0 . (32) which for Q˜ 1 can be expanded as: √2 | i ≪ 3E 8zJ˜ Q˜2 3E Q˜3+ 39E Q˜4+... (31) A solution with negative Q˜ indicates a disk-like ne- s ex s s − 3 − 2 16 matic state; the microscopic wave function is (cid:18) (cid:19) The cubic term leads to a first order phase transition 1 (n ψ )2 isnituthateiomneianncfilaesldsicaaplpnreomxiamtiactiloiqnu,iwdhcircyhstiaslssi.m17ilar to the 2ψη†ψη† − α2α† |0i. (33) (cid:18) (cid:19) In FIG.1 the Q˜-dependence of energy is plotted for variousη˜inthevicinityofaquantumcriticalpoint(mean at Q˜ = 1. For n = e the wave functions in field). For η˜< 0.985, the energy has only one minimum Eq.32,33be−co2me 1 ψ ψ 0 anzd 1 ψ ψ +ψ ψ 0 re- at Q˜ =0 and correspondingly the ground state is a spin √2 z† z†| i 2 x† x† y† y† | i spectively. singlet Mott state. When 0.985 < η˜ < 1.0, in addition (cid:0) (cid:1) to the global minimum at Q˜ = 0, there appears a local We have also tried a five-parameter variational ap- minimum at Q˜ > 0, which represents a spin nematic proach,takingintoaccountthefullon-siteHilbertspace. metastable state. When η˜>1.0 the solution with Q˜ >0 In a slightly different representation we write the trial becomes a global minimum and the solution at Q˜ =0 is wave function as: metastable; consequently the ground state is a nematic Mott state. For η˜ > 9, the solution at Q˜ = 0 becomes Ψ = (c xx +c yy +c zz 8 | i xx| ik yy| ik zz| ik unstable ; but an additional local minimum appears at k Q˜ <0 whichweinterpretas anew metastablestate (not Y+c xy +c xz +c yz ) (34) xy k xz k yz k shown in FIG. 1). | i | i | i The evolution of ground states as η˜ is varied, is sum- marized in FIG. 2. As is clearly visible, the phase tran- Here |ααik = √12ψk†,αψk†,α|0i (no summation) and psihtiaosne-itsraanwsiteiaoknly(η˜fir=st1o.0r)deisreoqnuea.lTtohe1j.ump in Q˜ at the |pαrβesiskio=n fψork†,tαhψek†e,βn|e0rig.y:This results in the following ex- 2 5 The Hamiltonian in Eq. 7 can be mapped to a Con- strained Quantum Rotor Model (CQR), describing the E = Es 4(c2xx+c2yy+c2zz −cxxcyy−cxxczz−cyyczz) dynamics of two unit vectors (n,eiχ) on a two-sphere +6(cid:2)(c2xy+c2xz +c2yz) and a unit circle: −zJ˜ex 6(c4xx+c4yy+(cid:3)c4zz)+4(c4xy+c4xz+c4yz) +4(c2xx(cid:2)c2yy+c2xxc2zz +c2yyc2zz) = t n n cos(χ χ )+ E Sˆ2+E ρˆ2 ρˆ µ +12(c2 c2 +c2 c2 +c2 c2 +c2 c2 HCQR − k· l k− l s k c k− k xx xy xx xz yy xy yy yz Xhkli Xk +c2 c2 +c2 c2 ) (41) zz xz zz yz t = Nt˜. The CQR-model has been introduced to study +8(c2xyc2xz+c2xyc2yz+c2xzc2yz) spin-one bosons in a few previous works and we refer +8√2(c +c +c )c c c to those papers for detailed discussions.5,6,16 For Mott xx yy zz xy xz yz states the effective Hamiltonian can be found as: +4(c2 c2 +c2 c2 +c2 c2 ) xx yz yy xz zz xy +8(c c c2 +c c c2 +c c c2 ) . (35) xx yy xy xx zz xz yy zz xz The conclusions are almost the same and s(cid:3)ummarized =E S2 J (n n)2; J = t2 , (42) below: H s k− ex k· l ex 2Ec i) For η˜ < 0.985. the only minimum is at c = c = Xk Xhkli xx yy c = 1 , c =0 for α=β. zz √3 αβ 6 ii) At η˜=0.985 additional local minima appear. iii) At η˜=1 a first order phase transition takes place. and we define η = zJex. Es iv) For 8 > η˜>1, the global minimum is at Q˜ >0, but 9 Ingeneral,wechoosethe on-sitetrialwavefunctionto the Q˜ =0-solution remains to be a local minimum. be: v) At η˜ = 9 the solution at c = c = c = 1 8 xx yy zz √3 becomes unstable. vi) However,the disk-like Q˜ <0-solutionappears in this σ case as a saddle point. ψ(n )=C exp (n n )2 . (43) k σ k 0 2 · h i B. Large N limit: An even number of particles per site C is a normalization constant. When σ 0 this yields σ → an isotropic state Y (n ), which indicates a spin singlet 00 k Foralargenumberofparticlespersite,itisconvenient state. When σ + , n is localized on the two-sphere k → ∞ to introduce the following coherent state representation: inthevicinityofn ,representingarod-likenematicstate 0 and when σ , n lies in a plane perpendicular k → −∞ to n , corresponding to a disk-like spin nematic state. n,χ = 1 N+δN exp( imχ) nαψα† m 0 (36) More0over this wave function has the following property: | i √2δNm=XN−δN − (cid:0)2(m−(cid:1)1)!| i ψpa(r−tinckle)s=peψr(snitke).1,6as is required for an even number of p where the director n is again a unit vector on S2 given Choosing n =e this gives: by (cosφsinθ,sinφsinθ,cosθ). In this representation 0 z ∂ ρˆ = i (37) ∂χ σ k ψ(φ ,θ )=C exp cos2θ . (44) k k σ k ∂ 2 Sˆ = in (38) h i × ∂n 1 ∂ ∂ 1 ∂2 Sˆ2 = sinθ + (39) − sinθ∂θ ∂θ sin2θ∂φ2 (cid:20) (cid:18) (cid:19) (cid:21) The expectationvalue ofthe Hamiltonianinthis state 1 Qˆ = N n n δ (40) is: αβ α β αβ − 3 (cid:18) (cid:19) 6 3 1 3eσ σ 12e2σσ 4eσ√π σ (3+2σ)Erfi σ +π(3+4(σ+σ2))Erfi2 σ E =E σ+ | | zJ − | | | | | | σ s −4 − 2 2√πEprfi σ!− ex p 8πσ2Erpfi2 σ p ! | | | | p p Q˜ E/Es 0.015 0.6 0.01 0.4 0.005 σ 0.5 1 1.5 2 2.5 3 0.2 -0.005 η -0.01 5 10 15 20 -0.015 -0.2 FIG.3: Energy(inunitsofEs)asafunctionofσ forvarious FIG.4: (Coloronline)Nematicorderparameterasafunction η (N =2k≫1). From top to bottom are curvesfor η=9.9, ofηforN =2k(≫1). Alongtheblacklinesaregroundstates; 9.96, 10, 10.0965, 10.3. alongthered(light)linesaremetastablestates. (Seealsothe caption of FIG. 2.) in which Erfi[x] is the complex error function defined by parameter can be calculated as: Erf[ix]/i. In a series expansion for σ 1, the result is: ≪ σ Qˆ σ Q˜ = h | αβ| i zJ 2 8 Qˆ − 3ex + 15Es− 675zJex σ2 h∞1| αβ3|∞i 3eσ (cid:18) (cid:19) = + (46) 4 32 −2 − 4σ 2 π σ Erfi σ + Es zJex σ3 | | | | 315 − 14175 (cid:18) (cid:19) When Q˜ is small, we obtain apn expressipon of energy in 8 32 + Es+ zJex σ4+o(σ5) (45) terms of Q˜: −4725 165375 (cid:18) (cid:19) zJ 15 2 75 1275 E = ex+ E zJ Q˜2 E Q˜3+ E Q˜4 The energy as a function of σ at different η is plot- Q˜ − 3 2 s− 3 ex −14 s 98 s (cid:18) (cid:19) ted in FIG. 3, which is qualitatively the same as FIG. (47) 1 for two particles per site. When η < 9.96, the en- The jump in Q˜ at the phase transition is equal to 0.323. ergy as a function of σ has only one (global) minimum, The evolution of ground state wave functions and re- which corresponds to a spin singlet ground state. When sults on quantum phase transitions are summarized in η > 9.96, in addition to the global minimum, there ap- figure4,wherethenematicorderparameterisplottedas pears a local minimum at σ > 0. At η = η = 10.0965, c a function of η. As stated before, these results are only these two minima become degenerate,signifying a phase valid in high dimensional lattices, where fluctuations in transition. At η > η , the solution at σ = 0 becomes a c orderedstatesaresmall. Fordetailedcalculationsoffluc- localminimumindicatingametastablespinsingletstate, tuations we refer to appendix B. whereas the global minimum at σ > 0 corresponds to a nematic ground state. As η further increases, the solu- tion at σ = 0 becomes unstable and a local minimum C. Large N limit: An odd number of particles per occurs at σ < 0, while the global minimum remains at site σ > 0. Following discussions on Eqs. 32,33 we inter- pret the σ < 0 solution as a metastable disk-like spin At last, we also present results for an odd number of nematics. atomspersite. Themaindifferencebetweenthiscaseand ForthetrialwavefunctioninEq. 44thenematicorder thecaseforanevennumberofparticlespersiteisthatat 7 zero hopping limit in the former case there is always an O unpaired atom at each site. Consequently in the mean 17.5 field approximation, we only find nematic Mott insulat- 15 ing phases. As in the case for even numbers of particles persite,weexpectthisapproximationtobevalidinhigh 12.5 dimensionallatticesbutfailinlowdimensions,especially 10 in one-dimensional lattices where long wave length fluc- 7.5 tuations are substantial. Here we restrict ourselves to high dimensional lattices only. 5 ForlargeN atrialwavefunctionwhichinterpolatesbe- 2.5 tween spin singlet states (dimerized) and nematic states η can be introduced as: 5 10 15 20 Ψ ( n ) = C(O,σ)[O(n n )(n n )+(n n )] odd { k} k· 0 l· 0 k· l FIG. 5: The valueof O as a function of η (N =2k+1≫1). hYklip exp[σ((n n )2+(n n )2)]. (48) k 0 l 0 × · · σ kl denotes that the summation should be taken over p h i parallelyorderedpairsofnearestneighborskandlcover- ingthelattice. C(O,σ) isanormalizationconstant. The 8 solution with O = 0,σ = 0 corresponds to a dimerized valence bond crystal state; and solutions with O =0, or 6 6 σ =0 represent nematic states. 6 It is straightforward, but tedious to compute the en- 4 ergyofthesestates. Minimizingitwithrespecttovarious values of η for d = 3 gives the results shown in FIG. 5 2 and 6. No phase transitions are found in the mean field approximation; and ground states break both rotational η and translational symmetries.18 5 10 15 20 At very small η, the on-site Hilbert space is truncated intotheoneforaspin-oneparticle.6 ThereducedHamil- FIG. 6: The valueof σ as a function of η (N =2k+1≫1). tonian in the truncated space is a Bilinear-Biquadratic model for spin-1 lattices Thus, we expect that fluctuations play a very important =J [cosθS S +sinθ(S S )2],S2 =2; (49) Hb.b. k· l k· l k role in these transitions and the full theory on quantum Xhkli phase transitions remains to be discovered. On the other hand, we have estimated fluctuations in θ in general varies between 3π/4 and π/2. We there- − − different regimes of the parameter space. We found that for expect ground states at small η limit should still ex- fluctuationsareindeedsmallawayfromthecriticalpoint, hibit nematic order (i.e. O =0). 6 ateithersmallhoppingorlargehoppinglimitforaneven It is worth emphasizing that conclusions about small number of particles per site. At the small hopping limit, η limit arrived here are only valid in high dimensional fluctuationsareproportionaltoη,whileatthelargehop- bipartite lattices. In low dimensional lattices, states of pinglimittheycanbeestimatedtobeproportionalto 1 correlatedatomsinthislimitwerediscussedrecentlyand √η ground states could be rotationally invariant dimerized- (see Appendix B). valence-bond crystals.6 For an odd number of particles per site, fluctuations are small only at large hopping limit and are significant at small hopping limit. The later fact implies a large IV. CONCLUSIONS degeneracy of Mott states at zero hopping limit which was emphasized in the discussions on low dimensional Wehavestudiedthemicroscopicwavefunctionsofspin Mottstates. Thephysicsinthislimitremainstobefully nematic and spinsingletMott states. Bothdisk-like and understood. rod-like spin nematic states were investigated. We also Inthecontextofantiferromagnets,spinnematicstates have analyzed quantum phase transitions between spin have also been proposed.19,20,21 Collective excitations in singlet Mott insulating states and nematic Mott insulat- atomic nematic states should be similar to those studied ing states. We show that in the mean field approxima- in previous works; we present some brief discussions on tion, the phase transitions are weakly first order ones. thissubjectinAppendixBandreferto19,20,21 fordetails. 8 V. ACKNOWLEDGMENT APPENDIX B: SPIN FLUCTUATIONS IN MOTT STATES ThisworkissupportedbythefoundationFOM,inthe Nonlinear dynamics and spin fluctuations in conden- Netherlands under contract 00CCSPP10, 02SIC25 and sates of spin-one bosons were discussed in a previous NWO-MK projectruimte 00PR1929;MS is alsopartially work.16 Here we carry out a similar discussion for Mott supported by a grant from Utrecht University. states. Following the Hamiltonian Note added: At a later stage of this work, we =E Sˆ2 J (n n)2 (B1) received a copy of the manuscript by A. Imambekov, M. H s k− ex k· l Lukin and E. Demler where similar results have been Xk Xhkli obtained. we derive the following equation of motion for the direc- tor n : k APPENDIX A: AN ALTERNATIVE dnk =2E Sˆ n (B2) DESCRIPTION dt s k× k Fluctuations when η is small Alternatively, one can also carry out the calculations For η =zJ /E = 0 and an even number of particles ex s in section III, using the following operator: per site, the ground state is the product state: Qˆ = SˆαSˆα′ 1δ SˆγSˆγ Ψη=0 = Y00(nk) (B3) 2,αα′ αα′ − 3 k Y = −ǫαβγǫα′β′γ′ψβ†ψβ†′ψγψγ′ +δαα′ψη†ψη When 0 < η ≪ 1, the ground state wave function can 1 also be obtained by a perturbation theory; the leading −ψα†′ψα− 3δαα′Sˆ2tot (A1) term is Defining the more conventionalorder parameter19 Ψ(0) J (n n )2 Ψ(0) Ψ(1) = h lm|− ex hkli k· l | 00i (B4) Q˜ = hQˆ2,αα′i , (A2) 0 l6=X0,m EP0(00)−El(m0) 2 Qˆ h 2,αα′iref InourcaseΨ(lm0) =Ylm withleven,andEl(0) =l(l+1)Es. A direct calculation yields weobtainthefollowingresultsforthetrialwavefunction in Eq. 28: 2 η Ψ(1) = Y (n )Y (n ) Y (n ). Q˜2 = 32sin2θ 0 45z XhijimX=−2 2m i 2,−m j kY6=i,j 00 k r (B5) E = 6Es 23Q˜2 resTualtksiningitnhtiosalicmcoitu,nthY00|Qˆαβ|Y00i=0,wefinddesired r 2 2zJ˜ 2√2 2(1 2Q˜ )Q˜ + 2Q˜ hQˆk,αβi=0. (B6) ex 2 2 2 −3 s 3 − 3 3 r r r To characterize fluctuations, we study the following correlation function Qˆk,ααQˆk′,αα . Calculations of this 16 2 h i = 2√6E zJ˜ Q˜ correlationfunction in the state given in Eq. B5 yield s ex 2 − 3 r3 ! 2η 16√42zJ Q˜3/2+ 28zJ˜ Q˜2+O(Q˜5/2) hQˆk,ααQˆk′,ααi= 45δ(kk′,hkli)× − √437 ex 2 9 ex 2 2 2 Y Qˆ Y Y Qˆ Y +h.c. . 2m αα 00 2, m αα 00 whichleadtothesameconclusionsasinsectionIII.How- h | | ih − | | i ever, in terms of the order parameter defined in Eq. A2, mX=−2(cid:16) (cid:17) the rod-like and disk-like structures shown in FIG.2 and δ(kk′, kl ) is unity if k′ and k sites are two neighboring h i FIG.4 are less obvious. sites as kl and otherwise is zero. The last expression h i can be calculated explicitly, In the case of a large number of particles per site, the order parameter introduced here has the same expecta- 8 tion value as the operator in Eq. 40. Y2m Qˆαα Y00 Y2, m Qˆαα Y00 +h.c. = (B7) h | | ih − | | i 45 (cid:16) (cid:17) 9 Clearly at small η, fluctuations are small. To obtain eigenmodes, we perform a Fourier transfor- Fluctuations when η is large mation (setting the lattice spacing to be unity), Againwe considerthe casefor anevennumber of par- ticles per site. In the limit of η , all directors n k → ∞ point in the direction of e . For a finite but large η we z 1 1 introduce Πˆk,α = Πˆq,αeiπk·q,Ck,α = Cq,αeiπk·q, √V √V T q T q X X (B13) nk =ez 1−Ck2x−Ck2y +Ckxex+Ckyey (B8) twohtehree VfoTlloiswtinhge Htoatmalilntuomnibaner of lattice sites. This leads q where C , α=x,y are much less than unity. kα Following discussions in section IIIB, we obtain the qπ = [E Πˆ2 +zJ sin2 | | C2 ]. (B14) following commutators, H s q,α ex 2 q,α q,α X [Sˆ ,C ] iδ ,[Sˆ ,C ] iδ (B9) ky k′x k,k′ kx k′y k,k′ ≈ ≈− Following a standard calculation, fluctuations in this whichdefinetwosetsofharmonicoscillators. Introducing limit are: Πˆ =Sˆ ,Πˆ = Sˆ (B10) y x x y − 1 2 1 1 C2 = C 2 = . the effective Hamiltonian becomes hXα k,αi VThXq,α| q,α| i √ηVT |qX|<qc sin|q|π/2 = [E Πˆ2 +J (C C )2]. (B11) (B15) H s k,α ex k,α− l,α Themomentumcut-offq ingeneraldependsontheshort αX=x,y Xk Xhkli distance behavior of ourcmodel and for simplicity we set And it as one. In high dimensional lattices, the sum in Eq. B15 is convergent; and we see the fluctuations are also [Πˆ ,C ]=iδ δ . (B12) small at the large η limit. kα k′β k,k′ αβ 1 M.Greiner,O.Mandel,T.Esslinger, T.W.H¨anschandI. 14 J. P. Burke Jr., C.H. Greene and J.L. Bohn, Phys. Rev. Bloch, Nature415, 39 (2002). Lett. 81, 3355 (1998). 2 M.Greiner,O.Mandel,T.W.H¨anschandI.Bloch,Nature 15 C.V. Ciobanu, S.K. Yip and T. L. Ho, Phys. Rev.A. 61, 419, 51 (2002). 033607 (2002). 3 M. P. A Fisher, P. B. Weichman, G. Grinstein and D.S. 16 FeiZhou,Phys.Rev.Lett.87, 080401-1 (2001); FeiZhou, Fisher, Phys.Rev.B. 40, 546 (1989). Int.Jour. Mod. Phys. B 17, 2643 (2003). 4 D. Jaksch, C. Bruder, J. I. Cirac, C.W. Gardiner and P. 17 P.G.deGennes,Thephysicsofliquidcrystals,OxfordUni- Zoller, Phys. Rev.Lett. 81, 3108 (1998). versity Press, 1974. 5 E. Demler and F. Zhou, Phys. Rev. Lett. 88, 163001-1 18 At thetime of submission, we noticed that nematicstates (2002) for one particle per site were also addressed very recently 6 F. Zhou, Europhys. Lett 63, 505 (2003); F. Zhou and M. in S.K. Yip, cond-mat/0306018. The 3D results there are Snoek, Annalsof Physics 308, 692 (2003). consistentwiththemeanfieldresultsinSectionIIIC.How- 7 S. Tsuchiya, S. Kurihara and T. Kimura, ever, most one-dimensional states proposed in that paper cond-mat/0209676 breaktherotationalsymmetryandareunlikelytobecan- 8 T. L. Ho,Phys. Rev.Lett. 81, 742 (1998). didates for true ground states. Some general discussions 9 C.K. Law, H. Pu and N.P. Bigelow, Phys. Rev. Lett. 81, on one-dimensional Mott states of spin-one bosons can be 5257 (1998). found in6. 10 T. Ohmiand K.Machida, Journal of thePhysical Society 19 A.F. Andreev and I.A. Grischuk, Zh. Exp. Teor. Fiz. 87, of Japan 76, 1822 (1998). 467 (1984)[Sov. Phys.JETP 60(2), 267 (1984)]. 11 Y. Castin and C. Herzog, cond-mat/0012040 20 A.V.Chubukov,J.Phys.Cond.Mat.2,1593(1990);Phys. 12 D.M. Stamper-Kurn, M. R. Andrews, A.P. Chikkatur, S. Rev.B 43, 3337 (1991). Inouye, H.-J. Miesner, J. Stenger and W. Ketterle, Phys. 21 P. Chandra, P. Coleman and A. I. Larkin, J. Phys. Cond. Rev.Lett. 80, 2027 (1998). Mat. 2, 7933 (1990). 13 J. Stenger,S.Inouye,D.M.Stamper-Kurn,H.-J.Miesner, A.P. Chikkaturand W. Ketterle, Nature 396, 345 (1998).