Charge Order of Strongly Bounded Electron Pairs on the Triangular Lattice: the Zero-Bandwidth Limit of the Extended Hubbard Model with Strong Onsite Attraction

The inhomogeneous distribution of charged particles is known as charge-order phenomena. Recent experimental studies of (quasi-) two-dimensional materials found different ordered patterns on the triangular lattice. In this work, the preliminary studies of the extended Hubbard model (with the intersite Coulomb density-density repulsion) in the atomic limit on the triangular lattice are presented. The model is investigated within the mean-field approach that treats the onsite term exactly. The finite temperature phase diagrams including metastable phases and some characteristics of the model are determined for the limit of strong onsite attraction. It is shown that in the limit considered all transitions between stable phases are discontinuous (for fixed chemical potential), which is totally different from the model on the alternate lattices. Continuous transitions between metastable phases are also found.


Introduction
The charge-order phenomenon is associated with inhomogeneous distribution of charged particles on a lattice [1,2]. The particles can form different patterns depending, among others, on relative values of interactions in the system as well as on a type of the lattice on which they interact. In this paper, the paradigmatic model for strongly correlated fermion (electron) system is studied on the triangular lattice. The motivation of this work is recent experimental [3][4][5][6][7][8][9][10][11] and theoretical [12][13][14][15][16][17][18][19] studies of quasi-two-dimensional systems, e.g., NbSe 2 or organic conductors, where various charge-ordered patterns have been observed on such a type of the lattice. The extended Hubbard model in the zero-bandwidth (or equivalently, atomic) limit analyzed in this work (e.g., Refs. [20][21][22][23][24][25][26][27][28][29][30][31][32][33]) has the following form wheren i =n i↑ +n i↓ is a operator of particle number at site i,n iσ =ĉ + iσĉ iσ is a operator of number of particles with spin σ (σ ∈ {↑, ↓}) at site i, andĉ i ĉ + i is a operator of annihilation (creation) of a fermion (electron) with spin σ at site i. i, j denotes the summation over the nearestneighbor sites, independently. U denotes the Hubbard onsite interaction, which is assumed to be attractive (i.e., U < 0), W is the intersite density-density Coulomb repulsion (i.e., W > 0), whereas μ is the chemical potential. z denotes the coordination number (the number of nearest neighbors).
The model is analyzed using the variational approach that treats the onsite term exactly and the intersite term within the standard mean-field approximation (including only the Hartree terms) and in the thermodynamic limit (in the grand canonical ensemble) [26][27][28][29][30]. This approach is an exact one in the limit of large coordination number (z → +∞) [26,[34][35][36], but even that it can be a good starting point for the triangular lattice with z = 6 ( Fig. 1). For the ground state of the model on this lattice, it is an exact theory. For finite temperatures, it obviously overestimates the critical temperatures, at which the long-range order vanishes.
The model is studied on the triangular lattice. This lattice can be divided into three equivalent sublattices labeled by α = A, B, C (cf. Fig. 1) (the tri-sublattice assumption). Various homogeneous phases occurring on phase diagrams Fig. 1 The diagram of the triangular lattice. Three equivalent sublattices (α = A, B, C) are denoted by different symbols. Symbol shapes schematically show a type of the charge order considered in the present paper (e.g., n A = n B = n C ) of the model are characterized by average occupations of the sublattices: n α = 3 L i∈α n i , where L is a number of lattice sites (within the approach used n α = n i for each i ∈ α).

Numerical Results for Large Onsite Attraction
For the limit studied (i.e., U/W = −10), only two different phases are found: (i) a nonordered (NO) phase with n A = n B = n C and (ii) a charge-ordered (CO) phase, where the concentrations in two of three sublattices are equal to each other and the concentration in the third sublattice is different (in a general case, one has three possibilities: n A = n B = n C , n B = n A = n C , or n C = n A = n B ). The concentration which occurs only in one of the sublattices is denoted as n 1 (e.g., n 1 = n A ), whereas the other one (occurring in two of the sublattices) is indicated as n 2 (e.g., n 2 = n B = n C ). One can imagine the existence of the ordered phase where all three n α 's are different (i.e., n A = n B n B = n C , and n C = n A ; cf., e.g., Ref. [33]), but such a phase has not been found in the studied range of model parameters. To distinguish the CO and NO phases, charge-order parameter Δ defined as Δ = (1/2)(n 1 − n 2 ) is introduced. In the CO phase, Δ = 0, whereas in the NO phase, Δ = 0. Notice that such defined Δ is different from order parameter Δ Q for checker-board-ordered phases occurring on alternate lattices, which is defined as the difference between the concentrations in two interpenetrating sublattices [27][28][29][30]. The expresions for Δ and n α for the triangular lattice can be derived from Ref. [28], where the equations for a general case were derived within a site-dependent approach. denote the non-ordered and two different charge-ordered phases. Solid lines denote discontinuous CO1-CO2 (atμ = 0) and order-disorder (CO1-NO and CO2-NO) transitions. Dashed line denotes continuous CO1-NO and CO2-NO transitions between metastable phases. The label in bold (italics in brackets) denotes that in the region such a phase is stable (metastable, respectively). The regions, where the NO phase is stable and the CO1 or CO2 phase is metastable, are very narrow (located above the order-disorder line) and not visible (but schematically labeled) in the figure

The Phase Diagrams
The phase diagram of the model for the fixed chemical potential (as a function ofμ = μ − U/2 − W ) and for U/W = −10 is shown in Fig. 2. First, the occurrence of stable phases, i.e., phases with the lowest grand canonical potential, is discussed. On the diagram, two distinguishable The phase diagram for W > 0 and U/W = −10 as a function of particle concentration n. Solid lines restrict the ranges of PS1:NO/CO1, PS2:CO2/NO, and PS3:CO1/CO2 states occurrence. In the PS regions, one or two homogeneous phases can be metastable. In PS1 and PS2 regions, the metastable phases are not labeled (details in the text). Dotted line denotes discontinuous CO1-CO2 transition between metastable phases. Other denotations as in Fig. 2 CO phases are present: (i) a CO1 phase with Δ > 0 stable forμ < 0 and (ii) a CO2 phase with Δ < 0 occurring for forμ > 0. These two regions of the CO phase are separated by the first-order CO1-CO2 transition Fig. 4 Dependencies of order parameters: a particle concentrations in sublattices: n 1 = n A and n 2 = n B = n C , b total particle concentration n, and c charge-order parameter Δ = 1 2 (n 1 − n 2 ) as a function of chemical potentialμ for k B T /W = 0.2, W > 0, and U/W = −10. Solid and dash-dotted lines correspond to stable solutions, whereas dotted and dashed lines correspond to metastable solutions. The metastable solutions for the CO1 and CO2 phases inside the NO region are not shown (cf. also Fig. 5). Labels as in Fig. 2 atμ = 0. The order-disorder transitions, i.e., NO-CO1 and CO2-NO transitions, are also discontinuous ones. The maximal temperature of the transitions is k B T M /W ≈ 0.5 and it is located atμ = 0. For |μ|/W > 1, a re-entrant behavior is found. Each of these discontinuous transitions is associated with an occurrence of the metastable phases in their neighborhood. Thus, as a consequence of the CO1-CO2 transition, the CO1 (CO2) phase is metastable in the CO2 (CO1, respectively) phase region of stability (cf. also Fig. 4). Similarly, next to the NO-CO1 (CO2-NO) line in the region of the NO phase stability, the CO1 (CO2, respectively) phase is metastable for some define range of model parameters (very narrow region, not visible in the figure), whereas the NO phase is metastable in the CO1 (CO2, respectively) region (cf. also Figs. 5 and 6). It turns out that inside the CO1 (CO2) region the continuous NO-CO2 (CO1-NO, respectively) transition between metastable phases occurs (dashed line in Fig. 2).   Figure 3 shows the phase diagram of the system considered for U/W = −10 as a function of total electron concentration n = (1/3)(n 1 + 2n 2 ). Due to the fact that the CO1-CO2, NO-CO1, and CO2-NO transitions are first order for fixed chemical potential, three regions of phase separations, PS1:NO/CO1, PS2:CO2/NO, and PS3:CO1/CO2, exist on the diagram. In these regions, phase-separated states have the lowest free energy and two homogeneous phases can coexist in define range of n. In particular, in the PS3 region, two charge-ordered phases, i.e., the CO1 and CO2 phases, coexist. The fraction of the system occupied by each phase depends on n (for details, see, e.g., Refs. [37][38][39]).
In the regions of the phase separations, one or two homogeneous phases can be metastable. In the whole region of the PS3 state occurrence, both CO1 and CO2 phases are metastable (they have higher free energy than the PS3 state, which free energy f PS3 can be calculated by Maxwell's construction [37][38][39]). But for n < 1, the free energy of the CO1 phase is lower than that of the CO2 phase (f PS3 < f CO1 < f CO2 ), whereas for n > 1, the CO1 phase has higher free energy than the CO2 phase (f PS3 < f CO2 < f CO1 ). At n = 1, a first-order transition between the CO1 and CO2 (metastable) phases occurs for any T < T M .
The situation is more complex in a case of regions of the PS1 and PS2 states occurrence. For fixed temperature k B T /W , there are two possible sequences of phase transitions with increasing n (for n < 1). The first one occurring at low temperatures (for T < T C , where k B T C /W ≈ 0.33) is NO-PS1(NO)-PS1(NO;CO1)-PS1(CO1;NO)-PS1(CO1;CO2)-CO1(CO2). The denotation, e.g., PS1(NO;CO1), denotes that in the PS1 region, two phases, NO and CO1, are metastable with f PS1 < f NO < f CO1 , whereas, e.g., PS1(NO) indicates that the NO phase is the only metastable solution with f NO > f PS1 . At the begging, the NO phase is only metastable phase in the PS1 region. Next, inside the PS1 region, the metastable CO1 solution appears [at n m (T )]. With increasing n, a first-order NO-CO1 transition between metastable phases occurs at n t (T ) (before it: f PS1 < f NO < f CO1 , whereas after it: f PS1 < f CO1 < f NO , cf. Fig. 5), analogously to the CO1-CO2 transition inside the PS3 region discussed previously. After that, but still in the PS1 region, the second-order NO-CO2 transition between metastable phases occurs at n c (T ) and f PS1 < f CO1 < f CO2 for n > n c (T ). The second possible sequence occurring for T C < T < T M is NO-PS1(NO)-PS1(NO;CO1)-PS1(CO1;NO)-CO1(NO)-CO1(CO2). In such a case, the metastable CO2 phase does not occur in the PS1 region and the secondorder NO-CO2 transition between metastable phases occurs in the CO1 region, not in the PS1 region as previously (cf. Fig. 6). At T = T C , the NO-PS1(NO)-PS1(NO;CO1)-PS1(CO1;NO)-CO1(CO2) sequence is present. Note that in Fig. 3, two lines inside the region of the PS1 state occurrence are not denoted: a line, where the CO1 solution appears as the second metastable phase [PS1(NO)-PS1(NO;CO1) boundary, n m (T )], and a line of the first order NO-CO1 transition between metastable phases [PS1(NO;CO1)-PS1(CO1;NO) boundary, n t (T )].
Due to the particle-hole symmetry of the model discussion for n > 1 is similar, one should change PS1 to PS2, exchange CO1 and CO2 with each other, and discuss the sequences with decreasing n (or equivalently change n into 2−n andμ into −μ). As it is shown before, near n = 1, a CO1(CO2)-PS3(CO1;CO2)-PS3(CO2;CO1)-CO2(CO1) sequence occurs for any T < T M .

Dependencies of the Order Parameters
To give deeper insight into properties of the system, Figs. 4, 5, and 6 present dependencies of order parameters as a function of chemical potential for fixed k B T /W . In Fig. 4a, concentrations n 1 and n 2 are presented, whereas in Fig. 4b and c, dependencies of total electron concentration n and charge-order parameter Δ are shown, respectively. It is clearly seen that atμ = 0, the first-order CO1-CO2 transition occurs, namely Δ and n change discontinuously from Δ − > 0 and n − < 1 (forμ = 0 − ) to Δ + = −Δ − < 0 and n + = 2 − n − > 1 (forμ = 0 + ). The metastable solutions exist on both sides of the transition. Notice that for fixed n between n − and n + , the PS3:CO1/CO2 state has the lowest free energy, but in this region, both homogeneous phases are metastable. The comparison of free energies of both metastable phases reveals that indeed for n = 1, there is the discontinuous CO1-CO2 transition between metastable phases. Now, behavior of the system in the neighborhood of the discontinuous order-disorder transition is considered. Only a case ofμ < 0 (n < 1) is discussed here, but discussion of a case ofμ > 0 (n > 1) is analogous due to the particle-hole symmetry (as indicated at the end of the previous section). At the NO-CO1 transition for fixedμ concentration, n changes discontinuously from n NO (in the NO phase) to n CO1 (in the CO1 phase), and thus for n NO < n < n CO1 , the PS1:NO/CO1 state occurs for fixed n. In Figs. 5 and 6, one can notice that the CO1 phase is the metastable on the left side from a point of the first-order NO-CO1 transition. One can also see that in the range between n NO and n CO1 , only the NO phase is metastable (for n < n m ) or two (CO1 and NO for n m < n < Min{n c , n CO1 }, or, CO1 and CO2 only for T < T C and n c < n < n CO1 ) homogeneous phases are metastable. The location n t of the first-order NO-CO1 transition between metastable phases for fixed n can be determined by the comparison of their free energies (not indicated in the figures, and n t > n m ). Moreover, in all figures, the continuous NO-CO2 transition between metastable phases occurring at n c is presented. For fixed n and for T < T C , the NO-CO2 transition between metastable phases occurs inside the PS1 region, i.e., n NO < n c < n CO1 (cf. Fig. 5a). Moreover, n c is also located in the region, where the CO1 phase is the metastable phase with lower free energy (i.e., n m < n t < n c and f PS1 < f CO1 < f NO = f CO2 at n c , precisely). In contrary, for T C < T < T M , the NO-CO2 transition between metastable phases occurs inside the CO1 region, i.e., n c > n CO1 (cf. Fig. 6a).
At the end of this section, one should stress that for T C < T < T M , the following relations between all critical concentrations are satisfied: n NO < n m < n t < n CO1 < n c , whereas for T < T C , one gets n NO < n m < n t < n c < n CO1 . To be strict, at low temperatures (smaller than k B T /W ≈ 0.1), n NO < n c < n m < n CO1 could be also possible, which would correspond to NO-PS(NO)-PS(CO2)-PS(CO1;CO2)-CO1(CO2) sequence of transitions (f CO1 < f CO2 for any T and n < 1). This is the only possible sequence not discussed previously. However, its existence cannot be unambiguously adjudged numerically due to very small difference |n c − n m |, which is of the order of a few thousandths. Nevertheless, for T → 0, all above critical concentrations goes to zero.
Note that generally homogeneous phases can be metastable on both sides of a PS region (such as for the PS3 region above or the PS region in the model of Ref. [40]) or any homogeneous phase cannot be metastable in some part of a PS region (cf. the model with longer-range interactions on alternate lattice [29]). In the present work, none of these cases has been found for the PS1 and PS2 regions.

Summary and Final Remarks
In this work, the extended Hubbard model in the atomic limit on the triangular lattice for strong onsite attraction was discussed in the context of charge-order phenomena and the analysis including metastable phases has been performed. The phase diagram of the model including metastable phases was derived for U/W = −10. All transitions between stable phases for fixed chemical potential were found to be discontinuous and thus phase separations occur for fixed concentration. It is different from a case of bipartite (alternate) lattices (with only nearest-neighbor interaction), where only a continuous order-disorder transition occurs for U < 0 [27][28][29][30].
Notice that the discussed model in the limit of U → −∞ is equivalent with the antiferromagnetic S = 1/2 Ising model on the triangular lattice [41][42][43]. The results obtained in these works do not predict long-range order at zero field (i.e.,μ = 0) [41]; however, longer-range interactions [43] or weak interactions between different (triangular) layers occurring in realistic systems could stabilize such an order (cf. also Ref. [33]).
Increasing U should decrease the temperatures of the order-disorder transitions in similar way as for alternate lattices, at least for U < 0 [27,30]. The analyses of this issue will be presented in the future work.
Additionally, one can note that the mean-field results for the atomic limit of the extended Hubbard model with attractive W < 0 are the same for both two-sublattice and tri-sublattice assumptions. In such a case, two different nonordered phases exist with the discontinuous first-order transition atμ = 0, and thus for fixed n, the electron-droplet state (phase separation NO/NO) exists (cf. Refs. [26,29,30] for large attractive U and Ref. [41] for the ferromagnetic Ising model).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.