RAPID COMMUNICATIONSPHYSICAL REVIEW A 96, 061603(R) (2017)

Quantum phases of two-component bosons with spin-orbit coupling in optical lattices
Daisuke Yamamoto,1 I. B. Spielman,2 and C. A. R. Sá de Melo2,312Department of Physics and Mathematics, Aoyama-Gakuin University, Sagamihara, Kanagawa 252-5258, Japan
Joint Quantum Institute, National Institute of Standards and Technology, and University of Maryland, Gaithersburg, Maryland 20899, USA
3
School of Physics, Georgia Institute of Technology, Atlanta, Georgia 30332, USA(Received 31 August 2017; published 26 December 2017)
Ultracold bosons in optical lattices are one of the few systems where bosonic matter is known to exhibit
strong correlations. Here we push the frontier of our understanding of interacting bosons in optical lattices
by adding synthetic spin-orbit coupling, and show that new kinds of density and chiral orders develop.
The competition between the optical lattice period and the spin-orbit coupling length—which can be madecomparable in experiments—along with the spin hybridization induced by a transverse field (i.e., Rabi coupling)
and interparticle interactions create a rich variety of quantum phases including uniform, nonuniform, andphase-separated superfluids, as well as Mott insulators. The spontaneous symmetry-breaking phenomena at
the transitions between them are explained by a two-order-parameter Ginzburg-Landau model with multiparticle
umklapp processes. Finally, in order to characterize each phase, we calculated their experimentally measurable
crystal momentum distributions.
DOI: 10.1103/PhysRevA.96.061603

The physics of spin-orbit coupling (SOC), which links thespin and momentum degrees of freedom in quantum particles,
is ubiquitous in nature, ranging from the microscopic worldof atoms, such as hydrogen, to macroscopic solid materials,
such as semiconductors. Recently, the effects of SOC havebeen explored in condensed-matter physics in connection
with topological insulators [1], as well as with topologicalsuperconductors [2], and superconductors without inversion
symmetry [3]. In these naturally occurring systems, it isvery difficult to control the magnitude of SOC and yet moredifficult to study correlated bosons. However, it is now possible
to create controllable artificial SOC for trapped ultracold
fermionic and bosonic atoms [4–9], the physics of which
was recently analyzed theoretically in the continuum limit[4,10–13] and one-dimensional (1D) optical lattices [14,15]. In
these cases weak-coupling effects were considered; however,
one of the emerging frontiers in this broad area of physics is
the interplay of the spin-orbit and lattice characteristic lengths,
which can be made comparable in optical lattice systems,
where additional contributions from a Zeeman field and strong
local interactions play an important role [16].
In this Rapid Communication, we first obtain the groundstate phase diagrams for two-component (↑,↓) bosons in
the presence of artificial SOC, an effective Zeeman field(created from Rabi coupling and detuning), and local interactions. With zero detuning, we identify four phases: uniform,
nonuniform, and phase-separated superfluids, along with Mott
insulating phases, depending on interactions. Secondly, wedevelop a Ginzburg-Landau theory for further characterizingthese phases. Lastly, we calculate their crystal momentum
distributions, which can be compared with experiments.
To describe the quantum phases of two-component bosons
with SOC, we begin by introducing the independent particle
Hamiltonian

Ĥ0 =

 †
(b̂k↑
k

†

b̂k↓ )k↑ − μh̄/22469-9926/2017/96(6)/061603(5)h̄/2k↓ − μ

 
b̂k↑
b̂k↓

(1)

in momentum space. Here, ks = −2t[cos(kx + skT ) +cos ky + cos kz ] + s h̄δ/2 for a three-dimensional (3D) opticallattice and kT = (kT ,0,0) is the SOC momentum. The lengthscale 2π/kT is of the order of the optical lattice spacing a,†chosen to be one. The operator b̂ks describes a creation ofs ∈ {↑,↓} ≡ {+,−} boson with momentum k. The chemical
potential μ tunes the average particle density ρ = ρ↑ + ρ↓ ≡†ks b̂ ks b̂ ks /M with M being the number of lattice sites.
In cold-atom experiments, the effective Zeeman energy  · F̂with  = (,0,δ) and F̂ being the total angular momentum
operator for spin-1/2 has two parts: spin flips through the
Rabi frequency  and a Zeeman shift via the detuning δ. The
Hamiltonian above can be engineered in the laboratory either
through Raman processes [4,5,17] or via radio-frequency chips[18,19].
The diagonalization of Ĥ0 gives two energy branches
Ek± = (k↑ + k↓ − 2μ ± (k↑ − k↓ )2 + (h̄)2 )/2.
For δ = 0 and small h̄/t, the lower branch Ek− has twodegenerate minima at kx ≈ ±kT and ky = kz = 0. The two
minima approach as  is increased, and eventually theycollapse into a single minimum at k = 0 when h̄/t 
4 sin kT tan kT . This double-minimum structure, the introduction of a new length scale 1/kT , and the interactions between
particlesĤint =1 † †Uss b̂ks b̂k+qs b̂k−qs b̂ks
2
ss

(2)

kq

provide additional contributions that are absent in the standard
spinless Bose-Hubbard system [20]. In this work, we explore
the special case where the same spin repulsions Uss are nearly
identical (U↑↑ ≈ U↓↓ = U ), but the opposite spin repulsion isdifferent from U , that is, U  U↑↓ = U↓↑  0. For instance,in the case of a mixture of the mF = 0 (↓) and mF = −1 (↑)
states from the F = 1 manifold of 87 Rb, these repulsions arenearly identical (U↑↑ ≈ U↓↓ ≈ U↑↓ ) [21].

061603-1©2017 American Physical Society

RAPID COMMUNICATIONSYAMAMOTO, SPIELMAN, AND SÁ DE MELOPHYSICAL REVIEW A 96, 061603(R) (2017)

We first investigate the regime of weak repulsive interactions (U √tρ). The bosonic fields b̂ks √can be written
as b̂ks = q Mψqs δk=(q,0,0) + âks , where Mψqs and âks
describe the Bose-Einstein condensate (BEC) with momentum
k = (q,0,0) and the residual bosons outside the condensate,

E0
∗
(ψq↑
=
M
q

∗
ψq↓
)


k↑ − μh̄/2respectively. Considering the single and double minima feafor
tures of Ek− within the first Brillouin zone, we allow 
multiple BECs with different momenta and take the sum q
to be over the set of possible momenta {q} along the (kx ,0,0)
direction. The energy per site of the condensates is   

U
h̄/2
ψq↑
∗
∗
∗
∗
ψq1 s ψq2 s ψq3 s ψq4 s + U↑↓ ψq1 ↑ ψq2 ↓ ψq3 ↓ ψq4 ↑ ,+k↓ − μ ψq↓
2
{qi }

s

(3)

where the sum {qi } is over momenta qi satisfying momentum
conservation q1 + q2 = q3 + q4 [mod 2π ].
After minimization of Eq. (3) with respect to ψqs and {q},we find four different ground states as shown in Fig. 1(a)for the weak-coupling regime with parameters U = t/ρ,U↑↓ = 0.9U , and kT = 0.2π . In the superfluid phases (SF± ),the set of BEC momenta {q} consists of a single value (q̄ > 0
in SF+ and −q̄ < 0 in SF− ) since the detuning δ tilts the
single-particle spectrum and lifts the degeneracy of the double
minima in Ek− . In these “single-q” states, the particle density
is uniform, while the phase of the condensate spatially varies
with pitch vector (±q̄,0,0). In the striped superfluid (ST) phasefor relatively small h̄/t, a BEC is formed with two different
momenta −q̄1 and q̄2 due to a double-minimum dispersion
in Ek− . The interference of these two momenta leads to anonuniform density profile along the x direction, resulting
in a stripe pattern. Moreover, the scattering process under
momentum conservation q1 + q2 = q3 + q4 with q3 = q4 =−q̄1 and q2 = q̄2 (or vice versa) gives rise to a higher harmoniccomponent with q1 = −2q̄1 − q̄2 (or q1 = q̄1 + 2q̄2 ). Similar
processes generate higher harmonics with interval q̄1 + q̄2 ,thus making the set {q} have a large number of different
momenta −q̄1 + n(q̄1 + q̄2 ), where n is an integer. Higher
harmonic generation is argued not to be qualitatively importantin the continuum limit [13] or in weak 1D optical lattices [14],but in our case of deep 3D optical lattices this mechanism plays
a crucial role in the stabilization of stripe phases commensuratewith the underlying optical lattice, as discussed below.
When h̄/t is large, the SF+ and SF− phases are continuously connected at δ = 0 through the conventional superfluid(SF0 ) with zero-momentum BEC. However, when h̄/t has
intermediate values, a direct first-order transition from SF+
to SF− takes place, and thus the spin population differenceρ↓ − ρ↑ exhibits a sudden jump from positive to negative.
Therefore, in the experimental situation where ρ↑ = ρ↓ , thesystem is unstable against spatial phase separation (PS) of
spin-down-rich SF+ and spin-up-rich SF− states.
The quadratic part of the Hamiltonian in terms of âks , ĤB = † (2)k â k H k â k , is a generalized Bogoliubov Hamiltonian [22]and includes quantum fluctuations outside the condensate
perturbatively. We diagonalize ĤB numerically [22,23], andobtain the spectrum of elementary excitations. In Fig. 1(b)we show typical excitation spectra of the SF+ states. In
the SF+ (SF− ) phase, the excitation energy of a rotonlike
minimum at a finite quasimomentum approaches zero as δ isdecreased (increased), and the softening of the rotonlike mode

induces the ST transition, similar to the standard superfluidsupersolid transition [24]. The momentum of the rotonlike
excitations largely determines the characteristic reciprocal
vector q̄1 + q̄2 of the ST state. Furthermore, the SF± -ST
transition can also be first order as indicated by the redsolid line shown in Fig. 1(a). In this case, the energy gapof rotonlike excitations jumps discontinuously to zero at theSF± /ST boundary, playing the role of an order parameter for
the first-order phase transition. The phase boundary where
the roton minimum jumps discontinuously corresponds to thelocation of the Maxwell construction, where the free energy of
the two phases have the same value. This is a generic featureof first-order (discontinuous) phase transitions. The roton
minimum has hysteretic behavior and remains finite in theregion of metastability of the hysteresis curve before jumping
to zero when sweeping from the SF± to ST phases. Upon areverse sweep from ST to SF± the roton mininum remains atzero below the phase boundary and then jumps discontinuously
back to a finite value at the end of the metastability region. Suchbehavior is similar to those encountered in continuum modelsof spin-orbit-coupled Bose-Einstein condensates [25,26].
The weak-coupling phase diagram shown in Fig. 1(a)reveals ground states which are very similar to those in thecontinuum limit [4,10–13], where the band structure due to theoptical lattice is not important. However, the phase diagramof SOC momentum kT /π versus h̄/t at ρ↓ = ρ↑ , shownin Fig. 1(c), illustrates the remarkable competition between
the reciprocal vector of the underlying optical lattice andcharacteristic vector q̄1 + q̄2 of the ST phase. In the spinsymmetric case (ρ↑ = ρ↓ ), the two wave vectors q̄1 and q̄2
are equal, that is, q̄1 = q̄2 ≡ q̄ leading to q̄1 + q̄2 = 2q̄. The
phase diagram in the range of kT = π to 2π is exactly the sameas that of Fig. 1(c) since the lattice Hamiltonian Ĥ0 + Ĥint isinvariant under the gauge tranformation b̂ks → b̂k+(π,0,0)s , aseasily verified by direct substitution.
When kT is nearly commensurate to the lattice reciprocal
wave number 2π , such as kT ≈ π/4, 2π/3, and π/2, thepitch vector q̄ of the ST state spontaneously takes anexact commensurate value over a finite range of kT . Higher
harmonic generation due to interactions and to umklapp
process q1 + q2 − q3 − q4 = 2π n with nonzero integer nfavor commensurate phases, since their energy is lower
than those of incommensurate ones. As a result, the curveof q̄/π versus kT /π exhibits multiple plateaux in the ST
phase [Fig. 1(d)]. In particular, when kT ≈ π/2, BEC occurs with only two momenta ±q̄ = ±π/2 since all the

061603-2

RAPID COMMUNICATIONSQUANTUM PHASES OF TWO-COMPONENT BOSONS WITH . . .

PHYSICAL REVIEW A 96, 061603(R) (2017)

FIG. 2. Ground-state phase diagrams in the (t/U,μ/U ) plane,
obtained by the Gutzwiller self-consistent calculations for different
values of h̄/t. We set the other parameters as U↑↓ = 0.9U , kT =0.2π , and ρ↑ = ρ↓ (δ = 0).
FIG. 1. Ground-state properties in the weak-coupling regime with
U/t = 1/ρ and U↑↓ = 0.9U . (a) Phase diagram of detuning δ versus
Rabi frequency  for kT = 0.2π . The thick red (thin black) curves
denote first- (second-) order transitions and the black dots indicate
multicritical points. In the δ > 0 (δ < 0) region to the left side of thedash-dotted line, the SF− (SF+ ) exists only as a metastable state. (b)Rotonlike softening in the elementary excitations for quasimomentum
k = (kx ,0,0) and h̄/t = 0.4. We set h̄δ/t = 0.4 (in SF+ ) for the
dotted lines and h̄δ/t = 0.06 (at the SF+ -ST boundary) for the solidlines. (c) The kT dependence of the ground state when ρ↑ = ρ↓ (δ =0). The yellow and darker green regions limited by the black-dashedand red lines are the CSF and period-locked ST phases illustrated in(e). (d) The plateaux in q̄ of the CSF and period-locked ST phases asa function of kT for h̄/t = 1.0. The dashed lines denote the width of
dominant plateaux with commensurate wave number q̄. (e) Density(the size of dots) and chiral (the direction of arrows) patterns in thecommensurate phases.

higher-harmonics momenta are reduced to ±π/2 due tothe Brillouin zone periodicity. In this special case whereq̄/π = 1/2, the interference of the two momenta does notlead to striped density pattern, but to Z2 chiral symmetry
breaking. This state is analogous to the chiral superfluid (CSF)
state, which has been discussed in Bose-Hubbard ladders[27–30]. In the present case, the 3D lattices for the two
spin components and the Rabi couplings play the role ofrails and rungs, respectively, of a synthetic “two-leg ladder”
in four (three spatial plus one extra spin) dimensions as
illustrated in Fig. 1(e). For other commensurate ST phases,
where q̄/π takes an irreducible fraction ζ /η with ζ and η
being integers, the superfluid phases break Zη symmetry, but
preserve a stripe pattern in the atom density. The stabilizationof these commensurate phases is a specific feature of spin-

orbit-coupled systems in optical lattices with interactions andare completely absent in interacting continuum systems.
In addition to the interplay between different length andmomentum scales discussed above, another particular featureof lattice systems is the existence of Mott insulator (MI) phases
induced by strong interactions and commensurate particle
fillings. To describe the Mott physics in the presence of SOC
and Zeeman fields, we employ the Gutzwiller variational
method [22]. Under the assumption that the ground state isgiven by a direct product state in real space, the Hamiltonianis mapped into inequivalent single-site problems connected
via mean fields ψis ≡ b̂is , where b̂is is the spin-s boson
operator at lattice site i. To deal with the nonuniformity of
the ST phases, we solve a set of self-consistency equations
with 2000 mean fields ψis along the x direction for each spin
and thus the momentum resolution is δkx ∼ 0.001π , while
the y and z directions are assumed to be uniform [22]. In the
ST phase, the inhomogeneous state is a result of the lengthscale introduced by the SOC, while in the absence of SOC a
new length scale leading to a supersolid state appears due tolong-range interactions [31].
Figure 2 shows phase diagrams in the μ/U -t/U plane forseveral values of  in the spin-symmetric case ρ↑ = ρ↓ (δ = 0).
In Fig. 2(a), where there is no spin hybridization ( = 0),
the phase boundaries of the MI lobes are identical to thosein the absence of SOC [32] since the gauge transformationb̂ks → b̂k+skT s eliminates the momentum transfer kT fromthe problem. The even-filling Mott transitions become firstorder in a two-component Bose-Hubbard model for large
intercomponent repulsions (for example, U↑↓  0.68U when
ρ = 2) [32–35]. In the superfluid phase outside the Mott lobes
for kT = 0, the spin-down and spin-up bosons independently

061603-3

RAPID COMMUNICATIONSYAMAMOTO, SPIELMAN, AND SÁ DE MELOPHYSICAL REVIEW A 96, 061603(R) (2017)

FIG. 3. Nonuniform superfluid-insulator transitions at ρ = 2 for(a) U↑↓ = 0.9U and (b) U↑↓ = 0.2U . We set kT = 0.2π and ρ↑ = ρ↓(δ = 0). The vertical dashed line in (a) marks h̄/t = 0.72. The
enlarged view of the region indicated by the dashed box in (b) isshown in (c). The fourth-order Ginzburg-Landau coefficients and thevalue of q̄ along the MI transition line of (c) are plotted in (d).

form the SF+ state with q̄ = kT and SF− with −q̄ = −kT ,respectively.
When the Rabi frequency  is nonvanishing, the two spincomponents mix, forming a nonuniform ST state with two
opposite momenta −q̄ and q̄ and their associated higher
harmonics. Figure 2(b) shows that the transition from the
odd-filling MI to the ST phase occurs via an intermediate SF0
state. A direct transition to the ST state occurs only for verysmall h̄/t (not shown: h̄/t  0.04 for ρ = 1). As seen in
Figs. 2(c) and 2(d), when the value of h̄/t is increased, theSF0 phase also emerges near the tip of the even-filling MI
lobes, and eventually joins other SF0 regions. The SF+ andSF− states only phase separate for small fillings ρ  1 and avery narrow region around the ρ = 2 MI lobe for large h̄/t.
To see the interplay between local correlations and spin
mixing, we plot in Figs. 3(a)–3(c) phase diagrams of U/tversus h̄/t for fixed density ρ = 2. Moreover, the natureof the superfluid-insulator transition when ρ↑ = ρ↓ can be
described by the Ginzburg-Landau energyEGL= ξ (k)
M

2
I +

2
II

+

1

2

4
I +

4
II

+

2

2
I

2
II

(4)

up to fourth order of the order parameters I = | q̄ | andII = | −q̄ | for the BEC with k = (±q̄,0,0). Note that thehigher harmonics are negligible in the vicinity of the transition.
The value of q̄ is determined so that the function ξ (k) attains
its minimum value −μ̄ at k = (±q̄,0,0). When μ̄ > 0, thebosons condense at q̄ and/or −q̄ with q̄ = 0, or simply at q̄ = 0.
For 1 < 2 , the minimization of Eq. (4) gives | q̄ | = 0 and| −q̄ | = 0 (or vice versa), and thus the Z2 symmetry related
to q̄ or −q̄ is broken. In this case, the transition from MI
to PS takes place. On the other hand, the condition 1 > 2
gives | q̄ | = | −q̄ | ≡ = 0, resulting in the transition to
the ST or CSF phase. When q̄/π is an irreducible fractionζ /η, the relative phase φ = arg( q̄ / −q̄ ) is determined by
the minimization of an additional η-particle umklapp process,FIG. 4. The crystal momentum distributions b̂†ks b̂ks  [k =(kx ,0,0)] of the four different states along the line of h̄/t = 0.72 in
Fig. 3(a) (at U/t = 0.5, 1.5, 30, and 44). The contribution from SF+(SF− ) in the PS phase is plotted by the solid (dashed) lines.

+ ( ∗−q̄ )η ( q̄ )η ) ∝ cos ηφ, which still has
η-fold degeneracy. Thus the ST transition is associated with
U (1) × Zη symmetry breaking about the global and relative
phases of ±q̄ .
The coefficients ξ (q), 1 , 2 , and η are related to themicroscopic system parameters in Ĥ0 + Ĥint by performing a
perturbative expansion based on a direct-product MI state. For
the specific relations, see the Supplemental Material [22]. Weshow in Fig. 3(d) the values of 1 and 2 along the line that
separates the MI phase from the others as seen in Fig. 3(c). Notethat if 1 < 0 for 1 < 2 or 1 + 2 < 0 for 1 > 2 , thecondensates have a negative compressibility, and the transition
becomes first order.
To assist in the experimental identification of these quantum
phases, Fig. 4 shows the crystal momentum distribution†
b̂ks b̂ks , which does not include the effects of Wannierfunctions, but can be easily extracted from standard momen†
tum distribution measurements. We evaluate b̂ks b̂ks  via the
Bogoliubov Hamiltonian ĤB for the PS and ST states atrelatively weak interactions, and via a generalized HolsteinPrimakoff approach based on the Gutzwiller variational state
for the SF0 and MI states in the strongly coupled regime [22].
Since the PS state consists of independent domains of SF+ andSF− , we plot the simple average of the two contributions.
As seen in Fig. 4, the momentum distribution of the PS state
exhibits two independent peaks around k = (q̄,0,0) and k =(−q̄,0,0), while the ST state shows additional peaks due to thehigher harmonics. The SF0 state exhibits a peak around k = 0
as in the case of a standard uniform superfluid state, although
the reflectional symmetry with respect to kx → −kx is absent
for each spin component. In the MI state, only a broad peakis observed at the momenta where the condensation occurs in
the neighboring superfluid state. The stark differences between
these crystal momentum distributions also enable the direct
imaging of the different phases present in inhomogeneous
trapped systems.
In summary, we investigated the quantum phases of twocomponent bosons in optical lattices as a function of spinorbit coupling, Rabi frequencies, and interactions. In phase
diagrams at zero detuning, we identified four different regionsη ((061603-4∗ ηq̄ ) (−q̄ )ηRAPID COMMUNICATIONSQUANTUM PHASES OF TWO-COMPONENT BOSONS WITH . . .

occupied by uniform, nonuniform, and phase-separated superfluids or Mott insulators. Finally, we characterized thesephases by calculating their crystal momentum distributions,which can be easily measured experimentally.

PHYSICAL REVIEW A 96, 061603(R) (2017)

We thank N. E. Lundblad and D. Trypogeorgos for a carefulreading of the manuscript. D.Y. acknowledges the support of
CREST, JST Grant No. JPMJCR1673, and of KAKENHI

from the Japan Society for the Promotion of Science, Grant
No. 26800200. I.S.B. acknowledges the support of AFOSRs
Quantum Matter MURI, NIST, and the NSF through the
PFC at the JQI. CARSdM acknowledges the support of
JQI and NIST via its visitors program, the Galileo Galilei
Institute for Theoretical Physics via a Simons Fellowship,
and the Aspen Center for Physics via NSF Grant No.
PHY1607611.

[1] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045
(2010).
[2] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
[3] Non-Centrosymmetric Superconductors, edited by E. Bauer and
M. Sigrist (Springer-Verlag, Berlin, 2012).
[4] Y.-J. Lin, K. Jiménez-García, and I. B. Spielman, Nature(London) 471, 83 (2011).
[5] J.-Y. Zhang, S.-C. Ji, Z. Chen, L. Zhang, Z.-D. Du, B. Yan, G.-S.
Pan, B. Zhao, Y.-J. Deng, H. Zhai, S. Chen, and J.-W. Pan, Phys.
Rev. Lett. 109, 115301 (2012).
[6] P. Wang, Z.-Q. Yu, Z. Fu, J. Miao, L. Huang, S. Chai, H. Zhai,
and J. Zhang, Phys. Rev. Lett. 109, 095301 (2012).
[7] L. W. Cheuk, A. T. Sommer, Z. Hadzibabic, T. Yefsah, W.
S. Bakr, and M. W. Zwierlein, Phys. Rev. Lett. 109, 095302(2012).
[8] H. Zhai, Int. J. Mod. Phys. B 26, 1230001 (2012); V. Galitski
and I. B. Spielman, Nature (London) 494, 49 (2013).
[9] J.-R. Li, Jeongwon Lee, W. Huang, S. Burchesky, B. Shteynas,
F. Ç. Top, A. O. Jamison, and W. Ketterle, Nature (London) 543,
91 (2017).
[10] C. J. Wu, I. M. Shem, and X. F. Zhou, Chin. Phys. Lett 28,
097102 (2011).
[11] T. L. Ho and S. Zhang, Phys. Rev. Lett. 107, 150403 (2011).
[12] Y. Li, L. P. Pitaevskii, and S. Stringari, Phys. Rev. Lett. 108,
225301 (2012).
[13] Y. Li, G. I. Martone, L. P. Pitaevskii, and S. Stringari, Phys. Rev.
Lett. 110, 235302 (2013).
[14] Z. Chen and Z. Liang, Phys. Rev. A 93, 013601 (2016).
[15] G. I. Martone, T. Ozawa, C. Qu, and S. Stringari, Phys. Rev. A
94, 043629 (2016).
[16] C. Hickey and A. Paramekanti, Phys. Rev. Lett. 113, 265302(2014).
[17] J. Higbie and D. M. Stamper-Kurn, Phys. Rev. Lett. 88, 090401(2002).

[18] N. Goldman, I. Satija, P. Nikolic, A. Bermudez, M. A. MartinDelgado, M. Lewenstein, and I. B. Spielman, Phys. Rev. Lett.
105, 255302 (2010).
[19] B. M. Anderson, I. B. Spielman, and G. Juzeliūnas, Phys. Rev.
Lett. 111, 125301 (2013).
[20] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch,
Nature (London) 415, 39 (2002).
[21] A. Widera, F. Gerbier, S. Fölling, T. Gericke, O. Mandel, and
I. Bloch, New J. Phys. 8, 152 (2006).
[22] See Supplemental Material at http://link.aps.org/supplemental/
10.1103/PhysRevA.96.061603 for technical details of calculations.
[23] J. H. P. Colpa, Physica A 93, 327 (1978).
[24] S.-C. Ji, L. Zhang, X.-T. Xu, Z. Wu, Y. Deng, S. Chen, and J.-W.
Pan, Phys. Rev. Lett. 114, 105301 (2015).
[25] G. I. Martone, Y. Li, L. P. Pitaevskii, and S. Stringari, Phys. Rev.
A 86, 063621 (2012).
[26] R. Liao, O. Fialko, J. Brand, and U. Zülicke, Phys. Rev. A 92,
043633 (2015).
[27] E. Orignac and T. Giamarchi, Phys. Rev. B 64, 144515 (2001).
[28] A. Dhar, M. Maji, T. Mishra, R. V. Pai, S. Mukerjee, and A.
Paramekanti, Phys. Rev. A 85, 041602(R) (2012).
[29] M. Atala, M. Aidelsburger, M. Lohse, J. T. Barreiro, B. Paredes,
and I. Bloch, Nat. Phys. 10, 588 (2014).
[30] S. Greschner, M. Piraud, F. Heidrich-Meisner, I. P. McCulloch,
U. Schollwöck, and T. Vekua, Phys. Rev. A 94, 063628 (2016).
[31] V. W. Scarola, E. Demler, and S. Das Sarma, Phys. Rev. A 73,
051601(R) (2006).
[32] D. Yamamoto, T. Ozaki, C. A. R. Sá de Melo, and I. Danshita,
Phys. Rev. A 88, 033624 (2013).
[33] A. Kuklov, N. Prokof’ev, and B. Svistunov, Phys. Rev. Lett. 92,
050402 (2004).
[34] P. Chen and M. F. Yang, Phys. Rev. B 82, 180510(R) (2010).
[35] T. Ozaki, I. Danshita, and T. Nikuni, arXiv:1210.1370v1.

061603-5