Formation of charge and spin ordering in strongly correlated electron systems

In this review we present results of our theoretical study of charge and spin ordering in strongly correlated electron systems obtained within various generalizations of the Falicov-Kimball model. The primary goal of this study was to identify crucial interactions that lead to the stabilization of various types of charge ordering in these systems, like the axial striped ordering, diagonal striped ordering, phase-separated ordering, phase-segregated ordering, etc. Among the major interactions that come into account, we have examined the effect of local Coulomb interaction between localized and itinerant electrons, long-range and correlated hopping of itinerant electrons, long-range Coulomb interaction between localized and itinerant electrons, local Coulomb interaction between itinerant electrons, local Coulomb interaction between localized electrons, spin-dependent interaction between localized and itinerant electrons, both for zero and nonzero temperatures, as well as for doped and undoped systems. Finally, the relevance of resultant solutions for a description of rare-earth and transition-metal compounds is discussed.


Introduction
The problem of inhomogeneous charge and magnetic ordering in strongly interacting electron systems is certainly one of the most intensively studied problems of the contemporary solid state physics. The reason is that the inhomogeneous charge ordering (e.g., the striped phases) were experimentally observed in many rare-earth and transition-metal compounds [1][2][3][4], like La 1.6 Nd 0.4 Sr x CuO 4 , YBa 2 Cu 3 O 6+x , Bi 2 Sr 2 Cu 2 O 8+x , La 1.5 Sr 0.5 NiO 4 , Na x CoO 2 , some of which exhibit a high temperature superconductivity. This phenomenon was most frequently studied in the literature within the Hubbard and t − J model [5][6][7][8][9][10][11][12]. Theoretical studies based on these models pointed to two possible mechanisms of forming the inhomogeneous charge ordering in these materials: (i) Kivelson, Emery, and co-workers [7] proposed that strongly correlated systems have a natural tendency toward phase separation and the inhomogeneous spatial charge ordering arises from a competition between this tendency to phase separate and the long-range Coulomb interaction which does not allow the electron density to stray too far from its average; and (ii) White and Scalapino [9][10][11][12] proposed that the stripe order arises from a competition between kinetic and exchange energies in a doped antiferromagnet which does not require long-range Coulomb forces to stabilize the stripes.
Moreover, shortly after the introduction of the spinless Falicov-Kimball model in 1986 by Kennedy and Lieb [13] and Brandt and Schmidt [14] it was found that this probably the simplest model of correlated electrons on a lattice is capable of describing various types of charge ordering including the periodic as well as phase-separated and segregated configurations [15][16][17][18][19][20][21][22][23][24][25][26][27][28]. This opened a new route of studying this important phenomenon. The model describes a system of itinerant and localized particles on the lattice, which interact only through the local Coulomb interaction. The Hamiltonian of the spinless Falicov-Kimball model can be written in the form where c + i (c i ) are the creation (anihilation) operators of the itinerant spinless particles at site i, t ij describes the hopping probability from site i to j and w i is the occupation number of the localized particles taking the value 1 or 0 according to whether the site i is occupied or unoccupied by the localized particle.
Accordingly, as the itinerant and localized particles are interpreted, we get different interpretations of the model. Historically, the first interpretation of the model: itinerant particles = electrons with spin up and localized particles = electrons with spin down, was already used in the original work by Hubbard [29] as the approximative solution of the original Hubbard model, in which one type of particles, e.g., electrons with spin down, are immobilized. Therefore, this approximation is sometimes referred to as the static Hubbard model. The second interpretation: itinerant particles = itinerant electrons and localized particles = immobile ions, comes from Kennedy and Lieb [13] and represents a very simple model of crystallization in solids. With respect to the fact that the Falicov-Kimball model was originally proposed to describe valence transitions and metal-insulator transitions in rare-earth compounds, one of the most frequent interpretations is to consider the conduction d electrons instead of itinerant particles and the valence f electrons instead of localized particles.
The greatest advantage of the spinless Falicov-Kimball model is its relative simplicity, which makes the model more accessible to analytical and numerical studies compared with e.g., the Hubbard or periodic Anderson model. Two diametrically different ways have been used to solve the spinless Falicov-Kimball model. The first direction represents analytical studies in the limit of infinite dimensions (D = ∞) and other analytical and numerical studies in reduced dimensions (D = 1 and D = 2).
Regarding the works devoted to the study of the Falicov-Kimball model in the infinite dimensional limit, it is worth noting that there are two good reasons to study the model in this limit. The first reason is that for D → ∞ the self-energy becomes local [30], which simplifies the analytical calculations to the extent that they permit to solve many problems exactly. The second reason for studying the model in the limit of infinite dimension is the fact that some physical quantities calculated for D = ∞ reproduce the three-dimensional results better than the one-dimensional solutions. A milestone in this direction is the work by Brandt and Mielsch [31], in which the exact solution of the spinless Falicov-Kimball model is presented for the symmetric case (the number of d electrons = the number of f electrons = L/2, where L is the number of lattice sites). The main result of this study was a determination of the critical transition temperature for a transition from the high temperature disordered phase to the low-temperature ordered (checkerboard) phase. In subsequent articles, Brandt et al. [31][32][33][34][35] built the foundations of physics of the Falicov-Kimball model in the limit of infinite dimensions, on which many other theoretical physicists built both the spinless and spin-one-half Falicov-Kimball model [36][37][38][39][40][41][42][43][44][45][46]. The results of these studies are summarized in an excellent review article by Freericks and Zlatic [47] devoted to exact solutions of the Falicov-Kimball model in the dynamic mean field theory.
Regarding the second direction, the analytical solutions of the Falicov-Kimball model in the limit of reduced dimensions, it should be noted that despite the huge efforts of theorists and relative simplicity of the model, so far only very few exact results for the ground state and thermodynamics of the spinless model have been obtained. In addition to the aforementioned evidence of the longrange arrangement at low temperatures and dimensions D 2, the following has been proven: (i) the absence of spontaneous hybridization at finite temperatures [48], (ii) the phase separation and periodic arrangement in the limit of strong Coulomb interactions and D = 1 [49], (iii) the phase separation in the two-dimensional model for selected values of electron concentrations and sufficiently large Coulomb interactions [21,27,28,50], (iv) the phase separation in the limit of strong Coulomb interactions for all dimensions [51,52].
Under these circumstances, numerical calculations seem to be very valuable. Although they do not always provide definitive answers to the questions searched, they are capable of pointing out the fundamental trends in the system, and thus help to understand the physics of the Falicov-Kimball model. Basically, there are two main directions of the numerical study of the Falicov-Kimball model. The first is an exact diagonalization of the model Hamiltonian on finite lattices over a complete set of accessible configurations of localized particles. As the number of configurations increases as 2 L , this method is severely limited by the size of clusters, which are eligible for the numerical study (L ∼ 36-40). In the one dimensional case such cluster sizes are sufficient to extrapolate the results obtained to the thermodynamic limit (L → ∞), yielding definite results on the ground state or thermodynamics of the model, at least for a certain area of model parameters (e.g., strong interactions). In the two-dimensional and three-dimensional cases, however, such cluster sizes are insufficient and, therefore, it is necessary to use approximate methods that somehow reduce the number of investigated configurations. This procedure was first used in the work by Freericks and Falicov [18] in the study of one-dimensional phase diagram of the spinless Falicov-Kimball model for selected concentrations of localized particles. Since the same formalism with minor modifications was used later by many other authors in the study of the ground states of the one-dimensional and two-dimensional models, we briefly summarize the main points of this algorithm. Consider some configuration w = {w 1 , w 2 , . . . , w L }, which corresponds to a specific distribution of particles localized on the lattice consisting of L sites, while the classical variable w i = 1, 0 indicates whether or not the site i is occupied by a localized particle. The energy of the system with N d itinerant electrons corresponding to the selected configuration is then given by where n(E, w) is the density of states corresponding to the configuration w and E F is the Fermi 42701-3 energy, which can be determined from the condition The density of states n(E, w) in the aforementioned expression is given by where G j (E) are the local Green's functions for which Freericks and Falicov [18] found the following expression where The problem is thus reduced to the calculation of infinite fractions ∆ ± j [53]. For the concentration of localized particles equal to n f = 0.5, Freericks and Falicov [18] investigated ten periodic configurations of ions (their list is given in table 1), which represent all different physical configurations Table 1. The complete list of all physically different periodic configurations of ions with the period p < 9 and n f = 1/2 [18]. k configuration period k configuration period 1 10 6 11101000 2 1100 7 11100100 3 111000 8 11011000 4 110100 9 11010100 5 11110000 10 11010010 of localized particles with the period less than 9 and a segregated phase (an incoherent mixture of completely empty and fully occupied lattice) with a density of states n seg (E) = (1 − n f )n empty (E) + n f n fully (E) .
Using the numerical analysis they obtained (for the case n f = 1/2) the coherent groundstate phase diagram (see figure 1) of the model that exhibits some general features: (a) The alternating phase w 1 = {1010 . . . 10} is the ground state at n d = 1/2 for all values of U as stated by previous investigations [13]. (b) The phase diagram tends to be simplified as the interaction strength increases indicating that many-body effects stabilize the system (this is a consequence of the segregation principle). (c) There is a trend for phases that disappear from the phase diagram as U increases to reappear as phase islands at even larger values of U (e.g., the configuration No. 3). (d) Phase islands of configurations not present at U = 0 may be formed at larger values of U (e.g., the configuration No. 8). (e) Some configurations are not the ground state for any value of U or n d (e.g., the configurations {11101000 . . .} and {11010100 . . .} do not appear).
A similar method was used later on by Gruber et al. [22] for the study of the one-dimensional phase diagram of the spinless Falicov-Kimball model for the so-called neutral case, where the concentration of itinerant electrons n d is equal to the concentration of localized electrons/ions n f . Unlike the previous work [18], Gruber et al. did not focus on studying the model only for selected values of n f , but they studied the comprehensive phase diagram of the model in the n f − U plane. Similarly to the previous case, the set of input configurations is not complete, since only periodic configurations with small periods and mixtures of these periodic configurations with an empty lattice have been taken into account 1 . The main result of their numerical and analytical studies is that the ground states of the spinless Falicov-Kimball model (n f 0.5) are either the most homogeneous distributions of localized particles (for U and n f large enough) or mixtures of periodic configurations and the empty lattice (U and n f small). Since the Fermi level for the mixtures of periodic configurations and the empty lattice lies in the conduction band, while for the most homogeneous configurations lie in the energy gap, the boundary between these two domains, is a boundary of metal-insulator transitions, which may be induced either by the Coulomb interaction or by changing the concentration of localized electrons.  [18]. The set of ground-state configurations consists of the segregated configuration and the periodic configurations with the smallest periods (see table 1).
For the two-dimensional case, this method was generalized by Watson and Lemanski [23] and later it was used by Lemanski, Freericks and Banach to study the ground-state phase diagram of the model in D = 2 [54]. In this case, the set of input configurations includes all periodic configurations with a unit cell having a smaller number of sites than the selected critical value N c and all possible mixtures thereof. The most interesting result of these studies was the observation of axial and diagonal striped phases of localized particles (see figure 2), suggesting that in the system of itinerant and localized particles, a sufficient mechanism leading to the formation of inhomogeneous charge arrangement is the local Coulomb interaction between these two electron subsystems. This is a much simpler mechanism of forming inhomogeneous charge stripes than the one considered earlier within the Hubbard model [5,6], respectively, within the t − J model [7][8][9][10][11][12]. Due to the incomplete input set of investigated configurations, an open question remains whether these results persist if the complete set of configurations will be considered, and especially, whether they persist for the more realistic three-dimensional case.  [54]. The typical ground-state configurations detected in the phase diagram. The large dots correspond to the occupied sites and the small dots correspond to the vacant sites.
One of the major shortcomings of the Falicov-Kimball model is that it does not include any spin interactions between electrons, and, therefore, it is not capable of describing the magnetic superstructures, which in many of rare-earth and transition-metal compounds coexist with charge ordering. This phenomenon was observed, for example, not only in nickelates [55][56][57], manganates [58], cobaltates [59,60], but also in materials exhibiting high-temperature superconductivity [1,61,62]. At present, there are still intensive discussions on possible mechanisms of the formation of inhomogeneous charge and spin ordering and its relation to the physical properties of the systems, e.g., high-temperature superconductivity. The easiest way of introducing spin interactions in the system of itinerant and localized electrons is to bind them by the Ising interaction. This idea was first used by Lemanski [63], who found that turning on the Ising interaction between itinerant and localized electrons leads to the stabilization of different types of charge and spin arrangement, including the axial and diagonal striped phases (see figure 3). Moreover, a number of simple rules of formation of various sorts of ground-state phases have been presented in reference [64]. Since these results were obtained on a restricted set of configurations, an open challenge for further theoretical studies was whether the nature of the ground state would remain unchanged after taking into account the complete set of configurations. For this reason we have decided to perform a systematic study (within an exact diagonalization and well-controlled approximate method described in the next section) of the ground state of the spinless as well as generalized spin-1/2 Falicov-Kimball model, in order to find the fundamental mechanisms of the formation of inhomogeneous charge and spin ordering in strongly correlated systems. Apart from the above mentioned local Coulomb interaction between d and f electrons, we have also investigated the effect of nonlocal Coulomb interaction [65,66], the correlated hopping [67][68][69][70][71], the lattice geometry [72,73], the dimension of the system [74,75], the anisotropic spin-dependent interaction [76][77][78] and the Hubbard interaction [76]. In what follows, we state a brief overview of the main results we have reached in our numerical studies.

Methods
To study the ground-state properties of the model Hamiltonians based on the spinless/ spin-1/2 Falicov-Kimball model we have used the method of exact diagonalization on finite clusters, where diagonalizations are performed over all possible distributions of localized particles, as well as the approximate method developed by us, in which the acceptance of configuration is realized using the principle of reducing the total energy of the system. Periodic boundary conditions are used in all the examined cases, since the fastest convergence of numerical results to the thermodynamic limit is observed for this type of boundary conditions.

Exact diagonalization technique
Although in the next major steps, the exact diagonalization method (EDM) will be illustrated for the spinless Falicov-Kimball model, the applicability of the method is much broader and with minor modifications it can be directly extended to the spin-1/2 Falicov-Kimball model, as well as the spin-1/2 Falicov-Kimball model with the Ising interaction between localized and itinerant electrons. The method is flexible with regard to the changes of the hopping matrix t ij and so it may be used, without any additional numerical complications, to study the effects of long-range and correlated hopping of electrons on the ground-state properties of the model.
Hereinafter we will use solely the interpretation of the spinless Falicov-Kimball model in which the itinerant particles are d electrons and the localized particles are f electrons from localized 4f or 5f states of rare-earth ions. Then, spinless Falicov-Kimball model can be written in the form where w i (w i = 1, 0) now describes the occupancy of the f orbital at lattice site i. It is important to note that for any distribution of f electrons w = {w 1 , w 2 , . . . , w L }, the Hamiltonian (8) is a single particle Hamiltonian in the representation of the second quantization where h ij (w) = t ij + U w i δ ij . Thus, the solution of the model (8) reduces to the problem of determining the spectra of matrix h(w) for different distributions of f electrons on the lattice of 42701-7 the size L. Since the problem is analytically solvable only for special types of configurations (e.g., the periodic configurations with the smallest periods), the only way to exactly solve this problem is to use the numerical diagonalization on finite clusters. Then, a fundamental task is to find a distribution of f electrons, for which the system has the lowest energy. The numerical algorithm for finding the configuration w 0 , which minimizes the energy of the system is as follows: (i) Having w, U and t ij fixed, we find all eigenvalues λ k of h(w). (ii) For a given N f = i w i we determine the ground-state energy E(w) = N d k=1 λ k of a particular f -electron configuration w by filling in the lowest N d one-electron levels. (iii) We find w 0 (examining all possible distributions of localized electrons), for which E(w, U ) has a minimum. Repeating this procedure for different values of model parameters one can immediately study the ground-state phase diagrams of the model.
Such exact calculations can be performed at present up to L ∼ 36-40 sites, which in some cases (the one-dimensional case and strong Coulomb interactions) is sufficient for an extrapolation of the results to the thermodynamic limit. In general, however, such cluster sizes are insufficient to obtain reliable conclusions about the behaviour of macroscopic systems in higher dimensions. Under these circumstances, the only way seems to be to use approximate methods. When selecting an appropriate approximate method one should have in mind the fact that charge and spin ordering as well as valence and metal-insulator transitions are very sensitive to the type of approximation used [79,80], and thus their description can only be successful within the approximations that introduce only small simplifications of the model system. Instead of searching for an appropriate method among the existing approximations, comparing them and excluding the least accurate candidates, we decided to develop a new numerical method, which would be sewn from the beginning on the Falicov-Kimball model, while retaining some degree of variability due to possible generalizations of this model.

Approximate method based on the reduction of the total energy
The natural starting point in building a new approximate method (AM) seemed to us to be the method of exact diagonalization. As stated above, within this method the single particle Hamiltonian h(w) is exactly diagonalized over all possible (2 L ) distributions of localized particles in order to find the only configuration w 0 , which minimizes the total energy of the system. This procedure is necessary, but not efficient. Much more efficient than passing through the complete set of configurations could be regulating the choice of configurations from the initial configuration w to the final configuration w 0 , under some criterion that would significantly reduce the number of configurations that should be examined. The most natural criterion seems to be a reduction of the total energy in a sequence of configurations from w to w 0 . This is the basic idea of our AM, the algorithm of which can be described as follows [74]: (i) Choose a trial configuration w = {w 1 , w 2 , . . . , w L }. (ii) Having w, U fixed, find all eigenvalues λ k of h(w). (iii) For a given N f = i w i determine the ground-state energy E(w) = N d k=1 λ k of a particular f -electron configuration w by filling in the lowest N d one-electron levels. (iv) Generate a new configuration w by moving a randomly chosen electron to a new position which is also chosen as random. (v) Calculate the ground-state energy E(w ). If E(w ) < E(w), the new configuration is accepted, otherwise w is rejected. Then, the steps (ii)-(v) are repeated until the convergence (for a given U ) is reached.
Of course, one can move instead of one electron (in step (iv)) simultaneously two or more electrons, thereby improving the convergence of method. Indeed, the tests that we have performed for a wide range of model parameters showed that the subsequent implementation of the method, in which 1 < p < p max electrons (p should be chosen at random) are moved to new positions, better overcomes the local minima of the ground-state energy. As usual, we have performed calculations with p max = N f . The main advantage of this implementation is that in any iteration step the system has a chance to lower its energy (even if it is in a local minimum), thereby the problem of local minima is strongly reduced (in principle, the method becomes exact if the number of iteration steps goes to infinity). On the other hand, a disadvantage of this selection is that the method converges slower than for p max = 2 and p max = 3. To speed up the convergence of the method (for p max = N f ) and still to hold its advantage we generate instead the random number p (in step (iv)) the pseudo-random number p that probability of choosing decreases (according to the power law) with increasing p. Such a modification considerably improves the convergence of the method.
Apart from the number of the moved electrons, the method was also tested on the optimum length of the iteration cycle M . It is obvious that if M → ∞, the method is exact. Unfortunately, with respect to the time factor, such a choice is not possible and we must consider only a finite number of iteration steps. The test process was realized on different finite clusters in the one-, two-and three-dimensional cases. We have observed that it is very convenient to divide the whole iteration process into several smaller independent cycles (from 5 to 10), among which the groundstate configuration with the lowest energy is selected. This also minimizes the problem of local minima. In table 2 we compare our numerical results on the one-dimensional cluster of L = 60 Table 2. The difference of the ground-state energies of the one-dimensional spinless Falicov-Kimball model calculated exactly and by our numerical method (∆ = |Eexact − Emet.|) for three different Coulomb interactions (U = 2, U = 4 and U = 8) for M = 200, 400 and 600 iterations. It should be noted that the exact results were obtained by the EDM on a set of the most homogeneous configurations, which are the ground states of the one-dimensional Falicov-Kimball model for all n f and U > 1.2. Therefore, the comparative tests were made for U = 2, U = 4 and U = 8. As seen in table 2 already relatively small number of iterations (M = 600) was sufficient to obtain the exact ground states for every N f on a sufficiently large cluster of L = 60. The same test was repeated on a twice larger lattice, where at 500 iterations, there were still observed small differences in the ground-state energy. These differences ranged in the order of about 10 −4 − 10 −9 , and their number decreased with the increasing number of iteration steps M . For M = 2000, we have found full consistency between our and exact results.

The effect of local Coulomb interaction
One-dimensional case We have started the study of ground-state properties of the spinless Falicov-Kimball model, in the limit U 1 [79]. This case was chosen for the reason that the analytical calculations showed [49] that the ground states of the model (at sufficiently large U ) can be only the most homogeneous distributions, when f electrons are as far apart as possible, taking into account the periodic boundary conditions. The knowledge of the ground states in the limit of strong interactions make possible a direct comparison between our results obtained on finite clusters with the results obtained in the thermodynamic limit (L → ∞) and at the same time it permits to specify more precisely the area of stability of these configurations. Numerical calculations were performed using the EDM, which permits to find the ground state of the model on the finite cluster of size L for any model parameters: n d , n f and U . Although n d and n f can be considered as independent parameters, we have bound them with the condition n f + n d = 1, since we are also interested in valence transitions, i.e., transitions induced by migration of valence f electrons to the conduction band. To minimize the finite-size effects, the model was studied on finite clusters from L = 4 to L = 24 sites. We have chosen the values of the Coulomb interaction from 1 to 10 with a unit step. The main result of our theoretical study was the finding that the most homogeneous configurations are the ground states of the spinless Falicov-Kimball model not only in the limit of strong Coulomb interactions, but for all the investigated values of U > 1, and for all f -electron concentrations. Since the Fermi level for the most homogeneous configurations always lies within the energy gap [22], all ground states for U > 1 are insulating, and thus the spinless Falicov-Kimball model is not capable of describing the metal-insulator transitions in this limit.
For this reason we have turned our attention to the case U < 1. Using the same method, we have performed a systematic study of the model for the selected value of Coulomb interaction (U = 0.6) and all even lattices from L = 16 to L = 48 lattice sites [80]. The main result of our study was the finding that for every L there exists a critical value of the f -electron occupation number N c below which the ground states are no longer the most homogeneous configurations, but the phaseseparated configurations that may be formally presented as an incoherent mixture of a configuration w and the empty lattice (w&w 0 ). A complete list of such configurations together with critical values N c is summarized in table 3. Thus, in accordance with the results by Gruber et al. [22] obtained on an incomplete set of configurations (the periodic configurations and the mixtures of periodic configurations with the empty lattice), we have found that the ground states of the spinless Falicov-Kimball model for U < 1 may be incoherent mixtures of the type w&w 0 , with a small difference, and namely, that the configuration w is not necessarily periodic.
To reveal the effect of Coulomb interaction U on the formation of charge ordering in these phase-separated configurations, we have performed an exhaustive numerical study of the model for a wide range of Coulomb interactions (U = 0.01, 0.02, . . . , 2) using our new AM that permits to treat much larger clusters. Thereby the finite size effects are considerably reduced. In figure 4 we present numerical results obtained for the ground-state phase diagram in the n f − U plane on finite clusters of L = 120 sites. We have found that the phase separation takes place for all Coulomb interaction U 1.2 for both small (n f < 0.25) and large (n f > 0.75) electron filling (the n f − U phase diagram is symmetric around n f = 0.5 line, and thus here we present only the results for n f 0.5) and that this domain, except a few isolated points/lines (denoted by * ) is continuous. This is an obvious difference between our phase diagram and one obtained by Gruber et al. [22] for periodic configurations and mixtures of these configurations with empty lattices, where large islands of the most homogeneous phase are observed in the phase-separated region. The second important difference is the existence of a narrow intermediate region (in our phase diagram) between the phase-separated configurations with regular (quasi-regular) distributions of f electrons (within w) and the most homogeneous domain (MHD). In this region (the medium gray area) the ground states are mixtures of an empty lattice and the aperiodic configurations w with two-molecule distributions in the middle of w and atomic distributions (a single occupied site) at the beginning and at the end of w (F 1 + F 2 ). This fact is clearly demonstrated in figure 5, where all ground-state configurations for U from 0.01 to 1.2 are displayed for two selected values of f -electron concentration n f = 1/20 and n f = 1/8. Analysing our numerical results we have found the following trends in the system: (i) In the weak coupling and low concentration limit there is an obvious tendency to form phasesegregated configurations (S) or mixtures of regularly (quasi-regularly) distributed n-molecules with the empty lattice (the phases F n with n = 3 and 4). (ii) With increasing U and n f , large n-molecules split into smaller ones, but their regular distribution persists.  From the viewpoint of rare-earth compounds where the role of itinerant particles is played by the d electrons and the role of immobile particles is played by the f electrons localized on the energy level E f (the term E f N f should be added to the model Hamiltonian (8)) it is interesting to transform the n f − U phase diagram into E f − U coordinates since the cuts of the E f − U phase diagram in the n f direction represent the valence transition at a given U . Since there is a  direct parametrization between E f and the external pressure p [81], the n f (E f ) behaviour is described (at least qualitatively) by the pressure induced valence changes in rare-earth compounds. The results of our numerical calculations for the E f − U phase diagram are summarized in figure 6. One can see that the phases with the largest area of stability in the E f − U phase diagram correspond to the periodic configurations with the smallest periods (p < 9) and the rational felectron concentrations. The number of phases with the relevant width is strongly reduced with increasing U , and thus only a few relevant phases (with p 5) form the basic structure of the phase diagram in the strong coupling limit. A detailed analysis of the model performed for U = 10 (L = 240 and L = 420) showed that some periodic phases with larger periods also persist in the strong-coupling regime, but their width is considerably smaller. A complete set of phases (with width w D > 10 −10 ) that have been determined numerically as the ground states of the model for U = 10 is shown in table 4. The phases with the smallest periods also persist in the weak coupling limit, but with decreasing U they are gradually suppressed by the periodic configuration with p 10 for E f < E c f (U ) and by phase-separated configurations for E f > E c f (U ). The corresponding picture of valence transition based on this phase diagram is displayed in figure 7. It is seen that the valence transitions have a staircase structure, where different phases correspond to different areas of stability. The largest stability regions correspond to the most homogeneous configurations with the smallest periods (n f =1/2, 1/3, 1/4, . . . ). These phases form a primary structure of the valence transition. Its characteristic feature is that it does not change with an increasing lattice size and therefore it can be used to represent the behaviour of macroscopic systems. The remaining phases form the secondary structure that unlike the previous one is very sensitive to the lattice size. This secondary structure is observed only for small values of the Coulomb interaction U and with an increasing U it rapidly disappears. Consequently, the valence transitions for intermediate values of Coulomb interactions consist of only a few valence steps, whose number is further reduced with an increasing U . For example, for 5 < U < 10 the valence transition consists of only four relevant transitions, and namely, from n f = 1 to n f = 2/3, from n f = 2/3 to n f = 1/2, from n f = 1/2 to n f = 1/3 and from n f = 1/3 to n f = 0 and for U > 10 there are even two relevant valence transitions: the first from n f = 1 to n f = 1/2 and the second from n f = 1/2 to n f = 0. Thus, we can conclude that the spinless Falicov-Kimball model is capable of describing two basic types of valence transitions, and namely, the transition from the integer-valence ground state into the inhomogeneous intermediate-valence state and transition from one inhomogeneous intermediate-valence state into another inhomogeneous intermediate-valence state. Moreover, our numerical results confirmed that the crucial role in the mechanism of valence transitions is played by the Coulomb interaction between itinerant and localized electrons.

Two-dimensional case
The combination of EDM and AM has been also used for the study of the charge ordering in a two-dimensional spinless Falicov-Kimball model [74]. In this case we were forced to limit ourselves only to the area of intermediate and strong Coulomb interactions (U > 1), whereas in the opposite limit U < 1 lattice effects were still strong, even on lattices with L = 400 sites.  Similarly to the one-dimensional case, these configurations are metallic, and thus the boundary between the phase-separated domain and the rest of the phase diagram is the boundary of metal-insulator transitions. For U 1, where the lattice effects are negligible, we have specified Figure 10. The ground-state configurations of the two-dimensional Falicov-Kimball model obtained for L = 16 × 16 and U = 2 (N f L/4).
the region of stability of this domain very precisely (figure 9). A surprising finding was that in the two-dimensional case the stability region of phase-separated domain shifts very significantly to higher values of U (U ∼ 3.3). This result is very important from the point of view of possible applications of the model for a description of metal-insulator transitions in rare-earth and transition-metal compounds. It is generally assumed that in these materials the values of the local Coulomb interaction U are much larger than the hopping integrals t ij , and, therefore, for a correct description of metal-insulator transitions in real materials one should consider the limit of U > t rather than U < t. Going with N f from zero to L/2 we Figure  We have observed the same picture for larger values of Coulomb interactions. The larger values of U only modify the stability regions of some phases, but no new configuration types appear. In particular, for U = 4, the region of phase segregation/separation is fully suppressed and the region of regular/quasi-regular distributions extends up to n f → 0.

Three-dimensional case
In principle, the same procedure as the one used in D = 1 and D = 2 can be also used in D = 3. Due to the numerical complexity of the problem in D = 3 we have performed numerical calculations only for selected values of the Coulomb interaction U representing the typical behaviour of the model at the weak (U = 1), intermediate (U = 2) and strong (U = 8) interactions [75]. In order to reveal the finite-size effects, numerical calculations were made on two different clusters of 4 × 4 × 4 and 6 × 6 × 6 sites. A direct comparison of numerical results obtained on 4 × 4 × 4 and 6 × 6 × 6 clusters showed that the ground-state configurations fall into several different categories whose stability regions are practically independent of L. For this reason we present here only the results obtained on L = 6 × 6 × 6 ( figure 12 and figure 13). The largest number of configuration types is observed in the weak-coupling limit. Going with N f from zero to half-filling (N f = L/2) we have observed the following configuration types for U = 1. At low f -electron concentrations, the ground-states are the phase-segregated configurations (f electrons clump together while the remaining part of lattice is free of f electrons). Typical examples of ground-state configurations from this region are depicted in figure 12 (left hand panel).
Above the region of phase segregation, we have observed the region of stripe formation (N f = 10, . . . , 20). In this region the f electrons form the one-dimensional charge lines (stripes) that can be perpendicular or parallel. This result shows that the crucial mechanism leading to the formation of stripes in strongly correlated systems should be a competition between the kinetic and short-range Coulomb interaction.
Going with N f to higher values, the stripes vanish and again appear at N f = 26, though in a fully different distribution. While at smaller values of N f the stripes have been distributed inhomogeneously (only over one half of the lattice), the stripes in the region N f = 26, . . . , 31 are distributed regularly.
Above this region a new type of configurations starts to develop. We call them diagonal charge planes with an incomplete chessboard structure, since the f electrons prefer to occupy the diagonal  planes with slope 1, and within these planes they form a chessboard structure. This region is relatively broad and extends up to N f ∼ 50. Then there follows the region in which the chessboard structure starts to develop. As illustrated in figure 13, the f electrons begin to preferably occupy the sites of sublattice A, leaving the sublattice B free of f electrons. Furthermore, the configurations that can be considered as mixtures of previous configuration types are also observed in this region. However, with increasing N f the configurations of chessboard type become dominant. Analysing these configurations we have found that the transition to the purely chessboard configuration is realized through several steps. The first step, the formation of the chessboard structure has been illustrated in figure 12. The second step is shown in figure 13. It is seen that the chessboard structure is fully developed in some regions (planes) that are separated by planes with an incomplete developed chessboard structure. Such type of distribution is replaced for larger values of N f by a new type of distributions (step three), where both regions with complete and incomplete chessboard structure have a three-dimensional character.

We have observed the same picture for intermediate values of Coulomb interactions (U = 2).
Larger values of U only slightly modify the stability regions of some phases, but no new configuration types appear. In particular, the domain of phase segregation, as well as the domain of stripe formation are reduced while the domain of diagonal planes with chessboard structure increases. This trend is also observed for larger values of U . In the strong coupling limit (U = 8) the phase segregated and striped phases are absent and the region of stability of the diagonal planes extends to relatively small values of N f ∼ 20. Below this value a homogeneous distribution of f electrons is observed.
For the same reasons as discussed above for the two-dimensional case we have investigated the stability of phase-separated (metallic) domain in the three-dimensional case and found (see figure 14) that the phase-separated region (and thus the metal-insulator transitions) extends up to U ∼ 5.5, which is a much larger value than in the two-dimensional case [82].

The effect of long-range hopping
Since the model including the electron hopping solely to the nearest neighbours may seem at first glance a very crude approximation, in order to have a more realistic description of electron processes in rare-earth compounds, we have generalized this model by taking into account the transitions to other neighbours [83,84]. Basically, there were two possible ways of performing such a generalization. The first way was to assign independent transition amplitudes for the first (t 1 ), second (t 2 ), third (t 3 ), . . . nearest neighbours, while the second way was to describe the electron hopping by a one-parametric formula with power decaying transition amplitudes t ij ∼ q |i−j| , where q 1. From the practical point of view, the second method is more suitable because it does not expand the model parameter space and has a clearer physical meaning, since the atomic wave functions have also the power law decay with an increasing distance. For this reason, for a description of electron hopping in the generalized model we have chosen the long-range hopping with power decreasing amplitudes. Explicit expressions of matrix elements t ij for the case of periodic boundary conditions in the one-dimensional case have the form: from which it follows that the case q → 0 approximates the nearest neighbours hopping (t ij = −t), while the q = 1 case corresponds to an unconstrained hopping, when the model is solvable exactly [85]. The results of our numerical simulations obtained within the AM are summarized in figure 15 in the form of n f − U diagrams. Figure 15 presents the one-dimensional ground-state

42701-20
states are configurations that can be considered as mixtures of the empty configuration and the alternating configuration w a (N f ) = {1010 . . . 10} with N f = 1, 2, . . . , L/2. The second domain located above U c and n f > 0.5 consists of two subdomains G 1 and G 2 , where the ground states are configurations of the following types: Analysing the energy gaps of ground-state configurations in all the above mentioned domains we have found (see figure 17) that only domains PSD 1 and PSD 2 are metallic, while all the remaining domains, including PSD 3 are insulating. Thus, in accordance with the q = 0 case, only the phase boundary between the small phase-separated domains PSD 1 and PSD 2 and the most homogeneous domain is a boundary of the metal-insulator transitions induced by Coulomb interaction (the felectron concentration).

The effect of nonlocal Coulomb interactions
One of the shortcomings of the basic variant of the spinless Falicov-Kimball model is that it neglects all nonlocal interactions between electrons, which immediately evokes the question of possible instability of numerical solutions discussed above with respect to the case when some of these nonlocal interactions are turned on. To answer this question we have examined the effects of two nonlocal interactions, and namely, the correlated hopping and the nearest-neighbour Coulomb interaction between d a f electrons

The effect of correlated hopping
Let us first describe the effect of the first term. This term is in the literature usually referred to as a term of correlated hopping, since it can be interpreted as a single-particle Hamiltonian describing the hopping of d electrons between the neighbouring d orbitals with an amplitude that explicitly depends on the occupancy of f orbitals. The selection of this term was motivated by earlier works [86][87][88], which showed its importance in describing the properties of strongly correlated systems, e.g., the superconducting state [89]. We have focused our attention on examining the effects of this term on charge ordering and valence and metal-insulator transitions [67]. A fundamental result of our numerical study is presented in figure 18, where we have displayed the phase diagram of the generalized model in the t − U plane for the half-filled band case n f = n d = 0.5. One can see that already very small values of the correlated hopping term lead to the instability of alternating phase w 1 = {1010 . . . 10} (which is for t = 0 the ground state of the model for all values of U > 0). For t < 0, this phase transforms to the alternating phase with a double period w 2 = {1100 . . . 1100} and for t > 0 on w 2 (U < 1), respectively, the segregated configuration w 3 = {11 . . . 100 . . . 0} (U > 1). Since the ground states corresponding to alternating phases w 1 and w 2 are insulating, while the ground state corresponding to w 3 is metallic, the phase boundary between w 1 and w 3 as well as between w 2 and w 3 is the boundary of metal-insulator transitions induced by the correlated hopping term. Similar instabilities were observed outside the symmetric case n f = n d = 0.5, which ultimately led to a completely different picture of valence and metal-insulator transitions for t = 0. Nonzero values of t reduce the total width of valence transitions (figure 18) as well as the width of stairs and above some critical value, the phase transition becomes continuous, initially only in certain areas (e.g., n f < 0.5 for t = 0.6), and finally, in the whole area (t 0.9). The parts of valence transition with staircase structure 42701-21 correspond to the most homogeneous (insulating) phases, while the continuous parts correspond to the segregated (metallic) phases. The border points between these two phases are, therefore, critical points of pressure (p ∼ E f ) induced metal-insulator transitions.
We have performed the same study in the two-dimensional case, where we have used our AM. We have found [68] that all characteristics discussed above, including the phase diagram for the half-filled band ( figure 19) as well as the picture of metal-insulator transitions remain unchanged Figure 19. Left: The t -U phase diagram of the two-dimensional Falicov-Kimball model with correlated hopping at half-filling (E f = 0, n f = n d = 0.5). The inset shows the t -U phase diagram for n f =1/4 and U > 2. Right: The ground-state configurations for n f =1/2 (w1, w2 and w3) and n f =1/4 (w4) [68].
in the two-dimensional case. A new and very interesting result is the observation of the charge striped ordering of the type w 2 in the half-filled band case. Since the mechanism of formation of inhomogeneous striped ordering in strongly correlated electron systems is a very intensively discussed topic in recent years, this result is also valuable from the cognitive perspective, because it opens up a new way to the study of this certainly interesting phenomenon. Figure 20. Representative ground-state configurations that form the basic structure of the phase diagrams in the n f −t plane: (a1) the chessboard phase, (a2, a3) the perturbed chessboard phases, (b1-b3) the diagonal stripes and perturbed diagonal stripes, (c1-c3) n-molecular phases, (d1-d7) the axial stripes and perturbed axial stripes, (e1, e2) the mixed phases and (f1, f2) the segregated and phase-separated configurations. The large dots correspond to occupied sites and the small dots correspond to vacant sites [70].
For this reason we have performed the exhaustive numerical studies of the model outside the half-filled band case. The primary goal of these studies was to identify all possible types of charge ordering induced by the term of correlated hopping. To fulfill this goal we have performed an exhaustive numerical study of the two-dimensional spinless Falicov-Kimball model for a wide range of Coulomb interactions U and correlated hopping t [70]. For each selected U and t the ground-state configurations for N f = 0, 1, . . . , L are calculated using our AM. To minimize the finite-size effects the same procedure is repeated on several different clusters. Of course, such a procedure demands in practice a considerable amount of CPU time, which imposes severe restrictions on the size of clusters that can be studied using this method (L = 4 × 4, 6 × 6, 8 × 8, 10 × 10, 12 × 12). Fortunately, we have found that the main features of the phase diagrams hold on all the examined lattices and thus can be used satisfactorily to represent the behaviour of macroscopic systems. In particular,

42701-23
we have found that for each L there is a finite number of basic types of ground-state configurations that form the basic structure of the phase diagram. This structure depends only very weakly on the size of clusters and covers practically the whole area of the phase diagram in the n f − t plane. Let us start a discussion of these phase diagrams (for U = 0.5, 1 and 2) with a description of configuration types that form their basic structure (see figure 20). (a) The chessboard phase (n f = 1/2). The f electrons occupy the A sublattice of the bipartite lattice and the B sublattice is empty (a1). The perturbed chessboard phases (n f close 1/2), denote the chessboard structure decorated by two-dimensional patterns of occupied or empty sites (e.g., a2, a3). (b) The diagonal stripes and perturbed diagonal stripes. The mentioned phases could be divided into three principal categories that are represented by examples b1, b2 and b3. (c) n-molecular phases (c1-c3). (d) The axial stripes and perturbed axial stripes. This group consists of several subgroups. In the first subgroup, ground states are configurations that can be constructed by spaced lines of occupations (vacancies) aligned with the lattice axes, into the chessboard structure (d1, d2). The second subgroup includes perpendicular axial stripes (d3, d4). The third subgroup is formed by simple axial and perturbed axial stripes (d5, d6). The last subgroup consists of n-molecular axial stripes (d7). diagrams of the Falicov-Kimball model with correlated hopping are presented for weak, intermediate and strong interactions. A direct comparison of these results reveals one general trend, and namely that the structure of phase diagrams is gradually simplified with increasing U and becomes very simple in the strong-coupling limit. In this case the basic structure of the phase diagram (in the n f −t plane) is formed by large segregated domains and several horizontal (band) domains cor-responding to the separated phases, n-molecular phases, diagonal and perturbed diagonal stripes, axial stripes (the first subgroup discussed above), mixed phases and finally perturbed chessboard configurations. Small deviations from the horizontal structure are observed for n-molecular phases for small (n f 0.2) and large (n f 0.7) f -electron concentrations. Comparing the phase diagram with conventional Falicov-Kimball model (t = 0) it is seen that the positive values of t do not essentially effect the ground states up to some critical value t c (n f ). However, at t c (n f ) the system exhibits a steep transition to the segregated phase that is stable for all n f and t > t c . A slightly different picture occurs for negative values of t . In this case the correlated hopping term induces new regions of axial stripes (the type d1, d2). The strong effect of negative t is apparent for felectron concentrations close to 1, where the diagonal configuration type (the type b) gradually disappears, while the segregated region is stabilized. As was discussed above, the phase diagrams become more complicated when the Coulomb interaction decreases. The simple band structure observed in the strong-coupling limit persists only for f -electron concentrations close to half-filling and it is suppressed gradually by axial-stripe configurations (the type d3-d7) for both positive and negative values of t . It should be noted that these axial stripe configurations have an arrangement principally different from the axial stripes observed in the conventional Falicov-Kimball model (d1, d2). The appearance of new types of axial stripes is one of the most interesting effects of correlated hopping on the ground-state properties of the two-dimensional Falicov-Kimball model. At the same time, this result positively answers the question whether the correlated hopping term can or cannot stabilize the stripe phases. Our results show that already relatively small values of t (positive as well as negative) stabilize this inhomogeneous charge ordering. Moreover, it was found (see figure 21) that the capability of correlated hopping to generate stripe ordering increases with a decreasing Coulomb interaction between localized and itinerant electrons. This opens up a new route towards the understanding of the nature of stripe formations in strongly correlated electron systems.

The effect of nearest-neighbour Coulomb interaction
Let us now discuss the effects of another nonlocal interaction term, i.e., the nearest-neighbour Coulomb interaction between the localized f and itinerant d electrons, that is of the same order as the term of correlated hopping. To study the effect of nonlocal Coulomb interaction U non on groundstate properties of the one-and two-dimensional Falicov-Kimball model, we have performed an exhaustive numerical study of the model for weak (U = 0.5), intermediate (U = 2) and strong (U = 8) on-site Coulomb interactions and for a wide range of nonlocal Coulomb interaction U non [65].
To determine the ground states of the model we used the EDM (up to L = 36 lattice sites) in combination with AM (up to L = 120 sites). We have started our study with the one dimensional case and U large (U = 8), which are relatively simple for a description.
Firstly, we have studied the ground-state phase diagram of the model in the n f − u non plane (u non = U non /U changes from 0 to 1 with step 0.001). To reveal the finite-size effects on the ground states of the model, we have performed numerical calculations for three different clusters of L = 12, 24 and 30 sites. The numerical results obtained for L = 30 are displayed in figure 22 3 . Our numerical results clearly demonstrate that already relatively small changes of u non can produce large changes in the ground-state f -electron distributions. Indeed, we have found that already values of U non around 60 times smaller than U (u c non = 0.015) are capable of fully destroying the most homogeneous distributions of f -electrons (that are the ground states of the conventional Falicov-Kimball model in the strong coupling limit for all n f ). It is interesting that only two configuration types are stabilized above the critical value of nonlocal interaction u c non . The first configuration type (w h ) is formed by the homogeneous distribution of [1100] and [10] clusters. A complete set of these configurations that are stable only in a very narrow region of u non (for 1/3 < n f < 2/3) are listed in figure 22 (the first panel below). The second type of configurations determined in the phase diagram of the Falicov-Kimball model with nonlocal Coulomb interaction in the strong coupling are the segregated configurations that are preferred as a ground state for all n f started at u non = 0.016.
In figure 22 we present the ground-state phase diagrams of the extended Falicov-Kimball model for intermediate (U = 2) and weak (U = 0.5) interactions. One can see that the main feature of the phase diagram found for U = 8, the nonlocal Coulomb interaction induced transition from the regular f -electron distributions to the phase-segregated distributions, also holds for smaller values of U , although the phase boundaries of different ground-state configuration types are now not so obvious due to the finite-size effects. A detailed analysis of the phase diagram for intermediate couplings showed that besides the most homogeneous configurations only three other configuration types enter the ground-state phase diagram, and namely: the segregated configurations (•), the weakly perturbed segregated configurations ( ) and w h (N f ) distributions (+). Contrary to U = 8, w h (N f ) configurations occur rarely at n f = 1/2, while weakly perturbed segregated configurations ( ) are observed on relatively large u non intervals. As shown in figure 22, in the weak interaction limit, the configurations w h (+) fully vanish and the set of ground-state configurations is much richer. Apart from the most homogeneous configurations (·) and the segregated configurations (•) we have determined a number of phase-separated configurations (denoted by ). A complete list of these configurations is given in figure 22.  We have performed the same calculations in two dimensions. To minimize the finite-size effects, the numerical calculations have been performed on three different clusters of 4 × 4, 6 × 6 and 8 × 8 sites. On the 4×4 cluster, the calculations have been performed by the EDM, and on larger clusters our AM was used.
Similarly to the one-dimension we have started our two-dimensional studies with U = 8. The ground-state phase diagram (calculated for 8 × 8) is shown in figure 23 (the first panel). Comparing this phase diagram with its one-dimensional counterpart one can find obvious similarities. In both cases the basic structure of the phase diagram only very weakly depends on L and consists of only three configuration types. Again one can see that the relatively small values of nonlocal interaction lead to changes of ground-state configurations from the regular distributions (·) to the segregated arrangements (•) straightforwardly, or through some n-molecular distributions (+) (usually diag-  figure 24 (the first row). While in the one-dimension the area of these configurations is stable only in isolated points of u non (u non = 0.015) and for 2/3 > n f > 1/3, in two dimensions these phases also persist for smaller n f and on the wider u non interval. On the other hand, it is interesting to note that the critical value of u c non = 0.018 above which all ground-state configurations are only the segregated phases is almost identical to 1D case (u c non = 0.016). Similarities between phase diagrams of 1D and 2D U=8, cases, can be also observed for intermediate and weak interactions. For both cases one can see a transition from ground states corresponding to u non =0 to phase-segregated distributions, due to the nonlocal Coulomb interaction. For U = 2, the obtained results show that the ground-state phase diagram consists of three different configuration types, denoted as regular distributions (·), phase-segregated distributions (•) and other types (+) including many n-molecular distributions (usually arranged to the "ladders" or to the blocks). Typical examples of these distributions (for L = 8 × 8) are displayed in figure 24 (the second row). Similarly, in the weak interaction limit the nonlocal Coulomb interaction u non prefers only a few types of ground-state configurations. We again observed regular distributions (·), phase segregated distributions (•) and some specific arrangements (+) discussed below. As was shown in figure 23 the critical value of u non , above which phase-segregated configurations are ground-states (for all n f ), shifts to higher values of u non . On the other hand, the critical value of u non , where the ground states of conventional Falicov-Kimball model are changed into the other ones, are observed already for u non ∼ 0.05 (for L = 8 × 8). Between these two boundaries one can find various f -electron distributions. In particular, there exist regular n-molecules, n-molecular "ladders", mixtures of chessboard structure and phase-segregated distributions, some periodic structures as well as the stripe formations (see figure 24, the third row). This clearly shows that nonlocal interaction can stabilize various types of inhomogeneous charge ordering in strongly correlated electron systems.

The effect of lattice geometry
There exists a large group of rare-earth and transition-metal compounds (e.g., GdI 2 , Na x CoO 2 , etc.) in which atoms instead of a square or cubic lattice decorate a triangular lattice (see figure 25) and they exhibit a number of anomalous physical characteristics [1,55,58,59,90,91]. From this point of view it is interesting to perform the same numerical study for a two-dimensional spinless Falicov-Kimball model on a triangular lattice. To reveal the effects of the lattice geometry on the ground-state properties of the Falicov-Kimball model we have started with the half-filled band case for which the nature of the ground state, its structural and energetic properties are quite understandable on the square lattice [73]. In this case the localized f electrons fill up one of two sublattices of the square lattice (the checkerboard structure) and the corresponding ground state is insulating for all U > 0. Thus, for finite interaction strength U there is no correlation-induced phase or metal-insulator transition at half-filling. The numerical calculations that we have performed in the half-filled band case on the triangular lattice revealed a completely different behaviour of the model for nonzero U in comparison with the square variant. Indeed, with increasing U we have observed a sequence of two correlation-  The stability regions of all the above described phases are displayed in figure 28. It is seen that the phase diagram of the triangular Falicov-Kimball model still keeps the band structure similarly to its square equivalent [70], with dominant area corresponding to the regular and quasi-regular

42701-30
Formation of charge and spin ordering in strongly correlated electron systems Having a complete set of ground-state configurations we have tried to construct the picture of valence transitions in different interaction limits. The resultant behaviours obtained for intermediate and strong values of Coulomb interactions are depicted in figure 29. Comparing these behaviours with their square counterparts [70] one can find significant differences. Indeed, while the valence transitions for square Falicov-Kimball model is symmetric with the largest step at n f = 1/2, the valence transitions on the triangular lattice are asymmetric and without any step at n f = 1/2. Instead of a significant step at n f = 1/2 there are now two large steps at n f = 1/3 and n f = 2/3. The origin of this different behaviour relates with the behaviour of the total energy E g (see figure 30), which instead of the total minimum at n f = 1/2 (the square Falicov-Kimball model) exhibits the total minimum at n f = 2/3 with the linear behaviour between n f = 1/3 and n f = 2/3. For this reason the valence transition for all finite clusters has a staircase structure, that follows the sequence n f = 1 → n f = 2/3 → n f = 1/3 to n f = 0 and it is practically independent of L for intermediate and strong interactions. For U → ∞, the total energy E g = 0 (for all f -electron concentrations), the steps at n f = 2/3 and n f = 1/3 vanish and the valence transition  is discontinuous from n f = 1 to n f = 0. Although in the strong coupling limit there exists a special type of phase separation, the dependence of the energy gap ∆ at the Fermi level on the f -electron concentration exhibits an insulating behaviour in the whole region. One can see (figure 31) that outside the interval n f ∈ (1/3, 2/3) the energy gap ∆ depends on the cluster size only very weakly and always has a finite value indicating an insulating behaviour. Inside the mentioned interval the energy gap slightly falls down and even depends on L, but the insulating behaviour is evident. The detailed analysis performed in the weak coupling limit at selected values of n f (n f = 2/3 and 8/9) showed (see figure 32), that there exists a critical value of U (U c ∼ 1 for n f = 2/3 and U c ∼ 3 for n f = 8/9) below which the Falicov-Kimball model on the triangular lattice also exhibits a metallic behaviour. From this point of view, the principal question occurring in the literature, namely, whether the systems with triangular lattice are capable of exhibiting the Mott-Hubbard transition, has been answered positively.

Spin-1/2 Falicov-Kimball model without the Ising interaction
The transition from spinless to spin version of the Falicov-Kimball model is formally trivial, and it is sufficient to add the spin variable to the creation and annihilation operators of itinerant and localized electrons: It should be noted, however, that from the physical point of view the spin Falicov-Kimball model describes a fully different physical reality. While in the spinless model, all states with double occupancy are projected out (the Coulomb interaction between d electrons with opposite spins U dd as well as between f electrons with opposite spins U f f are infinitely large) in the spin model, such states are permitted (U dd = 0, U f f = 0). The total omission of Coulomb interactions between the itinerant and localized electrons with opposite spins is, however, a too crude simplification of physical reality in real systems, and therefore as the first step of our study we have generalized the model Hamiltonian (15) by the term: describing the Coulomb repulsion of two f electrons with oppositely oriented spins localized at the same position, and studied its effect on the valence and metal-insulator transitions [92]. Since (16) does not violate the commutativity of f + iσ f iσ with the total Hamiltonian of the system, one can again replace f + iσ f iσ by the classical variable w iσ = 0, 1 and use for the study of the generalized spin-1/2 Falicov-Kimball model the same procedures and methods as for the Furthermore, we have found that the transition from w p to w s is realized through the following steps The last observation is very important for the extrapolation of small-cluster exact-diagonalization calculations since it allows us to avoid technical difficulties associated with a large number of configurations and consequently to study much larger systems. Figure 33 presents numerical results for critical interaction strengths U f f1 and U f f2 as functions of the f -electron occupation number n f obtained for U = 3 and L = 500. It is seen that there is a relative large region of U f f values where the configurations with a nonzero number of local f pairs are the ground states. The fact that the f electrons form the local f pairs, in spite of a relatively large repulsive interaction U f f , indicates that there is an attractive interaction that is capable of overcoming this direct repulsion. One of the most important results for the spin-1/2 Falicov-Kimball model is that the interaction of the localized f electrons with the itinerant d-band electrons leads to an effective on-site attraction between the localized f electrons. It is interesting to study whether this feature changes the picture of valence and metal-insulator transitions found for the spinless version of the model. The numerical only the most homogeneous configurations have been considered, and the MF curve represents the mean-field result obtained using the exact density of states [92].
results for E f dependence of n f (calculated for configurations of type (17) or (18)

273) there are discontinuous transitions from an integer-valence state n f = 1 into an inhomogeneous intermediate-valence state
f f the transitions are discontinuous from n f = 1 to n f = 0. They take place at E c = −U c f f independently of U f f . In the weak coupling limit, we have found strong finite-size effects and, therefore, we have turned our attention to the case U f f = ∞ that permits to reduce the total number of the investigated f -electron configurations from 4 L to 2 L , and thus to greatly increase the size of the clusters studied, which is very important for the correct analysis of structural properties of the model [93].

42701-34
Formation of charge and spin ordering in strongly correlated electron systems To reveal the basic structure of the phase diagram in the n f -U plane (E f = 0) we have performed an exhaustive study of the model on finite (even) clusters up to 36 sites. For fixed L, the numerical calculations have been performed along the lines discussed above with a step ∆N f = 2 and ∆U = 0.05. The results of numerical computations are summarized in figure 35. These results show that the phase diagram of the spin-1/2 Falicov-Kimball model consists of three main domains: the most homogeneous domain MHD (in figure 35 denoted as ·) and two phase separation domains PSD 1 (denoted as •, × and +) and PSD 2 (denoted as ). In the MHD the ground states are configurations in which the atomic or n-molecule clusters of f electrons are distributed in such a manner that the distances between two consecutive clusters are either d or d + 2. Furthermore, distribution of the distances of d and d + 2 has to be the most homogeneous. Two basic types of the ground state configurations that fill up practically the whole MHD are displayed in table 5. In the PSD 1 the ground states are configurations in which all f electrons are distributed only in one part of the lattice (w) while another (connecting) part of the lattice (w e ) is free of f electrons (phase separation). In accordance with Gruber et al. [22] we refer to such configurations as mixtures and denote them by w&w e . We have found three basic types of configurations w which form these mixtures: (i) aperiodic atomic configurations w a = {10 k1 10 k2 . . .  Table 5. Two basic types of the most homogeneous configurations that fill up practically the whole MHD for L = 24 [93].
100001000000100001000000 110000000000110000000000  6  100100001001000010010000 110000001100000011000000  8  100100100100100100100100 110000110000110000110000  10  110010011001001001100100 110011000011001100110000  12  110011001100110011001100 110011001100110011001100  14  111001110011100111001100  16  111100111100111100111100  18  111111001111110011111100  20 111111111100111111111100 Falicov-Kimball model [22] shows that this region roughly corresponds to a region  The second step in our numerical studies has been the extrapolation of small-cluster exact diagonalization results on large lattices. In figure 36 we present the ground state phase diagram of the spin-1/2 Falicov-Kimball model obtained for L = 100 on the extrapolated set of configurations that includes practically all possible types of the ground-state configurations found on finite lattices up to 36 sites. In particular, we have considered (i) the most homogeneous configurations w h of the type a and b (see table 5), (ii) all mixtures w a &w e with k i smaller than 6, (iii) all mixtures w b &w e with n i and k i smaller than 6, (iv) all mixtures w d &w f with periods smaller than 12, and (v) all segregated configurations. One can see that all fundamental features of the phase diagram found on small lattices hold on much larger lattices too. Of course, the phase boundaries of different regions corresponding to mixtures w a &w e (a), w b &w e (b), w c &w e (c) and w d &w f (d) are now more obvious. Since the mixtures w a &w e , w b &w e , w c &w e and w d &w f are metallic [22] one can expect that the phase boundary between the MHD and PSD is also the boundary of the correlation induced metal-insulator transition. To confirm this conjecture it is necessary to show that the most homogeneous configurations from the MHD are insulating, or, in other words, that there is a finite energy gap ∆ at the Fermi energy in the spectra of these configurations 4 . The numerical results for the most-homogeneous configurations of the type a and b that fill up practically the whole MHD are displayed in figure 37. The insulating character of these configurations is clearly demonstrated by the finite ∆ that exists for all nonzero f -electron concentrations n f and Coulomb interactions U . Thus we can conclude that the spin-1/2 Falicov-Kimball model undergoes (on the phase boundary between the MHD and PSD) the correlation induced metal-insulator transition that is accompanied by a discontinuous change of the energy gap ∆. This result is similar to what was found for the spinless Falicov-Kimball model by numerical [22] and analytical calculations in the strong [49] and weak-coupling limit [42].
Of course, this fact has to lead to a different picture of valence and metal-insulator transitions (figure 37) induced by pressure (increasing E f ). We have found that: (i) In the strong-coupling limit (U > 4), the model exhibits a pressure induced discontinuous insulator-metal transition from an integer-valence state (n f = 1) into another integer-valence state (n f = 0).

Spin-1/2 Falicov-Kimball model with the Ising interaction
Despite the unquestionable success of the spinless Falicov-Kimball model in describing the charge ordering in strongly correlated systems, this relatively simple model was not capable of accounting for all aspects of real experiments, and namely, that the charge superstructure is often accompanied by a magnetic superstructure. In order to more realistically describe electronic and spin processes in real materials, the original spin-1/2 Falicov-Kimball model has been extended by the spin-dependent interaction (of the Ising type) between localized and itinerant electrons [63,64,76] (19) and, in addition, the local Coulomb interaction between the f electrons in the limit of U f f → ∞ has been included. Thus, from the major interaction terms that come into account for the interacting d and f -electron subsystems, only the Hubbard type interaction between the spin-up and spin-down d electrons has been omitted in the Hamiltonian (19). In the work [63] Lemanski presented a simple justification for the omission of this term based on an intuitive argument: the longer time electrons occupy the same site, the more important becomes the interaction between them. According to this rule the interaction between the itinerant d electrons (U dd ) is smaller than the interaction between the localized f electrons (U f f ) as well as smaller than the spin-independent interaction between the localized and itinerant electrons U .
We have tried to verify the validity of these intuitive statements numerically and we have found that their validity is restricted only to the region of intermediate and strong Coulomb interactions U between the localized and itinerant electrons. Our numerical proof is based on the exact diagonalization of the model Hamiltonian (19) extended by H dd on small finite clusters. Such a study is numerically very exhaustive, since it is necessary for each configuration of localized electrons/spins to perform the Lanczos diagonalization [94] of the Hamiltonian over the full Hilbert space of the model and this process should be repeated L N f times. Of course, such a procedure demands in practice a considerable amount of CPU time, which imposes severe restrictions on the size of clusters that can be studied within the exact diagonalization calculations. For this reason we were able to exactly investigate only the clusters up to L = 12. Fortunately, it was found that in some parameter regimes the ground-state characteristics of the model are practically independent of L and thus already such small clusters can be satisfactorily used to represent the behaviour of macroscopic systems. In particular, we have studied the stability of the ground-state configuration

42701-38
Formation of charge and spin ordering in strongly correlated electron systems have been observed. In these cases the typical values of U c dd are of the order of 0.5 but for some N f even much smaller values were found. Thus, we can conclude that the Hubbard type interaction between the spin-up and spin-down d electrons can be neglected in the strong interaction limit between the localized and itinerant electrons (U 4). We have obtained the same conclusion for the spin-dependent interaction J > 0 [78]. For this reason, all subsequent calculations on the spin-1/2 Falicov-Kimball model with spin-dependent Coulomb interaction J between the f and d electrons have been done at U = 4.
We have started the study of the effect of anisotropic, spin-dependent interaction between localized and itinerant electrons on the ground-state properties of the model in the simplest case In general, the spin-dependent interaction J stabilizes the FP and PP phases, while the NP phase is gradually suppressed with increasing J. In figure 39 we present a complete set of groundstate configurations from the NP region. Among them one can find different types of periodic and nonperiodic configurations, but the most interesting examples represent configurations formed by antiparallel ferromagnetic domains, that convincingly illustrate the cooperative effects of spindependent interaction J between the localized and itinerant electrons.
As the next step, we have performed numerical studies of the model for the case N f = L. From the numerical point of view this case is considerably exact, since now we have to minimize the ground-state energy not only over all different spin configurations but also over all different f -electron distributions. This takes a considerable amount of CPU time and for this reason we were able to investigate exactly only the clusters up to L = 24. We have found that the spectrum of magnetic solutions that yields the model for the NP and PP phases is very rich. Indeed, for L = 24 and J = 0.5 we have found 140 different NP phases and 20 different PP phases that enter the N f − N d phase diagram. The main configuration types, with the largest stability regions are presented in figure 40 in the form of a skeleton phase diagram. From the NP phases, the largest stability region (denoted by I) corresponds to configurations of the type ↑ n ↓ n 0 L−2n . The second largest region (denoted by II) corresponds to NP configurations of the type [↑ 0 n ↓ 0 n ] k 0 L−2k(n+1) . Typical examples of the NP ground states from the central region of the phase diagram represent periodic configurations of the type ↑ n 0 m ↓ n 0 m (below the main diagonal) and configurations of the type In order to minimize the finite-size effects we have performed the same study using our AM on  pearing types of ground states in the N d − J phase diagram) are displayed in figure 42. Again one can see that the spectrum of magnetic solutions that yields the Falicov-Kimball model extended by spin-dependent interaction is very rich. In addition to the FP phase (that the stability region shifts to higher d-electron concentrations when J increases) there are various types of NP and PP structures like the antiparallel ferromagnetic domains (2-3), the axial magnetic stripes (4-7), the diagonal magnetic stripes (8)(9)(10)(11) and the perturbed diagonal magnetic stripes (12)(13)(14)(15). This again demonstrates strong effects of the spin-dependent interaction on the formation of magnetic superstructures in the extended Falicov-Kimball model and its importance for a correct description of correlated electron systems. Finally, we have performed the same numerical study for the case N f < L.
The typical examples of the resulting charge and spin ordering are displayed in figure 43. Among them one can find various types of phase-segregated (1), phase-separated (16) and n-molecular (3) configurations with FP, PP and NP ground states as well as various types of axial (9) and diagonal (5) magnetic/charge stripes. In general, we have observed that the system shows tendency towards phase segregation for small and large d-electron concentrations, while near the n d = 1 point the system prefers to form various types of axial and diagonal stripes.
For the most physically interesting cases N f + N d = L and N f + N d = 2L we have also performed a detailed analysis of charge and magnetic ordering. We have started our study with the case N f + N d = 2L, which is slightly simpler for a description. As shown in figure 44 the ground states for N f + N d = 2L are antiferromagnetic (AF) with alternating pattern, where the electrons (for n f < 1/2) or holes (for n f > 1/2) form axial distributions. Our calculations showed that these inhomogeneous stripe distributions are stable for large Coulomb interactions, while decreasing U  The second difference is that although with increasing f -electron concentration the ground states are the AF, these AF arrangements are formed by F ordered clusters (domains). In addition, for n f = 1/4 and n f = 1/3 (L 144) a new type of stripes (known as the ladders) occurs. And finally, a detailed analyses showed that there exists a critical f -electron concentration n c f ∼ 1/4 below which the ground-states are phase separated.

Ground-state properties of Na x CoO 2
The Na x CoO 2 system, which forms the basis for a quasi-two-dimensional transition metal oxide superconductor (T c = 4.5 K) when hydrated [95] shows a wide variety of unexplained behaviours in the accessible range 0 < x < 1. It consists of alternate stacks of electronically active triangular CoO 2 layers (with edge sharing CoO 6 octahedra) separated by Na layers that act not only as spacers, leading to electronic two-dimensionality ( figure 47), but also as charge reservoirs [96][97][98]. With increasing Na content x, the phase diagram nonhydrated Na x CoO 2 system exhibits a succession of ground states [98] changing from a paramagnetic metal at x = 0.3 through a chargeordered insulator at x = 0.5, a "Curie-Weiss metal" for x = 0.75, and finally a magnetically ordered state for x > 0. 75. At x = 0.5, detailed electron and neutron diffraction measurements [99,100] revealed an ordering of the Na ions along one crystallographic direction which decorates the chains of Co ions with different amounts of charge [100]. It is suggested that there are strong correlations between the Na ions and electrons (holes) from the CoO 2 layers, even though they occupy separate layers, and that a charge ordering in the CoO 2 layers is induced by sodium ion ordering. Unfortunately, numerical calculations performed at the Density Function Theory (DFT) and DFT+U level by Li et al. [101] did not confirm this conjecture. In the spirit of these results we have tried to describe the conducting properties and a formation of charge ordering in the CoO 2 layers within a relatively simple model that takes into account both the interplane interactions between the ordered Na ions and electrons (holes) form the CoO 2 planes as well as the intraplane interaction between electrons within CoO 2 layers [72]. The system is modelled by the spinless Falicov-Kimball Hamiltonian that automatically projects out the states with two holes on the 3d (d z 2 ) orbitals. Thus, our starting Hamiltonian can be written as a sum of three terms: The first term of (21) describes hopping of 3d-electrons on the triangle network of Co 4+ ions. These intersite hopping transitions are described by the matrix elements t ij which are −t if i and j are the nearest neighbours and zero otherwise. The second term describes the interplane interaction between the Na ions and the d electrons from the CoO 2 layers, where ω i represents the static potential of the first, second and third nearest Na ions to site i. The strengths of the Coulomb interactions corresponding to the first (U 1 ), second (U 2 ) and third (U 3 ) nearest neighbours are considered as the model parameters. The third term represents the nearest-neighbour Coulomb repulsion between two d electrons on the Co triangle network.
We have started our numerical study with the case x = 0.5 and V = 0. As mentioned above, at x = 0.5, the ordered Na ions form one-dimensional zig-zag chains and so, one can easily reconstruct the static potential ω i for any selected values of U 1 , U 2 and U 3 . Having ω i , the Hamiltonian (21) can   U 1 the charge gap at the Fermi level opens, and the system undergoes the metal-insulator transition at U 1 ∼ 5.8. This confirms a supposition that the interplane interactions between the ordered Na ions and the d electrons from Co layers play a crucial role in the stabilization of the insulating state in the Na 0.5 CoO 2 material.
From this point of view it is interesting to ask if these interactions could also stabilize the charge ordering within the Co layers. To answer this question we have calculated the on-site occupation where |ψ G is the ground state of Hamiltonian (21) obtained directly from the exact numerical diagonalization solution. The typical results for n i within the insulating phase are shown in figure 49 for a finite cluster of L = 8 × 8 sites.
Our results clearly demonstrate that an inhomogeneous charge ordering, of the stripe-like form, develops within the Co layers. For the sites below (above) the Na ions the occupation number n i is equal to 0.93 (the stripes a), while for the remaining sites n i = 0.07 (the stripes b). Thus, we arrive at the conclusion that interplane interactions between the ordered Na ions and the d electrons from the CoO 2 layers not only stabilize the insulating ground state in Na 0.5 CoO 2 , but they are also capable of describing the formation of inhomogeneous charge ordering within the CoO 2 layers.
Since the Coulomb interactions between d electrons inside CoO 2 planes are not negligible in comparison with interplane interactions, it was necessary to examine what happens if these interactions are switched on (V = 0). We have solved this task in a manner analogous to the case V = 0, i.e., by a direct calculation of the d-electron density of states and the d-electron on-site occupation n i , with only one technical complication, namely, that the diagonalization of the full Hamiltonian is performed by the Lanczos method, because for V = 0, the Hamiltonian (21) is no longer a single particle. The results of our numerical calculations for n i obtained for four different values of the intraplane interaction V are summarized in table 6. One can see that for nonzero interplane interactions, the nearest-neighbour intraplane interaction does not destroy the charge ordering induced by Na ions, but on the contrary, increasing V further stabilizes the inhomogeneous charge ordering. This effect is very strong and for sufficiently large V , practically all d electrons are localized on sites below (above) Na ions (the stripes a). Taking into account the fact that n i = 0 (in our notation) corresponds to Co 4+ and n i = 1 to Co 3+ , these results are consistent with the model proposed recently by Choy [102] for a description of the sign change of the Hall-coefficient (as a function of temperature), in which the rows of Co 4+ ions alternate with rows of Co 3+ ions. This indicates that our model, in spite of its relative simplicity, may still contain the relevant physics of cobaltate systems.
We have used the same procedure for the case of x = 0.5. Our main motivation in this case was to answer the question whether our simple model is capable of describing correctly the conducting properties of Na x CoO 2 for x < 0.5 as well as for x > 0.5. In figure 50 we present numerical results for d-electron density of states for three selected experimentally confirmed Na structures with x = 0.3, x = 0.71 a x = 0.75, representing the typical examples from the region of x < 0.5 and x > 0.5 (the experimentally proposed Na ion distributions are shown in figure 50, together with our results for n i ).
It is seen that in all three investigated cases the Fermi level lies inside the conduction band, and thus all three cases correspond to a metallic state which is perfectly consistent with experimental measurements [98].

Magnetization processes in rare-earth tetraborides
A spin system is frustrated when all local interactions between spin pairs cannot be satisfied at the same time. Frustration can arise from the competing interactions or/and from a particular geometry of the lattice, as seen in the triangular lattice. The Shastry-Sutherland lattice was considered more than 20 years ago by Shastry and Sutherland [103] as an interesting example of a frustrated quantum spin system with an exact ground state. It can be described as a square lattice with antiferromagnetic couplings J between nearest neighbours and additional antiferromagnetic couplings J between next-nearest neighbours in every second square (see figure 51). This lattice attracted much attention after its experimental realization in the SrCu 2 (BO 3 ) 2 compound [104]. saturated magnetization m s ) in this material [104] stimulated further theoretical and experimental studies of the Shastry-Sutherland lattice [105,106].
Similar phenomena of magnetization plateaus is also observed in the rare-earth tetraborid TmB 4 [107]. Since fully polarized state can be reached for experimentally accessible magnetic fields, this compound permits exploration of its complete magnetization process. It was found that the magnetization diagram of TmB 4 consists of magnetization plateaus located at small fractional values of m/m s =1/9, 1/8, 1/7 of the saturated magnetization, followed by the major magnetization plateau located at m/m s = 1/2. Due to strong crystal field effects, the effective spin model for TmB 4 has been suggested to be described by the spin-1/2 Shastry-Sutherland model under strong Ising (or easy-axis) anisotropy in a magnetic field h [107] where S z i = ±1/2 denotes the z-component of a spin-1/2 degree of freedom on site i of a square lattice and J, J are the antiferromagnetic exchange couplings between all nearest neighbour bonds (J) and next-nearest neighbour bonds in every second square (J ), as indicated in figure 51 (left hand panel).
Numerical simulations obtained within the Monte-Carlo [108] and tensor renormalization group methods [109] on large systems showed that the Ising model on the Shastry-Sutherland lattice exhibits, in the presence of the magnetic field, a magnetization plateau only at 1/3 of the saturated magnetization. The existence of the magnetization plateau at only 1/3 of the saturated magnetization and its absence at 1/2 indicates that it is necessary to go beyond the classical Ising limit to reach the correct description of the magnetization process in TmB 4 and other rare-earth tetraborides. The first attempt of this kind has been made by Meng and Wessel [108] who studied the spin-1/2 easy-axis Heisenberg model on the Shastry-Sutherland lattice with ferromagnetic transverse spin exchange using quantum Monte-Carlo and degenerate perturbation theory. Besides the magnetization plateau at 1/3 of the saturated magnetization they found a further plateau at 1/2, which persists only in the quantum regime.
Quite recently we have proposed an alternative model [113] of stabilizing the magnetization plateaus in the rare-earth tetraborides based on the fact that these materials, in contrast to SrCu 2 (BO 3 ) 2 , are metallic. Thus, for a correct description of ground-state properties of rare-earth tetraborides one should take into account both spin and electron subsystems as well as the coupling between them. Supposing that electron and spin subsystems interact only via the spin dependent

42701-49
Ising interaction J z , the Hamiltonian of the system can be written as The model described by (23) is a straightforward extension of the spin-1/2 Falicov-Kimball model with anisotropic spin-dependent interaction discussed in detail in [63]. The only differences are that we consider here a direct spin interaction (of the Ising type) between the localized spins and that the underlying lattice is of the Shastry-Sutherland type.
To examine the magnetization curve corresponding to the model Hamiltonian we have used our AM. Using this method we have performed exhaustive numerical studies of the model (23) for a wide range of model parameters h, J z , t (the hopping integral between the nearest-neighbours), t (the hopping integral between the next-nearest neighbours) and J/J = 1 selected based on the experimental measurements [107]. To exclude the finite-size effects, the numerical calculations have been performed for several different Shastry-Sutherland clusters consisting of L = 8 × 8, 10 × 10 and 12 × 12 sites. The most important result obtained from these calculations is that the switching on of J z and t leads to a stabilization of the uniform ground-state spin arrangement consisting of parallel antiferromagnetic bands separated by ferromagnetic stripes. A complete list of the ground-state spin arrangements (for 0 < m sp /m sp s < 1) that are stable on finite intervals of magnetic field values is depicted in figure 52. The second important observation is that the width w of the antiferromagnetic bands cannot be arbitrary, but fulfills severe restrictions. Indeed, we have found that with the exception of the case m sp /m sp s = 1/2, in all the remaining cases the permitted width of the antiferromagnetic band is only w or w + 2, where w is an even number. This fact is very important from the numerical point of view, since it allows us to perform the numerical calculations on much larger clusters with the extrapolated set of configurations of the above described type. The resulting magnetization curves obtained on the extrapolated set of ground-state spin configurations consisting of parallel antiferromagnetic bands of width w (w and w + 2) separated by ferromagnetic stripes are shown in figure 53 for the selected values of model parameters that represent the typical behaviour of the model. One can see that the switching on of the spin-dependent interaction J z and taking into account the electron hopping t on the Thus, our numerical results show that besides the pure spin mechanism (e.g., the easy-axis Heisenberg model on the Shastry-Sutherland lattice [108]) of stabilization the magnetization plateaus in rare-earth tetraborides, there also exists an alternative mechanism based on the coexistence of electron and spin subsystems that are present in these materials. From this point of view it is interesting to compare in detail the ground states obtained within these two different approaches. For m sp /m sp s = 1/3 our results are identical to the ones obtained within the Ising [107,109] as well as easy-axis Heisenberg [108,114] model on the Shastry-Sutherland lattice. The accordance between our and the easy-axis Heisenberg solution [108] is surprisingly found for m sp /m sp s = 1/2 as well. In this case both approaches predict the sequence of parallel antiferromagnetic and ferromagnetic stripes. For m sp /m sp s = 1/5 our results postulate a new type of spin ordering. Finally, it should be noted that more exhaustive studies of the model performed on much larger lattices have revealed the existence of magnetization plateaus at m sp /m sp s = 1/7, 1/9 and 1/11 (in accordance with experimental measurements in TmB 4 ) but the stability regions of these phases are much narrower in comparison with the ones of 1/2, 1/3 and 1/5 plateau phases.

Doping-induced valence changes in rare-earth compounds
The substitution of one type of atoms/ions by another type changes electronic relations in a given material which permits to directly study manifestations of electron correlations on physical properties of the system. From the perspective of the Falicov-Kimball model, which was originally introduced to describe valence and metal-insulator transitions in rare-earth compounds, it is therefore particularly interesting to examine what predictions the model is capable to provide for the case when we replace the rare-earth ions by other ions. Since previous theoretical works have shown [80,115] that the Falicov-Kimball model can provide good qualitative predictions for the transport, magnetic and spectroscopic properties of some mixed valence compounds (e.g., SmB 6 , SmS, etc.), we have focused on studying the effects of doping just for this type of materials. For SmB 6 it has been found, for example, that the substitution of Sm by nonmagnetic divalent ions (e.g., Sr 2+ , Yb 2+ ) increases the average samarium valence (the average occupancy of f orbitals decreases), while the substitution of Sm by nonmagnetic trivalent ions (e.g., Y 3+ , La 3+ ) produces the opposite effect [116,117]. In our previous papers [118,119] we have examined theoretically both cases: (i) the substitution of rare-earth Sm ions by non-magnetic trivalent ions (e.g. Y 3+ ) which introduce conduction electrons into the d-conduction band (one electron per dopant) and (ii) the substitution of rare-earth ions by non-magnetic divalent ions (e.g. Sr 2+ ), which play a dilution role and reduce the number of conduction electrons in the d-conduction band (no additional electrons are introduced to the system). We have started our study with the first case which is slightly simpler. In this case the  Figure 55. The average occupancy of f orbitals as a function of x for divalent and trivalent dopants. The one-dimensional case [118]. total number of electrons N = N f + N d is equal to the total number of lattice sites L (the case of one electron per rare-earth atom is considered) and does not depend on the number of dopants X. Consequently, the ground-state energy E g (N f ) = L−N f k=1 λ k is also independent of X and can be calculated directly using our AM for arbitrary N f from the interval [0, L − X]. The results of numerical calculations for the one dimensional lattices of L = 120 and L = 240 sites are summarized in figure 54 (a) for representative model parameters (E f = 0, U = 0.5). It is seen that finite-size effects are negligible and thus these results can be satisfactorily used to represent the behavior of macroscopic systems. From the N f dependence of the ground-state energy one can easily deduce (we remember that N f = 0, 1, . . . , L − X) that E g has a minimum at N 0 f = L/2 for N f L/2 and at N 0 f = L − X for N f < L/2. Thus, the average occupancy of the f -electron orbitals can be finally written as The situation is more complicated for the substitution of rare-earth ions by non-magnetic divalent ions that yield no additional electrons in the d-electron conduction band. In this case the total number of electrons N = L−X as well as the number of d-electrons N d = L−X −N f depend on the number of dopants X and thus the ground-state energy E g (N f , X) = N d k=1 λ k has to be calculated many times for the selected values of X. The numerical results obtained in the one-dimensional case for several selected values of n = N/L are presented in figure 54 (b). The finite-size effects are small again and so the results can be satisfactorily extrapolated to the thermodynamic limit (L → ∞). Comparing these results with the previous case (doping by trivalent ions) one can see a fully different behavior. For small values of N = L − X, the ground-state energy E g (N f , X) has s minimum at N 0 f = 0 and then monotonously increases. This type of behavior holds up to the critical value N c = L − X c ∼ L/2. Above this value E g reaches a minimum at N 0 f > 0 and with increasing N (decreasing X) the position of this minimum shifts to L/2. The resultant behavior  figure 55 we have displayed the dependence of n f on x for doping by trivalent ions. One can see that the substitution of rare-earth ions by trivalent and divalent ions produces altogether different effects. While in the first case the average occupancy of f -orbitals increases, in the second case n f decreases. These results are qualitatively in agreement with experimental measurements [116] of the average samarium valence v = 3 − n f performed on Sm 1−x M x B 6 (M=Y 3+ , La 3+ , Sr 2+ , Yb 2+ ) despite the fact that our calculations have been done at T = 0, while the experiments have been done at room temperatures. However, it is not expected that the increasing temperature could change dramatically this picture. From the theoretical investigation of the spinless Falicov-Kimball model at non-zero temperatures it is known [120] that the effect of temperature is such that a finite temperature only smears the behavior found for T = 0. Thus, one can expect that for non-zero temperatures the accordance of theoretical and experimental results could be even better.  [119]. The experimental results have been obtained for T=300 K (reference [116]).
To verify this conjecture, we have explicitly calculated the dependence of the f -state occupancy n f on the dopant concentration x for nonzero-temperatures using small-cluster exact-diagonalization technique at finite temperatures (this method is described in detail in reference [120]). The results of our numerical calculations are summarized in figure 56 for the selected values of temperature τ . It is seen, at a glance, that the spinless Falicov-Kimball model despite its simplicity can provide a qualitatively correct description of the effects of doping on thermodynamic properties of rare-earth systems for both divalent and trivalent ions. In addition, optimizing the n f (x) behaviour with respect to τ and E f one can obtain a nice quantitative correspondence between the theoretical and experimental results (see figure 57). Since no significant effects of the increasing dimension D (D = 2, 3) on the behaviour of n f (x) has been observed for both zero and non-zero temperatures [118,119], our results can be satisfactorily used to describe the behaviour of real three dimensional materials.

Stability of charge and spin ordering at finite temperatures
Examples of charge and spin ordering discussed in previous sections concerned exclusively the case of T = 0. Real experiments are, however, always performed at finite temperatures, and therefore getting an answer to the question about the stability of zero temperature solutions at finite temperatures seems to be the task of fundamental importance. A positive answer to this question is available at least for the chessboard charge ordering that is the ground state of the conventional Falicov-Kimball model at half-filling for all Coulomb interactions U > 0. For this case, there exists an exact proof [13] of the existence of a phase transition from the low-temperature ordered phase (the chessboard phase) to the high-temperature disordered phase at finite critical temperature τ c (for dimensions D 2) which strongly depends on the local Coulomb interaction U . In addition, the numerical simulations within the grand-canonical Monte-Carlo showed that the phase transitions are of the first order for small and intermediate values of the Coulomb interaction U and of the second order for strong interactions [121,122]. In our paper [123] we have extended the numerical study of the temperature induced phase transitions to the case of phase segregated and striped phases. Moreover, we have considered a more general situation where H FKM is the conventional spinless Falicov-Kimball model (8) and H t is the term of the correlated hopping (13). As was discussed in section 3.3.1, all three above mentioned phases (the chessboard phase, the segregated phase and the axial striped phase) are the ground states of the model Hamiltonian (25) at the symmetric band point, where both N f and N d are fixed to L/2 and, therefore, all numerical calculations at nonzero temperatures are done exclusively in the canonical ensemble. In this formalism, the partition function and the internal energy corresponding to the model Hamiltonian (25) can be written as: In the next step, the summation over all f and d distributions is replaced by the Monte-Carlo summation with the statistical weight e −E/τ /Z.
To identify the transition temperatures from low-temperature ordered phases to high-temperature disordered phase and the type of the phase transition, we have numerically calculated the specific heat C = ( E 2 − E 2 )/(Lτ 2 ), the thermal average of the f -electron occupation w s = w f and the energy distribution P (E). The numerical calculations are performed exclusively at U = 0.5, since the ground-state phase diagram exhibits a richer spectrum of solutions in the weak and intermediate coupling regions in comparison to the strong coupling limit. To verify the capability of our method to describe the phase transitions at finite temperatures we have started with the conventional two-dimensional Falicov-Kimball model (t = 0) at half-filling. As was mentioned above, the physical picture of temperature-induced phase transitions within this relatively simple model is quite understandable at present. For all finite Coulomb interaction U > 0, the ground state of the model is the chessboard phase that persists up to critical temperature τ c (U ), where the system undergoes a phase transition to the homogeneous phase. The phase transition is of the first order for U < 1 and of the second order for U > 1 [121]. Our numerical results obtained within the canonical Monte-Carlo method for C, w s and P (E) fully confirm this picture (see figure 58) [123]. The specific heat curves exhibit a sharp low-temperature peak at τ c ∼ 0.028 that is obviously connected with the phase transition from the chessboard phase to the homogeneous phase, as can be seen from the behaviour of the average f -electron occupation w s for temperatures slightly lower or slightly higher  than τ c . Moreover, the energy distribution function P (E) exhibits an apparent two-peak structure near the critical point τ c (it can be considered as a superposition of two Gaussians), which in accordance with the theory of Challa, Landau and Binder [124] points on the first order phase transition at τ c . Let us now discuss how this picture is changed when the correlated hopping term is added. Firstly, we have examined the case of small values of |t | for which the ground state of the model is still the chessboard phase [68]. The typical examples of C, w s and P (E) from the positive and negative region of t are displayed in figure 59 and in figure 60 for t = −0.3 and t = 0.3. One can see that the correlated hopping term (in the limit of small |t |) does not qualitatively change the picture of temperature induced phase transitions found for t = 0. For both, positive and negative t , there is the first order phase transition from the low-temperature ordered phase to the hightemperature disordered phase, similarly to t = 0, and the only difference between these cases is that the correlated hopping term reduces slightly the critical temperature τ c of the phase transition.
Therefore, in the next step we have turned our attention to the physically much less explored type of configurations, namely, the axial striped configurations that are ground states of the Falicov-Kimball model for the intermediate values of t (|t | ∼ 0.5). Note, that for the axial striped phase even the fundamental question concerning the temperature stability of this phase has remained unanswered so far. This is due to the fact that it is very difficult to find this phase in a pure form. For example, in the conventional Falicov-Kimball model (t = 0), the axial striped phases are stable for a relatively wide range of model parameters [54], but solely in mixtures with other phases (e.g., the empty configuration). In addition, strong finite-size effects have been observed on the stability of these mixtures and therefore it is practically impossible to make any conclusions concerning their stability at finite temperatures from the numerical calculations on finite clusters. However, in the Falicov-Kimball model with correlated hopping, the axial striped phase exists in a pure form for a wide range of model parameters t and U , the finite-size effects on the stability of this phase at τ = 0 are negligible, and so the corresponding numerical study of the temperature stability of the axial striped phase can be performed straightforwardly.
In figure 61 and figure 62 we present our canonical Monte-Carlo results for C, w s and P (E) obtained for two different values of t (t = 0.5 and t = 0.55) from the region where the groundstate of the model is just the axial striped phase. Again, the specific heat curves exhibit a sharp low-temperature peak, the existence of which indicates a phase transition from the axial striped phase to the homogeneous phase. This was independently verified by calculating the average felectron occupation w s and the energy distribution P (E) near the transition point τ c , which clearly  demonstrates the presence of the first order phase transition at τ c . Since the critical temperature τ c of the phase transition for both values of t shifts to smaller values with increasing L, we have performed a detailed finite-size scaling analysis of the τ c (L) dependence to exclude the possibility 42701-58 of τ c vanishing in the thermodynamic limit L → ∞. The resultant τ c (L) dependencies are plotted as insets in figure 61 and in figure 62. It is clearly seen that the critical temperatures τ c for both t = 0.5 and t = 0.55 also persist in the thermodynamic limit, meaning that the axial striped phase remains stable at finite temperatures as well. Moreover, our numerical results show that critical temperatures for an axial striped phase are considerably higher in comparison with critical temperatures for a chessboard phase. We have observed the same behaviour for negative values of t (t = −0.7). However, the critical temperature in this case was only slightly larger than the one corresponding to t = 0.
With increasing t , the half-filled Falicov-Kimball model with correlated hopping exhibits (at τ = 0) a phase transition from the axial striped phase to the segregated phase [68] that takes place at t ∼ 0.6. Since both the chessboard phase and the axial striped phase are insulating while the segregated phase is metallic [68], one can expect a fully different thermodynamic behaviour of the model for the last case. To verify this conjecture we have performed exhaustive numerical studies of the temperature dependence of C, w s and P (E) for t = 1. This study is also important from the point of view that the thermodynamics of the metallic phase has been so far examined only in a few cases [125,126], while for the insulating phase (usually the chessboard phase) there is a number of analytical and numerical results [13,84,127]. The results of our numerical calculations obtained for the specific heat C are shown in figure 63. To reveal the finite-size effects, the calculations for C have been done on several different clusters of L=6 × 6, 8 × 8, 10 × 10, 12 × 12 and 16 × 16 sites. We have found that the specific heat curves, in the low-temperature region, strongly depends on the cluster sizes, and, therefore, a very careful analysis has to be performed to find the correct behaviour of the model in the thermodynamic limit L → ∞. On small finite clusters (L = 6 × 6 and L = 8 × 8), the specific heat exhibits only one-peak structure in the low-temperature region (τ ∼ 0.15). With the increasing cluster size L, an additional peak is stabilized at slightly higher temperatures (τ ∼ 0.23), while the first peak is gradually suppressed and probably fully disappears in the thermodynamic limit. The behaviour of the average f -electron occupation shows (see figure 63) that the second peak in the specific heat corresponds to the phase transition from the low-temperature ordered (segregated) phase to the high-temperature disordered phase.
The nature of this phase transition is, however, different in comparison to the previous cases. While the energy distribution function P (E) is double peaked for the chessboard and the axial striped phase near the transition temperature τ c (the first order phase transition), P (E) exhibits a single-peak structure for the segregated phase, which points to the second order phase transition at τ c . Comparing the thermodynamic behaviour of the model in the chessboard, axial striped and segregated regions one can find two other important differences, namely, (i) the critical temperature of the second order phase transition is approximately ten times higher than critical temperatures of the first order phase transitions, and (ii) the specific heat (in the low-temperature region) exponentially decreases for the chessboard and axial striped phase, while in the segregated phase the specific heat C(τ ) seems to show a linear behaviour indicating the Fermi-liquid behaviour for τ < 0.08 (see the inset in figure 63 (a)). The observation of the linear contribution to the specific heat in the low-temperature region (τ < 0.08) is consistent with the behaviour of the average f -electron occupation in this region (see figure 63 (c)). One can see that despite the increasing  [128].

42701-60
temperature (from 0 to 0.08) the f -electrons preferably occupy only one half of the lattice leaving another part empty. Due to the on-site Coulomb interaction between the f and d electrons, the itinerant d electrons preferably occupy the empty part of the lattice, where they can move as free particles yielding the linear contribution to the specific heat. We have performed the same numerical study for the spin-1/2 Falicov-Kimball model extended by the spin-dependent Coulomb interaction J between the localized f and itinerant d electrons as well as the on-site Coulomb interaction (U f f ) between the localized f -electrons [128]. The coexistence of charge and spin degrees of freedom for such type of a system evokes the question concerning a possible sequence of two phase transitions corresponding to the breaking of different types of charge/spin ordering. We have tested the possibility of such a scenario independently for two configuration types displayed in figure 64, which are the ground states of the generalized Falicov-Kimball model for U = 2, J = 1, U f f = 4, E f = 0 and N = 2L or U = 4, J = 0.8, U f f = 8, E f = −2 and N = 3L/2. To identify the type of phase transition, we have used the the specific heat and the structure factor of the charge and spin ordering defined by [128]: and The temperature behaviours of these quantities are displayed in figure 65. It is seen that in the first case both types of ordering are simultaneously destroyed at the same critical temperature τ c ∼ 0.05, while in the second case the spin ordering disappears already at very low temperatures (τ c < 0.01), and the charge ordering persists up to τ c ∼ 0.25 5 .

Conclusion
In this review we have presented the results of our theoretical study of charge and spin ordering in strongly correlated electron systems obtained within various generalizations of the Falicov-Kimball model. The primary goal of this study was to identify crucial interactions that lead to the stabilization of various types of charge ordering in these systems such as the axial striped ordering, diagonal striped ordering, phase separated ordering, phase segregated ordering, etc. Among the major interactions that come into account, we have examined the effect of local Coulomb interaction between localized and itinerant electrons, long-range and correlated hopping of itinerant electrons, long-range Coulomb interaction between localized and itinerant electrons, local Coulomb interaction between itinerant electrons, local Coulomb interaction between localized electrons, spin-dependent interaction between localized and itinerant electrons, both for zero and nonzero temperatures, as well as for doped and undoped systems.
We have started our study of charge ordering in the strongly correlated systems with the spinless Falicov-Kimball model. First, we have focused on a description of the ground-state properties of the model at T = 0. Since the previous theoretical studies have shown that the ground-state properties of the Falicov-Kimball model are very sensitive to the type of approximation used, to study the charge ordering, valence and metal-insulator transitions in the Falicov-Kimball model we have used the EDM on finite clusters with subsequent extrapolation of the results to the thermodynamic limit (infinitely large system). Our results showed that the spinless Falicov-Kimball model is capable of describing homogeneous as well as phase-separated distributions of f electrons and thereby discontinuous valence and metal-insulator transitions induced by pressure. The major driving interaction of these transitions is the Coulomb interaction between itinerant and localized electrons.
In order to elaborate the most realistic description of charge ordering and valence and metalinsulator transitions in real materials (rare-earth and transition-metal compounds), we have generalized the original type of electron hopping (only to the nearest neighbours) to a much more realistic type of hopping (the long-range hopping with power decreasing hopping amplitudes) and in addition we have taken into account the term of correlated hopping, which modifies the hopping amplitudes of itinerant electrons from one site to another according to the occupancy of these sites by localized electrons. We have found that the term of correlated hopping has a significant effect on the dynamics of valence and metal-insulator transitions in the spinless Falicov-Kimball model and should, therefore, be taken into account for the correct description of these phenomena. Moreover, we have shown that the correlated hopping term leads to the stabilization of the axial charge stripes, suggesting an alternate and a much simpler mechanism (than the one which was previously considered within the Hubbard and t − J model) of forming an inhomogeneous charge ordering in strongly correlated systems. Regarding the effect of long-range hopping term on metal-insulator transitions, we have found that it is not very significant, at least for q 0.3, which physically turns out to be the most interesting case.
The next step towards a comprehensive description of cooperative phenomena in the Falicov-Kimball model was to study the model for dimensions D = 2 and D = 3. To fulfill this task, it was necessary to develop a new numerical method, since the EDM is capable of analysing only clusters up to L ∼ 40 sites, which does not offer any possibility to extrapolate the results to the thermodynamic limit (L = ∞). This goal has been fully fulfilled. Our numerical method (based on a modification of an exact diagonalization algorithm) is very accurate, even for cluster sizes of several hundred sites (this has been tested in the limiting cases where the exact solutions are known). Moreover, it is very flexible, and so it permits to study various generalizations of the model. The study of the Falicov-Kimball model by this method showed that the basic picture of charge ordering and valence and metal-insulator transitions remains unchanged for higher dimensions. The only important difference is that the stability area of metallic phases shifts to higher values of the Coulomb interaction U . This fact is very important from the practical point of view, since real materials are just in this limit.
The fact that we were able to analyse the model for cluster sizes of several hundred lattice sites allowed us to study the formation of charge ordering in two-and three-dimensional systems for arbitrary values of U , which was very valuable, since the previous results were limited either to the area of U 1 or U 1, or were limited to the type of the investigated configurations (e.g., the periodic configurations with the small periods and their mixtures). We have found that the spinless Falicov-Kimball model is capable (in D = 2 and D = 3) of describing a very wide range of charge orderings involving periodic charge ordering, phase-segregated and phase-separated ordering as well as axial and diagonal striped ordering. In the the three-dimensional case, our results represent the first and so far the only attempt to systematically describe the formation of inhomogeneous charge ordering in the spinless Falicov-Kimball model, which was to a large extent possible due to a new elaborated numerical method.
Moreover, the inclusion of spins and the spin interactions of the Hubbard type between itinerant electrons and of the Ising type between itinerant and localized electrons clearly demonstrated the enormous potential of the model in describing different types of charge and spin superstructures, including the segregated phase and phase-separated orderings, as well as axial and diagonal striped (band) orderings with ferro-, ferri-, or antiferromagnetic ground state. Furthermore, the spin model in the presence of interband Coulomb interaction between localized and itinerant electrons maintains advantages of a spinless model, namely, the capability of describing valence and metal-insulator transitions induced by pressure, temperature and alloying.