Charge-Order on the Triangular Lattice: A Mean-Field Study for the Lattice S = 1/2 Fermionic Gas

The adsorbed atoms exhibit tendency to occupy a triangular lattice formed by periodic potential of the underlying crystal surface. Such a lattice is formed by, e.g., a single layer of graphane or the graphite surfaces as well as (111) surface of face-cubic center crystals. In the present work, an extension of the lattice gas model to S=1/2 fermionic particles on the two-dimensional triangular (hexagonal) lattice is analyzed. In such a model, each lattice site can be occupied not by only one particle, but by two particles, which interact with each other by onsite U and intersite W1 and W2 (nearest and next-nearest-neighbor, respectively) density-density interaction. The investigated hamiltonian has a form of the extended Hubbard model in the atomic limit (i.e., the zero-bandwidth limit). In the analysis of the phase diagrams and thermodynamic properties of this model with repulsive W1>0, the variational approach is used, which treats the onsite interaction term exactly and the intersite interactions within the mean-field approximation. The ground state (T=0) diagram for W2≤0 as well as finite temperature (T>0) phase diagrams for W2=0 are presented. Two different types of charge order within 3×3 unit cell can occur. At T=0, for W2=0 phase separated states are degenerated with homogeneous phases (but T>0 removes this degeneration), whereas attractive W2<0 stabilizes phase separation at incommensurate fillings. For U/W1<0 and U/W1>1/2 only the phase with two different concentrations occurs (together with two different phase separated states occurring), whereas for small repulsive 0


Introduction
It is a well known fact that the classical lattice gas model is useful phenomenological model for various phenomena. It has been studied in the context of experimental studies of adsobed gas layers on crystaline substrates (cf., for example pioneering works [1][2][3][4]). For instance, the adsorbed atoms exhibit tendency to occupy a triangular lattice formed by periodic potential of the underlying crystal surface. This lattice is shown in Figure 1a. Such a lattice is formed by, e.g., a single layer of graphane or the graphite surface [i.e., the honeycomb lattice; (0001) hexagonal closed-packed (hcp) surface], and (111) facecentered cubic (fcc) surface. Atoms from (111) fcc surface are organized in the triangular lattice, whereas the triangular lattice is a dual lattice for the honeycomb lattice [5]. Note also that arrangements of atoms on (110) base-centered cubic (bcc) surface as well as on (111) bcc surface (if one neglects the interactions associated with other layers under surface) are quite close to the triangular lattice. One should mention that the triangular lattice and the honeycomb lattice are two examples of two-dimensional hexagonal Bravais lattices. Formally, the triangular lattice is a hexagonal lattice with a one-site basis, whereas the honeycomb lattice is a hexagonal lattice with a two-site basis. The classical lattice gas model is equivalent with the S = 1/2 Ising model in the external field [1,[6][7][8][9] (the results for this model on the triangular lattice will be discussed in more details in Section 2).
In the present work, an extension of the lattice gas model to S = 1/2 fermionic particles is analyzed. Such a model has a form of the atomic limit of the extended Hubbard model [10], cf. (1). In this model, each lattice site can be occupied not by only one particle as in the model discussed in previous paragraph, but also by two particles. In addition to long-range (i.e., intersite) interactions between fermions, the particles located at the same site can also interact with each other via onsite Hubbard U interaction. For a description of the interacting fermionic particles on the lattice, the single-orbital extended Hubbard model with intersite density-density interactions has been used widely [11][12][13][14][15][16][17][18][19]. It is one of the simplest model capturing the interplay between the Mott localization (onsite interactions) and the charge-order phenomenon [15][16][17][18][19][20][21][22][23][24][25]. However, in some systems the inclusion of other interactions and orbitals is necessary [11][12][13][14]26,27]. This work can be palced among recent theoretical and experimental studies of adsorption of various atoms on (i) the (0001) hcp surface of the graphite [28][29][30][31][32] or of other materials [33,34] and (ii) the (111) fcc surface of metals and semimetals [34][35][36][37][38][39]. Although in the mentioned works the adsorbed particles on surface are rather classical and the analysis of classical lattice gas can give some predictions, taking into account of the quantum properties of adsorbed particles is necessary for, e.g., a description of experiments with He 4 and He 3 [28,40,41]. Moreover, there is plethora of recent experimental and theoretical studies of quasi-two-dimensional systems, e.g., Na x CoO 2 [42], NbSe 2 [43][44][45][46][47], TiSe 2 [48], TaSe 2 [49], VSe 2 [50], TaS 2 [51], and other transition metals dichalcogenides [52] as well as organic conductors [53,54], where various charge-ordered patterns have been observed on the triangular lattice. However, for such phenomena the atomic limit of the model studies is less reliable and more realistic description includes also electron hoping term as in the extended Hubbard model [15][16][17][18][19] or coupling with phonons as in the Holstein-Hubbard model [55]. In such cases, results obtained for atomic limit can be treated as a benchmark for models including the itinerant properties of fermionic particles.
The present work is organized as follows. In Section 2 the model and the methods (together with the most important equations) are presented. Section 3 is devoted to the discussion of ground state phase diagrams of the model with non-zero next-nearest neighbor interactions. Next, the finite temperature properties of the model with only the nearest-neighbor interactions are presented in Section 4. Finally, the most important conclusions and supplementary discussion are included in Section 5.

The Model and the Method
The extended Hubbard model in the zero-bandwidth limit (i.e., in the atomic limit) with interactions restricted to the second neighbors (or, equivalently, the next-nearest neighbors) can be expressed as: wheren i = ∑ σniσ ,n iσ =ĉ † iσĉ iσ , andĉ † iσ (ĉ iσ ) denotes the creation (annihilation) operator of an electron with spin σ at the site i. U is the onsite density interaction, W 1 and W 2 are the intersite density-density interactions between the nearest neighbors (NN) and the next-nearest neighbors (NNNs), respectively. z 1 and z 2 are numbers of NN and NNNs, respectively. µ is the chemical potential determining the total concentration n of electrons in the system by the relation n = (1/L) ∑ i n i , where 0 ≤ n ≤ 2 and L is the total number of lattice sites. In this work phase diagrams emerging from this model are inspected. The analyses are performed in the grand canonical ensemble.
In [73] the model with W 2 = 0 was investigated on the triangular lattice at half-filling by using a classical Monte Carlo method, and a critical phase, characterized by algebraic decay of the charge correlation function, belonging to the universality class of the twodimensional XY model with a Z 6 anisotropy was found in the intermediate-temperature regime. Some preliminary results for model (1) on the triangular lattice and for large attractive U < 0 and W 2 = 0 within the mean-field approximation were presented in [74].
The model in the limit U → −∞ is equivalent with the S = 1/2 Ising model with antiferromagnetic (ferromagnetic) J n interactions if W n interaction in model (1) are repulsive, i.e., W n > 0 (attractive, i.e., W n < 0, respectively). The relation between interaction parameters in both models is very simple, namely J n = −W n . There is plethora of the results obtained for the Ising model on the triangular lattice. One should mention the following works (not assuming a comprehensive review): (a) exact solution in the absence of the external field H, i.e., for H = 0 (only with NN interactions, at arbitrary temperature) [5,[75][76][77][78][79]; (b) for the model with NNN interactions included: ground state exact results [2], Bethe-Peierls approximation [1], Monte Carlo simulation both for H = 0 [80] and H = 0 [3] (and other methods, e.g., [81,82]); (c) exact ground state results for the model with up to 3rd nearest-neighbor interactions for both H = 0 case [83] and H = 0 case [4,84]. The most important information arising from these analyses is that only for W 2 ≤ 0 (and arbitrary W 1 ) one can expect that consideration of √ 3 × √ 3 unit cells (i.e., tri-subblatice orderings) is enough to find all ordered states (particle arrangements) in the model. The reason is that the range of W 2 interaction is larger than the size of the unit cell. Thus, this is the point for that the present analysis of the model including only √ 3 × √ 3 unit cell orderings with restriction to W 2 ≤ 0 is justified. One should not expect occurrence of any other phases beyond the tri-sublattice assumption in the studied range of the model parameters.
Please note that for W 2 > 0 it is necessary to consider a larger unit cell to find the true phase diagram of the model even in the U → −∞ limit (cf., e.g., [2,4,81]). This is a similar situation as for model (1) on the square lattice, where for W 2 > 0 and any U not only checker-board order occurs (the two-sublattice assumption), but also other different arrangements of particles are present (the four-sublattice assumption, e.g., various stripes orders) [67,68].

General Definitions of Phases Existing in the Investigated System
In the systems analyzed only three nonequivalent homogeneous phases can exist (within the tree-sublattice assumption used). They are determined by the relations between concentrations n α 's in each sublattice α (n α = (3/L) ∑ i∈α n i ), but a few equivalent solutions exist due to change of sublattice indexes. For intuitive understanding of rather complicated phase diagrams each pattern is marked with adequate abbreviation. The nonordered (NO) phase is defined by n A = n B = n C (all three n α 's are equal), the charge-ordered phase with two different concentrations in sublattices (DCO phase) is defined by n A = n B = n C , n B = n C = n A , or n C = n A = n B (two and only two out of three n α 's are equal, 3 equivalent solution), whereas in the charge-ordered phase with three different concentrations in sublattices (TCO phase) n A = n B , n B = n C , and n A = n C (all three n α 's are different, 6 equivalent solutions). All these phases are schematically illustrated in Figure 1b. These phases exist in several equivalent solutions due to the equivalence of three sublattices forming the triangular lattice. Each of these patterns can be realized in a few distinct forms depending on specific electron concentrations on each sublattice (cf. Tables 1 and 2 for T = 0). In addition, the degeneracy of the ground state solutions is contained in Table 1 (including charge and spin degrees of freedom).

Expressions for the Ground State
In the ground state (i.e., for T = 0), the grand canonical potential ω 0 per site of model (1) can be found as where contributions associated with the onsite interaction, the intersite interactions, and the chemical potential, respectively, has the following forms In the above expressions, concentrations n α at T = 0 take the values from {0, 1, 2} set (cf. also Table 1). Please note that the above equations are the exact expressions for ω 0 of model (1) on the triangular lattice.
The free energy per site of homogeneous phases at T = 0 within the mean-field approximation is obtained as where E W is expressed by (5). D occ = (1/L) ∑ i n i↑ni↓ denotes the double occupancy and this quantity is found to be exact, cf. Table 2. One should underline that above expression for f 0 is an approximate result for model (1) on the triangular lattice. Here, it is assumed that concentration n α are as defined in Table 2 and they are the same in each √ 3 × √ 3 unit cell in the system. Formally, it could be treated as exact one only if the numbers z n (n = 1, 2) goes to infinity.
The expressions presented in this subsection (for W 2 = 0) can be obtained as the T → 0 limit of the equations for T > 0 included in Section 2.3.

Expressions for Finite Temperatures
For finite temperatures (T > 0), the expressions given in [59] for the three-sublattice assumption takes the following forms (cf. also these in [68] given for the four-sublattice assumption). In approach used, the onsite U term is treated exactly and for the intersite W 1 term the mean-field approximation (2) is used. For a grand canonical potential ω (per lattice site) in the case of the lattice presented in Figure 1 one obtains where and µ α is a local chemical potential in α sublattice (α ∈ {A, B, C}) For electron concentration n α in each sublattice in arbitrary temperature T > 0 one gets The set of three Equations (11) for n A , n B , and n C determines the (homogeneous) phase occurring in the system for fixed model parameters U, W 1 , and µ. If n = (1/3)(n A + n B + n C ) is fixed, one has also set of three equations, but it is solved with respect to µ, n A , and n B (the third n α is obviously found as n C = 3n − n A − n B ).
The free energy f per site is derived as where ω and n α 's are expressed by (8)-(11).

Macroscopic Phase Separation
The free energy f PS of the (macroscopic) phase separated state (as a function of total electron concentration n; and at any temperature T ≥ 0) is calculated from where f ± (n ± ) are free energies of separating homogeneous phases with concentrations n ± . The factor before f ± (n ± ) is associated with a fraction of the system, which is occupied by the phase with concentration n ± . Such defined phase separated states can exist only for n fulfilling the condition n − < n < n + . For n ± only the homogeneous phase exists in the system (one homogeneous phase occupies the whole system). Concentrations n ± are simply determined at the ground state, whereas for T > 0 they can be found as concentrations at the first-order (discontinuous) boundary for fixed µ or by minimizing the free energy f PS [i.e., (13)] with respect to n + and n − (for n fixed). For more details of the so-called Maxwell's construction and macroscopic phase separations see, e.g., [59,68,85,86]. The interface energy between two separating phases is neglected here. Table 1. Homogeneous phases (z n → ∞, n = 1, 2) or √ 3 × √ 3 unit cells (the triangular lattice) at T = 0 (for fixed µ). Star " * " in superscript indicates that the phase is obtained by the particle-hole transformation (i.e., n α → 2 − n α ; the NO 1 and TCO phases are invariant under this transformation). In the brackets also an alternative name is given. The degeneration d c × d s of the unit cells (equal to the degeneration of the ground state for z n → ∞ limit) and degeneration D c × D s of the ground state phases constructed from the corresponding unit cells for the triangular lattice is given (with respect to charge and spin degrees of freedom). Table 2. Homogeneous phases at T = 0 (for fixed n) defined by n α 's and D occ . n s and n f define the range [n s , n f ] of n, where the phase is correctly defined. In the last column, the phase separated state degenerated with the homogeneous phase in range (n s , n f ) for W 2 = 0 is mentioned. Star " * " in superscript indicates that the phase is obtained by the particle-hole transformation (i.e., n α → 2 − n α ; TCO A , TCO * A , TCO B , and TCO * B phases are invariant under this transformation).
Phase n A n B n C D occ n s n f PS The ground state diagram for model (1) as a function of (shifted) chemical potential µ = µ − W 1 − W 2 is shown in Figure 2. The diagram is determined by comparison of the grand canonical potentials ω 0 's of all phases collected in Table 1 [cf. (3)]. It consists of several regions, where the NO phase occurs (3 regions: NO 0 , NO 1 and NO 2 ), the DCO phase occurs (6 regions: DCO 1 , DCO 2 , DCO 3 , DCO * 1 , DCO * 2 , and DCO * 3 ) and the TCO phase occurs (1 region).
All boundaries between the phases in Figure 2 are associated with a discontinuous change of at least one of the n α . The only boundaries associated with a discontinuous jump of two n α 's are: DCO 2 -DCO 3 (DCO * 2 -DCO * 3 ) and TCO-NO 1 . At the boundaries ω 0 's of the phases are the same. It means that both phases can coexist in the system provided that a formation of the interface between two phases does not require additional energy. For W 2 = 0, only the boundaries DCO 2 -DCO 3 (DCO * 2 -DCO * 3 ) and TCO-NO 1 have finite degeneracy (6 and 7, respectively, modulo spin degrees of freedom) and the interface between different types of √ 3 × √ 3 unit cells increases the energy of the system. Thus, the mentioned phases from neighboring regions cannot coexist at the boundaries. The other boundaries exhibit infinite degeneracy (it is larger than 3 · 2 L/3 modulo spin) and entropy per site in the thermodynamic limit is non-zero. It means that at these boundaries both types of unit cells from neighboring regions can mix with any ratio and the formation of the interface between two phases does not change energy of the system. However, some conditions for arrangement of the cells can exist. For example, the DCO 2 phase with (0, 0, 2) can mix with the DCO * 2 phase with (0, 2, 2) or (2, 0, 2), but not with the DCO * 2 phase with (2, 2, 0). Please note that it is also possible to mix all three unit cells: (0, 0, 2), (0, 2, 2), and (2, 0, 2). In such a case, (0, 2, 2) and (2, 0, 2) cells of the DCO * 2 phase cannot be located next to each other, i.e., they need to be separated by (0, 0, 2) unit cells of the DCO 2 phase. Thus, the degeneracy of the DCO 2 -DCO * 2 boundary is indeed larger than 3 · 2 L/3 modulo spin. This is so-called macroscopic degeneracy, cf. [68]). In such a case, we say that the microscopic phase separation occurs. For W 2 < 0 these degeneracies are removed and all boundaries exhibit finite degeneracy (neglecting spin degrees of freedom). In this case the phases cannot be mixed on a microscopic level.
. The regions are labeled by the names of the phases defined in Table 1 and numbers corresponding to concentrations in each sublattice n A , n B and n C .
Please also note that for W 2 = 0 as well as for W 2 < 0 inside the regions shown in Figure 2, the √ 3 × √ 3 unit cells of the same type with different orientation cannot mix. It denotes that orientation of one type of the unit cell determines the orientation of other unit cells (of the same type). Thus, the degeneracy of the state of the system is finite (modulo spin) and the system exhibits the long-range order at the ground state inside each region of Figure 2. This is different from the case of two dimensional square lattice, where inside some regions different unit cells (elementary blocks) of the same phase can mix with each other [67,68].
One should underline that the discussed above ground state results for fixed chemical potential are the exact results for model (1) on the triangular lattice. This is due to the fact that the model is equivalent with a classical spin model, namely the S = 1 Blume-Cappel model with two-fold degenerated value of S = 0 (or the S = 1 classical Blume-Cappel with temperature-dependent anizotropy without degeneration), cf. [10,60,63]. For such a model, the mean-field approximation is an exact theory at the ground state and fixed external magnetic field (which corresponds to the fixed chemical potential in the model investigated).

Analysis for Fixed Particle Concentration n
The ground state diagram as a function of particle concentration n is shown in Figure 3. The rectangular regions are labeled by the abbreviations of homogeneous phases (cf. Table 2). At commensurate filling, i.e., i/3 (i = 0, 1, 2, 3, 4, 5, 6; but only on the vertical boundaries indicated in Figure 3) the homogeneous phase occurs, which can be found in Table 1 and Figure 2. On the horizontal boundaries the phases from both neighboring regions have the same energies.
For W 2 = 0 phase separated states (mentioned in the last column of Table 2) are degenerated with the corresponding homogeneous phases inside all regions of the phase diagram. This degeneracy can be removed in finite temperatures and in some regions the phase separated states can be stable at T > 0 (such regions are indicated by slantwise patter in Figure 3, cf. also Section 4). E.g., for W 2 = 0, the TCO phase can exist only in the range of 0 < U/W 1 < 1/2 at T = 0. For W 2 < 0 the phase separated states have lower energies and they occur on the phase diagram (inside the rectangular regions of Figure 3). Obviously, at commensurate filling and for any W 2 ≤ 0, the homogeneous states can only occur (i.e., solid vertical lines in Figure 3). Please note that the following boundaries between homogeneous states (obtained by comparing only energies of homogeneous phases): (i) the DCO A and DCO B phases, (ii) the DCO A and DCO C phases, and (iii) the TCO A and TCO B phases are located at U/W 1 = 0 (and these corresponding for n > 1; the dashed line in Figure 3). For W 2 < 0 these lines do not overlap with the boundaries between corresponding phase separated states at U/W 1 − |k| = 0 (or U/|W 2 | = 1), but in such a case the homogeneous states have higher energies than the phase separated states. In fact, the homogeneous states for W 2 < 0 are unstable (i.e., ∂µ/∂n < 0) inside the regions of Figure 1. For W 2 < 0 they are stable only for commensurate fillings (solid lines in Figure 3).
For the system on the square lattice the similar observation can be made (Figure 1 from [59])-compare HCO A -LCO A and HCO A -HCO B boundaries at U/W 1 = 0 with PS1 A -PS1 B and PS1 A -PS1 B boundaries at U/|W 2 | = 1, respectively. In [68] the boundaries between homogeneous phases for W 2 < 0 are not shown in Figure 3. Only boundaries between corresponding phase separated states are correctly presented in that figure for W 2 < 0. For U/W 1 > 0, the CBO A phase (corresponding to the HCO A phase from [59]) is not the phase with the lowest energy (among homogeneous phases) in any range of n (but for U/W 1 < 0 it has the lowest energy among all homogeneous states). However, the corresponding phase separated state NO 0 /CBO 2 (i.e., PS1 A from [59]) can occur for U/W 1 > 0 (and for U/|W 2 | < 1) as shown in Figure 3 of [68].
The vertical boundaries for homogeneous phases (i.e., the transitions with changing n) are associated with continuous changes of all n α 's and D occ , but the chemical potential µ (calculated as µ = ∂ f /∂n) changes discontinuously. Boundaries DCO A -DCO B , DCO A -DCO C , and TCO A -TCO B (and other transitions for fixed n at U/W 1 = 0) between homogeneous phases are associated with discontinuous change of only D occ . One should note that it is similar to transition between two checker-board ordered phases on the square lattice, namely CBO A -CBO B and CBO A -CBO C boundaries, cf. [68] (or the HCO A -LCO A and HCO A -HCO B boundaries, respectively, from [59]). At the other horizontal boundaries (i.e., transitions for fixed n at U/W 1 − |k| = 1/2 in Figure 3) two of n α 's and D occ change discontinuously. At commensurate fillings transitions with changing U/W 1 occur only at points indicated by squares in Figure 3.
All horizontal boundaries between phase separated states (which are stable for W 2 < 0) are connected with discontinuous changes of D occ . These boundaries located at U/W 1 − |k| = 0 are also associated to a discontinuous change of particle concentration in one of the domains.
The diagram presented in Figure 3 is constructed by the comparison of (free) energies of various homogeneous phases and phase separated states collected in Table 1. The energies of homogeneous phases are calculated from (7), whereas energies of phase separated states are calculated from (13). Please note that it is easy to calculate energies of f ± (n ± ) of separating homogeneous phase (with commensurate fillings) at the ground state by just taking µ = 0 in ω 0 's collected in Table 1. Obviously, one can also calculate energies of the phases collected in Table 2 at these fillings (from both neighboring regions). For example, the DCO B phase and the DCO C phase at n = 1/3 reduce to DCO 1 phase. Ground state phase diagram of the model as a function of particle concentration n for W 1 > 0 and W 2 ≤ 0 (|k| = |W 2 |/W 1 ). The regions are labeled by the names of the homogeneous phases (cf. Table 2). For W 2 = 0 all homogeneous phases are degenerated with macroscopic phase separated states indicated in the last column of Table 2. In regions filled by slantwise pattern the phase separated states occurs at infinitesimally T > 0 for W 2 = 0. For W 2 < 0 the phase separated states occur inside the regions, whereas at the vertical boundaries for commensurate filling the homogeneous states (defined in Table 1) still exist. The boundary at U/W 1 = 0 (schematically indicated by dashed green line) denotes the boundaries between homogeneous phases, which do not overlap with the boundaries between phase separated states for W 2 < 0. Squares denote transitions for fixed n between homogeneous phase at commensurate fillings.

Results for Finite Temperatures (W 1 > 0 and W 2 = 0)
One can distinguish four ranges of U interaction, where the system exhibits qualitatively different behavior, namely: (i) U/W 1 < 0, (ii) 0 < U/W 1 < (1/3) ln(2), (iii) (1/3) ln(2) < U/W 1 < 1/2, and (iv) U/W 1 > 1/2. In Figures 4-7, the exemplary finite temperature phase diagrams occurring in each of these ranges of onsite interaction are presented. All diagrams are found by investigation of the behavior of n α 's determined by (11) in the solution corresponding to the lowest grand canonical potential [Equation (8), when µ is fixed] or to the lowest free energy [Equations (12) and (13) if n is fixed]. The set of three nonlinear Equations (11) has usually several nonequivalent solutions and thus it is extremely important to find a solution, which has the minimal adequate thermodynamic potential. In Figure 8 the behavior of n α 's as a function of temperature or chemical potential is shown for some representative model parameters. Figure 9 presents the phase diagram of the system for half-filling.  For U/W 1 < 0 and U/W 1 > 1/2 the phase diagrams of the model are similar and the DCO phase is only ordered homogeneous one occurring on the diagrams. In the first range, there are two regions of ordered phase occurrence (cf. Figure 4 and [74]), whereas in the second case one can distinguish four regions of the DCO phase stability (cf. Figure 5). The NO-DCO transitions for fixed µ are discontinuous for any values of onsite interaction and chemical potential in discussed range of model parameters and thus phase separated state PS 1 :NO/DCO occurs in define ranges of n. In this state domains of the NO and the DCO phases coexist.
For U/W 1 < 0 the temperature of NO-DCO transition is maximal forμ = 0 (i.e., at half-filling)- Figure 4a. Its maximal value T M monotonously decreases with increasing of U from k B T M /W 1 = 1/2 for U → −∞ and at U = 0 it is equal to 1/4. This transition exhibits re-entrant behavior (for fixed |μ| > 1). At T = T M andμ = 0 and at only this point, this transition exhibits properties of a second order transition (cf. Figure 8a). In particular, with increasing T for fixedμ = 0 n α 's changes continuously at T M , but two equivalent solutions still exist for any T < T M (similarly as in the ferromagnetic Ising model at zero field [9]). Atμ = 0 and T < T M the discontinuous transition between two DCO phases occurs. In the DCO phase forμ < 0 (n < 1) [connecting with the DCO 1 (DCO A ) region at T = 0] the relation n A = n B < n C is fulfilled, whereas in the DCO phase forμ > 0 (n > 1) [connecting with the DCO * 1 (DCO * A ) region at T = 0] the relation n A < n B = n C occurs (n C can be larger than 1 for some temperatures), cf. also Figure 8g,h as well as [74]. Both discontinuous transitions for fixed chemical potential are associated with occurrence of phase separated states. On the diagrams obtained for fixed n three region of phase separated states occurs (Figure 4b). For W 2 = 0 the PS 1 :NO/DCO phase separated state occurs only for T > 0. For T → 0 the concentrations in both domains of the PS 1 state approach 0 (or 2), whereas for T → T M they approach to 1. Near n = 1 the PS 2 :DCO/DCO state is stable for 0 ≤ T < T M . In this state domains of two DCO phases (with different particle concentrations) coexist in the system.  Figure 6. Phase diagrams of the model for U/W 1 = 0.20 as a function of (a) chemical potential µ/W 1 and (b) particle concentration n (W 1 > 0, W 2 = 0). The boundary TCO-DCO is second order, the remaining are first order. Other denotations as in Figure 4. The diagrams are shown only for µ ≤ 0 and n ≤ 1, but they are symmetric with respect toμ = 0 and n = 1, respectively.  Figure 7. Phase diagrams of the model for U/W 1 = 0.35 as a function of (a) chemical potential µ/W 1 and (b) particle concentration n (W 1 > 0, W 2 = 0). The boundary TCO-DCO is second order, the remaining are first order. Other denotations as in Figure 4. The diagrams are shown only for µ ≤ 0 and n ≤ 1, but they are symmetric with respect toμ = 0 and n = 1, respectively.
For U/W 1 > 1/2 the diagrams are similar, but the double occupancy of sites is strongly reduced due to repulsive U ( Figure 5). Thus, their structure exhibits two lobs of the DCO phase occurrence in cotrary to the case of U/W 1 < 0, where a single lob of the DCO phase is present (as expected from previous studies of the model, cf. [10,59,60]). The maximal value k B T M /W 1 of NO-DCO transition occurs forμ/W 1 corresponding approximately quarter fillings (i.e., near n = 1/2 and n = 3/2). With increasing U it decreases and finally in the limit U → +∞ it reaches 1/8. At this point DCO-NO boundary exhibits features of continuous transition as discussed previously. In this range, the phase diagrams are (almost) symmetric with respect to these fillings (when one considers only one part of the diagram for 0 < n < 1 or for 1 < n < 2).
The most complex diagrams are obtained for 0 < U/W 1 < 1/2, where the TCO phase appears at T = 0 and for finite temperatures near half-filling. For 0 < U/W 1 < (1/3) ln(2) the region of the TCO phase is separated from the NO phase by the region of DCO phase, Figure 6a. The TCO-DCO transition is continuous (cf. Figure 8g,h for U/W 1 = 0.35) and its maximal temperature is located for half-filling (atμ = 0 or n = 1). At this point two firstorder NO-DCO and two second-order TCO-DCO boundaries merge (for fixed chemical potential). It is the only point for fixed U/W 1 in this range of model parameters, where a direct continuous transition from the TCO phase to the NO phase is possible (Figure 8b). The continuous TCO-DCO transition temperature can be also found as a solution of (11) and (A4) as discussed in Appendix A. Similarly as for U/W 1 < 0, the temperature of NO-DCO transition is maximal at half-filling. For fixed n, the narrow regions of PS 1 :NO/DCO states are present between the NO region and DCO regions. Please note that for T > 0 there is no signatures of the discontinuous DCO 1 -DCO 2 (DCO * 1 -DCO * 2 ) boundary occurring at T = 0. It is due to the fact that the discontinuous jumps of n α 's occurring for T = 0 at these boundaries are changed into continuous evolutions of sublattice concentrations at T > 0 and there is no criteria for distinction of these two DCO phases at finite temperatures (cf. also [59][60][61]). From the same reason, there is no boundary at T > 0 for fixed n associated to the DCO B -DCO C (DCO * B -DCO * C ) line occurring at T = 0 ( Figure 6b). However, strong reduction of one n α from the case where n α ≈ 2 to the case of n α ≈ 1 is visible (some kind of a smooth crossover inside the DCO region), cf. For (1/3) ln(2) < U/W 1 < 1/2, the maximum of the NO-DCO transition temperature is shifted towards larger |μ|/W 1 (or smaller |1 − n|). This is associated with forming of the two-lob structure of the diagram found for U/W 1 > 1/2. Inside the regions of the DCO phase occurrence discontinuous transitions between two DCO phases appear-See Figure 7a as well as Figure 8e,i. These new regions of the DCO phase at T > 0 [with n A < n B = n C (forμ < 0 or n < 1); cf. Figure  . From the analysis of (11) similarly as it was done in the case of the square lattice [10] (see also Appendix A) one obtains that the point, where the TCO-NO transition changes its order at half-filling, is k B T/W 1 = 1/6 and U/W 1 = (1/3) ln(2). For better overview of the system behavior, the phase diagram of the model for halffilling (μ = 0 or n = 1) is presented in Figure 9. The temperature of the order-disorder transition decreases with increasing U/W 1 . In low temperatures and for U/W 1 < 0, the DCO phases exist in the system (precisely, if µ is fixed-atμ = 0 the DCO-DCO discontinuous boundary occurs; whereas if n is fixed-the PS 2 :DCO/DCO state is stable at n = 1), cf. also Figure 4. For 0 < U/W 1 < 1/2 the TCO phase is stable below the order-disorder line, but for (1/3) ln(2) < U/W 1 < 1/2 and k B T/W 1 < 1/6 the TCO-NO phase transition is discontinuous (cf. also Figure 7). For U/W 1 < (1/3) ln(2) the order-disorder boundary presented in Figure 9 is a merging point of several boundaries as presented in Figures 4 and 6, and discussed previously. Thus, formally this order-disorder boundary for U/W 1 < (1/3) ln(2) occurring at half-filling is a line of some critical points of a higher order.
Please note that the order-disorder transition is discontinuous for any value of onsite interaction and chemical potential [excluding only the TCO-NO boundary for half-filling and 0 < U/W 1 < (1/3) ln (2)] in contrast to the case of two- [10,59,60] or four-sublattice [67,68] assumptions, where it can be continuous one for some range of model parameters). In [74] also metastable phases have been discussed in detail for the large onsite attraction limit and the triangular lattice.

Final Remarks
In this work, the mean-field approximation was used to investigate the atomic limit of extended Hubbard model [hamiltonian (1)] on the triangular lattice. The phase diagram was determined for the model with intersite repulsion between the nearest neighbors (W 1 > 0). The effects of attractive next-nearest-neighbor interaction (W 2 < 0) were discussed in the ground state. The most important findings of this work are that (i) two different arrangements of particles (i.e., two different charge-ordered phases: the DCO and TCO states) can occur in the system and (ii) attractive W 2 < 0 or finite T > 0 removes the degeneration between homogeneous phases and phase separated states occurring at T = 0 for W 2 = 0. It was shown that TCO phase is stable in intermediate range of onsite repulsion 0 < U/W 1 < 1/2 (for W 2 = 0). All transition from the ordered phases to the NO are discontinuous for fixed chemical potential (apart from TCO-NO boundary at half-filling for 0 < U/W 1 < (1/3) ln(2)) and the DCO-NO boundaries at single points corresponding to n = 1/2, 1, 3/2 as discussed in Section 4), thus the phase separated states occur on the phase diagram for fixed particle concentration.
One should stress that hamiltonian (1) is interesting not only from statistical point of view as a relatively simple toy model for phase transition investigations. Although it is oversimplified for quantitative description of bulk condensed matter systems, it can be useful in qualitative analysis of, e.g., experimental studies of adsorbed gas layers on crystalline substrates.
Notice that the mean-field decoupling of the intersite term is an approximation for purely two-dimensional model investigated, which overestimates the stability of ordered phases. For example, the order-disorder transition for the ferromagnetic Ising model is overestimated by the factor two (for the honeycomb, square and triangular lattices rigorous solution gives k B T c /|J| as 0.506, 0.568, 0.607, respectively, whereas the mean-field approximation gives k B T c /|J| = 1) [76]. Moreover, the results for the antiferromagnetic Ising model on the triangular lattice [the limit U → ±∞ of model (1)] do not predict long-range order at zero field [1,3,76] and T > 0 [corresponding to n = 1 or n = 1/2, 3/2, respectively, in the case of model (1)]. However, longer-range interactions [3] or weak interactions between adsorbed particles and the adsorbent material occurring in realistic systems could stabilize such an order (such systems are rather quasi-two-dimensional). It should be also mentioned that the charge Berezinskii-Kosterlitz-Thouless-like phase was found in the intermediate-temperature regime between the charge-ordered phase (with long-range order, coresponding to the TCO phase here) and disordered phases in the investigated model [73].
The recent progress in the field of optical lattices and a creation of the triangular lattice by laser trapping [89,90] could enable testing predictions of the present work. The fermionic gases in harmonic traps are fully controllable systems. Note also that the superconductivity in the twisted-bilayer graphene [91][92][93][94][95][96] is driven by the angle between the graphene layers.
It is associated with an occurrence of the Moiré pattern (the triangular lattice with very large supercell). Hetero-bilayer transition metals dichalcogenides system is the other field where this pattern appears [97,98]. This makes further studies of properties of different models on the triangular lattice desirable. Informed Consent Statement: Not applicable. The study does not involve humans.

Data Availability Statement:
The data presented in this study are available on request from the author. All data presented in Section 4 have been obtained by numerical solving of self-consistent equations given in Section 2.3.

Acknowledgments:
The author expresses his sincere thanks to J. Barański, R. Lemański, R. Micnas, P. Piekarz, and A. Ptok for very useful discussions on some issues raised in this work. The author also thanks R. Micnas and I. Ostrowska for careful reading of the manuscript.

Conflicts of Interest:
The author declares no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

NO
Non-ordered phase with the same concentrations in all three sublattices DCO Charge-ordered phase with two different concentrations in sublattices TCO Charge-ordered phase with three different concentrations in sublattices