Extended Falicov-Kimball model at weak onsite and intersite Coulomb interactions

We analyze in detail a behavior of the order parameter in the half-filled extended Falicov-Kimball model for small Coulomb interactions (both onsite $U$ and intersite $V$). The parameter is defined as the difference of localized electron concentrations in both sublattices of the Bethe lattice (in the limit of large coordination number). Using two methods, namely, the dynamic mean field theory and the Hartree-Fock approximation, we found the ranges of $U$ and $V$ for which the anomalous temperature dependence of the order parameter, characterized by the sharp reduction near $T \approx T_{\textrm{C}}/2$, occurs ($T_{\textrm{C}}$ is the temperature of the continuous order-disorder transition). In order to quantitatively describe this anomaly, we defined a function that measures the departure of the order parameter dependence from the standard mean-field $S=1/2$ Ising-like curve. We determined the $U$-dependent critical value $V_{\textrm{C}}$ of $V$ above which the anomaly disappears. Indicators of the anomalous behavior of the parameter dependence can be also observed in the temperature dependence of the specific heat.


Introduction
Correlations in electron systems sometimes cause surprising effects [1][2][3][4][5]. One of the effects is the anomalous behavior of the order parameter (Θ) with reduced temperature Θ = / c in the system described by the Falicov-Kimball model (FKM) for small values of the local Coulomb coupling at the half-filling. The parameter (Θ) is defined as the difference of localized electron concentrations in both sublattices of the alternate lattice (cf. also the equations in section 2.2). This anomalous behavior consists in a sudden drop in the value of (Θ) near Θ = 1/2 and its low, but non-zero, value in the range 1/2 < Θ < 1, before it reaches 0 at the order-disorder transition temperature at Θ = 1.
The anomalous behavior of (Θ) was first noticed by P.G.J. van Dongen and D. Vollhardt, who solved the half-filled FKM exactly on the infinite Bethe lattice in the limit → 0 using the dynamical mean field theory (DMFT) [6,7]. Then, accurate calculations also carried out using the DMFT showed that such a surprising course of (Θ) occurs not only in the limiting case → 0 but also for small finite values (roughly for < 0.01) for the Bethe lattice as well as for the hyper-cubic lattice in the limit of large coordination number ( → +∞; or, equivalently, in the limit of large dimensions → +∞) [8,9]. Moreover, it was confirmed that, for any finite , function (Θ) exhibits the square-root character when Θ → 1 [7][8][9]. A similar anomalous behavior of (Θ) was also obtained by means of numerical Monte Carlo calculations for the FKM on the square lattice ( = 2) [10].
On the other hand, it was reported in reference [7] that in the extended FKM (EFKM), i.e., after additional introduction of repulsion between electrons located on adjacent ions with the interaction constant , the shape of the order parameter curve for > 0 and comparable in magnitude takes the standard mean-field = 1/2 Ising-like form, i.e., it is completely different from that for = 0. But again, this outcome was found only in the limiting case of → 0. Therefore, the question arose: does this anomalous course of the function (Θ) occur only at the singular point = 0 and disappears at any positive ? If so, this anomaly could not occur in real systems, because there is always, perhaps even very small, a Coulomb repulsion of electrons located on the neighboring ions.
To clarify this issue, in this work we systematically study the behavior of (Θ) function in a range of small values of and couplings. Our goal was to check if there is a finite area in the parameter space ( , ) in which the anomalous course of (Θ) occurs, and if so, to determine this area. Using two methods, namely (i) the DMFT and (ii) the broken-symmetry mean-field Hartree-Fock approximation (HFA), we are able to find -dependent critical value c of coupling, which destroys the anomalous behavior of (Θ) occurring for < c .
The work is organized as follows. In section 2 we describe the model and the methods used in our analysis. In section 3 the results for (Θ) are presented and analyzed by means of the introduced measure of the anomaly of (Θ) curve. Some exemplary temperature dependencies of the specific heat are discussed in section 4. section 5 includes conclusions and final comments.

Half-filled extended Falicov-Kimball model
The Hamiltonian of the EFKM (cf. also references [6,7,[11][12][13][14]) on a lattice in the second quantization formalism has the following form where the denotions used here are the same as those used in references [7,[11][12][13][14]. Soˆ † ,↓ (ˆ, ↓ ) denotes creation (annihilation) of a fermion (electron) with spin ( ∈ {↓, ↑}) at lattice site andˆ, =ˆ † ,↓ˆ, ↓ . and denote onsite and intersite nearest-neighbor, respectively, density-density Coulomb interactions, is the coordination number, , is the local chemical potential for electrons with spin at site . In this work, we consider the case of half-filling, i.e., , = ( + 4 )/2. , denotes the sum over nearest-neighbor pairs and the prefactors in the first and the third term in equation (1) were chosen so that they yield a finite and non-vanishing contribution to the free energy per site in the limit → +∞.
As one of the simplest models of correlated electron systems, the standard spineless FKM (with = 0), as well as some of its extensions have been intensively studied . In particular, other than displayed in (1) extended versions of the FKM include, e.g., those that take into account correlated or extended electron hopping [24,45]. Reviews of the results obtained for the FKM and its various extensions can be found, e.g., in references [4,29,[46][47][48].
In this paper, we analyze the model on the Bethe lattice (Cayley tree) in the limit of large coordination number [2,7,49]. Originally, the lattice was demonstrated as an infinite connected cycle-free graph with the vertices all having the same coordination number [49]. The statistical mechanics of various lattice models on this graph are usually exactly solvable, due to its very specific topology [50]. Thus, we use the semi-elliptic density of states (DOS) in the HFA calculations, which is the DOS of non-interacting particles on the Bethe lattice with the coordination number → ∞ [2,4]. The explicit form for this DOS function (in → ∞ limit) is S-E ( ) = (2π 2 ) −1 √ 4 2 − 2 for | | 2 and S-E ( ) = 0 for | | > 2 , so in all these cases, the half-bandwidth is equal to 2 . However, this DOS can be also treated as an approximation for the DOS for the cubic lattice ( = 3). In the rest of the paper, we take as an energy unit ( = 1).

Mean-field approaches used
In the present work we consider the ordering on the alternate lattice, i.e., the lattice composed of two interpenetrating sublattices. Thus, one defines two mean-field order parameters (averages) and 1 as 43706-2 differences between average occupation of sublattices by localized and itinerant particles, respectively (cf. references [7,[11][12][13][14]37]). Namely, = ↑ − ↑ , 1 = ↓ − ↓ , where = ˆ, for any ∈ , where = , denotes the sublattice index ( is an average occupation of sublattice by particles with spin ). Due to the equivalence of these two sublattices, we restrict ourselves to the investigation of the solutions with > 0. The parameters are associated to charge polarization Δ Q and staggered magnetization Q by the following relations: Δ Q = (1/2) ( + 1 ) and Q = (1/2)( − 1 ), respectively. They are quantities directly determined in experiments usually.
Here, we treat intersite term within the Hartree-Fock approximation, restricting ourselves to the Hartree terms only, namely,ˆˆ . It is an exact approach for the intersite interaction in the limit of large dimensions [51]. For the onsite term, we use two complementary methods.
One is the DMFT method, which exactly captures the effects of quantum dynamics due to onsite interaction in the limit of large dimensions [2][3][4][5]51]. In the case of the FKM, the DMFT solution obtained by using analytical expresions, including those for the total free energy , was presented in references [7,9,37]. Its extention to the case of the EFKM with ≠ 0 was demonstrated in references [7,[11][12][13], from which we take the equations to obtain the results presented in this work. The equations are exact in the whole range of and parameters and temperature ( → +∞ limit), and due to the fact that they do not contain summation over Matsubara frequencies (as it occurs in reference [7]), they allow one to determine all significant physical quantities practically with any precision even at very low temperatures (see the comment to this in reference [12,13]).
The second method is the HFA used to term with restriction to only Hartree terms (done similarly as for -term above, but with = ) [1,3], which is an approximate method for this model (the treatment of the onsite interaction is not rigorous). However, there are some limits in which it gives rigorous results for the EFKM (e.g., in the ground state and for = 0), cf. references [11,14]. The equations derived within the HFA for the EFKM are presentend in [11,14,52].
One should also underline that all calculations presented in this work are performed for small and extremely small interactions as well as temperatures. This requires a really high numerical precision, particularly in numerical calculations of derivatives (e.g., specific heat) and integrals [e.g., parameter defined in equation (3)]. The results obtained within DMFT have much larger numerical errors than those obtained by the HFA due to the fact that they require minimization of free energy at each temperature [12,13], whereas the HFA computations are reduced to solving a set of two nonlinear equations [14]. Moreover, because the DMFT computations are more time-consuming, every (Θ) curve (figure 1) obtained within the DMFT and the HFA consists of about 40 and 5000 points, respectively.
Several (Θ) curves obtained for = 10 −2 , = 10 −3 , = 10 −4 , and = 10 −6 and for a few values of are presented in figure 1. They were obtained with the HFA and the DMFT (solid and dashed lines, respectively; for = 10 −6 only the HFA results are presented). It is clearly seen that the anomalous behavior present for = 0 vanishes with increasing , and for above some C the behavior of (Θ) curve is almost the same as the standard mean-field = 1/2 Ising-like dependence of ( C ; Θ).  (Θ) obtained within the DMFT is closer to ( C /2; Θ), cf. figure 1(c). It implies that even for such a weak interaction, the proper inclusion of onsite quantum dynamics (as done in the DMFT) is also important (for really small ). However, it does not mean that the HFA fails for small , because both these approaches give quite comparable results at finite temperatures, at least qualitatively. Let us also note that the function (Θ) exhibits the square-root character when Θ → 1 for any | | + ≠ 0 (from both the DMFT and the HFA results), as it is expected for the mean-field theories. For a quantitative description of the deviation of ( , ; Θ) curve from the standard mean-field = 1/2 Ising-like dependence of ( C ; Θ), we introduce the following measure ( , ) of this deviation: In the investigated range of the model parameters, the integrand function is always positive. It attains its maximum value for , → 0 [it is equal to lim , →0 ( , ) = . At larger , the deviation of ( , ; Θ) curve from ( C , Θ) approaches zero: ( , ) ≈ 0. One can identify some -dependent C , where ( , ) changes from ( , = 0) to 0. For = 0, the DMFT gives the values of ( , = 0) larger than the HFA. Due to the previously mentioned fact that (Θ) curve obtained within the DMFT can lie even above ( C , Θ), the DMFT calculations give ( , = 0) < 0 for larger than C ≈ 0.5 [8], cf. also figure 2(a).
More informative quantity to determine C is ( , ) = ( , )/ ( , = 0) (i.e., normalized by its = 0 value). In figure 3(a), the / -dependence of ( , ) is shown. The crossover region of the change of from 1 to 0 is clearly noticeable. The decay of is slightly faster within the DMFT. It turns out that the dependence of ( , ) for fixed can be quite well described by the single-parameter exponential-decay function ( ; ) ≡ ( ) where ≡ C is the -dependent fitting parameter. This parameter can be treated as an approximation of the effective value of intersite interaction, at which the low-anomaly of (Θ) curve is destroyed: Figure 3(b) presents a few fits for the data of obtained for several values of . One see that, indeed, the assumed fitting function can be considered as an approximation of -dependence of ( , ). The dependence of C as a function of is shown in figure 4. Note that C is an increasing function of , but C / does not exceed 1/2. In the ground state, a transition between two ordered phases has been found for = 2 , cf. phase diagrams presented in references [11][12][13][14]. For > 2 , the antiferromagnetism dominates over the charge order ( Q > Δ Q or 1 < 0), whereas < 2 the charge order is dominant ( Q < Δ Q or 1 < 0). At finite temperatures, due to the exact solution of the model, the location of the boundary is shifted for smaller than 2 and it is discontinuous for any and 43706-6 any temperature (with discontinuous changes of and 1 parameters) [12,13]. Thus, the anomaly of (Θ) dependence is destroyed inside the ordered phase with dominant antiferromagnetic order ( Q > Δ Q ).
Obviously, a better description of the -dependence of can be obtained, e.g., by two-parameter exponential-decay function particularly, for smaller [cf. dash-dotted lines in figure 3(b)]. The fitted parameters 1 and 2 are also presented in figure 4. From the data presented in figure 4, one conclude that the relation between C and (for small → 0) can be approximated by The proportional coefficient in the relation above is estimated to be between 0.25 and 0.45 and it cannot be precisely determined due to accumulating numerical errors (we are working here in the regime of very small , and C ). The linear fit to three points obtained for the smallest gives the value of the coefficient about ≈ 0.43. One should notice that we have also performed the calculation for the EFKM for = 0 (cf. also expresions in reference [52]). Additionally, we checked that, for = 0 and = 10 −9 (the smallest possible value in our calculations), according to equation (6), there are no signs of an unusual behavior of (Θ) and its dependence almost follows ( C ; Θ) ( ≈ 2 · 10 −4 in this case). For = 0 and large , one can analytically show that the behavior of (Θ) is the same as that of ( C ; Θ).
The FKM exhibits symmetry ↔ − , thus the results for < 0 are the same as for > 0 if = 0. In particular, ( = | |, = 0; Θ) = ( = −| |, = 0; Θ). This symmetry is destroyed by > 0, but for → 0, it is obviously restored and the differences between opposite signs of are negligible. However, it turns out that (Θ) curves for > 0 are located above (Θ) dependencies obtained for < 0 [i.e., ( = | |, ; Θ) > ( = −| |, ; Θ)]. Thus, C ( = | |) > C ( = −| |) and the anomalous behavior of (Θ) curve near Θ = 1/2 is destroyed faster for > 0. In that sense, the effect of on the anomalous behavior of (Θ) in the phase with dominant charge order (i.e., where Δ Q > Q ) is smaller than that in the phase with dominant antiferromagnetic ordering occurring for > 2 discussed previously. We do not present the results for < 0 in figures due to the fact that the differences are very small (below the thickness of lines in the figures) and they decrease with a decreasing | |.
Note that, in the considered range of small and , one can also distinguish two regions, where C behaves differently [7]. Namely, for ( → 0), the physics of the model is driven by  interaction and C ∝ 2 ln(1/ ). For larger , where interaction dominates, C ∝ /2. These two behaviours are clearly indicated in figure 5. Here, we plotted B C / ( B C normalized by ), where = 1 = (2π) −1 2 ln(1/ ) or = 2 = /2. These two different regions of C behaviour are clearly visible. The first one corresponds to almost constant B C / 1 (linearly decreasing B C / 2 ), whereas the second one is associated with constant B C / 2 (linear increase of B C / 1 ). The crossover region between these two regimes can be characterized by the parameter * being the solution of 1 = 2 , i.e., * ≈ ( 2 /π) ln(1/ ). The interesting observation implied from these plots is that when starts to decrease from its maximal value of 1, also B C / 1 starts to increase from its minimum value of 1. Moreover, when approaches its minimum value of 0, B C / 2 does not change and takes on its minimum value equal to 1. Thus, there is a connection between the C and anomaly of (Θ) dependence. The scaling of C found here [i.e., equation (6)] is in an agreement with these asymptotic expressions for C . However, the proportional factor for C in equation (6) could be different from that for * .

Specific heat
An anomalous dependence of (Θ) also implies an unusual dependence of the specific heat = − ( 2 / 2 ) as a function of reduced temperature Θ. The results for the FKM with small were presented in [9]. It was found (for = 0) that the function (Θ) changes its shape from the lambda-type for = 1, assuming two maxima in the intermediate range of 0.1 < < 0.3, and finally reaching only one maximum around Θ = 1/2 when < 0.01.
A similar evolution from the lambda-type behavior (for > C ) to a curve with one maximum around Θ = 1/2 for → 0 (and small ) is observed. Figure 6 presents some exemplary dependencies of specific heat (normalized by its maximum value max ) for two fixed values of and increasing values of . Similarly to the case of = 0, the maximum of occurs somewhat above Θ = 1/2, when there is a substantial decrease in the order parameter, although the maxima of and the derivative of (with respect to Θ) are sligthly displaced. With an increasing , it shifts towards higher Θ. The position of max calculated within the DMFT is located in slightly lower Θ than the one obtained within the HFA, which coincides with the presented results for (Θ) [cf., e.g., figures 1(a) and 6(a)]. The other drop in the specific heat always also occurs at the transition point Θ = 1, but it is relatively small (for small ) in comparison with max (due to very small values of (Θ) near C ) and it is not visible in figure 6. Note that in the nonordered phase, i.e., for Θ > 1, is also small ( / max is almost 0) due to small absolute temperatures just above C . In figure 6(b) the effects of numerical errors, which accumulate due to twofold derivation and very small changes in both free energy and , are visible (particularly for = 0), although the calculations were performed with a very high precision (in the figure, B-spline curves for the data points are plotted). This issue is absent for larger .

Conclusions and final remarks
In this paper we studied the anomalous behavior of the order parameter (Θ) versus the reduced temperature Θ in the EFKM for weak couplings and using the DMFT and the HFA approaches (in the large coordination number limit on the half-filled Bethe lattice). At this regime, both these methods are equivalent, therefore, the results are consistent with each other. But since it is easier to perform calculations using the HFA, most of the the data were obtained using this method. We focused mainly on finding the range of Coulomb interaction parameters in which there is an anomalous course of the order parameter as a function of temperature. For this purpose, we defined in formula (3) the parameter determining the degree of deviation of (Θ) from the standard mean-field = 1/2 Ising-like curve. This enabled us to determine the degree of this anomaly and to describe its evolution in a quantitative manner for any small and (their ratio is also arbitrary).