Nematic quantum liquid crystals of bosons in frustrated lattices

# Nematic quantum liquid crystals of bosons in frustrated lattices

Guanyu Zhu Departments of Physics and Astronomy, Northwestern University, Evanston, Illinois 60208, USA    Jens Koch Departments of Physics and Astronomy, Northwestern University, Evanston, Illinois 60208, USA    Ivar Martin Materials Science Division, Argonne National Laboratory, Argonne, Illinois 60439, USA
July 2, 2019
###### Abstract

The problem of interacting bosons in frustrated lattices is an intricate one due to the absence of a unique minimum in the single-particle dispersion where macroscopic number of bosons can condense. Here we consider a family of tight-binding 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 mean-field 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.Kt

## I 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 tight-binding 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 composite-fermion states of hard-core 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 time-reversal 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 single-particle 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 Bose-Hubbard 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 well-controlled projection onto the flat-band subspace Huber and Altman (2010); Zhitomirsky and Tsunetsugu (2004, 2005); Tovmasyan et al. (2013); Takayoshi et al. (2013) in the weak-interaction regime, , and yields an effective low-energy Hamiltonian applicable to a wide filling range. This is analogous to the lowest-Landau 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 mean-field treatment predicts transitions to a non-uniform 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 many-body 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 real-space correlation function. In addition to the standard delta-function 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 time-of-flight imaging in the context of ultra-cold 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 ultra-cold 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 single-particle 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 Wigner-crystal ground states below the close-packing filling of loop states. In Section V, we study the quantum phases beyond the close-packing filling. We apply a flat-band 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 mean-field 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 Bose-Hubbard model on the kagome lattice subject to gauge flux [Fig. 1] is described by the Hamiltonian

 H=∑⟨r,r′⟩(|t|eiArr′b†r′br+h.% c.)+U∑rb†rb†rbrbr, (1)

where creates a single boson on the site labeled . The first term is the tight-binding Hamiltonian determining the band structure of the non-interacting bosons. We denote the hopping amplitude by to stress that it is positive (frustrated) 111In ultracold-atom 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 nearest-neighbor 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 Single-particle eigenstates

### iii.1 Gapped flat band in the presence of flux

We first discuss the tight-binding band structure for . In the absence of external flux, the single-particle 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 single-particle spectrum takes the typical Hofstadter butterfly form [Fig. 1]. For non-integer , the lowest band 222The 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, time-reversal (TR) symmetry is intact and the model can be realized with real-valued 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 tight-binding kagome Hamiltonian:

1. The energy of the single-particle eigenstate is exactly -2.

2. The elementary flat-band eigenstates are single-particle 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 (non-elementary) flat-band state can be composed as a linear superposition of loop states.

3. “Flux quantization”: any loop eigenstate encloses an integer number of flux quanta,

 ϕL=to0.0ptto7.499886pt\hss$↺$\hss∑loopArr′=∑loopϕ∈2πN. (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 zero-flux 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 flux-quantization 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 positive-hopping bonds and are identical across decorated negative-hopping 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 on-site boson repulsion due to the Hubbard term

 V=∑rVr=∑rUb†rb†rbrbr. (3)

Since the interaction is local, we note that any many-body state of the form

 |ψ⟩=∏m∈AL†m|0⟩ (4)

with single-particle occupation of a set of non-overlapping LOs is an exact ground state of the interacting system for filling 333Note: 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 interaction-induced 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 15-fold 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 444Note that the Wigner crystal state, is first discovered in the spin context, and termed as valence-bond 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 valence-bond 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 hard-core loop gas with macroscopic degeneracy determined by all possible configurations of non-overlapping 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 flat-band 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 ground-state 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 Flat-band 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 low-energy effective Hamiltonian by adapting the approach by Huber and Altman Huber and Altman (2010), consisting of a projection onto the subspace spanned by flat-band 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 weak-interaction 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 negative-hopping bonds [Fig. 5(a)]. Due to the unit-cell 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 left-dimer states only containing hexagons from right-dimer 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 mutually-orthogonal Wannier orbitals (WOs). As usual, there is not a unique set of WOs and different choices can vary significantly in their real-space localization. Since we will ultimately employ local-decoupling mean-field theory, it is particularly important to obtain well-localized WOs 555 Note that it is possible to construct symmetric but less compact localized orbitals as our flat-band basis. However, the mean-field 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 Gram-Schmidt procedure (see Appendix C for details). The results for two adjacent WOs are depicted in Fig. 5(b). The major part of the real-valued 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 power-law decay of WO amplitudes in the 0-flux 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

 w†j≡∑rwj(r)b†r, (5)

where the Wannier function gives the amplitudes of the dimer-type WO centered at position of the effective triangular lattice [Fig. 5(a)] on each site of the underlying kagome lattice. The flat-band projection corresponds to the inverse transformation

 b†r→∑jw∗j(r)w†j, (6)

where the Wannier states of the dispersive bands have been dropped as an approximation. Upon projection and switching to the grand-canonical ensemble, the effective Hamiltonian takes the form

 H→Heff=∑j(−2|t|−μ)w†jwj+∑ijklIijklw†iw†jwkwl, (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 real-valued (since the constructed Wannier functions are real-valued 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 on-site repulsion term has the largest strength. The next sub-leading terms come from two other types of effective interaction, namely, density-density 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 density-density 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 on-site repulsion strength .

Within the low-density regime (i.e., in the effective triangular lattice), we therefore employ a hard-core 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 density-density repulsion and assisted hopping, the only remaining interaction type is ring-exchange [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 density-density repulsion and assisted hopping.

### v.2 Mean-field theory.

While density-density repulsion favors density-wave 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 mean-field theory (MFT). We adopt the Gutzwiller approach Rokhsar and Kotliar (1991) and employ a product ansatz consistent with the hard-core constraint

 |ψMF⟩=∏j(fj,0+fj,1w†j)|0⟩, (8)

which decouples sites on the effective triangular lattice of WOs. The mean-field ansatz naturally captures the nematic Wigner crystal phase since it is a product of single-particle states with occupation of non-overlapping LOs (in this case approximated by WOs). Above close packing, mean-field solutions continue to break the symmetry due to the anisotropic nature of the Wannier orbitals.

To describe states with density-wave order such as the nematic Wigner crystal, we must allow for the dependence of the mean-field amplitudes on the spatial index . To obtain mean-field solutions, we decouple the effective Hamiltonian, replacing density-density interaction and assisted-hopping terms by

 VDD→ ∑i,j′2Idijni⟨n% j⟩, VAH→ ∑i,j′∑kIahijk⟨w†i⟩⟨wj⟩nk + ∑i,j′∑kIahijk(w†i⟨wj⟩+h.c.)⟨nk⟩. (9)

(We have verified that inclusion of ring-exchange does not lead to significant changes.) With this, we obtain a mean-field Hamiltonian , where depends on the mean-field 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 re-calculate order parameters until reaching self-consistency (see Appendix E for details).

For a range of chemical potentials, we calculate results for the mean filling , density-wave 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 density-wave 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 hard-core 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 density-wave order parameter remains nonzero and decays slowly, overall suggesting a second-order transition to a nematic supersolid [Fig. 6(c)] in which a fraction of the bosons condense on interstitial sites between the Wigner-crystal structure. Further on at , the density-wave order abruptly drops to zero, accompanied by a sudden increase in the superfluid order . This indicates a sudden melting of the Wigner-crystal structure and a first-order 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 non-monotonic 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 time-of-flight (TOF) imaging in the context of cold-atom 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 time-of-flight experiments

The microscopic signature, i.e., the ground-state 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 single-particle density matrix:

 ⟨nq⟩=1Nsite∑r1,r2eiq⋅(r1−r2)⟨b†r1br2⟩. (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 non-overlapping 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

 ⟨nq⟩NWC=1Nsite∑j∈A∑α,α′eiq⋅(rα−rα′)⟨0|Ljb†αbα′L†j|0⟩, (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 single-particle 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 delta-function 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 mean-field 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 nearest-neighbor hopping along the major axis ( direction) is zero in the nematic superfluid phase (within mean-field approximation). Only successive hopping along other directions will contribute to superflow in the direction. On the other hand, the large effective nearest-neighbor 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. Single-particle localized loop states can be combined to construct many-body 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 flat-band projection using dimer-shaped Wannier orbitals and subsequent mean-field 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 time-of-flight 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 strong-interaction regime.

###### Acknowledgements.
We are indebted to Steven Girvin, Ashvin Vishwanath, Tigran Sedrakayan, Eliot Kapit, Hakan Türeci, Tzu-Chieh 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 PHY-1055993. 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. DE-AC02-06CH11357.

## 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 positive-hopping tight-binding lattice Hamiltonian in the presence of an additional gauge field has the form:

 Htb=∑⟨r,r′⟩|t|eiArr′b†r′br+H.c.≡∑⟨r,r′⟩Trr′, (12)

where is the hopping operator on the nearest-neighbor bond . Assuming that there exist single-particle wavefunctions of the loop eigenstate type, they can be expressed as

 |ψL⟩=∑r∈Lψrb†r|0⟩, (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

 Htb|ψL⟩=∑⟨r,r′⟩∈LTrr′|ψL⟩+∑⟨r,r′⟩∈L,l∉L[Trl+Tr′l]|ψL⟩. (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:

 ∑⟨r,r′⟩∈L|t|[eiArr′ψrb†r′+e−iArr′ψr′b†r]|0⟩ (15) =−2|t||ψL⟩≡∑⟨r,r′⟩∈L−|t|[ψr′b†r′+ψrb†r]|0⟩.

The above equation leads to the following relation between the amplitudes of neighboring sites,

 ψr′=−ψreiArr′. (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 0-flux 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 single-valued, we derive the flux quantization condition for a loop eigenstate (property 3), namely

 ϕL=∑⟨r,r′⟩↺   Arr′= 2πn, n∈N. (17)

Now we consider the cancellation of the second sum in Eq. (14) which leads to caging, and get

 [Trl+Tr′l]|ψL⟩ = |t|[eiArlψrb†l+e−iAr′lψr′b†l+H.c.]|0⟩=0. (18)

This in turn leads to . Combined with Eq. (16), we get

 ϕΔ=Arl+Alr′+Ar′r=0. (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 (gauge-equivalent to -flux in the negative hopping model).

To understand one additional feature of the loop state, now we consider the gauge-invariant current operator on the bond (from site to site ), namely,

 Jrr′=2|t|i[ b†r′breiArr′−H.c.]. (20)

Its expectation value is given by

 ⟨ψL|Jrr′|ψL⟩=2|t|i[ψ∗r′ψreiArr′−H.c.], (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 time-reversal (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:

 H=∑⟨ij⟩[J⊥ij(SxiSxj+SyiSyj)+JzijSziSzj]−h∑jSzj. (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:

 H=∑⟨ij⟩[12J⊥ij(S+iS−j+S−iS+j)+JzijSziSzj]−h∑jSzj. (23)

From above, one can see that the transverse interaction can be written as a flip-flop term, which in the case actually corresponds to hopping of hard-core bosons. The term induces nearest neighbor interactions between hard-core bosons. If the magnetic field is sufficiently large, the ground state is the magnon vacuum . For the flat-band 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 on-site 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 valence-bond crystal phase (equivalent to the Wigner crystal) formed by non-overlapping hexagon loop magons has been found in Ref. Schulenburg et al. (2002) (earlier than its boson analog). This valance-bond crystal phase, in the spin-1/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 valence-bond crystal phase is not limited to spin-1/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 valence-bond 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 density-density interaction between magnons is repulsive and hence stabilizes the valence-bond 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 sign-tunable spin model, a promising candidate is the nitrogen-vacancy center array Cai (2013), although in that case the spin-spin interaction is not restricted to nearest-neighbors but has a power-law decay due to its dipole-dipole nature.

Besides the nematic Wigner crystal (nematic valence-bond crystal) phase, the nematic supersolid or superfluid phases may also occur in the spin systems. However, due to the hard-core 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 weakly-interacting 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 0-flux 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 2-dimensional 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 flat-band degeneracy). The operators , which create the two types of loops within each unit cell, can be represented by the original lattice boson operators as:

 L†1,R= bA†R−bB†R+bC†R−bE†R−bE†R−a2+bF†R−a1 −bF†R−a1+a2+bB†R+a2−bA†R+2a2+bC†R+a2, L†2,R= bB†R+a1−bB†R+a1−a2−bC†R−bC†R+a2+bD†R−bD†R+2a2 +bE†R+bE†R+a2−bF†R+bF†R+a2. (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 6-site 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 translationally-invariant superposition (with a particular wavevector ) of the two types of loop states respectively, i.e.

 L†s,k=∑Re−ik⋅RL†s,R. (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 6-component Bloch vector , where we have and The two Bloch vectors got from the two chosen loop states are represented as

 →u1,k= (1−e−i2k⋅R2,−1+e−ik⋅R2,1+e−ik⋅R2,0, →u2,k= (0,e−ik⋅R1−e−ik⋅(R1−R2),−1−e−ik⋅R2, (26) 1−e−i2k⋅R2,1+e−ik⋅R2,−1−e−ik⋅R2)T.

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:

 w†s,R=1√Nsite∑keik⋅R~L†s,k≡∑R′,lwls,R(R′)bl†R′, (27)

where the Wannier wavefunction is given by

 wls,R(R′)=1√Nsite∑keik⋅Rβls,ke−ik⋅R′. (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 real-valued 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 Gram-Schmidt process, i.e.,

 |γ1,k⟩=|β2,k⟩−|β1,k⟩⟨β1,k|β2,k⟩||β2,k⟩−|β1,k⟩⟨β1,k|β2,k⟩|, (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 flat-band 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:

 |α1,k⟩ =1√2[|β1,k⟩+eiθk|γ2,k⟩], |α2,k⟩ =1√2[|β2,k⟩+eiθk|γ1,k⟩], (30)

where and are the mutually-orthogonal 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 k-space, which ensures that the generated Wannier function is exponentially localized Kohn (1973) and hence is more compact.

The Bloch states generated from the Gram-Schmidt 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 square-wave 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 hard-core constraint. The types of terms are following:

(1) On-site repulsion:

 Vonsite=∑jU′w†jw†jwjwj, (31)

where is the largest energy scale in the effective Hamiltonian.

(2) Density-density repulsion:

 VDD=∑(i,j)2Idijw†iwiw% †jwj=∑i,j′Idijw†iwiw†jwj, (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 pair-hopping (involving two different sites):

 VPHa=∑(i,j)Ipij(w†iw†iwjwj+H.c.)=∑i,j′Ipij%w†iw†iwjwj, (33)

where .

(3b) Off-site pair-hopping (involving three different sites):

 VPHb= ∑(i,j)∑k≠i,j[Ipijk(w†iw†j+w†jw†i)wkwk+H.c.] = ∑i,j,k′Ipijk[w†iw†jwkwk+H.c.], (34)

where .

(4a) Assisted-hopping (involving three different sites):

 VAHa= ∑(i,j)∑k≠i,jIahijk(w†iwjw†kwk+H.c.) (35) = ∑i,j,k′Iahijkw†iwjw†kwk. (36)

Four terms (and their H.c.) in the effective interaction correspond to the pair , namely . Thus, we get ().

(4b) Assisted-hopping (involving only two different sites):

 VAHb= ∑(i,j)[Iahijj(w†iwj+H.c.)w†jwj+Iahjii(w†iwj+H.c.)w†iwi] = ∑i,j′Iahijj(w†iwj+H.c.)w†jwj. (37)

Here, we have , due to the presence of two terms of each type, e.g. .

(5) Ring-exchange interaction:

 VRing=∑(i,j),(k,l)IRingij,kl(w†iw†jw