Extended Falicov-Kimball model: Hartree-Fock vs DMFT approach

The Hamiltonian of the EFKM (cf. also Refs. van Dongen and Vollhardt (1990); van Dongen (1992); Lemański et al. (2017); Kapcia et al. (2019, 2020)) has the following form

where $\hat{c}^{{\dagger}}_{i,\downarrow}$ ($\hat{c}_{i,\downarrow}$) denotes creation(annihilation) of fermion (electron) with spin $\sigma$ ($\sigma\in\{\downarrow,\uparrow\}$) at site $i$ and $\hat{n}_{i,\sigma}=\hat{c}^{{\dagger}}_{i,\downarrow}\hat{c}_{i,\downarrow}$. $U$ and $V$ denote on-site and intersite nearest-neighbor, respectively, density-density Coulomb interactions. $\sum_{\langle i,j\rangle}$ indicates summation over the nearest-neighbor pairs. $\mu_{\sigma}$ is the site-independent chemical potential for electrons with spin $\sigma$. In this work, we consider the case of half-filling, i.e., $\mu_{\sigma}=(U+4V)/2$ for both directions of the spin. The denotation used are the same as those used in Refs. van Dongen (1992); Lemański et al. (2017); Kapcia et al. (2019, 2020).

A review of properties of the standard Falicov-Kimball model (FKM) [called also as the spin-less Falicov-Kimball model, i.e., $V=0$ case of model (II.1), particularly in the infinite dimension limit] can be found, e.g., in Refs. Falicov and Kimball (1969); Kennedy and Lieb (1986); Lieb (1986); Brandt and Mielsch (1989, 1990, 1991); Brandt and Urbanek (1992); Freericks et al. (1999); Jędrzejewski and Lemański (2001); Chen et al. (2003); Freericks and Zlatić (2003); Freericks (2006); Hassan and Krishnamurthy (2007); Lemański and Ziegler (2014); Lemański (2016); Krawczyk and Lemański (2018); Žonda et al. (2019); Astleithner et al. (2020). One should note that other extensions of the Falicov-Kimball model such as an explicit local hybridization, a level splitting, various nonlocal Coulomb interaction, correlated and extended hoppings, or a consideration of a larger number of localized states are also possible, e.g., Refs. Brydon et al. (2005); Yadav et al. (2011); Farkašovský (2015); Hamada et al. (2017); Farkašovský (2019) and extensive lists of references, which can be found in the reviews (cf. Refs. Freericks and Zlatić (2003); Gruber and Macris (1996); Jędrzejewski and Lemański (2001)).

From the historical perspective, the standard FKM appears in Hubbard’s original work Hubbard (1963) and it was analyzed as an approximation of the full Hubbard model in Ref. Hubbard (1964). Next, it was proposed for a description of transition metal oxides Falicov and Kimball (1969) as well as a model for crystallization Kennedy and Lieb (1986); Lieb (1986). Moreover, the FKM can describe some anomalous properties of rare-earth compounds with an isostructure valence-charge transition, such as Yb_1-xY_xInCu₄ or EuNi₂(Si_1-xGe_x)₂ materials Zlatić and Freericks (2001, 2003). It can also successfully describe electron Raman scattering features in, e.g., SmB₆ and FeSi in the insulting phase Freericks and Devereaux (2001a, b); Freericks et al. (2001a). The FKM can be also applied to the pressure-induced isostructural metal-insulator transition in NiI₂ Pasternak et al. (1990); Freericks and Falicov (1992); Chen et al. (1993). Other example, where the FKM can give prediction on real systems, is the field of Josephson junctions (e.g., in Ta_xN) Freericks et al. (2001b, 2002, 2003). The model was used also for an explanation of behavior of colossal magneto-resistance materials Allub and Alascio (1997); Letfulov and Freericks (2001). However, one should underline that, in reality, not only onsite interactions occurs, thus the inclusion of intersite repulsion, as it is done in the EFKM [Eq. (II.1)], could give a better insight into the physics of real materials.

In this paper, we use mainly the semi-elliptic density of states, which is the DOS of non-interacting particles on the Bethe lattice with the coordination number $z\rightarrow\infty$ Georges et al. (1996); Freericks and Zlatić (2003). In addition, we use the gaussian DOS, which is specific one for the model of tight-binding electrons on the $d$-dimensional hypercubic lattice Metzner and Vollhardt (1989); Müller-Hartmann (1989); Georges et al. (1996), and the lorentzian DOS, which can be realized with a hopping matrix involving the long-range hopping Georges et al. (1992). For a comparison we also use the rectangular DOS. The explicit forms for the used DOSs are as follows: (i) the semi-elliptic DOS: $D_{S-E}(\varepsilon)=(2\pi t^{2})^{-1}\sqrt{4t^{2}-\varepsilon^{2}}$ for $|\varepsilon|\leq 2t$ and $D_{S-E}(\varepsilon)=0$ for $|\varepsilon|>2t$; (ii) the gaussian DOS: $D_{G}(\varepsilon)=\left(t\sqrt{2\pi}\right)^{-1}\exp\left[-\varepsilon^{2}/\left(2t^{2}\right)\right]$; (iii) the lorentzian DOS: $D_{L}(\varepsilon)=t\left[\pi\left(\varepsilon^{2}+t^{2}\right)\right]^{-1}$; (iv) the rectangular DOS: $D_{R}(\varepsilon)=1/(4t)$ for $|\varepsilon|\leq 2t$ and $D_{R}(\varepsilon)=0$ for $|\varepsilon|>2t$. In all these cases, the half-bandwidth is defined as $2t$. In the rest of the paper, we take $t$ as an energy unit.

Let us consider the model on a bipartite (alternate) lattice, i.e., on the lattice which can be divided into two sublattices (denoted by $\alpha=A,B$) in such a way that all nearest-neighbors of a site from one sublattice belong to the other sublattice.

The Hamiltonian (II.1) is treated within the standard broken symmetry mean-field Hartree-Fock approach Penn (1966); Micnas et al. (1990); Imada et al. (1998); Robaszkiewicz and Bułka (1999) using the Bogoliubov transformation Bogoljubov (1958); Valatin (1958) and restricting only to Hartree terms Müller-Hartmann (1989). Namely, we use the following decoupling of two-particle operators:

where $\langle\hat{A}\rangle$ denotes the average value of the operator $\hat{A}$ (in the thermodynamic meaning). Note that this decoupling is an exact one for intersite term (i.e., $i\neq j$) in the limit of large dimension Müller-Hartmann (1989). The interaction part of the Hamiltonian (II.1) including $U$-, $V$- and $\mu$- terms at the half-filling (i.e., $\sum_{i}\langle\hat{n}_{i,\sigma}\rangle=N/2$, $N$ – the number of the lattice sites) can be written in the form:

where $C=\frac{1}{4}\left\{-U\left(1+\Delta_{\uparrow}\Delta_{\downarrow}\right)+V\left[\left(\Delta_{\uparrow}+\Delta_{\downarrow}\right)^{2}-4\right]\right\}$, $\exp{\left(\mathbf{i}\vec{Q}\cdot\vec{R}_{i}\right)}=\pm 1$ if $i\in A,B$, respectively. $\vec{Q}$ is a half of the largest reciprocal lattice vector in the first Brillouin zone, $\vec{R}_{i}$ indicates the location of $i$-th site, and $W_{\sigma}=U\Delta_{\sigma}/2-V(\Delta_{\uparrow}+\Delta_{\downarrow})$. Parameters $\Delta_{\downarrow}$ and $\Delta_{\uparrow}$ are defined as differences between average occupation of sublattices by itinerant and localized electrons, respectively (cf. Refs. van Dongen and Vollhardt (1990); van Dongen (1992); Lemański et al. (2017); Kapcia et al. (2019, 2020)). Namely,

where $n^{\alpha}_{\sigma}=\langle\hat{n}_{i,\sigma}\rangle$ for any $i\in\alpha$, where $\alpha=A,B$ denotes the sublattice index. Now, the Hamiltonian (II.1) in the HFA can be written in terms of two sublattice operators in the reduced Brillouin zone (RBZ) as $\hat{H}_{MF}=\sum_{\vec{k},\sigma}\hat{\Phi}_{\vec{k},\sigma}^{\dagger}\mathbb{H}_{\vec{k},\sigma}\hat{\Phi}_{\vec{k},\sigma}+NC$, where $\hat{\Phi}_{\vec{k},\sigma}^{\dagger}=\left(\hat{c}^{\dagger}_{A,\vec{k},\sigma},\hat{c}^{\dagger}_{B,\vec{k},\sigma}\right)$ are the Nambu spinors, $\hat{c}^{\dagger}_{\alpha,\vec{k},\sigma}$ and $\hat{c}_{\alpha,\vec{k},\sigma}$ are fermion operators defined by the discrete Fourier transformation:

Matrix $\mathbb{H}_{\vec{k},\sigma}$ has the form

with $\epsilon_{\vec{k},\uparrow}=0$ and $\epsilon_{\vec{k},\downarrow}=(t/\sqrt{z})\sum_{m}\exp{(-\mathbf{i}\vec{k}\cdot\vec{\delta}_{m})}$, $\vec{\delta}_{m}$ defines the locations of the nearest-neighbor sites in the unit cell consisting of two lattice sites and the sum is done over all nearest neighbors. All summations over momentum $\vec{k}$ are performed over the RBZ (there is $N/2$ states in the RBZ for each sublattice per spin). The eigenvalues $E_{\vec{k},\sigma}^{\pm}$ of $\mathbb{H}_{\vec{k},\sigma}$ are given by $E_{\vec{k},\sigma}^{\pm}=\pm\sqrt{\epsilon_{\vec{k},\sigma}^{2}+W_{\sigma}^{2}}$. Thus, the full mean-field Hamiltonian can be expressed in the quasi-particle excitation as $\hat{H}_{MF}^{{}^{\prime}}=\sum_{\vec{k},\sigma,r=\pm}E_{\vec{k},\sigma}^{\pm}\gamma^{\dagger}_{\vec{k},\sigma,r}\gamma_{\vec{k},\sigma,r}+NC$, where $\gamma^{\dagger}_{\vec{k},\sigma,r}$ are creation operators of the fermionic quasiparticles. The grand canonical potential $\Omega$ (per site) of the system is defined as $\Omega=-1/(N\beta)\ln\{\textrm{Tr}[\exp(-\beta\hat{H}_{MF})]\}$. One gets $\Omega=-1/(2N\beta)\sum_{\vec{k},\sigma,r=\pm}\ln\left[1+\exp\left(-\beta E_{\vec{k},\sigma}^{r}\right)\right]+C$. The free energy $F=\Omega+\sum_{i,\sigma}\mu_{\sigma}\langle\hat{n}_{i,\sigma}\rangle$ in the half-filled case is derived as $F=-1/(N\beta)\sum_{\vec{k},\sigma}\ln\left[2\cosh\left(\beta E_{\vec{k},\sigma}^{+}/2\right)\right]+C+(U+4V)/2$. Finally, one obtains the free energy of the EFKM (per site) in the following form

where $A=-U\Delta_{\uparrow}/2+V\left(\Delta_{\uparrow}+\Delta_{\downarrow}\right)$, $B=-U\Delta_{\downarrow}+2V\left(\Delta_{\uparrow}+\Delta_{\downarrow}\right)$, and $\beta=1/(k_{B}T)$. $T$ denotes temperature and $k_{B}$ is the Boltzmann constant. $D(\varepsilon)$ is the non-interacting density of states (DOS), which an explicit form is dependent on the particular lattice on which model (II.1) is considered (as described in Sec. II.1).

After some straightforward transformations, one also gets the self-consistent equations for parameters $\Delta_{\downarrow}$ and $\Delta_{\uparrow}$ in the following form:

Value $\Delta_{F}(\varepsilon)=2A$ can be also interpreted as an energy gap at the chemical potential (the Fermi level at $T=0$) for Bogoliubov quasiparticles Lemański et al. (2017).

Let us also notice, that Eqs. (4)–(5) could be also obtained from the conditions $\partial F/\partial\Delta_{\sigma}=0$ (formally, with an exclusion of $U=0$ and $U=4V$). These are, however, only the necessary conditions for the extremum value of (II.2) with respect to $\Delta_{\downarrow}$ and $\Delta_{\uparrow}$. Thus, the solutions of (4)–(5) can correspond to a local minimum, a local maximum, or a point of inflection of $F(\Delta_{\downarrow},\Delta_{\uparrow})$. In addition, the number of minima can be larger than one (particularly for small $T>0$), so it is very important to find the solution which corresponds to the global minimum of (II.2).

Note that if pair $(\Delta_{\downarrow},\Delta_{\uparrow})$ is a solution of set (4)–(5) then pair $(-\Delta_{\downarrow},-\Delta_{\uparrow})$ is also a solution of the set [Eqs. (4)–(5) do not change their forms if one does substitution $\Delta_{\downarrow}\rightarrow-\Delta_{\downarrow}$ and $\Delta_{\uparrow}\rightarrow-\Delta_{\uparrow}$]. Also from (II.2), one gets that $F(U,V,\Delta_{\downarrow},\Delta_{\uparrow})=F(U,V,-\Delta_{\downarrow},-\Delta_{\uparrow})$ at any $\beta$. Thus, we can further restrict ourselves to find the solutions of (4)–(5) only with $\Delta_{\uparrow}\geq 0$ (parameter $\Delta_{\downarrow}$ can be of any sign). These properties are connected with an equivalence of both sublattices of an alternate lattice.

These two parameters can be connected with charge polarization $n_{Q}$ and staggered magnetization $m_{Q}$ by relation: $n_{Q}=\left(\Delta_{\uparrow}+\Delta_{\downarrow}\right)$ and $m_{Q}=\left(\Delta_{\uparrow}-\Delta_{\downarrow}\right)$, which create a different, but totally equivalent, set of parameters ($n_{Q}\geq 0$ and $m_{Q}\geq 0$, because of assumed $\Delta_{\uparrow}\geq 0$ and relation $\Delta_{\uparrow}\geq\Delta_{\downarrow}$ founded) Lemański et al. (2017); Kapcia et al. (2019, 2020). These quantities define various phases occurring in the system. A solution with $\Delta_{\downarrow}=\Delta_{\uparrow}=0$ ($n_{Q}=m_{Q}=0$) corresponds to the nonordered (NO) phase. In the ordered phases, $\Delta_{\uparrow}\neq 0$ or $\Delta_{\downarrow}\neq 0$ ($n_{Q}\neq 0$ or $m_{Q}\neq 0$). We distinguish two such phases ($\Delta_{\uparrow}>0$ assumed): (i) the CO phase, where charge order dominates, i.e., $\Delta_{\downarrow}>0$ ($n_{Q}>m_{Q}$) and (ii) the AF phase, where antiferromagnetic order is dominant, i.e., $\Delta_{\downarrow}<0$ ($n_{Q}<m_{Q}$). A very special case of $\Delta_{\uparrow}>0$ and $\Delta_{\downarrow}=0$ (i.e., $n_{Q}=m_{Q}>0$) occurring for $U=2V$ is discussed in detail in Sec. III.1. Note that, in both CO and AF phases, the long-range order breaks the same translation symmetry.

Please also note that the HFA for the intersite term restricted only to the Hartree terms is an exact approach in the limit of large dimensions Müller-Hartmann (1989), but for the onsite term in a general case the DMFT needs to be used Georges et al. (1996); Imada et al. (1998). However, the HFA can work correctly in the ground state for some particular models and phases (cf. Appendix A and Ref. Lemański et al. (2017)).

The diagram of the EFKM for $T>0$ is determined by finding all solutions of the set of Eqs. (4)–(5), comparing their free energies and checking if they correspond to the local minima of free energy $F$ [Eq. (II.2)] with respect to $\Delta_{\downarrow}$ and $\Delta_{\uparrow}$ parameters. The set usually has one solution (we restricted ourselves to the case of $\Delta_{\uparrow}\geq 0$) corresponding to the stable phase (i.e., free energy $F$ has a single minimum). Only in some restricted ranges it has two solutions. In such a case, $F$ has two local minima and the minimum with lower (higher) free energy corresponds to a stable (metastable) phase. The parameters $\Delta_{\downarrow}$ and $\Delta_{\uparrow}$ characterize a phase occurring in the system.

The general structure of the finite temperature phase diagram of the model obtained within the HFA for fixed $V/t$ is not dependent on particular value of $V/t\neq 0$. The exemplary phase diagram for $V/t=1.0$ is shown in Fig. 1(a) [for other values of $V/t$ see also Fig. 6]. It consists of three regions. For large temperatures the NO phase is stable. With decreasing temperature, the continuous order-disorder transition at $T_{c}$ occurs [which coincides with a solution of (8)]. For $U<2V$, the transition is to the CO phase (i.e., the phase, where charge-order dominates over antiferromagnetic order), whereas, for $U>2V$, the low-temperature phase is the AF phase (antiferromagnetism dominates). $T_{c}$ is minimal for $U=2V$ and it is equal to $k_{B}T_{c}=V/2$. It increases with increasing $|U-2V|$. For $U=2V$ and $T<T_{c}^{*}$ [$k_{B}T_{c}^{*}/t=0.357t$ for $V/t=1.0$; $\beta_{c}^{*}=(k_{B}T_{c}^{*})^{-1}$ is a solution of equation $I_{c}(\beta_{c}^{*})=1/V$, cf. Appendix B], a discontinuous (first order) transition occurs between two different ordered phase. The discontinuous boundary starts at $T=0$ and ends at isolated critical point (labeled by $I$-point) for $T=T_{c}^{*}$ (similarly as other found, e.g., for atomic limit of the model Micnas et al. (1984); Kapcia and Robaszkiewicz (2016)). In Fig. 1, also regions of occurrence of the metastable phases (near discontinuous transition) are determined. A range of $U/t$ where both phases (i.e., the CO and AF phases) coexist is $1.710<U/t<2.290$ at the ground state and vanishes at $T=T_{c}^{*}$. For $U=2V$ and $T_{c}^{*}<T<T_{c}$, there is no transition but only a smooth crossover between the CO and AF phases occurs through the point, where both antiferromagnetic and charge order parameters are the same, and none of them is dominant (i.e., $\Delta_{Q}=m_{Q}$). Although, at $U=2V$, the parameter $\Delta_{\downarrow}=0$, but $\Delta_{\uparrow}=\Delta_{Q}=m_{Q}>0$, thus the system is still in the ordered phase at $T_{c}^{*}<T<T_{c}$ (cf. also Fig. 4 and Appendix B). Please note that, in Fig. 1(a), the order-disorder transition line (which can be continuous as well as discontinuous dependently on $U/t$) and the discontinuous CO–AF boundary obtained within the DMFT are also shown by grey lines. They are taken from Ref. Kapcia et al. (2019). The lines of metal-insulator transitions are not indicated there.

In Fig. 1(b), we present a dependence of order-disorder transition temperature $T_{c}$ as a function of $\left(U-2V\right)/t$ for different values of $V/t$. One can notice that it is an increasing function of $|U-2V|$ with the minimal value of $k_{B}T_{c}(U=2V)=V/2$ [this results is easily obtained from Eq. (8)]. With increasing $V$ the lines $k_{B}T_{c}-V/2$ (as a function of $U-2V$) line are one above the other [i.e., $T_{c}(U-2V,V=V_{1})\leq T_{c}(U-2V,V=V_{2})$ if $V_{1}<V_{2}$]. For $V/t\rightarrow+\infty$ or $U\rightarrow\pm\infty$ (which is equivalent with $t\rightarrow 0$) the critical temperatures approaches $k_{B}T_{c}\rightarrow V/2+|U-2V|/4$ (cf. Sec. II.2.2). In particular, for $U=0$, $T_{c}$ approaches the results for the atomic limit of the model: $k_{B}T_{c}(U=0)=V$ Micnas et al. (1984); Kapcia and Robaszkiewicz (2016) (obviously, for $U\neq 0$ it does not coincide with rigorous results for the $t=0$ limit). One should also stress that the first order boundary is located at $U=2V$ for any $V/t$.

Figure 2(a) shows $U$-dependence of $T_{c}$ for $V/t=1.0$ and different DOSs (listed in Sec. II). One can see that $k_{B}T_{c}(U=2V)=V/2$ independently on the DOS used for calculations (Appendix B). Moreover, there is no significant qualitative difference between all the curves, however, they are different quantitatively. Namely, temperatures $T_{c}$ obtained for gaussian $D_{G}(\varepsilon)$ are slightly higher that those obtained for semi-elliptic $D_{S-E}(\varepsilon)$. The line of $T_{c}$ for $D_{R}(\varepsilon)$ is located below the curve obtained for the semi-elliptical DOS. The lowest critical temperatures $T_{c}$ are calculated for the lorentzian DOS. This sequence of $T_{c}$ curves occurs for any $V/t$. One should notice that, for $V\rightarrow+\infty$ or $U\rightarrow\pm\infty$, one gets $k_{B}T_{c}=V/2+|V/2-U/4|$ (for any symmetric DOS this is proven analytically in Sec. II.2.2), but this behavior is not visible in Fig. 2(a), which is obtained for relatively small $V/t$ and $U/t$. Note also that, for $U=2V$, the first order boundary for temperatures smaller than $T_{c}^{*}$ (which depends on the DOS) occurs as well as the smooth crossover region for $T_{c}^{*}<T<T_{c}$ (not shown in the figure). The dependence of $T_{c}^{*}$ temperature (which is located for $U=2V$) as a function of $V/t$ and different DOSs is shown in Fig. 2(b). It increases monotonously with increasing $V/t$ from $T_{c}^{*}/V=0$ (at $V/t\rightarrow 0$) to $k_{B}T_{c}^{*}/V=k_{B}T_{c}/V\rightarrow 0.5$ (at $V/t\rightarrow\infty$) for any DOS. However, one should note that the lines $T_{c}^{*}/V$ obtained for the lorentzian and rectangular DOSs cross at $V/t\approx 0.632$ and this obtained for $D_{L}(\varepsilon)$ is above that calculated for $D_{R}(\varepsilon)$.

Let us now discuss limitation of the HFA and compare the results derived within this approach with the results obtained by the DMFT. As we said previously, the HFA gives rigorous results at the ground state (Sec. II.2.1 and Appendix A, cf. also Ref. Lemański et al. (2017)). For finite temperatures the situations is more complex. In Fig. 6, one can see the comparison of the order-disorder transition temperature obtained within the HFA and the DMFT (for various $V/t$). Notice that, for values of $V/t$ larger than approximately $0.7$, the DMFT predicts a discontinuous order-disorder transition. This behavior is not captured within the HFA (cf. also Fig. 1). Nevertheless, for small $U$ both approaches give very comparable results in finite temperatures. For $U=0$, they give exactly the same results, because in this limit both approaches are equivalent Müller-Hartmann (1989); Kapcia and Majewska-Albrzykowska (2020).

Please also note, that the HFA results do not agree with the results of the DMFT, which predicts that CO–AF boundary is a standard first order transition (for $T\neq 0$) associated with the discontinuity of $\Delta_{\uparrow}$ Kapcia et al. (2019, 2020). As it is shown above, at $T_{c}>T>T^{*}_{c}$, the HFA finds a smooth crossover between these two ordered phases (cf. Sec. III.1).

As one can expect the HFA fails in the description of metal-insulator transition. Within the HFA all ordered phases are insulators with the gap $\Delta_{F}(\varepsilon)=2A$ ($\Delta_{F}(\varepsilon)\neq 0$ if $\Delta_{\uparrow}\neq 0$), whereas the DMFT finds several metallic phases with the long-range order Lemański and Ziegler (2014); Kapcia et al. (2019, 2020).

Finally, let us discuss the dependence of $\Delta_{\uparrow}$ as a function of reduced temperature $T/T_{c}$ for really small interactions parameters. It was found that for small $U/t\rightarrow 0^{+}$ and $V/t=0$ the temperature dependence of $\Delta_{\uparrow}$ is quite unusual, namely, it resembles mean-field solution of Eq. (10), but with $T_{c}^{Is}=T_{c}/2$ van Dongen and Vollhardt (1990); van Dongen (1992); Chen et al. (2003); Krawczyk and Lemański (2018). In Fig. 7(a), we show the comparison of the DMFT results with those obtained within the HFA for small values of $U/t$ (for the Bethe lattice). As one can see the HFA values are higher than those of the DMFT, but the similarity is noticeable. In Fig. 7(b), we also presents the dependence of normalized $\Delta_{\downarrow}/\Delta_{\downarrow}(T=0)$ for $V/t=0$ and this quantity almost follows behavior of $\Delta_{\uparrow}$ ($\Delta_{\uparrow}(T=0)=1$ for any $U$ and $V$).