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

**Konrad Jerzy Kapcia**

Faculty of Physics, Adam Mickiewicz University in Pozna ´n, ulica Uniwersytetu Pozna ´nskiego 2, PL-61614 Pozna ´n, Poland; konrad.kapcia@amu.edu.pl; Tel.: +48-61-829-5051

**Abstract:** 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 *W*<sup>1</sup> and *W*<sup>2</sup> (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 *W*<sup>1</sup> > 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 *W*<sup>2</sup> ≤ 0 as well as finite temperature (*T* > 0) phase diagrams for *W*<sup>2</sup> = 0 are presented. Two different types of charge order within <sup>√</sup> 3 × √ 3 unit cell can occur. At *T* = 0, for *W*<sup>2</sup> = 0 phase separated states are degenerated with homogeneous phases (but *T* > 0 removes this degeneration), whereas attractive *W*<sup>2</sup> < 0 stabilizes phase separation at incommensurate fillings. For *U*/*W*<sup>1</sup> < 0 and *U*/*W*<sup>1</sup> > 1/2 only the phase with two different concentrations occurs (together with two different phase separated states occurring), whereas for small repulsive 0 < *U*/*W*<sup>1</sup> < 1/2 the other ordered phase also appears (with tree different concentrations in sublattices). The qualitative differences with the model considered on hypercubic lattices are also discussed.

**Keywords:** charge order; triangular lattice; extended Hubbard model; atomic limit; mean-field theory; phase diagram; longer-range interactions; thermodynamic properties; fermionic lattice gas; adsorption on the surface

**1. 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–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

**Citation:** Kapcia, K.J. Charge-Order on the Triangular Lattice: A Mean-Field Study for the Lattice *S* = 1/2 Fermionic Gas. *Nanomaterials* **2021**, *11*, 1181. https://doi.org/ 10.3390/nano11051181

Academic Editor: Orion Ciftja

Received: 25 March 2021 Accepted: 29 April 2021 Published: 30 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

model is equivalent with the *S* = 1/2 Ising model in the external field [1,6–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–19]. It is one of the simplest model capturing the interplay between the Mott localization (onsite interactions) and the charge-order phenomenon [15–25]. However, in some systems the inclusion of other interactions and orbitals is necessary [11–14,26,27].

**Figure 1.** (**a**) The schema of the triangular lattice on which the extended Hubbard model in the atomic limit is studied in the present work. The lattice is divided into three equivalent sublattices (*α* = *A*, *B*, *C*) denoted by different symbols. The dash-dotted line denotes the boundaries of <sup>√</sup> 3 × √ 3 unit cell. By solid and dashed lines all nearest neighbors and all next-nearest neighbors of a chosen site from sublattice *A* are indicated, respectively. (**b**) There different types of particle arrangements in <sup>√</sup> 3 × √ 3 unit cells (i.e., the tri-sublattice assumption) corresponding to NO, DCO, TCO phases (as labeled). Symbol shapes on each panel correspond to respective concentrations at the lattice sites.

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–32] or of other materials [33,34] and (ii) the (111) fcc surface of metals and semimetals [34–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<sup>4</sup> and He<sup>3</sup> [28,40,41]. Moreover, there is plethora of recent experimental and theoretical studies of quasi-two-dimensional systems, e.g., Na*x*CoO<sup>2</sup> [42], NbSe<sup>2</sup> [43–47], TiSe<sup>2</sup> [48], TaSe<sup>2</sup> [49], VSe<sup>2</sup> [50], TaS<sup>2</sup> [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–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.

#### **2. 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:

$$\hat{H}\_{\perp} = -\mathcal{U}\sum\_{\mathbf{i}} \hbar\_{\mathbf{i}\uparrow} \hbar\_{\mathbf{i}\downarrow} + \frac{1}{2} \frac{\mathcal{W}\_{\mathbf{1}}}{z\_{1}} \sum\_{\langle i,j\rangle\_{1}} \hbar\_{\mathbf{i}} \hbar\_{\mathbf{j}} + \frac{1}{2} \frac{\mathcal{W}\_{\mathbf{2}}}{z\_{2}} \sum\_{\langle i,j\rangle\_{2}} \hbar\_{\mathbf{i}} \hbar\_{\mathbf{j}} - \mu \sum\_{\mathbf{i}} \hbar\_{\mathbf{i}\prime} \tag{1}$$

where *n*ˆ*<sup>i</sup>* = ∑*<sup>σ</sup> n*ˆ*iσ*, *n*ˆ*i<sup>σ</sup>* = *c*ˆ † *iσ c*ˆ*iσ*, and *c*ˆ † *iσ* (*c*ˆ*iσ*) denotes the creation (annihilation) operator of an electron with spin *σ* at the site *i*. *U* is the onsite density interaction, *W*<sup>1</sup> and *W*<sup>2</sup> are the intersite density-density interactions between the nearest neighbors (NN) and the next-nearest neighbors (NNNs), respectively. *z*<sup>1</sup> and *z*<sup>2</sup> 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*) ∑*<sup>i</sup>* h*n*ˆ*i*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 this work the mean-field decoupling of the intersite term is used in the following form

$$
\hbar \hbar \mathfrak{H}\_{\hat{l}} = \langle \hbar\_{\hat{l}} \rangle \hbar\_{\hat{l}} + \hbar\_{i} \langle \hbar\_{\hat{l}} \rangle - \langle \hbar\_{i} \rangle \langle \hbar\_{\hat{l}} \rangle,\tag{2}
$$

which is an exact treatment only in the limit of large coordination number (*z<sup>n</sup>* → ∞; or limit of infinite dimensions) [10–13,56–58]. Thus, for the two-dimensional triangular lattice (with *z*<sup>1</sup> = *z*<sup>2</sup> = 6) it is an approximation in the general case. It should be underlined that the treatment of the onsite term is rigorous in the present work. Please note that that the interactions *U*, *W*<sup>1</sup> and *W*<sup>2</sup> should be treated as effective parameters for fermionic particles including all possible contributions and renormalizations originating from other (sub-)systems.

Model (1) for *W*<sup>2</sup> 6= 0 has been intensively studied on the hypercubic lattices (see, e.g., [59–68] and references therein). Also the case of two-dimensional square lattice was investigated in detail for *W*<sup>2</sup> = 0 [61–64] as well as for *W*<sup>2</sup> 6= 0 [65–69]. There are also rigorous results for one-dimensional chain for *W*<sup>2</sup> = 0 [70,71] and *W*<sup>2</sup> 6= 0 [72].

In [73] the model with *W*<sup>2</sup> = 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<sup>6</sup> 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*<sup>2</sup> = 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<sup>n</sup>* interactions if *W<sup>n</sup>* interaction in model (1) are repulsive, i.e., *W<sup>n</sup>* > 0 (attractive, i.e., *W<sup>n</sup>* < 0, respectively). The relation between interaction parameters in both models is very simple, namely *J<sup>n</sup>* = −*Wn*. 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–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* 6= 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* 6= 0 case [4,84]. The most important information arising from these analyses is that only for *<sup>W</sup>*<sup>2</sup> <sup>≤</sup> <sup>0</sup> (and arbitrary *<sup>W</sup>*1) one can expect that consideration of <sup>√</sup> 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*<sup>2</sup> 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*<sup>2</sup> ≤ 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*<sup>2</sup> > 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*<sup>2</sup> > 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].

#### *2.1. 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 *<sup>n</sup>α*'s in each sublattice *<sup>α</sup>* (*n<sup>α</sup>* = (3/*L*) <sup>∑</sup>*i*∈*<sup>α</sup>* h*n*ˆ*i*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<sup>A</sup>* = *n<sup>B</sup>* = *n<sup>C</sup>* (all three *nα*'s are equal), the charge-ordered phase with two different concentrations in sublattices (DCO phase) is defined by *n<sup>A</sup>* = *n<sup>B</sup>* 6= *nC*, *n<sup>B</sup>* = *n<sup>C</sup>* 6= *nA*, or *n<sup>C</sup>* = *n<sup>A</sup>* 6= *n<sup>B</sup>* (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<sup>A</sup>* 6= *nB*, *n<sup>B</sup>* 6= *nC*, and *n<sup>A</sup>* 6= *n<sup>C</sup>* (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).

#### *2.2. Expressions for the Ground State*

In the ground state (i.e., for *T* = 0), the grand canonical potential *ω*<sup>0</sup> per site of model (1) can be found as

$$
\omega\_0 = \langle \hat{H} \rangle / L = E\_D + E\_W + E\_{\mu \nu} \tag{3}
$$

where contributions associated with the onsite interaction, the intersite interactions, and the chemical potential, respectively, has the following forms

$$E\_D = \underbrace{\frac{U}{6}}\_{\text{(\\_}} [n\_A(n\_A - 1) + n\_B(n\_B - 1) + n\_C(n\_C - 1)]} \tag{4}$$

$$E\_W = -\frac{1}{6}\mathcal{W}\_1(n\_A n\_B + n\_B n\_\mathbb{C} + n\_\mathbb{C} n\_A) + \frac{1}{6}\mathcal{W}\_2(n\_A^2 + n\_B^2 + n\_\mathbb{C}^2),\tag{5}$$

$$E\_{\mu} \quad = \quad -\frac{1}{3}\mu(n\_A + n\_B + n\_C). \tag{6}$$

In the above expressions, concentrations *n<sup>α</sup>* 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 *ω*<sup>0</sup> 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

$$f\_0 = \langle \hat{H} + \mu \sum\_i \hbar\_i \rangle / L = \mathsf{U} \mathsf{D}\_{\mathsf{occ}} + \mathsf{E}\_{\mathsf{W}} \tag{7}$$

where *E<sup>W</sup>* is expressed by (5). *Docc* = (1/*L*) ∑*<sup>i</sup>* h*n*ˆ*i*↑*n*ˆ*i*<sup>↓</sup> i denotes the double occupancy and this quantity is found to be exact, cf. Table 2. One should underline that above expression for *f*<sup>0</sup> is an approximate result for model (1) on the triangular lattice. Here, it is assumed

that concentration *<sup>n</sup><sup>α</sup>* are as defined in Table <sup>2</sup> and they are the same in each <sup>√</sup> 3 × √ 3 unit cell in the system. Formally, it could be treated as exact one only if the numbers *z<sup>n</sup>* (*n* = 1, 2) goes to infinity.

The expressions presented in this subsection (for *W*<sup>2</sup> = 0) can be obtained as the *T* → 0 limit of the equations for *T* > 0 included in Section 2.3.

#### *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*<sup>1</sup> 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

$$
\omega = -\frac{1}{6} \sum\_{\mathfrak{a}} \Phi\_{\mathfrak{a}} n\_{\mathfrak{a}} - \frac{1}{3\mathfrak{P}} \sum\_{\mathfrak{a}} (\ln Z\_{\mathfrak{a}}).\tag{8}
$$

where *β* = 1/(*kBT*) is inverted temperature, coefficients Φ*<sup>α</sup>* are defined as Φ*<sup>α</sup>* = *µ* − *µα*,

$$Z\_{\mathfrak{u}} = 1 + 2 \exp(\beta \mu\_a) + \exp\left[\beta (2\mu\_a - \mathcal{U})\right],\tag{9}$$

and *µ<sup>α</sup>* is a local chemical potential in *α* sublattice (*α* ∈ {*A*, *B*, *C*})

$$
\mu\_A = \mu - \frac{1}{2} \mathcal{W}\_1 (n\_B + n\_\subset), \; \mu\_B = \mu - \frac{1}{2} \mathcal{W}\_1 (n\_A + n\_\subset), \; \mu\_\subset = \mu - \frac{1}{2} \mathcal{W}\_1 (n\_A + n\_B). \tag{10}
$$

For electron concentration *n<sup>α</sup>* in each sublattice in arbitrary temperature *T* > 0 one gets

$$n\_{\mathfrak{a}} = \frac{2}{Z\_{\mathfrak{a}}} \{ \exp \left( \beta \mu\_{\mathfrak{a}} \right) + \exp \left[ \beta (2 \mu\_{\mathfrak{a}} - \mathsf{U}) \right] \} \quad \text{(for } \mathfrak{a} \in \{ A, B, \mathbb{C} \} \text{)}.\tag{11}$$

The set of three Equations (11) for *nA*, *nB*, and *n<sup>C</sup>* determines the (homogeneous) phase occurring in the system for fixed model parameters *U*, *W*1, and *µ*. If *n* = (1/3)(*n<sup>A</sup>* + *n<sup>B</sup>* + *nC*) is fixed, one has also set of three equations, but it is solved with respect to *µ*, *nA*, and *n<sup>B</sup>* (the third *n<sup>α</sup>* is obviously found as *n<sup>C</sup>* = 3*n* − *n<sup>A</sup>* − *nB*).

The free energy *f* per site is derived as

$$f = \omega + \frac{1}{3}\mu (n\_A + n\_B + n\_{\mathbb{C}}) \,. \tag{12}$$

where *ω* and *nα*'s are expressed by (8)–(11).

#### *2.4. Macroscopic Phase Separation*

The free energy *fPS* of the (macroscopic) phase separated state (as a function of total electron concentration *n*; and at any temperature *T* ≥ 0) is calculated from

$$f\_{PS}(n) = \frac{n - n\_{-}}{n\_{+} - n\_{-}} f\_{+}(n\_{+}) + \frac{n\_{+} - n\_{-}}{n\_{+} - n\_{-}} f\_{-}(n\_{-}),\tag{13}$$

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*<sup>−</sup> < *n* < *n*+. For *n*<sup>±</sup> only the homogeneous phase exists in the system (one homogeneous phase occupies the whole system). Concentrations *n*<sup>±</sup> 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 *fPS* [i.e., (13)] with respect to *n*<sup>+</sup> and *n*<sup>−</sup> (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<sup>n</sup>* <sup>→</sup> <sup>∞</sup>, *<sup>n</sup>* <sup>=</sup> 1, 2) or <sup>√</sup> 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<sup>α</sup>* → 2 − *nα*; the NO<sup>1</sup> and TCO phases are invariant under this transformation). In the brackets also an alternative name is given. The degeneration *d<sup>c</sup>* × *d<sup>s</sup>* of the unit cells (equal to the degeneration of the ground state for *z<sup>n</sup>* → ∞ limit) and degeneration *D<sup>c</sup>* × *D<sup>s</sup>* 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 *Docc*. *n<sup>s</sup>* and *n<sup>f</sup>* define the range [*ns*, *n<sup>f</sup>* ] of *n*, where the phase is correctly defined. In the last column, the phase separated state degenerated with the homogeneous phase in range (*ns*, *n<sup>f</sup>* ) for *W*<sup>2</sup> = 0 is mentioned. Star "∗" in superscript indicates that the phase is obtained by the particle-hole transformation (i.e., *n<sup>α</sup>* → 2 − *nα*; TCOA, TCO∗ A , TCOB, and TCO∗ B phases are invariant under this transformation).


### **3. Results for the Ground State (***W***<sup>1</sup>** > **0 and** *W***<sup>2</sup>** ≤ **0)**

*3.1. Analysis for Fixed Chemical Potential µ*

The ground state diagram for model (1) as a function of (shifted) chemical potential *µ*¯ = *µ* − *W*<sup>1</sup> − *W*<sup>2</sup> 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: NO0, NO<sup>1</sup> and NO2), the DCO phase occurs (6 regions: DCO1, DCO2, DCO3, DCO<sup>∗</sup> 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: DCO2–DCO<sup>3</sup> (DCO<sup>∗</sup> 2 –DCO∗ 3 ) and TCO–NO1. 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*<sup>2</sup> = 0, only the boundaries DCO2–DCO<sup>3</sup> (DCO<sup>∗</sup> 2 –DCO∗ 3 ) and TCO–NO<sup>1</sup> have finite degeneracy (6 and 7, respectively, modulo spin degrees of freedom) and the interface between different types of <sup>√</sup> 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 *<sup>L</sup>*/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<sup>2</sup> 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<sup>2</sup> phase. Thus, the degeneracy of the DCO2–DCO<sup>∗</sup> 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*<sup>2</sup> < 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.

**Figure 2.** Ground state phase diagram of the model on the triangular lattice as a function of shifted chemical potential *µ*¯ = *µ* − *W*<sup>1</sup> − *W*<sup>2</sup> for *W*<sup>1</sup> > 0 and *W*<sup>2</sup> ≤ 0 (|*k*| = |*W*2|/*W*<sup>1</sup> ). The regions are labeled by the names of the phases defined in Table 1 and numbers corresponding to concentrations in each sublattice *nA*, *n<sup>B</sup>* and *nC*.

Please also note that for *W*<sup>2</sup> = 0 as well as for *W*<sup>2</sup> < 0 inside the regions shown in Figure 2, the <sup>√</sup> 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).

#### *3.2. 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*<sup>2</sup> = 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*<sup>2</sup> = 0, the TCO phase can exist only in the range of 0 < *U*/*W*<sup>1</sup> < 1/2 at *T* 6= 0. For *W*<sup>2</sup> < 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*<sup>2</sup> ≤ 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<sup>A</sup> and DCO<sup>B</sup> phases, (ii) the DCO<sup>A</sup> and DCO<sup>C</sup> phases, and (iii) the TCO<sup>A</sup> and TCO<sup>B</sup> phases are located at *U*/*W*<sup>1</sup> = 0 (and these corresponding for *n* > 1; the dashed line in Figure 3). For *W*<sup>2</sup> < 0 these lines do not overlap with the boundaries between corresponding phase separated states at *U*/*W*<sup>1</sup> − |*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*<sup>2</sup> < 0 are unstable (i.e., *∂µ*/*∂n* < 0) inside the regions of Figure 1. For *W*<sup>2</sup> < 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 HCOA–LCO<sup>A</sup> and HCOA–HCO<sup>B</sup> boundaries at *U*/*W*<sup>1</sup> = 0 with PS1A–PS1<sup>B</sup> and PS1A–PS1<sup>B</sup> boundaries at *U*/|*W*2| = 1, respectively. In [68] the boundaries between homogeneous phases for *W*<sup>2</sup> < 0 are not shown in Figure 3. Only boundaries between corresponding phase separated states are correctly presented in that figure for *W*<sup>2</sup> < 0. For *U*/*W*<sup>1</sup> > 0, the CBO<sup>A</sup> phase (corresponding to the HCO<sup>A</sup> phase from [59]) is not the phase with the lowest energy (among homogeneous phases) in any range of *n* (but for *U*/*W*<sup>1</sup> < 0 it has the lowest energy among all homogeneous states). However, the corresponding phase separated state NO0/CBO<sup>2</sup> (i.e., PS1<sup>A</sup> from [59]) can occur for *U*/*W*<sup>1</sup> > 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 *Docc*, but the chemical potential *µ* (calculated as *µ* = *∂ f* /*∂n*) changes discontinuously. Boundaries DCOA–DCOB, DCOA–DCOC, and TCOA–TCO<sup>B</sup> (and other transitions for fixed *n* at *U*/*W*<sup>1</sup> = 0) between homogeneous phases are associated with discontinuous change of only *Docc*. One should note that it is similar to transition between two checker-board ordered phases on the square lattice, namely CBOA–CBO<sup>B</sup> and CBOA–CBO<sup>C</sup> boundaries, cf. [68] (or the HCOA–LCO<sup>A</sup> and HCOA–HCO<sup>B</sup> boundaries, respectively, from [59]). At the other horizontal boundaries (i.e., transitions for fixed *n* at *U*/*W*<sup>1</sup> − |*k*| = 1/2 in Figure 3) two of *nα*'s and *Docc* change discontinuously. At commensurate fillings transitions with changing *U*/*W*<sup>1</sup> occur only at points indicated by squares in Figure 3.

All horizontal boundaries between phase separated states (which are stable for *W*<sup>2</sup> < 0) are connected with discontinuous changes of *Docc*. These boundaries located at *U*/*W*<sup>1</sup> − |*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<sup>B</sup> phase and the DCO<sup>C</sup> phase at *n* = 1/3 reduce to DCO<sup>1</sup> phase.

**Figure 3.** Ground state phase diagram of the model as a function of particle concentration *n* for *W*<sup>1</sup> > 0 and *W*<sup>2</sup> ≤ 0 (|*k*| = |*W*2|/*W*<sup>1</sup> ). The regions are labeled by the names of the homogeneous phases (cf. Table 2). For *W*<sup>2</sup> = 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*<sup>2</sup> = 0. For *W*<sup>2</sup> < 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*<sup>1</sup> = 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*<sup>2</sup> < 0. Squares denote transitions for fixed *n* between homogeneous phase at commensurate fillings.

#### **4. Results for Finite Temperatures (***W***<sup>1</sup>** > **0 and** *W***<sup>2</sup>** = **0)**

One can distinguish four ranges of *U* interaction, where the system exhibits qualitatively different behavior, namely: (i) *U*/*W*<sup>1</sup> < 0, (ii) 0 < *U*/*W*<sup>1</sup> < (1/3)ln(2), (iii) (1/3)ln(2) < *U*/*W*<sup>1</sup> < 1/2, and (iv) *U*/*W*<sup>1</sup> > 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.

**Figure 4.** Phase diagrams of the model for *U*/*W*<sup>1</sup> = −1.00 as a function of (**a**) chemical potential *µ*¯/*W*<sup>1</sup> and (**b**) particle concentration *n* (*W*<sup>1</sup> > 0, *W*<sup>2</sup> = 0). All transitions are first order and regions of phase separated state (PS<sup>1</sup> :NO/DCO and PS2:DCO/DCO) occurrence are present on panel (**b**). NO and DCO denote homogeneous phases defined in Figure 1b.

**Figure 5.** Phase diagrams of the model for *U*/*W*<sup>1</sup> = 0.75 as a function of (**a**) chemical potential *µ*¯/*W*<sup>1</sup> and (**b**) particle concentration *n* (*W*<sup>1</sup> > 0, *W*<sup>2</sup> = 0). All transitions are first order. Other denotations as in Figure 4.

For *U*/*W*<sup>1</sup> < 0 and *U*/*W*<sup>1</sup> > 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 PS1:NO/DCO occurs in define ranges of *n*. In this state domains of the NO and the DCO phases coexist.

For *U*/*W*<sup>1</sup> < 0 the temperature of NO–DCO transition is maximal for *µ*¯ = 0 (i.e., at half-filling)–Figure 4a. Its maximal value *T<sup>M</sup>* monotonously decreases with increasing of *U* from *kBTM*/*W*<sup>1</sup> = 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<sup>M</sup>* 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 *TM*, but two equivalent solutions still exist for any *T* < *T<sup>M</sup>* (similarly as in the ferromagnetic Ising model at zero field [9]). At *µ*¯ = 0 and *T* < *T<sup>M</sup>* the discontinuous transition between two DCO phases occurs. In the DCO phase for *µ*¯ < 0 (*n* < 1) [connecting with the DCO<sup>1</sup> (DCOA) region at *T* = 0] the relation *n<sup>A</sup>* = *n<sup>B</sup>* < *n<sup>C</sup>* 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<sup>A</sup>* < *n<sup>B</sup>* = *n<sup>C</sup>* occurs (*n<sup>C</sup>* 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*<sup>2</sup> = 0 the PS1:NO/DCO phase separated state occurs only for *T* > 0. For *T* → 0 the concentrations in both domains of the PS<sup>1</sup> state approach 0 (or 2), whereas for *T* → *T<sup>M</sup>* they approach to 1. Near *n* = 1 the PS2:DCO/DCO state is stable for 0 ≤ *T* < *TM*. 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*<sup>1</sup> = 0.20 as a function of (**a**) chemical potential *µ*¯/*W*<sup>1</sup> and (**b**) particle concentration *n* (*W*<sup>1</sup> > 0, *W*<sup>2</sup> = 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*<sup>1</sup> = 0.35 as a function of (**a**) chemical potential *µ*¯/*W*<sup>1</sup> and (**b**) particle concentration *n* (*W*<sup>1</sup> > 0, *W*<sup>2</sup> = 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*<sup>1</sup> > 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*<sup>1</sup> < 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 *kBTM*/*W*<sup>1</sup> of NO–DCO transition occurs for *µ*¯/*W*<sup>1</sup> 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

1

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*<sup>1</sup> < 1/2, where the TCO phase appears at *T* = 0 and for finite temperatures near half-filling. For 0 < *U*/*W*<sup>1</sup> < (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*<sup>1</sup> = 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*<sup>1</sup> 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*<sup>1</sup> < 0, the temperature of NO– DCO transition is maximal at half-filling. For fixed *n*, the narrow regions of PS1: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 DCO1–DCO<sup>2</sup> (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–61]). From the same reason, there is no boundary at *T* > 0 for fixed *n* associated to the DCOB–DCO<sup>C</sup> (DCO∗ B –DCO∗ C ) line occurring at *T* = 0 (Figure 6b). However, strong reduction of one *n<sup>α</sup>* from the case where *n<sup>α</sup>* ≈ 2 to the case of *n<sup>α</sup>* ≈ 1 is visible (some kind of a smooth crossover inside the DCO region), cf. Figure 8f–h for *U*/*W*<sup>1</sup> = 0.35.

**Figure 8.** Dependencies of particle concentrations *nα*'s in the sublattices (red dotted, blue dashed, and green dot-dashed lines) as a function of *kBT*/*W*<sup>1</sup> [(**a**–**f**)] and *µ*¯/*W*<sup>1</sup> [(**g**–**h**)] for *W*<sup>2</sup> = 0. Black solid lines denote total particle concentration *n* = (*n<sup>A</sup>* + *n<sup>B</sup>* + *nC*)/3. They are obtained for: (**a**) *U*/*W*<sup>1</sup> = −1.00, *µ*¯/*W*<sup>1</sup> = 0; (**b**) *U*/*W*<sup>1</sup> = 0.20, *µ*¯/*W*<sup>1</sup> = 0; (**c**) *U*/*W*<sup>1</sup> = 0.35, *µ*¯/*W*<sup>1</sup> = 0; (**d**) *U*/*W*<sup>1</sup> = 0.35, *µ*¯/*W*<sup>1</sup> = −0.15; (**e**) *U*/*W*<sup>1</sup> = 0.35, *µ*¯/*W*<sup>1</sup> = −0.5; (**f**) *U*/*W*<sup>1</sup> = 0.35, *µ*¯/*W*<sup>1</sup> = −0.8; (**g**) *U*/*W*<sup>1</sup> = 0.35, *kBT*/*W*<sup>1</sup> = 0.40; (**h**) *U*/*W*<sup>1</sup> = 0.35, *kBT*/*W*<sup>1</sup> = 0.80; (**i**) *U*/*W*<sup>1</sup> = 0.35, *kBT*/*W*<sup>1</sup> = 0.11. Vertical solid and dashed lines indicate points of continuous and discontinuous transitions, respectively.

1

For (1/3)ln(2) < *U*/*W*<sup>1</sup> < 1/2, the maximum of the NO–DCO transition temperature is shifted towards larger |*µ*¯|/*W*<sup>1</sup> (or smaller |1 − *n*|). This is associated with forming of the two-lob structure of the diagram found for *U*/*W*<sup>1</sup> > 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<sup>A</sup>* < *n<sup>B</sup>* = *n<sup>C</sup>* (for *µ*¯ < 0 or *n* < 1); cf. Figure 8e,i] are connected with the DCO<sup>3</sup> and DCO∗ 3 regions occurring at the ground state. The boundaries DCO–DCO weakly dependent on *µ*¯ are associated with occurrence of phase separated PS2:DCO/DCO states (at high temperatures) in some ranges of *n*, cf. Figure 7b. The other DCO–DCO transitions (which are almost temperature-independent) are not connected with phase separated states. Also the first-order TCO–NO line is present near half-filling, cf. Figure 8d. One should underline that all four lines (three first-order boundaries: DCO–NO, DCO–DCO, TCO–NO and the second-order TCO–DCO boundary) merge at single point with numeric accuracy. However, it cannot be excluded that the DCO–NO and TCO–DCO boundaries connect with the temperature-independent line in slightly different points, what could result in, e.g., the TCO–DCO–NO sequence of transition with increasing temperature for small range of chemical potential *µ*¯. All of these almost temperature-independent boundaries (i.e., the DCO–DCO and the TCO–NO lines) are located at temperature, which decreases with increasing *U*/*W*<sup>1</sup> and approaches 0 at *U*/*W*<sup>1</sup> = 1/2 [i.e., they connect with the DCO2–DCO<sup>3</sup> (DCO∗ 2 –DCO∗ 3 ) and TCO–NO<sup>1</sup> boundaries at *T* = 0 for fixed *µ* or with the TCO–DCO<sup>D</sup> (TCO∗–DCO∗ D ) lines at *T* = 0 for fixed *n*]. 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 *kBT*/*W*<sup>1</sup> = 1/6 and *U*/*W*<sup>1</sup> = (1/3)ln(2).

**Figure 9.** Phase diagram of the model for half-filling (*µ*¯ = 0 or *n* = 1) as a function of onsite interaction *U*/*W*<sup>1</sup> (*W*<sup>1</sup> > 0 and *W*<sup>2</sup> = 0). The order-disorder boundary for *U*/*W*<sup>1</sup> < (1/3)ln(2) is a line consisting of some higher-order critical points as discussed in the text.

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*<sup>1</sup> < 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 PS2:DCO/DCO state is stable at *n* = 1), cf. also Figure 4. For 0 < *U*/*W*<sup>1</sup> < 1/2 the TCO phase is stable below the order-disorder line, but for (1/3)ln(2) < *U*/*W*<sup>1</sup> < 1/2 and *kBT*/*W*<sup>1</sup> < 1/6 the TCO– NO phase transition is discontinuous (cf. also Figure 7). For *U*/*W*<sup>1</sup> < (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*<sup>1</sup> < (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*<sup>1</sup> < (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.

#### **5. 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*<sup>1</sup> > 0). The effects of attractive next-nearest-neighbor interaction (*W*<sup>2</sup> < 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*<sup>2</sup> < 0 or finite *T* > 0 removes the degeneration between homogeneous phases and phase separated states occurring at *T* = 0 for *W*<sup>2</sup> = 0. It was shown that TCO phase is stable in intermediate range of onsite repulsion 0 < *U*/*W*<sup>1</sup> < 1/2 (for *W*<sup>2</sup> = 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*<sup>1</sup> < (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.

Additionally, one notes that the mean-field results for model (1) with attractive *W*<sup>1</sup> < 0 and *W*<sup>2</sup> ≤ 0 are the same for both two-sublattice and tri-sublattice assumptions. In such a case, three different nonordered phases exist with the discontinuous first-order transition between them (at *µ*¯ = 0 for *U* < 0 or for |*µ*¯| 6= 0 for *U*/(|*W*1| + |*W*2|) > 1), and thus for fixed *n*, several so-called electron-droplet states (phase separations NO/NO) exist (cf. [60,68,87,88], particularly Figure 2 of [60]).

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 *kBTc*/|*J*| as 0.506, 0.568, 0.607, respectively, whereas the mean-field approximation gives *kBTc*/|*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–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.

**Author Contributions:** Conceptualization, K.J.K.; methodology, K.J.K.; software, K.J.K.; validation, K.J.K.; formal analysis, K.J.K.; investigation, K.J.K.; resources, K.J.K.; data curation, K.J.K.; writing—original draft preparation, K.J.K.; writing—review and editing, K.J.K.; visualization, K.J.K.; supervision, K.J.K.; project administration, K.J.K.; funding acquisition, K.J.K. The author has read and agreed to the published version of the manuscript.

**Funding:** The support from the National Science Centre (NCN, Poland) under Grant SONATINA 1 no. UMO-2017/24/C/ST3/00276 is acknowledged. Founding in the frame of a scholarship of the Minister of Science and Higher Education (Poland) for outstanding young scientists (2019 edition, no. 821/STYP/14/2019) is also appreciated.

**Institutional Review Board Statement:** Not applicable. The study does not involve humans or animals.

**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´nski, R. Lema´nski, 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:


#### **Appendix A. Analytic Expressions for Continuous Transition Temperatures**

Equations (11) can be written in a different form, namely *n<sup>α</sup>* = *fα*, where *f<sup>α</sup>* ≡ 2(*t<sup>α</sup>* + *t* 2 *α a*)/(1 + 2*t<sup>α</sup>* + *t* 2 *α a*), *t<sup>α</sup>* ≡ exp(*βµα*) and *a* ≡ exp(−*βU*). One can define ∆ ≡ (*n<sup>A</sup>* − *nB*)/2 and *χ* ≡ (*n<sup>B</sup>* − *nC*)/2. From (10) one gets:

$$
\mu\_A \quad = \quad \mu - W\_1 \Big[ n - \frac{1}{3} (2\Delta + \chi) \Big] \,\tag{A1}
$$

$$
\mu\_B \quad = \quad \mu - \mathcal{W}\_1 \Big[ n + \frac{1}{3} (\Delta - \chi) \Big], \tag{A2}
$$

$$
\mu\_{\mathbb{C}} \quad = \quad \mu - \mathcal{W}\_1 \Big[ n + \frac{1}{3} (\Delta + 2\chi) \Big]. \tag{A3}
$$

Taking the limit *χ* → 0 of both sides of the equation (*f<sup>B</sup>* − *fC*)/(2*χ*) = 1 (using de l'Hospital theorem) one gets (*g<sup>B</sup>* − *gC*)/2 = 1, where *g<sup>α</sup>* ≡ *∂ fα ∂χ* = *∂ fα ∂tα ∂tα ∂µα ∂µα ∂χ* . One easily finds that *∂ fα*/*∂t<sup>α</sup>* = 2(1 + 2*tαa* + *t* 2 *α a*)/(1 + 2*t<sup>α</sup>* + *t* 2 *α a*) 2 , *∂tα*/*∂µ<sup>α</sup>* = *βt<sup>α</sup>* as well as *∂µA*/*∂χ* = −*W*1/3, *∂µB*/*∂χ* = *W*1/3, *∂µC*/*∂χ* = −2*W*1/3. Finally, the equation determining temperature *T<sup>c</sup>* of a continuous transition (at which *n<sup>B</sup>* → *nC*) has the form

$$\frac{1}{\beta\_{\text{c}}W\_{1}} = \frac{\left(1 + 2t\_{\text{BC}}\overline{a} + t\_{\text{BC}}^{2}\overline{a}\right)t\_{\text{BC}}}{\left(1 + 2t\_{\text{BC}} + t\_{\text{BC}}^{2}\overline{a}\right)^{2}},\tag{A4}$$

where *tBC* ≡ exp(*βcµBC*), *µBC* = *µ* − *W*1(*n* + ∆/3) (in the considered limit *µ<sup>B</sup>* = *µ<sup>C</sup>* and *n<sup>B</sup>* = *nC*), *a*¯ ≡ exp(−*βcU*), *β<sup>c</sup>* ≡ 1/(*kBTc*). Concentrations *n<sup>A</sup>* and *nBC* ≡ *n<sup>B</sup>* = *n<sup>C</sup>* are calculated from (11) for *β<sup>c</sup>* self-consistently. Thus, for fixed *µ* (or *n*) one has a set of three equation which is solved with respect to *βc*, *n* (or *µ*) and ∆.

The solutions of (A4) and (11) with ∆ 6= 0 (i.e., *n<sup>A</sup>* 6= *nB*) correspond to the TCO–DCO boundaries. Such determined temperatures coincide with those found from the analysis of (11) and (8) or (12) and presented in Figures 6 and 7, what supports the findings that the TCO–DCO boundaries are indeed continuous.

The solutions of (A4) and (11) with ∆ = 0 (i.e., *n<sup>A</sup>* = *nB*) correspond to the continuous DCO–NO boundaries. On the diagrams presented in Section 4 such solutions for *T<sup>c</sup>* are located inside the regions of the DCO phase occurrence (and they correspond to the transitions between metastable phases [74] or to vanishing of the NO metastable solution, cf. [88,99]). In the present case of model (1) studied, they coincide with the DCO–NO transitions presented in Figures 4–7 only at *T* = 0 (i.e., for *n* = 0, 2 as well as for *n* = 1 and *U*/*W*<sup>1</sup> > 1/2; or corresponding *µ*¯) and at *T* = *T<sup>M</sup>* (i.e., maximal temperature of the DCO–NO transition, occurring for *U*/*W*<sup>1</sup> < (1/3)ln(2) and *n* = 1 or *µ*¯ = 0, as well as for *U*/*W*<sup>1</sup> > 1/2 and *n* ≈ 1/2, 3/2 or corresponding *µ*¯; for (1/3)ln(2) < *U*/*W*<sup>1</sup> < 1/2 it is located for some intermediate concentrations 1/2 < *n* < 1 and 1 < *n* < 3/2). For ∆ = 0, (A4) and (11) give the following results: (i) for *U* → −∞: *kBTc*/*W*<sup>1</sup> = *n*(2 − *n*)/2; (ii) for *U* = 0: *kBTc*/*W*<sup>1</sup> = *n*(2 − *n*)/4; and (iii) for *U* → +∞: *kBTc*/*W*<sup>1</sup> = *n*(1 − *n*)/2 (if *n* < 1) and *kBTc*/*W*<sup>1</sup> = (2 − *n*)(*n* − 1)/2 (if *n* > 1). Please note that such determined *T<sup>c</sup>* for ∆ = 0 is two times smaller than corresponding continuous transitions for the model considered on the hypercubic lattice within the mean-field aprroximation for the intersite term (for the same *U*/*W*<sup>1</sup> and *n*) [10,59,60].

#### **References**


#### **Short Biography of Author**

**Konrad Jerzy Kapcia** is an assistant professor (adjunct) at Condensed Matter Theory Division of Faculty of Physics of Adam Mickiewicz University (AMU) in Pozna´n (Poland). He graduated in physics (M.Sc.) in 2009 and defended his dissertation (PhD) in theoretical physics in 2014 (both at AMU). Next, he worked in the Institute of Physics of the Polish Academy of Sciences (PAS) in Warsaw and in the Henryk Niewodnicza ´nski Institute of Nuclear Physics of PAS in Kraków. He was a visiting scientist at SISSA in Trieste (Italy) and CFEL-DESY in Hamburg (Germany). He is interested in theoretical physics, mainly in the field of strongly correlated systems (including cold gases on optical lattices) and superconductivity as well as in ab initio modeling of real materials (electron and phonon properties). He was the head and principal investigator of several projects financed by National Science Centre (Poland) and a recipient of a few scholarships of Ministry of Science and Higher Education (Poland). He is a member of the Polish Physical Society, the Institute of Physics (IoP, UK), and the American Physical Society.
