Extended Falicov-Kimball model: Exact solution for finite temperatures

TL;DR

This paper provides an exact finite-temperature solution for the extended Falicov-Kimball model in large dimensions, revealing phase transitions, phase diagrams, and the effects of onsite and intersite interactions on electronic states.

Contribution

It offers the first exact finite-temperature analysis of the extended Falicov-Kimball model including $U$ and $V$, detailing phase diagrams and insulator-metal transitions.

Findings

Identification of temperature-dependent density of states (DOS)

Discovery of insulator-metal transition temperature $T_{MI}$

Detection of eight different ordered phases and their transitions

Abstract

The extended Falicov-Kimball model is analyzed exactly for finite temperatures ($T\geq0$) in the limit of large dimensions. Onsite and intersite density-density interactions $U$ and $V$ are included in the model. Using the dynamical mean field theory formalism on the Bethe lattice we find rigorously the temperature dependent density of states (DOS) at half-filling. At $T=0$ the system is ordered to form the checkerboard pattern and the DOS has the gap $\Delta(\varepsilon_F) > 0$ at the Fermi level, if only $U\neq 0$ or $V\neq 0$. If $U <0$ or $U > 2V$, two additional subbands develop inside the principal energy gap. They become wider with increasing $T$ and at a certain $U$- and $V$-dependent temperature $T_{MI}$ they join with each other at $\varepsilon_F$. Since above $T_{MI}$ the DOS is positive at $\varepsilon_F$, we interpret $T_{MI}$ as the transformation temperature from insulator to metal. Moreover, we show that if $V\lesssim 0.54$ then $T_{MI}=0$ at two quasi-quantum critical points $U_{cr}^{\pm}$ (one positive and the other negative), whereas for $V\gtrsim 0.54$ there is only one negative $U^-_{cr}$. Having calculated the temperature dependent DOS we study thermodynamic properties of the system starting from its free energy and then we construct the phase diagrams in the variables $T$ and $U$ for a few values of $V$. Our calculations give that inclusion of the intersite coupling $V$ causes the finite temperature phase diagrams to become asymmetric with respect to a change of sign of $U$. On these phase diagrams we detected stability regions of eight different kinds of ordered phases, where both charge-order and antiferromagnetism coexists (five of them are insulating and three are conducting) and three different nonordered phases (two of them are insulating and one is conducting). Moreover, both continuous and discontinuous transitions between various phases were found.

Equations

$$\hat{H}={\hat{H}}_t+{\hat{H}}_U+{\hat{H}}_V+{\hat{H}}_{\mu },$$

$${\hat{H}}_t$$

$${\hat{H}}_V$$

$$G^{+}(z)=\frac{z+v+\frac{1}{2}Ud-G^{-}(z)}{[z+v+\frac{1}{2}U-G^{-}(z)][z+v-\frac{1}{2}U-G^{-}(z)]}$$

$$G^{-}(z)=\frac{z-v-\frac{1}{2}Ud-G^{+}(z)}{[z-v+\frac{1}{2}U-G^{+}(z)][z-v-\frac{1}{2}U-G^{+}(z)]},$$

$$d=n_{\uparrow }^{+}-n_{\uparrow }^{-},$$

$$d_1=n_{\downarrow }^{+}-n_{\downarrow }^{-}.$$

$${\Delta }_Q=\frac{1}{2}\left(d+d_1\right)\text{and}{m}_Q=\frac{1}{2}\left(d-d_1\right),$$

$$G^{\pm }(z)=\frac{4z^2-A^2-\sqrt{(4z^2-A^2)(4z^2-A^2-16)}}{4(2z\pm A)},$$

$$W(G^{+})$$

$${\rho }^{\pm }(U,V,d,d_1;\epsilon )=-\frac{1}{\pi }\text{Im}\left[G^{\pm }(U,V,d,d_1;\epsilon +i0)\right].$$

$$n_{\downarrow }^{\pm }={\int }_{-\infty }^{{\epsilon }_F}\frac{{\rho }^{\pm }(U,V,d,d_1;\epsilon )}{1+\exp \left[\epsilon /\left(k_{B}T\right)\right]}d\epsilon$$

$$\frac{dW(G^{+})}{d{G}^{+}}$$

$$\rho(\varepsilon_{F})=\left\{\begin{array}[]{r}\frac{1}{\pi}\quad\text{if }U=0\text{ and }(d+d_{1})V=0\\ \\ \quad\frac{1}{\pi}\left|\frac{\textrm{Im}\sqrt{w(U,V,d,d_{1})}}{8[(d+d_{1})V]^{2}-2U^{2}}\right|\quad\text{otherwise},\end{array}\right.$$

$$w(U,V,d,d_1)=U^2(U^4-4U^2+4d^2)$$

$$+8(U^2-2)dU(d+d_1)V-8(U^4-2U^2-2)[(d+d_1)V{]}^2$$

$$-32dU[(d+d_1)V{]}^3+16U^2[(d+d_1)V{]}^4.$$

$$\Delta(\varepsilon_{F})_{T\rightarrow 0^{+}}=\left\{\begin{array}[]{ll}\left|2V(1+d_{1})-U\right|&\quad\text{if }0\leq U<2V\\ &\\ \left|\left(\sqrt{1+4U^{2}-4(1+d_{1})UV+4(1+d_{1})^{2}U^{2}V^{2}}-U^{2}-1\right)/U\right|&\quad\text{if }U<0\quad\text{or}\quad 2V<U.\end{array}\right.$$

$$\Delta(\varepsilon_{F};U=2V)_{T\rightarrow 0^{+}}=\left\{\begin{array}[]{ll}Ud_{1}&\quad\text{if }d_{1}>0\\ &\\ \left|\left(\sqrt{1+4U^{2}-2(1+d_{1})U^{2}+(1+d_{1})^{2}U^{4}}-U^{2}-1\right)/U\right|&\quad\text{if }d_{1}<0.\end{array}\right.$$

$$F=F_{el}+F_{ions}+U/4+V,$$

$$F_{el}$$

$$F_{ions}$$

$$a_0$$

$$a_1$$

$$a_2$$

$$a_3$$

$$a_4$$

$$a_5$$

Full text

††preprint: Submitted to: PHYSICAL REVIEW B

I Introduction

Correlated electron systems exhibit many diverse and interesting properties, as e.g. charge and magnetic ordering, superconductivity, mixed valence, metal-insulator phase transition, etc. Micnas et al. (1990); Imada et al. (1998); Yoshimi et al. (2012); Frandsen et al. (2014); Comin et al. (2015); da Silva Neto et al. (2015); Cai et al. (2016); Hsu et al. (2016); Pelc et al. (2016); Park et al. (2017); Novello et al. (2017). Unfortunately, it is very difficult to find complete solutions even for simple models that describe these systems, therefore many issues are not yet explained satisfactorily Nagaoka (1966); Lieb and Wu (1968); Lieb (1986); Georges et al. (1996); Freericks and Zlatić (2003).

One of a few methods of reliable studying of strongly correlated fermion systems is the dynamic mean field theory (DMFT) van Dongen and Vollhardt (1990); Georges et al. (1996); Freericks and Zlatić (2003); Freericks (2006); Kotliar et al. (2006), which is the exact theory in the limit of high dimensions ($D\rightarrow+\infty$) or, equivalently, of large coordination number. And one of a few models for which this method can be used to achieve accurate results in the thermodynamic limit is the Falicov-Kimball model (FKM) Falicov and Kimball (1969); Brandt and Mielsch (1989, 1990, 1991); Brandt and Urbanek (1992); Freericks and Lemański (2000); Zlatić et al. (2001); Freericks and Zlatić (2003); Chen et al. (2003); Hassan and Krishnamurthy (2007); Lemański and Ziegler (2014); Lemański (2016), sometimes referred to as the simplified Hubbard model. Recently, the FKM has been also investigated by a cluster extension of the DMFT Haldar et al. (2016, 2017a, 2017b). The simplest version of the FKM describes spinless electrons interacting with localized ions via only the local (on-site) Coulomb coupling $U$.

So far, the FKM has been used to describe various effects, such as crystal formation, mixed valence, metal-insulator phase transition, and so on, e.g., Refs. Falicov and Kimball (1969); Georges et al. (1996); Imada et al. (1998); Freericks and Zlatić (2003). In particular, in Refs. Lemański and Ziegler (2014); Lemański (2016) there were analyzed exactly properties of this model (in $D\rightarrow\infty$) related to the order-disorder phase transition caused by a rise of temperature $T$ and the associated to it insulator-conductor transition. The results reported in Lemański and Ziegler (2014) show that the rising $T$ causes an evolution of the density of states (DOS) consisting, e.g., in formation of additional bands within the main energy gap. It turns out that above a certain $U$-dependent temperature $T_{MI}(U)$, $0\leq T_{MI}(U)\leq T_{OD}(U)$, where $T_{OD}(U)$ is the order-disorder transition temperature, the DOS at the Fermi level $\rho(\epsilon_{F})$ becomes positive (in the ordered phase). Moreover, it was fix there the value $U_{cr}$, which points out the quasi-quantum critical point, for which $T_{MI}(U_{cr})=0$, what means that if $U=U_{cr}$, then $\rho(\epsilon_{F})>0$ for any $T>0$ (the details are also discussed in the subsequent sections of this work).

The results obtained for the simplest FKM have prompted us to investigate the extended FKM (EFKM), which, in addition to the local interactions, also includes nonlocal couplings represented by the Coulomb’s repulsion force $V$ between electrons located on neighboring sites of the crystal lattice. As it is well known, in some systems this effect can be quite significant and sometimes it can lead even to a change in nature of the metal-insulator phase transition Lemański et al. (2017); Amaricci et al. (2010); Kapcia et al. (2017a).

The effects of Coulomb’s interactions between electrons located on neighboring lattice sites have already been intensively investigated for the extended Hubbard model, e.g., Refs. Amaricci et al. (2010); Ayral et al. (2012, 2013); Huang et al. (2014); Giovannetti et al. (2015); Kapcia et al. (2017a); Terletska et al. (2017); Ayral et al. (2017); Kapcia et al. (2017b); Terletska et al. (2018) and references therein, while for the EFKM only very few results have been reported so far van Dongen (1992); van Dongen and Leinung (1997); Gajek and Lemański (2004); Čenčariková et al. (2008); Hamada et al. (2017); Lemański et al. (2017). As far as we know, exact results for the EFKM were only reported in Refs. van Dongen (1992); van Dongen and Leinung (1997); Lemański et al. (2017) for $D\rightarrow\infty$, but Refs. van Dongen (1992); van Dongen and Leinung (1997) only refer to the limiting cases of $U\rightarrow 0$ and $U\rightarrow\infty$ and they did not attain the lowest temperatures.

In our previous work Lemański et al. (2017) we obtained the exact solution of the EFKM, but only in the ground state and discussed a few properties at infinitesimally small temperatures. In this work we report our exact results also for the EFKM, but those obtained at finite temperatures. This paper can be also viewed as an extension of the work Lemański and Ziegler (2014), where the simplest version of the FKM was investigated at finite temperatures. Our findings end up with the resulting phase diagram of the EFKM at finite temperatures for a wide range of interaction parameters $U$ and $V$.

Here, we need to emphasize that the inclusion of repulsion between electrons located on neighboring sites significantly increases the level of difficulty in studying the system. This is because in the ordered phase all calculated physical quantities are then expressed not only explicitly by order parameter $d$, as it is in the case of the simplest version of the FKM, but also by additional $d$-dependent parameter $d_{1}(d)$, which for given $d$ has to be determined from the self-consistent equation (all details and the precise definitions of parameters $d$ and $d_{1}$ are given in the next paragraphs). In the ground state the task is relatively simple because then $d=1$, and for given values of $U$ and $V$ parameter $d_{1}$ has a certain fixed value. But at finite temperature it becomes challenging, as then both $d$ and $d_{1}$ depend on $T$.

The rest of the paper is organized as follows. In Sec. II the model considered is introduced (Sec. II.1), and the equations for Green’s functions are determined with dynamical mean field theory (Sec. II.2). Section III is devoted to a discussion of analytical (Sec. III.1) and numerical solutions at $T>0$, including phase diagrams of the model (Sec. III.2). Finally, in Sec. IV the results of this work are summarized and the conclusions are provided.

II Model and method

II.1 Extended Falicov-Kimball model

Here, we study the same Hamiltonian $\hat{H}$ as was used in Ref. van Dongen (1992) and then in Ref. Lemański and Ziegler (2014) on the Bethe lattice. It includes electrons’ kinetic energy term $\hat{H}_{t}$, Coulomb interaction terms (onsite $\hat{H}_{U}$ and intersite $\hat{H}_{V}$) and also $\hat{H}_{\mu}$ term representing an influence of the chemical potential. Thus, the considered Hamiltonian has the following form:

[TABLE]

where

[TABLE]

with $Z$ being the coordination number. $\hat{n}_{i\sigma}=\hat{c}^{+}_{i\sigma}\hat{c}^{\ }_{i\sigma}$ is the occupation number and $\hat{c}^{+}_{i\sigma}$ ($\hat{c}^{\ }_{i\sigma}$) denotes the creation (annihilation) operator of an electron with spin $\sigma=\uparrow,\downarrow$. Electrons with spin $\sigma=\uparrow$ ($\downarrow$) are localized (itinerant), that is why here we call them ions (electrons), respectively. The prefactors in $\hat{H}_{t}$ and $\hat{H}_{V}$ have been chosen such that they yield a finite and non-vanishing contribution to the free energy per site in the limit $Z\rightarrow\infty$. $\left\langle i,j\right\rangle$ denotes the sum over nearest-neighbor pairs. At half-filling, i.e., for $n=1$ ($n=\tfrac{1}{L}\sum_{i,\sigma}\left\langle\hat{n}_{i\sigma}\right\rangle$, $L$ is the number of lattice sites) the chemical potential $\mu$ for the both types of electron is given by $\mu\equiv\mu_{\sigma}=\tfrac{1}{2}U+2V$ (and $n_{\sigma}=1/2$ for both $\sigma=\uparrow,\downarrow$, $n_{\sigma}=\tfrac{1}{L}\sum_{i}\left\langle\hat{n}_{i\sigma}\right\rangle$) van Dongen (1992).

Note that in this work the model is analyzed on the Bethe lattice, which is alternate one, i.e., it can be divided into two equivalent sublattices. We take $t$ as an energy unit, i.e., $t=1$, and basically interaction couplings $U$ and $V$, temperature $k_{B}T$ ($k_{B}$ denotes the Boltzmann constant), gap at the Fermi energy $\Delta(\varepsilon_{F})$, and energies $\varepsilon$ are given in units of $t$. Nevertheless, for a clarity, $t$ will be explicitly given in some expressions (similarly like in Ref. Lemański et al. (2017)). We also assume that $V$ interaction is repulsive, i.e., $V>0$.

II.2 Dynamical mean field theory

The dynamical mean field theory (DMFT) enables exact studies of the correlated electron systems, including EFKM, in the high-dimension limit van Dongen and Vollhardt (1990); Georges et al. (1996); Freericks and Zlatić (2003); Freericks (2006); Kotliar et al. (2006); Lemański and Ziegler (2014); Lemański et al. (2017). Moreover, it was proven that the nonlocal interaction term $V$ can be treated at the Hartree level because the exchange (Fock) and the correlation energy due to the intersite term are negligible in that limit Müller-Hartmann (1989a, b); Metzner and Vollhardt (1989); Metzner (1989).

The basic quantity calculated within the DMFT is the retarded Green’s function $G(z)$, which is defined for the complex $z$ with $\textrm{Im}(z)>0$. Due to the fact that we are dealing with the system composed of two sublattices, we need to determine two Green’s functions $G^{+}$ and $G^{-}$ separately for “$+$” and “$-$” sublattice. Here, we use the Green’s functions derived by van Dongen for the EFKM on the Bethe lattice in the limit of large dimension van Dongen (1992). The formulas [for $t=1$ and at half-filling (i.e., for $\mu=U/2+2V$)] have the following forms:

[TABLE]

where $v=V(d+d_{1})$. In such an approach we have two parameters $d$ and $d_{1}$, which need to be determined self-consistently. $d$ stands for the order parameter (because it is obtained by minimalization of the free energy), which is equal to the difference in mean occupancies of the localized electrons on sublattices ’$+$’ and ’$-$’:

[TABLE]

whereas $d_{1}$ is the difference of the mean occupancies of the itinerant electrons on the both sublattices:

[TABLE]

In the above definitions $n_{\sigma}^{\alpha}:=\langle\hat{n}_{i\sigma}\rangle$ for any $i\in\alpha$, where $\alpha=+,-$ denotes the sublattice index. In fact, $d$ and $d_{1}$ are not independent quantities, as for a given temperature $T$ and for a given parameter $d$ the value of $d_{1}$ can be determined unambiguously (excluding the case of coexistence of two phases at the first-order transition points, as it is discussed further). However, $d$ needs to be found from the condition for a minimum of the free energy.

Notice that due to the equivalence of two sublattices in the Bethe lattice the solution with parameters $d$ and $d_{1}$ is equivalent to the solution in which the parameters have the opposite signs (i.e., in which they are equal to $-d$ and $-d_{1}$, respectively). As a consequence, in the rest of the paper, we consider the solutions with $d\geq 0$ only. With such a choice, as it will be shown further, solutions with both signs of $d_{1}$ can be found depending on values of the interactions and the temperature.

The mentioned parameters can be associated to charge polarization $\Delta_{Q}$ and staggered magnetization $m_{Q}$ by the following relations:

[TABLE]

which create a different, but equivalent, set of parameters. For $d>0$ one gets that $\Delta_{Q}>0$ and $m_{Q}>0$ because $d>|d_{1}|$ for any finite value of $U$ or $V$ (for $d=0$ one always gets $d_{1}=0$, hence also $\Delta_{Q}=0$ and $m_{Q}=0$). For calculations presented in this work it is not relevant, which combination [$(d,d_{1})$ or $(\Delta_{Q},m_{Q})$] is used. Note also that any of $d$, $d_{1}$, $\Delta_{Q}$, and $m_{Q}$ reflects breaking of the system in the ordered phases mentioned below. Additionally, a use of $d$ and $d_{1}$ ensures a simple correspondence to previous works on the FKM van Dongen and Vollhardt (1990); Freericks et al. (1999); Chen et al. (2003); Maśka and Czajka (2006); Hassan and Krishnamurthy (2007); Lemański and Ziegler (2014) and the EFKM van Dongen (1992); Lemański et al. (2017).

At $T=0$, i.e., for $d=1$, the solutions of the set of equations (2) can be written in a simple analytical form, because then the system reduces to the quadratic equation for $G^{+}(z)$ or $G^{-}(z)$. The formulas for $G^{\pm}(z)$ at the ground state have the following form:

[TABLE]

where $A=2V(1+d_{1})-U$. Their analyses are reported in Ref. Lemański et al. (2017). However, finding solutions of (2) at arbitrary temperature $T>0$ is equivalent to finding of roots of the polynomial: of the third rank in the disordered phase (when $d=0$ and $d_{1}=0$, $G^{+}=G^{-}$, and it is independent on $V$) Lemański and Ziegler (2014); Hubbard (1964); Velický et al. (1968); van Dongen and Leinung (1997) and of the fifth rank in a general case, when $0<d<1$. So, even though we have no simple analytical formulas for $G^{\pm}(z)$ when $T>0$, we are able to get very precise numerical values on $G^{+}(z)$ and $G^{-}(z)$ for any $T$ and any values of $U$ and $V$.

For $V=0$, the formulas on $G^{+}$ or $G^{-}$ and their analysis in the whole temperature region are provided in Ref. Lemański and Ziegler (2014). In a general case of $V\geq 0$ one needs to solve the following fifth rank polynomial equation on $G^{+}$ [if we know $G^{+}$ then we can find $G^{-}$ from Eq. (2)]:

[TABLE]

Coefficients $a_{0}$, $a_{1}$, $a_{2}$, $a_{3}$, $a_{4}$, and $a_{5}$ are functions of $z$, $U$, $V$, $d$ and $d_{1}$. Since the explicit expressions for these coefficients are rather lengthy, we present them into Appendix A.

As can be seen from formula (2), the entire temperature dependence of the Green’s functions, and thus also some other characteristics of the system that are expressed by them, such as, e.g., the DOS and the energy gap at the Fermi level, comes merely from the temperature dependence of parameters $d(T)$ and $d_{1}(T)$. Of course, all other thermodynamic characteristics of the system also depend on these functions, in addition to the explicit dependence on $T$. Therefore, our primary task is to determine $d(T)$ and $d_{1}(T)$.

The procedure of finding $d(T)$ and $d_{1}(T)$ is as follows. First, having determined $G^{+}$and $G^{-}$ we calculate the DOS functions $\rho^{+}$ and $\rho^{-}$ from the standard formulas:

[TABLE]

Then, for given parameters $U$, $V$, and $T$ we solve the self-consistent equation for $d_{1}$ [Eq. (4)], from which we get the dependence of $d_{1}$ as a function of $d$. Concentrations $n_{\downarrow}^{\pm}$ appearing in Eq. (4) are calculated from the expression

[TABLE]

(in our case the Fermi level is located at $\varepsilon_{F}=0$).

Next, we insert function $d_{1}(T,d)$ into the expression for free energy $F$ given by Eqs. (15), (16), and (III.1.3), and finally, we minimize $F$ over $d$, from where we get $d(T)$ and then also $d_{1}(T)$ from Eqs. (4) and (9).

III Results

The model predicts the existence of a variety of ordered and nonordered (NO) phases. The full phase diagram of the model is quite complex. In particular, the diagram exhibits five distinguishable regions in which an ordered insulator occur (COI1, COI2, and COI3 with $d_{1}>0$; AFI1 and AFI2 with $d_{1}<0$). Three distinguishable regions of ordered metal are also present (COM1 and COM2 with $d_{1}>0$; AFM with $d_{1}<0$). Note that in these ordered phases both a charge order and an antiferromagnetic order exist simultaneously, i.e., $\Delta_{Q}\neq 0$ and $m_{Q}\neq 0$ [excluding $U\rightarrow-\infty$ and $V\rightarrow+\infty$ limits ($U<2V$), where $m_{Q}\rightarrow 0$ as well as $U\rightarrow+\infty$, where $\Delta_{Q}\rightarrow 0$]. Moreover, the nonordered metal (NOM) is found to be stable is some range of the model parameters. This NOM is not the Fermi liquid phase due to a huge scattering rate that is maintained even when approaching the ground state (for details for the FKM see Refs. Si et al. (1992); van Dongen and Leinung (1997); Freericks and Zlatić (2003)). Finally, two separated ranges of the model parameters are determined, where the nonordered insulator (denoted here as NOI) occurs. In the NOM and the NOI phases $d=0$ and $d_{1}=0$ (or $\Delta_{Q}=0$ and $m_{Q}=0$, equivalently). One should be also aware of the fact that the NOI of the FKM has a quite different nature than the Mott insulator described by the HM Freericks and Zlatić (2003); Philipp et al. (2017).

In this paper all phases with $d_{1}>0$, i.e., the phases with dominant charge order ($\Delta_{Q}>m_{Q}$), are named as charge-ordered (CO) ones and all phases with $d_{1}<0$, i.e., the phases where antiferromagnetism dominates ($m_{Q}>\Delta_{Q}$), are named as antiferromagnetic (AF) ones (note that we assumed that $d\geq 0$, cf. Sec. II.2). The lower index (which can be $1$, $2$, and $3$) labels the phases with the same properties (details are given further in the text). We also introduced the denotations NOI${}_{\textrm{A}}$ and NOI${}_{\textrm{B}}$ to distinguish two region of the NOI occurrence for $U<-2$ and $U>2$, respectively. In the following sections we characterize the properties of each phase and construct the full phase diagram of the model.

III.1 Exact formulas for finite temperatures

III.1.1 Density of states

In order to be able to perform accurate calculations, it is extremely important to precisely determine the boundaries of energy bands. In our case, this is possible thanks to the polynomial form of the equation for the Green function, as it is expressed in Eq. (7). Indeed, the simultaneous solution of Eq. (7) and the equation $dW(G^{+})/dG^{+}=0$, which can be written explicitly as

[TABLE]

allows us (at fixed values $U$, $V$, $d$ and $d_{1}$) to obtain energies at which the edges of energy bands occur. Coefficients $a_{1}$, $a_{2}$, $a_{3}$, $a_{4}$, and $a_{5}$ are defined in Appendix A.

The exact formula on the density of states (DOS) at the Fermi level $\rho(\varepsilon_{F})=\left(\rho^{+}(\varepsilon_{F})+\rho^{-}(\varepsilon_{F})\right)/2$ is expressed in Eq. (11) below through parameters $U$, $V$, $d$ and $d_{1}$.

[TABLE]

where

[TABLE]

It applies throughout the entire temperature range, but only through parameters $d(T)$ and $d_{1}(T)$. Consequently, in the disordered phase with $d=0$ and $d_{1}=0$ it does not change with temperature. Notice also that, in a general case of any value of $\varepsilon$, it is not possible to do determine expression for $\rho(\varepsilon)$ analogous to (11).

In Figs. 1 and 2 DOS functions $\rho^{+}(\varepsilon)$ and $\rho^{-}(\varepsilon)$ [$\rho^{\pm}(\varepsilon)\equiv\rho^{\pm}(U,V,d,d_{1};\varepsilon)$ evaluated from Eq. (8)] in all found distinguishable ordered phases are presented for several representative values of the model parameters. The phases presented in Fig. 2 can occur only if $V>0$, whereas those shown in Fig. 1 are present on the phase diagram also for $V=0$ (cf. also Refs. Lemański and Ziegler (2014); Lemański (2016)). Note that $\rho^{+}(\varepsilon_{F})=\rho^{-}(\varepsilon_{F})=0$ in all insulating phases and the gap at the Fermi level $\Delta(\varepsilon_{F})$ is finite (but small, and thus not clearly visible in the figures). At the half-filling, it turns out that $\rho^{+}(\varepsilon)=\rho^{-}(-\varepsilon)$ and the total (resultant) DOS $\rho(\varepsilon)=(1/2)\left[\rho^{+}(\varepsilon)+\rho^{-}(\varepsilon)\right]$ is symmetric, but $\rho^{+}(\varepsilon)$ and $\rho^{-}(\varepsilon)$, when considered separately, are not symmetric. The obvious observation is that in phases with $d_{1}>0$ the $\rho^{+}$ has the larger weight than $\rho^{-}$ in the main band below the Fermi level (the lower main band) due to the fact that $n_{\downarrow}^{+}>n_{\downarrow}^{-}$. Namely, $\rho^{+}$ has lower weight in its upper main band than in the lower main band. Likewise, $\rho^{-}$ has greater weight in its upper main band than in the lower main band. In all phases with $d_{1}<0$ the relation between weights in main bands is reversed.

The behavior of the itinerant electrons in the phases collected in Fig. 1 is associated with the specific features of the subbands, which arise inside the principal energy gap between the main bands (the region indicated by the gray shadow in Fig. 1). These subbands can be present only at $T>0$ (i.e., for $d<1$), whereas at $T=0$ (i.e., for $d=1$) the DOS exhibits only two main bands separated by the principal gap at the Fermi level (cf. Fig. 1 of Ref. Lemański et al. (2017)). In particular, the characteristic feature of the COI1 which distinguishes it from the COI2 is that $\rho^{+}$ ($\rho^{-}$) has lower (higher) weight in the subband lying below $\varepsilon_{F}$ than in the subband located above $\varepsilon_{F}$. In contrast, in the COI2 $\rho^{+}$ ($\rho^{-}$) has lower (higher) weight in the subband above $\varepsilon_{F}$ than in the subband situated below $\varepsilon_{F}$. The relation between weights in the subbands of $\rho^{+}$ (or $\rho_{-}$) in the AFI1 is the same as in the COI2, whereas those in the AFI2 are such as in the COI1. One can also say that in the COI1 and the AFI1 the fillings of the main bands and of the subbands are inverted for both $\rho^{+}$ and $\rho^{-}$, i.e., the weights in $\rho^{+}$ ($\rho^{-}$) of the subbands below and above the Fermi level are in the opposite relation than the weights in the main lower and upper bands of $\rho^{+}$ ($\rho^{-}$).

The metallic behavior of the COM1 and the AFM is a result of merging of the subbands near the Fermi level, as it is visible in Fig. 1, but the precursors of the subbands can be still visible in the DOSs of these phases (at least in the regions of their occurrence in the neighborhood of boundaries with the COIs or the AFIs, respectively, on the phase diagram). In these metallic phases the continuous change from a case of inverted subbands (in the COI1 and the AFI1, in the sense discussed above) to a case of “noninverted” subbands (in the COI2 and the AFI2) occurs. In higher temperatures, the central subband with nonzero weight at $\varepsilon_{F}$ can merge with the main bands (cf. also Fig. 2 of Ref. Lemański and Ziegler (2014)).

The DOSs in the COI3 and the COM2 are shown in Fig. 2. In contrast to the DOSs discussed above, the structure of the DOS in the COI3 consists only of two main bands and does not exhibit any subbands (also at $T>0$). The metallic behavior of the COM2 is associated with closing of the gap between the main bands.

To be rigorous, one should also add that for large $|U|$ in the COI1 and the AFI1 two additional subbands can appear in the DOS, one below the lower main band and one above the upper main band. Its precursors are visible in Fig. 1 (cf. also Fig. 3 of Ref. Lemański and Ziegler (2014)). In Fig. 1 of Ref. Lemański and Ziegler (2014) the DOSs in the non-ordered phases, i.e., in the NOM and in the NOIs, are also presented.

III.1.2 Energy gap at the Fermi level

It appears, that the energy gap at the Fermi level $\Delta(\varepsilon_{F})$ is a continuous function at $T=0$ only within the interval $0\leq U\leq 2V$, but it is discontinuous both for $U<0$ and for $2V<U$. The continuity of $\Delta(\varepsilon_{F})$ at $T=0$ for $0\leq U\leq 2V$ is due to the fact, that no subbands are formed inside the principal energy gap when the temperature is raised above $T=0$ (the COI3, cf. Sec. III.1.1). Indeed, at $T=0$ one has $\Delta(\varepsilon_{F})_{T=0}=|2V(1+d_{1})-U|$ Lemański et al. (2017), but in the limit of $T\rightarrow 0^{+}$ we get the following formulas:

[TABLE]

And for the special case of $U=2V$, when two ordered phases with two different energy gaps coexist at $T=0$ (a point of a discontinuous transition Lemański et al. (2017)), are two solutions:

[TABLE]

In the above formulas [i.e., Eqs. (13) and (14)] parameter $d_{1}$ is taken at $T=0$. The behavior of $\Delta(\varepsilon_{F})$ for a few representative values of $V$ is illustrated in Fig. 3(a). Figure 3(b) shows difference $D_{\Delta}=\Delta(\varepsilon_{F})_{T\rightarrow 0^{+}}-\Delta(\varepsilon_{F})_{T=0}$ as a function of $U$ for the same values of $V$.

It turns out that if $U$ lies outside of the region $0\leq U\leq 2V$, then $\Delta(\varepsilon_{F})_{T\rightarrow 0^{+}}$ is always smaller than $\Delta(\varepsilon_{F})_{T=0}$ and it attains zero at the quasi-quantum critical points $U_{cr}^{\pm}$. There are two such points for $V\lesssim 0.54$, one positive ($U_{cr}^{+}$) and another negative ($U_{cr}^{-}$), but there is only one negative $U^{-}_{cr}$ for $V\gtrsim 0.54$ (cf. also Ref. Lemański et al. (2017)). At the first order phase transition point at $T=0$, i.e., for $U=2V$, $\Delta(\varepsilon_{F})_{T\rightarrow 0^{+}}$ exhibits discontinuous jump, but $\Delta(\varepsilon_{F})_{T=0}$ is continuous. For large positive $U$ (i.e., $U>U_{cr}^{+}$ if $V\lesssim 0.54$ and $U>2V$ if $V\gtrsim 0.54$, in the AFI1) or large negative $U$ ($U<U^{-}_{cr}$, in the COI1), $\Delta(\varepsilon_{F})_{T\rightarrow 0^{+}}$ as a function of $U$ is a monotonic function of $U$ increasing with $|U|$. Then difference $D_{\Delta}$ is obviously zero for $0\leq U\leq 2V$ (in the COI3) and it exhibits discontinuity at $U=2V$. It is a monotonous function of $U$ for $U<0$ and $U>2V$ with a change of slope at $U_{cr}^{\pm}$ (if $U_{cr}^{+}$ exists).

In the general case of finite $T$ it is possible to precisely determine the energy gap based on the knowledge of the edges of energy bands, but before that we need to determine the values of $d(T)$ and $d_{1}(T)$ numerically.

III.1.3 Free energy

Total free energy $F$ of the system per site is given by the following expression

[TABLE]

where

[TABLE]

and

[TABLE]

where $\rho(U,V,d,d_{1};\varepsilon)=\left(\rho^{+}(\varepsilon)+\rho^{-}(\varepsilon)\right)/2$ and $\rho^{\pm}(\varepsilon)\equiv\rho^{\pm}(U,V,d,d_{1};\varepsilon)$ is expressed by Eq. (8).

Free energy $F$ calculated using formulas (15), (16) and (III.1.3) has a fundamental meaning here, because only after its minimization with respect to order parameter $d$ we obtain $d(T)$ and $d_{1}(T)$ [from Eqs. (4) and (9)], that enter into all physical characteristics of the system. And on the basis of such determined free energy we constructed finite temperature phase diagrams which we present in the next section. In these diagrams, in the phases marked as insulators (conductors) the condition $\rho_{F}=0$ ($\rho_{F}\neq 0$) is met. Note that the term ’insulator’ is used to characterize the DOS with a gap at the Fermi level. It is clear that, strictly speaking, such a system would be insulating only at $T=0$.

We still need to mention that we derived formula (16) by generalizing the expression given for the FKM in Ref. Freericks et al. (1999). Then we have verified that this is equivalent to the expression presented by Van Dongen in van Dongen (1992), but has a simpler form than that in van Dongen (1992) and allow for greater precision of calculations in the whole range of the model parameters and of temperature, which is especially important near $T=0$. This is a consequence of the fact that in the formulation of this work one does not do summation over Matsubara frequencies. Such a summation is difficult to perform numerically because the tails of the Green’s functions are vanishing slowly with the frequencies. Certainly, both approaches are formally equivalent, as it was shown in Ref. Shvaika and Freericks (2003), but in practice the method we use in this work enables us to calculate all relevant physical quantities with any precision and at any temperature, what is not the case when the summations over Matsubara frequencies are performed (cf. Refs. Hassan and Krishnamurthy (2007); Lemański and Ziegler (2014) for the FKM).

III.2 Phase diagrams

In this subsection we will focus on the evolution of the phase diagram of the model with increasing intersite interaction $V$. On the diagram there are a variety of transitions, thus we discuss firstly an evolution of the order-disorder transitions (Sec. III.2.1). Next, we focus on transitions between metallic and insulating ordered phases at lower temperatures (Secs. III.2.2 and III.2.3). Finally, some dependencies of important quantities at finite temperatures are presented (Sec. III.2.4).

At the beginning let us briefly discuss the diagram for model with $V=0$. The detailed study of the phase diagram of the FKM is contained in Refs. Lemański and Ziegler (2014); Lemański (2016). For $V=0$ the phase diagram is symmetric with respect to $U=0$. However, for $U<0$ quantities $d$ and $d_{1}$ have the same signs (and the charge order dominates over antiferromagnetism), whereas for $U>0$ they have opposite signs (and the antiferromagnetic order is dominant; note that we assumed that $d\geq 0$). To be precise, for $V=0$ the solution with $d$ and $d_{1}$ for $U$ corresponds to the solution with $d$ and $-d_{1}$ occurring for $-U$.

For $V=0$ a continuous order-disorder transitions between the ordered and the nonordered phases occurs when $T$ is raised. At temperatures above the transition, two nonordered phases can exist: the NOM for $|U|<2$ and the nonordered insulator for $|U|>2$ (NOI${}_{\textrm{A}}$ for $U<-2$ and NOI${}_{\textrm{B}}$ for $U>2$). With an increase of $|U|$, the continuous NOM-NOI${}_{\textrm{A}}$ (NOM-NOI${}_{\textrm{B}}$) transformation occurs at $U=-2$ ($U=2$), respectively and it is not dependent on temperature. To be precise, this is the type of metal-insulator transition predicted by Mott (e.g., Refs. Mott (1990); Montorsi (1992); Gebhard (1997)), but unlike Mott’s prediction that the transition would be generically discontinuous Müller-Hartmann (1989b); van Dongen and Leinung (1997); Freericks and Zlatić (2003). The $U$-dependence of the order-disorder transition temperature $T_{OD}$ is shown in Fig. 4 (the first solid line from the bottom is for $V=0$). $T_{OD}$ is a nonmonotonous function of $U$. It increases with increasing $|U|$ starting from zero at $U=0$ to reach its maximal value of $k_{B}T_{OD}\approx 0.113$ at $|U|\approx 2.61$. With further increase of $|U|$ it decreases to zero for $U\rightarrow+\infty$.

At low temperatures (below $T_{OD}$ line) the model exhibits six long-range-ordered phase, in which the orders (charge and antiferromagnetic) can coexist with the metallic or insulating behavior [two ordered metals: COM1 and AFM (for $U<0$ and $U>0$, respectively); or four charge-ordered insulators: COI1 and COI2 (for $U<0$) as well as AFI1 and AFI2 (for $U>0$)]. The COM1 (AFM) can exists only for $0<|U|<2$ and $T>0$. The regions of their occurrence divides the regions of the ordered insulator occurrence into two separated areas [the COI2 (AFI2) for $0<|U|<\sqrt{2}$ and the COI1 (AFI1) for $|U|>\sqrt{2}$]. The COI1 and the COI2 for $U<0$ (the AFI1 and the AFI2 for $U>0$) are differentiated by behavior of a capacity of the main lower band and a capacity of the lower subband, lying inside the main energy gap just below the Fermi level (cf. Ref. Lemański (2016)). At $U^{-}_{cr}(V=0)=-\sqrt{2}$ ($U^{+}_{cr}(V=0)=\sqrt{2}$) the COM1 (AFM, respectively) exist at any infinitesimally small (but finite) $T>0$. There are so called quasi-critical points at $T=0$ Lemański et al. (2017).

For $V\neq 0$ the phase diagram loses its symmetry with respect to $U=0$ and it needs to be discussed for negative as well as positive onsite interactions. Nevertheless, for $V\neq 0$ at temperatures above the order-disorder transition, one can distinguish three regions of the nonordered phases: two nonordered insulators: NOI${}_{\textrm{A}}$ for $U<-2$ and NOI${}_{\textrm{B}}$ for $U>2$, as well as the NOM for $-2<U<2$. It turns out that the boundaries between nonordered phases do not depend either on $T$, or on $V$. This is a consequence of the fact that $G^{\pm}(z)$ are dependent on $V$ only through the term $v=V(d+d_{1})$ [cf. Eq. (2)] and $d=d_{1}=0$ in the NOM, the NOI${}_{\textrm{A}}$, and the NOI${}_{\textrm{B}}$.

At the phase diagram for $V>0$, two new regions appear for $U\geq 0$ (cf. Figs. 5, 6, and 7). One is that of the COI3 extending from the ground state (in the range $0\leq U<2V$ at $T=0$). Another one is a region of the COM2, which appears for $0<U<U_{c}$ at finite temperatures between regions of the COI1 and NOM and is separated from that of the COM1. $U_{c}$ is the value of $U$ interaction for the bicritical point, where two second-order lines: COM2-NOM and AFM-NOM and one first-order line: COM2-AFM merge together (it is explicitly denoted in Figs. 6, and 7). Notice that $U_{c}<2V$. The bicritical point is found for intersite repulsion smaller than $V\approx 0.7$. The detailed discussion of the evolution of the phase diagram with increasing $V$ is contained below.

III.2.1 Order-disorder transition

Figure 4 presents temperatures $T_{OD}$ of the order-disorder transition for several values of $V$ obtained within an assumption that the transition is continuous. The diagram is determined by comparing energies of the ordered phases (with $d=0.04$ and both signs of $d_{1}$ determined self-consistently) and the NO phase. Decreasing $T$ or $U$ (starting from the region of nonordered phase) with a small step the first point, when the free energy of the ordered phase is smaller, determines the continuous boundary. In Fig. 4 the line of points is also indicated, where parameter $d_{1}$ in the ordered phase changes its sign ($d_{1}>0$ on the left and $d_{1}<0$ on the right). As we will show below for large $V$ (larger than $V\approx 0.7$), the order-disorder transition is also first-order one in some range of intermediate values of $U$ (temperatures of such a transition are not shown in Fig. 4). Notice that for $V\gtrsim 1.5$ temperature $T_{OD}(U)$ for the continuous transition exhibits reentrant behavior (e.g., the boundary for $V=1.6$ is $Z$-shaped near $U\approx 1.9$, it is slightly visible in Fig. 4, cf. also Fig. 9).

For small $V\neq 0$ the order-disorder transition is indeed continuous one for any $U$, but the boundary line is no longer symmetric with respect to $U=0$, however, it retains its two maxima for both signs of $U$ (Figs. 4 and 5). The local minimum still exists between these two maxima, but it is located at $U>2V$. In the range $0<|U|<2$ the transition occurs between two metallic phases: for $-2<U<0$ between the COM1 and the NOM, for $0<U<U_{c}$ between the COM2 and the NOM, and for $U_{c}<U<2$ between the AFM and the NOM (Figs. 5 and 6). For larger $|U|$, i.e., for $|U|>2$, the transition line separates the regions of insulating phases: for $U<-2$ — the COI1 and the NOI${}_{\textrm{A}}$, whereas for $U>2$ — the AFI1 and the NOI${}_{\textrm{B}}$ (e.g., Fig. 5)). At the single point $U=0$ the transition is directly from the COI3 to the NOM (in the neighborhood of the transition point the COM1 and the COM2 are also stable).

For larger values of $V$, the maximum of transition temperature $T_{OD}(U)$ for $U>0$ disappears, but the other one for $U<0$ is present for any $V$ (Figs. 4 and 7). The disappearance of the maximum of $T_{OD}(U)$ for repulsive $U$ occurs at $V\approx 0.6$. In addition, the transition in some range of $U$ changes its order into first one (for the intersite repulsion larger than $V\approx 0.7$, Fig. 8) with the transition temperature higher than $T_{OD}$ found with an assumption of the continuous transition. This first-order transition is directly from the COI3 to the NOM, without passing through the region of the COM2. With further increase of $V$ (for $V$ larger than $V\approx 1.05$) also discontinuous COI3–NOI${}_{\textrm{B}}$ transition appears with simultaneous disappearance of the region of the AFM (Fig. 9). When $V\rightarrow+\infty$, the results approaches to the CO-NO transition line for the atomic limit of the model ($n=1$, $D\rightarrow+\infty$) with an occurrence of the first-order transition for $4/3\ln(2)<U/V<2$ and $2/3>k_{B}T/V>0$ (the region of the AFI1 vanishes in this limit) Micnas et al. (1984); Kapcia and Robaszkiewicz (2016, 2012). For any finite $V$, the second-order transition for $U>0$ with increasing $T$ from the ordered phases with $d_{1}>0$ is only from the COM2 and can be to the NOM (as shown for $V=1.6$) or even to the NOI${}_{\textrm{B}}$ (for larger $V$, not shown in the figures). Note also that the reentrant behavior of $T_{OD}$ for the second-order transition has been also found in the atomic limit of the model Micnas et al. (1984); Kapcia and Robaszkiewicz (2012). Moreover, the second-order transition temperature $T_{OD}$, which is lower than the true $T_{OD}$ of the first-order transition, was identified as the boundary of the NO phase metastability inside the CO phase region Kapcia and Robaszkiewicz (2012); Kapcia (2014).

Finally, let us discuss the limit of large $|U|$, i.e., $U\rightarrow\mp\infty$. The results for $U\rightarrow-\infty$, i.e., $k_{B}T_{OD}\rightarrow 2V$ (where the COI1-NOI${}_{\textrm{A}}$ transition is present) are in an agreement with the results for the atomic limit of the model Micnas et al. (1984); Kapcia and Robaszkiewicz (2016, 2012). Notice that the order-disorder line exhibits local maximum for some $U<0$ and for any finite $V$. The local maximum moves to $U\rightarrow-\infty$ if $V\rightarrow+\infty$. In this limit the transition temperature decreases monotonously with increasing of $U$. The local maximum for $U>0$ exists only for $V$ smaller than $V\approx 0.6$. For large positive $U$ the transition temperature behaves as $k_{B}T_{OD}\approx 1/(2U)$ and for $U\rightarrow+\infty$ it decreases to $T_{OD}\rightarrow 0$ (the AFI1-NOI${}_{\textrm{B}}$ transition). It is a result of an equivalence of the EFKM with the Ising model for large $U>0$ with $J\approx t^{2}/(2U)$ in this limit (here $t$ is explicitly given) and $k_{B}T_{OD}\approx 1/(2U)$ ($k_{B}T/J=1$ for the Ising model).

III.2.2 Discontinuous transitions between ordered phases

In this section we focus on the first-order (discontinuous) transitions between various ordered phase associated with discontinuous change of $d$. In the model investigated in this work such transitions occur only between phases with different signs of $d_{1}$, i.e, only between the CO and AF phases. However, the criterion for the transition is equality of the free energies of both phases (cf. Sec. III.2.4). As it was mentioned previously, for $V=0$ the regions of the COI2 and the AFI2 are connected only by single point at $U=0$ and $T=0$. For $V\neq 0$ three new discontinuous transitions appear on the phase diagram: COI3-AFI2, COI3-AFM, COM2-AFM. They are mentioned in the sequence consistent with an increase of temperature (cf. Figs. 5 and 6; in Fig. 5 this behavior is slightly visible). Notice that the discussed here first-order transitions can be found only for $U_{c}<U<U_{cr}^{+}$ with an increase of temperature. The first-order boundary on the phase diagram is a decreasing function of $U$. With increasing $V$ the COI3-AFI2 line shrinks and at $V\approx 0.54$ it totally vanishes (accompanied with the disappearance of the AFI2 region, Fig. 7). The further increase of $V$ results in an appearance of COI3-AFI1 transition. Finally, the first-order line evolves into direct phase transitions from COI3 to the nonordered phases (the NOM and the NOI${}_{\textrm{B}}$) as described before in Sec. III.2.1 (Figs. 8 and 9). Notice that the discontinuous COI3-AFI1 boundary line is still present on the diagram and it totally vanishes only in $V\rightarrow+\infty$ limit.

III.2.3 Continuous metal-insulator transformations

On the phase diagram of the model also several continuous metal-insulator transformations between the ordered phases were found, namely,

• (i)

for $U<0$: COI1-COM1 and COI2-COM1;

• (ii)

for $U>0$: AFI1-AFM, AFI2-AFM, and COI3-COM2.

As one can notice, they occur only between phases with the same sign of $d_{1}$. However, the physics of the COI3-COM2 transformation is different than the other ones. Namely, it is associated with a disappearance of the energy gap between main (lower and upper) bands, whereas the rest of the transformations are connected to closing the gap between subbands inside the principal energy gap (cf. Sec. III.1.1).

Here, we use the term “transformation” for the change between an insulator [with $\rho(\varepsilon_{F})=0$] and conductor [where $\rho(\varepsilon_{F})>0$], because this change is not accompanied by a discontinuity of the first or second derivative of free energy, and according to the convention it cannot be classified as a phase transition of the first or second kind.

For $U<0$ the evolution of the boundaries with $V$ is not very complicated. The region of the COM1 extends from the ground state (precisely it occurs at any $T_{OD}>T>0$ for $U=U^{-}_{cr}$) due to the fact of existence of the quasi-critical point at $T=0$ and $U=U^{-}_{cr}(V)<0$ for any $V>0$ Lemański et al. (2017). For small $V$, the temperature of the COI1-COM1 transformation decreases with $U$ from its maximal value at $U=-2$ (which is equal to $T_{OD}(U=-2)$) to zero at $U_{cr}^{-}$ (Fig. 5). With an increase of $V$ the COI1-COM1 boundary changes its slope at $V\approx 0.268$ (in that point $U^{-}_{cr}(V\approx 0.268)=-2$), and for larger $V$ the temperature of the transformation is an increasing function of $U$ (Fig. 10). Nevertheless, the area of the COM1 separates the regions of the COI1 and the COI2 for any $V$. The COI1-COM1 transformation temperature is a nonmonotonous function of $U$, increasing from zero at $U_{cr}^{-}$ to its maximal value and next it decreases to the value equal $T_{OD}(U=0)$.

The situation for $U>0$ is more complex. For any $V\neq 0$ the region of the COM2 phase appears for $U_{c}>U>0$ and $T>0$ (cf., e.g., Fig. 5). The COI3-COM2 boundary line decreases with $U$. The quasi-critical point at $T=0$ exist for $V\lesssim 0.54$ Lemański et al. (2017) and thus the AFM can occur at any $T>0$. The AFI2-AFM boundary is a nonmonotonous function of $U$. The region of the AFI2 disappears at $V\approx 0.54$ (cf. Figs. 6 and 7) as a result of that the first-order COI3-AFI2 transition line and the line of quasi-critical points merge at $T=0$. The AFI1-AFM boundary is increasing function of $U$. The region of the AFM moves towards higher temperatures and shrinks with increasing $V$. With further increase of $V$ the region of the AFM occurrence disappears (cf. Figs. 8 and 9).

For the discussion of the phase diagram to be complete one needs to mention also the COI2-COI3 transformation. It is a continuous one between two insulating phases with positive $d_{1}$. It occurs at $U=0$ for any $V>0$ and its location is not dependent on $k_{B}T$ (for $V=0$ the boundary is reduced to a point at $T=0$). It is associated with vanishing of the subband structure inside the main energy gap of the charge-ordered insulator, as it was described in Sec. III.1.1 [the lower (upper) main band and the lower (upper, respectively) subband merge together at $U=0$]. At the transformation, no discontinuities of $d$ and $d_{1}$, as well as $\rho(\varepsilon_{F})$ and $\Delta(\varepsilon_{F})$, are found. Similarly as for previously mentioned in this section metal-insulator transformations, at COI2-COI3 boundary the derivatives of $F$ do also not exhibit standard behavior expected from second-order transitions.

Please note that all transformations discussed in this section occur between phases with the same order and are associated with continuous changes of $d$ and $d_{1}$ parameters. They are defined by the behavior of quantity $\rho(\varepsilon_{F})$, but at the transformations no discontinuities $\rho(\varepsilon_{F})$ and $\Delta(\varepsilon_{F})$ are found. Thus, they are not phase transition in the usual sense (e.g., there is no kink of free energy or entropy at the transition point) similarly to the metal-insulator transition in the FKM Lemański and Ziegler (2014). Nevertheless, one should underline that order-disorder transitions (both continuous and discontinuous) found in Sec. III.2.1 as well as discontinuous transitions between various ordered phases collected in Sec. III.2.2 are conventional phase transitions. The details are included in Sec. III.2.4, where dependencies of thermodynamic parameters are shown for a few exemplary boundaries discussed previously.

III.2.4 Changes of thermodynamic quantities at the phase boundaries

Let us start this section from revisiting of the FKM [i.e., Eq. (1) with $V=0$]. As we wrote at the beginning of Sec. III, in all ordered phases both charge-order and antiferromagnetism coexist. It is not an obvious fact, because the previous works on the EFM van Dongen and Vollhardt (1990); Freericks et al. (1999); Chen et al. (2003); Maśka and Czajka (2006); Hassan and Krishnamurthy (2007); Lemański and Ziegler (2014) concentrated on the analysis of the behavior of a difference of concentration of immobile particles in both sublattices, i.e., parameter $d$ in this work [cf. Eq. (3)] and parameter $d_{1}$ was not determined. For $V=0$ free energy $F$ does not depend on parameter $d_{1}$, but it can be calculated from Eqs. (4) and (9). Parameters $d$ and $d_{1}$ are presented as a function of temperature for $U=1.0$ on the left panel of Fig. 11 (cf. also Figs. 5 and 6 of Ref. Lemański and Ziegler (2014)). This corresponds to AFI2-AFM-NOM sequence of continuous transitions. The right panel of Fig. 11 presents the parameters as a function of $U$ for fixed $k_{B}T=0.06$ corresponding to NOM-AFM-AFI1-NOI${}_{\textrm{B}}$ sequence (cf. also Fig. 7 of Ref. Lemański and Ziegler (2014)). Using relations (5), one can calculate charge polarization $\Delta_{Q}$ and staggered magnetization $m_{Q}$, which are also plotted in Fig. 11. All four discussed parameters change continuously at the phase boundaries and they vanish to zero at the order-disorder transition points.

It is clearly seen that, even for $V=0$, both $\Delta_{Q}$ and $m_{Q}$ are nonzero (what is equivalent to $d\neq|d_{1}|$ and $d,d_{1}\neq 0$). This implies that both charge order and antiferromagnetism exist simultaneously in the ordered phases of the FKM. Please also note that for $U>0$ and $V=0$ parameter $d_{1}<0$ and thus $\Delta_{Q}<m_{Q}$ (the antiferromagnetic order dominates). For $V=0$ the model exhibits a symmetry and one can obtain the results for $-U$ simply by transforming $d_{1}\rightarrow-d_{1}$, $\Delta_{Q}\rightarrow m_{Q}$, and $m_{Q}\rightarrow\Delta_{Q}$ together with changing $U\rightarrow-U$ ($d$ does not change under the transformation). For $U<0$ parameter $d_{1}>0$ and thus $\Delta_{Q}>m_{Q}$ (the charge-order dominates over antiferromagnetism).

For completing our discussion on the phase diagram of the EFKM we present a few quantities such as parameters $d$ and $d_{1}$, DOS at the Fermi level $\rho(\varepsilon_{F})$, and energy gap at Fermi level $\Delta(\varepsilon_{F})$ as well as free energy $F$ per site and entropy $S=-\tfrac{\partial F}{\partial T}$ per site at finite temperatures for representative sets of the model parameters.

In Fig. 12 the temperature dependencies of them are presented for $U=1.0$ and $U=2.5$ (and for $V=0.1$). It is clearly seen that the both order-disorder transitions: AFM-NOM and AFI1-NOI${}_{\textrm{B}}$ are the second-order transitions with the standard behavior of order parameter $d$, which continuously decreases with increasing temperature and goes to zero at the transition point. The behavior of parameter $|d_{1}|$ is similar (notice that $d_{1}<0$, i.e., $\Delta_{Q}<m_{Q}$, in the phases discussed). It should be also noted that, obviously, both $\Delta_{Q}$ and $m_{Q}$ [as a linear combination of $d$ and $d_{1}$, Eq. (5)] also vanish continuously at the transition temperature. Gap $\Delta(\varepsilon_{F})$ in the both insulating phases (i.e., AFI2 and AFI1) decreases with an increase of $k_{B}T$. One can distinguish two regions with different slopes, but the boundary between them cannot be undoubtedly determined. It is a kind of smooth crossover between these different behaviors. At the AFI2-AFM transformation $\Delta(\varepsilon_{F})$ goes continuously to zero, whereas at the AFI1-NOI${}_{\textrm{B}}$ transition it gets the value of the gap in the NOI${}_{\textrm{B}}$. $\rho(\varepsilon_{F})$ increases with $k_{B}T$ inside the AFM from zero at the AFI2-AFM boundary to the value of the $\rho(\varepsilon_{F})$ in the NOM. At NOM it does not change with further increase of temperature. Notice that at the AFI2-AFM transformation all discussed parameters change continuously, thus this boundary indeed is a continuous one (just like another metal-insulator) between phases with the same signs of $d_{1}$, cf. Sec. III.2.3).

Moreover, from the bottom panels of Fig. 12 it is clearly seen that at the transition temperatures of the both order-disorder transitions $F$ and $S$ are continuous, but the different slope of $S(T)$ in ordered and nonordered phases is visible ($S$ has a kink at the transition point). It is associated with the discontinuity of specific heat ($c=\tfrac{1}{T}\tfrac{\partial S}{\partial T}$) as it should be for the second-order transition with continuous change of the order parameter (a discontinuity of the second derivative of free energy). In contrast, at the AFI2-AFM transformation there is no such behavior and thus one can conclude that the continuous metal-insulator transformation is not a phase transition in the usual sense [$d$ (as well as $d_{1}$), free energy $F$, entropy $S$, and specific heat are continuous there] (cf. also Ref. Lemański and Ziegler (2014)).

Figure 13 presents the behavior of mentioned quantities as a function of onsite interaction $U$ in the neighborhood of two discontinuous transitions: COI3-AFM and COI3-NOM. One can notice characteristic features of the first-order transitions. Namely, in their neighborhood there is a region of a coexistence of the two phases (indicated in the figure by the gray shadow), where both solutions can coexist. In the coexistence region one of them has a higher free energy (the metastable phase) than the other (the stable phase) (cf. the bottom panels of Fig. 13, where the free energies $F$ of both solutions are presented). The value of $U$ at which energies of the both phases are the same is the transition point (denoted by vertical dashed-dotted lines in the figure). Moreover, in these panels also the different slope of $F(U)$ in each phase from both sides of the boundary is also visible. It indicates that first derivative $\tfrac{\partial F}{\partial U}$ is discontinuous at the transition point as expected for a first-order transition. These two first-order transitions are associated with a discontinuous change of all discussed quantities at the transition point, in particular, of parameters $d$ and $d_{1}$ (cf. the top panels of Fig. 2). Note also that in Fig. 13 the dependencies of mentioned quantities in the metastable phases (i.e., in the phases, which energies are not the lowest one for given set of model parameters) are also shown (dotted lines).

The COI3-AFM transition is an example of the transition between two ordered phases, at which $d_{1}$ changes not only discontinuously, but it even changes its sign (i.e., $d_{1}>0$ in the COI3 and $d_{1}<0$ in the AFM; $V=0.6$, $k_{B}T=0.085$, cf. Sec. III.2.2). At the transition the discontinuous change of parameter $d$ is not so high but it is also clearly visible (note that it decreases with increasing $U$). Obviously, $\Delta_{Q}$ and $m_{Q}$ also exhibit discontinuous jump at the transition point. $\Delta(\varepsilon_{F})$ is decreasing functions of $U$ in the COI3, whereas $\rho(\varepsilon_{F})$ slightly decreases with increasing $U$ in the AFM.

Finally, the discontinuous COI3-NOM transition is an example of the order-disorder transition, which can occur for $V$ larger than $V\approx 0.7$. Parameters $d$ and $d_{1}$ as well as $\Delta(\varepsilon_{F})$ decrease with increasing $U$ in the COI3 and, at the transition point, they exhibit abrupt drop to zero in the NOM. It is worth to notice that $\rho(\varepsilon_{F})$ in the NOM is decreasing function of $U$ due to the enhanced correlation between electrons.

Note that a coexistence of both charge-order and antiferromagnetic order in ordered phases (for finite $U$) is a specific feature of systems, where the hopping integrals for electron species are different (cf. also Refs. Dao et al. (2012); Winograd et al. (2011, 2012); Sekania et al. (2017)). In the case of the EFKM considered here, spin-$\downarrow$ electrons can move so they are “more delocalized” than the spin-$\uparrow$ electrons (ions). Thus, if there is an order (inhomogeneous distribution) of electrons in the system, the difference of spin-$\uparrow$ electron concentrations in the sublattices (parameter $d$) needs to be larger than those of mobile electrons with spin-$\downarrow$ (parameter $|d_{1}|$). The argument is especially justified in the case of $V=0$, i.e., in the case of the FKM. This is totally different from the extended Hubbard model with $V>0$ at $n=1$, where a phase only with charge order ($\Delta_{Q}\neq$ and $m_{Q}=0$) or only with magnetic order ($m_{Q}>0$ and $\Delta_{Q}=$) exists, e.g., Refs. Hirsch (1984); Lin and Hirsch (1986); Lin et al. (2000). In particular, for $U>0$ the half-filled Hubbard model ($V=0$) exhibits only the antiferromagnetic order Penn (1966); Georges et al. (1996).

IV Conclusions and final remarks

In this work we investigated the extended Falicov-Kimball model [Eq. (1)] at half-filling (i.e., $n=1$ or, equivalently, $n_{\uparrow}=n_{\downarrow}=1/2$) on the Bethe lattice within an approach which captures properly the effects of the local electron correlations. Both onsite $U$ and intersite $V$ terms are treated in the consistent approach to give the rigorous results in the limit of large dimensions. In this limit the intersite term reduces to the Hartree contribution. The main achievements of the research contained in this work are as follows:

Derivation of the exact expressions (part of them are analytical) for the DOS, the energy gap at the Fermi level (for insulators) and the free energy (expressed by parameters $d$ and $d_{1}$) for the extended Falicov-Kimball model in the limit of $D\rightarrow+\infty$.

• 2)

Construction and analysis of the phase diagrams of the model (they appeared to be quite rich) obtained within the rigorous method for the whole range of interaction parameters $U$ and $V$ as well as temperature $T$.

Once again, it appears that the relatively simple model of correlated electrons, when it is solved exactly, provides quite complicated phase diagrams with many ordered phases detected on them. The differences between some of these phases are related with additional parameters, that may not be detected using approximate calculations. Perhaps this is why sometimes phase transitions with hidden order parameters are reported in scientific literature. On the other hand, the exact solutions permit to distinguish very precisely various ordered phases and to determine precisely of regions of their stability in the space of their model parameters. In our case, for example, such “hidden” parameters that enable to distinguish some ordered insulating phases are weights of subbands of the DOS located inside the principal energy gap.

It is worth emphasizing here that a thorough examination of this system in the whole range of interaction parameters $U$ and $V$ at any $T$ was possible due to obtaining precise expressions not based on the summation over the Matsubara frequencies. It allowed us to investigate, e.g., a nonanalyticity of the ground state and determine the quasi-critical points at $T=0$. Moreover, by performing integration on the real axis, the numerical results could be obtained (in practice) with any precision. It would be not possible to attain this (in reasonable computing time) if the method of summation over Matsubara frequencies was applied due to the fact that the Green’s functions vanish very slowly on the imaginary axis [$G(i\omega_{n})\sim 1/\omega_{n}$].

Our results show that taking into account even a small Coulomb repulsive force $V$ between electrons located on neighboring sites is very important, because it causes a qualitative change in some properties of the system. For example, if $0\leq U<2V$ (and $d_{1}>0$, precisely in the COI3 region), then no additional subbands within finite temperatures arise inside the main energy gap. They only arise when $U<0$ (and $d_{1}>0$) or $U>2V$ (and $d_{1}<0$). Whereas for $V=0$ such additional subbands arise for any $U\neq 0$. The phase diagrams for $V>0$ at finite temperatures becomes asymmetric with the conversion of $U\rightarrow-U$ and much more complex than that with $V=0$. In particular, for $U>0$, the quasi-critical quantum point exists only when $V<V_{cr}\approx 0.54$, while for $U<0$ it exists for any value of $V$. A very interesting effect is also the existence of the discontinuous phase transition associated with a step change of $d_{1}$ parameter from a positive to a negative value (or vice versa). Moreover, various order-disorder and metal-insulator transitions were found in the model.

Although our results are exact only in the limit of $D\rightarrow\infty$, they are also useful for finite dimensions. Indeed, there is a qualitative similarity between the DOS of the cubic $D=3$ system and the DOS of the Bethe lattice in the limit $D\rightarrow\infty$ (in both these cases, the DOS close to its band edges has the square-root-type behavior). On the other hand, some Monte Carlo calculations performed for the FKM (see, e.g., Refs. Maśka and Czajka (2006); Žonda et al. (2009)) show that in the square $D=2$ systems subbands appear inside the principal energy gap in a similar way as it was observed for the Bethe lattice in the limit $D\rightarrow\infty$.

Let us stress here that the decomposition of the system into two sublattices and consideration only of such classes of solutions work for the half-filling case only, when the number of both localized and itinerant particles are equal to one-half of the number of lattice sites. This is the rigorous result, not an approximation, for any alternate lattice Kennedy and Lieb (1986); Lieb (1986). Out of this symmetry point, e.g., in the system with a doping, one deals with either phase separation (see, e.g., Refs. Freericks et al. (1999); Freericks and Lemański (2000)) or with higher periodic or even incommensurate phases, as it was shown, for example, for 2D systems at $T=0$ in Ref. Lemański et al. (2002).

Acknowledgements.

The authors express their sincere thanks to J. K. Freericks and M. M. Maśka for useful discussions on some issues raised in this work. K.J.K. acknowledges the support from the National Science Centre (NCN, Poland) under Grant No. UMO-2017/24/C/ST3/00276.

Appendix A Coefficients of the polynomial given in Eq. (7)

Here are the coefficients $a_{0}$, $a_{1}$, $a_{2}$, $a_{3}$, $a_{4}$, $a_{5}$ given in Eq. (7) that are obtained from the transformation of the system of Eqs. (2).

[TABLE]

As one can notice, the coefficients are indeed expressed explicitly by $z$, interactions $U$ and $V$, and parameters $d$ and $d_{1}$. It allows us to determine the properties of the investigated system with very high precision.