The composite operator method route to the 2D Hubbard model and the cuprates

In this review paper, we illustrate a possible route to obtain a reliable solution of the 2D Hubbard model and an explanation for some of the unconventional behaviours of underdoped high-$T_\text{c}$ cuprate superconductors within the framework of the composite operator method. The latter is described exhaustively in its fundamental philosophy, various ingredients and robust machinery to clarify the reasons behind its successful applications to many diverse strongly correlated systems, controversial phenomenologies and puzzling materials.

The composite operator method (COM) has been designed and devised and is still currently developed with the aim of providing an analytical (operatorial) method that would seamlessly account for the natural emergence in SCES of elementary excitations, that is quasi-particles, whose operatorial description can only be realized in terms of fields whose commutation relations are inherently non-canonical. Such a philosophy requires a complete rethinking, more than a mere rewriting, of the whole apparatus of the GF framework and of the equations of motion (EM) formalism, which are at the basis of the operatorial approach. As a matter of fact, the standard ingredients of the GF framework (spectral weights, dispersions, decay rates, . . . ) need to be partly reinterpreted in order to properly and effectively describe the system under analysis. Moreover, new concepts -and consequently a new terminology -emerge naturally when the EM formalism is applied to such non-canonical fields.
The presence in the Hamiltonian describing the system under analysis of strongly correlated termsterms not quadratic in the canonical fields {a λ , b µ , . . .} describing the (quasi-)particles building up the system (electrons, phonons, magnons, . . . ) and their internal degrees of freedom λ = {λ i } (spin, orbital, site, momentum, branch, . . . ) -immediately and naturally leads to the emergence, in the EM of the canonical fields, of products of the same canonical fields which is not possible to recast as linear combinations of the canonical fields themselves. Such products of canonical fields are what we define composite operators (CO): they denote the very essence of strong correlations and operatorially describe the possible excitations in the system. Such excitations embody the strong interplay of different degrees of freedom and/or of different families of canonical fields characterizing such systems. A complete basis for these new operators (a set of operators in terms of which is possible to express all CO emerging at the first order in the EM hierarchy as linear combinations) contains all possible products of the number . .} operators with the canonical fields themselves. For instance, the complete basis for a single site, single orbital fermionic system is c ↑ , c ↓ , c † ↓ c ↓ c ↑ , c † ↑ c ↑ c ↓ ; the latter two operators being the main components of the well-known Hubbard operators: η σ = nσc σ , ξ σ = (1 − nσ) c σ . In principle, at each order of the EM hierarchy, to find the complete basis to express all possibly emerging new CO as linear combinations, one more layer of number and mixing operators is necessary in front of each element of the previous complete basis. The procedure to construct a complete basis at each order of the EM hierarchy is heavily affected by the Pauli principle for purely fermionic systems, that is by the algebraic properties of the canonical fields. For instance, in the previous complete basis we have no c † ↑ c ↑ c ↑ because it is simply zero or we have that c † ↑ c ↓ c ↑ = −c † ↑ c ↑ c ↓ and no mixing operators seemed to be present while the fictitious absence is simply related to the fact that the charge and spin degrees of freedom on the same site are just indissolubly linked. This is so very true that such algebraic links are at the basis of a further fundamental ingredient of the COM [the Algebra Constraints (AC)] and unfortunately they are often overlooked leading to doomed-to-fail approximations.
Only in purely fermionic systems with a finite number of degrees of freedom, where the process of adding layers is limited again by the Pauli principle (e.g., nσnσc σ = nσc σ ), it is possible to close selfconsistently the system of EM for any possible Hamiltonian, that is for any possible strongly correlated term. On infinite systems, the possibility to close the hierarchy of EM is related to (i) the effective expression of the strongly correlated terms and (ii) if they interest those degrees of freedom which are truly infinite (e.g., sites and momenta for a bulk system). For instance, for the Hubbard model on the bulk, it is not possible to close the hierarchy of the EM because the degree of freedom in which the hopping term is diagonal, the momentum, is infinite and mixed by the Hubbard repulsion, which is instead diagonal in the infinite site space where the hopping term mixes. This is just the very magic of the Hubbard model! Then, when it is not possible to close the hierarchy of the EM, one can always choose where to stop being sure that the processes up to those defined by the CO taken into account will be correctly described in such an approximation. However, this is only one necessary step towards a complete solution because to compute the GF G and the correlators C of such CO under the effect of the Hamiltonian H it is necessary to calculate explicitly their weights and overlaps appearing in the normalization matrix I and the connections among them appearing in the matrix m [4,46]. This will allow one to obtain the eigenenergies E (eigenvalues of the energy matrix ε) and the spectral weights σ (from the matrix Ω of eigenvectors of ε) where ψ a belongs to the chosen set of CO (ψ 1 , . . . , ψ n ), which acts as the basis in the operatorial space, J a is the current of ψ a , δJ a is the residual current of ψ a defined by {δJ a , ψ † b } = 0, Σ ab is the residual self-energy, A ab is the spectral density and f F (ω) is the Fermi function. F is the Fourier transform from time to frequency domain. The latter expression on the right for G ab is exact if the hierarchy of EM closes on a finite number n of CO, otherwise it is the approximate expression of G in the n-pole approximation. This is not the very end of the story since I and m usually contains parameters to be computed. Some of these parameters are just correlators C ab of the basis that can be self-consistently computed through their expression in terms of the GF G (the fluctuation-dissipation theorem), though some others are higher-order correlators (correlators of CO appearing at higher orders in the hierarchy of the EM) and need to be computed in some way. COM uses the AC dictated by the algebra obeyed by the CO in the basis, which leads to relationships between correlators of the basis (e.g., ξ σ η † σ = 0), to fix such parameters achieving a two-fold result: (i) enforce in the solution the AC that are not automatically satisfied contrarily to what many people erroneously think [45,47] and (ii) compute all unknowns in the theory.
There is still something to be computed (or neglected): the residual self-energy Σ. According to the physical properties under analysis and the range of temperatures, dopings, and interactions that we want to explore, we have to choose an approximation to compute the residual self-energy. In the last years, we have been using the n−pole approximation , the asymptotic field approach [72][73][74] and the non-crossing approximation (NCA) [15,[75][76][77][78].

2D Hubbard model
The Hamiltonian of the single-orbital 2D Hubbard model reads as is the electronic field operator in spinorial notation and Heisenberg picture (i = (i, t i )). · and ⊗ stand for the inner (scalar) and the outer products, respectively, in spin space. i is a vector of the two-dimensional square Bravais lattice, n σ (i) = c † σ (i) c σ (i) is the particle density operator for spin σ at site i, n (i) = σ n σ (i) = c † (i) · c (i) is the total particle density operator at site i, µ is the chemical potential, t is the hopping integral and the energy unit hereafter, U is the Coulomb on-site repulsion and α ij is the projector on the nearest-neighbor sites with Fourier transform . For any operator ψ (i), we use the notation ψ κ (i) = j κ ij ψ (j, t i ) where κ ij can be any function of the two sites i and j.

Basis and equations of motion
Following the COM prescription [15,45,46,91], we have chosen a basic field and, in particular, we have selected the following composite triplet field operator ψ , the Hubbard operators, satisfy the following EM: where the higher-order composite field π (i) The use of c s (i) will highlight the most relevant physics as we do expect spin fluctuations to be the most relevant ones. For further details, such as the EM satisfied by the field c s (i), see reference [46].

Normalization I matrix
In the paramagnetic and homogeneous case, the entries of the normalization matrix I(k) = F k {ψ (i, t) , ψ † (j, t)} have the following expressions where n = n (i) is the filling, χ κ s = 1 3 n κ k (i) n k (i) is the spin-spin correlation function at distances determined by the projector κ and f s = 1 3 c † (i) · σ k · c α (i) n α k (i) is a higher-order (up to three different sites are involved) spin-spin correlation function. We have also introduced the following definitions, which are based on those related to the correlation functions of the fields of the basis: where no summation over σ is intended. β (k) and η (k) are the projectors onto the second-nearest-neighbor sites along the main diagonals and the main axes of the lattice, respectively.

m-matrix
In the same conditions, the entries of the matrix where ∆ = C α ξ ξ − C α ηη is the difference between upper and lower intra-Hubbard-subband contributions to the kinetic energy and p = 1 We can avoid cumbersome and somewhat meaningless calculations (they will lead to the appearance of many unknown higher-order correlation functions) by restricting m 33 (k) just to the local and the nearest-neighbor terms:m 0 33 andm α 33 . Given the overall choice of cutting cubic harmonics higher than the nearest-neighbor ones, for the sake of consistency, we also neglected the β (k) and η (k) terms in I 33 (k). We checked that this latter simplification does not lead to any appreciable difference: within the explored paramagnetic solution, χ β s and χ η s have not very significative values. By checking systematically all operatorial relations existing among the fields of the basis, we can recognize the following AC where D = n ↑ (i) n ↓ (i) is the double occupancy. These relations lead to the following very relevant ones: On the other hand, we can compute χ α 0 , χ α s , χ α p and f s by operatorial projection, which is equivalent to the well-established one-loop approximation [91][92][93] for same-time correlations functions As a matter of fact, the energy matrix ε(k) = m(k)I −1 (k) is sure to have real eigenvalues if the normalization matrix I (k) is semi-positive.1 Then, the presence of χ α s and f s in the normalization matrix I (k) imposes a special care in evaluating their values. Accordingly, we have decided to avoid using AC to fix them, and to fix χ α 0 and χ α p for the sake of consistency. AC, in the attempt to preserve the operatorial relations they stem from, can lead to values of the unknowns slightly off their physical bounds in the spirit of using them as mere parameters to achieve the ultimate task of satisfying the operatorial algebra at the level of averages.
In figure 1, we report the behaviour of the chemical potential µ and of the double occupancy D as functions of the filling n for U = 1, 2 and 4 and T = 1/6. The COM 3-pole solution [COM(3p)] is clearly in very good agreement for all values of U reported in the whole range of filling n with the 12 × 12-site qMC [94] and 2-site DCA [95]  1The product of two symmetric matrices has real eigenvalues if one of the two is semi-positive, that is, it has positive or null eigenvalues.

Theory
In this case, we choose as basic field just the composite doublet field operator ψ (i) = ξ † (i) , η † (i) . By considering the two-time thermodynamic GF [96][97][98], we define the retarded GF G(i, j) = R[ψ(i)ψ † ( j)] = θ(t i − t j ) {ψ(i), ψ † ( j)} that satisfies the following equation where the derivative operator Λ(∂ i ) = i ∂ ∂t i − ε(−i∇ i ) and the propagator G 0 (i, j) is defined by the equation Λ(∂ i )G 0 (i, j) = iδ(t i − t j )I(i, j). By introducing the Fourier transform, equation (4.1) can be formally solved as where the self-energy Σ(k, ω) has the expression Σ(k, ω) = B irr (k, ω) is the Dyson equation for composite fields and represents the starting point for a perturbative calculation in terms of the propagator G 0 (k, ω) = 1 ω−ε(k) I(k), which corresponds to the 2-pole approximation and can be easily obtained by the calculations in the previous sections.
The calculation of the self-energy Σ(k, ω) requires the calculation of the higher-order propagator B(k, ω). We shall compute this quantity by using the non-crossing approximation (NCA). By neglecting the pair term c(i)c †α (i)c(i), the source J(i) can be written as J(i, t) = j a(i, j, t)ψ(j, t) where .

(4.4)
In order to calculate the retarded function F(i, j), first we use the spectral theorem to express where C(i − j, ω ) is the correlation function C(i, j) = σ µ δn µ (i)c α (i)c †α ( j)δn λ ( j)σ λ . Next, we use the NCA and approximate σ µ δn µ (i)c α (i)c †α ( j)δn λ ( j)σ λ ≈ δn µ (i)δn µ ( j) c α (i)c †α ( j) . By means of this decoupling and using again the spectral theorem, we finally have where G cc (k, ω) = 2 a,b=1 G ab (k, ω) is the retarded electronic GF and χ(k, ω) = µ F R[δn µ (i)δn µ ( j)] is the total charge and spin dynamical susceptibility. The result (4.6) shows that the self-energy calculation requires the knowledge of the bosonic propagator. Remarkably, up to this point, the system of equations for the GF and the anomalous self-energy is similar to the one derived in the two-particle self-consistent approach (TPSC) [99,100], the DMFT+Σ approach [101][102][103][104][105] and a Mori-like approach by Plakida and coworkers [41][42][43]. It would be the way to compute the dynamical spin and charge susceptibilities to be completely different because, instead of relying on a phenomenological model and neglecting the charge susceptibility like these approaches do, we will use a self-consistent two-pole approximation [51].
Finally, the electronic GF G(k, ω) is computed through the following self-consistency scheme: we first compute G 0 (k, ω) and χ µ (k, ω) in the two-pole approximation, then Σ(k, ω) and consequently G(k, ω). Finally, we check how much the fermionic parameters (µ, ∆, and p) changed and decide whether to stop or to continue by computing new χ µ (k, ω) and Σ(k, ω) after G(k, ω).

Results
We intend to characterize some electronic properties by computing the spectral function A(k, ω) = − 1 π Im[G cc (k, ω)] and the density of states per spin N(ω) = 1 . We also investigate the electronic self-energy Σ cc (k, ω), defined through the equation where 0 (k) = −µ − 4tα(k) is the non-interacting dispersion, and introduce the quantity r(k) = 0 (k) + Σ cc (k, ω = 0) that determines the Fermi surface locus in the momentum space, r(k) = 0, in a Fermi liquid. The actual Fermi surface (or its relic in a non-Fermi-liquid) is given by the relative maxima of A(k, ω = 0), which takes into account, on equal footing, both Σ cc (k, ω) and Σ cc (k, ω), and is directly related, within the sudden approximation and forgetting any selection rules, to the effective ARPES measurements.

Spectral function and dispersion
The electronic dispersion, or better its relic in a strongly correlated system, can be obtained through the relative maxima of A(k, ω). In figure 2, we show the latter, in scale of grays (red is for the above-scale 33701-7 In the regions where Σ (k, ω) is instead finite, A(k, ω) assumes very low values, very difficult to be detected by ARPES, which, accordingly, would report only the red areas in the picture. This is crucial to understand the experimental evidences on the Fermi surface in the underdoped regime and to explain ARPES findings together with those of quantum oscillations measurements.
The dispersion loses significance close to M where A(k, ω) loses weight as Σ (k, ω) increases: correspondingly, χ 3 (k, ω) is strongly peaked at M due to the strong antiferromagnetic correlations in the system. The bandwidth reduces from 4t to values of the order J = 4t 2 /U = 0.5t, as expected for the dispersion of few holes in a strong antiferromagnetic background. The characteristics of the dispersion are also compatible with this scenario, such as the sequence of minima and maxima, actually driven by the doubling of the Brillouin zone induced by the strong antiferromagnetic correlations, as well as the dynamical formation of a t diagonal hopping indicated by the pronounced warping of the dispersion along the X → Y direction. Moreover, remarkably, the dispersion at X(Y ) coming from both Γ and M is substantially flat: this feature is in very good agreement both with qMC computations (see [106] and references therein) and with ARPES experiments [107], which detect a similar behaviour for overdoped systems. It is worth noting that such flatness at X(Y ) is present for larger dopings too (n = 0.7, 0.78 and 0.85), as shown in figure 3 of reference [15], and that the extension of the plateau increases upon decreasing the doping. It is now clear that the dispersion warping along X → Y , displayed in figure 2, is also responsible for two maxima in N (ω): one related to the van-Hove singularities at X and Y and one to the dispersion maximum close to S. The entity of the dip between these two maxima just depends on the number of available well-defined (red) states in k present between these two values of ω. This will determine the formation of a more or less pronounced pseudogap in N (ω), as we will discuss after the analysis of the region close to M.
Definitely, the absence of spectral weight close to M and, in particular and surprisingly, at ω = 0 (i.e., on the Fermi surface, in contradiction with the Fermi-liquid scenario), is the most relevant result of this study: it will determine almost all relevant and anomalous features of the single-particle properties. Last, but not least, quite remarkable is also the presence of kinks in the dispersion in both the nodal (Γ → M) and the antinodal (X → Γ) directions, as highlighted by the dark green guidelines, in qualitative agreement with some ARPES experiments [108]. Remarkably, similar results for the singleparticle excitation spectrum were obtained in the self-consistent projection operator method [109,110], the operator projection method [111][112][113] and a Mori-like approach by Plakida and coworkers [42,43].

Spectral function and Fermi surface
Considering A(k, ω = 0), we attain the closest concept to Fermi surface for a strongly correlated system. In figure 3, we show A(k, ω = 0) as a function of k in a quarter of the Brillouin zone for U = 8 and (left) n = 0.85 and T = 0.01 and (right) n = 0.92 and T = 0.02. According to the interpretation of the ARPES measurements in the sudden approximation [107,108], the Fermi surface can be defined as the locus in k space of the relative maxima of A(k, ω = 0). Such a definition leads to a possible explanation of ARPES measurements, but also allows one to go beyond them, with their finite resolution and sensitivity, aiming at a reconciliation with other types of measurements, in particular quantum oscillations ones, which seemingly report results in disagreement, also up to dichotomy in some cases, with ARPES.
The evolution of the Fermi surface topology on increasing n is the consequence of two Lifshitz transitions [114], at n 0.82 and n 0.9. The first Lifshitz transition (n 0.82, not shown) is due to the crossing of the chemical potential by the van Hove singularity: the exact doping can be easily identified by looking at the chemical potential behaviour as a function of n (not shown). The second transition (n 0.9, not shown) is characterized by the formation of a hole pocket in the proximity of S. Changes in the Fermisurface topology were also found in other studies on the cuprates [37,42,115]. A very detailed analysis for the Lifshitz transitions in the t-J model was performed by Korshunov and Ovchinnikov [116]: they obtained the model parameters from ab initio calculations for cuprates [117] and found a large overdoped Fermi surface for n < 0.76, two concentric Fermi surfaces for 0.76 < n < 0.85 and a small hole Fermi surface around S for n > 0.85. Consequently, in reference [118], Ovchinnikov et al. characterized better these two quantum phase transitions, obtaining a logarithmic singularity for N (ω = 0) at n = 0.85 and a stepwise feature at n = 0.76.
Then, in figure 3, we report results for n = 0.85 and n = 0.92 as representatives of two non-trivial topologies of the Fermi surface. We can clearly distinguish two arcs that for n = 0.92 somehow join. Given the current sensitivities, only the arc with the larger intensities, among the two, may be visible to ARPES: at n = 0.92, the region close to S is the only one with an appreciable signal. The less-intense arc, reported in reference [42] too, is the relic of a shadow band, as clearly seen in figure 2, and thus never changes its curvature, in contrast to what happens to the other arc, subject to the crossing of the van Hove singularity (n 0.82, not shown).
Three additional ingredients can help us better understand the evolution with doping of the Fermi surface: (i) the n(k) = 0.5 locus (solid line), i.e., the Fermi surface if the system would be non-interacting; (ii) the r(k) = 0 locus (dashed line), i.e., the Fermi surface if the system would be a Fermi liquid or a state close to it conceptually; (iii) the values (grey lines and labels) of Σ cc (k, ω = 0). Through the combined analysis of these three ingredients with the positions and intensities of the relative maxima of A(k, ω = 0), we can understand better what these latter imply and classify the behaviour of the system on changing filling. At high dopings (not shown), the positions of the two arcs match exactly r(k) = 0 lines: this fact corroborates our definition of Fermi surface, making it versatile and valid beyond the Brillouin zone. From [15].

33701-9
Fermi liquid picture without contradicting this latter. Increasing the filling, we find the first topological transition from a Fermi surface closed around Γ (hole like in cuprates language) to a Fermi surface closed around M (electron like in cuprates language) at n 0.82 (not shown), where the van Hove singularity is at ω = 0 [see figure 3 (left-hand panel) for n = 0.85, a close doping]. Close to the antinodal points (X and Y ), a net discrepancy between the position of the relative maxima of A(k, ω = 0) and the line n(k) = 0.5 arises on decreasing doping. This feature does not only allow the topological transition, absent for the n(k) = 0.5 locus that reaches the anti-diagonal (X → Y ) at n = 1 satisfying the Luttinger theorem, but it also accounts for the broadening of the relative maxima of A(k, ω = 0) close to the anti-nodal points. The broadening is due to the small, but finite, value of Σ cc (k, ω = 0) in those regions, signaling the net enhancement of the correlation strength, and then the impossibility to consider the system in this regime as a conventional non-(or weakly-) interacting one within a Fermi-liquid scenario or its ordinary extensions for ordered phases. The emergence of such features only in well defined regions in k space is of great interest and goes beyond the problem under analysis (cuprates).
As the most important result, at n = 0.92, the relative maxima of A(k, ω = 0) deviate also from the r(k) = 0 line, at least partially, leading to a completely new scenario: r(k) = 0 defines a pocket, while the peaks of A(k, ω = 0) feature the very same pocket together with quite well-defined wings closing, with one half of the pocket, a kind of relic of a large Fermi surface. This is the second and most surprising topological transition for the Fermi surface; in fact, the two arcs, clearly visible for all other dopings, join and do not close just a pocket, as expected from the conventional theory for an antiferromagnet: they develop a fully independent branch. The actual Fermi surface is neither a pocket nor a large Fermi surface. This very unexpected result can be connected to the dichotomy between the experiments (e.g., ARPES) pointing to a small and the ones (e.g., quantum oscillations) pointing to a large Fermi surface.
The pocket too is definitely unconventional. In fact, there are two distinct halves of the pocket: one with very high intensity pinned at S (the only possibly detectable by ARPES) and one with very low intensity (detectable only by some quantum oscillations experiments). This is our interpretation for the Fermi arcs, seen in many ARPES experiments [108,119] and unaccountable for any ordinary theory based on the Fermi liquid picture, though modified by including an incipient spin or charge ordering. Clearly, looking only at the Fermi arc (as necessarily in ARPES), the Fermi surface looks ill defined, not enclosing a definite region of k space, but having access also to the other half of the pocket, such problem is greatly alleviated. In our understanding, the antiferromagnetic fluctuations are so strong to destroy the quasi-particle coherence in that region of k space, as similarly reported in the DMFT+Σ approach [101][102][103][104][105] and in a Mori-like approach by Plakida and collaborators [42,43]. The Fermi arc is pinned at S: this pinning of the center of mass of the ARPES-visible Fermi arc has been obtained also in ARPES experiments [120].

Density of states and pseudogap
The other main feature in underdoped cuprates is a substantial depletion in the electronic N (ω): the pseudogap. In figure 4 (left-hand panel), we show N(ω) for U = 8 and four couples of values of filling and temperature: n = 0.7 and T = 0.01, n = 0.78 and T = 0.01, n = 0.85 and T = 0.01, and n = 0.92 and T = 0.02, in the frequency region close to the chemical potential. As a reference, we also show, in figure 4 (right-hand panel), A(k, ω ∼ 0) at k = S = (π/2, π/2), k = S which is where the phantom half of the pocket touches the diagonal Γ → M (i.e., where the dispersion cuts the diagonal Γ → M closer to M), and k = X for U = 8, n = 0.92 and T = 0.02. As apparent in figure 4 (left-hand panel), N (ω) has two maxima separated by a dip, which plays the role of pseudogap. Its presence is due to the dispersion warping along the X → Y direction (see figure 2), which induces the two maxima [one due to the van-Hove singularity at X and one due to the dispersion maximum close to S -see figure 4 (right-hand panel)] and the loss of states, in this frequency window, in the region in k close to M. In order to realize the weight loss due to the finite value of Σ cc (k, ω = 0) in the region in k close to M, we can look at the impressive difference between the values of A(S, ω = 0) and A(S, ω = 0) [figure 4 (right-hand panel)]. Increasing the filling, there is a net spectral weight transfer between the two maxima; in particular, from the dispersion top close to S to the antinodal point X, where the van Hove singularity resides. At the lowest doping (n = 0.92), a well formed pseudogap can be seen for ω < 0 and will clearly influence the properties of the system. For n = 0.92, we do not find any divergence of Σ cc (k, ω = 0), in contrast to what is stated in reference [121] where this feature is presented as the definite reason for the pseudogap formation. In our study, the pseudogap is just the result of the weight transfer from the single-particle density of states to the two-particle one, related to the (antiferro)magnetic excitations occurring in the system on decreasing doping at low T. A similar doping behaviour of the pseudogap has been found by the DMFT+Σ approach [101][102][103][104][105], a Mori-like approach by Plakida and collaborators [42,43] and the cluster perturbation theory [100,122,123].

Conclusions
In this review paper, written in honor of the 80 th anniversary of Professor Ihor Stasyuk, we have illustrated the composite operator method by means of a prototypical example: the 2D Hubbard model and its application to the description of high-T c cuprate superconductors. The method has been reported first for a general case in order to appreciate its philosophy and its capability to be applied to any strongly correlated system in a controlled fashion.