Mott Insulator-like Bose-Einstein Condensation in a Tight-Binding System of Interacting Bosons with a Flat Band

We propose a new class of tight-binding systems of interacting bosons with a flat band, which are exactly solvable in the sense that one can explicitly write down the unique ground state. The ground state is expressed in terms of local creation operators, and apparently resembles that of a Mott insulator. Based on an exact representation in terms of a classical loop-gas model, we conjecture that the ground state may exhibit quasi Bose-Einstein condensation (BEC) or genuine BEC in dimensions two and three or higher, respectively, still keeping Mott insulator-like character. Our Monte Carlo simulation of the loop-gas model strongly supports this conjecture, i.e., the ground state undergoes a Kosterlitz-Thouless transition and exhibits quasi BEC in two dimensions.


I. INTRODUCTION
Tight-binding models of interacting particles with a flat band, i.e., a set of highly degenerate single-particle energy eigenstates, have been studied intensively over the decades. Flat-band systems do not only serve as idealized models of materials with a narrow band, but also provide a theoretical playground for investigating various collective phenomena arising from the interplay between particle motion and interactions. This is because the effect of interactions is magnified due to the flatness of the band. Such an approach was fruitful in the study of the origin of ferrimagnetism [1] and ferromagnetism [2][3][4][5] in the Hubbard model. See [6] for a review. For the formation of a Wigner crystal and the effect of the change in the density in bosonic systems with a flat band, see [7][8][9][10][11]. The recent proposal that the nearly flat band in twisted bilayer graphene supports superconductivity is intriguing [12].
In this paper we propose a new class of tight-binding systems of interacting bosons with a flat lowest band. The model is based on the construction in [13] of flat band Hubbard models, and can be regarded as a spinless version of the models studied in [14]. We write down the ground state of the model explicitly as in (5), and prove that it is the unique ground state. The expression suggests that the ground state is essentially different from that of a non-interacting system, and resembles that of a Mott insulator.
In spite of the simple expression, the property of the ground state is nontrivial and rich. By examining representations of the norm and the correlation functions in terms of a classical loop-gas model, we conjecture that the ground states may exhibit quasi off-diagonal longrange order (ODLRO) in two dimensions, and genuine ODLRO in three or higher dimensions. We present some results of Monte Carlo simulation of the loop-gas model, which strongly indicates that the two-dimensional model undergoes a Kosterlitz-Thouless (KT) transition and exhibits quasi ODLRO in its ground states.
We note that a class of states (without parent Hamiltonians) very similar to ours was proposed and examined in [15]. It was found that these states do not exhibit off-diagonal (quasi) long-range order.
It was found in a two-component system of bosons that one component may exhibit Bose-Einstein condensation (BEC) while the other is in the Mott insulating state [16]. This is different from our Mott insulator-like ground state (5), which consists only of the states created by theb † u operators. Note also that our ground state in the BEC phase, although solid-like, is not a supersolid since there is no spontaneous breakdown of translation symmetry [17][18][19]. It is a challenging problem to design a similar exactly solvable model that has supersolid ground states.
We stress that our ground states maintain Mottinsulating-like nature even when they exhibit (quasi) ODLRO. This is most clearly seen in the anomalously small particle number fluctuation observed in a specific setting. See section VI, in particular, Figure 8. This does not mean, however, that the ground states describe genuine Mott insulators. See section V.
It would be exciting if our exactly solvable model provides an example of a novel exotic phase of matter where Mott-insulator like nature and (quasi) ODLRO coexists. Although we are not able to give a definite conclusion at the moment, we find it rather likely that (unfortu- nately) our Mott-insulator like condensate is smoothly connected to ordinary Bose-Einstein condensate realized, e.g., in non-interacting systems. See section VI. We nevertheless believe that it is of essential importance that an entirely new class of exactly solvable models of strongly interacting bosons has been discovered. We hope that the present work opens a new direction in the research of quantum many-body systems.

II. THE MODEL AND THE EXACT GROUND STATE
Let (E, B) be a finite lattice, where E is the set of sites and B is the set of bonds. A bond is an unoriented segment that connects two distinct sites in E. We allow multiple bonds to connect the same pair of sites. At the center of each bond b ∈ B, we take a new site and denote it as u b . We let I be the set of all sites u b with b ∈ B, and consider the decorated lattice Λ = E ∪ I. We typically choose E to be the set of sites of the d-dimensional hypercubic lattice with periodic boundary conditions, and take p distinct bonds connecting a pair of neighboring sites. See Fig. 1 for the resulting decorated lattices.
We shall define a tight-binding model of bosons on Λ. We denote byâ † r andâ r the creation and annihilation operators, respectively, of a boson at site r ∈ Λ. They satisfy the canonical commutation relations [â r ,â s ] = 0 and [â r ,â † s ] = δ r,s for any r, s ∈ Λ. The number operator is defined asn r =â † râr . We denote by |Φ vac the state without any bosons, i.e., the unique normalized state such thatâ r |Φ vac = 0 for any r ∈ Λ. For each x ∈ E, we defined where ζ > 0 is a model parameter, and N (x) is the set of sites in I that are on the bonds connected to x. We consider the HamiltonianĤ =Ĥ hop +Ĥ int with ǫ(k) The single particle energy eigenvalue ǫ(k) as a function of the wave number vector k = (k1, k2) for the models defined on the square lattice (E , B) as in Figure 1. Here kj (j = 1, 2) runs between −π and π. There are flat bands (formed by theb states) with zero energy and a dispersive band (formed by thed states) with width 4pt. The band gap is equal to ζ 2 t. See Appendix A, in particular, (A10).

the hopping Hamiltonian
where t > 0, and the interaction Hamitlonian where U > 0. Our model is characterized by the three parameters ζ, t, and U . Note thatĤ hop is rewritten in the standard form asĤ hop = r,s∈Λ t r,sâ † râ s , where the amplitude t r,s describes hopping between the nearest and some of the next nearest neighbor sites, and on-site potentials. Note also thatĤ int describes on-site repulsive interaction only on sites in I.
The hopping HamiltonianĤ hop describes a tightbinding model with a flat band. To see this we definê for u ∈ I where x, y ∈ E are the sites connected by the bond corresponding to u. One can easily verify the orthogonality [d x ,b † u ] = 0 for any x ∈ E and u ∈ I. This means thatĤ hopb † u |Φ vac = 0 for any u ∈ I. Sincê H hop ≥ 0 and the statesb † u |Φ vac with u ∈ I are linearly independent, we see thatĤ hop has |I| independent singleparticle ground states (with zero energy). See Figure 2 and Appendix A for more about the single-particle energy eigenvalues.
Let us state our first theorem, which characterizes the ground state of the model. Theorem 1 .-Consider the above model with particle number N = |I|. For any ζ > 0, t > 0, and U > 0, the ground state ofĤ is unique, has vanishing energy, and is written as Note that, in the ground state (5), exactly one particle is associated with adjacent three sites including u ∈ I. This suggests that the particles are distributed almost uniformly over the lattice, and the state resembles the ground state of a "solid" or, more precisely, a Mott insulator. But the property of the ground state turns out to be much richer than a simple Mott insulator. In fact we shall argue in section V that the ground state does not correspond to a genuine Mott insulator when it exhibits (quasi) ODLRO.
The proof of Theorem 1 is an easy application of the technique developed in [4,5,13] for the (fermionic) Hubbard model.
Proof of Theorem 1 .-It follows from the definitions thatĤ|Φ GS = 0. SinceĤ ≥ 0, this proves that |Φ GS is a ground state. We only need to prove that it is the unique ground state.
We first note that any single-particle state can be expressed by a linear combination of the operatorsb † u with u ∈ I andd † x with x ∈ E. This follows by observing that b † u andd † x are all linearly independent, and there are exactly |I| + |E| = |Λ| operators. We thus see that any N particle state is written as a linear combination of the basis states |Φ n = { u∈I (b † u ) nu }{ x∈E (d † x ) nx }|Φ vac , where the "occupation number" n = (n r ) r∈Λ satisfies r∈Λ n r = N . Our proof is based on the standard argument for frustration free Hamiltonians. Let |Ψ be an arbitrary state with N = |I| particles such thatĤ|Ψ = 0, and expand it as |Ψ = n α n |Φ n . Noting thatd † xdx ≥ 0 and n u (n u − 1) = (â † u ) 2 (â u ) 2 ≥ 0, we see thatd † xd x |Ψ = 0 and (â † u ) 2 (â u ) 2 |Ψ = 0 for all x ∈ E and u ∈ I. These relations further imply thatd x |Ψ = 0 for all x ∈ E and (â u ) 2 |Ψ = 0 for all u ∈ I. The first condition implies that α n = 0 whenever n x = 0 for some x ∈ E. Thus the ground state contains only theb-states. Then the second condition implies that α n = 0 whenever n u ≥ 2 for some u ∈ I. This means that |Ψ is a constant multiple of |Φ GS .

III. GROUND STATE PHASE TRANSITION
We first derive an exact loop-gas model representation of the ground state. Further details of the derivation are given in Appendix B. Note that in this quantumclassical correspondence, unlike in the standard correspondence via path-integral, a d-dimensional quantum state is mapped to a d-dimensional classical system. This is a peculiar point about our model, and is analogous to the loop-gas representations of the Affleck-Kennedy-Lieb-Tasaki (AKLT) model [20] and the Kitaev model [21]. Note first that theb-operators satisfy the commutation relations where Here u ≈ v indicates that the bonds corresponding to u and v connect the same pair of sites, and u ∼ v indicates that the bonds for u and v share a single common site. By repeatedly using (6), one finds that the normalization factor Φ GS |Φ GS is represented as [15] where the sum is over all possible sets L = {ℓ 1 , . . . , ℓ n } with n = 0, 1, 2, . . . of oriented loops. See Fig. 3. By an oriented loop of length m, we mean a sequence ℓ = (u 1 , . . . , u m ) of m distinct sites in I such that u j ≈ u j+1 or u j ∼ u j+1 for j = 1, . . . , m, where we set u m+1 = u 1 .
Here we identify the new sequence obtained by the shift u j → u j+1 for j = 1, . . . , m with the original sequence. This means that, for u 1 , u 2 ∈ I such that u 1 ∼ u 2 , there is a unique loop, ℓ = (u 1 , u 2 ), of length two, while for We also assume that, in any set L = {ℓ 1 , . . . , ℓ n }, no loops share a common site. To take into account the factor 2 in (6), we interpret u, v ∈ I with u ≈ v as being connected via two distinct paths, and properly over-count loops according to this interpretation. Finally we wrote |L| = n j=1 |ℓ j |, where |ℓ| denotes the length of a loop ℓ. See Appendix B for details.
We can derive a similar representation for correlation functions. Noting thatâ u |Φ GS = ζ √ β( w =ub † w )|Φ vac for u ∈ I, we find that where L is again summed over sets of oriented loops, and ω is summed over all self-avoiding walks connecting u to v, i.e., a sequence ω = (u 0 , . . . , u m ) of m + 1 distinct sites in I such that u 0 = u, u m = v, and u j ≈ u j+1 or u j ∼ u j+1 for j = 0, . . . , m − 1. The length of the walk is defined as |ω| = m. The prime in the sum indicates that ω and loops in L do not share common sites. By combining (8) and (9), we see that the off-diagonal correlation is represented as (10) The correlation function â † râs GS with general r, s ∈ Λ has a similar (but slightly more complicated) representation, and should behave almost similarly as â † vâu GS , especially when the distance between r and s is large.
The above summations over loops and walks are in general nontrivial and hardly evaluated explicitly. When the basic lattice (E, B) is a chain with bonds connecting nearest neighbor sites, i.e., d = p = 1 in the notation of section II, one can evaluate the summations by using, e.g., the transfer matrix method to show for a long enough chain that where C = ζ 2 β/ 1 + 4β 2 . The correlation always decays exponentially, and the ground state is disordered. Let us turn to the models in higher dimensions. When β is small (or ζ is large), we can easily prove that the corresponding loop-gas model is in the disordered phase, where configurations with only small loops are dominant. This corresponds to a disordered ground state that describes a Mott insulator. To see this, we relax the constraint L ∩ ω = ∅ in (9) to get This implies where Ω u,v (n) is the total number of self-avoiding walks of length n that connect u and v, and dist(u, v) is the minimum length of such walks. For u ∈ I, let s u and d u be the numbers of v such that u ∼ v and u ≈ v, respectively. We let ν be a constant such that s u + 2d u ≤ ν +1 for any u ∈ I. Then we have Ω u,v (n) ≤ (ν +1)ν n−2 .
(Proof: There are at most (ν+1) choices for the first step, and at most ν choices for each of the following n − 2 steps. There is no choice in the final step because the walk ends at v.) This bound, with (13), implies that the correlation decays exponentially for sufficiently small β, or, equivalently, sufficiently large ζ.
When the dimension is larger than one and β is sufficiently large, on the other hand, it is expected that the loop-gas model is in the percolating phase where a macroscopically large loop (or a walk) appears. The transition can be characterized by the divergence of the "static structure factor" If the critical value of β for the transition is less than 1/2, which is the upper bound for β, the ground state undergoes a phase transition. Let us first examine the standard large-d (or meanfield) approximation, in which one fixesβ = νβ and lets d ↑ ∞ (and hence ν ↑ ∞). It is well-known that, for β < 1, one can neglect contributions from loops in this limit. See Appendix C. Then S is evaluated for anyβ < 1 as where ω is summed over all self-avoiding walks (with an arbitrary length) that starts from u. We here noted that the number of walks is almost equal to ν n . The structure factor S diverges as β approaches the critical value ν −1 .
The existence of a phase transition is also suggested by a more careful examination of the loop-gas model. Recall that, in (8) or (9), every loop of length three or more is summed exactly twice with different orientations. This is equivalent to considering unoriented loops, but with an extra factor 2 for each loop. We thus see that our loop-gas model resembles that obtained from the hightemperature expansion of an O(2) symmetric classical ferromagnetic spin system at a finite temperature. The O(2) symmetry corresponds to the U(1) phase symmetry in the original quantum system. From this analogy we conjecture that in dimensions three or higher the ground state |Φ GS may exhibit BEC (where â † vâu GS has longrange order) while in two dimensions it may exhibit a quasi BEC (where â † vâu GS shows power law decay), both at sufficiently large β.
To prove the existence of BEC is extremely difficult, if not impossible, since we need to show the breakdown of a continuous symmetry. See [22] and references therein for rare cases where the existence of BEC can be established rigorously by using the method of reflection positivity. (A readable account of the reflection positivity method can be found in [6].) In the present model, we have a rigorous result only for the system defined on a tree, where one can apply standard techniques (see, e.g., [23,24]). See Appendix D.

IV. NUMERICAL EVIDENCE OF A KOSTERLITZ-THOULESS TRANSITION
To see whether the ground state exhibits a phase transition, we carried out a Monte Carlo simulation of the corresponding loop-gas model by using the worm algorithm [25]. We focus on the two-dimensional models constructed from the L × L square lattice with bond number p = 1 and 2 and with periodic boundary conditions. See Fig. 1.
A central quantity that we measured is the static structure factor S defined in (15). When the correlation â † vâu GS becomes long-ranged, the structure factor should diverge as a function of L. Specifically, S = O(L d ) if the system possesses ODLRO, and S = O(L d−η ) if it has quasi-ODLRO characterized by the correlation-decay exponent η, while S = O(1) in the disordered phase. In particular, η = 1/4 at the KT transition point with the multiplicative logarithmic correction of (log L) 1/8 [26].
To capture more specific features of the KT transition, we also measured the helicity modulus, which is proportional to the super-fluid density. It can be computed as the Monte Carlo average of the squared total winding number of all loops [15,27], where the summation is over all loops in the loop-gas configuration L and w l is the winding number of the loop l around the horizontal direction of the lattice. As we pass the KT transition point entering the quasi-ODLRO phase, this quantity is expected to show a discontinuous jump from zero to the universal value 2/π, and keep increasing afterwards. Let us start from the simplest model with p = 1. The top panel of Fig. 4 shows the β dependence of the static structure factor S for varying system sizes up to L = 64. Clearly we do not observe any singularity in the region β ≤ 1/2, where the classical-quantum correspondence is valid. The data, however, seem to indicate that the system, as a classical loop-gas model, has a phase transition beyond β = 1/2.
The bottom panel of Fig. 4 shows the β dependence of the helicity modulus Υ. It can be seen from the figure that Υ goes from zero to some finite value greater than the expected jump, indicating that a KT transition takes place near β ∼ 1. This is in contrast to the case studied in [15]. We shall recall, however, that the loop-gas model with β ∼ 1 does not correspond to a quantum mechanical ground state.
Having observed that the model with p = 1 does not exhibit a transition (in the range β < 1/2 that is relevant for us), let us focus on the model with p = 2. Figure  5 shows the β dependence of the structure factor (top) and the helicity modulus (bottom) for varying system sizes up to L = 64 in the loop-gas model corresponding to p = 2. Clearly the static structure factor S shows a diverging behavior for large but not too large values of β. The helicity modulus Υ also shows an increase from zero to some finite value greater than the expected universal jump π/2, strongly suggesting the presence of a KT transition in the physically meaningful region β < 1/2.
Further evidence of a transition can be obtained from the size-dependence of the structure factor, shown in Fig. 6. For small values of β, before the slope of the curve reaches 7/4, it is not a straight line but bends and starts converging to a finite value, as seen for β = 0.15. As we increase β, the curve becomes a straight line with slope approximately equal to 7/4 (β = 0.20 in the figure). Further increasing of β brings the slope of the straight line greater than 7/4 (β = 0.5 in the figure). To estimate the slope of the straight part of the curve more quantitatively, we analyzed the data with or without the expected logarithmic correction, i.e., fit the form L 2−η or (log L) 1/8 L 2−η to the data regarding the exponent η as a fitting parameter, although according to [28] it is not clear whether including such a logarithmic correction would improve the estimates or not. The inset shows the result of the fitting with or without the logarithmic factor. From this result, we may conclude that the KT transition takes place at β = 0.20(1), which is consistent with the universal jump in the bottom panel of Fig. 5 considering the expected slow (∼ (log L) −1 ) convergence of Υ at the transition point.
To get an additional evidence of a KT transition, we examined this size-dependence of the helicity modulus near the transition point in more detail. As shown in Fig.7, we plot Υ as a function of 1/ log L around the critical value of β. At β = 0.2, the value close to the critical point estimated from the size dependence of the structure factor, the data is consistent with the linear convergence in 1/ log L to the universal jump, as predicted in [29] as a characteristic size-dependence for the KT trnasition. This may be taken as another evidence suggesting the KT nature of the transition.

V. THE MODEL WITH OTHER PARTICLE NUMBERS
Let us briefly discuss the properties of the models with different particle numbers, and related issue about the charge gap.
When the number of particles, N , is less than |I|, we see that the space of the ground states is spanned by |Φ S = ( u∈Sb † u )|Φ vac where S is an arbitrary subset of I such that |S| = N . The ground states show macroscopic degeneracy, which should be immediately lifted when a generic infinitesimal perturbation is added to the Hamiltonian.
The case with N > |I| is difficult, and we have almost no exact results. We can nevertheless construct a class of exact energy eigenstates, which may be regarded as new examples of quantum many-body scars [30,31]. Assume that the lattice (E, B) is connected and bipartite. To be precise we say that (E, B) is bipartite if there is a decomposition E = E + ∪ E − such that two sites x, y may be connected by a bond in B only when x ∈ E + , y ∈ E − or x ∈ E − , y ∈ E + . Let us definê Then, for any subset S ⊂ I and n = 1, 2 . . ., it is easily shown that the state which has |S| + n particles, is an energy eigenstate with eigenvalue E n = ntζ 2 . See Appendix E. Let us denote by E GS N the ground state energy of the model with N particles, and define the charge gap (or the jump in chemical potential) as It is believed that the charge gap provides a simple criterion for conductivity in the sense that the ground state is insulating if ∆ N is positive and of order 1 [32].
Since we have E GS N = 0 for N ≤ |I| in the present model, we see that E GS |I|+1 is nothing but the charge gap ∆ |I| , which is directly relevant to the property of the model with |I| particles. We conjecture that, for β < β c , the charge gap ∆ |I| is strictly positive and our ground state describes a Mott insulator, while, for β > β c where the ground state exhibits (quasi) BEC, ∆ |I| vanishes in the infinite volume limit according to the theorem in [33]. Therefore the ground state, although Mott insulator-like, is not a genuine Mott insulator.

VI. DISCUSSION
We proposed a new class of exactly solvable models of interacting bosons with a flat band, and argued that the Mott insulator-like ground states may exhibit (quasi) BEC. The conjecture is supported by the strong numerical evidence that the two-dimensional model exhibits a KT transition. The properties of the three-dimensional models remain to be investigated.
We believe it important that a new exactly solvable model that exhibits (or, that is conjectured to exhibit) nontrivial condensation phenomena has been discovered. It is also interesting to investigate the possibility of similar models of electrons, which should exhibit superconductivity.
It may be counterintuitive that our exact ground state (5), which consists of bosons almost localized at each u ∈ I, exhibits ODLRO. One should note however that the operatorb † u creates a coherent superposition of the three states in which a particle is at x, y, and u. The coherence "propagates" in the system thus generating offdiagonal correlation, which may be short-ranged or longranged [15]. At least mathematically, the situation is parallel to that for the long-range Néel order in the exact valence-bond ground states of the AKLT model in high dimensions [20,24,34].
Indeed this point is related to the fact that the states proposed in [15] only have short-ranged off-diagonal correlation while our ground state on a suitable twodimensional lattice exhibits quasi ODLRO. The basic difference is not of qualitative but of quantitative nature, namely, the commutator [b u ,b † v ] of neighboring sites, which give the basic parameter β, and the coordination number ν + 1 can be larger in our models compared with that in [15].
It is worth recalling that a two-dimensional system of interacting bosons at zero temperature generically exhibits genuine ODLRO rather than quasi ODLRO. We believe that the present model exhibits only quasi ODLRO because of its peculiar ground state structure (5) expressed only in terms of local bosonic operatorŝ b † u . This is consistent with the fact that the ground state of the present two-dimensional model is represented in terms of a classical statistical mechanical model (i.e., the loop-gas model) in two dimensions, while a ground state in two-dimensional quantum system generally corresponds to a three-dimensional classical system. We conjecture that the observed quasi ODLRO will be immediately elevated into genuine ODLRO when the model is perturbed. See, e.g., [35] for the discussion of a similar behavior in the ground state correlation function of the one-dimensional AKLT model.
We have repeatedly stressed that our exact ground state (5) has a Mott-insulator like character since it is generated by local operatorsb † u . Its peculiar property is most clearly seen if one considers a special lattice (E, B) obtained by connecting two arbitrary standard lattices by a single bond. See Fig. 8. In the ground state of the model defined on the corresponding decorated lattice Λ, the number of particles in each of the two sublattices is almost constant. (To be precise, it fluctuates only by one.) Note that this is also true when we properly choose the lattice so that the ground state exhibits (quasi) ODLRO. Such a ground state with (quasi) ODLRO and vanishing particle number fluctuation is quite exotic since (quasi) ODLRO is usually accompanied by large density fluctuation. We should not, however, jump to the conclusion that a novel exotic phase of matter has been discovered. It is possible that the zero fluctuation is a singular property of the exactly solvable model, and normal large fluctuation is recovered once the model is perturbed. For the moment we believe that this (less exciting) scenario is plausible. We nevertheless stress that the discovery of models in which anomalously small particle number fluctuation and (quasi) ODLRO may coexist indicates that strongly interacting systems of bosons may exhibit unexpectedly rich behavior.  Appendix A: Single-particle energy eigenvalues Let us examine the band structure determined by the hopping Hamiltonian (2). We first recall that the stateŝ b † u |Φ vac with u ∈ I andd † x |Φ vac with x ∈ E form a basis of the single-particle Hilbert space. Since the stateb † u |Φ vac has zero energy, we focus on states spanned byd † x |Φ vac . Here, for simplicity, we assume that any pair x, y ∈ E is either connected by p distinct bonds in B or not connected at all. We write x ∼ y when x and y are connected by p bonds, and denote by E(x) the set of y ∈ E such that x ∼ y. Recalling the definitiond we see that where z x = |E(x)| is the coordination number (i.e., the number of neighboring sites) of the original lattice (E, B). Note that |N (x)| = pz x . The commutation relation, along with the definition (2) ofĤ hop , implies Consider and arbitrary state spanned byd operators where ϕ x ∈ C. By using the commutation relation (A3) andĤ hop |Φ vac = 0, we see that where we switched the roles of x and y to get the final expression. Thus the Schrödinger equation which is nothing but the standard tight-binding Schrödinger equation with hopping pt and on-site potential (ζ 2 +pz x )t. Now we suppose that (E, B) is the d-dimensional L × · · · × L hypercubic lattice with periodic boundary conditions, i.e., and x ∼ y when |x − y| = 1. The coordination number is z x = 2d for all x. Let us take the standard plane wave ϕ (k) x = e ik·x with k · x = d j=1 k j x j and k ∈ K, where the set of wave number vectors is Substituting ϕ (k) x into (A7), one readily confirms that this is an energy eigenstate with eigenvalue See Figure 2 for the case with d = 2.
Since there isb u but nob † u , we must have [b u ,b † u1 ] with u ∼ u 1 to have a nonzero contribution. This implies that we also need [b u1 ,b † u2 ] with u 1 ∼ u 2 , and so on. This is terminated only when we have We have obtained the contribution from the random walk. The contribution from loops can be derived as in the above.
We now turn to a general class of models where we need to take into account the factor 2 that appears in the commutation relation [b u ,b † v ] = 2β for u, v ∈ I such that u ≈ v. This can be of course done by redefining the weights for loops and walks, but there is a more elegant way of incorporating the factor. We use the same definitions for the weights, but regard in general that sites u and v with u ≈ v are connected by two distinct paths.
To see that this works, it (almost) suffices to consider the simplest model constructed from the lattice E = {x, y} with only two sites where x and y are connected by two distinct bonds. We then let u and v be the sites at the center of these bonds. We thus have u ≈ v. In this case we see from (B1) that The factor 4 is reproduced as in Fig. 9. The configurations of loops in these general models can be much more complicated than Fig. 3 in the main text. In the model with d = 2 and p = 2 depicted in Fig. 1 (b), for example, the loops are defined on a lattice similar to Fig. 3, but each white site should be replaced by a pair of white sites connected by two distinct paths, and a path connecting two white sites should be replaced by four paths connecting two pairs of white sites.
x y u v FIG. 9. The graphs corresponding to the normalization factor ΦGS|ΦGS in the simplest model with only two u, v ∈ I such that u ≈ v. 10. The first three generations of the tree with branching number three. The site at the root is denoted as o.

Appendix C: Large-d approximation
In the large-d approximation, we can neglect the possibility that a random walk accidentally intersect its trajectory. This means that the number of length n walks is given by ν n . We used this estimate in (16). Note that we here do not make distinction between ν and ν + 1.
Let us see why we can neglect the contribution from loops. Let Ω(n) be the number of loops that contain a given site in I. Since there are at most (ν + 1) choices for the first step, at most n choices for each of the following n − 2 steps, and no choices in the final step, we have Ω(n) ≤ (ν + 1)ν n−2 . (C1) Thus the total contribution from all the loops containing a site is bounded from above by This vanishes as ν ↑ ∞ provided that the sum converges. (Remark: To be precise, to fix a site and sum over all the loops containing it is not a proper way of evaluating the summations in (10). But it gives a correct order estimate in terms of ν.) Appendix D: Bose-Einstein condensation in the model on a tree We shall study the model defined on a tree, and show that the ground state exhibits Bose-Einstein condensation for sufficiently small ζ.
Let (E n , B n ) be the regular n-generation tree with branching number three, as depicted in Fig. 10. As in the main text, we define the corresponding set of sites I n , and consider the model of interacting boson on E n ∪ I n . Our goal is to show that the ground state exhibits spontaneous symmetry breaking associated with BEC. We note that the model on the tree with branching number two does not exhibit a phase transition.
A standard method to test for the existence of spontaneous symmetry breaking is to impose boundary conditions that explicitly favor certain order, and see if the effect of the boundary conditions survives in the infinite volume limit. In the case of ferromagnetic spin systems, this is done by enforcing spins at the boundary to point in a certain fixed direction. In the case of BEC, the corresponding procedure is to replaceb † u at the boundary by α + γb † u with some nonzero α, γ ∈ C and then to examine the expectation value of the annihilation operator deep inside the tree. In the ferromagnetic Ising model, for example, it is known that the same procedure exactly recovers the result of the Bethe approximation [23]. See also [24] for a treatment of a quantum spin state.
To be precise, let ∂I n be the set of sites at the boundary in I n . See Fig. 11. We then consider the ground state with plus boundary conditions defined as Note that we have chosen α = γ = 1 for simplicity. We are interested in the expectation value especially in its limiting value as n ↑ ∞, where o denotes the site at the root of the tree (E n , B n ). As in the main text, we develop loop-gas representations for Φ + GS |Φ + GS and Φ + GS |â o |Φ + GS . Reflecting the special geometry of I n , the representations contain loops of length two, three, or four, but not larger. Apart from these loops the representations contain random walks that starts from a site in ∂I n and ends in another site in ∂I n . Note that the symmetry breaking boundary terms play the roles of sources and sinks of the walks. Of course the loops and the walks should satisfy the site-avoiding conditions. See Fig. 12 (a). These are all contributions for the representation for Φ + GS |Φ + GS , and we have The representation of Φ + GS |â o |Φ + GS must contain a random walk that starts from the site o ′ , the root of I n , and ends at a site in ∂I n . See Fig. 12 Following the standard procedure for models on a tree [23,24], we shall evaluate these sums by using exact recursion relations. We define four sums X n , Y n , Z + n , and Z − n of loops and walks as in (D3) and (D4) with different condition on the site o ′ at the root of I n . In X n we sum over all configurations where no loop or walk touching o ′ . In Y n we sum over all configurations where there are two segments (which are part of a loop or a walk) touching o ′ . In Z + n (resp. Z − n ) we sum over all configurations where there is exactly one segment (which is a part of a walk) coming into (resp. going out of) o ′ . Note that we have X 0 = 1, Y 0 = 0, and Z ± 0 = 1. We see that Φ + GS |Φ + GS = X n + Y n , Φ + GS |â o |Φ + GS = √ β Z − n , and hence where we defined y n = Y n X n , z n = Z − n X n . (D6) x ∈ E − , y ∈ E + . For simplicity we shall assume that a pair of sites x, y ∈ E is either connected by p bonds in B or not connected at all. But the following result about exact energy eigenstates is valid without this restriction. As in section V, we defineD Then we shall prove, for any subset S ⊂ I and n = 1, 2, . . ., that the state satisfiesĤ |Φ S,n = ntζ 2 |Φ S,n .
Thus |Φ S,n is an exact energy eigenstate with particle number N = |S| + n and energy E S,n = ntζ 2 . Note that we can expressD 0 in terms of thed-operators aŝ Then, by using (A3), we find We note in passing that tζ 2 is the minimum energy among single-particle energy eigenstates that are orthogonal to the flat band (i.e., the states generated by theb † operators), and that the state generated byD † 0 is the unique energy eigenstate with energy tζ 2 . This fact is obvious from the dispersion relation (A10) for the models on the hypercubic lattice, but can be proved in general.
Note that by construction the operatorD † 0 does not containâ † u for any u ∈ I. This immediately implies that It is clear from (E5) and (E6) thatĤ hop |Φ S,n = ntζ 2 |Φ S,n andĤ int |Φ S,n = 0. We thus find that |Φ S,n is an eigenstate ofĤ =Ĥ hop +Ĥ int with eigenvalue E S,n = ntζ 2 .