1 Introduction↩︎

A 2D triangular lattice is typical for a number of organic conductors [1]–[3], transition-metal oxides and dichalcogenides, can be formed by adsorbed helium atoms on a surface, and can describe the moiré lattices [4]–[6]. The latter is an interesting platform to investigate various strongly-correlated and frustration-induced phenomena, since the interaction parameters and carrier itineracy can be controlled by changing a twist angle and the choice of two layers out of rich family of 2D materials [7]–[9]. Another useful platform to investigate the triangular lattice experimentally are the ultracold atomic gases on optical lattices [10]–[17]. The systems with triangular-lattice structure are found to host various phenomena, such as superconductivity [18]–[25], variety of charge and magnetic orderings [26]–[35], topological states [36], [37]. Among them, the charge ordering, e.g., (generalized) Wigner crystals [38]–[41], charge-transfer insulators, charge-density waves [42]–[45], or charge glasses, attracts researchers interest due to its interplay with superconductivity, as well as, possible applications in new devices, such as involving pyroelectric or ferroelectric materials [46]. Triangular lattices are commonly recognized to be perspective for searching the exotic charge orders, e.g., a pinball-liquid (PL) order [47]–[52]. The PL phase consists of the lattice sites that are insulating for the charge carriers (pins) but surrounded with the lattice sites where the charge carriers are itinerant (balls), and thus the PL shares some properties of supersolids [53].

The model for the investigation of charge-ordering phenomenon is the extended Hubbard model (EHM) [54]–[62]. It has been investigated with the mean-field approximation (MFA) both in the context of organic conductors for the quarter-filling or 3/4-filling [63], [64], and in the context of the moiré lattices for a series of a few fillings [40], [65], [66]. The methods beyond the MFA, such as the dynamical mean-field theory [67]–[69], are actively used to investigate the triangular-lattice materials with a specific charge concentration as well [40], [49]. Nevertheless, to understand the full picture of various charge orders in the triangular lattice, the investigation in the grand-canonical ensemble with arbitrary charge concentration is required. It is especially relevant since a particular charge concentration may appear in a phase-separation region: this can be captured as a first-order transition within the fixed-chemical-potential approach, while the fixed-concentration approach might yield an artificial result without consideration of phase separated states explicitly, cf. e.g., [70], [71]. Here we take advantage on the investigation of the EHM itself, encompassing both existing and theoretical materials with arbitrary repulsive electron interaction, in search of attractive properties and the mechanisms behind them. It is worth to highlight that a simple Hubbard model provided significant insights into condensed-matter physics and described various phenomena despite its extreme simplicity and artificiality, while the EHM is its most obvious continuation. Meanwhile, the successes in the ultracold atomic gases in optical lattices have also renewed the interest in Hubbard-like models, since the physics of such quantum simulators is not just modeled with these models, but described by them.

The mean-field study of the triangular-lattice EHM in the grand-canonical ensemble can provide significant insights into the problem despite ignoring the correlation effects. One should remember that the Hartree-Fock MFA overestimates critical temperatures and particularly the temperature range of stability of long-range orders, but it can, nevertheless, give a qualitative description of the system in the ground state in certain interaction parameter ranges [54], [70], [71]. Besides using it as a benchmark for further investigations, a number of unusual phenomena can already be found within the MFA, which makes them both easy to analyze and distinguish with strongly-correlated phenomena [72], [73]. The great advantage here (besides the required computational and time resources) is provided by the opportunity to analyze band structures in details, which cannot be done within dynamical mean-field theory that we consider to be the next step of this research (to take into account local correlations). Moreover, the non-correlated phases are found within the dynamical mean-field theory when the intersite interaction prevails over the onsite interaction [74], while some other phases and phase transitions can have similar qualitative features as withing the MFA. It is advantageous to use MFA results as a reference point in future studies together with the atomic limit results [75], [76].

Here, we present the solution of the triangular-lattice EHM within the MFA focusing on the zero-temperature systems without a magnetic order and utilizing the \(\sqrt{3}\times\sqrt{3}\) hexagonal supercells. The found phase diagram consists of the large variety of phase transitions and is more complex than one could assume in the absence of strong-correlation effects. Hence, the results included in this work are i) required as a starting point for the further investigations that would yield even more features in the phase diagram together with being more computationally expensive and suffering from convergence problems; and ii) useful as a tool to interpret the future results while being an extensive standalone research itself. Both the pinball-liquid phase and the (often ignored) strong particle-hole asymmetry are found and analyzed.

In the paper, after the discussion of the method of investigation (Sec. 2), we present the known band structures of the non-interacting triangular and honeycomb lattices, required for the further discussion (Sec. 3), and present the found mean-field phase diagram of the triangular lattice (Sec. 4) with a description of the found phases and phase transitions (Secs. 4.1-4.3). Sec. 4.5 is devoted to the brief comparison of the mean-field and the atomic-limit results. The main finding are summarized in Sec. 5.

2 Model and Method↩︎

The extended Hubbard model in a grand-canonical ensemble of electrons is represented by the Hamiltonian [54]–[62]: \[\begin{align} \label{ehm} \hat{H} & = & - t \sum_{\left< i,j \right>, \sigma} \left( \hat{c}^\dagger_{i\sigma} \hat{c}_{j\sigma} + \text{H.c.} \right) + U\sum_i \hat{n}_{i\uparrow}\hat{n}_{i\downarrow} \\ & + & V \sum_{\left< i,j \right>} \hat{n}_i \hat{n}_j - \mu\hat{N}, \nonumber \end{align}\tag{1}\] where \(t\), \(\mu\), \(U\), and \(V\) are a hopping amplitude, a chemical potential, an onsite and an intersite nearest-neighbor (NN) interaction parameters, respectively. These parameters are effective, meaning they can include not only Coulomb repulsion but other interactions and renormalizations such as those that involve phonons. The \(i\) and \(\sigma\) are site and spin indices while the summation over \(\left<i,j\right>\) means a summation over NN pairs without repeating (i.e., if the term with \(ij\) is present in a sum, the term with \(ji\) is not). The \(\hat{c}^\dagger_{i\sigma}\) and \(\hat{c}_{i\sigma}\) are creation and annihilation operators, \(\hat{n}_{i\sigma} = \hat{c}^\dagger_{i\sigma}\hat{c}_{i\sigma}\) is an occupation number operator, \(\hat{n}_i = \hat{n}_{i\uparrow} + \hat{n}_{i\downarrow}\), and \(\hat{N} = \sum_{i} \hat{n}_i\).

The mean-field approximation \[\hat{n}_{i\sigma}\hat{n}_{j\sigma'} = n_{j\sigma'} \hat{n}_{i\sigma} + n_{i\sigma} \hat{n}_{j\sigma'} - n_{i\sigma} n_{j\sigma'}\] (\(n_{i\sigma} = \left< \hat{n}_{i\sigma} \right>\)), and the Fourier transform to a reciprocal space turn the Hamiltonian into a sum of independent terms. In particular, for the triangular lattice and \(\sqrt{3}\times\sqrt{3}\) supercell \[\begin{align} \label{hamiltonian} \hat{H} &=& \sum_{\mathbf{k}\sigma} \left[ \varepsilon_\mathbf{k} (\hat{c}^\dagger_{1\mathbf{k}\sigma}\hat{c}_{2\mathbf{k}\sigma} + \hat{c}^\dagger_{3\mathbf{k}\sigma}\hat{c}_{1\mathbf{k}\sigma} + \hat{c}^\dagger_{2\mathbf{k}\sigma}\hat{c}_{3\mathbf{k}\sigma}) +\text{H.c.}+ \right. \nonumber\\ &+& \left. \sum_\alpha \epsilon_{\alpha\sigma} \hat{n}_{\alpha\mathbf{k}\sigma}\right]+C=\sum_{\mathbf{k}\sigma} \hat{H}_{\mathbf{k}\sigma} + C, \end{align}\tag{2}\] where \(\alpha = 1, 2, 3\) is a sublattice index, \(\mathbf{k}\) is a reciprocal-space vector, the constant term \[C = - \frac{L}{3} \sum_\alpha \left( U n_{\alpha\uparrow}n_{\alpha\downarrow} + \frac{zV}{2} n_{\bar\alpha} n_{\bar{\bar\alpha}} \right),\] (\(\bar\alpha\) and \(\bar{\bar\alpha}\) are sublattice indices different from \(\alpha\) and from each other), \(L\) is the number of lattice sites, \(n_\alpha = \frac{3}{L}\sum_{\sigma\mathbf{k}} \left< \hat{n}_{\alpha\mathbf{k}\sigma} \right>\), \(z=6\) is a coordination number, \[\label{en-sbl} \epsilon_{\alpha\sigma} = U n_{\alpha\bar\sigma} + \frac{zV}{2} (n_{\bar\alpha} + n_{\bar{\bar\alpha}}) -\mu,\tag{3}\] (\(\bar\sigma\) is the spin index different from \(\sigma\)) and \[\varepsilon_\mathbf{k} = - t \left( e^{i\mathbf{k}\mathbf{r}_1} + e^{i\mathbf{k}\mathbf{r}_2} + e^{i\mathbf{k}\mathbf{r}_3} \right)\] with: \[\mathbf{r}_1 = \left(-\frac{1}{3},-\frac{2}{3}\right), \; \mathbf{r}_2 = \left(-\frac{1}{3},\frac{1}{3}\right), \; \mathbf{r}_3 = \left(\frac{2}{3},\frac{1}{3}\right),\] provided that the vectors \(\mathbf{k}\) are written in the basis of a cell that is reciprocal to the \(\sqrt{3}\times\sqrt{3}\) supercell.

Each term \(\hat{H}_{\mathbf{k}\sigma}\) in the Hamiltonian (2 ) can be written as \[\hat{H}_{\mathbf{k}\sigma} = \hat{I}_{\mathbf{k}_1\uparrow} \otimes \hat{I}_{\mathbf{k}_1\downarrow} \otimes \hat{I}_{\mathbf{k}_2\uparrow} \otimes \ldots \otimes \hat{\mathbf{H}}_{\mathbf{k}\sigma} \otimes \ldots \otimes \hat{I}_{\mathbf{k}_{L/3}\downarrow},\] where \(\hat{\mathbf{H}}_{\mathbf{k}\sigma}\) can be represented by in a form of a block-diagonal \(8 \times 8\) matrix: \[\label{hammatrix} \hat{\mathbf{H}}_{\mathbf{k}\sigma} = \begin{pmatrix} 0 & \mathbb{0}_{1\times 3} & \mathbb{0}_{1\times3} & 0 \\ \mathbb{0}_{3\times1} & \check{\mathbb{H}}_{1\mathbf{k}\sigma} & \mathbb{0}_{3\times3} & \mathbb{0}_{3\times1}\\ \mathbb{0}_{3\times1} & \mathbb{0}_{3\times3} & \check{\mathbb{H}}_{2\mathbf{k}\sigma} & \mathbb{0}_{3\times1}\\ 0 & \mathbb{0}_{1\times3} & \mathbb{0}_{1\times3} & \epsilon_1 + \epsilon_2 + \epsilon_3 \end{pmatrix}_{\mathbf{k}\sigma},\tag{4}\] Here, \(\mathbb{0}_{m \times n}\) denotes a block of \(m \times n\) size with all elements \(0\), whereas \(\check{\mathbb{H}}_{1\mathbf{k}\sigma}\) and \(\check{\mathbb{H}}_{2\mathbf{k}\sigma}\) are the following \(3 \times 3\) matrices: \[\label{eq46bandstructure} \check{\mathbb{H}}_{1\mathbf{k}\sigma} \equiv \begin{pmatrix} \epsilon_{1\sigma} & \varepsilon_\mathbf{k} & \varepsilon_\mathbf{k}^*\\ \varepsilon_\mathbf{k}^* & \epsilon_{2\sigma} & \varepsilon_\mathbf{k}\\ \varepsilon_\mathbf{k} & \varepsilon_\mathbf{k}^* & \epsilon_{3\sigma} \end{pmatrix},\tag{5}\] \[\check{\mathbb{H}}_{2\mathbf{k}\sigma} \equiv \begin{pmatrix} \epsilon_{1\sigma}+\epsilon_{2\sigma} & \varepsilon_\mathbf{k} & -\varepsilon_\mathbf{k}^*\\ \varepsilon_\mathbf{k}^* & \epsilon_{1\sigma}+\epsilon_{3\sigma} & \varepsilon_\mathbf{k}\\ - \varepsilon_\mathbf{k} & \varepsilon_\mathbf{k}^* & \epsilon_{2\sigma}+\epsilon_{3\sigma} \end{pmatrix}.\] The operator matrix (4 ) is built from the matrix elements with respect to the basis states \[\begin{align} \left| 0_{1\mathbf{k}\sigma} 0_{2\mathbf{k}\sigma} 0_{3\mathbf{k}\sigma} \right>, \\\left| 1_{1\mathbf{k}\sigma} 0_{2\mathbf{k}\sigma} 0_{3\mathbf{k}\sigma} \right>, && \left| 0_{1\mathbf{k}\sigma} 1_{2\mathbf{k}\sigma} 0_{3\mathbf{k}\sigma} \right>, && \left| 0_{1\mathbf{k}\sigma} 0_{2\mathbf{k}\sigma} 1_{3\mathbf{k}\sigma} \right>,\\ \begin{align} \left| 1_{1\mathbf{k}\sigma} 1_{2\mathbf{k}\sigma} 0_{3\mathbf{k}\sigma} \right>, && \left| 1_{1\mathbf{k}\sigma} 0_{2\mathbf{k}\sigma} 1_{3\mathbf{k}\sigma} \right>, && \left| 0_{1\mathbf{k}\sigma} 1_{2\mathbf{k}\sigma} 1_{3\mathbf{k}\sigma} \right>, \end{align} \\ \left| 1_{1\mathbf{k}\sigma} 1_{2\mathbf{k}\sigma} 1_{3\mathbf{k}\sigma} \right>, \end{align}\] and the subscript \(\mathbf{k}\sigma\) under the matrix (4 ) means that it is associated with only \(\mathbf{k}\sigma\) single-particle states.

Thus, the model can be easily solved on a fine grid of \(\mathbf{k}\)-vectors when the parameters \(V/t\), \(U/t\), \(\mu/t\), and \(n_{\alpha\sigma}\) are provided. To distinguish stable and metastable phases the grand potential is calculated as \[\label{eq46omega} \frac{\Omega}{L} = \frac{1}{L} \left( C - \frac{1}{\beta}\ln\mathcal{Z} \right) = \frac{1}{L} \left( C - \frac{1}{\beta} \sum_{\mathbf{k}\sigma} \ln\mathcal{Z}_{\mathbf{k}\sigma} \right),\tag{6}\] where \(\mathcal{Z} = \prod_{\mathbf{k}\sigma}\mathcal{Z}_{\mathbf{k}\sigma}\) is the partition function of the system, \(\mathcal{Z}_{\mathbf{k}\sigma}=\sum_m e^{-\beta E^{(m)}_{\mathbf{k}\sigma}}\), \(E^{(m)}_{\mathbf{k}\sigma}\) is the \(m\)-th eigenvalue of \(\hat{\mathbf{H}}_{\mathbf{k}\sigma}\) (i.e., of the matrix (4 )), and \(\beta = T^{-1}\) is the inverse temperature.

The solution of the model gives the occupation numbers: \[\label{eq46occup} n_{\alpha\sigma} = \frac{3}{L} \sum_\mathbf{k} \frac{\sum_m n_{\alpha\mathbf{k}\sigma}^{(m)} e^{-\beta E^{(m)}_{\mathbf{k}\sigma}}}{\mathcal{Z}_{\mathbf{k}\sigma}},\tag{7}\] where \[\begin{align} n_{\alpha\mathbf{k}\sigma}^{(m)} & = & \left|a_{\alpha\mathbf{k}\sigma}^{(m)}\right|^2 + \left|a_{\alpha\mathbf{k}\sigma,\bar\alpha\mathbf{k}\sigma}^{(m)}\right|^2 \\ & + & \left|a_{\alpha\mathbf{k}\sigma,\bar{\bar\alpha}\mathbf{k}\sigma}^{(m)}\right|^2 + \left|a_{1\mathbf{k}\sigma,2\mathbf{k}\sigma,3\mathbf{k}\sigma}^{(m)}\right|^2, \nonumber \end{align}\] and \(a^{(m)}\) are components of the \(m\)th eigenvector of the matrix (4 ).

In this research, we focus on the charge-order phenomenon neglecting possible spin-order formation. Hence, the equations are simplified such that \(n_{\alpha\sigma} = n_\alpha/2\) and \(\epsilon_\alpha \equiv \epsilon_{\alpha\sigma}\).

Since the \(n_\alpha\) take role of both input and output quantities of the model, it can be solved self-consistently. In this research, our point of interest is the zero-temperature phase diagram. However, due to the lack of convergence, the calculations are performed for \(T = 10^{-3}\cdot 4.5t\) which stabilizes the algorithm. Still, rather strong mixing is used: when the Fermi level is in a proximity of a singularity of a spectral function, only \(0.005-0.07\) fraction of a new solution is used for the next iteration. A strict criterion of convergence is used: \(10^{-8}\) for \(n_\alpha\) together with \(10^{-8}\cdot 4.5t\) for \(\Omega/L\). The grid of \(\mathbf{k}\)-points contained \(96 \times 96\) points (817 irreducible points; the \(\mathbf{k}\)-dependent quantities exhibit the symmetry of the non-symmetry-broken triangular lattice (wallpaper group p6m), despite the fact that the supercell has a lower symmetry (p3m1)).

For the sake of better comparison with other lattice models, the quantities are expressed in the units of a half-bandwidth of a non-interacting triangular lattice \(D=4.5t\). For the same reason, the intersite interaction \(V\) is expressed in the units of \(D/z = 0.75t\). In the plots, the displaced chemical potential \(\bar\mu = \mu - U/2 - zV\) is used.

We denote the total concentration \(n = \sum_\alpha n_\alpha / 3\). For the analysis of discontinuous phase transitions, the contributions to the grand potential \(\Omega\) (all per lattice site \(L\)) are defined as the thermal averages of the corresponding terms in the Hamiltonian (1 ): kinetic energy (\(t\)-term contribution), on-site interaction energy (\(U\)-term contribution), intersite-interaction energy (\(V\)-term contribution), potential energy (\(V\)- and \(U\)-term contributions together), and chemical energy (\(\mu\)-term contribution). At finite temperatures \(\Omega\) also has a contribution \(-TS\) (\(S\) is an entropy of the system) which is zero in this research.

One of the advantages of MFA is the ability to plot and analyze a band structure \(\omega(\mathbf{k})\). The Hamiltonian (2 ) can be written in the form of \[\hat{H} = \sum_{\mathbf{k}\sigma} \begin{pmatrix} \hat{c}^\dagger_{1\mathbf{k}\sigma} & \hat{c}^\dagger_{2\mathbf{k}\sigma} & \hat{c}^\dagger_{3\mathbf{k}\sigma} \end{pmatrix} \check{\mathbb{H}}_{1\mathbf{k}\sigma} \begin{pmatrix} \hat{c}_{1\mathbf{k}\sigma} \\ \hat{c}_{2\mathbf{k}\sigma} \\ \hat{c}_{3\mathbf{k}\sigma} \end{pmatrix} + C,\] and the band structure is found as eigenvalues of the matrix \(\check{\mathbb{H}}_{1\mathbf{k}\sigma}\) (i.e., the matrix (5 )). Note that \(\mathcal{Z}_{\mathbf{k}\sigma}=\prod_m \left( 1 + e^{-\beta \lambda^{(m)}_{\mathbf{k}\sigma}} \right)\), where \(\lambda^{(m)}_{\mathbf{k}\sigma}\) is the \(m\)-th eigenvalue of \(\check{\mathbb{H}}_{1\mathbf{k}\sigma}\). Such a formulation is in the language of non-interacting quasiparticles, it gives exactly the same results as the reasoning leading to the equations (6 ) and (7 ).

The spectral function is calculated as \[A_\alpha(\omega + i\eta) = -\frac{1}{\pi}\frac{3}{L}\sum_\mathbf{k} \operatorname{Im} G_{\alpha\alpha\mathbf{k}}(\omega + i\eta),\] where \(\eta=0.001D-0.005D\), and \(G_{\alpha\alpha\mathbf{k}}(z)\) are diagonal elements of a \(3 \times 3\) matrix \(\left( z\check{I} - \check{\mathbb{H}}_{1\mathbf{k}\sigma} \right)^{-1}\). When the spectral function is calculated with such a small \(\eta\), the grid of \(396 \times 396\) \(\mathbf{k}\)-points or more is used to make plots of spectral function clear.

The calculations are implemented in the Python language with the use of Matplotlib library for visualization.

3 Non-Interacting Limit↩︎

For the reference point and the following analysis, the band structures and the densities of states (DOSs) of non-interacting triangular and honeycomb lattices are presented in Fig. 1. Both the no-sublattice case (Fig. 1a) and the 3-sublattice case with the bands folded in accordance with \(\sqrt{3}\times\sqrt{3}\) supercell (Fig. 1b) are shown.

The DOS (per spin) of a triangular lattice is [77]: \[\rho_\text{tri}(\omega) = \frac{K(z_1/z_0)}{\pi^2t\sqrt{z_0}},\] where: \(K(m)\) is the complete elliptic integral of the first kind, \[\begin{align} & & K(m) = \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1 - m\sin^2\theta}}, \\ & & z_0 = \begin{cases} 3 + 2\sqrt{3 - \omega/t} - \left(\frac{\omega/t}{2}\right)^2, &\text{for } 2 \leq \omega/t \leq 3 \\ 4\sqrt{3 - \omega/t}, &\text{for } -6 \leq \omega/t \leq 2 \end{cases}, \nonumber \\ & & z_1 = \begin{cases} 4\sqrt{3 - \omega/t}, &\text{for } 2 \leq \omega/t \leq 3 \\ 3 + 2\sqrt{3 - \omega/t} - \left(\frac{\omega/t}{2}\right)^2, &\text{for } -6 \leq \omega/t \leq 2 \end{cases}. \nonumber \end{align}\] It is asymmetric, exists from \(-6t\) to \(3t\), and has a van Hove singularity at \(2t\). The zero-temperature half-filling would correspond to a Fermi level \(\omega_\text{F} \approx 0.8347t\).

The DOS (per spin) of a honeycomb lattice is [78], [79]: \[\rho_\text{hcb}(\omega) = 2 |\omega| \rho_\text{tri}(3 - \omega^2).\] It is symmetric around the Dirac cone at \(\omega=0t\), exists from \(-3t\) to \(3t\), and has van Hove singularities at \(\pm t\) (Fig. 1c).

4 Phase Diagram↩︎

The phase diagram is shown in Fig. 2 while the structures of the found charge orders are schematically presented in Fig. 3. In the model (1 ), the lattice is fully unoccupied (I\(_{000}\)) for \(\bar\mu < -(zV + U/2 + 6t)\) and fully occupied (I\(_{222}\)) for \(\bar\mu > zV + U/2 + 3t\). The corresponding regions on the phase diagrams are indicated by gray color in Fig. 2.

The metallic phase without charge order (CO) is denoted as M\(_{AAA}\) and marked by the red color. Its band structure and spectral function fully repeat the non-interacting ones (Fig. 1b) except the Fermi level is shifted by: \[\omega_\text{shift} = \left(zV+\frac{U}{2}\right)(n - 1) - \bar\mu,\] (cf. Eq. (3 ) when all \(n_{\alpha\sigma}\) are equal). It is possible to have convergence of an algorithm for the non-charge-ordered phase throughout the whole phase diagram. However, the further it goes into the charge-ordered regions of the phase diagram, the easier it can be destabilized toward CO phases by introduction of small deviations in the input \(n_\alpha\), and the larger the difference between its grand potential and the grand potential of the CO phases.

The blue and yellow regions of the phase diagram in Fig. 2 mark the CO phases with the sublattice occupation numbers \(n_1=n_A\), \(n_2=n_3=n_B\) (\(ABB\) phases, Fig. 3c) and \(n_1=n_2=n_A\), \(n_3=n_B\) (\(AAB\) phases, Fig. 3d), respectively, with \(n_A > n_B\). These phases have degeneracy of 3, e.g., the possible sublattice occupation numbers \(n_\alpha\) of the yellow region (\(AAB\) region) are \((n_A, n_A, n_B)\), \((n_A, n_B, n_A)\), and \((n_B, n_A, n_A)\).

The yellow region (subsection 4.1) encompasses two CO metallic phases (M\(_{AAB}^1\) and M\(_{AAB}^2\)) and an insulating one (I\(_{AAB}\)). The total concentration in the latter is always \(4/3\), while in the atomic-limit case it has the sublattice occupation numbers \((2,2,0)\). The M\(_{AAB}^1\) phase is in fact a so-called pinball liquid (PL) phase, i.e. the electrons cannot hop through one of the sublattices which takes role of pins, while the honeycomb lattice formed by the other two sublattices is metallic. The phase loses its PL properties as it approaches a non-ordered phase (the dotted line).

The blue region in Fig. 2 (subsection 4.2) encompasses two CO metallic phases (M\(_{ABB}^1\) and M\(_{ABB}^2\)) and an insulating one (I\(_{ABB}\), \(n=2/3\)) which in the atomic-limit case has the sublattice occupation numbers \((2,0,0)\).

When \(U > 0D\), CO metallic phases with the most broken symmetry occur, i.e., all three occupation numbers in the sublattices are different from each other (M\(_{ABC}\), white region, Fig. 3b, subsection 4.4). The two phases can be identified in this region: an \(ABB\)-like M\(_{ABC}\) phase and an \(AAB\)-like M\(_{ABC}\) phase. The degeneracy of the phases is \(6\).

4.1 AAB Region↩︎

Two continuous phase transitions and the band structures of the \(AAB\) phases (yellow-region phases in Fig. 2) are shown in Fig. 4. For the whole band-structure evolution between \(\bar\mu = 1.50D\) and \(\bar\mu = 2.60D\) for \(U = 2.00D\) and \(zV = 3.50D\) see the file M1_U2.00_V3.50_AAB.mp4 ¹. Since sublattices in this region are occupied by two different occupation numbers only, the useful order parameter is a polarization \[\Delta = \frac{n_A - n_B}{2}.\] The spectral weight \(A_A(\omega) = 2A_1(\omega) = 2A_2(\omega)\) and \(A_B(\omega) = A_3(\omega)\).

When not in the proximity of the M\(_{AAA}\) or I\(_{222}\) phases (see subsection 4.3), the AAB phases can be described as a separate triangular lattice formed by the sublattice with the occupation \(n_B\) (triangular-\(B\) lattice), a honeycomb lattice represented by the sublattices with the occupation \(n_A\) (honeycomb-\(A\) lattice), and a certain amount of hybridization between them (cf. the non-interacting band structures in Fig. 1a and 1c). The hopping between sites of the triangular-\(B\) lattice takes place by means of this hybridization with the honeycomb-\(A\) lattice and hance has a modified amplitude: the opposite sign and a smaller absolute value. As the intersite interaction \(V\) increases, the hybridization decreases. This turns the triangular-like band into a dispersionless atomic level, while the lower bands become indistinguishable from the non-interacting honeycomb band structure (Fig. 1c) with a shifted Fermi level. For a band-structure evolution with the increase of \(V\) see the file M2_U2.00_mu2.00_AAB.mp4 [81]. The above is valid for all \({AAB}\) phases.

Noticeably, a position of a Dirac cone at the K-point of the honeycomb-like bands and a position of an extremum at the K-point of the triangular-like band are always fixed at \[\begin{align} \omega_A = \omega_\text{shift} - \frac{\omega_\Delta}{3}, &&\quad \omega_B = \omega_\text{shift} + \frac{2\omega_\Delta}{3}, \end{align}\] respectively, where the distance between them: \[\omega_\Delta = (zV - U)\Delta.\] The other way to write these expressions is \[\label{bandlocation} \omega_\alpha = \omega_\text{shift} + (zV-U)\frac{n - n_\alpha}{2}\tag{8}\] (in this subsection, \(\omega_A = \omega_1 = \omega_2\) and \(\omega_B = \omega_3\)).

To get more insights on \(AAB\) phases, the lower bands (honeycomb-like bands) are shown in Fig. 5 together with \(\pm 3t\) and \(\pm t\) lines where the non-interacting honeycomb band structure has edges and van Hove singularities. It turns out that the upper edge and the lower singularity are always fixed at \(\omega_A + 3t\) and \(\omega_A - t\), respectively. It justifies the lack of dependency on the interaction \(V\) of a left boundary of the I\(_{AAB}\) phase in Fig. 2. Consequently, the band gap of the I\(_{AAB}\) phase: \[\omega_\text{gap} = \omega_\Delta - 3t = (zV - U)\Delta - 3t\] (asymptotically approaches \(zV - U - 3t\) for a large \(V\)).

On the energy scale the hybridization with the triangular-\(B\) lattice is concentrated around the lower edge and the upper singularity of the honeycomb-like bands. It is clear from the non-zero spectral weight of the triangular-\(B\) lattice. Locations of these points are not fixed to \(-3t\) and \(t\), but approach these values in the large-\(V\) limit where the hybridization disappears (see the file M2_U2.00_mu2.00_AAB.mp4 [81]).

Since there is no spectral weight of the triangular-\(B\) lattice close to the maximum of the lower bands, the M\(_{ABC}^1\) phase is a pinball liquid. Close to the non-charge-ordered phase, the upper and lower bands start to merge and the phase is not a PL anymore. This transition is continuous and very smooth (see for example Fig. 7d and 7e [81]). To give an approximate boundary of the PL phase, the dotted line in Fig. 2 shows the points where the maximum of the lower bands (located at the \(\Gamma\)-point) and the minimum of the upper band (located at the K-point) have the same energy. From the expressions above it is clear that this line exactly corresponds to \(\omega_\Delta = 3t\).

For \(U=0D\) or close to \(0D\), the I\(_{AAB}\)-M\(_{AAB}^2\) transition is sharp and gets discontinuous for \(\bar\mu \gtrsim 3D\) (see for example Fig. S2a [81]). It can be justified by the proximity of the transition to the fully occupied phase (I\(_{222}\)). The total concentration in I\(_{AAB}\) and I\(_{222}\) phases is always \(4/3\) and \(2\), respectively, while the transition between M\(_{AAB}^2\) and I\(_{222}\) phases is continuous (subsection 4.3). As a consequence the sharp change in occupation numbers occurs at the I\(_{AAB}\)-M\(_{AAB}^2\) transition. When moving from the CO metallic phase, this discontinuous transition corresponds to the release of the potential energy at the expense of the chemical energy.

It is worth noting, since the M\(_{AAB}^2\) phase exists for any large values of \(V\) (when \(\bar\mu\) is also large), and in the large-\(V\) limit the hybridization between the honeycomb-\(A\) and the triangular-\(B\) lattices disappears, it eventually becomes an insulator: the Fermi level is located in the atomic-like peak formed by the sublattice with the occupation \(n_B\) (and both \(n\) and \(n_B\) are not integer!), while the hopping with neighboring sites is restricted. This situation is rather unphysical, and both an inclusion of a next-nearest-neighbor hopping, and an inclusion of the Mott physics can resolve it.

4.2 ABB Region↩︎

The band structures of the \(ABB\) phases (blue region in Fig. 2) and phase transitions between them are shown in Fig. 6. For the whole band structure evolution during these transitions see also the file M3_U2.00_V3.50_ABB.mp4 [81]. Here, the spectral weight \(A_A(\omega) = A_1(\omega)\) and \(A_B(\omega)= 2A_2(\omega) = 2A_3(\omega)\).

In the \(ABB\) region, honeycomb-like bands of a lattice formed by the sublattices with the occupation \(n_B\) (honeycomb-\(B\) lattice) are located above a triangular-like band of a lattice formed by the sublattice with the occupation \(n_A\) (triangular-\(A\)). It is clear, that the upper bands look less like the honeycomb bands (Fig. 1c) for \(zV=3.50D\) in comparison to the \(AAB\) lower bands for the same \(V\). It is justified by the fact that the 1st and the 2nd bands in Fig. 1b can be transformed to the honeycomb bands (Fig. 1c) with less effort than the 2nd and the 3rd bands. When \(V\) increases, the bands change accordingly to form non-interacting honeycomb-like and triangular-like (eventually, an atomic-like) bands while the hybridization between them disappears: see the file M4_U2.00_mu-1.80_ABB.mp4 [81].

Positions of the K-point extremum of the triangular-like band and the Dirac cone of the honeycomb-like bands are still fixed to the energies defined by Eq. (8 ) (\(\omega_A=\omega_1\) and \(\omega_B=\omega_2=\omega_3\)) or, in terms of the polarization, \[\begin{align} \omega_A = \omega_\text{shift} - \frac{2\omega_\Delta}{3}, &&\quad \omega_B = \omega_\text{shift} + \frac{\omega_\Delta}{3}. \end{align}\]

Similar to the \(AAB\) phases, the upper edge and the lower singularity of the honeycomb-like bands are also fixed to \(3t\) and \(-t\) in contrast to the lower edge and the upper singularity that approach \(-3t\) and \(t\) only in the large-\(V\) limit when the hybridization disappears. Since, the lower edge is now adjacent to a band gap, the \(ABB\) phases are noticeably different from the \(AAB\) phases. The right boundary of the I\(_{ABB}\) phase in Fig. 2 is dependent on \(V\). The expression of the lower limit only of the I\(_{ABB}\)-phase band gap can be written: \[\omega_\text{gap} \ge (zV - U)\Delta - 3t.\] When \(V\) increases, the \(\omega_\text{gap}\) asymptotically approaches the right-side expression.

Since the hybridization between the honeycomb-\(B\) and triangular-\(A\) lattices is concentrated around the minimum of the honeycomb-like bands, exactly where the Fermi level of the M\(_{ABB}^1\) phase is located, the phase is not a PL. In the limit of large \(V\) the hybridization disappears; hence, in this limit, the M\(_{ABB}^1\) phase would eventually become a PL, but the border of such a transition is vague. Similarly to the M\(_{AAB}^2\), the M\(_{ABB}^2\) phase in large-\(V\) limit is an insulator with non-integer occupations which can be viewed as an artifact of a model or its mean-field treating.

In contrast to the \(AAB\) region, the M\(_{ABB}^2\) phase does not exist for \(U=0D\) or close to \(0D\). The discontinuous transition to the non-charge-ordered phase (subsection 4.3) takes place directly from I\(_{ABB}\) phase.

4.3 Symmetry Breaking from a Non-Charge-Ordered Phase↩︎

Throughout most of the range of \(\mu\) the transitions from the metallic non-charge-ordered phase (M\(_{AAA}\), red region in Fig. 2) to the CO phases (M\(_{ABB}^1\), M\(_{ABB}^2\), I\(_{ABB}\), M\(_{AAB}^1\), M\(_{AAB}^2\), I\(_{AAB}\)) are discontinuous (see for example Fig. 7 [81]). Noticeably, these transitions take place for \(zV \gtrsim U + D\). The discontinuous symmetry breaking corresponds to the release of intersite-interaction energy at the expense of both kinetic and on-site-interaction energy. Additionally, the chemical energy decreases with transitions to the M\(_{ABB}^2\) and I\(_{ABB}\) phases, and increases with transitions to the M\(_{AAB}^2\) and I\(_{ABB}\) phases. The CO phases are metastable inside the M\(_{AAA}\) region for a rather small range of parameters: on the \(zV\) axis the span of this range is up to \(0.1D\) for the most noticeable hysteresis (I\(_{ABB}\)-M\(_{AAA}\) transition for \(U=0D\) and large hole doping).

The analysis of band structures shows that the CO phases get unstable when the Fermi level approaches certain valleys or saddle points. Moreover, there are ranges of \(\mu\) where the symmetry-breaking transitions to the CO metallic phases are continuous (see continuous evolutions in the files M*_SymBreak.mp4 [81]).

4.3.1 M\(_{ABB}^1\) and I\(_{ABB}\) Phases↩︎

The CO M\(_{ABB}^1\) and I\(_{ABB}\) phases become unstable when the Fermi level approaches (from above) the maximum of the lower (triangular-like) band located at the K-point (see their phase diagrams in Fig. 6 for the reference).

For the continuous transition to the M\(_{AAA}\) phase, the maximum of the triangular-like band should merge with the Dirac point of the honeycomb-like bands at the K-point. The Fermi level cannot be above the Dirac point because it destabilizes the M\(_{ABB}^1\) phase towards an \(AAB\) phase. However, the Fermi level stays between the maximum of the triangular-like band and the Dirac point close to the triple phase-transition point of the M\(_{AAA}\), M\(_{ABB}^1\), and M\(_{AAB}^1\) phases. There, the continuous symmetry breaking to the M\(_{ABB}^1\) phase takes place. It is visualized in the file M5_U2.00_mu-0.85_ABB_SymBreak.mp4 [81].

Since the Fermi level appears where all three bands of the non-charge-ordered phase are degenerate at the K-point, the \(\omega_\text{shift} = 0t\) around the continuous transition, i.e. \[\label{triplepoint} \bar\mu = \left(zV + \frac{U}{2}\right) (n - 1),\tag{9}\] where \(n \approx 0.7973\) is the non-interacting triangular-lattice occupation that corresponds to \(\omega_\text{F} = 0t\). Noticeably, around the triple point the \(zV \approx U+D\) which brings us to the expression: \[\frac{\bar\mu}{D} \approx 0.10 - 0.30\frac{zV}{D}\] for the triple point.

4.3.2 M\(_{AAB}^1\) Phase↩︎

Similarly, the M\(_{AAB}^1\) phase (but not I\(_{AAB}\)) becomes unstable towards the M\(_{AAA}\) phase when the Fermi level approaches (from below) the minimum of the upper (triangular-like) band at the K-point (see its phase diagram in Fig. 4 for the reference). The same as for the M\(_{ABB}^1\) phase, the continuous phase transition (merging the extremum of the triangular-like band and the Dirac point at the K-point) happens around the same triple point of the phase diagram where Eq. (9 ) is satisfied with the same \(n\). The continuous transition is visualized in the file M6_U2.00_mu-0.80_AAB_SymBreak.mp4 [81].

4.3.3 M\(_{ABB}^2\) Phase↩︎

For \(U\) close to \(0D\), there is no transition between the M\(_{AAA}\) and M\(_{ABB}^2\) phases. For a larger \(U\), the CO phase becomes unstable towards the M\(_{AAA}\) phase in two ways: 1) when the Fermi level approaches (from above) the singularity (saddle point) of the lower (triangular-like) band at the M-point; 2) when the Fermi level approaches (from below) the minimum of the upper (honeycomb-like) bands which in this range of parameters is at the M-point.

The continuous transition between the M\(_{AAA}\) and M\(_{ABB}^2\) phase happens in the small range of \(\mu\) when the Fermi level stays in between of the mentioned saddle point and minimum at the M-point. It is visualized in the file M7_U2.00_mu-1.80_ABB_SymBreak.mp4 [81]. Since the Fermi level appears where the two bands of the M\(_{AAA}\) phase are degenerate at the M-point, the \(\omega_\text{shift} = t\) around the continuous transition, i.e. \[\bar\mu + t = \left(zV + \frac{U}{2}\right) (n - 1),\] where \(n \approx 0.6078\) is the non-interacting triangular-lattice occupation when \(\omega_\text{F} = -t\).

4.3.4 M\(_{AAB}^2\) and I\(_{AAB}\) Phases↩︎

The M\(_{AAB}^2\) and I\(_{AAB}\) phases become unstable towards the M\(_{AAA}\) phase when the Fermi level approaches (from above) the maximum of the lower (honeycomb-like) bands at the \(\Gamma\)-point.

In contrast to the M\(_{ABB}^2\), the M\(_{AAB}^2\) phase cannot continuously transform to the M\(_{AAA}\) phase but has an infinite range of \(\mu\) (particularly, \(\bar\mu \gtrsim 3(U/2+D)\)) where it continuously transforms to the fully occupied phase (I\(_{222}\)). The continuous transition is visualized in the file M8_U0.00_mu3.10_AAB_SymBreak.mp4 [81].

At this continuous transition the \(\omega_\text{shift} = -3t\) with \(n=2\) which corresponds to the top of the non-interacting triangular-lattice bands: \[\bar\mu - 3t = zV + \frac{U}{2}.\]

4.4 AAB-ABB Transition and ABC Region↩︎

The transition between the M\(_{ABB}^1\) and M\(_{AAB}^1\) phases (between the yellow and blue regions) is discontinuous (see for example Fig. 8a [81]). When moving from the \(ABB\) phase, it happens with the release of chemical energy at the expense of intersite-interaction energy.

For \(U>0D\) an intermediate region of M\(_{ABC}\) phases (\(n_1=n_A\), \(n_2=n_B\), \(n_3=n_C\), \(n_A>n_B>n_C\), white region) opens between the \(ABB\) and \(AAB\) regions (see for example Fig. S3d [81]). For the most range of parameters the transition between the \(ABC\) region and the M\(_{ABB}^1\) phase on one side and the M\(_{AAB}^1\) phase on the other side is continuous. Note that the less-symmetry-broken phases (\(ABB\) and \(AAB\)) can have the self-consistent solution of the Hamiltonian (2 ) inside the \(ABC\) region despite other indicators that the phase transition is continuous.

The \(ABC\) region can be characterized by identifying two phases: the \(ABB\)-like M\(_{ABC}\) phase and the \(AAB\)-like M\(_{ABC}\) phase with the discontinuous transition between them. Moreover, there is a range of \(V\) where the \(AAB\)-like M\(_{ABC}\) phase discontinuously transforms directly into the M\(_{ABB}^1\) phase (see Fig. 2). The same as for the M\(_{ABB}^1\)-M\(_{AAB}^1\) transition, these discontinuous transitions take place with the release of chemical energy at the expense of intersite-interaction energy when moving from the \(ABB\)-like side.

There is a range of parameters where the \(ABC\) phase can be viewed as both \(ABB\)-like and \(AAB\)-like M\(_{ABC}\) phase (in Fig. 2 the region around \(zV=4.4D\) for \(U=2D\)). The visualization of two continuous transitions between the \(ABB\), \(ABC\), and \(AAB\) regions for \(U=2D\) and \(V=4.4D\) is presented in the file M9_U2.00_V4.40 _ABB-ABC-AAB.mp4 [81].

The band structure of the M\(_{ABC}\) phase consists of \(3\) bands: a lower triangular-like which spectral weight primarily comes from the \(A\)-sublattice, an upper triangular-like which spectral weight primarily comes from the \(C\)-sublattice, and an intermediate band which spectral weight primarily comes from the \(B\)-sublattice. The Fermi level is located inside the intermediate band. All three bands have an extremum at the K-point that is fixed to the energy defined from Eq. (8 ) (\(\omega_A=\omega_1\), \(\omega_B=\omega_2\), and \(\omega_C=\omega_3\)). In the large-\(V\) limit the three bands are rather three atomic levels: a fully occupied level, a fully unoccupied level, and a level where the Fermi level is located.

4.5 Correspondence to the Atomic-Limit Phase Diagram↩︎

For a large \(V\) the atomic-limit (\(D=0\)) solution of the triangular-lattice EHM (exact for the ground state) predicts \(7\) phases that are denoted according to the three occupation numbers: 000, 100, 200, 210, 220, 221, and 222 [75]. The phase transitions take place at the following lines: \[\begin{align} \label{atomic-limit} \bar\mu &= &\mp (zV + U/2),\;\text{for 000-100 and 222-221}, \nonumber\\ \bar\mu &= &\mp (zV - U/2),\;\text{for 100-200 and 221-220}, \quad \\ \bar\mu & = &\mp U/2, \quad \;\;\qquad \text{for 200-210 and 220-210}, \nonumber \end{align}\tag{10}\] where the ‘\(\mp\)’ corresponds to the first or second boundary mentioned above, respectively. It is clear from the formulas that the phases 100, 210 and 221 do not exist for \(U=0D\).

There is a correspondence between the atomic-limit and mean-field solutions at a large \(V\). Besides the obvious correspondence for the fully-unoccupied and fully-occupied phases, the atomic-limit 100 phase can be associated with the mean-field M\(_{ABB}^2\) phase, the 200 phase with the I\(_{ABB}\) and M\(_{ABB}^1\) phases together, the 210 phase with M\(_{ABC}\) phase, and correspondingly for the other side of the phase diagram.

Thus, the disappearance of the M\(_{ABC}\) and M\(_{ABB}^2\) phases for \(U=0D\) could be predicted from the atomic-limit solution. Despite this correspondence, the mean-field phase M\(_{AAB}^2\) exists even for \(U=0D\) and large \(V\), which makes it an unexpected effect of non-zero itineracy.

To compare the atomic-limit and the mean-field solutions, one can put \(U=2D\) in Eqs. (10 ) (despite the fact that \(D=0\) in the atomic limit): \[\begin{align} \bar\mu & = &\mp (zV + D), \;\text{for 000-100 (`--') and 222-221 (`+')}, \nonumber\\ \bar\mu & = &\mp (zV - D), \;\text{for 100-200 (`--') and 221-220 (`+')}, \nonumber\\ \bar\mu & = &\mp D, \qquad \;\;\quad \text{for 200-210 (`--') and 220-210 (`+')}. \nonumber \end{align}\] The mean-field boundaries of M\(_{ABC}\) phase at large \(V\) and \(U=2D\) indeed slowly approach \(\bar\mu = \mp D\) (Fig. 2). The same is valid for the boundaries of the I\(_{000}\) and I\(_{222}\) phases since the terms \(-6t\) and \(3t\) become negligible in comparison to the large \(zV\). The left boundary of the I\(_{ABB}\) phase coincide well with the left boundary of the atomic-limit 200 phase even at the small values of \(V\), and the same is valid for the 220 and I\(_{AAB}\) phases.

It is reasonable to expect that primarily those phases that are associated with non-zero \(U\) (M\(_{ABC}\), M\(_{ABB}^2\), M\(_{AAB}^2\)) will be modified when taking into account the Mott physics (local correlations). Note also, that in the atomic limit there are also 111 (Mott insulator), 110, and 211 phases for small \(V\) (particularly, for \(zV \leq 2U\)).

5 Summary↩︎

The full zero-temperature phase diagram of the extended Hubbard model on the triangular lattice in the all-encompassing range of chemical potential values and repulsive onsite and nearest-neighbor electron interactions is presented and analyzed. Despite the limitations of the mean-field approximation, restriction to the \(\sqrt{3}\times\sqrt{3}\) supercells and the absence of magnetic order, the large variety of features are found. Those are the diverse and numerous phase transitions, including the continuous and discontinuous symmetry breakings; the pinball liquid phase; the strong particle-hole asymmetry manifested in every phase transition with a charge-ordered phase. Both the found phases and the found phase transitions are accompanied by the extensive band-structure analysis that provided clarity and the understanding of the phase diagram.

It is worth noting, that the quarter-filling for electrons (\(n=0.5\)) and holes (\(n=1.5\)) that is common for some organic conductors occurs in the M\(_{AAA}\), M\(_{ABB}^2\), and M\(_{AAB}^2\) phases, which are not pinball-liquid phases within the MFA. Nevertheless, for a large enough \(U\), we expect that the M\(_{ABB}^2\), M\(_{AAB}^2\), and M\(_{ABC}\) phases should experience prominent Mott localization at one of the sublattices which is consistent with the correspondence to the exact atomic-limit results. This localization may contribute to formation of a pinball liquid [49]. Moreover, for a region \(zV \lesssim 2U\) we except the appearance of a non-ordered Mott insulator and the charge-ordered phases that correspond to the 110 and 211 phases of an atomic-limit phase diagram [75]. Additionally, when considered, magnetic orders might be found particularly for \(zV < U/2\) [54], [67], but detailed analysis of the magnetism in the system is out of the scope of this work.

It is evident that the future investigation beyond the mean-field approximation can benefit from our analysis. Additionally, our research can easily be expanded to include next-nearest-neighbor interactions and finite temperatures. The more work is also required to include magnetic order and to consider the orderings that require larger supercells.

In this Supplemental Material we present additional results which clarify the findings included in the main text and extend the discussion, in particular, concerning:

• the dependencies of various observables near phase transitions to supplement the discussion of the phase diagram (see Sec. 6 including the figures and multimedia description)

• the full list of all transitions found in the system in the used approach (Sec. 7).

6 The Dependencies of Order parameters and Grand-potential Contributions↩︎

Figures 7, 8, and 9 show the order parameters and the contributions in the grand potential as a function of \(V\), \(\bar{\mu}\), and \(U\), respectively. The metastability of the phases (hysteresis) is also shown with dotted lines. In some cases the spectral weight at the Fermi level (\(A(i\eta)\)) is shown on two different scales due to its proximity to a spectral-weight singularity in the \(M_{ABB}^1\) and \(M_{ABC}\) phases. In particular, the figures contain the results for the following model parameters:

• Fig. 7 shows \(V\)-dependencies for \(U=2.00 D\) and various values of \(\bar{\mu}\): (a) \(\bar{\mu} =- 3.00 D\), (b) \(\bar{\mu} = 3.00 D\), (c) \(\bar{\mu} =- 1.25 D\), (d) \(\bar{\mu} =- 0.30 D\), and (e) \(\bar{\mu} = 1.00 D\);

• Fig. 8 shows \(\bar\mu\)-dependencies for (a) \(U=0.00 D\) and \(zV = 3.50 D\), (b) \(U=0.00 D\) and \(zV = 4.00 D\), and (c) \(U=2.00 D\) and \(zV = 7.00 D\);

• Fig. 9 shows \(U\)-dependencies for (a) \(\bar{\mu} =- 3.00 D\) and \(zV = 4.00 D\), (b) \(\bar{\mu} = 3.00 D\) and \(zV = 4.00 D\), (c) \(\bar{\mu} =- 0.30 D\) and \(zV = 5.50 D\), and (d) \(\bar{\mu} = 0.00 D\) and \(zV = 9.00 D\).

In addition, we present nine multimedia movies that show the evolution of the band structure and density of states:

1. M1_U2.00_V3.50_AAB.mp4
   — evolution between \(\bar\mu = 1.50D\) and \(\bar\mu = 2.60D\) for \(U = 2.00D\) and \(zV = 3.50D\) (transitions M\(_{AAB}^1\)-I\(_{AAB}\) and I\(_{AAB}\)-M\(_{AAB}^2\))

2. M2_U2.00_mu2.00_AAB.mp4
   — evolution between \(zV = 3.50D\) and \(zV = 9.00D\) for \(U = 2.00D\) and \(\bar{\mu} = 2.00D\)

3. M3_U2.00_V3.50_ABB.mp4
   — evolution between \(\bar\mu = -2.38D\) and \(\bar\mu = -1.20D\) for \(U = 2.00D\) and \(zV = 3.50D\) (transitions M\(_{ABB}^2\)-I\(_{ABB}\) and I\(_{ABB}\)-M\(_{ABB}^1\))

4. M4_U2.00_mu-1.80_ABB.mp4
   — evolution between \(zV = 3.50D\) and \(zV = 9.00D\) for \(U = 2.00D\) and \(\bar{\mu} = -1.80D\)

5. M5_U2.00_mu-0.85_ABB_SymBreak.mp4
   — evolution between \(zV = 2.95D\) and \(zV = 3.45D\) for \(U = 2.00D\) and \(\bar{\mu} = -0.85D\) (transition M\(_{AAA}\)-M\(_{ABB}^1\))

6. M6_U2.00_mu-0.80_AAB_SymBreak.mp4
   — evolution between \(zV = 2.95D\) and \(zV = 3.45D\) for \(U = 2.00D\) and \(\bar{\mu} = -0.80D\) (transition M\(_{AAA}\)-M\(_{AAB}^1\))

7. M7_U2.00_mu-1.80_ABB_SymBreak.mp4
   — evolution between \(zV = 2.95D\) and \(zV = 3.45D\) for \(U = 2.00D\) and \(\bar{\mu} = -1.80D\) (transitions M\(_{AAA}\)-M\(_{ABB}^2\) and M\(_{ABB}^2\)-I\(_{ABB}\))

8. M8_U0.00_mu3.10_AAB_SymBreak.mp4
   — evolution between \(zV = 2.35D\) and \(zV = 2.85D\) for \(U = 0.00D\) and \(\bar{\mu} = 3.10D\) (transition I\(_{222}\)-M\(_{AAB}^2\))

9. M9_U2.00_V4.40_ABB-ABC-AAB.mp4
   — evolution between \(\bar{\mu} = -0.70D\) and \(\bar{\mu} = -0.10D\) for \(U = 2.00D\) and \(zV = 4.40D\) (transitions M\(_{ABB}^1\)-M\(_{ABC}\) and M\(_{ABC}\)-M\(_{AAB}^1\))

7 Transitions Found in the System↩︎

The list of 24 phase transitions:

1. \(ABB\)-metal-insulator transitions (blue region):

   1. M\(_{ABB}^1\) to I\(_{ABB}\) (continuous);

   2. M\(_{ABB}^2\) to I\(_{ABB}\) (continuous);

2. \(AAB\)-metal-insulator transitions (yellow region):

   1. M\(_{AAB}^1\) to I\(_{AAB}\) (continuous);

   2. M\(_{AAB}^2\) to I\(_{AAB}\) (continuous);

   3. M\(_{AAB}^2\) to I\(_{AAB}\) (discontinuous), for large \(\bar\mu\) but small \(U\);

3. M\(_{AAA}\) to \(ABB\) transitions (red to blue):

   1. M\(_{AAA}\) to M\(_{ABB}^1\) (discontinuous);

   2. M\(_{AAA}\) to M\(_{ABB}^1\) (continuous), a small region in the very proximity of M\(_{ABB}^1\)-M\(_{AAB}^1\) transition;

   3. M\(_{AAA}\) to M\(_{ABB}^2\) (discontinuous), not for small \(U\);

   4. M\(_{AAA}\) to M\(_{ABB}^2\) (continuous), a small region, not for small \(U\);

   5. M\(_{AAA}\) to I\(_{ABB}\) (discontinuous);

4. M\(_{AAA}\) to \(AAB\) transitions (red to yellow):

   1. M\(_{AAA}\) to M\(_{AAB}^1\) (discontinuous);

   2. M\(_{AAA}\) to M\(_{AAB}^1\) (continuous), a small region in the very proximity of M\(_{ABB}^1\)-M\(_{AAB}^1\) transition;

   3. M\(_{AAA}\) to M\(_{AAB}^2\) (discontinuous);

   4. M\(_{AAA}\) to I\(_{AAB}\) (discontinuous);

5. M\(_{ABB}^1\)-M\(_{AAB}^1\) transition (blue to yellow) (discontinuous);

6. \(ABC\) region (blue-white-yellow):

   1. M\(_{ABC}\) to M\(_{ABB}^1\) (continuous);

   2. M\(_{ABC}\) to M\(_{AAB}^1\) (continuous);

   3. inside the M\(_{ABC}\) the discontinuous transition from \({ABB}\)-like M\(_{ABC}\) phase to \({AAB}\)-like M\(_{ABC}\) phase. As shown on phase diagram, there is a region where this transition disappears (the line is terminated)

   4. directly from \({AAB}\)-like M\(_{ABC}\) phase to M\(_{ABB}^1\) phase (discontinuous);

7. The M\(_{AAB}^1\) phase is not a PL as it gets closer to M\(_{ABB}^1\)-M\(_{AAA}\)-M\(_{AAB}^1\) triple point. The approximate border between M\(_{AAB}^1\)-PL and M\(_{AAB}^1\)-non-PL is the dotted line;

8. In the limit of large \(V\), the M\(_{ABB}^1\) phase becomes a PL because the hybridization between triangular and honeycomb lattices disappears, but the border of this phase is hard to define;

9. Additionally, there is one transition from M\(_{AAA}\) to the fully unoccupied phase, and two transitions to the fully occupied phase: from M\(_{AAA}\) and directly from M\(_{AAB}^2\) phase.

References↩︎