Charge-order on the triangular lattice: Effects of next-nearest-neighbor attraction in finite temperatures

The extended Hubbard model in the atomic limit, which is equivalent to lattice $S=1/2$ fermionic gas, is considered on the triangular lattice. The model includes onsite Hubbard $U$ interaction and both nearest-neighbor ($W_{1}$) and next-nearest-neighbor ($W_{2}$) density-density intersite interactions. The variational approach treating the $U$ term exactly and the $W_l$ terms in the mean-field approximation is used to investigate thermodynamics of the model and to find its finite temperature ($T>0$) phase diagrams (as a function of particle concentration) for $W_{1}>0$ and $W_{2}<0$. Two different types of charge-order (i.e., DCO and TCO phases) within $\sqrt{3} \times \sqrt{3}$ unit cells as well as the nonordered (NO) phase occur on the diagram. Moreover, several kinds of phase-separated (PS) states (NO/DCO, DCO/DCO, DCO/TCO, and TCO/TCO) are found to be stable for fixed concentration. Attractive $W_{2}<0$ stabilizes PS states at $T=0$ and it extends the regions of their occurrence at $T>0$. The evolution of the diagrams with increasing of $|W_{2}|/W_{1}$ is investigated. It is found that some of the PS states are stable only at $T>0$. Two different critical values of $|W_{2}|/W_{1}$ are determined for the PS states, in which two ordered phases of the same type (i.e., two domains of the DCO or TCO phase) coexist.


I. INTRODUCTION
The classical lattice gas model (equivalent with the S = 1/2 Ising model) is useful effective model for description of adsorbed particles on crystalline substrates (cf. pioneering works on the triangular lattice [1][2][3][4][5][6][7]). In the case of a graphine surface or a single layer of graphene as well as (111) face-centered cubic surface, the periodic potential of the underlying crystal surface forms a triangular lattice, which can be occupied by adsorbed atoms, e.g., [8][9][10][11][12]. Although the adsorbed particles are rather classical, taking into account the quantum properties is necessary for a description of helium atoms adsorption [13][14][15][16]. Thus, in the present work, an extension of the classical lattice gas model to S = 1/2 fermionic particles is analyzed. Particular attention is taken for effects of the next-nearest-neighbor attraction on the phase diagrams at finite temperatures.
i,j l denotes the summation over lth neighbors independently (l = 1, 2). U is the onsite interaction, whereas W 1 and W 2 denote the intersite interactions between the nearest neighbors (NNs) and the next-nearest neighbors (NNNs), respectively. z 1 and z 2 are numbers of NNs and NNNs, respectively (z 1 = z 2 = 6 for the triangular lattice). The chemical potential µ is related with the total concentration n of particles in the system through n = (1/L) i n i , where 0 ≤ n ≤ 2 and L is the total number of lattice sites.
Within the variational approach treating onsite U term exactly and intersite W l terms in the mean-field approximation, i.e., n inj = n i n j +n i n j − n i n j , arXiv:2111.02699v1 [cond-mat.str-el] 4 Nov 2021 this model was intensively investigated on the hypercubic lattices, e.g., for W 2 = 0 [21][22][23] and W 2 = 0 [24][25][26][27]. Moreover, rigorous results for one-dimensional chain were found in [28,29] as well as the model on two-dimensional square lattice was analyzed by various methods [27,[30][31][32][33][34]. On the triangular lattice [within the variational approach with decoupling (2)], the evolution of metastable phases in the model was determined for U < 0 and W 2 = 0 [35], whereas the full phase diagram (for all temperatures, at T = 0 and T > 0) for W 2 = 0 was obtained in [36]. In [36] the effects of W 2 < 0 only at the ground state were also discussed.
In the present work, the effects of the attractive NNN interaction (i.e., W 2 < 0) are investigated at T > 0 within the variational approach mentioned, for details cf. [22-25, 27, 35, 36]. In particular, the evolution of phase diagrams (for fixed U/W T , where W T = W 1 − 2W 2 , W 1 > 0, and W 2 ≤ 0) with increasing k = |W 2 |/W 1 is presented and an emergence of novel phase separation states (not occurring for W 2 = 0 or T = 0) is noticed. Because the triangular lattice can be divided into three equivalent sublattices only orderings within √ 3 × √ 3 unit cell are considered (the three-sublattice assumption), which is justified in the range of model parameters considered.
For T > 0, the expressions given in [24] for the triangular lattice and W 2 = 0 take the following forms (cf. also these in [27,36]). For a grand canonical potential ω (per lattice site) one obtains where β = 1/(k B T ) is inverted temperature, coefficients ϕ α are defined as ϕ α = µ − µ α , and µ α is a local chemical potential in α sublattice (α ∈ {A, B, C}) defined as Here, α and α denote two other sublattices than α (and α = α ). Particle concentration n α = (3/L) i∈α n i in each sublattice α for arbitrary T > 0 is expressed by Three equations for n α determine the solution for a (homogeneous) phase occurring in the system for fixed U , W 1 , W 2 , and µ. If n = (n A + n B + n C )/3 is fixed, the set is solved with respect to µ, n A , and n B (the third concentration is found as n C = 3n − n A − n B ). This set has usually several solution, thus it is extremely important to find the solution corresponding to the lowest ω (if µ is fixed) or free energy f = ω + µn (if n is fixed). For fixed n, the phase-separated (PS) states can also occur, which free energy is determined by the Maxwell's construction (macroscopic phase separation), e.g., [23,[36][37][38][39][40].  It was shown that in the model (within the approximation used) the following phases can occur: (i) the nonordered (NO) phase with n A = n B = n C , (ii) the charge-ordered phase with two different concentration in sublattices (the DCO phase, e.g., n A = n B = n C and other cyclic permutations, 3 equivalent solutions), and (iii) the charge-ordered phase with n A = n B , n B = n C , and n A = n C (the TCO phase; three different concentrations in sublattices, 6 equivalent solutions) [36,41]. Moreover, for W 2 = 0, two PS states were found in some ranges of n: (i) PS1:NO/DCO, where the NO and the DCO phases coexist and (ii) PS2:DCO/DCO, where two different DCO phases coexist [35,36]. Due to the particle-hole symmetry of model (1) the phase diagram is symmetric with respect to half-filling, i.e., n = 1 or [17,22,35,36].
The simplest phase diagram of the model is for U < 0. For W 2 = 0, the diagram for fixed µ consists of two regions of the DCO phase and one region of the NO phase. These regions are separated by two kinds of firstorder (discontinuous) boundaries: (i) the DCO-NO line with its maximum temperature T M at half-filling and (ii) the DCO-DCO line atμ = 0 extending from T = 0 to T = T M (it isμ-independent) [35,36] (cf. also Fig. 1). As a result, on the diagram as a function of n one finds regions of two PS states occurrence: (i) the PS1 state oc- curring in narrow range of n for 0 < T < T M (it does not exist at T = 0, Fig. 2) and (ii) the PS2 state, which is stable for 0 ≤ T < T M (with concentrations in coexisting domains as n − = 2/3 and n + = 4/3 at T = 0; Fig. 3), respectively.
Nonzero W 2 < 0 extends the regions of the PS states occurrence. For W 2 < 0, the PS1 state is stabilized at T = 0 with n − = 0 and n + = 2/3. For 0 < k < k c1 (where k c1 ≈ 3/20) and above some T * c1 (which is U/W Tand k-dependent), the PS2 state appears in a narrow region (shadowed regions in Fig. 2). It vanishes continuously at k = k c1 . With increasing T for fixed n (higher concentrations), at T * c1 the PS1-PS2 transition occurs, which is associated with the change of the phase in the domain of lower concentration (from the NO to the DCO phase) with simultaneous a discontinuous change of concentration n − in this domain. For lower concentrations n, there is also a transition between two different PS1 states (the PS1-PS1 transition). At T * c1 the discontinuous change of concentration in the DCO phase domain of the PS1 state occurs. This behavior is associated with a new first-order DCO-DCO boundary inside the DCO region [ending at a bicritical-end (also called as isolated-critical) point, cf. [22,24,25,37]], which is present on the diagram for fixedμ. This is schematically shown only in the inset of Fig. 1, where three first-order lines merge in the triple point located at T * c1 . T * c1 increases with k and for k → k c1 the bicritical-end point goes to the DCO-DCO line. On the k B T /W T -μ/W 1 diagrams the DCO region shrinks with increasing k and simultaneously the re-entrant feature of the DCO-DCO line is destroyed for large k. However, the evolution of the DCO-DCO boundary is nonmonotonous (these for k > k c1 can cross each other, approximately, in the range of −1.0 <μ/W 1 < −0.9, Fig. 1). For k > k c1 , the PS2 region is absent (the triple and the bicritical-end points do not exist) and only the PS1 state occurs on the diagram for n < 2/3 (the PS1 region separates the DCO and the NO regions for all T < T M ).
The evolution of the PS2 region boundaries in the range 2/3 < n < 4/3 is shown in Fig. 3. Increasing k enlarges the PS2 region existing for 0 ≤ T < T M with simultaneous change of boundary curvatures.
The above discussed behaviors are generic for all U < |W 2 | (or U/W 1 − k < 0). For U/W T > 1 the critical behaviors are similar to those discussed above, but the structure of phase diagrams exhibits two lobs of the DCO phase occurrence [36] and the maximal temperature T M for the DCO-DCO transition is located forμ corresponding to n = 1/2 and n = 3/2 (cf. also [22,24,25]).
For U > |W 2 | and U/W T < 1/2 (0 < U/W 1 −k < 1/2), the TCO phase appears on the phase diagram (Fig. 4). For k = 0, the TCO phase is stable near half-filling (in a range of 2/3 < n < 4/3 at T = 0) and the TCO-DCO transition is continuous (second order) [36], cf. also Appendix A. For k > 0 the following qualitative changes occur, which are shown in Fig. 4 for U/W T = 0.2. The evolution of the DCO-NO boundary for fixedμ is the same as for the case discussed previously. With increasing k (for 0 < k < k c1 ) the triple point appears [with the DCO-DCO boundary at T > T * c1 , not shown in Fig.  4(a); cf., inset of Fig. 1] and the PS1 and PS2 states are present in define ranges of n, whereas for k > k c1 only the PS1 state occurs (but now, at T = 0, the concentrations in domains of the PS1 state are n − = 0 and n + = 1/3). For smaller |μ| other discontinuous DCO-DCO bound-ary extending from T = 0 exists [almost straight lines in Fig. 4(a)], which results in the PS2 state occurrence in 1/3 < n < 2/3 range (this behavior is not present for k = 0).
The attractive W 2 = 0 affects also on the DCO-TCO boundary. It changes its order for small T and for T < T * c2 the DCO-TCO transition is first order (for fixed µ). For T < T * c2 the PS3:DCO/TCO state, which is a coexistence of the DCO and the TCO phases, is stable (at T = 0 in the range of 2/3 < n < 1). For 0 < k < k c2 (where k c2 ≈ 3/10) and T > T * c2 , the PS4:TCO/TCO state appears in a narrow region [ Fig. 4(b)], which shrinks with increasing k and vanishes continuously at k = k c2 . This behavior is connected with the occurrence of a first-order TCO-TCO transition at T > T *

c2
[not shown in Fig. 4(a)], which ends at a bicriticalend point [the discontinuous DCO-TCO, the continuous DCO-TCO, and discontinuous TCO-TCO lines merge at a critical-end point, schematically shown only in the inset of Fig. 4(a)]. Note that at T * c2 the PS3-PS4 transition occurs for fixed n, which is associated with the TCO-DCO transition in one domain. T * c2 increases with k and for k = k c2 the bicritical-end and the criticalend points merges into one critical point of higher order (cf. [37]). For k > k c2 , the discontinuous TCO-TCO line is no longer present on the diagram and a tricritical point appears, at which the DCO-TCO boundary changes its order. As a result, only the PS3 state occurs for 2/3 < n < 1 separating the DCO and the TCO regions below the tricritical point. The discussed behavior is similar to those occurring at the boundaries between checker-board charge-ordered phase and the NO phase for k H c = 3/5 for model (1) considered on the hypercubic lattices [24,25]. Note also that for U/W T = 0.2 and k > 1/3 one gets that U/|W 2 | < 1 and the TCO phase is no longer present on the diagram, which is the regime discussed at the beginning of this section (cf. Figs. 1-3). Thus, the tricritical point for the PS3 state exists for U/W T = 0.2 only in the range k c2 < k < 1/3.
One should note that, for (1/3) ln(2) < U/W T < 1/2, the phase diagram at high temperatures near half-filling is different than that shown in Fig. 4 (these regions are schematically indicated by gray rectangles). For example, for k = 0, the maximum of the DCO-NO transition (for fixedμ) is not located at half-filling, discontinuous DCO-DCO transitions appears at T > 0 (that results in new regions of the PS1 and PS2 states stability for fixed n), the direct discontinuous TCO-NO transitions is present in define range ofμ/W 1 and n [36]. The detailed analysis of these issues in the presence of W 2 = 0 is beyond the scope of this work. However, the main results for the evolution of the DCO-NO boundaries (with vanishing of the triple point at k = k c1 and associated PS1-PS2 transition) and of the TCO-DCO boundary (with a change of the critical-end point into the tricritical point at k = k c2 and associated PS3-PS4 transition) are still valid (cf. also [22,24,25]). shows schematically the structure of the diagram with critical-end and bicritical-end points for 0 < k < kc2, where a discontinuous TCO-TCO line appears at T > T * c2 . Squares, diamonds, circles, and triangles indicate triple, bicritical-end, criticalend, and higher-order critical points. Triple and critical-end points correspond to three or two different concentrations, respectively. Above critical-end or higher-order critical points, or for k = 0, the DCO-TCO transition is second order. All other boundaries are first order. Not all bicritical-end points, not all DCO-DCO lines, neither no TCO-TCO lines are shown on panel (a).

III. CONCLUSIONS
In this work, the extended Hubbard model in the atomic limit [Eq. (1) with W 1 > 0 and W 2 ≤ 0] on the triangular lattice was investigated. In particular, the effects of next-nearest-neighbor attraction W 2 were analyzed in detail. Increasing k = |W 2 |/W 1 > 0 affects the boundaries by increasing discontinuity of n at the transitions (for fixedμ) and extends the regions of phase separation occurrence for fixed n (at T = 0 they are stable for any incommensurate fillings, i.e., n = 2i/3, i = 1, 2, 3 for U/W 1 − k < 0 and n = j/3, j = 1, . . . , 6 for U/W 1 − k > 0). Two different critical values of k exist, below which at T > 0 the first-order transition between two phases of the same type occurs (and corresponding phase separated states for define ranges of concentrations are stable), namely: (i) k c1 ≈ 3/20 for the DCO-NO line and emerging DCO-DCO line and the PS2 region at T > T * c1 , (ii) k c2 ≈ 3/10 for the TCO-DCO boundaries and emerging TCO-TCO line and the PS4 region at T > T * c2 . For k > 0 also the PS2 state is stable for some T ≥ 0 inside the range of 1/3 < n < 2/3 (for U/W 2 > 1 and U/W T < 1/2).
Note that, in the U → −∞ limit, model (1) reduces to the well-known S = 1/2 Ising model. In the general case, model (1) can be mapped onto the S = 1 Blume-Capel model in the field with an effective temperaturedependent single-ion anisotropy [22,32].
Decoupling (2) of the intersite terms is exact only in z l → +∞ limit [42,43]. Thus, it is an approximation for the triangular lattice in the general case. It overestimates the critical temperatures and stability regions of ordered phases and could not properly describe the purely two-dimensional system investigated, where Berezinskii-Kosterlitz-Thouless-like state exists [41]. However, it is rigorous theory at the ground state for the model considered. The longer-range interaction or small interactions between two-dimensional layers can stabilize long-range order [6,44].
Hamiltonian (1) is a relatively simple toy model and it is oversimplified in many aspects for description of real materials. However, it can be treated as a benchmark for various approximate approaches for models with finite intersite hopping. One should underline here, that in the case of nonzero hopping, for U < 0 (and any W l ) or U > 0 and W l < 0, various superconducting states could appear and the stability ranges of the charge-ordered phases might be reduced (cf., e.g., the results for hypercubic lattices [17,[45][46][47][48][49]). Moreover, an occurrence of the Moiré pattern (which is the triangular lattice with a very large supercell) in the twisted-bilayer graphene (associated with emergence of superconductivity) [50,51] and hetero-bilayer transition metal dichalcogenides [52,53] makes further studies of various models on the triangular lattice worthwhile. Also ultra-cold atomic gases on the triangular lattice created by laser trapping [20,[54][55][56] are systems, which could enable testing of some theoretical predictions of this work.

Acknowledgments
The author thanks R. Micnas for very useful discussions on some issues raised in this work. The support from the National Science Centre (Poland) under Grant SONATINA 1 no. UMO-2017/24/C/ST3/00276 is ac-knowledged. Founding in the frame of a scholarship of the Minister of Science and Higher Education (Poland) for outstanding young scientists (no. 821/STYP/14/2019) is also appreciated.