Evolution of the spectrum of the Hubbard model with filling

The diagram technique for the one-band Hubbard model is formulated for the case of moderate to strong Hubbard repulsion. The expansion in powers of the hopping constant is expressed in terms of site cumulants of electron creation and annihilation operators. For Green's function an equation of the Larkin type is derived and solved in a one-loop approximation for the case of two dimensions and nearest-neighbor hopping. With decreasing the electron concentration in addition to the four bands observed at half-filling a narrow band arises near the Fermi level. The dispersion of the new band, its bandwidth and the variation with filling are close to those of the spin-polaron band in the t-J model.


Introduction
Systems with strong electron correlations which, in particular, cuprate perovskites belong to are characterized by a Coulomb repulsion that is comparable to or larger than hopping constants. One of the simplest models for the description of such systems is the Hubbard model [ 1,2,3] with the Hamiltonian where t nm is the hopping constants, the operator a † nσ creates an electron on the site n with the spin projection σ = ±1, U is the on-site Coulomb repulsion and the electron number operator n nσ = a † nσ a nσ . In the case of strong electron correlations, U ≥ |t nm |, it is reasonable to use a perturbation expansion in powers of the hopping constants for investigating Hamiltonian (1). Apparently the first such expansion was considered in reference [ 4]. The further development of this approach was given in references [ 5,6,7,8,9] where the diagram technique for Hubbard operators was developed and used for investigating the Mott transition, the magnetic phase diagram and the superconducting transition in the Hubbard model.
The rules of the diagram technique for Hubbard operators are rather intricate. Besides, these rules and the graphical representation of the expansion vary depending on the choice of the operator precedence. The diagram technique suggested in references [ 10,11,12,13] is free from these defects. In this approach the power expansion for Green's function of electron operators a nσ and a † nσ is considered and the terms of the expansion are expressed as site cumulants of these operators. Based on this diagram technique the equations of the Larkin type [ 15] for Green's function were derived [ 10,12,13].
However, the application of this approach runs into problems. In particular, at half-filling the spectral weight obtained after a resummation of diagrams appears to be negative near frequencies ω d = ±U/2 [ 13]. This drawback is connected with divergencies in cumulants at these frequencies [ 14]. As can be seen from formulas given below, all higher-order cumulants have such divergencies at ω d with sign-changing residues, which are expected to compensate the negative spectral weight in the entire series. On the other hand, at frequencies neighboring to ω d cumulants are regular. If a selected subset of diagrams is expected to give a correct estimate of the entire series for these frequencies the values at ω d can be corrected by an interpolation using results for the regular regions. This procedure was applied in reference [ 14] for the case of two dimensions, half-filling, nearest-neighbor hopping and with the use of the one-loop approximation. The spectrum was shown to consist of four bands. These band structure and the calculated shapes of the electron spectral function are close to those obtained in the Monte-Carlo [ 16,17,18], cluster perturbation [ 19] and the two-particle self-consistent [ 20] calculations. The four-band structure of the spectrum of the Hubbard model was also considered in reference [ 21].
In the approach of reference [ 14] the mentioned four-band structure of the spectrum at half-filling has its origin in the regions of large damping which separate the low-and high-frequency bands. In the major part of the Brillouin zone both these types of the bands owe their origin to the same terms of the irreducible part, i.e. to the same interaction processes. This brings up the questions: How this spectrum is changed with filling and how the quasiparticle peak, which determines the photoemission leading edge, arises in the spectrum? As expected, this peak differs in nature from other spectral features. In the present work it is shown that on certain deviation from half-filling an additional narrow band arises near the Fermi level on the background of the four above-mentioned bands. The dispersion of the new band, its bandwidth and the variation with filling are close to those of the spin-polaron band in the t-J model [ 22,23].
In the following section the perturbation expansion for the electron Green's function is formulated in the form convenient for calculations and Larkin's equation is derived. In Sec. III the equations of the previous section are used for calculating the spectral function and the obtained results are discussed. Concluding remarks are presented in Sec. IV.

Diagram technique
We consider Green's function where the angular brackets denote the statistical averaging with the Hamiltonian H = H − µ nσ n nσ , µ is the chemical potential, T is the time-ordering operator which arranges other operators from right to left in ascending order of times τ , a nσ (τ ) = exp(Hτ )a nσ exp(−Hτ ) andā nσ (τ ) = exp(Hτ )a † nσ exp(−Hτ ). Choosing as the unperturbed Hamiltonian and the perturbation, respectively, and using the known expansion [ 24] for the evolution operator we get where β = T −1 is the inverse temperature, the subscript "0" near the angular bracket indicates that the averaging and time dependencies of operators are determined with the Hamiltonian H 0 . The subscript "c" indicates that terms which split into two or more disconnected averages have to be dropped out.
The Hamiltonian H 0 is diagonal in the site representation. Therefore the average in the right-hand side of equation (4) splits into averages belonging to separate lattice sites. To be nonzero these latter averages have to contain equal number of creation and annihilation operators. Let us consider the average which appears in the first order: Tā n ′ σ (τ ′ )a nσ (τ )ā n ′ 1 σ 1 (τ 1 )a n 1 σ 1 (τ 1 ) (for short here and below subscripts "0" and "c" are omitted). For this average to be nonvanishing operators have to belong either to the same site or to two different sites, The multiplier 1 − δ nn 1 ensures that n = n 1 in the two last terms (the case n = n 1 is taken into account by the first term in the right-hand side). After the rearrangement of the terms we find [ 25] of the first and second order, respectively. All operators in cumulants (5) belong to the same lattice site. Due to the translational symmetry of H 0 in equation (3) the cumulants do not depend on the site index which is therefore omitted in equation (5). Averages which appear in higher-order terms of equation (4) can be transformed in the same manner. The average in the k-th order term which contains k + 1 creation and annihilation operators is represented by the sum of all possible products of cumulants with the sum of orders equal to k + 1. All possible distributions of operators between the cumulants in these products have to be taken into account. The sign of a term in this sum is determined by the number of permutations of fermion operators, which bring the sequence of operators in the initial average to that in the term. Actually the above statements determine rules of the diagram technique. Additionally we have to take into account the presence of topologically equivalent terms -terms which differ only by permutation of operators H 1 (τ i ) in equation (4). Since these terms are equal, in the expansion only one of them can be taken into account with the prefactor ν = j/k! where j is the number of topologically equivalent terms. Following reference [ 11] in diagrams we denote a cumulant by a circle and the hopping constant t nn ′ by a line directed from n ′ to n. The external operatorsā n ′ σ (τ ′ ) and a nσ (τ ) are denoted by directed lines leaving from and entering into a cumulant. The order of a cumulant is equal to a number of incoming or outgoing lines. Summations and integrations over the internal indices n i , n ′ i , σ i and τ i are carried out. Since site indices of operators included in a cumulant coincide, some site summations disappear. Also some summations over σ i get lost, because in any cumulant spin indices of creation and annihilation operators have to match. Taking into account the multiplier (−1) k in equation (4), the sign of the diagram is equal to (−1) l where l is the number of loops formed by hopping lines. Figure 1 demonstrates connected diagrams of the first four orders of the power expansion (4) with their signs and prefactors. Here the thick line with arrow in the left-hand side of the equation is the total Green's function. Notice that if we set t nn = 0, contributions of the diagrams (b), (d)-(f), (i)-(n), (p), (t), and (u) vanish. However, below a renormalized hopping parameter will be introduced which is nonzero for coinciding site indices and therefore the mentioned diagrams are retained in figure 1.
All diagrams can be separated into two categories -those which can be divided into two parts by cutting some hopping line and those which cannot be divided in this way [ 8,15]. These latter diagrams are referred to as irreducible diagrams. In figure 1 the diagrams (c), (e), (f), (h), (j), (k), (n), (p), and (r)-(v) belong to the first category, while the others are from the second category. The former diagrams can be constructed from the irreducible diagrams by connecting them with hopping lines. Notice that a prefactor ν of some composite diagram is the product of prefactors of its irreducible parts [cf., e.g., the irreducible diagrams (a), (d) and the composed diagrams (j), (k) in figure 1]. If we denote the sum of all irreducible diagrams as K(n ′ τ ′ , nτ ), the equation for Green's function can be written as This equation has the form of the Larkin equation [ 15]. The partial summation can be carried out in the hopping lines of cumulants by inserting the irreducible diagrams into these lines. In doing so the hopping constant t nn ′ in the respective formulas is substituted by For the diagram (b) in figure 1 inserting irreducible diagrams into the hopping line leads to the diagrams (g), (m), and (q). Due to the translation and time invariance of the problem the quantities in equations (6) and (7) depend only on differences of site indices and times. The use of the Fourier transformation where ω l = (2l + 1)πT is the Matsubara frequency, simplifies significantly these equations: In the approximation used below the total collection of irreducible diagrams K(k, iω l ) is substituted by the sum of the two simplest diagrams (a) and (b) appearing in the first two orders of the perturbation theory. Due to the form of the latter diagram this approximation is referred to as the one-loop approximation. In the diagram (b) the hopping line is renormalized in accordance with equation (9). Thus, where K 1 (iω l ) and are the Fourier transforms of cumulants (5), N is the number of sites and we set k t k = 0. Notice that in this approximation K does not depend on momentum. Now we need to calculate the cumulants in equation (10). To do this it is convenient to introduce the Hubbard operators X ij n = |in jn| where |in are eigenvectors of site Hamiltonians forming H 0 , equation (3). For each site there are four such states: the empty state |0n with the energy E 0 = 0, the two degenerate singly occupied states |σn with the energy E 1 = −µ and the doubly occupied state |2n with the energy E 2 = U − 2µ. The Hubbard operators are connected by the relations with the creation and annihilation operators. The commutation relations for the Hubbard operators are easily derived from their definition. Using equation (11) the first cumulant in equation (5) can be computed straightforwardly: where Z 0 = e −βE 0 + 2e −βEσ + e −βE 2 is the partition function, E ij = E i − E j . As indicated in [ 8,10,13], if K(k, iω l ) is approximated by this cumulant the result corresponds to the Hubbard-I approximation [ 2]. To calculate K 2 it is convenient to use Wick's theorem for Hubbard operators [ 5,6,7,8]: where the averaging and time dependencies of operators are determined by the Hamiltonian H 0 , α is the index combining the state and site indices of the Hubbard operator. If X α is a fermion operator (X 0σ , X σ2 and their conjugates), P k is the number of permutation with other fermion operators which is necessary to transfer the operator X α from its position in the left-hand side of equation (13) to the position in the right-hand side. In this case where i and j are the state indices of X α . If X α is a boson operator (X 00 , X 22 , X σσ ′ , X 02 , and X 20 ), P k = 0 and In equation (13), [X α k , X α ] ± denotes an anticommutator when both operators are of fermion type and a commutator in other cases. The substitution of equation (11) in Tā σ (τ ′ )a σ (τ )ā σ 1 (τ ′ 1 )a σ 1 (τ 1 ) in K 2 , equation (5), leads to six nonvanishing averages of Hubbard operators (such averages are nonzero if the numbers of the operators X 0σ and X σ0 coincide in them, and the same is true for the pair X σ2 and X 2σ ). Applying Wick's theorem (13) to these averages the number of operators in them is sequentially decreased until only timeindependent operators are left. For H 0 in equation (3) these are X 00 , X σσ ′ , and X 22 . Their averages are easily calculated. As a result, after some algebra we find 0 U e −βE 0 g 0σ (iω l )g 0σ (iω l 1 )g 02 (iω l + iω l 1 ) g 0σ (iω l ) + g 0σ (iω l 1 ) +e −βE 2 g σ2 (iω l )g σ2 (iω l 1 )g 02 (iω l + iω l 1 ) g σ2 (iω l ) + g σ2 (iω l 1 ) where g ij (iω l ) = (iω l + E ij ) −1 is the Fourier transform of functions (14) and (15). Equations (12) and (16) can be significantly simplified for the case of principal interest U ≫ T . In this case if µ satisfies the conditions λ ≫ T , the exponent exp(−βE 1 ) is much larger than exp(−βE 0 ) and exp(−βE 2 ). By passing to real frequencies we can ascertain that terms in (16) with the two latter multipliers contain the same peculiarities as other terms. Therefore terms with these multipliers can be omitted in the equations and we get Let us turn to real frequencies by substituting iω l with z = ω + iη where η is a small positive constant which affords an artificial broadening. Results given in the next section were calculated with G(k, z) in the right-hand side of equation (10) taken from the Hubbard-I approximation. As mentioned, this Green's function is obtained if K(kω) in equation (8) is approximated by K 1 from equation (18) which gives Below the two-dimensional square lattice is considered. It is supposed that only the hopping constants between nearest neighbor sites t are nonzero which gives t k = 2t[cos(k x ) + cos(k y )] where the intersite distance is taken as the unit of length. Due to the electron-hole symmetry of Hamiltonian (1) the consideration can be restricted to the range of the chemical potentials µ ≤ U 2 .  (10), (18) and (19). The change to real frequencies carried out in the previous section converts the Matsubara function (2) into the retarded Green's function [ 24]. It is an analytic function in the upper half-plane which requires that ℑK(ω) be negative. As seen from figure 2, this condition is violated at ω d = −µ and U − µ. This difficulty of the considered approximation was indicated in reference [ 13]. The problem is connected with divergencies at these frequencies introduced by functions g 0σ (ω) and g σ2 (ω) in the above formulas. As can be seen from the procedure of calculating the cumulants in the previous section, these functions and divergencies with signchanging residues appear in all orders of the perturbation expansion (4). In the entire series the divergencies are expected to compensate each other so that the resulting ℑK(ω) is negative everywhere. However, in the considered subset of terms such compensation does not occur. Nevertheless, as seen from figure 2, at frequencies neighboring to ω d cumulants are regular and, if the used subset of diagrams is expected to give a correct estimate of the entire series for these frequencies, the value of ℑK(ω) at the singular frequencies can be reconstructed using an interpolation and its values in the regular region. Examples of such interpolation are given in figure 2.

Spectral function
As seen from figure 2a, at half-filling, µ = U 2 , ℑK(ω) has two broad minima. With decreasing the chemical potential from this value these minima shift with respect to the Fermi level without a noticeable change of their shapes until the Fermi level enters one of the minima at µ ≈ 0.17U. As this takes place, two new sharp minima arise near frequencies −µ and U − µ on the background of the above-mentioned broad minima. The appearance of the broad features in figure 2 is connected with the third term in σ 1 K 2 in equation (18), while the sharp minima are related to the second term in this formula. Its contribution to K(ω), equation (10), grows rapidly when the Fermi level enters the broad minimum.
The function K(z) has to be analytic in the upper half-plane also and therefore its real part can be calculated from its imaginary part using the Kramers-Kronig relations. We use this way with the interpolated ℑK(ω) to avoid the influence of the divergencies on ℜK(ω). However, the use of the interpolation overrates somewhat values of |ℑK(ω)| which leads to the overestimation of the tails in the real part. To correct this defect the interpolated K(ω) is scaled so that in the far tails its real part coincides with the values obtained from equation (10).
The spectral function obtained in this way for momenta along the symmetry lines of the square Brillouin zone is shown in figure 3. The shapes of the spectral function in figure 3a are nearly the same as at half-filling -with decreasing µ from 0.5U to approximately 0.17U these curves shift with respect to the Fermi level without perceptible changes in their shapes. As indicated in reference [ 14], four bands can be distinguished in these  spectra. For parameters of figure 3a these bands are located near frequencies −4|t|, |t|, 4|t|, and 9|t| (see figure 3b). For the major part of the Brillouin zone the peaks forming the bands arise at frequencies which satisfy the equation 1 − t k ℜK(ω) = 0 and fall into the region of a small damping |ℑK(ω)| [see equation (20)]. As seen from figure 2, such regions of small damping are located between and on the outside of the two broad minima in ℑK(ω). This is the reason of the existence of the four well separated bands -two of them are located between the minima of ℑK(ω), while two others are on the outside of these minima. Broader maxima of A(kω) for momenta near the boundary of the magnetic Brillouin zone are of different nature -since t k is small for such momenta, the resonant denominator in equation 20 does not vanish and the shape of the spectral functions is determined by ℑK(ω) in the numerator of this formula.
More substantial changes in A(kω) occur for µ ≤ 0.17U when the Fermi level enters one of the broad minima in ℑK(ω). As seen from figure 3c, in addition to the mentioned four bands there appear sharp dispersive features near ω = −µ and U −µ for momenta in the vicinity of the boundary of the magnetic Brillouin zone. It is clear that these changes in the spectral function are connected with the changes in ℑK(ω) shown in figure 2. The peaks near −µ are more intensive than those near U − µ and are located in the nearest vicinity of the Fermi level. For µ ≈ 0.17U the width of the band formed by the former peaks is comparable to the superexchange constant J = 4t 2 /U of the effective Heisenberg model which describes magnetic excitations in the limit U ≫ |t|. This indicates obviously the participation of the spin excitations in the formation of these band states. The bandwidth decreases with further reduction of the electron concentration. As this takes place, the peak intensities first grow and then saturate. The maximum energies of the band are located near the boundary of the magnetic Brillouin zone. For parameters of figure 3c for these momenta the band touches the Fermi level and the corresponding peaks disappear above it (see figure 3d).
These properties of the band resembles those of the spin-polaron band in the t-J model. This latter band is also located near the Fermi level, has the similar dispersion and the bandwidth, which decreases with the reduction of the electron concentration (with the increase of the hole doping counted from half-filling) [ 23]. However, there is one essential difference in the behavior of these bands in the Hubbard and t-J models. In the former model the narrow band near the Fermi level appears at a certain deficiency of electrons, while the above-mentioned four bands exist in the entire considered range of electron concentrations. In the t-J model the situation is opposite -the spin-polaron band exist in the wide range of hole concentrations 0 ≤ x ≤ 0.17, while the wider band -an analog of the four-band structure -starts to form at x ≈ 0.06 [ 23].
The value of the chemical potential for which the Fermi level enters the minimum of ℑK(ω) and the narrow band begins to form depends on the ratio t/U. For example, for t = −U/4 this happens at µ ≈ 0.27U.
There is some difficulty in the comparison of our results with the data of Monte-Carlo calculations and cluster theories. It is connected with the fact mentioned in reference [ 14]: the one-loop approximation overestimates the spectral weights of the two internal bands of the four-band structure near the momenta (0, 0) and (π, π). As a consequence the electron concentration calculated from the formula  [ 18]) it can be concluded that these concentrations have to be approximately 0.95 and 0.85, respectively. With this in mind we find that the spectral functions and dispersions in figure 3 are close to the Monte-Carlo spectra (cf. with figures 10 and 11 in [ 18]). In these latter spectra for some deviation from half-filling a weakly dispersive feature is also observed near the Fermi level. However, this feature has low intensity and is lost at the foot of a more intensive maximum on approaching the boundary of the magnetic Brillouin zone. These differences may be connected with the comparative high temperature T = 0.33|t| used in the Monte-Carlo simulations. A similar weakly dispersive feature near the Fermi level was obtained also in quantum cluster theories, however for much smaller deviations from half-filling (cf. with figure 2c in [ 27] and figure 32 in [ 28]).

Conclusion
The considered diagram technique is very promising for a generalization to manyband Hubbard models for which energy parameters of the one-site parts of the Hamiltonians exceed or at least are comparable to the intersite parameters. The expansion in powers of these latter parameters can be expressed in terms of cumulants in the same manner as discussed above. Now there are distinct cumulants which belong to different site states characterized by dissimilar parameters of the repulsion and the level energy. These cumulants are described by formulas similar to equations (12) and (16). For example, for the Emery model [ 29] which describes oxygen 2p σ and copper 3d x 2 −y 2 orbitals of Cu-O planes in high-T c superconductors there are two types of cumulants corresponding to these states. Diagrams of the lowest orders, e.g., for Green's function on copper sites resemble those shown in figure 1 where "oxygen" cumulants are included in hopping lines. Equations of the type of (8) and partial summations similar to equation (9) can be derived also in this case.
In summary, the diagram technique for the one-band Hubbard model was formulated for the case of moderate to strong Hubbard repulsion. The expansion in powers of the hopping constant is expressed in terms of site cumulants of electron creation and annihilation operators. For Green's function the equation of the Larkin type was derived and solved for the case of two dimensions and nearest-neighbor hopping. With decreasing the electron concentration in addition to the four bands observed at half-filling a narrow band arises near the Fermi level for momenta near the boundary of the magnetic Brillouin zone. On the occurrence the width of the band is comparable to the superexchange constant J = 4t 2 /U which indicates the participation of the spin excitations in the band formation. The bandwidth decreases with decreasing the electron concentration. The maximum energies of the band are located near the boundary of the magnetic Brillouin zone. For some deviation from half-filling in these points the band touch the Fermi level. By these properties the band resembles the spin-polaron band of the t-J model.