Spin and charge density waves in the Lieb lattice
SSpin and charge density waves in the Lieb lattice
J. D. Gouveia, R. G. Dias
Departamento de F´ısica, I3N, Universidade de Aveiro, Campus de Santiago, Portugal (Dated: May 8, 2015)We study the mean-field phase diagram of the two-dimensional (2D) Hubbard model in the Lieblattice allowing for spin and charge density waves. Previous studies of this diagram have shown thatthe mean-field magnetization surprisingly deviates from the value predicted by Lieb’s theorem [1] asthe on-site repulsive Coulomb interaction ( U ) becomes smaller [2]. Here, we show that in order forLieb’s theorem to be satisfied, a more complex mean-field approach should be followed in the caseof bipartite lattices or other lattices whose unit cells contain more than two types of atoms. In thecase of the Lieb lattice, we show that, by allowing the system to modulate the magnetization andcharge density between sublattices, the difference in the absolute values of the magnetization of thesublattices, m Lieb , at half-filling, saturates at the exact value 1 / U , as predictedby Lieb. Additionally, Lieb’s relation, m Lieb = 1 /
2, is verified approximately for large U , in the n ∈ [2 / , /
3] range. This range includes not only the ferromagnetic region of the phase diagram ofthe Lieb lattice (see Ref. 2), but also the adjacent spiral regions. In fact, in this lattice, below orat half-filling, m Lieb is simply the filling of the quasi-flat bands in the mean-field energy dispersionboth for large and small U . I. INTRODUCTION
Despite intense research in the last few decades, the2D Hubbard model in the square lattice has remained anopen theoretical problem in the field of the strong corre-lated systems [3, 4]. Although it is known that at half-filling, the spin dynamics of the 2D Hubbard model isdescribed by the Heisenberg antiferromagnetic exchangeterm [5], there is no consensus regarding the ground statemagnetic phase diagram of this model. In fact, even atthe mean-field (MF) level, depending on the magneticphases allowed, different authors obtain different dia-grams for the square lattice [6]. The traditional order-ings are ferromagnetism, antiferromagnetism and para-magnetism [7–11]. Later, spin spiral phases, a general-ization of the previous three, were introduced [12]. TheMF phase diagram became even more complex with theconsideration of spatial phase separation [13–15].The Hubbard model in decorated 2D lattices has alsobeen extensively studied, motivated by the search formetallic (flat-band) ferromagnetism. These decoratedlattices fall into three categories: the Lieb [1], Mielke[16] and Tasaki lattices [17]. All of these lattices share acommon feature: the presence of flat bands in the energydispersion relation. In the particular case of the Lieb’slattices, the flat bands result from the topology of the lat-tice, while in the case of Mielke and Tasaki lattices, theflat bands reflect longer-range transfer integrals in thesystem. One of the most representative examples of dec-orated lattices is the Lieb lattice, which can be obtainedfrom the 2D square lattice, for example, by inserting anextra atom between every two nearest-neighbours (seeFig. 1a). Each unit cell (shaded rectangle in Fig. 1a) has3 atoms, one of each kind: A, B and C, whose relativeoccupation is depicted in Fig. 1b, in the limit of no in-teractions. Fig. 1c shows the energy bands of the Lieblattice in this limit. Examples of materials whose struc-ture resembles the Lieb lattice include La − x Sr x CuO and YBa Cu O , two well-known high- T c superconduc-tors with weakly coupled CuO planes [18, 19].A theorem by Lieb [1] states that, in the particularcase of bipartite lattices (i.e., lattices with two sublat-tices, A and B, such that each site on sublattice A hasits nearest neighbors on sublattice B, and vice versa),the ground state is ferromagnetic at half-filling ( n = 1,or one electron per lattice site), as long as the number ofatoms of each sublattice is different. However, for exam-ple in the case of the Lieb lattice (a line-centered squarelattice [20]), this ground state should be identified withferrimagnetism [21]. In fact, although each sublattice isindeed ferromagnetic, there is antiferromagnetic orderingbetween every pair of nearest neighbours [2] (see Fig. 1a).The magnetic phase diagram of the Hubbard modelin the Lieb lattice was recently studied by us [2]. Weshowed that the mean-field magnetization per unit cellat half-filling, m Lieb = ( | m B | + | m C | − | m A | ) /
2, assum-ing the particle density ( n ) of the tight-binding limit (seeFig. 1b) and the same magnetization ( m ) in the wholelattice, surprisingly deviates from the value predicted byLieb’s theorem [1] as the on-site repulsive Coulomb inter-action ( U ) becomes smaller, although these two assump-tions are common in mean-field studies [2, 22–25]. Lieb’stheorem predicts that the magnetization per unit cell, m Lieb , is 1 / U at half-filling. Fig. 2 shows boththe mean-field phase diagram of the Hubbard model inthe Lieb lattice, and the value of m Lieb , using the mean-field results from Ref. 2. For the results to agree withLieb’s theorem, one should have m Lieb = 1 / U .Although in the strong coupling limit ( U (cid:29) t ), the mean-field result satisfies Lieb’s theorem, it is far from correctnear the tight-binding limit ( U = 0).In this manuscript, we study the magnetic phase di-agram of the Lieb lattice allowing for different averageoccupations ( n A , n B , and n C ) and magnetization am-plitudes ( m A , m B , and m C ) in each sublattice. Wefind that with these new considerations, Lieb’s relation, a r X i v : . [ c ond - m a t . s t r- e l ] M a y B C A (a) n n A , n B = n C n A n B = n C (b) x y ε k (c) FIG. 1. (a) The Lieb lattice is a line-centered square lattice,comprising three sublattices, A, B, and C and having oneatom of each type in a unit cell. The circles represent atomicnuclei and the arrows represent spins. At half-filling, the Lieblattice is ferromagnetic within each sublattice but antiferro-magnetic overall. Moreover, at half-filling, the total spin perunit cell is 1 /
2, as predicted by Lieb [1]. (b) Plot of the tight-binding ( U = 0) particle density of each sublattice of the Lieblattice, A, B or C, as a function of the total particle density.(c) Plot of the tight-binding dispersion relation ε ( k x , k y ) ofthe Lieb lattice. The flat band is made up entirely of B andC orbitals. m Lieb = 1 /
2, is satisfied for any U at half-filling, andsatisfied approximately for large U , in the n ∈ [2 / , / m Lieb is no longer the unit cell magnetization,but gives the difference in the absolute values of the mag-netization of the sublattices). An important point in thisresult is that finite m Lieb reflects the existence of quasi-flat bands in the mean-field energy dispersion. Thesequasi-flat bands are present, not only for small U , butalso for large U . In fact, below or at half-filling, m Lieb isapproximately the filling of the flat bands both for largeand small U .Our mean-field approach follows that of Bach andPoelchau [26] (see also Refs. 13 and 27). The formal-ism with further mathematical details can be found inRef. 28.The organization of this paper is as follows. We beginby presenting some key results for the Hubbard modeland mean-field. Secondly, we revisit the tight-bindinglimit of the Lieb lattice. We then proceed to adding elec-tronic interactions to the Hamiltonian and calculatingits mean-field counterpart. Finally, we show our results,discuss their meaningfulness, compare them to analyti- cal calculations and take conclusions. In Appendix A, we U m L i e b (a) (b) FIG. 2. (a) Mean-field magnetic phase diagram of the Hub-bard model in the Lieb lattice, obtained in Ref. 2 and (b)difference in the absolute values of the magnetization of thesublattices, m Lieb = ( | m B | + | m C | − | m A | ) /
2, of the Lieb lat-tice at half-filling, using the mean-field results from the samereference, where the on-site magnetization is assumed to bethe same on all sublattices. briefly outline the derivation which leads to the results insection II, and in Appendix B, we explain in more detailour method of calculating saddle points, which we usedto obtain the results in section V.
II. MEAN-FIELD METHOD FOR THEHUBBARD MODEL
In this section, we adapt a key result presented as The-orem 4.14 in Ref. 28. This derivation was first done byLieb and his colaborators [27] and simplified by Bach andPoelchau [26]; in Appendix A we show an adaptation ofthe derivation in Ref. 26. Alternative approaches can befound in Refs. 25 and 29. The result consists of thefollowing abridged derivation.We begin with the Hubbard Hamiltonian, given by H = t (cid:88) (cid:104) x,y (cid:105) ,σ c † x,σ c y,σ + U (cid:88) x ˆ n x, ↑ ˆ n x, ↓ . (1)Here, t is the hopping parameter between nearest neigh-bours and c † x,σ ( c y,σ ) is the creation (annihilation) oper-ator of an electron on site x ( y ) with spin σ = ↑ , ↓ . Theletters x and y denote lattice sites, (cid:104) x, y (cid:105) stands for near-est neighbors and U is the on-site repulsive. The totalnumber of particles is N .With the intent of finding the mean-field Helmholtzfree energy, F HF , associated to this Hamiltonian (notethat we use F because we work with fixed number ofparticles; if we worked with fixed chemical potential, wewould use the grand canonical potential instead), we firstreplace the interaction term, ˆ n x, ↑ ˆ n x, ↓ , with the Hartreeand Fock terms,ˆ n x, ↑ ˆ n x, ↓ → ˆ n x, ↑ (cid:104) ˆ n x, ↓ (cid:105) + (cid:104) ˆ n x, ↑ (cid:105) ˆ n x, ↓ − (cid:104) ˆ n x, ↑ (cid:105)(cid:104) ˆ n x, ↓ (cid:105) − c † x, ↑ c x, ↓ (cid:104) c † x, ↓ c x, ↑ (cid:105) − (cid:104) c † x, ↑ c x, ↓ (cid:105) c † x, ↓ c x, ↑ + (cid:104) c † x, ↑ c x, ↓ (cid:105)(cid:104) c † x ↓ c x, ↑ (cid:105) . (2)Replacing the averages by the mean-field parameters (cid:126)m and n (see Appendix A for details), it follows that themean-field Hamiltonian, H ( (cid:126)m, n ), corresponding to the Hamiltonian in Eq. 1 is H ( (cid:126)m, n ) = t (cid:88) (cid:104) x,y (cid:105) ,σ c † x,σ c y,σ + U (cid:88) x (cid:20) (cid:0) (cid:126)m x − n x (cid:1) + 12 (cid:16) n x ˆ n x − (cid:126)m x · ˆ (cid:126)s x (cid:17)(cid:21) , (3)where ˆ (cid:126)s x and ˆ n x are the spin and electron density op-erators at the site x , and (cid:126)m x and n x are the respectivemean-field parameters.The Helmholtz free energy F ( (cid:126)m, n ) is calculated from H ( (cid:126)m, n ) using the partition function, Z ( (cid:126)m, n ), F ( (cid:126)m, n ) = − β ln Z ( (cid:126)m, n ) = − β ln (cid:16) Tr (cid:16) e − βH ( (cid:126)m,n ) (cid:17)(cid:17) = − β Tr (cid:0) ln (cid:0) e − βh (cid:1)(cid:1) + U (cid:88) x (cid:0) (cid:126)m x − n x (cid:1) , (4)where h xσyσ (cid:48) = tδ σσ (cid:48) + U n x δ σσ (cid:48) − (cid:126)m x · (cid:126)σ σσ (cid:48) ) δ xy , (5)Here, (cid:126)σ σσ (cid:48) is the vector of Pauli matrices.The important result is that in the Hubbard model, theminimum of the Helmholtz free energy, F HF , correspondsto a saddle point of its mean-field counterpart, F ( (cid:126)m, n ), F HF = min (cid:126)m max n F ( (cid:126)m, n ) . (6)See Appendix A for a detailed discussion of the aboverelation.Computing the partial derivatives of F ( (cid:126)m, n ) with re-spect to (cid:126)m and n and setting them equal to zero, we findthe self-consistency relations n x = (cid:104) ˆ n x (cid:105) and (cid:126)m x = (cid:104) ˆ (cid:126)s x (cid:105) .In the case of the Hubbard model, solving the Hartree-Fock equations self-consistently is actually equivalent tofinding a saddle point of the mean-field energy, F ( (cid:126)m, n ).If one fixes the particle density, n on all sites of the lat-tice, then the mean-field calculation is reduced to findinga minimum of the mean-field energy with respect to thespin density, (cid:126)m . Note that these results are for finitetemperatures, but they are still valid in the limit T → n and (cid:126)m to extremize F ( (cid:126)m, n ). To accomplishthis, we use the generalized HF theory, which turns theoriginal minimization problem for the Helmholtz free en-ergy F into a saddle-point problem for the mean-fieldHelmholtz free energy F ( (cid:126)m, n ) (see Appendices). III. THE LIEB LATTICE IN THETIGHT-BINDING LIMIT
The Lieb lattice is a square lattice, with a quarter ofits atoms removed in a regular pattern. Introducing adifferent creation operator in each sublattice, A † , B † , and C † , the tight-binding term of the Hamiltonian of such amodel, H t , is given by [30] t L x (cid:80) x =1 L y (cid:80) y =1 (cid:2) ( A † x,y B x,y + A † x,y C x,y + H.c.)+( A † x,y B x,y − + A † x,y C x − ,y + H.c.) (cid:3) . (7) L x ( L y ) is the number of unit cells along the x ( y ) direc-tion. The hopping terms in the first line are intra-unitcell and the remaining are inter-unit cell. Its eigenvaluesoriginate three energy bands, one of which is flat. Thedispersion relation for periodic boundary conditions is ε ± = ± t (cid:114) cos k x k y , (8)for the two non-flat energy bands, where k α = 2 πn α /L α with n α = 0 , , · · · , L α and α ∈ { x, y } . The flat band is L x × L y -fold degenerate with zero energy. These threeenergy bands are shown in Fig. 1c. The three branchesintersect at the point ( k x , k y ) = ( π, π ). Expanding thedispersion relation in Eq. 8 around this momentum, wefind the Dirac cones ε = t ( k x + k y ). The flat band isbuilt up from B- and C-type orbitals, while the lower andupper bands involve all three lattices A, B, and C. Thislack of uniformity in the distribution of the sublattices inthe energy bands justifies the difference in the occupationnumbers of the sublattices presented in Fig. 1b. IV. INTERACTIONS AND MEAN-FIELD
In this section, we add interactions to the tight-bindingHamiltonian of the Lieb lattice, H t , and reduce the quar-tic dependance of the resulting Hamiltonian on the cre-ation and destruction operators to a quadratic one, usingthe mean-field approximation [12].The key difference between our approach and previousapproaches is that we allow sublattices A, B, and C tohave a different average occupation number each, whilekeeping the total number of particles of the system, N ,fixed [on each point of the ( n, U ) phase diagram]. We be-gin by defining an average particle density on each sub-lattice, n A = n + δ A n B = n + δ B n C = n + δ C , (9)along with the number of particles on each sublattice.For sublattice A, this would be N A = n A L , where L = L x L y is the number of unit cells (which in turn equalsthe number of sites A). To keep the number of particlesequal to N , we apply the restriction N A + N B + N C = N ,which is equivalent to δ A + δ B + δ C = 0 . (10)Setting δ A = δ B = δ C = 0, we would obtain a latticewith its particles evenly distributed, which is what hap-pens in the usual 2D square lattice: all sites have thesame particle density. Aside from total particle numberconservation (Eq. 10) and motivated by the symmetry ofthe lattice, we impose is that δ B = δ C . This gives the im-portant relation δ B = δ C = − δ A / ⇒ n B = (3 n − n A ),which leaves us with one unknown with respect to which F ( (cid:126)m, n ) needs to be maximized. We chose to maximizewith respect to δ A .In our case, we work at T = 0, so that the mean-field free energy of the system can be found by summingthe mean-field energies of the lowest levels that the N particles can occupy (this is the usual Fermi sea). Oneach point ( n, U ), the total energy of the system, E HF ,is obtained by adding the lowest N eigenvalues of themean-field Hamiltonian H HF [2, 23, 24, 31], H HF = (cid:18) H t ( (cid:126)k ) + H δ H m H † m H t ( (cid:126)k + 2 (cid:126)q ) + H δ (cid:19) , (11)and then adding the diagonal terms U L m A + m B + m C − ( n + δ A ) − ( n + δ B ) − ( n + δ C ) ) . (12)The smaller matrices that compose the Hamiltonian H HF are H t ( (cid:126)k ) = t (1 + e ik y ) t (1 + e ik x ) t (1 + e − ik y ) 0 0 t (1 + e − ik x ) 0 0 , (13) H δ = U n + δ A n + δ B
00 0 n + δ C , (14)and H m = − U m A m B e iq y
00 0 m C e iq x . (15)The matrix H HF above is written in the basis { A (cid:126)k , B (cid:126)k , C (cid:126)k , A (cid:126)k +2 (cid:126)q , B (cid:126)k +2 (cid:126)q , C (cid:126)k +2 (cid:126)q } , where the vector (cid:126)q =( q x , q y ) defines the spin orientation in the system, as inthe works by Dzierzawa [24] and Singh [31]. In this paper,we assume that the spin spiral wavenumber (cid:126)q remains thesame as in Ref. 2, even though we allow the system tohave spin and charge modulation. H t ( (cid:126)k ) is the matrixthat corresponds to the tight-binding term of the Hamil-tonian. All other terms have correspondence with theinteraction terms of the Hamiltonian in Eq. 3. Namely, U (cid:80) x (cid:126)m x → UL ( m A + m B + m C ) , U (cid:80) x n x → UL (( n + δ A ) + ( n + δ B ) + ( n + δ C ) ) , U (cid:80) x n x ˆ n x → U diag( n + δ A , n + δ B , n + δ C ) = H δ , − U (cid:80) x (cid:126)m · ˆ (cid:126)s x → − U diag( m A , m B e iq y , m C e iq x ) = H m . (16)The extra imaginary coefficients involving components of (cid:126)q arise from coupling sites on unit cells other than the celllabelled as ( x, y ). From this point forward, we consider t = 1, so that U is given in units of t . V. RESULTS AND DISCUSSION
Our results consist of the values of m A , m B and n A which, for each pair ( n, U ) ∈ [0 , × [0 , E HF , the mean-field energy of theLieb lattice (see Appendix B for a more detailed ex-planation on how to find these saddle points). Fromthese three quantities, we can calculate m C = m B , and n B = n C = (3 n − n A ). We do not impose different occu-pation or magnetization, we simply let the system choosethe values which lead to a saddle point of the mean-fieldenergy. Before performing the calculations for the Lieblattice, we tested the algorithm for the Hubbard modelin a square lattice and found that the occupations of allfour sublattices (A, B, C, and D) were the same, i. e., δ A = δ B = δ C = δ D = 0, while the system chose to havetwo different magnetizations, m A = m D and m B = m C ,reproducing our results in Ref. 23.The results for the Lieb lattice are in Fig. 3. In thefollowing subsections, we discuss each region of interestin more detail. In particular, we study the low U andhigh U regions separately, and finally the near-half-fillingregion ( n ≈ A (a) B (b) A (c) B (d) FIG. 3. Raw numerical results of our mean-field approach, namely for (a) n A , (b) n B , (c) m A and (d) m B . The bold red lineshighlight the U = 0 and U = 20 edges of the plots, and the bold blue line at the center of each plot allows a clearer visualizationof each plot at half-filling ( n = 1 / A. Results near the tight-binding limit ( U → ) As expected, for low U , the average occupations ofthe sublattices, n A and n B = n C (Figs. 3a and 3b, re-spectively), approach those of the tight-binding limit, de-scribed in section III and plotted in Fig. 1b. As for themagnetization amplitudes, sublattice A is paramagnetic( m A = 0) for any n and sublattice B displays a behavioursimilar to that of m in Ref. 2 for low U . In other words, m B is finite between n = 2 / n = 4 / U perturbation, the system becomes aLieb lattice without electron-electron interactions anddisplays the electronic density in Fig. 1b. The energy dispersion relation comprises three bands, as in Fig. 1c,which can be doubly occupied with no additional energycost (because of the absence of U ). Keep in mind that thedispersive bands comprise A-, B-, and C-type orbitals,while the flat band comprises only B- and C-type orbitals.The introduction of the repulsive first-order perturbationmodifies the total energy of the system by shifting the en-ergy bands and adding the diagonal terms of Eq. 12. Thedispersive bands are shifted by + Uδ A and the flat band isshifted by − Uδ A . Additionally, the flat band (which usedto have zero energy) splits into two bands, separated byan amount proportional to U m B .At zero filling, the energy bands are in their original( U = 0) position, because all δ and m are set to zero(having at least one finite m would lead to higher energydue to the term in Eq. 12). As we insert electrons inthe system, they occupy the lowest states in the lowerdispersive band, with n A = 2 n B = 2 n C . Due to theperturbation, this slowly causes the flat band to shift tolower energy and the dispersive bands to shift to higherenergy. Note that before the filling n ≈ /
3, the systemis able to have the lowest energy by displaying paramag-netism ( m A = m B = 0), because, to first order, the onlydependance of the energy on m is in the diagonal term ofEq. 12. This dependance is maintained for any n , in thecase of m A . At n ≈ /
3, the lower dispersive band is al-most full and, due to the flat bands having shifted by thesmall amount − Uδ A , electrons begin to occupy the flatbands, rather than the dispersive bands. From here on,up to n = 1, the magnetization amplitude of sublattice B( m B ) increases, in order to separate the flat bands intotwo, and decrease the energy of the lower flat band, whichis the one being filled. Due to the high density of statesin the flat bands (which are related to B and C atomsonly), newly added electrons choose to occupy sublatticesB and C, until both flat bands are filled, which occurs at n ≈ /
3. Between n = 1 and n = 4 /
3, the magneti-zation m B decreases again, because the lower flat bandis full and electrons are now occupying the higher flatband, which has energy proportional to m B . For filling n ∈ [1 , n ∈ [0 ,
1] region.
B. Results in the strong coupling limit ( U (cid:29) t ) In this subsection, we discuss our results for high U ,which we can assume to be nearly identical to those at U → ∞ , for two reasons. Firstly, through inspectionof the plots in Fig. 3, one readily realizes that the be-haviour at U = 20 is approximately the same as, say, U = 15, and therefore should not change with an in- crease in U . Moreover, if a certain value of U is enoughto impose ferromagnetism (and therefore, for n <
1, theenergy bands are singly occupied with spins in the samedirection), higher values of U will have no additional ef-fect. Secondly, the analytical results for U → ∞ whichwe present in the following paragraphs are in agreementwith the numerical results in Fig. 3 for U = 20.One important first remark is that, as can be seen fromFigs. 3a and 3b, the value of n A or n B for U = 20along the line n ∈ [0 , U = 0 alongthe line n ∈ [0 , U , the inequivalence of sublattices is imposedsolely by the tight-binding terms of the Hamiltonian (seeEq. 11). Therefore, the behaviour of n A and n B is thesame as in the tight-binding limit, albeit with all spinsequally aligned. In fact, at U = 0, the tight-bindingbands become doubly occupied without any additionalenergy cost, while for high U this cost is so high that alltight-binding states (thus, all sublattices) become singlyoccupied before double occupancies are created.To study the U → ∞ limit from a perturbation theorypoint of view, we begin by setting t = 0 in the Hamilto-nian in Eq. 11 and taking the result as the unperturbedHamiltonian. Its eigenvalues are U ( n A ± m A ) , U ( n B ± m B ) , U ( n C ± m C ) . (17)These are six flat bands, with L states each. Positive m and negative m give the same set of eigenvalues, solet us assume positive m , with no loss of generality. At n <
1, electrons occupy the three lowest energy bands: U ( n A − m A ), U ( n B − m B ) and U ( n C − m C ), so that thetotal mean-field energy of the system is given by E U →∞ = U L m A + m B + m C − n A − n B − n C ) + U (cid:34)(cid:88) N A ( n A − m A ) + (cid:88) N B ( n B − m B ) + (cid:88) N C ( n C − m C ) (cid:35) , (18)where we have reintroduced the diagonal terms of Eq. 12.This expression can be simplified using the symmetriesmentioned above: (i) δ A + δ B + δ C = 0, (ii) n B = n C and (iii) m B = m C . Using these three relations and performing the summations up to some fixed (cid:101) N A = L (cid:101) n A ,the total energy for large U and n < E U →∞ = U L m A + 2 m B − n A − n + 3 nn A ) + U L (cid:20)(cid:101) n A ( n A − m A ) + (3 n − (cid:101) n A ) (cid:18)
12 (3 n − n A ) − m B (cid:19)(cid:21) . (19)We find the ground state energy by taking (cid:126) ∇ E U →∞ = (cid:126)
0, where the derivatives are taken with respect to m A , m B and n A . The result is the self-consistency n A = (cid:101) n A andthe relations m A = n A and m B = n B . These two rela-tions hold true for U = 20, as can be realized by compar-ing n A with m A , and n B with m B in Fig. 3. Going backto E U →∞ and replacing n A = m A and n B = m B , we findthat the three bands in Eq. 17 become degenerate withzero energy. Note that setting n A = m A and n B = m B does not lead to a minimum of E U →∞ , but to a saddlepoint, as expected. Again, Hartree-Fock mean-field the-ory is not about finding minima, but rather about find-ing self-consistency, as explained in section II. Finally, weremark that having n A = m A and n B = m B is a conse-quence of a ferromagnetic ground state with the bandsbeing filled with only one spin direction.The next step in the perturbation analysis is to intro-duce the hopping terms of the Hamiltonian of Eq. 11, asthe perturbation. This perturbation lifts the degeneracyof the three zero-energy bands. Up to first order, one ofthem, comprising B- and C-type orbitals, retains zero en-ergy, while each of the other two bands are pushed to pos-itive or negative energy, proportionally to t , and compriseorbitals of the three types, A, B, and C. This new energyband configuration mimics that of the tight-binding limit.This justifies the relative filling of the sublattices at high U (Figs. 3a and 3b), and therefore their magnetization m A = n A and m B = n B . We now have a lower bandinvolving A, B, and C orbitals, an intermediate flat bandbuilt up from B and C orbitals, and a higher energy bandinvolving all three types of atoms. It follows that, as webegin inserting electrons in the system, sublattice A fillsat two times the rate of the B/C sublattices. At n = 1 / − t ) is fulland n A = 1 / n B = 2 × /
8. Between n = 1 / n = 2 /
3, the flat band is filled up using only sublat-tices B and C, therefore n A does not change. Finally, for n ∈ [2 / , t isfilled. For filling n ∈ [1 , n ∈ [0 ,
1] region.
C. Results near half-filling ( n ≈ ) At half-filling, the ground state is a (cid:126)q = ( π, π ) phase.Lieb’s theorem [1] states that in bipartite systems, suchas our Lieb lattice, the ground state at half-filling is fer-rimagnetic [21]. In our case, this ferrimagnetic orderingis ferromagnetic within each sublattice, and antiferro-magnetic between every two nearest-neighbour sites (seeFig. 1a). Although Lieb’s theorem does not provide in-formation on the spin per site, it does state that in oursystem the quantity (2 m B − m A ) should be equal to1 /
2. This served as one of the motivations for this work,as previous studies of the Lieb lattice using mean-field[2] failed to yield the value 1 / U . One of thereasons for this was that the magnetization amplitudewas the same on all sites, A, B, and C, and the relativeoccupation of the three sublattices in the tight-binding Lieb (a) / / / / ( m B - m A ) / ≈ / m B - m A = m B - m A = m A ≈ / ; ∆ m B ≈ ∆ n / m A ≈ / ; ∆ m B ≈ - ∆ n / m A ≈ ; ∆ m B ≈ ∆ n / m A ≈ ; ∆ m B ≈ - ∆ n / ( m B - m A ) / ≈ / m Lieb =(2m B -m A )/2=1/2 n U (b) FIG. 4. (a) Plot of the difference in the absolute values ofthe magnetization of the sublattices, m Lieb = (2 m B − m A ),as a function of n and U , using the data from the plots inFigs. 3c and 3d. The bold red lines highlight the U = 0 and U = 20 edges of the plot, and the bold blue line at the centerallows for easier visualization of the behaviour of (2 m B − m A ) at half-filling. According to a theorem by Lieb [1], thevalue of (2 m B − m A ) at half-filling is 1 /
2, a value which theplot shows to have been achieved by our mean-field approach.(b) Phase diagram illustrating the relative behavior of themagnetization of sublattices A and B, as a function of n and U . limit was used for finite U .One can plot the function (2 m B − m A ), using ourresults for m A and m B , in Figs. 3c and 3d. Such aplot can be found in Fig. 4a. At exactly n = 1 (thebold blue line at the center of the plot), our general-ized Hartree-Fock approach succeeds in yielding the re-sult (2 m B − m A ) = , thus verifying Lieb’s theorem.This leads to the conclusion that our mean-field study,allowing for modulation of m and n in the Lieb lattice,produces more accurate results than imposing the samemagnetization and electronic density in the whole lat-tice. Moreover, this plot indicates that Lieb’s theorem isverified approximately in the wider n ∈ [2 / , /
3] range.The phase diagram given by Fig. 4b was constructedusing the relative behavior of the magnetization of sub-lattices A and B. The dashed lines indicate the bound-aries where the behavior of the magnetization of eachsublattice changes. The magnetization m Lieb is exactly1 / n = 1). For small U ,within the interval 2 / < n < /
3, only sublattice B hasfinite magnetization. For large U , in the n ∈ [2 / , / n A = n B = n C = 1. Plotting n A and n B as a function of U at fixed n = 1, using our results for n A and n B (shown in Figs.3a and 3b, respectively), we checked that we obtainedthe lines n A = 1 and n B = 1. As discussed in SectionV B of this paper, the behaviour of n A and n B at large U consists of two compacted copies of the behaviour atzero U . This holds true for any bipartite system. Conse-quently, for large U , not only are there not charge densitymodulations at half-filling, but also there are no chargemodulations at n = 1 / n = 3 / U can be interpreted as a consequence of the modifica-tions of the mean-field energy dispersion as U increases.More precisely, for small U , one has a single quasi-flatband (twice degenerate with energy ε ∼ U , one has two non-degenerate one-particle quasi-flat bands well separated in energy ( ε ∼ ε ∼ U ).Note that the flat bands for large U are present in theexact solution of the Hubbard model in the subspaceof eigenstates associated with saturated ferromagnetism.As one can conclude from Fig. 4a, the difference in theabsolute values of the magnetizations of the sublatticesgrows only when a flat band is being filled and in fact,we can write m Lieb = 12 (2 m B − m A ) ≈ filling of flat bands , (20) for filling n ≤ , both for large and small U . For n > n ≤ VI. CONCLUSIONS
In summary, we have studied the Lieb lattice usinga mean-field approach and allowing for charge and spindensity modulation. Although theory about the cor-respondence between Hartree-Fock self-consistency andsaddle points of the mean-field energy is relatively old(20 years old), to the extent of our knowledge, this isthe first time it is in fact applied to a system wherecharge modulation is known to occur. We have foundthat, in the limits of low interaction ( U →
0) and veryhigh interaction ( U → ∞ ) results agree with what onewould expect. Namely, the relative occupation of sublat-tices A and B of the bipartite Lieb lattice (where sublat-tice A comprises the atoms with four nearest neighbours,and B denotes the remaining atoms) in the tight-bindinglimit coincides with the results in the literature (for in-stance, in Ref. 20). We have also found that the profile ofthe relative occupation of the sublattices in both strong-coupling and tight-binding are analogous, and one can beinferred if the other is known. The argument is relativelysimple. On the one hand, in the U = 0 limit, the energydispersion is governed by the tight-binding terms of theHamiltonian only and the energy bands in the case ofthe Lieb lattice are as depicted in Fig. 1c. On the otherhand, in the limit U (cid:29) t , the energy bands are separatedby a very high energy gap (of the order of U ), the lowerbands corresponding to no double occupancies and thehigher bands to double occupancies. The gradual fillingof the system is done by singly filling all sites and onlythen jumping to the higher bands and doubly occupyingall sites. Each one of these two filling regimes followsthe relative occupation of sublattices that occurs in thetight-binding limit.At half-filling, our numerical mean-field results ver-ify two important exact results for the Hubbard model.Firstly, one theorem [1] states that in bipartite latticeswith more B-type atoms than A-type atoms, the to-tal spin per unit cell at half-filling is equal to m Lieb =( | B | − | A | ) /
2, where | x | denotes the number of x -typeatoms a unit cell. In the case of our bipartite Lieb lat-tice, we get (2 − / /
2. Our results are plotted inFig. 4. At exactly half-filling we obtained the expectedvalue 1 /
2. Additionally, for large U , in the n ∈ [2 / , / m Lieb ≈ /
2. Interestingly, thisrange includes not only the ferromagnetic region of thephase diagram of Fig. 2a, but also the adjacent spiralregions. Secondly, another theorem [32] states that inthe bipartite Hubbard model, there are no charge den-sity modulations at half-filling, that is, at half-filling allsublattices are equally occupied and half-filled. Our nu-merically calculated relative occupations of sublattices Aand B of a bipartite Lieb lattice, shown in Figs. 3a and3b, are in agreement with this theorem. In addition, wefound that, for large U , not only are there no chargedensity modulations at half-filling, but also there are nocharge modulations at n = 1 / n = 3 / U , the difference in the absolute values of the sub-lattice magnetizations ( m Lieb ) grows or decreases onlywhen a flat band is being filled and furthermore, for n ≤ m Lieb is given approximately by the filling of theflat bands. Note that m Lieb is the unit cell magnetizationin the case of the (cid:126)q = ( π, π ), the ferrimagnetic phase ofthe Lieb lattice found at half-filling.We suggest that much of the above discussion is valid inthe case of other bipartite lattices with flat bands in theenergy dispersion. In the case of non-bipartite latticeswith more than two types of atoms in the unit cell, theanalysis is more complex, because, for instance, an an-tiferromagnet configuration between sublattices may notbe commensurate with the unit cell. Bipartite latticesof various geometries can be realized by manipulatingquantum dot arrays [33] or cold atoms in optical lattices[34].
APPENDIX A: MIN-MAX THEOREM FOR THEHUBBARD MODEL
In this appendix, we discuss how the mean-fieldmethod should be applied to the Hubbard model if parti-cle density is not fixed and in particular, we justify Eq. 6.Our approach follows the method presented in Ref. 26.We start by presenting the usual mean-field approach andthen we explain, using the method presented in Ref. 26,why this approach is somewhat misleading.
A. The usual mean-field method
We begin with the Hubbard Hamiltonian, given by H = t (cid:88) (cid:104) x,y (cid:105) ,σ c † x,σ c y,σ + U (cid:88) x ˆ n x, ↑ ˆ n x, ↓ (21)We then replace the interaction term with the Hartreeand Fock terms (see Eq. 2). These terms are obtainedby considering that each fermionic operator only devi-ates slightly from its mean value, so that in a product ofoperators, we can neglect quadratic terms in these devi-ations.The operators c † x, ↑ c x, ↑ and c † x, ↑ c x, ↓ can be identifiedwith the particle density operator n x and the s + x opera-tor. For consistence with the notation in reference [28],the spin operators in this appendix are c †↑ c ↓ = s + = ( s x + is y ) c †↓ c ↑ = s − = ( s x + is y ) c †↑ c ↑ − c †↓ c ↓ = s z . (22)The vector ˆ (cid:126)s x is the spin density operator on site x . Thepeculiarity of this definition of s + and s − is the factor 1 / s x and s y instead.Replacing this in the interaction term of the Hamiltonianin Eq. 21 gives H HF = t (cid:88) (cid:104) x,y (cid:105) ,σ c † x,σ c y,σ + U (cid:88) x (cid:20) (cid:16) (cid:104) ˆ (cid:126)s x (cid:105) − (cid:104) ˆ n x (cid:105) (cid:17) + 12 (cid:16) (cid:104) ˆ n x (cid:105) ˆ n x − (cid:104) ˆ (cid:126)s x (cid:105) · ˆ (cid:126)s x (cid:17)(cid:21) , (23)We now replace the averages by the mean-field parameters (cid:126)m x and n x . These parameters can be identified with themean values of the magnetization and particle density, respectively, upon extremization of the mean-field free energy,i.e., when the self-consistency equations are satisfied. The Hartree-Fock Hamiltonian becomes H ( (cid:126)m, n ) = (cid:88) x,y,σ,σ (cid:48) h xσyσ (cid:48) c † x,σ c y,σ (cid:48) + U (cid:88) x (cid:0) (cid:126)m x − n x (cid:1) , (24)where h xσyσ (cid:48) = tδ σσ (cid:48) + U n x δ σσ (cid:48) − (cid:126)m x · (cid:126)σ σσ (cid:48) ) δ xy . (25)Here, (cid:126)σ σσ (cid:48) is the vector of Pauli matrices. The function F ( (cid:126)m, n ) is calculated from H ( (cid:126)m, n ) using the partitionfunction, F ( (cid:126)m, n ) = − β ln Z ( (cid:126)m, n ) = − β ln (cid:0) Tr (cid:0) e − βH ( (cid:126)m,n ) (cid:1)(cid:1) = − β Tr (cid:0) ln (cid:0) e − βh (cid:1)(cid:1) + U (cid:80) x (cid:0) (cid:126)m x − n x (cid:1) . (26)At this stage, one usually finds the minimum free en- ergy, F HF , by minimizing F ( (cid:126)m, n ) with respect to the0mean-field parameters and this would lead to the usualself-consistency equations [35]. Clearly, this works forfixed particle density, but in the previous expressionwe allow for variable particle density and the respectivequadratic term has a negative coefficient. If one imposeda minimization with respect to the parameter n x , con-vergence would not be achieved in a numerical approach,unless one limits the possible values of n x to a certaininterval, in which case the result of the numerical mini-mization would lie at the boundary of this interval. Thisreflects the fact that one should not minimize with re-spect to the parameter n x , but instead maximize, as weexplain in the next subsection. B. A different perspective for the mean-fieldmethod
In this subsection, we present the mean-field approachwhich should be applied when one takes into account thepossibility of non-uniform particle density in a lattice.For a rigorous proof see, for instance, Refs. 26–28.Our goal is to know the thermal equilibrium state ofthe system which minimizes the free energy. These statesare defined in terms of density matrices. Having an ex-act free energy would require having an exact partitionfunction, which in turn would require an exact diagonal-ization of the Hubbard model. The Hartree-Fock methodprovides an approximation to the exact equilibrium state,in terms of many-body states of non-interacting parti-cles, replacing the quartic terms of the Hamiltonian byone-particle potentials (which are adjusted to provide thebest possible approximation). The free energy obtainedin the Hartree-Fock method provides an upper boundto the exact free energy (a consequence of the varia-tional theorem). Since the particles are independent, theHartree-Fock state can be written as a one-particle den-sity matrix, γ ij = (cid:104) c † i c j (cid:105) [26]. The objective of the mean-field method is indeed to find the minimum free energyin the set of free energies associated with the possiblestates of N independent particles (or equivalently, asso-ciated with the possible one-particle density matrices), F HF = min γ F ( γ ) = min γ [ E ( γ ) − β S ( γ )] . (27)where β is the einverse temperature. In the previousexpression, no mean field parameters are present andthe free energy is determined for each γ using the exact Hamiltonian (the mean field approximation is associatedwith the state, not with the Hamiltonian). However, cal-culating min γ F ( γ ) going through all γ is not pratical,so one introduces mean-field parameters. In the nextparagraphs, we introduce these parameters, following aderivation in Ref. 26. For now, let us assume zero tem-perature.In the case of the Hubbard Hamiltonian, the Hartree-Fock energy functional, E ( γ ), can be written as (Ref. 26Lemma 2) E ( γ ) = Tr[ T γ ] + U (cid:88) x [ (cid:104) ˆ n x (cid:105) − (cid:104) ˆ (cid:126)s x (cid:105) ] , (28)where x labels the lattice sites and T is the matrix whoseelements t xy are the transition amplitudes of an electronto move from site x to site y or vice versa. The averagesof the electronic and spin densities at site x are (cid:104) ˆ n x (cid:105) = Tr (cid:104) ˆ I · γ (cid:105) , (cid:104) ˆ (cid:126)s x (cid:105) = Tr (cid:104) (ˆ I ⊗ ˆ (cid:126)σ ) · γ (cid:105) , (29)respectively.The first step to introduce the variational parametersis to use the simple fact that x ≥ xy − y for all x, y ∈ R n , (30)where equality holds if and only if x = y . This implies x = max y (2 xy − y ) , (31)and therefore, we can write (cid:104) ˆ n x (cid:105) = max n x (cid:8) (cid:104) ˆ n x (cid:105) n x − n x (cid:9) (32)and − (cid:104) ˆ (cid:126)s x (cid:105) = min (cid:126)m x (cid:110) (cid:126)m x − (cid:104) ˆ (cid:126)s x (cid:105) · (cid:126)m x (cid:111) . (33)where n x and (cid:126)m x are, at this point, an arbitrary constantand vector. Inserting this into Eq. 28, the energy for agiven state γ at zero temperature assumes the form E ( γ ) = min (cid:126)m max n { E ( n, (cid:126)m, γ ) } , (34)where E ( n, (cid:126)m, γ ) is given by E ( n, (cid:126)m, γ ) = Tr (cid:34)(cid:32) T − U (cid:88) all sites (cid:104) n x · ˆ I − (cid:126)m x · ˆ (cid:126)σ (cid:105)(cid:33) γ (cid:35) + U (cid:88) all sites (cid:0) (cid:126)m x − n x (cid:1) , (35)The aim of the mean-field method is to obtainmin γ E ( γ ) = min γ min (cid:126)m max n { E ( n, (cid:126)m, γ ) } . (36) The crucial point of the mean-field method is the possi-1bility of exchanging the order of the extremization, whichis implicitly done in standard mean-field. In fact, if par-ticle density is fixed, one hasmin γ E ( γ ) = min γ min (cid:126)m { E ( (cid:126)m, γ ) } = min (cid:126)m min γ { E ( (cid:126)m, γ ) } , (37)since the order of minimization is irrelevant, and one re-covers the usual mean-field picture, in which one has tominimize the energy with respect to the mean-field pa-rameter (cid:126)m .On the other hand, if the particle density n x is allowedto actually depend on x , in Ref. 26 it was shown thatindeedmin γ min (cid:126)m max n { E ( n, (cid:126)m, γ ) } = min (cid:126)m max n min γ { E ( n, (cid:126)m, γ ) } . (38)Here, we present a simple physical interpretation for thisresult.First, one should note that Eq. 35 for given (cid:126)m and n corresponds to the energy of independent particles sub-jected to one-particle potentials which depend on thegiven (cid:126)m and n . Looking at the left-hand side of Eq. 38,after one has minimized and maximized E ( n, (cid:126)m, γ ) withrespect to (cid:126)m and n respectively, the remaining minimiza-tion (with respect to γ ) will lead, at zero temperature, toa filled Fermi sea associated to these one-particle poten-tials. These one-particle potentials might, for instance,generate imbalance between the number of up and downspins and lead to finite magnetization.Second, Eqs. 32 and 33, associated with minimizationand maximization with respect to (cid:126)m and n , imply theself-consistency equations (cid:104) ˆ n x (cid:105) = n x , (cid:104) ˆ (cid:126)s x (cid:105) = (cid:126)m x , (39)and therefore, computing the left-hand side of Eq. 38can be interpreted as the following instruction: ”Amongthe set of one-particle density matrices that satisfy theself-consistency equations, find the one with minimumenergy”. As we said above, this minimum energy cor-responds to a filled Fermi sea at zero temperature, soin the set of states corresponding to filled Fermi seas,there is one that satisfies the self-consistency equations.This helps us understand the right-hand side of Eq. 38in the following way: ”Among the set of all filled Fermiseas, find (the energy of) the state which satisfies theself-consistency equations”.This argument can be generalized for finite tempera-ture [26], taking into account the entropy contribution.In this case, instead of filled Fermi seas, one has the one-particle density matrices that minimize the free energyin the case of independent particles. If instead of fixednumber of particles, one imposes a fixed chemical poten-tial (in which case the grand-canonical potencial wouldreplace the free energy), this one-particle density matrixwould become the Fermi-Dirac distribution function. APPENDIX B: HOW TO EXTREMIZE E HF In this appendix, we explain in more detail the algo-rithm we used for finding saddle points.Our first approach to extremize E HF was perhaps themost intuitive non-brute force one. It consisted in max-imizing E HF with respect to δ A , δ B , and δ C using ourresults for q x , q y , m A and m B of the Hubbard model ina square lattice [23]. The lattice size was also kept at100 × E HF . We used δ A . We note that the minimum of E HF with respect to δ A was found to be negative infinity (us-ing the results in Ref. 23 for m A and m B as startingpoints). Fixing q x , q y , m A and m B means that the mag-netic phases are kept the same as our previous ones, onlythe occupation numbers in the sublattices change. How-ever, finding saddle points by starting from a minimumwith respect to one direction (the m direction) and thenfinding the maximum with respect to an orthogonal di-rection (the n direction) turned out to diverge on mostpoints of the phase diagram. In fact, the suggestion thatsaddle points of a function can be found by starting ata random point in the function, minimizing the functionwith respect to the minimizer variables, and then usingthose new points to maximize the function with respectto the maximizer variables, is false. As a matter of fact,this statement remains false even in a more generalizedcase. One might think that by successively minimizingand maximizing the function, one would eventually reacha saddle point. This is also not necessarily true. Whentrying to use this method to extremize E HF , one findsthat, in most points of the diagram, the value of E HF re-sulting from the last maximization and the value of E HF resulting from the last minimization differ by orders ofmagnitude comparable to those of E HF itself. Addition-ally, the values of m A , m B , δ A and δ D which extremize E HF (or so one thought) do not converge on each suc-cessive iteration, they alternate between several possibleresults.In order to better understand why this method mayfail to find saddle points, one can consider a much simplerexample, the function f ( x, y ) = sin( x ) cos( y ). This func-tion is periodic, oscillates between -1 and 1, and has in-finitely many minima, maxima and saddle points (Fig. 5).Let us now assume that we want to find the saddle pointsof f ( x, y ) using the method of minimizing and maximiz-ing several times alternating between the two, using x asminimizer and y as maximizer. If we start with a non-saddle point, the algorithm fails to find a saddle point andinstead, after 2 ou 3 steps, alternates between absoluteminima and maxima. Even if we start at a saddle point,successive iterations will move away from it and back tojumping between maxima and minima (see Fig. 5a). Noinitial point will allow this algorithm to converge to asaddle point, no matter how close it is to it. This is pre-cisely what happened with the case of E HF . Depending2 x y −5 0 5−6−4−20246 (a) x y −5 0 5−6−4−20246 (b) x y −2 0 2 4 6−2−1012345 (c) FIG. 5. (a,b) Contour plots of the function function f ( x, y ) = sin( x ) cos( y ) in the range x, y ∈ [ − π, π ]. The maxima are atthe center of the red circles, the minima are at the center of the blue circles, and the saddle points are at the intersections ofthe green lines. The straight bold lines overlapped with the contour plots are the path taken by an algorithm which attemptsto find a saddle point of f ( x, y ), by starting at point ( − . , π − .
5) (red asterisk), and (a) successively minimizing in the x direction and maximizing in the y direction, and (b) successively minimizing in the direction which makes an angle of π/ x axis and maximizing in a direction which makes an angle of π/
10 with the y axis. (c) Contour plot of the function f ( x, y ) = 2 x + 6 xy − y − x in the range ( x, y ) ∈ [ − , × [ − , f ( x, y ) withrespect to the x direction and the large (red) dot represents the maximum of these minima, which is a good approximation ofa saddle point. on the pair ( n, U ) in question, E HF may or not behavein a way that allows saddle points to be found using thismethod.One possible alternative approach to finding saddlepoints of f ( x, y ) is to attempt to extremize the func-tion in a different direction (or rotating the axes, whichis another way to see it). Instead of using x and y asvariables, we can use two auxiliary variables which arelinear combinations of x and y but still orthogonal. Forinstance, we can alternate between the following two: • replace x → x + r cos θ , and y → y + r sin θ , andminimize with respect to r ; • replace x → x + r sin θ , and y → y − r cos θ , andmaximize with respect to r .Setting θ = 0 would reduce this to simply minimizingwith respect to x and maximizing with respect to y , asdescribed above. Due to the symmetry of the simplefunction we used, we should ideally use θ = π/ π/ π/
4, saddlepoints may not be found anymore. Using a ”wrong” anglewill make the algorithm diverge even if we start veryclose to a saddle point, even if the starting point is asaddle point itself. The algorithm will often be circlingthe saddle point (see Fig. 5b).The conclusion is that alternating between maximiza-tion and minimization is too sensitive to the initial pointto be used to find saddle points. In the case of the sim-ple highly-symmetric function f ( x, y ), the saddle pointis located at the center of the squares or rectangles thatare drawn when connecting the points given by the algo-rithm, but this was shown to not be the case with E HF . Moreover, this ideal angle depends on the pair ( n, U )which is being studied, and some angles may not evenproduce closed polygons, but rather diverge to infinity ina certain direction.The working alternative is failproof in the sense thatit will always converge to one saddle point. It consists ofcalculating the value of f ( x, y ) for many points in an areawhich is known to contain at least one saddle point, mak-ing a list of the maximum value of the function for each x and then finding the minimum of all the maxima thatwere found. If f ( x, y ) is continuous we will end up on asaddle point. In the case of E HF , we have not two vari-ables but three: m A , m B and δ A . In order to achieve anacceptable precision, it would be reasonable to calculate,say, 100 values of E HF in each direction ( m A , m B , and δ A ), for a total of 10 values for each pair ( n, U ) in ourphase diagram, which could take a long time and wouldhave a precision of about two decimal places. A moreefficient alternative is to divide each direction into fewerparts, say 7 or 8, finding an initial approximation to thesaddle point, and then repeating this a dozen times con-sidering a smaller hypercube centered on the new point.Assuming each direction is divided into 7 parts (i.e. wecalculate 8 function values in each direction) and the sideof the hypercube is halved with each of 12 iterations, wenow calculate a total of 7 × ≈ E HF perpair ( n, U ). This was the procedure we used, which has aprecision of around 1 / f ( x, y ) = 2 x + 6 xy − y − x and successfully found the saddle point (3 , ACKNOWLEDGEMENTS
R. G. Dias acknowledges the financial support from thePortuguese Science and Technology Foundation (FCT) through the program PEst-C/CTM/LA0025/2013. J.D. Gouveia acknowledges the financial support from thePortuguese Science and Technology Foundation (FCT)through the grant SFRH/BD/73057/2010. [1] Elliott H. Lieb. Two Theorems on the Hubbard model.
Physical Review Letters , 62:1201–1204, 1989.[2] J. D. Gouveia and R. G. Dias. Magnetic phase diagramof the Hubbard model in the Lieb lattice.
Journal ofMagnetism and Magnetic Materials , 382:312–317, 2015.[3] J. G. Bednorz and K. A. M¨uller. Possible high T c super-conductivity in the Ba-La-Cu-O system. Zeitschrift furPhysik B Condensed Matter , 64:189–193, June 1986.[4] M. A. Kastner, R. J. Birgeneau, G. Shirane, and Y. En-doh. Magnetic, transport, and optical properties ofmonolayer copper oxides.
Rev. Mod. Phys. , 70:897–928,Jul 1998.[5] R. G. Dias and J. M. B. Lopes dos Santos. Simplerepresentation of the eigenstates of the U → ∞ one-dimensional Hubbard model. Journal de Physique I ,2:1889–1897, October 1992.[6] Michael P. Marder.
Condensed Matter Physics . JohnWiley and Sons, 2000.[7] David R. Penn. Stability theory of the magnetic phasesfor a simple model of the transition metals.
Phys. Rev. ,142:350–365, Feb 1966.[8] J. Dorantes-Davila, J. L. Moran-Lopez, and M Avignon.Ground-state solutions of the Hubbard model.
Phys. Rev.B , 27:575–577, 1983.[9] E. Kaxiras and E. Manousakis. Ground state of thestrong-coupling Hubbard Hamiltonian: A numerical di-agonalization study.
Phys. Rev. B , 37:656–659, 1988.[10] S. N. Coppersmith and Clare C. Yu. Phase diagram of theHubbard model: A variational wave-function approach.
Phys. Rev. B , 39:11 464 – 11 474, 1989.[11] A. Richter, G. R¨opke, and F. Goedsche. Functional In-tegral Approach for the Hubbard Model with ArbitraryElectron Density.
Physica Status Solidi B Basic Research ,88:189–198, July 1978.[12] Sanjoy Sarker, C. Jayaprakash, H. R. Krishnamurthy,and Wolfgang Wenzel. Spiral states in the square-latticeHubbard model.
Physical Review B , 43:8775–8778, Apr1991.[13] Edwin Langmann and Mats Wallin. Mean Field Mag-netic Phase Diagrams for the Two Dimensional t − t (cid:48) − U Hubbard Model.
Journal of Statistical Physics , 127:825–840, 2007.[14] P. A. Igoshev, M. A. Timirgazin, A. A. Katanin, A. K.Arzhnikov, and V. Yu. Irkhin. Incommensurate mag-netic order and phase separation in the two-dimensionalHubbard model with nearest- and next-nearest-neighborhopping.
Phys. Rev. B , 81:094407, 2010.[15] W. Schumacher. On incommensurate phases in the mag-netic phase diagram of the hubbard model.
Physica Sta-tus Solidi (b) , 119:235–238, 1983.[16] A. Mielke. Exact ground states for the Hubbard modelon the Kagome lattice.
Journal of Physics A , 25:4335,1992.[17] Hal Tasaki. Ferromagnetism in the Hubbard Models withDegenerate Single-Electron Ground States.
Physical Re-view Letters , 69:1608–1612, 1992. [18] V. J. Emery. Theory of high- T c superconductivity inoxides. Physical Review Letters , 58:2794–4797, 1987.[19] R. T. Scalettar, D. J. Scalapino, R. L. Sugar, and S. R.White. Antiferromagnetic, charge-transfer and pairingcorrelations in the three-band Hubbard model.
PhysicalReview B , 44:770–781, 1991.[20] Hu Wang, Shun-Li Yu, and Jian-Xin Li. Spin fluctuationsand unconventional pairing on the Lieb lattice.
PhysicsLetters A , 378:3360–3365, 2014.[21] Andreas Mielke and Hal Tasaki. Ferromagnetism inthe Hubbard Model.
Communications in MathematicalPhysics , 158:341–371, 1993.[22] K. Noda, K. Inaba, and M. Yamashita. Flat-band ferro-magnetism in the multilayer Lieb optical lattice.
Phys.Rev. A , 90:043624, 2014.[23] J. D. Gouveia and R. G. Dias. Spiral ferrimagnetic phasesin the two-dimensional Hubbard model.
Solid State Com-munications , 185:21–24, 2014.[24] M. Dzierzawa. Hartree-Fock theory of spiral magneticorder in the 2-d Hubbard model.
Z. Phys. B , 86:49–52,1992.[25] E. Langmann and M. Wallin. Restricted path integralapproach to the doped Hubbard model.
Europhysics Let-ters , 37 (3):219–224, 1997.[26] V. Bach and J. Poelchau. Hartree-Fock Gibbs states forthe Hubbard model.
Markov Processes and Rel Fields ,2(1):225–240, 1996.[27] V. Bach, E. H. Lieb, and J. P. Solovej. GeneralizedHartree-Fock theory and the Hubbard model.
Journalof Statistical Physics , 76:3–89, 1994.[28] Jonas de Woul. A restricted Hartree-Fock study of the2D Hubbard model. Master’s thesis, Royal Institute ofTechnology, 2007.[29] Edwin Langmann and Mats Wallin. Mean-field ap-proach to antiferromagnetic domains in the doped Hub-bard model.
Physical Review B , 55:9439–9451, 1997.[30] M. Nita, B. Ostahie, and A. Aldea. Spectral and trans-port properties of the two-dimensional Lieb lattice.
Phys-ical Review B , 87:125428, 2013.[31] A. Singh, Z. Tesanovic, and H. H. Kim. Instability of thespiral state of the doped Hubbard model.
Pramana - J.Phys. , 38:211–217, 1992.[32] Elliott H. Lieb, Michael Loss, and Robert J. McCann.Uniform density theorem for the Hubbard model.
Journalof Mathematica Physics , 34:891–898, 1993.[33] Hiroyuki Tamura, Kenji Shiraishi, and HideakiTakayanagi. Ferromagnetism in Semiconductor DotArray.
Jpn. J. Appl. Phys. , 39:L241, 2000.[34] N. Goldman, D. F. Urban, and D. Bercioux. Topologi-cal phases for fermionic cold atoms on the Lieb lattice.
Physical Review A , 83:063601, 2011.[35] H. Bruus and K. Flensberg.