Bogoliubov theory of interacting bosons on a lattice in a synthetic magnetic field
Stephen Powell, Ryan Barnett, Rajdeep Sensarma, Sankar Das Sarma
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] S e p Bogoliubov theory of interacting bosons on a lattice in a synthetic magnetic field
Stephen Powell, Ryan Barnett, Rajdeep Sensarma, and Sankar Das Sarma
Joint Quantum Institute and Condensed Matter Theory Center,Department of Physics, University of Maryland, College Park, MD 20742, USA
We consider theoretically the problem of an artificial gauge potential applied to a cold atomicsystem of interacting neutral bosons in a tight-binding optical lattice. Using the Bose-Hubbardmodel, we show that an effective magnetic field leads to superfluid phases with simultaneous spatialorder, which we analyze using Bogoliubov theory. This gives a consistent expansion in terms ofquantum and thermal fluctuations, in which the lowest order gives a Gross-Pitaevskii equationdetermining the condensate configuration. We apply an analysis based on the magnetic symmetrygroup to show how the spatial structure of this configuration depends on commensuration betweenthe magnetic field and the lattice. Higher orders describe the quasiparticle excitations, whosespectrum combines the intricacy the Hofstadter butterfly with the characteristic features of thesuperfluid phase. We use the depletion of the condensate to determine the range of validity of ourapproximations and also to find an estimate for the onset of the Mott insulator phase. Our theoryprovides concrete experimental predictions, for both time-of-flight imagery and Bragg spectroscopy.
PACS numbers:
I. INTRODUCTION
One of the main goals of experiments with cold atomshas been to extend to new contexts physical effects thatare familiar from the study of condensed matter. Fromthis point of view, a particularly important area has beenthe study of vortices in superfluid systems, previouslyconsidered in the context of type-II superconductors andsuperfluid helium. When a trapped superfluid comprisedof weakly interacting bosonic atoms is caused to rotate,it forms a lattice of vortices, a remarkably clear demon-stration of the quantization of circulation.
In a corotating reference frame, the (neutral) atoms ex-perience a Coriolis force of the same form as the Lorentzforce on charged particles in a magnetic field. This mag-netic analogy has inspired considerable experimental ef-fort to achieve effective magnetic fields large enough toreach the quantum Hall regime, and corresponding the-oretical work to extend the theory of the quantum Halleffect to the context of trapped bosons. Much recent effort has been directed towards combin-ing effective magnetic fields and optical lattices, both toincrease the stability of the experimental systems, andbecause of interesting physical effects expected in thepresence of a lattice. The first experiments with rotatinglattices used masks to produce parallel beams, whose sub-sequent interference formed the optical lattice potential.Mechanical rotation of the mask caused the interferencepattern to rotate, and a density of vortices comparable tothe density of lattice sites was achieved. Mechanical in-stabilities limited the lattice strength, however, restrict-ing to the regime where the vortices are weakly pinned. More recent experiments have replaced the mask withan acousto-optic modulator, allowing for considerablydeeper lattices and lower temperatures. For a sufficientlydeep lattice, a single-band approximation becomes valid,and the Bose-Hubbard model gives a good descriptionof the physics. Experiments with Raman lasers allow for effective magnetic fields without rotation, instead us-ing coupling between internal atomic states to imprintthe required geometric phases. A static optical latticecan be applied to such an arrangement, giving an effec-tive magnetic field within the laboratory frame. Variousother proposals have been made to induce phases directlywithin the lattice, some of which have the advantageof producing a perfectly commensurate flux density.A major goal of these experiments is, as noted above,to reach the quantum Hall regime, in which the ef-fective flux density (in units of the flux quantum φ =2 π ~ /Q where Q is the effective particle charge) is com-parable to particle density. Continuum systems in thisregime exhibit a sequence of incompressible phases, theinteger and fractional quantum Hall states, and, in thepresence of a lattice, closely related physics has been pre-dicted in certain areas of the phase diagram. Systems of bosons also support compressible superflu-ids, which are more closely related to the phases in theabsence of a magnetic field, and are likely competitorswith the quantum Hall states in the phase diagram. De-tailed study of these phases is beneficial for the searchfor quantum Hall physics, both to calibrate experimentsand to understand this competition, but their nontrivialproperties and phenomena, especially in comparison withconventional superfluids, make them worthy of consider-able interest in their own right.These properties are inherited from the remarkablestructure of the corresponding noninteracting problem,a single particle moving on a tight-binding lattice in thepresence of a uniform magnetic field.
The densityof states exhibits a fractal structure known as the ‘Hofs-tadter butterfly’ (see Section II C, and in particular Fig-ure 2), with the spectrum depending sensitively on α ,the magnetic flux per plaquette of the lattice (measuredin units of φ ). For rational α = p/q (with p and q co-prime), there are q bands, and each state is q -fold degen-erate. This degeneracy can be understood in terms of the‘magnetic symmetry group’, which takes into accountthe modification of the symmetries inherent in making aparticular choice for the gauge potential.The present work introduces a theoretical approachthat starts with the noninteracting Hofstadter spectrumand treats the superfluid using Bogoliubov theory. Thisgives a mean-field description of the superfluid phase ofthe Bose-Hubbard model, in which the interplay betweenthe ‘magnetic’ vorticity and the lattice potential (analo-gous to pinning effects for a weak lattice ) leads to real-space density modulations (‘supersolidity’). The theoryfurthermore provides a consistent expansion in terms ofboth thermal and quantum fluctuations, and describesthe quasiparticle spectrum above the condensate, whichcombines features of the Hofstadter butterfly and theGoldstone mode characteristic of superfluid order. Wepresent detailed calculations within the superfluid phase,based on the microscopic Hamiltonian, which include pre-dictions for both time-of-flight and Bragg spectroscopymeasurements. A brief outline of our methods and re-sults has been presented elsewhere. Bogoliubov theory has previously been applied to un-frustrated superfluids in optical lattices, and also morerecently to a proposed ‘staggered-flux’ model, whichshares some features of the present analysis. Ðurić andLee have applied a spin-wave analysis to the systemconsidered here, using a real-space perspective that pro-vides results consistent with ours. A related analysis hasbeen applied to fermions, which similarly show spatialorder in the paired superfluid phase.The real-space structure of the condensate also bearsmany similarities to other systems with phase coherencein a magnetic field. In particular, superconducting lat-tices with applied fields support states with current pat-terns similar to those discussed in Section III below, and experiments on these systems have clearly demon-strated the effect of the Hofstadter spectrum on the su-perconducting transition. Frustrated Josephson junc-tion arrays, where the charging energy can becomecomparable to the Cooper-pair tunneling amplitude, pro-vide a close analogue of the present system, and imple-mentations using cold atoms have been proposed.
Previous work on this system has also considered theMott insulator that should exist for strong interactionsand commensurate density.
This phase, whichis very similar to its analogue in the absence of a mag-netic field, is favored when the Hubbard U interactionsuppresses number fluctuations and eliminates the co-herence between neighboring sites of lattice. The tran-sition between this phase and the superfluid has beenstudied using the Gutzwiller ansatz and effectivefield-theory methods, with the main effect being an en-hancement of the insulator phase due to the frustration of the hopping.In this work, we restrict consideration to uniform mag-netic fields. By modulating the effective field in realspace, it is possible to produce states with nontrivialtopological properties. These include topological in- sulators and metals, which display similar physics tothe single-particle states underlying the phenomena de-scribed here. Recent work has also addressed systems ofinteracting bosons in spatially varying magnetic fields. Although the framework that we present has broadapplicability, our specific model includes several simpli-fications. First, we use a single-band Hubbard model,incorporating only the lowest band of the optical-latticepotential. (One effect of the magnetic field is to split thetight-binding spectrum into multiple ‘Hofstadter bands’,and these are included exactly in our approach.) Whilethis is certainly not a valid approximation for the shallowlattices of earlier experiments, the single-band limit isapproached by more recent work with rotating lattices, and should be readily achievable in combination withRaman-induced gauge potentials. We also assume a spatially uniform system with ex-actly rational α . The main effect of the finite trap size isto wash out the small-scale fractal structure of the Hof-stadter butterfly, so we expect most of our resultsto be valid for any α . The experimental consequencesof nonuniformity and incommensurate magnetic flux willbe addressed in more detail in future work. A. Outline
Our approach is as follows: In Section II, we use a sym-metry analysis to find the full spectrum of single-particlestates, described by bosonic annihilation operators a k ℓγ ,defined below in Eq. (18). We then apply the Bogoliubovansatz, a k ℓγ = A ℓγ (2 π ) δ ( k ) + ˜ a k ℓγ , (1)where A ℓγ is a set of c-numbers giving the condensateorder parameter, and the operators ˜ a k ℓγ describe theresidual bosons outside the condensate. An expansion interms of these operators is then developed, which in phys-ical terms is an expansion in fluctuations around mean-field theory. In Section III, we address the lowest-orderterm, which involves only the constants A ℓγ and leadsto a time-independent Gross-Pitaevskii equation for thecondensate configuration.Higher-order terms, describing the spectrum of low-energy excitations and their interactions, are treated inSection IV. In this work, we truncate the expansionat quadratic order, leading to a theory of noninteract-ing Bogoliubov quasiparticles. For this approximation tobe reasonable, one requires a sufficiently low density ofquasiparticles, and we determine the regime of validityby calculating the condensate depletion.In Section V, we describe the consequences for experi-ments, giving predictions for both time-of-flight imageryand Bragg spectroscopy. We conclude with discussion inSection VI. Some details of the Bogoliubov transforma-tion and the symmetry analysis of the interacting spec-tra are given in Appendices A and B. In Appendix C weprovide more details and analytic results for the simplestcase, α = . B. Hamiltonian
The single-band Bose-Hubbard model in the presenceof a static magnetic field can be written as H = − t X h ij i (cid:16) e iΦ ij b † j b i + h.c. (cid:17) + U X j n j ( n j − , (2)where b j and n j = b † j b j are the annihilation and num-ber operators respectively for site j . We treat a two-dimensional square lattice, and assume hopping withamplitude t only between nearest neighbors, denoted h ij i . The second term gives an on-site interaction withstrength U ; the approach described here can straightfor-wardly be extended to incorporate more complicated in-teractions, including between particles on different sites.(Note that we treat electrically neutral particles in thepresence of a synthetic magnetic field, and we will as-sume that there are no long-range interactions.)The phase Φ ij on the directed link from site i to site j , to be denoted i → j , can be expressed in terms of themagnetic vector potential A as Φ ij = Q Z x j x i d r · A ( r ) . (3)The integral is to be taken along the straight-line pathbetween the positions x i and x j of the two sites; one has Φ ij = − Φ ji . While the individual phases depend on thechoice of gauge, the lattice curl X ij (cid:9)(cid:3) Φ ij = Q Z (cid:3) d r · ∇ × A ( r ) (4)is gauge-independent, where the sum is over links enclos-ing a given plaquette (cid:3) in a counterclockwise sense.The integral in Eq. (4) is simply the magnetic fluxthrough the plaquette, given by | B | a , assuming a uni-form magnetic field B = ∇ × A perpendicular to a squarelattice with spacing a . The dimensionless flux per plaque-tte is defined by α = | B | a /φ , and completely specifiesthe effect of the magnetic field. We will henceforth useunits where ~ = a = 1 , and take the effective charge onthe bosons as Q = +1 , so the flux quantum is simply φ = 2 π , and Eq. (4) becomes X ij (cid:9)(cid:3) Φ ij = 2 πα . (5)It is convenient theoretically, and also for recent ex-periments with effective gauge potentials, to use theLandau gauge, in which the vector potential is parallelto one of the square lattice axes. We take the vector po-tential along ˆ y , which is conventional in the theoretical xy H L RI x T x T y Π Α
FIG. 1: An illustration of the vector potential Φ ij in the Lan-dau gauge, and the symmetries defined in Section II A. Thenumber of arrowheads on each link of the square lattice givesthe value of Φ ij on the directed link i → j , in units of πα .The configuration is one possible choice obeying Eq. (5): go-ing around any square plaquette in a counterclockwise sense,the number of forward arrows minus backward arrows is +1 .In this choice of gauge, Φ ij vanishes on all links in the x di-rection. Any choice of gauge reduces the full symmetry ofthe square lattice, and the translation and rotation opera-tions shown are only symmetries when accompanied by phasefactors (see Section II A). Reflections must be combined withtime-reversal operations to preserve the direction of the ap-plied magnetic field. literature, but it should noted that in the recent experi-ments of Lin et al., it is instead aligned with ˆ x . Ourchoice of gauge is illustrated in Figure 1.With this choice, one has e iΦ j,j +ˆ x = 1 and e iΦ j,j +ˆ y =e π i αx j , where (for example) j + ˆ x denotes the site adja-cent to j in the +ˆ x direction. The kinetic term H t thenbecomes H t = − t X j (cid:2) b † j +ˆ x b j + b † j − ˆ x b j + ω x j b † j +ˆ y b j + ω − x j b † j − ˆ y b j (cid:3) , (6)where ω = e π i α . Note that the phases act to ‘frustrate’the hopping, so that for noninteger α it is not possible tominimize the kinetic energy on every link of the latticesimultaneously. The transformation to an alternative gauge is imple-mented by applying a spatially-varying phase rotationto b j . For example, in the symmetric gauge, which ismore natural for a rotating lattice, one defines opera-tors ˜ b j = ω − x j y j / b j , so that the gauge field becomes e iΦ j,j +ˆ x = ω − y j / and e iΦ j,j +ˆ y = ω x j / . The numberoperator n j = b † j b j is invariant under any such gaugetransformation, as required for a physically observablequantity, and so the interaction energy H U is also invari-ant.As noted above, we assume precisely rational α = p/q (with p and q coprime), so that ω q = 1 , and the unit cellof H is q × sites. (In the symmetric gauge, the unit cellis considerably larger, containing q × q sites.) II. SINGLE-PARTICLE SPECTRUM
We begin by describing in detail the spectrum of thesingle-particle kinetic term H t , given in Eq. (6). Whilethese results are well established (and have beenfor many decades), we will present them here in somedetail, in order to introduce the concepts and formalismthat will be central to our subsequent analysis.The general structure of the spectrum for any α = p/q can be determined by considering the “magnetic symme-try group” (MSG; also known as the “projective symme-try group”). This arises because any choice of gaugenecessarily reduces the physical symmetry of the latticeand leads to a Hamiltonian that does not commute withthe standard spatial symmetry operators. For exam-ple, while the (uniform) magnetic field is invariant un-der translation by a single lattice site in any direction, H t is manifestly asymmetric under translations in the x direction.One can nonetheless define a group of operators inone-to-one correspondence with the physical symmetries,that obey the multiplication table apart from phase fac-tors, and that commute with H t (and indeed the fullHamiltonian H ). We define this group by specifying op-erators corresponding to the elementary translation, ro-tation, and reflection transformations from which the fullgroup can be constructed. A. Magnetic symmetry group
We first define the elementary translation operators T x and T y , illustrated in Figure 1, in terms of their commu-tation with the annihilation operator b j . With our choiceof the Landau gauge, the Hamiltonian is symmetric un-der translations in the y direction, and so one can define T y by T y b j = b j +ˆ y T y . (7)In contrast, the x -dependent phase factors in Eq. (6) im-ply that a pure translation does not commute with H t ,and we instead define T x by T x b j = b j +ˆ x T x ω − y j . (8)The gauge field has periodicity q in the x direction, sothe combination T qx obeys T qx b j = b j + q ˆ x T qx . The multiplication relations of the MSG are equal upto phase factors to those of the ordinary spatial group.With total particle number N = P j n j , one finds T x T y = T y T x ω N , (9)by induction, starting from the (totally symmetric) vac-uum state. While the phase factors associated with in-dividual operators are dependent on the choice of gauge,this relation is gauge independent.Besides translations, it is useful to consider rotationand reflection operations. We define the unitary operator R giving a rotation by ◦ counterclockwise about thesite at the origin (see Figure 1): R b j = b R j ω x j y j R , (10)where j → R j under the rotation ( x R j = − y j , y R j = x j ).The phase factor again ensures that [ R , H t ] = 0 .For reflection operators, the situation is somewhat dif-ferent, since the magnetic field explicitly breaks chiralityand time-reversal symmetry. While reflection reverseschirality and so is not a symmetry of the Hamiltonian, acombination of reflection and time reversal restores theappropriate sense of circulation and remains a symme-try in the presence of the field. One can therefore de-fine antiunitary operators corresponding to such com-binations; we define I x and I y for the transformationsobeying x I x j = − x j and y I y j = − y j respectively. Nophase factors are required in these cases.It is also useful to define the combinations P = R = I x I y and I xy = RI y . The former gives inversion aboutthe origin, x → − x , which in two dimensions is a properrotation and hence represented by a unitary operator.The antiunitary operator I xy gives reflection in the line y = x .Note that the interaction term H U is a function only ofthe gauge-invariant combination n j = b † j b j and so is un-affected by the phase factors included in expressions suchas Eq. (8). Any interaction term with the full symmetryof the lattice, such as the explicit example in Eq. (2), istherefore invariant under the magnetic symmetry group. B. Momentum-space operators
While the unit cell of the Hamiltonian H t contains q sites and hence gives a reduced Brillouin zone, it is usefulto construct momentum space operators based on thefull lattice Brillouin zone B L : − π ≤ k x , k y < π . Themomentum-space annihilation operator b k is defined by b k = X j e − i x j · k b j , (11)so that the commutator is given by [ b k , b † k ′ ] = (2 π ) δ ([ k − k ′ ] B L ) . (12)Here and throughout, we use the notation [ k ] B to denotethe momentum k reduced to the Brillouin zone B by theaddition of an appropriate lattice vector. We also use the shorthand notation [ x ] q = x mod q .Written in momentum space, the Hamiltonian H t is H t = − t Z k ∈ B L d k (2 π ) (cid:18) k x b † k b k + e − i k y b † [ k + X ] B L b k + e i k y b † [ k − X ] B L b k (cid:19) , (13)where X = 2 πα ˆ x . The enlarged unit cell ( q × sites) of the Hamiltonian allows mixing between momentum statesthat coincide when reduced to the magnetic Brillouin zone B M : − π ≤ k y < π , − π/q ≤ k x < π/q .The operators b k commute with T y , T y b k = e i k y b k T y , (14)but the phase factor in the definition of T x causes it to mix momenta, T x b k = e i k x b [ k + Y ] B L T x , (15)where Y = 2 πα ˆ y . Since T x commutes with H t , this implies degeneracies in the single-particle spectrum betweenpoints separated by Y . To make this transparent, it is convenient to define the doubly reduced Brillouin zone B N : − π/q ≤ k x , k y < π/q .Any point within B L can be specified as [ k + ℓ X + n Y ] B L , where k ∈ B N and n and ℓ are integers from to q − ,so H t can be rewritten as H t = Z k ∈ B N d k (2 π ) q − X ℓ =0 q − X n,n ′ =0 b † [ k + n X + ℓ Y ] B L H nn ′ ([ k + ℓ Y ] B L ) b [ k + n ′ X + ℓ Y ] B L , (16)where the q × q matrix H ( k ) has elements (for q > ) H nn ′ ( k ) = − t × e +i k x ω n + e − i k x ω − n if n ′ = n e +i k y if n ′ = [ n + 1] q e − i k y if n ′ = [ n − q otherwise. (17)(If q = 2 , [ n +1] q = [ n − q and H = H = − t cos k y .)Noting that H nn ′ ([ k + ℓ Y ] B L ) = ω ( n − n ′ ) ℓ H nn ′ ( k ) , wediagonalize H t by writing b [ k + n X + ℓ Y ] B L = ω − nℓ X γ ψ γn ( k ) a k ℓγ , (18)where a k ℓγ is the annihilation operator for a single-particle state labeled by band index γ ∈ { , . . . , q } . Foreach k ∈ B N , ψ γn ( k ) is an eigenvector of H ( k ) witheigenvalue ǫ γ ( k ) . The q bands can be understood fromthe ‘folding’ of the Brillouin zone due to the reducedtranslation symmetry of H .It should be noted that the annihilation operator a k ℓγ has momentum k when referred to B N , but momentum k + ℓ Y in the magnetic Brillouin zone B M . For each k ∈ B N , the eigenvectors ψ γn ( k ) corresponding to dif-ferent bands are orthogonal, so the operators a k ℓγ obeycanonical commutation relations, [ a k ℓγ a † k ′ ℓ ′ γ ′ ] = (2 π ) δ ([ k − k ′ ] B N ) δ ℓℓ ′ δ γγ ′ . (19) The single-particle Hamiltonian can finally be rewrit-ten as H t = Z k ∈ B N d k (2 π ) q − X ℓ =0 X γ ǫ γ ( k ) a † k ℓγ a k ℓγ . (20)The single-particle energy ǫ γ ( k ) is independent of ℓ , soevery state is (at least) q -fold degenerate. The band la-bels γ can be arranged so that ǫ γ ( k ) ≤ ǫ γ +1 ( k ) for every k and ǫ γ ( k ) is a continuous function of k .Using Eqs. (14) and (15), one finds the effect of theoperators T y and T x as T y a k ℓγ = e i k y ω ℓ a k ℓγ T y (21) T x a k ℓγ = e i k x a k [ ℓ +1] q γ T x . (22)The operator T x therefore transforms one degeneratesingle-particle state into another, and it is this symmetrythat enforces the degeneracy.Determining the effects of the rotation and reflectionoperators R , I x , and I y is somewhat more involved. Considering first R , its commutation with H t impliesthat R a k ℓγ R † can be written in terms of a ( R k ) ℓ ′ γ in thesame band γ . With an appropriate choice of the arbi-trary phase of the eigenvectors ψ γn at the points k and R k , one finds R a k ℓγ = 1 √ q e − i φ γ X ℓ ′ ω ℓℓ ′ a ( R k ) ℓ ′ γ R , (23)where the phase φ γ is independent of k . The requirementthat R = 1 implies that φ γ is a multiple of π/ for allbands γ . One similarly finds I x a k ℓγ = a ( I y k )[ − ℓ ] q γ I x (24) I y a k ℓγ = e − φ γ a ( I x k ) ℓγ I y (25) P a k ℓγ = e − φ γ a ( − k )[ − ℓ ] q γ P . (26)Note that e − φ γ = ± is real. Similar expressions, suchas ψ γn ( R k ) = 1 √ q e − i φ γ X n ′ ω nn ′ ψ γn ′ ( k ) , (27)relate the eigenvectors ψ γn ( k ) at symmetry-equivalentmomenta k .Applying these operators to the Hamiltonian in theform of Eq. (20) immediately shows that the single-particle dispersion ǫ γ ( k ) is symmetric under the corre-sponding transformations of the momentum k . For ex-ample, ǫ γ ( R k ) = ǫ γ ( k ) , which implies that the dispersionhas the full four-fold rotation symmetry of the lattice, de-spite the reduced symmetry of H t . C. Spectrum
In summary, for α = p/q , the single-particle spectrumconsists of q bands, labeled by γ , with each state q -folddegenerate. These degenerate states have energy ǫ γ ( k ) and momentum [ k + ℓ Y ] B M , with ℓ ∈ { , . . . , q − } .Figure 2 shows the ‘Hofstadter butterfly’, a plot of theallowed single-particle energies ǫ γ (for any momentum k ) as a function of the flux α . The plot has a fractalstructure that is sensitively dependent on α , and forclarity only rational α = p/q with q ≤ have beenincluded. For each α = p/q , points mark the top andbottom of each of the q bands.As can be seen in Figure 2, most of the bands areseparated by nonzero gaps. The only exceptions are thetwo central bands for q even, which touch exactly at thepoint of zero energy. These occur at k = for q aninteger multiple of , and at the corner of B N , k X = πq ˆ x + πq ˆ y , otherwise. For both, the spectrum has a linear‘Dirac-cone’ dispersion near the degeneracy point. In all other cases, including the lowest band for all α ,the dispersion is quadratic near its minimum. The mini-mum of the lowest band always occurs at k = , as canbe shown using the Perron-Frobenius theorem, and theeffective mass near this point can be found by perturba-tion theory for small k . In all cases, the coefficients areequal in the x and y directions (i.e., the effective mass Α- Ε Γ (cid:144) t FIG. 2: The Hofstadter butterfly, a plot of the single-particle energies ǫ γ as a function of flux α . The butterflyhas a fractal structure, but only a finite set of α = p/q canbe plotted; for clarity, we restrict to q ≤ . Points mark thetop and bottom of each of the q bands, which become increas-ingly narrow as q becomes larger. In the limit q → ∞ with p fixed, or equivalently α ≪ , the low-lying bands become theLandau levels of the continuum. tensor is proportional to the unit matrix), a straightfor-ward consequence of the symmetry R . Figure 3 shows thelowest band of the noninteracting dispersion, ǫ ( k ) , for α = and . In both cases, the dispersion is quadraticand isotropic near the top and bottom of the band. D. Interactions
The operators a k ℓγ defined in Section II B are chosen todiagonalize the kinetic energy operator H t . To incorpo-rate the effects of the interaction term H U , this must alsobe expressed in terms of these operators. This involvesthe straightforward process of substituting Eqs. (11) and(18) into H U , and can be performed for any choice ofinteraction. Our explicit calculations are for the on-siteHubbard interaction in Eq. (2), appropriate to bosons ina deep optical lattice.A general quartic interaction can be written in termsof the operators a k ℓγ as H U = Z k ··· k X ℓ ··· ℓ X γ ··· γ u a † k ℓ γ a † k ℓ γ a k ℓ γ a k ℓ γ ,(28)where the coefficient u is a function of the four sets ofindices k , ℓ , and γ . It can be chosen symmetric underexchange of the first two ( ↔ ) or last two sets ( ↔ ),and the requirement that H U be hermitian implies that u (3 , , ,
2) = u ∗ (1 , , , .A translation-invariant interaction conserves momen-tum, so the coupling coefficient u is nonzero only if (cid:2) ( k + k − k − k ) + ( ℓ + ℓ − ℓ − ℓ ) Y (cid:3) B M = . (29) - - - - - G MX - Π Π -Π- Π Π Π G MX - Π Π -Π- Π Π Π - - - - - FIG. 3: Contour plot of the lowest band of the noninteractingsingle-particle spectrum t ǫ ( k ) , for α = pq = (left) and (right). In both cases, the spectrum is plotted in the magneticBrillouin zone B M , in which − πq ≤ k x < πq and − π ≤ k y < π .Points separated by momentum Y = 2 πα ˆ y are degenerate;identifying these gives the doubly reduced Brillouin zone B N , − πq ≤ k x , k y < πq , indicated by the horizontal dashed lines.Note that the lowest band is considerably narrower in the case α = , as is also evident in Figure 2. Note that this allows for umklapp processes where the netmomentum is zero only when reduced to B N . Factoringout (2 π ) δ ([ k + k − k − k ] B N ) gives ¯ u : ¯ u ( { k } , { ℓ } , { γ } ) = U δ { ℓ } X n ··· n δ { n } ψ ∗ γ n ( k ) ψ ∗ γ n ( k ) ψ γ n ( k ) ψ γ n ( k ) ω n ℓ + n ℓ − n ℓ − n ℓ , (30)where δ { ℓ } and δ { n } denote Kronecker deltas enforcing Eq. (29) and a similar constraint on n ··· . The complicatedstructure of the single-particle states gives ¯ u a nontrivial dependence on the momenta k of the interacting particles,despite the choice of a purely on-site interaction.Note that ¯ u is, up to a phase factor that is only nontrivial for umklapp processes, only dependent on two of the ℓ ’s: ¯ u ( ℓ , ℓ , ℓ , ℓ ) = e − i p ¯ pℓ ( k + k − k − k ) . ˆ x ¯ u (0 , [ ℓ − ℓ ] q , [ ℓ − ℓ ] q , [ ℓ − ℓ ] q ) , (31)where ¯ p is the modulo- q reciprocal of p , the integer such that [ p ¯ p ] q = 1 and < ¯ p < q . This identity, and themomentum-conservation constraint of Eq. (29) are consequences of the symmetry of the interactions under T x and T y . Further constraints on the coefficients u result from the symmetry properties of the eigenvectors ψ γn ( k ) underrotations and reflections. For example, requiring that P commutes with H U and using Eq. (26) gives e − P i =1 φ γi u ( k . . . k , ℓ . . . ℓ , γ . . . γ ) = u ( − k . . . − k , [ − ℓ ] q . . . [ − ℓ ] q , γ . . . γ ) . (32)This implies that u is odd in momentum, and hence van-ishes for k ··· = , for certain combinations of ℓ ··· and γ ··· . A detailed discussion of the restrictions imposedby symmetries has been presented by Balents et al. ina different context. III. MEAN-FIELD THEORY
Having described the spectrum of noninteracting par-ticles, we now turn to the effects of interactions on themany-body physics. As described in Section I A, our ap-proach will be based on the Bogoliubov theory, using theansatz Eq. (1) and performing an expansion in powersof the fluctuation operators. The Bogoliubov ansatz canbe viewed as a statement about correlation functions inthe superfluid phase, with the first term in Eq. (1) giv-ing the one-point correlation function, the condensate or-der parameter h a k ℓγ i = A ℓγ (2 π ) δ ( k ) . The zeroth-orderterm in the Bogoliubov expansion is given by neglectinghigher-order connected correlation functions.This mean-field theory can be viewed as a specialcase of that derived by using the Gutzwiller ansatz,which assumes a state Q j | ψ j i that is factorizable in realspace. In general, one allows | ψ j i to be an arbi-trary state within the on-site manifold, but our ansatzassumes a bosonic coherent state and is appropriate onlywithin the superfluid phase.Substituting the mean-field ansatz into the Hamilto-nian, written as in Eqs. (20) and (28), gives the energydensity h = X ℓ,γ A ∗ ℓγ A ℓγ [ ǫ γ ( ) − µ ]+ X { ℓ } , { γ } ¯ u ( { } , { ℓ } , { γ } ) A ∗ ℓ γ A ∗ ℓ γ A ℓ γ A ℓ γ , (33)where the interaction strength ¯ u is evaluated with allmomenta equal to zero. [This expression has been di-vided by a factor of (2 π ) δ ( ) , corresponding physicallyto the system volume.] The mean-field condensate con-figuration can be found by minimizing h with respectto A ℓγ . The resulting equation can be viewed as a time-independent Gross-Pitaevskii equation for the conden-sate wavefunction in momentum space.The corresponding real-space wavefunction can befound using Eqs. (11) and (18) and is given by h b j i = X ℓn ω nx j + ℓy j − nℓ X γ ψ γn ( ) A ℓγ . (34)This is in general a function of [ x ] q and [ y ] q , and so givesa q × q site unit cell in real space. To this order, the par-ticle density is simply given by h n j i = |h b j i| . Note thatthe presence of A ℓγ for nonzero ℓ implies that the con-densate contains components for k = ℓ Y = and hencethat spatial symmetry is broken. Within this mean-fieldtheory, this spatial order develops simultaneously withthe breaking of phase-rotation symmetry, and is a sim-ple consequence of the degeneracy in ℓ . (In fact, thefinite-temperature transition in two dimensions is of theBerezinskii-Kosterlitz-Thouless type, and so this partic-ular result is not necessarily reliable.)The configuration of currents within the superfluidphase can be calculated using the gauge-invariant cur-rent operator for the link i → j , J ij = i t e iΦ ij b † j b i + h.c. . (35)Figure 4 shows the currents hJ ij i in the mean-field con-densate configurations for α = q with ≤ q ≤ ; theyhave the same q × q unit cell as the condensate wave-function. For the larger values of q , particularly α = , FIG. 4: Example mean-field condensate configurations for α = (top left), (top right), (bottom left), and (bot-tom right). The blue arrows show the direction of the current hJ ij i on each link i → j of the lattice, and their lengths indi-cate the magnitude (the length scale is not consistent betweendifferent values of α ). The black points show the positions ofthe lattice sites i and have area proportional to the density h n i i . In each case, the condensate reduces the spatial symme-try of the square lattice, and is one member of a discrete setof degenerate configurations related by the action of the bro-ken symmetries. The degeneracies and residual symmetriesare listed in Table I. The quantitative details, but not thesymmetries, depend on the interaction strength U ; the plotsshow the case U ≪ t . these resemble Abrikosov lattices: the plaquettes withlow density and high current circulation can be viewedas containing vortices. The symmetry properties of theseconfigurations are listed in Table I.An analysis of the patterns that are allowed for generalinteractions and various values of q has been given byBalents et al., who considered the same problem in adifferent context. In the present case, it is valid to assumepurely on-site interactions, allowing the ordered states tobe determined unambiguously.Minimization of h with respect to A ℓγ is equivalentto minimizing with respect to the real-space condensatewavefunction or vortex configuration. The latter perspec-tive is more appropriate in the continuum, whereas herethe lattice potential provides a strong pinning potentialthat simplifies the momentum-space approach.Our ansatz for the condensate configuration, Eq. (1),also involves bands γ other than the lowest. Occupa-tion of higher bands costs kinetic energy, increasing thefirst term of Eq. (33), and so is disfavored when inter- α A ℓ degeneracy symmetries (1 i) 2 T y T x , T x R , I xy (1 ω ω ) 6 T y T x , R , I xy (cid:16) p ω/ − i p ω/ (cid:17) R , T y T x I y (1 ω ω ∗ ω ∗ ω ) 10 T y T x , R TABLE I: Properties of the mean-field condensate configu-rations shown in Figure 4, including their degeneracies andunbroken symmetries. The vectors in the column labeled A ℓ are configurations minimizing h in the limit of weak inter-actions, U/t → , where the condensate is restricted to thelowest band, γ = 1 . The residual symmetry group is given byproducts of powers of the operators listed, along with T qx and T qy , which are always preserved by the ansatz of Eq. (1). Thesymmetries T x , T y , and R are illustrated in Figure 1, and thecombination I xy = I x R gives a reflection in the line y = x . actions are very weak. For stronger interactions, the en-ergy is reduced by smoothing out density fluctuations,which requires incorporating higher bands into the con-densate. This competition between kinetic and potentialenergy also allows for first-order transitions between dif-ferent local minima of h as the interaction strength ormean density varies. We have not found any examplesfor q ≤ , however, and their observation in experimentswould anyway likely require considerable enhancementsin stability and cooling.The mean-field energy h is symmetric under the sametransformations of A ℓγ as the full Hamiltonian is un-der transformations of a ℓγ , as discussed in more de-tail in Appendix B. As noted above, certain symme-tries are spontaneously broken by the condensate con-figuration, and the corresponding operators transform agiven A ℓγ into a symmetry-equivalent degenerate con-figuration. The degeneracies of the patterns shown inFigure 4 are listed in Table I. The number of degener-ate configurations is in every case a multiple of q , as weprove in Appendix B. The degeneracy in the orderingpatterns allows for the possibility of real-space domainformation, which would not affect time-of-flight imagesand would likely require more sophisticated in situ probesto confirm. It should be noted that the ansatz of Eq. (1) implicitlyexcludes ordered states with larger unit cells than q × q sites. (Previous work using a real-space approach hassuggested that this may happen for α = .) It is straight-forward to include such states (with larger but finite unitcells) at the mean-field level, by allowing nonzero conden-sate amplitude at a discrete set of momenta [ k ] B N = .This complicates somewhat the analysis that follows, and we will not treat this possibility further.The condensate configuration A ℓγ also determines theoccupation numbers in momentum space, and so canbe used to predict the result of a time-of-flight expan-sion measurement, as discussed in detail in Section V A.Briefly, the terms in the Hamiltonian Eq. (13) mixingmomenta differing by X imply Bragg peaks at pointscorresponding to momenta n X , while the nonzero con-densate amplitude A ℓγ for ℓ = 0 gives further peaks at ℓ Y . (It should be recalled that our axes are reversedfrom those of Lin et al. ) IV. BOGOLIUBOV THEORY
The mean-field theory of Section III results from us-ing the Bogoliubov ansatz of Eq. (1) and keeping only thelowest-order term in an expansion in terms of the fluctua-tion operators ˜ a k ℓγ . To improve upon this theory and de-termine the spectrum for single-particle excitations abovethe condensate, we consider in this section the followingorder in the expansion. The terms containing a singleoperator vanish when the mean-field energy density h isminimized, and so we next treat the quadratic terms.The quadratic part of the Hamiltonian H (2) can beconveniently expressed in matrix form, by combining cre-ation and annihilation operators into a column vector, α k ℓγ = ˜ a k ℓγ ˜ a † ( − k ) ℓγ ! . (36)One can then write H (2) as H (2) = 12 Z k ∈ B N d k (2 π ) X ℓℓ ′ ,γγ ′ α † k ℓγ M ℓγ,ℓ ′ γ ′ ( k ) α k ℓ ′ γ ′ + H (2) c , (37)where H (2) c is a term that contains no operators andcomes from a commutator.This expression can straightforwardly be generalizedto allow for other choices of single-particle basis, by re-placing ℓ and γ by a single generic index λ . The matrix M ( k ) can then be written as M λλ ′ ( k ) = [ ǫ λλ ′ ( k ) − µδ λλ ′ ] + B λλ ′ ( k ) , (38)where ǫ λλ ′ , the generalization of ǫ γ δ γγ ′ δ ℓℓ ′ , is not diago-nal in the general case, and B λλ ′ ( k ) = X λ λ u ( k , , k , ; λ, λ , λ ′ , λ ) A ∗ λ A λ u ( k , − k , , ; λ, λ ′ , λ , λ ) A λ A λ u ∗ ( k , − k , , ; λ ′ , λ, λ , λ ) A ∗ λ A ∗ λ u ( − k , , − k , ; λ ′ , λ , λ, λ ) A ∗ λ A λ ! . (39)0Similarly to h in Eq. (33), M ( k ) has contributionsfrom both the kinetic and potential energy. The lattercan be viewed as self-energy terms for the quasiparticlesdue to scattering with bosons in the condensate. They in-clude ‘anomalous’ processes in which a pair of condensedparticles scatter from each other into an excited state andthe reverse process where they return to the condensate.These result in the off-diagonal elements in Eq. (39),giving terms in H (2) that do not conserve the numberof ˜ a k ℓγ quanta. The standard Bogoliubov theory forzero magnetic field is recovered by taking q = 1 , in whichcase the ℓ and γ indices are redundant.It is useful to consider M ℓγ,ℓ ′ γ ′ ( k ) as a q × q ma-trix (for each k ). Using the properties of u given afterEq. (28), one can show that this matrix is hermitian.The zero-momentum limit M ( ) is the Hessian of themean-field term h (with respect to variations in A ℓγ andits conjugate) and so is a nonnegative-definite matrix.The single vanishing eigenvalue corresponds to the bro-ken U(1) symmetry of h . For nonzero k , all eigenvaluesof M ( k ) are strictly positive. A. Bogoliubov quasiparticles
To find the spectrum of quasiparticles, one must definea new set of annihilation and creation operators in termsof which H (2) is diagonal. Momenta (referred to B N ) arenot mixed in Eq. (37), so the new operators are labeledby k , but since the condensate breaks symmetry under T y , ℓ is no longer a good quantum number. We there-fore define annihilation operators for these modes as d k ζ ,where ζ ∈ { , , . . . , q } . (In certain cases, there are un-broken translation symmetries, as shown in Table I. Thestates can then be labeled by the eigenvalues of the cor-responding operators, as discussed in Appendix B.)In order to preserve the bosonic commutation rela-tions, [ d k ζ , d † k ′ ζ ′ ] = (2 π ) δ ( k − k ′ ) δ ζζ ′ , the transforma-tion between α k ℓγ and d k ζ must be symplectic; detailsare given in Appendix A. In terms of the new operators,the quadratic part of the Hamiltonian is given by H (2) = Z k ∈ B N d k (2 π ) X ζ ξ k ζ (cid:20) d † k ζ d k ζ − (2 π ) δ ( ) X ℓγ Y ζ k ℓγ ∗ Y ζ k ℓγ (cid:21) , (40)where d k ζ = X ℓγ h X ζ k ℓγ ∗ ˜ a k ℓγ − Y ζ k ℓγ ∗ ˜ a † ( − k ) ℓγ i . (41)To this order, the system is therefore described by non-interacting Bogoliubov quasiparticles with annihilationoperators d k ζ and energies ξ k ζ > . Because of the off-diagonal elements in M ( k ) , the quasiparticles are super-positions of particles and holes, with X ζ k ℓγ and Y ζ k ℓγ re-spectively giving these components. G M X G G MX Π q Π (cid:144) q FIG. 5: Quasiparticle dispersion (solid lines) and noninter-acting single-particle dispersion (dashed), both in units ofhopping t , for α = pq = . An analytic expression for thespectrum for this case is given in Eq. (C7) of Appendix C.The dispersions are plotted along a path in the reduced Bril-louin zone B N shown in the left inset. In the interacting case, U = 4 t , the mean density is ρ = 1 , and the real-space con-figuration is as shown in the right inset (see also Figure 4).In both cases there are q = 4 modes, including, in the in-teracting case, one Goldstone mode with linear dispersion.For U = 0 , the modes are q -fold degenerate and have beenshifted vertically by an arbitrary choice of chemical potential.At k = π ˆ x + π ˆ y , the corner X of B N , the modes meet ata point, with a linear dispersion. Such a ‘Dirac cone’ occurswhenever q is even. Other notable features of the interactingspectrum include a twofold degeneracy along the line from Mto X, and the unshifted modes (relative to the noninteractingdispersion) from X to Γ . The former can be understood as aKramers degeneracy due to the antiunitary symmetry under T x I y , as discussed in Appendix B. For vanishing interactions, the second term in Eq. (38)is absent and the quasiparticle spectrum ξ k ζ is identicalto the single-particle spectrum ǫ γ ( k ) described in Sec-tion II C. With nonzero interactions, the q -fold degen-eracy within each band is split and, for generic k , thespectrum consists of q distinct modes. The quasiparticledispersion for α = and are shown in Figures 5 and 6respectively, along with the noninteracting single-particlespectrum. (For clarity, only out of the q = 9 modesare shown in the latter case.) For general q , the diago-nalization of M ( k ) must be performed numerically, butthe simplest case is analytically tractable, and is treatedin detail in Appendix C.As the figures show, the lowest energy approaches zeroin the limit k → , giving the Goldstone mode thatresults from broken U(1) symmetry. For small | k | this‘phonon’ has a linear dispersion, and can be described interms of long-wavelength fluctuations of the condensatephase. A low-energy theory of the mode can be foundby allowing gradual deviations from the mean-field valueof the phase and its conjugate density, as described inSection C 2.1 X' G M X G G M XX' Π q Π (cid:144) q FIG. 6: Quasiparticle dispersion (solid lines) and noninter-acting single-particle dispersion (dashed), for α = pq = . Inthe interacting case, U = 2 t , the mean density is ρ = 1 . Inboth cases, there are q = 9 modes, of which only the lowest are shown. The dispersions are plotted along a path in thereduced Brillouin zone B N shown in the left inset. Becausethe condensate configuration breaks the symmetry under R ,rotation by π , the quasiparticle dispersions along the linesfrom Γ to X and from Γ to X ′ are different. The interact-ing spectrum includes several (unavoided) level crossings, forexample in the lower band near X ′ , a result of the unbrokentranslation symmetry T y T x , as discussed in Appendix B. The phase velocity c of the Goldstone mode, whichis expressed in terms of the spectrum at k = in Sec-tion A 3, is in many cases independent of direction, in-cluding both α = and . This isotropy is a straight-forward result of the symmetry in both cases under I xy = RI y , as noted in Table I. This property, andother features of the spectra shown in Figures 5 and 6,are discussed in Appendix B.The second term in Eq. (40) includes H (2) c fromEq. (37) and represents the change in the zero-point en-ergy associated with the superfluid state. This term,coming from quantum fluctuations, is accompanied atnonzero temperature T by a contribution from thermallyexcited Bogoliubov quasiparticles, leading to a free en-ergy per site of ∆ f = Z k ∈ B N d k (2 π ) X ζ (cid:20) T log(1 − e − ξ k ζ /T ) − ξ k ζ | Y ζ k ℓγ | (cid:21) ,(42)where we use units such that k B = 1 . (Within the mean-field theory of Section III, all particles are in the conden-sate, so the entropy vanishes and the free energy is givenby h .) These contributions in principle allow a configu-ration with a higher mean-field energy h to be selectedbecause of its enhanced fluctuations and hence lower freeenergy h + ∆ f . It should be noted, however, that thedegeneracy of the symmetry-equivalent condensate con-figurations discussed in Section III cannot be lifted by ∆ f . The calculated spectra lead to important experimentalpredictions, as discussed below in Section V. Occupa-tion of the quasiparticle modes, due to both thermal andquantum fluctuations, gives the structure of time-of-flightimages away from the Bragg peaks mentioned previously,and spectroscopic methods should be able to measure themode dispersions ξ k ζ directly (see Section V B). B. Condensate depletion
The ansatz of Eq. (1) and the expansion in powersof operators is in principle exact, with the higher-orderterms leading to interactions between the Bogoliubovquasiparticles. Here, the series is truncated at quadraticorder, an approximation that is valid provided that thequasiparticles remain at sufficiently low density for theirinteractions to be neglected.This criterion can be quantified by calculating the de-pletion of the condensate, equal to the quasiparticle con-tribution to the total particle number. This is found byexpressing the number operator n j is terms of the quasi-particle operators d k ζ , summing over sites j , and takingthe ensemble average. The mean particle density is thengiven by ρ = X ℓγ | A ℓγ | + Z k ∈ B N d k (2 π ) X ζℓγ (cid:26) | X ζ k ℓγ | n B ( ξ k ζ )+ | Y ζ k ℓγ | [1 + n B ( ξ k ζ )] (cid:27) , (43)where n B ( ξ ) = (e ξ/T − − is the Bose-Einstein distribu-tion function. The first term in Eq. (43) is the condensatedensity, and is simply the spatial average of the mean-field density calculated in Section III, while the secondterm gives the average density of particles outside thecondensate. The relative magnitude of these two termsgives a measure of the significance of fluctuations, and weuse the ratio of the second to the first as our definitionof the depletion.Figure 7 shows the depletion for α = , as a functionof density, interaction strength, and (in the inset) tem-perature. It is small deep within the superfluid phaseand increases to roughly for the largest values of U and T shown. Neglecting cubic and quartic termswithin the Bogoliubov theory relies on the assumptionof small depletion, and so the conclusions presented hereare only qualitatively applicable for larger values of U and T . (The order of magnitude is consistent with thespin-wave analysis of Ðurić and Lee. )For zero temperature, n B ( ξ >
0) = 0 and the onlyfluctuation contribution is from the second term withinthe braces [zero-point fluctuations; compare Eq. (42)]. Inthis case, the integrand diverges as v | k | − for small | k | ,where v is a coefficient in the expansion of Y ζ k ℓγ for small k (see Appendix A 3). The integral is therefore finite inthis case.2 Ut Tt FIG. 7: (Color online) Condensate depletion for α = asa function of interactions U/t (main figure) and temperature
T /t (inset), where t is the hopping strength. In the main fig-ure, T = 0 and the densities are ρ = 1 (top curve), (middle),and (bottom), while in the inset, ρ = 1 and U/t = 2 . Thedepletion is smallest, and hence the approximation best, deepin the superfluid phase, with weak interactions, high density,and low temperature. For
T > , small-momentum cutoffs of k = 0 . (solid line) and k = 0 . (dashed line), in latticeunits, have been used to remove the logarithmic divergence ofthe depletion integral. For
T > , the Bose-Einstein distribution function be-comes T ( c | k | ) − for small | k | . This results in a logarith-mically divergent integral, an instance of the Mermin-Wagner-Hohenberg theorem, which states that, intwo dimensions for nonzero temperature, the continuousphase symmetry cannot be broken. In an infinite two-dimensional system, there is no true condensate and sothe ‘depletion’ is complete.In the presence of an external trapping potential, how-ever, nonzero temperature condensation is possible evenin two dimensions. This can be captured in a crude wayby applying a small-momentum cutoff k on the integralover k , with k ≃ R − eff , where R eff is the effective radiusof the system in the trap (in units of the lattice spacing).If the two-dimensional plane is embedded within a deeplattice in the z direction, then hopping in this transversedirection can also stabilize the condensate. An appro-priate momentum cutoff is then given by k ≈ √ m ∗ t ⊥ ,where t ⊥ is the transverse hopping matrix element, andhence the energy scale over which the system appearsthree-dimensional, and m ∗ is the effective mass at theminimum of the lowest band in the single-particle dis-persion.In either case, the depletion integral is finite, with alogarithmic dependence on k of ρ log = v T πc log k . (44)In the inset of Figure 7, the depletion is shown at nonzerotemperature, using two different values of k . The dif-ference between the two curves is well approximated for most values by Eq. (44). (The exact difference involvesother terms that are not singular at k = 0 .) To deter-mine the cutoff used in the plot, we have assumed thatthe finite system size will be the most important effect,and taken R eff ≃ – lattice sites. The depletion calculation also provides a rough esti-mate for the boundary of the superfluid phase, at thepoint where the depletion reaches , although theapproximation of independent quasiparticles is proba-bly not valid at this point. For α = , ρ = 1 and T = 0 , this gives an estimate of ( t/U ) c = 0 . , inreasonable agreement with the value of ( t/U ) c = 0 . (at the tip of the ρ = 1 Mott lobe) found using theGutzwiller ansatz.
It should be noted that the lat-ter approach, which neglects fluctuations within the Mottinsulator, generally underestimates ( t/U ) c . V. EXPERIMENTAL PREDICTIONS
The Bogoliubov theory that we have presented for thesuperfluid phase provides several concrete predictions forexperiments, most notably for time-of-flight images andBragg spectroscopy.
A. Time-of-flight images
As noted above in Section III, the enlarged unit cell ofthe condensate has important consequences for time-of-flight images. The corresponding reduction of the Bril-louin zone leads to additional Bragg peaks that give aclear indication of the formation of spatial order in thecondensate. The intensity away from these peaks is de-termined by bosons excited to states with [ k ] B N = bythermal and quantum fluctuations.In a time-of-flight measurement, the trapping poten-tial confining the atoms within the lattice is suddenlyswitched off, causing a rapid expansion. After a fixedperiod of the time, the density profile of the cloud isdetermined, for example by illuminating the atoms andmeasuring the transmitted intensity. If the interactionsbetween the atoms during the expansion are sufficientlyweak, then it can be treated as ballistic, and we will as-sume that this is the case throughout. In the absenceof a magnetic field, the density profile after a fixed timeof flight measures the original momentum distribution inthe trap. The same is true with a field, apart from somemodifications that we discuss in the following.The time-of-flight images depend on certain details ofthe experiment and, in particular, the means used toproduce the effective magnetic field. In the case of arotating system, the momentum in the stationary (labo-ratory) frame is equal, up to a possible global rotation,to the symmetric-gauge canonical momentum in the ro-tating frame.With a Raman-induced gauge field, the results dependon whether the Raman beams remain after release; we3assume that they are suddenly switched off simultane-ously with the trap, such as in the experiments of Lin etal. The trajectory of an atom is determined by the mo-mentum immediately after the gauge field is switched off,which is equal, using the sudden approximation, to theLandau-gauge canonical momentum before switch-off. Within the approximation of a ballistic expansion, thetime-of-flight image shows the continuum momentum oc-cupation N ( k ) , defined by N ( k ) = h ˜Ψ † ( k ) ˜Ψ( k ) i , (45)where ˜Ψ( k ) is the (continuum) momentum-space annihi-lation operator. The real-space operator Ψ( r ) can, af-ter projection to the lowest Bloch band, be expressed interms of the lattice operator b j using the Wannier func-tion W j ( r ) , Ψ( r ) = X j W j ( r ) b j . (46)In the presence of a magnetic field, this expression cannotgenerally be written as a convolution.Combing Eqs. (45) and (46) and using Eq. (11) to ex-press b j in terms of b k gives N ( k ) = Z k , k ∈ B L F ∗ ( k , k ) h b † k b k i F ( k , k ) (47)(where, for brevity, the standard integration measure forboth integrals has been omitted). The kernel of this dou-ble integral transform, analogous to a matrix similaritytransformation, is given by F ( k , k ′ ) = X j e i k ′ · x j Z d r e − i k · r W j ( r ) , (48)where r is integrated over all space. (The Wannier func-tion W j restricts the integral to the neighborhood of thetwo-dimensional plane.)The kernel F depends on experimental details, includ-ing the optical lattice parameters and the effective gaugepotential, while theoretical analysis based on the Bose-Hubbard model leads to predictions for the correlationfunction h b † k b k i . We will first outline the form of F appropriate to experiments using rotation and Raman-induced gauge fields, before giving our results for thecorrelation function based on the Bogoliubov theory.In the absence of a magnetic field, the Wannier func-tion at site j is a function only of r − x j , and so can bewritten in the form W (0) j ( r ) = w (0) ( r − x j ) . Shifting theintegration variable r in Eq. (48) allows the sum to beevaluated, giving F (0) ( k , k ′ ) = (2 π ) δ ([ k − k ′ ] B L ) ˜ w (0) ( k ) , (49)where ˜ w (0) is the Fourier transform of w (0) . This leadsto the simple result N (0) ( k ) = | ˜ w (0) ( k ) | D b † [ k ] B L b [ k ] B L E , (50) so the time-of-flight image gives the lattice-momentumdistribution, with an overall envelope given by the Wan-nier function. With nonzero gauge potential A , one can instead ex-press the Wannier function as W j ( r ) = w ( r − x j ) exp " i Z rx j d r ′ · A ( r ′ ) , (51)where the integral is taken along a straight-line path, asin Eq. (3). In this case, the kernel F is no longer diagonalin k and k ′ , and will depend on the appropriate choiceof gauge.In the Landau gauge (with A parallel to ˆ y ), one canwrite A L ( r ) = B × r x ˆ x , and the kernel is given by F L ( k , k ′ ) = Z d r e − i k · r w ( r )e i | B | r x r y × (2 π ) δ ([ k − k ′ + B × r y ˆ y ] B L ) . (52)In the symmetric gauge, A S ( r ) = − B × r , leading to F S ( k , k ′ ) = Z d r e − i k · r w ( r ) × (2 π ) δ ([ k − k ′ + 14 B × r ] B L ) . (53)It should be noted that, in both cases, w ( r ) differs fromthe Wannier function w (0) ( r ) in the absence of a mag-netic field.For our purposes, it is sufficient to note that the ker-nels both give contributions to N ( k ) from a range oflattice momenta [ k ] B L + δ k . The scale is determined by | δ k | . | B | d = 2 παd , where d is the characteristic size ofthe Wannier function w . This point-spreading effect canbe understood as resulting from the position-dependentimpulse A ( r ) imparted to the atoms when the gauge po-tential is switched off.We now discuss the form of the correlation function h b † k b k i , which also depends on the gauge. As in previoussections, we will focus on the Landau gauge, appropriatefor experiments with Raman-induced gauge potentials.The symmetry under T qx and T qy implies that h b † k b k i vanishes unless [ k − k ] B N = . We therefore define (2 π ) δ ( k − k ′ ) N nℓ,n ′ ℓ ′ ( k ) = h b † k + n X + ℓ Y b k ′ + n ′ X + ℓ ′ Y i ,(54)for k , k ′ ∈ B N . (The formally infinite factor when k = k ′ corresponds physically to system volume.)The dominant contribution to N nℓ,n ′ ℓ ′ ( k ) is a delta-function peak at k = , coming from the first term inEq. (1), N nℓ,n ′ ℓ ′ ( k ) = (2 π ) δ ( k ) X γγ ′ A ∗ ℓγ ψ ∗ γn ( k ) A ℓ ′ γ ′ ψ γ ′ n ′ ( k ) .(55)The corresponding peaks in h b † k b k i , at momenta suchthat [ k ] B N = [ k ] B N = , are separated by multiples of4 X and Y . If they are to be resolved in time-of-flight im-ages, we require that their separation, π/q , be greaterthan the point-spread παd of the kernel F . The con-dition is then simply that the Wannier function w ( r ) bewell localized compared to the lattice spacing.If this condition is satisfied, then the time-of-flight in-tensity N ( k ) consists of sharp peaks near each of thereciprocal lattice vectors n X + ℓ Y of the reduced Bril-louin zone B N . These extra Bragg peaks in fact resultfrom two separate physical effects.The first is the enlargement of the unit cell of theHamiltonian to q × sites, as a result of the phasesappearing the hopping term H t . States with momen-tum differing by X are therefore mixed at the single-particle level, leading to additional Bragg peaks at mo-menta π ˆ x n/q even in the absence of interactions. Theobservation of such peaks in an experiment is a clear signthat the flux per plaquette is at (or sufficiently close to)a rational value.The second effect occurs only in the presence of inter-actions and is due to the spontaneous breaking of spatialsymmetry in the superfluid. As discussed in Section III,the condensate contains contributions from the q degen-erate minima of the single-particle dispersion, and there-fore enlarges the unit cell to q × q sites. Peaks at momenta π ˆ y ℓ/q are clear indications of the formation of such anordered state.Importantly, the kernel for the Landau gauge, F L ( k , k ′ ) in Eq. (52), does not change the y componentof the momentum, and so the point-spreading effect isentirely in the x direction. The second class of Braggpeaks, resulting from interaction effects, are therefore notaffected, increasing the likelihood that they can be ob-served in experiment. Note that this separation does notapply in the symmetric gauge, appropriate for the case ofrotation, and furthermore that the spacing of the Braggpeaks is reduced, as a result of the q × q unit cell of H t ,making their observation considerably more challenging.Finally, it should be noted that in many cases, includ-ing α = and , the condensate has equal amplitude(but not phase) for all values of ℓ . Bragg peaks withthe same value of n therefore have the same intensity,apart from the envelope coming from the on-site Wannierwavefunction. This is in contrast to the case of strictlyvanishing interactions, when any distribution of particlesbetween the q minima of the single-particle dispersionhas equal probability. B. Spectroscopy
Developments in spectroscopic measurements forultracold atomic systems have allowed experi-mental access to dynamic correlation functions withinthese systems. We consider two such techniques,Bragg spectroscopy and lattice-modulation spectroscopy, and describe the information regard-ing the quasiparticle spectrum that can be determinedfrom both.Bragg spectroscopy involves applying a weak pe-riodic perturbation of the form cos( K · x − Ω t ) to the sys-tem, using two laser beams at an angle and with frequen-cies differing by Ω . One then measures, usually throughtime-of-flight imaging, the total momentum or energy im-parted to the system. The response is given by the dy-namical structure factor S ( K , Ω) , the density correlationfunction in momentum and frequency space.Lattice-modulation spectroscopy involves oscillat-ing the lattice depth at frequency Ω ; using time-of-flightimaging to determine the imparted energy then gives S ( , Ω) . It is also possible to measure S ( K , Ω) at cer-tain high-symmetry points K in the lattice Brillouin zone B L by the application of lattices with enlarged periods.(Each point in B L corresponds to a point in B N in away that depends on q .) Even if only the point K = is accessible, the presence of multiple Hofstadter bandsshould be clear, and the splitting of the Goldstone modefrom the rest of the first band is also measurable.In either case, the coupling to the perturbation canbe expressed in terms of the momentum-space densityoperator, ρ ( K ) = Z k ∈ B L d k (2 π ) b † k − K b k = X j e − i x j · K n j , (56)which, being a function only of n j , is gauge-invariant.The dynamic structure factor is given in spectral repre-sentation by S ( K , Ω) = 1 Z X Ψ , Ψ e − E Ψ1 /T δ ( E Ψ − E Ψ − Ω) × |h Ψ | ρ ( K ) | Ψ i| , (57)where Z = P Ψ e − E Ψ /T is the partition function, and | Ψ , i are eigenstates of the Hamiltonian H with energy E Ψ , .While S ( K , Ω) is given by a four-point correlationfunction, it can be factorized into two-point functionswithin the quadratic Bogoliubov theory. In the con-densed phase, and assuming depletion is not too large,the dominant contribution to the integral in Eq. (56) infact comes from the points where either [ k ] B N = or [ k − K ] B N = . (For [ K ] B N = these cases coincide,and there is an extra term in ρ ( K ) which, however, con-tributes only at Ω = 0 .) The structure factor is thereforegiven by a two-point correlation function multiplied bythe condensate density.Within this approximation, the density operator canbe expanded in terms of the operators d k ζ and d † k ζ :5 ρ ([ K + N X + L Y ] B L ) = X ζ n r ζNL ( K ) d K ζ + ω − y L [ r ζNL ( K )] ∗ d †− K ζ o , (58)where r ζNL ( k ) = q − X n,l =0 X γγ ′ h ω − nℓ +( n − N )( ℓ − L ) A ∗ [ ℓ − L ] q γ ψ ∗ γ [ n − N ] q ( ) ψ γ ′ n ( k ) X ζ k ℓγ ′ + ω nℓ − ( n + N )( ℓ + L ) A [ ℓ + L ] q γ ψ γ [ n + N ] q ( ) ψ ∗ γ ′ n ( − k ) Y ζ k ℓγ ′ i . (59)A subdominant third term involving two d k ζ operators has been dropped from Eq. (58).The dynamic structure factor defined in Eq. (57) is calculated using the eigenstates of H , which, at the level ofthe quadratic approximation of Section IV, are eigenstates of the occupation numbers of each Bogoliubov mode. Thematrix elements of Eq. (58) between any pair of states can be expressed in terms of these occupation numbers, giving S ([ K + N X + L Y ] B L , Ω) = X ζ | r ζNL ( K ) | n δ (Ω + ξ K ζ ) n B ( ξ K ζ ) + δ (Ω − ξ K ζ )[1 + n B ( ξ K ζ )] o . (60)The structure factor at frequency Ω therefore has reso-nances at each quasiparticle mode ζ , allowing the quasi-particle spectrum to be measured directly. VI. DISCUSSION
We have studied the effect of a synthetic magnetic fieldon the superfluid phase of bosons in a lattice. Our the-oretical approach is based on Bogoliubov theory, whichdetermines the condensate configuration and allows inter-actions to be taken into account within an expansion interms of fluctuations. We predict broken spatial symme-try in the condensed phase, leading to qualitative changescompared to the Hofstadter spectrum for noninteractingparticles.This analysis leads to several clear predictions thatshould be testable in experiment. The density modu-lations in the superfluid phase, illustrated in Figure 4,may be directly measurable using recently developed real-space imaging techniques.
Our order-of-magnitudeestimate for the extent of superfluidity given in Sec-tion IV B, which is in agreement with independent the-oretical approaches, can be tested in experimentsanalogous to those performed in the absence of a mag-netic field. Predictions for spectroscopic measurementshave been detailed in Section V B.Our approach also provides predictions for time-of-flight imaging, the most well-established technique incold-atom experiments. As described in Section V A, wepredict extra Bragg peaks due to the spatial symmetrybreaking. In experiments using Raman-induced gaugefields, these result from two distinct physical effects. Thegradient in the applied synthetic vector potential (dueto a gradient in the physical magnetic field in the ex-periments of Lin et al. ) breaks translation symmetry explicitly, leading to an extra set of Bragg peaks in thedirection of the gradient. By contrast, symmetry un-der translation in the perpendicular direction is brokenspontaneously when the bosons condense, and this leadsto further Bragg peaks, in the direction of propagation ofthe applied Raman lasers ( ˆ x in our convention, but ˆ y inthe experiments). The appearance of these latter peaksis therefore a clear signature of many-body effects.Among the approximations made in the present work isthe assumption that thermal equilibrium can be reachedon the time scale of the experiments. Previous studies ofclosely related systems have shown that the process ofvortex formation can exhibit hysteresis, and experimentswith effective gauge potentials exhibit a considerable de-pendence of the vortex density on hold times. In themodel considered here, two-body scattering is sufficientto populate modes of nonzero ℓ and hence generate non-trivial spatial structures, but further work is required toprovide quantitative estimates of the rate for these pro-cesses.We have also neglected the influence of higher latticebands and hopping between pairs of sites other than near-est neighbors. Neither is expected to have qualitativeeffects on our conclusions, as long as the magnetic sym-metries described in Section II A are preserved. (This iscertainly the case with a synthetic magnetic field due torotation or Raman lasers, but not necessarily so whenhopping phases are induced by other methods. ) Asalready noted in Section I B, weak interactions betweenbosons on different sites will also have only quantitativeeffects.As discussed in Section III, the broken spatial sym-metry implies the existence of multiple degenerate con-figurations in the superfluid phase. This allows for thepossible formation of real-space domains, especially onshorter time scales, upon which the effect of the external6trapping potential is likely to be important.Besides the simplifications inherent in our startingmodel, our analysis has made the approximation of trun-cating the Bogoliubov expansion at quadratic order,neglecting interactions between quasiparticles. Conse-quences of these interactions include finite quasiparticlelifetimes and also the possibility of spectrum terminationat the point where decay into the two-particle continuumis allowed by kinematics. These are likely to have impli-cations for spectroscopy experiments; an understandingof these is left for future work. Acknowledgments
We thank Ian Spielman, Trey Porto, Chris Foot, RossWilliams, and Sarah Al-Assam for helpful discussions.This work is supported by JQI-NSF-PFC, ARO-DARPA-OLE, and Atomtronics-ARO-MURI.
Appendix A: Bogoliubov transformation
In this Appendix, we will show that to diagonalize H (2) , given in Eq. (37), one must find the eigenvaluesand -vectors of the matrix η M k for each k ∈ B N , where η ℓγ,ℓ ′ γ ′ = δ ℓℓ ′ δ γγ ′ − ! , (A1)and the matrix product is taken treating both η and M k as q × q matrices. While η M k is not hermitian, itcan be shown that, since M k is nonnegative-definite,the eigenvalues of η M k are all real. Furthermore, for k = , when M k is positive -definite, the eigenvalues areall nonzero and come in pairs with equal magnitude andopposite sign. In this case, the q positive eigenvalues ξ k ζ and the corresponding eigenvectors V ζ k , V ζ k ℓγ = X ζ k ℓγ Y ζ k ℓγ ! , (A2)defined by η M k V ζ k = ξ k ζ V ζ k , describe the Bogoliubovquasiparticles.Section A 3 treats separately the special case of k = ,allowing us to develop a series expansion for the proper-ties near this point.
1. Inversion symmetry
The Bogoliubov transformation, which mixes annihila-tion operators at momentum k with creation operatorsat − k , requires the existence of an inversion symmetryin the condensed phase. Because of broken translationalsymmetry in the presence of a condensate, it is not neces-sarily the case that P , inversion about the origin x = , remains a good symmetry. Instead, define the operator P y = T y y P representing inversion about the real-spacepoint x = y ˆ y (a lattice site if y is an integer or thecenter of a bond if y is a half integer).Using Eqs. (21) and (26), P y can be shown to obey P y a k ℓγ = e − y k y e − φ γ ω − y ℓ a ( − k )[ − ℓ ] q γ P y , (A3)so the condensate is invariant under this transformationif A ℓγ obeys A ℓγ = A [ − ℓ ] q γ e − φ γ ω − y ℓ . (A4)The following assumes that the condensate has an inver-sion point, and hence there exists some value of y forwhich this relation holds.Corresponding to this inversion symmetry, and analo-gous to the matrix η , define the q × q matrix π , π ℓγ,ℓ ′ γ ′ = δ [ ℓ + ℓ ′ ] q , δ γγ ′ e − φ γ ω − y ℓ ω y ℓ ! . (A5)It is also convenient to define γ , γ ℓγ,ℓ ′ γ ′ = δ ℓℓ ′ δ γγ ′ ! , (A6)which has the effect of exchanging creation and annihila-tion operators: γα k = ( α †− k ) T . (The matrices η , π , and γ are all both hermitian and unitary, and obey ηπ = πη , πγ = γπ ∗ , and γη = − ηγ .)Under the assumption that A ℓγ obeys Eq. (A4),one can show that π M k π = M − k , and furthermore, γπ M k πγ = M ∗ k . These symmetries imply that for ev-ery eigenvector V ζ k of η M k with a positive eigenvalue ξ k ζ , there is a corresponding eigenvector W ζ k = πγ V ζ k ∗ with eigenvalue − ξ k ζ . The corresponding eigenvector of η M − k is W ζ − k = γ V ζ k ∗ , so that ξ − k ζ = ξ k ζ .
2. Nonzero momentum
For nonzero k , the matrix M k is positive-definite, andall eigenvalues of η M k are real and nonzero. The eigen-values then come in pairs of equal magnitude and op-posite sign, as claimed previously. It can furthermore beshown that one can normalize the q vectors V ζ k so that V ζ k † η V ζ ′ k = δ ζζ ′ , (A7)and hence W ζ k † η W ζ ′ k = − δ ζζ ′ and W ζ k † η V ζ ′ k = 0 .These orthonormality relations immediately lead tothe results used in Section IV A. First, they imply thatthe operators d k ζ = V ζ k † ηα k = − α †− k η W ζ − k obey thecanonical commutation relations, [ d k ζ , d † k ′ ζ ′ ] = (2 π ) δ ( k − k ′ ) δ ζζ ′ . (A8)7Second, they lead to the inverse expression α k = X ζ (cid:16) V ζ k d k ζ + W ζ k d †− k ζ (cid:17) , (A9)giving α † k M k α k = X ζ ξ k ζ (cid:16) d † k ζ d k ζ + d k ζ d † k ζ (cid:17) , (A10)from which Eq. (40) follows.
3. Near zero momentum
The invariance of the mean-field energy h underchanges of phase of A ℓγ leads to a vanishing eigenvalueof the matrix M , which is by assumption the only zeroeigenvalue (generically the case when no further contin-uous symmetries are broken). The corresponding eigen-vector is given by P ℓγ = 1 p | A | i A ℓγ − i A ∗ ℓγ ! , (A11)and is obviously also an eigenvector of η M with zeroeigenvalue. We choose the normalization and phase of P so that P † P = 1 and P = π P = γ P ∗ .It is convenient to define the vector Q satisfying η M Q = − i ν P (A12) Q † η P = i (A13) Q † P = 0 , (A14)where ν is a constant with dimensions of energy. Thevector Q is specified uniquely by these three equations,because M has only one zero eigenvalue, and so its in-verse can be defined in the subspace orthogonal to P .One can therefore write Q = − i ν M − η P , because η P isorthogonal to P , and so is Q by Eq. (A14). Eq. (A13)simply fixes the normalization of Q , and thus the con-stant ν , which can be shown to be positive and givenby ν − = P † η M − η P . (A15)The vectors P and Q are ‘conjugate’ in the sense thatthey obey Eq. (A13) along with P † η P = Q † η Q = 0 , (A16)so that operators constructed using the vectors P and Q obey the commutation relations of momentum andposition. These operators, however, apply only preciselyat the zero-measure point k = , and so are not directlyrelevant for the properties in the thermodynamic limit.Together with the eigenvectors V ζ and W ζ corre-sponding to nonzero eigenvalues, P and Q span the q -dimensional vector space, and can therefore be used as a basis for a perturbative description of the region near k = . For small k (along ˆ x , say), the matrix M k canbe expanded as M k x ˆ x ≃ M + k x M (1) + k x M (2) , (A17)and M (1 , can be treated as perturbations. The coef-ficient matrices can themselves be calculated using per-turbation theory for the wavefunctions ψ γn ( k ) . The re-sulting expressions are analytic functions of the momen-tum k , implying the important symmetry property that P † M (1) P = Q † M (1) Q = 0 .In the following, we restrict to cases where the conden-sate is sufficiently symmetric that the phonon velocity isisotropic. This requires a symmetry of A under either R or I xy (see Appendix B), and as seen in Table I, is alwaysthe case for q ≤ . (The extension to general symmetryis straightforward.)Because of the nonhermitian nature of the matrix η M k , standard Rayleigh-Schrödinger perturbation the-ory cannot be applied directly to this problem, and theresulting expressions for the eigenvalues and -vectors arenot analytic in k . Instead, the smallest eigenvalue ξ k islinear in | k | , ξ k = c | k | + O ( | k | ) , (A18)with phonon velocity c given by c ν = P † M (2) P − X ζ | P † M (1) V ζ | ξ ζ . (A19)As described in Section A 2, away from k = alleigenvectors have nonzero eigenvalues and can be nor-malized so that V ζ k † η V ζ k = − W ζ k † η W ζ k = 1 . As k approaches , the vectors V k and W k correspondingto the smallest eigenvalue both approach P , the uniquezero-eigenvector of M . Since this eigenvector satisfies P † η P = 0 , the normalization of both V k and W k mustdiverge as k → . To leading order, one finds V k = v | k | − / P + O ( | k | / ) , (A20)with coefficient v = p ν/ (2 c ) . The omitted higher-order terms maintain the appropriate normalization: P † η V k = v | k | / c/ν + O ( | k | / ) . Appendix B: Symmetries and degeneracies
In Section II A, the magnetic symmetry group (MSG)was introduced, and the operators corresponding to var-ious elementary operations were defined. In this Ap-pendix, we will present in detail the consequences ofthese symmetries for the condensate configurations andthe quasiparticle spectrum.We denote the full group of symmetries of the Hamil-tonian H as G , which includes the translations, rota-tions, and reflections considered in Section II A. As noted8in Section III, the mean-field condensate configurationbreaks a subset of G ; the subgroup of symmetries thatare preserved will be denoted H . Examples are given inTable I, which lists the symmetries preserved in the con-figurations illustrated in Figure 4. We will discuss thevarious properties of the spectrum, including symmetriesand degeneracies, that result from H .A general spatial transformation, such as an element of G , is represented by the (anti)unitary operator S , underwhich a k ℓγ transforms as S a k ℓγ = X ℓ ′ ,γ ′ s ℓγ,ℓ ′ γ ′ ( k ) a ( k ) ℓ ′ γ ′ S . (B1)It should be recalled that reflection operators must becombined with time reversal to give symmetries of theHamiltonian, and are hence represented by antiunitaryoperators. The coefficients in Eq. (B1), which can bewritten as a q × q matrix s k , are diagonal in γ (withpossible exceptions at points where two bands touch).Preservation of the commutation relations requires that s k be a unitary matrix, including for antiunitary S .Under the general operation S , the condensate config-uration A is mapped to s − ( ) A ( ∗ ) , with complex con-jugation if S is antiunitary. (The inverse matrix arisesfrom considering transformations of states versus oper-ators.) The mean-field energy h , defined in Eq. (33),is therefore symmetric under the same transformationsof A ℓγ as the full Hamiltonian is under transformationsof a ℓγ . A given configuration that minimizes h willin general have lower symmetry, however, being invari-ant only under (a group of mappings isomorphic to) H .As argued in Section III, the superfluid therefore breaksspatial symmetries as well as the U(1) phase symmetry.As usual, the broken symmetries, comprising the sub-set G \ H , imply the existence of multiple degenerate con-figurations. For example, h is invariant under A ℓγ → A [ ℓ − q γ and under A ℓγ → ω − ℓ A ℓγ , corresponding to T x and T y respectively [see Eqs. (21) and (22)]. These twomappings do not commute, and so there is no config-uration A that preserves both symmetries. This leadsto the conclusion that the number of distinct configu-rations that minimize h is always a multiple of q (bythe same argument that implies the degeneracy of thesingle-particle states with different ℓ ).The consequences of the symmetries H for the quasi-particle spectrum can be determined by considering thetransformations of the quadratic matrix M k . The single-particle contribution to M k , given by the first term inEq. (38), commutes with any operation S in the full sym-metry group G , while the second term results from inter-actions with the condensate, and so is symmetric onlyunder the elements of H .These preserved symmetries imply constraints on thematrix M k . In particular, defining s as in Eq. (B1)and taking S to be (anti)unitary, one can show that if s − A ( ∗ ) = σ A , then Σ k M ( ∗ ) k = M k Σ k , where Σ k = σ ∗ s † k σ s T − k ! . (B2)In the presence of a condensate that is symmetric under atransformation S , the quasiparticle energies are thereforeequal at k and k (including in the case where S is antiu-nitary). For antiunitary operators, it is convenient to usethe notation ˇ Σ for the operation of complex conjugationfollowed by multiplication by the matrix Σ .Applied to momenta near k = , this leads to con-straints on the phonon speed c , calculated in Section A 3.In all four cases listed in Table I, the symmetry is suf-ficient to have isotropic phonon speed, as assumed inEq. (A18). For α = and , the relevant symmetry isthe (antiunitary) reflection I xy , while for α = and ,it is the rotation R . (For α = , there is also symmetryunder rotation by π about a plaquette center, T x R .)In all cases, the spectrum is of course symmetric onlyunder, at most, the fourfold rotation symmetry of thesquare lattice. The symmetry under continuous rotationsof k in Eq. (A18) is a simple example of an emergent low-energy symmetry.In many cases (for example α = , , and ; see Fig-ure 4 and Table I), the condensate configuration pre-serves a nontrivial translation symmetry. In particular,suppose translation T R by a displacement R is unbroken,where R is not a lattice vector of the enlarged unit cell.(Because T x and T y do not commute, one must specifythe path to define T R precisely.) Since translation doesnot change k , one can define an operator using Eq. (B2)that commutes with M k . For each k , the modes ζ cantherefore be labeled according to their eigenvalue under Σ k . Modes with different eigenvalues are allowed to have(unavoided) crossings, as visible for example in Figure 6.
1. High-symmetry points
Points in momentum space separated by the reciprocallattice vectors π ˆ x and π ˆ y are physically equivalent, sothe full Brillouin zone B L has the topology of a torus.The same applies, with some modifications, to the doublyreduced Brillouin zone B N .With momentum shift operators defined by K x k = [ k + πq ˆ x ] B L and K y k = [ k + πq ˆ y ] B L , the matrix H nn ′ ( k ) ,defined in Eq. (17), obeys H nn ′ ( K x k ) = H [ n +¯ p ] q [ n ′ +¯ p ] q ( k ) (B3) H nn ′ ( K y k ) = ω − ¯ p ( n − n ′ ) H nn ′ ( k ) , (B4)so one can extend the definition of ψ γn ( k ) beyond B N by ψ γn ( K x k ) = ψ γ [ n +¯ p ] q ( k )e i θ xγ ( k ) (B5) ψ γn ( K y k ) = ω − ¯ pn ψ γn ( k )e i θ yγ ( k ) . (B6)9The phases θ x,yγ ( k ) , corresponding to the flux threadedthrough the holes of the torus, are arbitrary apart fromconstraints due to symmetry. These can be found by con-sidering the commutation relations of the operators K x,y with each other and with the symmetry operators, andby requiring that ψ γn ( k ) be continuous (apart from pos-sibly at degeneracy points). The definitions in Eqs. (B5)and (B6) are particularly useful on the boundary of B N ,where the number of constraints on θ x,yγ ( k ) is larger, andespecially at the high-symmetry points at the corner (X)and edge-center (M) of B N .The relations between the eigenvectors at k and K µ k lead to corresponding relations for M k , given by K µ k M k = M K µ k K µ k , where K xℓγ,ℓ ′ γ ′ ( k ) = δ ℓℓ ′ δ γγ ′ ω − ¯ pℓ e − i θ xγ ( k )
00 e i θ xγ ( k ) ! (B7) K yℓγ,ℓ ′ γ ′ ( k ) = δ γγ ′ e − i θ yγ ( k ) δ ℓ ′ , [ ℓ +¯ p ] q
00 e i θ yγ ( k ) δ ℓ ′ , [ ℓ − ¯ p ] q ! .(B8)These immediately imply that the mode energies areequal at points on opposite sides of B N .
2. Kramers degeneracy
For α = , there is a twofold degeneracy at every pointalong the line from M to X, as can be seen in Figure 5,and by symmetry at every point on the edge of B N . (Thesame degeneracy also occurs for α = .)This is in fact a Kramers degeneracy, and is a conse-quence of the symmetry under the glide reflection T x I y (for α = ; the corresponding symmetry for α = is T y T x I y ). Its action in momentum space is to map k to I x k , so that a point k = πq ˆ x + k y ˆ y , on the line from Mto X, is mapped to − πq ˆ x + k y ˆ y , on the opposite side of B N . The quadratic matrix M k at such a point thereforeobeys [ˇ Γ k , η M k ] = , where Γ k = K x I x k Σ T x I y k , (B9)with Σ T x I y k defined in Eq. (B2) and K x k in Eq. (B7).The antiunitary operator ˇ Γ k has the property ˇ Γ k = Γ k Γ ∗ k = − , (B10)which implies, by Kramers’ theorem, that all eigenval-ues of η M k are twofold degenerate. Note that, in thecase α = , this Kramers degeneracy exists despite theexplicit breaking of time-reversal symmetry by the ap-plied magnetic field. (For α = , ω = − is real, and theHamiltonian preserves time-reversal symmetry. The con-densate configuration nonetheless breaks this symmetry,as shown in Figure 4.) Appendix C: Analytics for α = The simplest nontrivial case is α = , and it is thenpossible to perform many of the calculations analytically.In this case, the matrix H ( k ) defined in Section II B isgiven by H ( k x ˆ x + k y ˆ y ) = − t cos k x cos k y cos k y − cos k x ! , (C1)with eigenvalues ǫ k = ± t q cos k x + cos k y . (C2)Note that the spectrum has a Dirac cone at the corner of B N , where | k x | = | k y | = π . At k = , the eigenvectorsare, for bands γ ∈ { , } , ψ γ ( ) = 1 q − γ √ − γ √ ! . (C3)The mean-field condensate configuration can be deter-mined by calculating h , given in Eq. (33), and minimiz-ing with respect to A ℓγ at fixed average density. In thiscase, however, the appropriate configuration is more eas-ily found by inspection. With the choice A ℓ, = i ℓ p ρ/ and A ℓ, = 0 , direct calculation shows that the real-spacewavefunction, given by Eq. (34), has uniform magnitude, h b j i = √ ρ exp (cid:26) i( − y j [( − x j π − π (cid:27) . (C4)This configuration, in which the condensate is restrictedto the lower band, therefore has uniform density of ρ particles per lattice site. Calculation of the currents usingEq. (35) gives configurations as illustrated in Figure 4,with the current on each link having magnitude |hJ i| = √ tρ . Replacing A ℓγ by its complex conjugate gives theequivalent configuration with currents reversed on eachlink.Since these configurations have uniform density, theysimultaneously minimize both terms in Eq. (33) globally,and are therefore global minima of h , at fixed density ρ . (They minimize the kinetic energy because they con-tain contributions only from the degenerate minima ofthe lowest band, and they minimize the potential energybecause they have uniform density.) The case α = isunique in this regard, with the condensate configurationminimizing both terms simultaneously and so insensitiveto the value of the interaction strength U . For q > , thedensity modulations can be reduced by including higherbands in the condensate configuration, and the extent towhich they contribute is determined by U/t .These two configurations are in fact the only minima(up to a redundant overall phase rotation), and providean example of the general result that there is always adiscrete set of degenerate minima, whose number is a0multiple of q . Either of the translation operators T x and T y relates one of the two configurations to the other, upto an overall phase (using the transformation of A ℓγ spec-ified in Appendix B).Both configurations are symmetric under T y T x , asnoted in Table I, with the first obeying X ℓℓ ′ (cid:0) s T y T x (cid:1) − ℓℓ ′ A ℓ ′ , = i A ℓ, , (C5)where s T y T x = − ! is the matrix defined by Eq. (B1)for the transformation T y T x (at k = and restricted to γ = 1 ). One can therefore construct the matrix Σ T y T x according to Eq. (B2); it has eigenvalues ± , allowing thequasiparticle modes to be labeled as even or odd under T y T x .
1. Quasiparticle dispersion
To find the dispersion, one must construct the × matrix M k given in Eq. (38). Rather than using thesingle-particle basis labeled by k , ℓ , and γ , it is somewhateasier to find analytic results by starting in the basis of k , ℓ , and n as in Eq. (16). While the first term in Eq. (38) isnot diagonal in this basis, the second has a considerablysimpler expression, as a result of the simple form of theon-site interaction in momentum space.After transforming to the basis of eigenvectors of Σ T y T x , the matrix M k splits into two × blocks, M ( ± ) k = U ρ + 2 t ( √ k x ) 2 t cos k y ± U ρ/ √ U ρ/ √ t cos k y U ρ + 2 t ( √ − cos k x ) U ρ/ √ ∓ U ρ/ √ ± U ρ/ √ U ρ/ √ U ρ + 2 t ( √ k x ) 2 t cos k y U ρ/ √ ∓ U ρ/ √ t cos k y U ρ + 2 t ( √ − cos k x ) , (C6)corresponding to the eigenvalues ± . It is then straightforward to calculate the eigenvalues of η M k as discussed inSection A 1. They come in pairs of opposite sign, and so their squares are given by the roots of a quadratic equation.The quasiparticle energies ξ k are finally given by ξ k = 8 t + 4 √ tU ρ + ǫ k ± q (32 t + 16 √ tU ρ + 2 U ρ ) ǫ k ± t U ρ cos k x cos k y , (C7)where the two choices of ± are independent, giving the q = 4 modes of the interacting dispersion. Taking thefirst sign as − and the second as + gives the Goldstonemode, which has ξ k = | k | p √ ρU t + O ( | k | ) near k = .Using Eq. (C7), one can confirm the twofold degener-acy along the line from M to X established in Section B 2.For these points, cos k x = 0 and the second choice of ± is redundant.
2. Amplitude-phase description of gapless mode
The small- | k | dispersion of the Goldstone mode can bederived by more elegant means if one restricts to long-wavelength fluctuations of the condensate configuration.Gradual variations of the real-space wavefunction can beparametrized by writing b j = h b j i e i ϑ j r ̺ j ρ , (C8)where ϑ j and ̺ j describe deviations in the phase andamplitude respectively, and have canonical commutationrelations [ ϑ i , ̺ j ] = i δ ij . We now rewrite the Hamiltonian H in terms of thesenew degrees of freedom. Assuming ̺ j ≪ ρ and that both ̺ j and ϑ j vary only over distances large compared to thelattice scale, one can expand to give H = h + tρ √ X h ij i ( ϑ i − ϑ j ) + U X j ̺ j + · · · . (C9)Note that the frustration in H t implies that the kineticenergy of each link is not separately minimized in themean-field configuration. Each link i → j therefore con-tributes a term linear in ϑ i − ϑ j , but their sum vanishes,since the mean-field configuration is a minimum of thetotal kinetic energy.Writing this equation in terms of the Fourier compo-nents of ϑ j and ̺ j , we obtain H = h + Z d k (2 π ) (cid:18) tρ √ | k | | ϑ k | + U | ̺ k | (cid:19) + · · · ,(C10)where the integral is restricted to small | k | by the as-sumption of slowly varying fluctuations. This takes theform of a harmonic oscillator for each momentum, so thedispersion is ξ k = 2 q tρ √ | k | × U , in agreement with the1result given in the previous section. M. Tinkham,
Introduction to Superconductivity , McGraw-Hill, New York (1996). R. J. Donnelly,
Quantized Vortices in Helium II , Cam-bridge University Press, Cambridge (1991). N. R. Cooper, Adv. Phys. , 539 (2008). S. Tung, V. Schweikhard, and E. A. Cornell, Phys. Rev.Lett. , 240402 (2006). J. W. Reijnders and R. A. Duine, Phys. Rev. Lett. ,060401 (2004); Phys. Rev. A , 063607 (2005). R. A. Williams, S. Al-Assam, and C. J. Foot, Phys. Rev.Lett. , 050404 (2010). M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S.Fisher, Phys. Rev. B , 546 (1989). D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P.Zoller, Phys. Rev. Lett. , 3108 (1998). I. B. Spielman, Phys. Rev. A , 063613 (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.-J. Lin, R. L. Compton, K. Jiménez-García, J. V. Porto,and I. B. Spielman, Nature , 628 (2009). R. Dum and M. Olshanii, Phys. Rev. Lett. , 1788 (1996). D. Jaksch and P. Zoller, New J. Phys. , 56 (2003). E. J. Mueller, Phys. Rev. A , 041603(R) (2004). A. S. Sørensen et al., Phys. Rev. Lett. , 086803 (2005). F. Gerbier and J. Dalibard, New J. Phys. , 033007(2010). M. Hafezi, A. S. Sørensen, E. Demler, and M. D. Lukin,Phys. Rev. A , 023613 (2007). R. N. Palmer, A. Klein, and D. Jaksch, Phys. Rev. A ,013609 (2008). G. Möller and N. R. Cooper, Phys. Rev. Lett. , 105303(2009). J. M. Luttinger, Phys. Rev. , 814 (1951). P. G. Harper, Proc. Phys. Soc., London, Sect. A , 874(1955). G. H. Wannier, Rev. Mod. Phys. , 645 (1962). J. Zak, Phys. Rev. , A1602; , A1607 (1964). D. Hofstadter, Phys. Rev. B , 2239 (1976). G. M. Obermair and H.-J. Schellnhuber, Phys. Rev. B ,5185 (1981); H.-J. Schellnhuber, G. M. Obermair, and A.Rauh, ibid. , 5191 (1981). D. S. Goldbaum and E. J. Mueller, Phys. Rev. A ,033629 (2008); Phys. Rev. A , 021602(R) (2009). S. Powell, R. Barnett, R. Sensarma, and S. Das Sarma,Phys. Rev. Lett. , 255303 (2010). A. M. Rey, K. Burnett, R. Roth, M. Edwards, C. J.Williams, and C. W. Clark, J. Phys. B: At. Mol. Opt.Phys. , 825 (2003). L.-K. Lim, C. Morais Smith, and A. Hemmerich, Phys.Rev. Lett. , 130402 (2008); Phys. Rev. A , 023404(2010). T. Ðurić and D. K. K. Lee, Phys. Rev. B , 014520 (2010). H. Zhai, R. O. Umucalılar, and M. Ö. Oktel, Phys. Rev.Lett. , 145301 (2010). S. Alexander, Phys. Rev. B , 1541 (1983). B. Pannetier, J. Chaussy, R. Rammal, and J. C. Villegier,Phys. Rev. Lett. , 1845 (1984). S. Teitel and C. Jayaprakash, Phys. Rev. Lett. , 1999(1983). T. C. Halsey, Phys. Rev. B , 5728 (1985). M. Polini, R. Fazio, A. H. MacDonald, and M. P. Tosi,Phys. Rev. Lett. , 010401 (2005). K. Kasamatsu, J. Low Temp. Phys. , 593 (2008); Phys.Rev. A , 021604(R) (2009). M. Ö. Oktel, M. Nita, and B. Tanatar, Phys. Rev. B ,045133 (2007). R. O. Umucalılar and M. Ö. Oktel, Phys. Rev. A ,055601 (2007). V. W. Scarola and S. Das Sarma, Phys. Rev. Lett. ,210403 (2007). D. S. Goldbaum and E. J. Mueller, Phys. Rev. A ,063625 (2009). E. Lundh, Europhys. Lett. , 10007 (2008). S. Sinha and K. Sengupta, arXiv:1003.0258v1 (unpub-lished). R. Moessner, Can. J. Phys. , 1283 (2001). T. D. Stanescu, V. Galitski, J. Y. Vaishnav, C. W. Clark,and S. D. Sarma, Phys. Rev. A , 053639 (2009); T. D.Stanescu, V. Galitski, and S. Das Sarma, Phys. Rev. A ,013608 (2010). K. Saha, K. Sengupta, and K. Ray, arXiv:1005.4476v1 (un-published). N. N. Bogoliubov, J. Phys. (USSR) , 23 (1947). A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinski,
Methods of Quantum Field Theory in Statistical Physics ,Dover, New York (1963). X.-G. Wen, Phys. Rev. B , 165113 (2002). L. Balents, L. Bartosch, A. Burkov, S. Sachdev, and K.Sengupta, Phys. Rev. B , 144508 (2005); Prog. Theor.Phys. Suppl. , 314 (2005). M. Kohmoto, Phys. Rev. B , 11943 (1989). W. S. Bakr, J. I. Gillen, A. Peng, S. Fölling, and M.Greiner, Nature , 74 (2009); W. S. Bakr, A. Peng, M.E. Tai, R. Ma, J. Simon, J. I. Gillen, S. Fölling, L. Pollet,M. Greiner, arXiv:1006.0754v1 (unpublished). J. F. Sherson, C. Weitenberg, M. Endres, M. Cheneau, I.Bloch, and S. Kuhr, Nature advance online publication,doi:10.1038/nature09378 (2010). J.-P. Blaizot and G. Ripka,
Quantum theory of finite sys-tems , MIT Press, Cambridge, Mass. (1986). N. D. Mermin and H. Wagner, Phys. Rev. Lett. , 1133(1966). P. C. Hohenberg, Phys. Rev. , 383 (1966). V. Bagnato and D. Kleppner, Phys. Rev. A , 7439 (1991) K. Jiménez-García, R. L. Compton, Y.-J. Lin, W.D. Phillips, J. V. Porto, and I. B. Spielman,arXiv:1003.1541v1 (unpublished). B. Capogrosso-Sansone et al., Phys. Rev. A , 015602(2008). E. Toth, A. M. Rey, and P. B. Blakie, Phys. Rev. A ,013627 (2008). One can easily check that this gives the correct result fora classical particle that starts at position r with zero ki-netic energy. It has zero dynamical momentum, and hencecanonical momentum A ( r ) . When the Raman beams are suddenly switched off, its canonical momentum is pre-served, so its final momentum (dynamical or canonical)is A ( r ) . This momentum can be viewed as the result ofan induced electric field − ˙ A acting during the period thatthe gauge potential is switched off. J. Stenger, S. Inouye, A. P. Chikkatur, D. M. Stamper-Kurn, D. E. Pritchard, and W. Ketterle, Phys. Rev. Lett. , 4569 (1999). D. M. Stamper-Kurn, A. P. Chikkatur, A. Görlitz, S. In-ouye, S. Gupta, D. E. Pritchard, and W. Ketterle, Phys.Rev. Lett. , 2876 (1999). C. Kollath, M. Köhl, and T. Giamarchi, Phys. Rev. A ,063602 (2007). J. T. Stewart, J. P. Gaebler, and D. S. Jin, Nature , 744 (2008). R. Sensarma, D. Pekker, M. D. Lukin, and E. Demler,Phys. Rev. Lett. , 035303 (2009). N. Strohmaier, D. Greif, R. Jördens, L. Tarruell, H. Moritz,T. Esslinger, R. Sensarma, D. Pekker, E. Altman, and E.Demler, Phys. Rev. Lett. , 080401 (2010). D. Clément, N. Fabbri, L. Fallani, C. Fort, and M. Ingus-cio, J. Low Temp. Phys. , 5 (2010). A. M. Rey, P. B. Blkier, G. Pupillo, C. J. Williams, andC. W. Clark, Phys. Rev. A , 023407 (2005). M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, andI. Bloch, Nature415