Nematic quantum liquid crystals of bosons in frustrated lattices
Abstract
The problem of interacting bosons in frustrated lattices is an intricate one due to the absence of a unique minimum in the singleparticle dispersion where macroscopic number of bosons can condense. Here we consider a family of tightbinding models with macroscopically degenerate lowest energy band, separated from other bands by a gap. We predict the formation of exotic states that spontaneously break rotational symmetry at relatively low filling. These states belong to three nematic phases: Wigner crystal, supersolid, and superfluid. The Wigner crystal phase is established exactly at low filling. Supersolid and superfluid phases, at larger filling, are obtained by making use of a projection onto the flat band, construction of an appropriate Wannier basis, and subsequent meanfield treatment. The nematic superfluid that we predict is uniform in real space but has an anisotropic momentum distribution, providing a novel scenario for Bose condensation with an additional nematic order. Our findings open up a promising direction of studying microscopic quantum liquid crystalline phases of bosons.
pacs:
67.85.d, 75.10.Jm, 67.80.K, 75.10.KtI Introduction
Interactions take the center stage when the kinetic energy of particles is quenched. The most prominent example is the fractional quantum Hall effect, where Coulomb interaction lifts the degeneracy of a partially filled Landau level and leads to formation of topological fractionalized states. There, kinetic energy of particles is quenched due to the magnetic field that localizes them within the magnetic length, leading to formation of flat bands – the Landau levels. Flat bands also appear in some tightbinding lattice models, as is the case for the kagome lattice. Hence, a question naturally arises – what kind of exotic states can emerge in such systems?
The problem of interacting bosons in frustrated lattices with a flat band has been a subject of significant theoretical interest Huber and Altman (2010); Wu et al. (2007); Bergman et al. (2008); Wang et al. (2011); You et al. (2012); Wang et al. (2012); Sedrakyan et al. (2014). For exactly degenerate single particle spectrum, there is no preferred momentum state for the Bose condensation to occur. Depending on the specific model and parameter regime, bosonic ground states in such frustrated lattices may include chiral compositefermion states of hardcore bosons Sedrakyan et al. (2014); Yang et al. (1993); Kumar et al. (2014) and chiral superfluid/Mott insulator states Zaletel et al. (2014) which spontaneously break timereversal symmetry, fractional Chern insulators Wang et al. (2011, 2012), and other exotic broken symmetry states Huber and Altman (2010); Bergman et al. (2008); Wu et al. (2007); Huse and Rutenberg (1992); You et al. (2012); Möller and Cooper (2010).
From the experimental standpoint, rapid advance of artificial condensed matter systems such as cold atoms and interacting photons in circuit QED system have enabled not only the realization of geometrically frustrated lattices but also lattices subject to synthetic gauge fields Jaksch and Zoller (2003); Dalibard et al. (2011); Goldman et al. (2014); Gerbier and Dalibard (2010); Lin et al. (2012); Aidelsburger et al. (2013); Ji et al. (2014); Jotzu et al. (2014); Koch et al. (2010); Nunnenkamp et al. (2011); Hafezi et al. (2013); Hafezi and Rabl (2012). These recent developments make the search for accurate theoretical approaches and concrete proposals timely.
In this paper, we focus on a modified kagome lattice model, constructed in a way to allow controlled treatment thanks to a spectral gap. This gap can be generated by inserting an additional gauge flux into each hexagon of the kagome lattice Katsura et al. (2010). We show that the singleparticle eigenstates comprising the gapped flat band can be chosen as localized loop states, which typically break the lattice point group symmetry. Previous works in the context of the BoseHubbard model have primarily considered the simple kagome lattice Sachdev (1992); Yang et al. (1993); Schulenburg et al. (2002); Huber and Altman (2010); Huse and Rutenberg (1992); Balents et al. (2002); You et al. (2012); Capponi et al. (2013); Punk et al. (2014); Kumar et al. (2014). There, the lowest band is flat but gapless since it is in contact with another band at the point Bergman et al. (2008); Huber and Altman (2010); You et al. (2012). This makes the analysis of the interacting problem quite subtle due to the ability of particles to leak easily into the higher band.
The introduced gap enables a wellcontrolled projection onto the flatband subspace Huber and Altman (2010); Zhitomirsky and Tsunetsugu (2004, 2005); Tovmasyan et al. (2013); Takayoshi et al. (2013) in the weakinteraction regime, , and yields an effective lowenergy Hamiltonian applicable to a wide filling range. This is analogous to the lowestLandau level projection employed in the fractional quantum Hall effect. Depending on the filling fraction of the lattice, we find three types of exotic nematic phases. At close packing of maximally compact loop states, a nematic Wigner crystal is the exact ground state of the system. In the specific case of flux and higher filling fraction, our meanfield treatment predicts transitions to a nonuniform nematic supersolid followed by a uniform nematic superfluid phase.
The nematic superfluid phase is quite unusual since it is not featureless but contains internal structure in its microscopic manybody wavefunction. The lattice rotational symmetry is spontaneously broken due to the anisotropic correlations among the loop orbitals. Such anisotropic internal structure is encoded in the momentum distribution, i.e. the Fourier transform of the realspace correlation function. In addition to the standard deltafunction peak, there is an anisotropic and squeezed continuous background in the momentum distribution. It clearly reveals a novel nematic Bose condensation and can be detected through timeofflight imaging in the context of ultracold atoms. In addition, we show that the nematicity can also manifest itself in macroscopic quantities, namely the anisotropic superfluid stiffness tensor and superflow, which can be probed with phase imprinting techniques in ultracold atom setups Isoshima et al. (2000).
From a broader perspective, the possibility of such microscopic liquid crystalline phases has been pointed out previously in the context of strongly correlated electronic materials Kivelson et al. (1998); Fernandes et al. (2014). While in this article, we focus on translationally invariant/periodic states, additional rich physics is associated with topological defects de Gennes and Porst (August 10, 1995) and warrants future studies.
The paper is organized as follows. Section II shows the interacting boson model we study throughout the paper. Section III discusses the singleparticle eigenstates of the model, including the gapped flat band structure and the presence of localized loop eigenstates in real space and their specific properties such as “flux quantization”, with the detailed derivation shown in Appendix A. In Section IV, we work out the exact Wignercrystal ground states below the closepacking filling of loop states. In Section V, we study the quantum phases beyond the closepacking filling. We apply a flatband projection based on the construction of mutually orthogonal and spatially compact Wannier states, with the details of construction shown in Appendix C. We then perform a subsequent meanfield analysis in the Wannier basis, which predicts the existence of nematic superfluid and supersolid phases. In Section VI, we show the novel signatures of the nematic superfluidity, namely the anisotropic momentum distribution and the anisotropic superflow. The detailed calculations of the momentum distribution can be found in Appendix F. We conclude our work and provide a brief outlook in Section VII.
Ii Model
The BoseHubbard model on the kagome lattice subject to gauge flux [Fig. 1] is described by the Hamiltonian
(1) 
where creates a single boson on the site labeled . The first term is the tightbinding Hamiltonian determining the band structure of the noninteracting bosons. We denote the hopping amplitude by to stress that it is positive (frustrated) ^{1}^{1}1In ultracoldatom experiments, the natural negative hopping can be turned positive by threading triangles of the kagome lattice with a flux.. The gauge potential is defined on each nearestneighbor bond determines the flux threading each plaquette in the lattice. The second term in captures the repulsive Hubbard interaction on each site with strength .
Iii Singleparticle eigenstates
iii.1 Gapped flat band in the presence of flux
We first discuss the tightbinding band structure for . In the absence of external flux, the singleparticle spectrum has a lowest flat band as shown in Fig. 2(a). However, the flat band touches the higher dispersive band at the point []. Once flux is inserted into each hexagon of the kagome lattice, the singleparticle spectrum takes the typical Hofstadter butterfly form [Fig. 1]. For noninteger , the lowest band ^{2}^{2}2The term “band” is applied loosely here. For rational flux, the unit cell is enlarged and the flat band decomposes into multiple ones. No such simple picture applies to irrational flux. remains flat but acquires a gap that reaches its maximum size of at flux [Fig. 2(a)]. At this point, timereversal (TR) symmetry is intact and the model can be realized with realvalued hopping of positive and negative sign. The energy and degeneracy of the flat band are independent of the flux. The latter is given by the number of hexagons in the lattice, ( denoting the number of sites).
iii.2 Localized loop eigenstates
The presence of the flat band is directly linked to the existence of degenerate eigenstates that form localized loops. The localization mechanism, called caging Vidal et al. (1998); Häusler (2015), is due to the destructive interference of the wavefunction amplitude anywhere outside the loop Bergman et al. (2008). In the kagome lattice with positive hopping, the flat band and localization persist as long as there is no flux through triangles.
To be systematic, we list three important properties of the single particle eigenstates in the flat band of the tightbinding kagome Hamiltonian:

The energy of the singleparticle eigenstate is exactly 2.

The elementary flatband eigenstates are singleparticle loop states that have equal probability on each involved site. The amplitudes on the adjacent sites outside the loop (on the outward and inward triangles) cancel due to destructive interference (caging). Any other (nonelementary) flatband state can be composed as a linear superposition of loop states.

“Flux quantization”: any loop eigenstate encloses an integer number of flux quanta,
(2) Here, the direction of the gauge potentials is chosen to be counterclockwise () around the loop.
Previous studies Bergman et al. (2008); Zhitomirsky and Tsunetsugu (2004, 2005); Huber and Altman (2010) have focused on the zeroflux case and identified the state with amplitudes of equal magnitude but alternating signs on the hexagon loop as maximally compact eigenstates [see Fig. 3(a)]. While the existence of flat band and localized eigenstates remains unharmed by the flux through hexagons, we find that the shapes of loop states must change. Specifically, in the spirit of the fluxquantization condition (property 3), maximally compact loop eigenstates encircling two and four hexagons in Figs. 3(b) and 3(c), respectively. Note that orientation and the shape of the maximally compact loop states are not generally unique for .
The main focus of this paper is the case of flux where the TR symmetry is intact and maximally compact loop states are dimers encircling two hexagons. By a convenient gauge choice on decorated bonds [Fig. 1(b)], all hopping elements are real and given by on regular and on decorated bonds. In this gauge, the amplitudes of loop eigenstates simply alternate in sign across positivehopping bonds and are identical across decorated negativehopping bonds. In the following, we consider the occupation of these states by multiple bosons and refer to the maximally compact loop states as loop orbitals (LOs).
Iv Exact nematic Wigner crystal ground state
We next turn to the interacting case, accounting for onsite boson repulsion due to the Hubbard term
(3) 
Since the interaction is local, we note that any manybody state of the form
(4) 
with singleparticle occupation of a set of nonoverlapping LOs is an exact ground state of the interacting system for filling ^{3}^{3}3Note: for filling below close packing, there will generally be a large number of degenerate ground states.. Here, the operator creates a single particle occupying the LO labeled by . Indeed, the above product state is an eigenstate with eigenenergy per particle and interaction does not contribute since double occupancy of sites is avoided.
Once the filling reaches close packing, the ground state becomes an incompressible Wigner crystal Wu et al. (2007). No additional particle can be placed on the lattice without incurring an interactioninduced energy increase due to unavoidable overlap. At the critical filling of close packing, bosons occupy maximally compact LOs while avoiding double occupation. As discussed above, maximally compact LOs may break the lattice point group symmetry (here, ), which directly leads to ground states with spontaneously broken lattice symmetry.
In general, the filling fraction for close packing depends on flux. In the flux case, maximally compact LOs are dimers and close packing occurs at . Due to the three possible orientations of a dimer [Fig. 1(b)] and the freedom to use one Wigner crystal representative [Fig. 1(d)] and produce four other inequivalent ones by translations to four neighboring hexagons, we predict that the ground state is overall 15fold degenerate. These ground states are nematic Wigner crystals. Here, nematicity refers to the emergence of dimers that break the lattice rotation symmetry. In this aspect, the flux case is dramatically different from the flux case studied before in the context of both the boson model and the antiferromagnetic Heisenberg spin model ^{4}^{4}4Note that the Wigner crystal state, is first discovered in the spin context, and termed as valencebond crystal. It is the exact ground state corresponding to the magnetization plateau Schulenburg et al. (2002) of the antiferromagnetic kagome Heisenberg model in an external magnetic field. See Appendix B for details., where the Wigner crystal Huber and Altman (2010) and the valencebond crystal Schulenburg et al. (2002) ground states do not exhibit any nematicity. We note that the nematic Wigner crystal state can also be found in the magnetization plateau of the corresponding spin model in the presence of both positive and negative XY interaction (see Appendix B).
We note that for and filling below , bosons form an infinitely compressible hardcore loop gas with macroscopic degeneracy determined by all possible configurations of nonoverlapping loops Zhitomirsky and Tsunetsugu (2004, 2005). One such configuration is depicted in Fig. 4(b), where an extra particle can be added to an unoccupied loop orbital without costing any interaction energy. Hence the chemical potential, i.e. energy cost per extra particle , is fixed to be the flatband energy and hence does not change with the filling, i.e. . Equivalently, we get , which means infinite compressibility. As shown in Fig. 4(b), by adding a local perturbation, certain loops can move freely to a nearby vacancy (dashed loops) and hence make a transition to another state with the same groundstate energy. The shown state (and other states which connect to this state by a local perturbation) breaks the lattice rotational symmetry and hence is also nematic. Finally, we mention that there are infinitely many incompressible glassy states below , which cannot be connected to other ground states by a local perturbation. We will leave the discussion of these glassy states to future works.
V Nematic superfluid and supersolid
v.1 Flatband projection and construction of Wannier orbitals.
In the following, we exclusively focus on the flux case. For filling above close packing, interaction cannot be avoided anymore and hence no exact solution in the above manner is possible. To make approximations, we derive a lowenergy effective Hamiltonian by adapting the approach by Huber and Altman Huber and Altman (2010), consisting of a projection onto the subspace spanned by flatband eigenstates. In our case of nonzero flux, however, we forego the more subtle situation of an ungapped band encountered in Ref. Huber and Altman (2010). In the presence of a gap and in the weakinteraction limit, boson occupation is to a good approximation limited to the flat band and the projection is appropriate unless the filling fraction becomes too large (details depend on the ratio ).
To facilitate the projection, we construct an orthonormal basis of the flat band. For flux, the unit cell is doubled and contains a left and right hexagon, and , which differ by the relative positions of negativehopping bonds [Fig. 5(a)]. Due to the unitcell doubling there are, strictly speaking, two degenerate flat bands. Accordingly, we choose two sets of maximally compact dimer LOs aligned in the direction [Fig. 5(a)] as our basis for the two degenerate flat bands. We distinguish leftdimer states only containing hexagons from rightdimer states only containing hexagons. Although these sets of LOs together form a basis of the two degenerate flat bands, not all basis states are mutually orthogonal. We thus need to determine appropriate superpositions of the dimer LOs to form a set of mutuallyorthogonal Wannier orbitals (WOs). As usual, there is not a unique set of WOs and different choices can vary significantly in their realspace localization. Since we will ultimately employ localdecoupling meanfield theory, it is particularly important to obtain welllocalized WOs ^{5}^{5}5 Note that it is possible to construct symmetric but less compact localized orbitals as our flatband basis. However, the meanfield ansatz with such orbitals has a much larger interaction energy cost due to the larger overlap and hence is not energetically favorable (see Appendix G for details)..
Our construction scheme for suitable WOs involves an important step of orthogonalizing the sets of left and right LOs by means of a symmetrized version of the GramSchmidt procedure (see Appendix C for details). The results for two adjacent WOs are depicted in Fig. 5(b). The major part of the realvalued WO amplitude is essentially concentrated on each original dimer [Fig. 5(a)]. From there, the amplitudes decrease rapidly (asymptotically in an exponential fashion). This is in contrast to the slower powerlaw decay of WO amplitudes in the 0flux case, which is caused by the touching of bands Huber and Altman (2010). The WOs we obtain respect translational symmetry (in terms of probability), TR symmetry, and preserve the mirror symmetry along their major axes, just as the original dimer LOs. They weakly break mirror symmetry along their minor axes.
We define the creation operator for occupation of these Wannier orbitals by
(5) 
where the Wannier function gives the amplitudes of the dimertype WO centered at position of the effective triangular lattice [Fig. 5(a)] on each site of the underlying kagome lattice. The flatband projection corresponds to the inverse transformation
(6) 
where the Wannier states of the dispersive bands have been dropped as an approximation. Upon projection and switching to the grandcanonical ensemble, the effective Hamiltonian takes the form
(7) 
where is the chemical potential. For convenience, we may define the shifted chemical potential which absorbs the energy constant of the flat band. The coefficients determine the strength of the effective interaction terms and involve overlaps of four Wannier functions centered on specific sites , , , and of the triangular lattice. Due to the localization of WOs, the interaction is short range and falls off rapidly with growing spatial distance between the four sites. We note that is translationally invariant and realvalued (since the constructed Wannier functions are realvalued themselves).
The distinct spatial configurations of the four dimer WOs labeled through give rise to different types of interaction terms. Whenever all four indices coincide, the contribution corresponds to an effective onsite repulsion with strength . Among the set of all effective interaction terms, this onsite repulsion term has the largest strength. The next subleading terms come from two other types of effective interaction, namely, densitydensity repulsion [Fig. 5(d)], and assisted hopping [Fig. 5(e)]. Here, denotes the Wannier number operator, and the primes on sums signal that those terms with coinciding summation indices are to be omitted. The strengths of densitydensity interaction and assisted hopping depend on the specific arrangement of the involved WO dimers. The largest contributing terms are and , and thus significantly smaller than the onsite repulsion strength .
Within the lowdensity regime (i.e., in the effective triangular lattice), we therefore employ a hardcore approximation which forbids double occupation of WOs Huber and Altman (2010). Within this approximation, interaction terms with repeated Wannier operators on the same site, or , drop out. This includes effective onsite interaction as well as pair hopping . Besides densitydensity repulsion and assisted hopping, the only remaining interaction type is ringexchange [Fig. 5(f)], in which the Wannier functions are centered on four different sites on the triangular lattice. We find that the maximum strength of ring exchange is , which is significantly weaker than both densitydensity repulsion and assisted hopping.
v.2 Meanfield theory.
While densitydensity repulsion favors densitywave order and formation of a Wigner crystal, assisted hopping may lead to melting and formation of a superfluid. In addition, this competition also allows for an intermediate supersolid phase in which both types of order are present. Here, we study the competition between different types of orders within meanfield theory (MFT). We adopt the Gutzwiller approach Rokhsar and Kotliar (1991) and employ a product ansatz consistent with the hardcore constraint
(8) 
which decouples sites on the effective triangular lattice of WOs. The meanfield ansatz naturally captures the nematic Wigner crystal phase since it is a product of singleparticle states with occupation of nonoverlapping LOs (in this case approximated by WOs). Above close packing, meanfield solutions continue to break the symmetry due to the anisotropic nature of the Wannier orbitals.
To describe states with densitywave order such as the nematic Wigner crystal, we must allow for the dependence of the meanfield amplitudes on the spatial index . To obtain meanfield solutions, we decouple the effective Hamiltonian, replacing densitydensity interaction and assistedhopping terms by
(9) 
(We have verified that inclusion of ringexchange does not lead to significant changes.) With this, we obtain a meanfield Hamiltonian , where depends on the meanfield order parameters and on each site of the triangular lattice. Starting from a random initial set of order parameters on a lattice of 200 sites with periodic boundary conditions, we repeatedly solve for the eigenstates and recalculate order parameters until reaching selfconsistency (see Appendix E for details).
For a range of chemical potentials, we calculate results for the mean filling , densitywave order parameter defined as the difference between maximum and average density taking into account the six surrounding sites, and the mean superfluid order parameter . The key results from this calculation are presented in Fig. 6. MFT reproduces the exact nematic Wigner crystal [Fig. 6(b)] for at close packing , showing maximum densitywave order and vanishing superfluidity . Below close packing, MFT produces a gradual change of average filling and superfluid order, which differs from the exact solution discussed above based on LOs. The exact solution exhibits a density plateau at containing the entire nematic Wigner crystal phase, and a vertical jump corresponding to the hardcore loop gas phase which is more appropriately represented in the canonical ensemble.
Above , superfluid order sets in and grows gradually while, at the same time, the densitywave order parameter remains nonzero and decays slowly, overall suggesting a secondorder transition to a nematic supersolid [Fig. 6(c)] in which a fraction of the bosons condense on interstitial sites between the Wignercrystal structure. Further on at , the densitywave order abruptly drops to zero, accompanied by a sudden increase in the superfluid order . This indicates a sudden melting of the Wignercrystal structure and a firstorder transition into a superfluid phase [Fig. 6(d)].
Based on our MFT, we predict that the superfluid phase is nematic since condensation of bosons is based on hopping among anisotropic dimer WOs. Within the superfluid phase, phase angles form stripes in which neighboring stripes differ by a phase difference. The nematic supersolid has similar phase stripes, the only difference being that sites with maximum density have an additional phase flip. Finally, we find a narrow region in which the nonmonotonic dependence of the density on the chemical potential suggests phase coexistence between the superfluid and supersolid.
Vi The signatures of nematic superfluidity and detection methods
The interesting aspect of the uniform superfluid phase is that it is nematic. Its internal structure, i.e., the correlation in the loop/Wannier orbitals, is encoded in the momentum distribution, which can be probed through timeofflight (TOF) imaging in the context of coldatom experiments. The nematicity can also be identified through macroscopic quantities, such as the superfluid stiffness tensor and the anisotropic superflow, which can be probed with phase imprinting technique Isoshima et al. (2000). In the following two sections, we discuss microscopic and macroscopic signatures, along with methods to detect them.
vi.1 Momentum distribution and timeofflight experiments
The microscopic signature, i.e., the groundstate momentum distribution , as mentioned above, can be directly measured experimentally through the time of flight images. It serves as a useful probe of the correlation properties of the ground state. The momentum distribution is equal to the Fourier transform of the singleparticle density matrix:
(10) 
Recall that for the nematic Wigner crystal (NWC) ground state, the wavefunction can be expressed as . Here, is the label of the loops and is the set of nonoverlapping loop states forming the crystal. Since the loops do not overlap, we have the correlation if and only if and are sites on the same loop. Thus, we get
(11) 
where represents the boson operator on the loop (). There are two types of loops: one encircles the left hexagons and the other encircles the right hexagons, as shown in Fig. 5(a). As we can see, the nonzero contributions for the singleparticle density matrix in Eq. (11) are , where “” is determined by the sign of the overlap of wavefunction amplitudes.
Now we discuss the momentum distribution of the nematic superfluid phase (see Appendix F for derivation). It is shown in Fig. 7. The continuous background originates from the correlation within the loop orbitals; apart from a prefactor it is identical to NWC. It is squeezed in the direction of the major axis (along which the loop is elongated). This continuous background encodes the internal correlation of the loop/Wannier orbitals. The deltafunction peaks (represented by white circles) originate from the Bose condensation. Only the four equivalent peaks on the boundary of the first Brillouin zone are shown in the plot. Such a momentum distribution of a uniform superfluid implies a novel scenario of Bose condensation, where the ground state is unstable against developing an additional nematic order which, spontaneously breaks the lattice rotational symmetry. Therefore in the corresponding TOF experiment, one expects that the sample prepared under similar conditions repeatedly would spontaneously pick up one of the three directions in which the image pattern is squeezed.
vi.2 Superfluid stiffness tensor and anisotropic superflow
Now we discuss the macroscopic signature of the nematic superfluid phase. We consider the superfluid stiffness tensor , where refer to the directions along and . Note that direction is special because it aligns with the major axis of the dimers. To study this quantity, we apply phase differences and across the boundaries of the finite sample (16 16), as shown in Fig. 8(a). Experimentally such phase difference can be achieved with the phase imprinting technique developed in the cold atom setup Isoshima et al. (2000).
The phase differences across the boundaries induce superflow in the corresponding directions and hence increase the kinetic energy. The contour plot in panel (b) shows the meanfield energy as a function of the phase differences, namely . The superfluid stiffness tensor corresponds to the curvature of the energy profile in the vicinity of the origin. The anisotropy of the energy contours suggests that the superfluid stiffness is also anisotropic. To see this more clearly, we make cuts along the  and axis (blue solid and red dashed line). We then show the energy profile along the two cuts, namely and , in panel (c). It is obvious that, in the vicinity of the origin, the curvature of the blue solid line is larger than that of the red dashed line, which means that is larger than . This suggests that the superfluid stiffness along the two directions is different. Now we consider the first derivative, , which is the current generated when applying a phase difference at direction . We can see that, not far away from the origin, is always larger than . This is consistent with the results in Appendix D, where we show that the effective nearestneighbor hopping along the major axis ( direction) is zero in the nematic superfluid phase (within meanfield approximation). Only successive hopping along other directions will contribute to superflow in the direction. On the other hand, the large effective nearestneighbor hopping in the other two directions ( and ) leads to larger superflow in those directions.
In sum, the anisotropy of the two macroscopic quantities, the superfluid stiffness and superflow, reveals breaking of discrete rotational symmetry and hence microscopic nematicity of the superfluid phase.
Vii Conclusion
We studied the emergence of nematic phases in a kagome lattice with a gapped flat band, obtained when a flux is threaded through each hexagon of the lattice. Singleparticle localized loop states can be combined to construct manybody eigenstates below a critical filling. This critical filling corresponds to close packing of nonoverlapping loop states and marks the formation of a nematic Wigner crystal ground state. For larger filling, the effective Hamiltonian based on flatband projection using dimershaped Wannier orbitals and subsequent meanfield treatment predict nematic supersolid and superfluid phases. The latter is a uniform quantum liquid with anisotropic internal structure, which is encoded in the momentum distribution and can be probed by timeofflight experiment. Interesting future directions include the study of phases at higher density, especially the possibility of a featureless Mott insulator Parameswaran et al. (2013) at 1/3 filling, resonating valence bond states and fermionized ground states in the stronginteraction regime.
Acknowledgements.
We are indebted to Steven Girvin, Ashvin Vishwanath, Tigran Sedrakayan, Eliot Kapit, Hakan Türeci, TzuChieh Wei, Anupam Garg, James Sauls, Murad Tovmasyan and Andy C. Y. Li for insightful discussions. Work performed at Northwestern University (G.Z. and J.K.) was supported by the NSF under Grant PHY1055993. Work performed at Argonne National Laboratory (I.M.) is supported by the U. S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DEAC0206CH11357.Appendix A Derivation of the properties of loop eigenstates
In this appendix, we show the derivation of some properties of loop eigenstates stated in Sec. III.2.
The flux quantization condition (property 3) can be derived from properties 1 and 2 as follows. The positivehopping tightbinding lattice Hamiltonian in the presence of an additional gauge field has the form:
(12) 
where is the hopping operator on the nearestneighbor bond . Assuming that there exist singleparticle wavefunctions of the loop eigenstate type, they can be expressed as
(13) 
where the summation index runs over all sites on the loop . When acting with the Hamiltonian on the loop eigenstates, we can split the expressions into two parts, namely
(14) 
The first sum includes hopping along the bonds on the loop , while the second sum corresponds to hopping from the bonds on the loop to the adjacent sites on the outward/inward triangles, as illustrated in Fig. 9 by red dashed lines. The cancellation of the probability amplitude outside the loop (caging, property 2) requires the second sum to be zero (see Appendix A), while the requirement of eigenenergy being (property 1) implies the first sum being equal to:
(15)  
The above equation leads to the following relation between the amplitudes of neighboring sites,
(16) 
That is, the wavefunction has equal probability on every site along the loop, and adjacent sites differ by a minus sign and an additional phase shift due to the gauge potential on the bond . Note that in the 0flux case (), Eq. (16) simplifies to alternating signs on the loop, including the hexagon loop state shown in Fig. 3(a). By applying Eq. (16) around the loop and requiring the probability amplitude to be singlevalued, we derive the flux quantization condition for a loop eigenstate (property 3), namely
(17) 
Now we consider the cancellation of the second sum in Eq. (14) which leads to caging, and get
(18) 
This in turn leads to . Combined with Eq. (16), we get
(19) 
Thus, we have just shown the necessary condition for the persistence of lowest flat band with energy , which we stated in the beginning of Sec. III.2, is that the flux threading all the outward/inward triangles must be zero (gaugeequivalent to flux in the negative hopping model).
To understand one additional feature of the loop state, now we consider the gaugeinvariant current operator on the bond (from site to site ), namely,
(20) 
Its expectation value is given by
(21) 
which equals zero after applying Eq. (16). Therefore, for any flux value , the loop states in the flat band carry no current, and thus do not break timereversal (TR) symmetry, even though the Hamiltonian itself breaks TR except at . From the butterfly spectrum in Fig. 1(b), we observe that the lowest flat band does not change as a function of , which leads to zero current due to the linear response formula of the current . The fact that there is no current can also be understood in another way, i.e. the Chern number for the lowest flat band is always zero. This is in contrast to the higher bands near [Fig. 1(b)], which are essentially Landau levels with nonzero Chern numbers etc..
Appendix B Connection to frustrated spin models
Although the current paper focuses on the interacting boson models, some of the results, such as the Wigner crystal phases, also apply for the frustrated spin models corresponding to the same type of frustrated lattice. Here, we consider the frustrated anisotropic Heisenberg model with spin in the presence of external magnetic field:
(22) 
Here, determines the interaction and determines the (Ising) interaction. The isotropic situation () corresponds to the usual Heisenberg model. The external magnetic field acts as the chemical potential of magnons and has nothing to do with the gauge flux which we discuss throughout the whole paper. To make the connection with the boson model more explicit, we rewrite it as:
(23) 
From above, one can see that the transverse interaction can be written as a flipflop term, which in the case actually corresponds to hopping of hardcore bosons. The term induces nearest neighbor interactions between hardcore bosons. If the magnetic field is sufficiently large, the ground state is the magnon vacuum . For the flatband hopping models, a single magnon on the loop is created by the operator out of the magnon vacuum , where represents the wavefunction amplitude on each site of the loop. Just like in the case of onsite interacting bosons, the magnon loop gas and the loop crystal are the eigenstates of the Hamiltonian (the term does not change that since the adjacent loops do not occupy neighboring sites.
Indeed, for the antiferromagnetic Heisenberg model (), a valencebond crystal phase (equivalent to the Wigner crystal) formed by nonoverlapping hexagon loop magons has been found in Ref. Schulenburg et al. (2002) (earlier than its boson analog). This valancebond crystal phase, in the spin1/2 case, corresponds to the magnetization plateau (m=0 corresponds to no polarization; corresponds to full polarization in the up/down direction), and is equivalent to the 1/9 state of the interacting boson model. We note that the valencebond crystal phase is not limited to spin1/2 case, but actually exists for arbitrary spin Schulenburg et al. (2002).
To implement the flux model described in the main text, one can choose different signs of the transverse coupling according to the hopping signs of the corresponding boson model. Thus, there will be a nematic valencebond crystal (at 13/15 magnetization plateau) in the spin model which corresponds to 1/15 nematic Wigner crystal in the boson model. To make sure the densitydensity interaction between magnons is repulsive and hence stabilizes the valencebond crystal as the ground state, it is preferable that the interaction is positive, namely . However, we caution that, if the anisotropy ratio is sufficiently small, the sign of the term does not matter since the term can be treated as a small perturbation. To experimentally realize such a signtunable spin model, a promising candidate is the nitrogenvacancy center array Cai (2013), although in that case the spinspin interaction is not restricted to nearestneighbors but has a powerlaw decay due to its dipoledipole nature.
Besides the nematic Wigner crystal (nematic valencebond crystal) phase, the nematic supersolid or superfluid phases may also occur in the spin systems. However, due to the hardcore nature of the spin system (especially for spin 1/2), the projection approach may break down at certain critical filling. Only below that critical filling, the predicted phase for the weaklyinteracting boson model is expected to apply.
Appendix C Construction of orthogonal Wannier Orbitals
For the flux case, the unit cell is doubled with respect to the 0flux case and hence includes 6 sites (labeled as A, B, C, D, E, and F), as shown in Fig. 10(a). Hence, the lowest flat bands are doubly degenerate, which means that there is an arbitrary choice to decompose the two flat bands, since for each one can arbitrarily choose two orthogonal Bloch vectors in the 2dimensional degenerate subspace. A sensible choice of basis has to be physically motivated and has to respect certain symmetries.
For the sake of convenience, we choose loops along the direction as our preferred flat band basis (total number equals the number of flatband degeneracy). The operators , which create the two types of loops within each unit cell, can be represented by the original lattice boson operators as:
(24) 
Here, the index labels the loops encircling left/right hexagons. The wavefunctions of the two types of loops are already shown in Fig. 3(a). Here, the lattice vectors labels the enlarged 6site unit cell and the lattice vectors and translate the cells in the two oblique directions [shown in Fig. 10(a)].
One can construct two classes of Bloch states by a translationallyinvariant superposition (with a particular wavevector ) of the two types of loop states respectively, i.e.
(25) 
Here, left/right label can also be thought as the band index and in this particular case labels the two degenerate flat bands. The generated state can be represented by a 6component Bloch vector , where we have and The two Bloch vectors got from the two chosen loop states are represented as
(26)  
However, the Bloch vector is not yet normalized. We call the normalized Bloch vectors , and define the normalized Bloch state as , where the redefined operator now becomes canonical bosonic operators satisfying the commutation relation . Thus, we get a set of orthonormal Bloch states for each of the two flat bands, and can be transformed into two sets of Wannier states as:
(27) 
where the Wannier wavefunction is given by
(28) 
Here, the Wannier wavefunction sits on the coordinate . The coordinate labels where the center of the wavefunction locates. We note that the more detailed notations, and , which we use here, are equivalent to the more compact notations we have used in the main text, namely, and . The direct correspondence is and . The two sets of wavefunctions are illustrated in Fig. 10(b), where the one encircles only the left hexagons in every unit cell and the one encircles only the right hexagons. The major part of the realvalued Wannier functions are essentially the two dimer loop states, which we start with. The amplitude tail spreads out and decays exponentially along the major axis of the loop which ensures orthogonalization. However, the two sets of Wannier functions are not mutually orthogonal to each other (for example, those neighboring ones will still have finite overlap) since the two sets of Bloch vectors are not mutually orthogonalized yet.
For each , one can orthogonalize the two Bloch vectors through the GramSchmidt process, i.e.,
(29) 
which generates a normalized orthogonal to . Similarly, one can get a normalized , which is orthogonal to . Therefore, one can choose either orthogonal pair as the flatband basis. However, from either choice, the acquired Wannier wavefunctions belonging to the two bands have completely different shapes and hence lose translational symmetry. To preserve translational symmetry and being closer to a dimer shape, we can make a symmetric superposition as:
(30) 
where and are the mutuallyorthogonal sets of Bloch states we choose. The additional free choice of phase factor will give rise to different Wannier states. The sensible choice of the phase factor makes sure that the Bloch vectors are analytically continuous in kspace, which ensures that the generated Wannier function is exponentially localized Kohn (1973) and hence is more compact.
The Bloch states generated from the GramSchmidt process, , are unfortunately not analytically continuous. For example, as shown in Fig. 10(c), the component of one of the Bloch vectors, , has a diagonal discontinuity cut in its real part. Same cut happens for its imaginary part and most of the other components of , and .
Thus, one has to take advantage of the additional phase factor to remove the discontinuity. In this particular case, a simple sign flip of every other strip in the space, which can be expressed as a squarewave function: , is able to remove the discontinuity [see Fig. 10(d)]. In addition, we employ an extra phase factor to ensure no breaking of TR symmetry and closeness in shape to the dimer loop state. Thus, our choice of phase factor for Eq. (C) is (a independent relative sign or phase factor does not affect the probability distribution of the Wannier functions). This particular choice yields the complete Wannier basis illustrated in Fig. 3(b).
Our Wannier orbitals (WOs) preserve the mirror symmetry (in terms of probability) with respect to its major axis, similar to the original loop orbitals (LOs) which they are based on. However, due to the additional phase factor we choose to keep the analytical continuity, the mirror symmetry along the minor axis is slightly broken. We can see that the lower part of the WO has slightly larger probability than the upper part. If we replace part of the phase factor with , the shape of the WO will be flipped with respect to the minor axis, namely the higher part will have larger probability. We also note that we do not claim that we have found the maximally compact WOs, even though the construction is based on the maximally compact LOs. In general, it should be possible to numerically/analytically determine such maximally compact WO, which also preserve both types of mirror symmetries. Thus, our current approach is just a simple mathematical construction, which aims to approximate the maximally compact WOs, since the shape we have acquired is not too far from the original LOs, which they are based on.
Finally, we note that, since we have successfully found a complete orthogonal Bloch or Wannier basis from superposition of the dimer loop states, we have explicitly shown the completeness of the loop states which is mentioned in property 2 of Sec. III.2.
Appendix D Summary of terms in the effective Hamiltonian
Here we classify all types of effective interaction , not limited by the hardcore constraint. The types of terms are following:
(1) Onsite repulsion:
(31) 
where is the largest energy scale in the effective Hamiltonian.
(2) Densitydensity repulsion:
(32) 
where means sum over pairs of sites (). Thus, and correspond to the same term and should not be double counted. Now we determine the coefficients of the effective interaction. Four terms in the effective interaction correspond to the pair , namely . Thus, we get .
(3a) Onsite pairhopping (involving two different sites):
(33) 
where .
(3b) Offsite pairhopping (involving three different sites):
(34) 
where .
(4a) Assistedhopping (involving three different sites):
(35)  
(36) 
Four terms (and their H.c.) in the effective interaction correspond to the pair , namely . Thus, we get ().
(4b) Assistedhopping (involving only two different sites):
(37) 
Here, we have , due to the presence of two terms of each type, e.g. .
(5) Ringexchange interaction: