# Composite Operator Method Analysis of the Underdoped Cuprates Puzzle

Authors: Adolfo Avella
Journal: Advances in Condensed Matter Physics
 (2014)

---

## Abstract

The microscopical analysis of the unconventional and puzzling physics of the underdoped cuprates, as carried out lately by means of the composite operator method (COM) applied to the 2D Hubbard model, is reviewed and systematized. The 2D Hubbard model has been adopted as it has been considered the minimal model capable of describing the most peculiar features of cuprates held responsible for their anomalous behavior. COM is designed to endorse, since its foundation, the systematic emergence in any SCS of new elementary excitations described by composite operators obeying noncanonical algebras. In this case (underdoped cuprates—2D Hubbard model), the residual interactions—beyond a 2-pole approximation—between the new elementary electronic excitations, dictated by the strong local Coulomb repulsion and well described by the two Hubbard composite operators, have been treated within the noncrossing approximation. Given this recipe and exploiting the few unknowns to enforce the Pauli principle content in the solution, it is possible to qualitatively describe some of the anomalous features of high-Tc cuprate superconductors such as large versus small Fermi surface dichotomy, Fermi surface deconstruction (appearance of Fermi arcs), nodal versus antinodal physics, pseudogap(s), and kinks in the electronic dispersion. The resulting scenario envisages a smooth crossover between an ordinary weakly interacting metal sustaining weak, short-range antiferromagnetic correlations in the overdoped regime to an unconventional poor metal characterized by very strong, long-but-finite-range antiferromagnetic correlations leading to momentum-selective non-Fermi liquid features as well as to the opening of a pseudogap and to the striking differences between the nodal and the antinodal dynamics in the underdoped regime.


---

## Body

## 1. Introduction

### 1.1. Composite Fields

One of the most intriguing challenges in modern condensed matter physics is the theoretical description of the anomalous behaviors experimentally observed in many novel materials. By anomalous behaviors we mean those not predicted by standard many-body theory, that is, behaviors in contradiction with the Fermi-liquid framework and diagrammatic expansions. The most relevant characteristic of such novel materials is the presence of so strong correlations among the electrons that classical schemes based on the band picture and the perturbation theory are definitely inapplicable. Accordingly, it is necessary to move from asingle-electron physics to amany-electron physics, where the dominant contributions come from the strong interactions among the electrons: usual schemes are simply inadequate and new concepts must be introduced.Theclassical techniques are based on the hypothesis that the interactions among the electrons are weak enough, or sufficiently well screened, to be properly taken into account within the framework of perturbative/diagrammatic methods. However, as many and many experimental and theoretical studies of highly correlated systems have shown, with more and more convincing evidence, all these methods are no more viable. The main concept that breaks down is the existence of the electrons as particles or quasiparticles with quite-well-defined properties. The presence of the interactions radically modifies the properties of the particles and, at a macroscopic level, what are observed are new particles (actually they are the only observable ones) with new peculiar properties entirely determined by the dynamics and by the boundary conditions (i.e., the phase under study, the external fields,
…). These new objects appear as the final result of the modifications imposed by the interactions on the original particles and contain, by the very beginning, the effects of correlations.On the basis of this evidence, one is induced to move the attention from the original fields to the new fields generated by the interactions. The operators describing these excitations, once they have been identified, can be written in terms of the original ones and are known as composite operators. The necessity of developing a formulation to treat composite operators as fundamental objects of the many-body problem in condensed matter physics has been deeply understood and systematically noticed since quite long time. Recent years have seen remarkable achievements in the development of a modern many-body theory in solid-state physics in the form of an assortment of techniques that may be termed composite particle methods. The foundations of these types of techniques may be traced back to the work of Bogoliubov [1] and later to that of Dancoff [2]. The work of Zwanzig [3], Mori [4–9], and Umezawa [10] definitely deserves to be mentioned too. Closely related to this work is that of Hubbard [11–13], Rowe [14], Roth [15], and Tserkovnikov [16, 17]. The slave boson method [18–20], the spectral density approach [21, 22], the diagram technique for Hubbard operators [23], the cumulant expansion based diagram technique [24], the generalized tight-binding method [25–27], self-consistent projection operator method [28, 29], operator projection method [30–32], and the composite operator method (COM) [33–35] are along the same lines. This large class of theories is very promising as it is based on the firm conviction that strong interactions call for an analysis in terms of new elementary fields embedding the greatest possible part of the correlations so permitting to overcome the problem of finding an appropriate expansion parameter. However, one price must be paid. In general, composite fields are neither Fermi nor Bose operators, since they do not satisfy canonical (anti)commutation relations, and their properties must be determined self-consistently. They can only be recognized as fermionic or bosonic operators according to the number and type of the constituting original particles. Accordingly, new techniques have to be developed in order to deal with such composite fields and to design diagrammatic schemes where the building blocks are the propagators of such composite fields: standard diagrammatic expansions and the Wick’s theorem are no more valid. The formulation of the Green’s function method itself must be revisited and new frameworks of calculations have to be devised.Following these ideas, we have been developing a systematic approach, the composite operator method (COM) [33–35], to study highly correlated systems. The formalism is based on two main ideas: (i) use of propagators of relevant composite operators as building blocks for any subsequent approximate calculations and (ii) use of algebra constraints to fix the representation of the relevant propagators in order to properly preserve algebraic and symmetry properties; these constraints will also determine the unknown parameters appearing in the formulation due to the noncanonical algebra satisfied by the composite operators. In the last fifteen years, COM has been applied to several models and materials: Hubbard [36–47], p-d [48–51], t-J [52], t-t
′-U [53–55], extended Hubbard (t-U-V) [56, 57], Kondo [58], Anderson [59, 60], two-orbital Hubbard [61–63], Ising [64, 65], J
1
-
J
2 [66–71], Hubbard-Kondo [72], Cuprates [73–78], and so forth.


### 1.2. Underdoped Cuprates

Cuprate superconductors [79] display a full range of anomalous features, mainly appearing in the underdoped region, in almost all experimentally measurable physical properties [80–84]. According to this, their microscopic description is still an open problem: non-Fermi-liquid response, quantum criticality, pseudogap formation, ill-defined Fermi surface, kinks in the electronic dispersion, and so forth still remain unexplained (or at least controversially debated) anomalous features [84–87]. In the last years, the attention of the community has been focusing on three main experimental facts [84]: the dramatic change in shape and nature of the Fermi surface between underdoped and overdoped regimes, the appearance of a psedudogap in the underdoped regime, and the striking differentiation between the physics at the nodes and at the antinodes in the pseudogap regime. The topological transition of the Fermi surface has been first detected by means of ARPES [82] and reflects the noteworthy differences between the quite-ordinary, large Fermi surface measured in the overdoped regime [88–90] and quite well described by LDA calculations [91] and quantum oscillations measurements [92], and the ill-defined Fermi arcs appearing in the underdoped regime [85, 93–98]. The enormous relevance of these experimental findings, not only for the microscopic comprehension of the high-T
c superconductivity phenomenology, but also for the drafting of a general microscopic theory for strongly correlated materials, called for many more measurements in order to explore all possible aspects of such extremely anomalous and peculiar behavior: plenty of quantum oscillations measurements in the underdoped regime [99–118], Hall effect measurements [100], Seebeck effect measurements [116], and heat capacity measurements [115]. The presence of a quite strong depletion in the electronic density of states, known as pseudogap [119], is well established thanks to ARPES [82], NMR [120], optical conductivity [121], and quantum oscillations [122] measurements. The microscopic origin of such a loss of single-particle electronic states is still unclear and the number of possible theoretical, as well as phenomenological, explanations has grown quite large in the last few years. As a matter of fact, this phenomenon affects any measurable properties and, accordingly, was the first to be detected in the underdoped regime granting to this latter the first evidences of its exceptionality with respect to the other regimes in the phase diagram. The plethora of theoretical scenarios present in the literature [85, 123–125], tentatively explaining few, some, or many of the anomalous features reported by the experiments on underdoped cuprates, can be coarsely divided into those not relying on any translational symmetry breaking [126–131] and those instead proposing that it should be some kind of charge and/or spin arrangements to be held responsible for the whole range of anomalous features. Among these latter theories, there are those focusing on the physics at the antinodal region and those focusing on the nodal region. We can account for proposals of (as regards the antinode): a collinear spin AF order [132], an AF quantum critical point [133, 134], and a 1D charge stripe order [132, 135] with the addition of a smectic phase [136]. Instead, at the node, we have a d-density wave [137], a more-or-less ordinary AF spin order [75–78, 138–142], and a nodal pocket from bilayer low-Q charge order, slowly fluctuacting [135, 143–145]. This latter proposal, which is among the newest on the table, relies on many experimental measurements: STM [146–150], neutron scattering [122], X-ray diffraction [151], NMR [152], RXS [153], and phonon softening [154–156]. A proposal regarding the emergence of a hidden Fermi liquid [157] is also worth mentioning.


### 1.3. 2D Hubbard Model: Approximation Methods

Since the very beginning [158], the two-dimensional Hubbard model [11] has been universally recognized as the minimal model capable of describing the Cu–O2 planes of cuprates superconductors. It certainly contains many of the key ingredients by construction: strong electronic correlations, competition between localization and itineracy, Mott physics, and low-energy spin excitations. Unfortunately, though being fundamental for benchmarking and fine tuning analytical theories, numerical approaches [159] cannot be of help to solve the puzzle of underdoped cuprates owing to their limited resolution in frequency and momentum. On the other hand, there are not so many analytical approaches capable of dealing with the quite complex aspects of underdoped cuprates phenomenology [9]. Among others, the two-particle self-consistent (TPSC) approach [86, 87, 160] has been the first completely microscopic approach to obtain results comparable with the experimental findings. Almost all other promising approaches available in the literature can be essentially divided into two classes. One class makes use of phenomenological expressions for the electronic self-energy and the electronic spin susceptibility [5–9, 161, 162]. The electronic self-energy is usually computed as the convolution of the electronic propagator and of the electronic spin susceptibility. Then, the electronic spin susceptibility is modeled phenomenologically parameterizing correlation length and damping as functions of doping and temperature according to the common belief that the electronic spin susceptibility should present a well-developed mode at M
=
(
π
,
π
) with a damping of Landau type. The DMFT
+
Σ approach [163–167] and the DGammaA approach [168, 169] should also be mentioned. All cluster-dynamical-mean-field-like theories (cluster-DMFT theories) [86, 87, 170, 171] (the cellular dynamical mean-field theory (C-DMFT) [172, 173], the dynamical cluster approximation (DCA) [174], and the cluster perturbation theory (CPT) [86, 87, 175, 176]) belong to the second class. The dynamical mean-field theory (DMFT) [177–181] cannot tackle the underdoped cuprates puzzle because its self-energy has the same identical value at each point on the Fermi surface without any possible differentiation between nodal and antinodal physics or visible and phantom portion of the Fermi surface. The cluster-DMFT theories instead can, in principle, deal with both coherent quasi-particles and marginal ones within the same Fermi surface. These theories usually self-consistently map the generic Hubbard problem to a few-site lattice Anderson problem and solve this latter by means of, mainly, numerical techniques. What really distinguishes one formulation from another, within this second class, is the procedure used to map the small cluster on the infinite lattice. Anyway, it is worth noticing that these approaches often rely on numerical methods in order to close their self-consistency cycles (with the abovementioned limitations in frequency and momentum resolutions and with the obvious difficulties in the physical interpretation of their results) and always face the emergence of a quite serious periodization problem since a cluster embedded in the lattice violates its periodicity. The Composite Operator Method (COM) [33–35] does not belong to any of these two classes of theoretical formulations and has the advantage to be completely microscopic, exclusively analytical, and fully self-consistent. COM recipe uses three main ingredients [33–35]:composite operators,algebra constraints, andresidual self-energy treatment. Composite operators are products of electronic operators and describe the new elementary excitations appearing in the system owing to strong correlations. According to the system under analysis [33–35], one has to choose a set of composite operators as operatorial basis and rewrite the electronic operators and the electronic Green’s function in terms of this basis. One should think of composite operators just as a more convenient starting point, with respect to electronic operators, for any mean-field-like approximation/perturbation scheme. Algebra constraints are relations among correlation functions dictated by the noncanonical operatorial algebra closed by the chosen operatorial basis [33–35]. Other ways to obtain algebra constraints rely on the symmetries enjoined by the Hamiltonian under study, the Ward-Takahashi identities, the hydrodynamics, and so forth, [33–35]. One should think of algebra constraints as a way to restrict the Fock space on which the chosen operatorial basis acts to the Fock space of physical electrons. Algebra constraints are used to compute unknown correlation functions appearing in the calculations. Interactions among the elements of the chosen operatorial basis are described by the residual self-energy, that is, the propagator of the residual term of the current after this latter has been projected on the chosen operatorial basis [33–35]. According to the physical properties under analysis and the range of temperatures, dopings, and interactions you want to explore, one has to choose an approximation to compute the residual self-energy. In the last years, we have been using the n-pole Approximation [36–43, 46, 48–57, 61–65, 72–74], the Asymptotic Field Approach [58–60], the NCA [44, 45, 75–78, 182], and the Two-Site Resolvent Approach [183, 184]. You should think of the residual self-energy as a measure in the frequency and momentum space of how much well defined are, as quasi-particles, your composite operators. It is really worth noticing that, although the description of some of the anomalous features of underdoped cuprates given by COM qualitatively coincides with those obtained by TPSC [86, 87] and by the two classes of formulations mentioned above, the results obtained by means of COM greatly differ from those obtained within the other methods as regards the evolution with doping of the dispersion and of the Fermi surface and, at the moment, no experimental result can tell which is the unique and distinctive choice nature made.


### 1.4. Outline

To study the underdoped cuprates modeled by the 2D Hubbard model (see Section2.1), we start from a basis of two composite operators (the two Hubbard operators) and formulate the Dyson equation (see Section 2.2) in terms of the 2-pole approximated Green’s function (see Section 2.3). According to this, the self-energy is the propagator of nonlocal composite operators describing the electronic field dressed by charge, spin, and pair fluctuations on the nearest-neighbor sites. Then, within the noncrossing approximation (NCA) [185], we obtain amicroscopic self-energy written in terms of the convolution of the electronic propagator and of the charge, spin, and pair susceptibilities (see Section 2.4). Finally, we close, fully analytically, the self-consistency cycle for the electronic propagator by computingmicroscopic susceptibilities within a 2-pole approximation (see Section 2.5). Our results (see Section 2.4) show that, within COM, the two-dimensional Hubbard model can describe some of the anomalous features experimentally observed in underdoped cuprates phenomenology. In particular, we show how Fermi arcs can develop out of a large Fermi surface (see Section 3.2), how pseudogap can show itself in the dispersion (see Section 3.1) and in the density of states (see Section 3.3), how non-Fermi-liquid features can become apparent in the momentum distribution function (see Section 3.4) and in the frequency and temperature dependences of the self-energy (see Section 3.5), how muchkinked the dispersion can get on varying doping (see Section 3.1), and why, or at least how, spin dynamics can be held responsible for all this (see Section 3.6). Finally (see Section 4), we summarize the current status of the scenario emerging by these theoretical findings and which are the perspectives.


## 1.1. Composite Fields

One of the most intriguing challenges in modern condensed matter physics is the theoretical description of the anomalous behaviors experimentally observed in many novel materials. By anomalous behaviors we mean those not predicted by standard many-body theory, that is, behaviors in contradiction with the Fermi-liquid framework and diagrammatic expansions. The most relevant characteristic of such novel materials is the presence of so strong correlations among the electrons that classical schemes based on the band picture and the perturbation theory are definitely inapplicable. Accordingly, it is necessary to move from asingle-electron physics to amany-electron physics, where the dominant contributions come from the strong interactions among the electrons: usual schemes are simply inadequate and new concepts must be introduced.Theclassical techniques are based on the hypothesis that the interactions among the electrons are weak enough, or sufficiently well screened, to be properly taken into account within the framework of perturbative/diagrammatic methods. However, as many and many experimental and theoretical studies of highly correlated systems have shown, with more and more convincing evidence, all these methods are no more viable. The main concept that breaks down is the existence of the electrons as particles or quasiparticles with quite-well-defined properties. The presence of the interactions radically modifies the properties of the particles and, at a macroscopic level, what are observed are new particles (actually they are the only observable ones) with new peculiar properties entirely determined by the dynamics and by the boundary conditions (i.e., the phase under study, the external fields,
…). These new objects appear as the final result of the modifications imposed by the interactions on the original particles and contain, by the very beginning, the effects of correlations.On the basis of this evidence, one is induced to move the attention from the original fields to the new fields generated by the interactions. The operators describing these excitations, once they have been identified, can be written in terms of the original ones and are known as composite operators. The necessity of developing a formulation to treat composite operators as fundamental objects of the many-body problem in condensed matter physics has been deeply understood and systematically noticed since quite long time. Recent years have seen remarkable achievements in the development of a modern many-body theory in solid-state physics in the form of an assortment of techniques that may be termed composite particle methods. The foundations of these types of techniques may be traced back to the work of Bogoliubov [1] and later to that of Dancoff [2]. The work of Zwanzig [3], Mori [4–9], and Umezawa [10] definitely deserves to be mentioned too. Closely related to this work is that of Hubbard [11–13], Rowe [14], Roth [15], and Tserkovnikov [16, 17]. The slave boson method [18–20], the spectral density approach [21, 22], the diagram technique for Hubbard operators [23], the cumulant expansion based diagram technique [24], the generalized tight-binding method [25–27], self-consistent projection operator method [28, 29], operator projection method [30–32], and the composite operator method (COM) [33–35] are along the same lines. This large class of theories is very promising as it is based on the firm conviction that strong interactions call for an analysis in terms of new elementary fields embedding the greatest possible part of the correlations so permitting to overcome the problem of finding an appropriate expansion parameter. However, one price must be paid. In general, composite fields are neither Fermi nor Bose operators, since they do not satisfy canonical (anti)commutation relations, and their properties must be determined self-consistently. They can only be recognized as fermionic or bosonic operators according to the number and type of the constituting original particles. Accordingly, new techniques have to be developed in order to deal with such composite fields and to design diagrammatic schemes where the building blocks are the propagators of such composite fields: standard diagrammatic expansions and the Wick’s theorem are no more valid. The formulation of the Green’s function method itself must be revisited and new frameworks of calculations have to be devised.Following these ideas, we have been developing a systematic approach, the composite operator method (COM) [33–35], to study highly correlated systems. The formalism is based on two main ideas: (i) use of propagators of relevant composite operators as building blocks for any subsequent approximate calculations and (ii) use of algebra constraints to fix the representation of the relevant propagators in order to properly preserve algebraic and symmetry properties; these constraints will also determine the unknown parameters appearing in the formulation due to the noncanonical algebra satisfied by the composite operators. In the last fifteen years, COM has been applied to several models and materials: Hubbard [36–47], p-d [48–51], t-J [52], t-t
′-U [53–55], extended Hubbard (t-U-V) [56, 57], Kondo [58], Anderson [59, 60], two-orbital Hubbard [61–63], Ising [64, 65], J
1
-
J
2 [66–71], Hubbard-Kondo [72], Cuprates [73–78], and so forth.


## 1.2. Underdoped Cuprates

Cuprate superconductors [79] display a full range of anomalous features, mainly appearing in the underdoped region, in almost all experimentally measurable physical properties [80–84]. According to this, their microscopic description is still an open problem: non-Fermi-liquid response, quantum criticality, pseudogap formation, ill-defined Fermi surface, kinks in the electronic dispersion, and so forth still remain unexplained (or at least controversially debated) anomalous features [84–87]. In the last years, the attention of the community has been focusing on three main experimental facts [84]: the dramatic change in shape and nature of the Fermi surface between underdoped and overdoped regimes, the appearance of a psedudogap in the underdoped regime, and the striking differentiation between the physics at the nodes and at the antinodes in the pseudogap regime. The topological transition of the Fermi surface has been first detected by means of ARPES [82] and reflects the noteworthy differences between the quite-ordinary, large Fermi surface measured in the overdoped regime [88–90] and quite well described by LDA calculations [91] and quantum oscillations measurements [92], and the ill-defined Fermi arcs appearing in the underdoped regime [85, 93–98]. The enormous relevance of these experimental findings, not only for the microscopic comprehension of the high-T
c superconductivity phenomenology, but also for the drafting of a general microscopic theory for strongly correlated materials, called for many more measurements in order to explore all possible aspects of such extremely anomalous and peculiar behavior: plenty of quantum oscillations measurements in the underdoped regime [99–118], Hall effect measurements [100], Seebeck effect measurements [116], and heat capacity measurements [115]. The presence of a quite strong depletion in the electronic density of states, known as pseudogap [119], is well established thanks to ARPES [82], NMR [120], optical conductivity [121], and quantum oscillations [122] measurements. The microscopic origin of such a loss of single-particle electronic states is still unclear and the number of possible theoretical, as well as phenomenological, explanations has grown quite large in the last few years. As a matter of fact, this phenomenon affects any measurable properties and, accordingly, was the first to be detected in the underdoped regime granting to this latter the first evidences of its exceptionality with respect to the other regimes in the phase diagram. The plethora of theoretical scenarios present in the literature [85, 123–125], tentatively explaining few, some, or many of the anomalous features reported by the experiments on underdoped cuprates, can be coarsely divided into those not relying on any translational symmetry breaking [126–131] and those instead proposing that it should be some kind of charge and/or spin arrangements to be held responsible for the whole range of anomalous features. Among these latter theories, there are those focusing on the physics at the antinodal region and those focusing on the nodal region. We can account for proposals of (as regards the antinode): a collinear spin AF order [132], an AF quantum critical point [133, 134], and a 1D charge stripe order [132, 135] with the addition of a smectic phase [136]. Instead, at the node, we have a d-density wave [137], a more-or-less ordinary AF spin order [75–78, 138–142], and a nodal pocket from bilayer low-Q charge order, slowly fluctuacting [135, 143–145]. This latter proposal, which is among the newest on the table, relies on many experimental measurements: STM [146–150], neutron scattering [122], X-ray diffraction [151], NMR [152], RXS [153], and phonon softening [154–156]. A proposal regarding the emergence of a hidden Fermi liquid [157] is also worth mentioning.


## 1.3. 2D Hubbard Model: Approximation Methods

Since the very beginning [158], the two-dimensional Hubbard model [11] has been universally recognized as the minimal model capable of describing the Cu–O2 planes of cuprates superconductors. It certainly contains many of the key ingredients by construction: strong electronic correlations, competition between localization and itineracy, Mott physics, and low-energy spin excitations. Unfortunately, though being fundamental for benchmarking and fine tuning analytical theories, numerical approaches [159] cannot be of help to solve the puzzle of underdoped cuprates owing to their limited resolution in frequency and momentum. On the other hand, there are not so many analytical approaches capable of dealing with the quite complex aspects of underdoped cuprates phenomenology [9]. Among others, the two-particle self-consistent (TPSC) approach [86, 87, 160] has been the first completely microscopic approach to obtain results comparable with the experimental findings. Almost all other promising approaches available in the literature can be essentially divided into two classes. One class makes use of phenomenological expressions for the electronic self-energy and the electronic spin susceptibility [5–9, 161, 162]. The electronic self-energy is usually computed as the convolution of the electronic propagator and of the electronic spin susceptibility. Then, the electronic spin susceptibility is modeled phenomenologically parameterizing correlation length and damping as functions of doping and temperature according to the common belief that the electronic spin susceptibility should present a well-developed mode at M
=
(
π
,
π
) with a damping of Landau type. The DMFT
+
Σ approach [163–167] and the DGammaA approach [168, 169] should also be mentioned. All cluster-dynamical-mean-field-like theories (cluster-DMFT theories) [86, 87, 170, 171] (the cellular dynamical mean-field theory (C-DMFT) [172, 173], the dynamical cluster approximation (DCA) [174], and the cluster perturbation theory (CPT) [86, 87, 175, 176]) belong to the second class. The dynamical mean-field theory (DMFT) [177–181] cannot tackle the underdoped cuprates puzzle because its self-energy has the same identical value at each point on the Fermi surface without any possible differentiation between nodal and antinodal physics or visible and phantom portion of the Fermi surface. The cluster-DMFT theories instead can, in principle, deal with both coherent quasi-particles and marginal ones within the same Fermi surface. These theories usually self-consistently map the generic Hubbard problem to a few-site lattice Anderson problem and solve this latter by means of, mainly, numerical techniques. What really distinguishes one formulation from another, within this second class, is the procedure used to map the small cluster on the infinite lattice. Anyway, it is worth noticing that these approaches often rely on numerical methods in order to close their self-consistency cycles (with the abovementioned limitations in frequency and momentum resolutions and with the obvious difficulties in the physical interpretation of their results) and always face the emergence of a quite serious periodization problem since a cluster embedded in the lattice violates its periodicity. The Composite Operator Method (COM) [33–35] does not belong to any of these two classes of theoretical formulations and has the advantage to be completely microscopic, exclusively analytical, and fully self-consistent. COM recipe uses three main ingredients [33–35]:composite operators,algebra constraints, andresidual self-energy treatment. Composite operators are products of electronic operators and describe the new elementary excitations appearing in the system owing to strong correlations. According to the system under analysis [33–35], one has to choose a set of composite operators as operatorial basis and rewrite the electronic operators and the electronic Green’s function in terms of this basis. One should think of composite operators just as a more convenient starting point, with respect to electronic operators, for any mean-field-like approximation/perturbation scheme. Algebra constraints are relations among correlation functions dictated by the noncanonical operatorial algebra closed by the chosen operatorial basis [33–35]. Other ways to obtain algebra constraints rely on the symmetries enjoined by the Hamiltonian under study, the Ward-Takahashi identities, the hydrodynamics, and so forth, [33–35]. One should think of algebra constraints as a way to restrict the Fock space on which the chosen operatorial basis acts to the Fock space of physical electrons. Algebra constraints are used to compute unknown correlation functions appearing in the calculations. Interactions among the elements of the chosen operatorial basis are described by the residual self-energy, that is, the propagator of the residual term of the current after this latter has been projected on the chosen operatorial basis [33–35]. According to the physical properties under analysis and the range of temperatures, dopings, and interactions you want to explore, one has to choose an approximation to compute the residual self-energy. In the last years, we have been using the n-pole Approximation [36–43, 46, 48–57, 61–65, 72–74], the Asymptotic Field Approach [58–60], the NCA [44, 45, 75–78, 182], and the Two-Site Resolvent Approach [183, 184]. You should think of the residual self-energy as a measure in the frequency and momentum space of how much well defined are, as quasi-particles, your composite operators. It is really worth noticing that, although the description of some of the anomalous features of underdoped cuprates given by COM qualitatively coincides with those obtained by TPSC [86, 87] and by the two classes of formulations mentioned above, the results obtained by means of COM greatly differ from those obtained within the other methods as regards the evolution with doping of the dispersion and of the Fermi surface and, at the moment, no experimental result can tell which is the unique and distinctive choice nature made.


## 1.4. Outline

To study the underdoped cuprates modeled by the 2D Hubbard model (see Section2.1), we start from a basis of two composite operators (the two Hubbard operators) and formulate the Dyson equation (see Section 2.2) in terms of the 2-pole approximated Green’s function (see Section 2.3). According to this, the self-energy is the propagator of nonlocal composite operators describing the electronic field dressed by charge, spin, and pair fluctuations on the nearest-neighbor sites. Then, within the noncrossing approximation (NCA) [185], we obtain amicroscopic self-energy written in terms of the convolution of the electronic propagator and of the charge, spin, and pair susceptibilities (see Section 2.4). Finally, we close, fully analytically, the self-consistency cycle for the electronic propagator by computingmicroscopic susceptibilities within a 2-pole approximation (see Section 2.5). Our results (see Section 2.4) show that, within COM, the two-dimensional Hubbard model can describe some of the anomalous features experimentally observed in underdoped cuprates phenomenology. In particular, we show how Fermi arcs can develop out of a large Fermi surface (see Section 3.2), how pseudogap can show itself in the dispersion (see Section 3.1) and in the density of states (see Section 3.3), how non-Fermi-liquid features can become apparent in the momentum distribution function (see Section 3.4) and in the frequency and temperature dependences of the self-energy (see Section 3.5), how muchkinked the dispersion can get on varying doping (see Section 3.1), and why, or at least how, spin dynamics can be held responsible for all this (see Section 3.6). Finally (see Section 4), we summarize the current status of the scenario emerging by these theoretical findings and which are the perspectives.


## 2. Framework

### 2.1. Hamiltonian

The Hamiltonian of the two-dimensional Hubbard model reads as(1)
H
=
∑
i
j
‍
-
μ
δ
i
j
-
4
t
α
i
j
c
†
i
c
j
+
U
∑
i
‍
n
↑
i
n
↓
i
,
where
(2)
c
i
=
c
↑
i
c
↓
i
is the electron field operator in spinorial notation and Heisenberg picture (i
=
(
i
,
t
i
)), i
=
R
i is a vector of the Bravais lattice, n
σ
(
i
)
=
c
σ
†
(
i
)
c
σ
(
i
) is the particle density operator for spin σ, n
(
i
)
=
∑
σ
n
σ
(
i
) is the total particle density operator, μ is the chemical potential, t is the hopping integral and the energy unit, U is the Coulomb on-site repulsion, and α
i
j is the projector on the nearest-neighbor sites:
(3)
α
i
j
=
1
N
∑
k
‍
e
i
k
·
R
i
-
R
j
α
k
,
α
k
=
1
2
cos
⁡
k
x
a
+
cos
⁡
k
y
a
,
where k runs over the first Brillouin zone, N is the number of sites, and a is the lattice constant.


### 2.2. Green’s Functions and Dyson Equation

Following COM prescriptions [33–35], we chose a basic field; in particular, we select the composite doublet field operator:
(4)
ψ
i
=
ξ
i
η
i
,
where η
(
i
)
=
n
(
i
)
c
(
i
) and ξ
(
i
)
=
c
(
i
)
-
η
(
i
) are the Hubbard operators describing the main subbands. This choice is guided by the hierarchy of the equations of motion and by the fact that ξ
(
i
) and η
(
i
) are eigenoperators of the interacting term in the Hamiltonian (1). The field ψ
(
i
) satisfies the Heisenberg equation:
(5)
i
∂
∂
t
ψ
i
=
J
i
=
-
μ
ξ
i
-
4
t
c
α
i
-
4
t
π
i
U
-
μ
η
i
+
4
t
π
i
,
where the higher-order composite field π
(
i
) is defined by
(6)
π
i
=
1
2
σ
μ
n
μ
i
c
α
i
+
c
i
c
†
α
i
c
i
with the following notation: n
μ
(
i
)
=
c
†
(
i
)
σ
μ
c
(
i
) is the particle (μ
=
0) and spin (μ
=
1
,
2
,
3) density operator and σ
μ
=
1
,
σ
→, σ
μ
=
-
1
,
σ
→, and σ
k  k
=
1
,
2
,
3 are the Pauli matrices. Hereafter, for any operator Φ
(
i
), we use the notation Φ
α
(
i
,
t
)
=
∑
j
‍
α
i
j
Φ
(
j
,
t
).It is always possible to decompose the sourceJ
(
i
) under the form:
(7)
J
i
=
ɛ
-
i
∇
ψ
i
+
δ
J
i
,
where the linear term represents the projection of the source on the basis ψ
(
i
) and is calculated by means of the equation:
(8)
δ
J
i
,
t
,
ψ
†
j
,
t
=
0
,
where ⋯ stands for the thermal average taken in the grand-canonical ensemble.This constraint assures that the residual currentδ
J
(
i
) contains all and only the physics orthogonal to the chosen basis ψ
(
i
). The action of the derivative operator ɛ
(
-
i
∇
) on ψ
(
i
) is defined in momentum space
(9)
ɛ
-
i
∇
ψ
i
=
ɛ
-
i
∇
1
N
∑
k
‍
e
i
k
·
R
i
ψ
k
,
t
=
1
N
∑
k
‍
e
i
k
·
R
i
ɛ
k
ψ
k
,
t
,
where ɛ
(
k
) is named energy matrix.The constraint (8) gives
(10)
m
k
=
ɛ
k
I
k
after defining the normalization matrix
(11)
I
i
,
j
=
ψ
i
,
t
,
ψ
†
j
,
t
=
1
N
∑
k
‍
e
i
k
·
R
i
-
R
j
I
k
,
and the m-matrix
(12)
m
i
,
j
=
J
i
,
t
,
ψ
†
j
,
t
=
1
N
∑
k
‍
e
i
k
·
R
i
-
R
j
m
k
.Since the components ofψ
(
i
) contain composite operators, the normalization matrix I
(
k
) is not the identity matrix and defines the spectral content of the excitations. In fact, the composite operator method has the advantage of describing crossover phenomena as the phenomena in which the weight of some operator is shifted to another one.By considering the two-time thermodynamic Green’s functions [186–188], let us define the retarded function
(13)
G
i
,
j
=
R
ψ
i
ψ
†
j
=
θ
t
i
-
t
j
ψ
i
,
ψ
†
j
.By means of the Heisenberg equation (5) and using the decomposition (7), the Green’s function G
(
i
,
j
) satisfies the equation
(14)
Λ
∂
i
G
i
,
j
Λ
†
∂
←
j
=
Λ
∂
i
G
0
i
,
j
Λ
†
∂
←
j
+
R
δ
J
i
δ
J
†
j
,
where the derivative operator Λ
(
∂
i
) is defined as
(15)
Λ
∂
i
=
i
∂
∂
t
i
-
ɛ
-
i
∇
i
,
and the propagator G
0
(
i
,
j
) is defined by the equation
(16)
Λ
∂
i
G
0
i
,
j
=
i
δ
t
i
-
t
j
I
i
,
j
.By introducing the Fourier transform(17)
G
i
,
j
=
1
N
∑
k
‍
i
2
π
∫
‍
d
ω
e
i
k
·
R
i
-
R
j
-
i
ω
t
i
-
t
j
G
k
,
ω
equation (14) in momentum space can be written as
(18)
G
k
,
ω
=
G
0
k
,
ω
+
G
0
k
,
ω
I
-
1
k
Σ
k
,
ω
G
k
,
ω
,
and can be formally solved as
(19)
G
k
,
ω
=
1
ω
-
ɛ
k
-
Σ
k
,
ω
I
k
,
where the self-energy Σ
(
k
,
ω
) has the expression
(20)
Σ
k
,
ω
=
B
irr
k
,
ω
I
-
1
k
with
(21)
B
k
,
ω
=
F
R
δ
J
i
δ
J
†
j
.The notationF denotes the Fourier transform and the subscript irr indicates that the irreducible part of the propagator B
(
k
,
ω
) is taken. Equation (18) is nothing else than the Dyson equation for composite fields and represents the starting point for a perturbative calculation in terms of the propagator G
0
(
k
,
ω
). This quantity will be calculated in the next section. Then, the attention will be given to the calculation of the self-energy Σ
(
k
,
ω
). It should be noted that the computation of the two quantities G
0
(
k
,
ω
) and Σ
(
k
,
ω
) is intimately related. The total weight of the self-energy corrections is bounded by the weight of the residual source operator δ
J
(
i
). According to this, it can be made smaller and smaller by increasing the components of the basis ψ
(
i
) (e.g., by including higher-order composite operators appearing in δ
J
(
i
)). The result of such a procedure will be the inclusion in the energy matrix of part of the self-energy as an expansion in terms of coupling constants multiplied by the weights of the newly included basis operators. In general, the enlargement of the basis leads to a new self-energy with a smaller total weight. However, it is necessary pointing out that this process can be quite cumbersome and the inclusion of fully momentum and frequency dependent self-energy corrections can be necessary to effectively take into account low-energy and virtual processes. According to this, one can choose a reasonable number of components for the basic set and then use another approximation method to evaluate the residual dynamical corrections.


### 2.3. Two-Pole Approximation

According to (16), the free propagator G
0
(
k
,
ω
) is determined by the following expression:
(22)
G
0
k
,
ω
=
1
ω
-
ɛ
k
I
k
.For a paramagnetic state, straightforward calculations give the following expressions for the normalizationI
(
k
) and energy ɛ
(
k
) matrices:
(23)
I
k
=
1
-
n
2
0
0
n
2
=
I
11
0
0
I
22
,
ɛ
11
k
=
-
μ
-
4
t
I
11
-
1
Δ
+
1
-
n
+
p
α
k
,
ɛ
12
k
=
4
t
I
22
-
1
Δ
+
p
-
I
22
α
k
,
ɛ
21
k
=
4
t
I
11
-
1
Δ
+
p
-
I
22
α
k
,
ɛ
22
k
=
U
-
μ
-
4
t
I
22
-
1
Δ
+
p
α
k
,
where n
=
n
i is the filling and
(24)
Δ
=
ξ
α
i
ξ
†
i
-
η
α
i
η
†
i
,
p
=
1
4
n
μ
α
i
n
μ
i
-
c
↑
i
c
↓
i
α
c
↓
†
i
c
↑
†
i
.
Then, (22) can be written in spectral form as
(25)
G
0
k
,
ω
=
∑
n
=
1
2
‍
σ
n
k
ω
-
E
n
k
+
i
δ
.The energy spectraE
n
k and the spectral functions σ
(
n
)
(
k
) are given by
(26)
E
1
k
=
R
k
+
Q
k
,
E
2
k
=
R
k
-
Q
k
,
σ
11
1
k
=
I
11
2
1
+
g
k
2
Q
k
,
σ
11
2
k
=
I
11
2
1
-
g
k
2
Q
k
,
σ
12
1
k
=
m
12
k
2
Q
k
,
σ
12
2
k
=
-
m
12
k
2
Q
k
,
σ
22
1
k
=
I
22
2
1
-
g
k
2
Q
k
,
σ
22
2
k
=
I
22
2
1
+
g
k
2
Q
k
,
where
(27)
R
k
=
-
μ
-
4
t
α
k
+
1
2
U
-
ɛ
12
k
2
I
11
,
Q
k
=
1
2
g
2
k
+
4
ɛ
12
2
k
I
22
I
11
,
g
k
=
-
U
+
1
-
n
I
11
ɛ
12
k
.The energy matrixɛ
(
k
) contains three parameters: μ, the chemical potential, Δ, the difference between upper and lower intrasubband contributions to kinetic energy, and p, a combination of the nearest-neighbor charge-charge, spin-spin, and pair-pair correlation functions. These parameters will be determined in a self-consistent way by means of algebra constraints in terms of the external parameters n, U, and T.


### 2.4. Noncrossing Approximation

The calculation of the self-energyΣ
(
k
,
ω
) requires the calculation of the higher-order propagator B
(
k
,
ω
) [cf. (21)]. We will compute this quantity by using the noncrossing approximation (NCA). By neglecting the pair term c
(
i
)
c
†
α
(
i
)
c
(
i
), the source J
(
i
) can be written as
(28)
J
i
,
t
=
∑
j
‍
a
i
,
j
,
t
ψ
j
,
t
,
where
(29)
a
11
i
,
j
,
t
=
-
μ
δ
i
j
-
4
t
α
i
j
-
2
t
σ
μ
n
μ
i
α
i
j
,
a
12
i
,
j
,
t
=
-
4
t
α
i
j
-
2
t
σ
μ
n
μ
i
α
i
j
,
a
21
i
,
j
,
t
=
2
t
σ
μ
n
μ
i
α
i
j
,
a
22
i
,
j
,
t
=
U
-
μ
δ
i
j
+
2
t
σ
μ
n
μ
i
α
i
j
.Then, for the calculation ofB
irr
(
i
,
j
)
=
R
δ
J
i
δ
J
†
j
irr, we approximate
(30)
δ
J
i
,
t
≈
∑
j
‍
a
i
,
j
,
t
-
a
i
,
j
,
t
ψ
j
,
t
.Therefore(31)
B
irr
i
,
j
=
4
t
2
F
i
,
j
1
-
σ
1
,
where we defined
(32)
F
i
,
j
=
R
σ
μ
δ
n
μ
i
c
α
i
c
†
α
j
δ
n
λ
j
σ
λ
with δ
n
μ
(
i
)
=
n
μ
(
i
)
-
n
μ
i. The self-energy (20) is written as
(33)
Σ
k
,
ω
=
4
t
2
F
k
,
ω
I
11
-
2
-
I
11
-
1
I
22
-
1
-
I
11
-
1
I
22
-
1
I
22
-
2
.In order to calculate the retarded functionF
(
i
,
j
), first we use the spectral theorem to express
(34)
F
i
,
j
=
i
2
π
∫
-
∞
+
∞
‍
d
ω
e
-
i
ω
(
t
i
-
t
j
)
×
1
2
π
∫
-
∞
+
∞
‍
d
ω
′
1
+
e
-
β
ω
′
ω
-
ω
′
+
i
ɛ
C
i
-
j
,
ω
′
,
where C
i
-
j
,
ω
′ is the correlation function
(35)
C
i
,
j
=
σ
μ
δ
n
μ
i
c
α
i
c
†
α
j
δ
n
λ
j
σ
λ
=
1
2
π
∫
‍
d
ω
e
-
i
ω
t
i
-
t
j
C
i
-
j
,
ω
.Next, we use the noncrossing approximation (NCA) and approximate(36)
σ
μ
δ
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(37)
F
k
,
ω
=
1
π
∫
-
∞
+
∞
‍
d
ω
′
1
ω
-
ω
′
+
i
δ
a
2
2
π
3
∫
‍
d
2
p
d
Ω
α
2
(
p
)
×
tanh
β
Ω
2
+
coth
⁡
β
ω
′
-
Ω
2
×
I
G
c
c
p
,
Ω
I
χ
k
-
p
,
ω
′
-
Ω
,
where G
c
c
(
k
,
ω
) is the retarded electronic Green’s function [cf. (13)]
(38)
G
c
c
k
,
ω
=
∑
a
,
b
=
1
2
‍
G
a
b
k
,
ω
,
(39)
χ
k
,
ω
=
∑
μ
‍
F
R
δ
n
μ
i
δ
n
μ
j
is the total charge and spin dynamical susceptibility. The result (37) shows that the calculation of the self-energy requires the knowledge of the bosonic propagator (39). This problem will be considered in the following section.It is worth noting that the NCA can also be applied to the casual propagators giving the same result. In general, the knowledge of the self-energy requires the calculation of the higher-order propagatorB
Q
(
i
,
j
)
=
Q
δ
J
i
δ
J
†
j, where Q can be R (retarded propagator) o
T (causal propagator). Then, typically we have to calculate propagator of the form
(40)
H
R
i
,
j
=
R
B
i
F
i
F
†
j
B
†
j
,
where F
(
i
) and B
(
i
) are fermionic and bosonic field operators, respectively. By means of the spectral representation we can write
(41)
H
R
k
,
ω
=
-
1
π
∫
-
∞
+
∞
‍
d
ω
′
1
ω
-
ω
′
+
i
δ
coth
⁡
β
ω
′
2
I
H
C
k
,
ω
′
,
where H
C
(
i
,
j
)
=
T
B
i
F
i
F
†
j
B
†
j is the causal propagator. In the NCA, we approximate
(42)
k
H
C
i
,
j
≈
f
C
i
,
j
b
C
i
,
j
,
f
C
i
,
j
=
T
F
i
F
†
j
,
b
C
i
,
j
=
T
B
i
B
†
j
.
Then, we can use the spectral representation to obtain
(43)
f
C
k
,
ω
=
-
1
π
∫
-
∞
+
∞
‍
d
ω
′
1
-
f
F
ω
′
ω
-
ω
′
+
i
δ
+
f
F
ω
′
ω
-
ω
′
-
i
δ
×
I
f
R
k
,
ω
′
,
b
C
k
,
ω
=
-
1
π
∫
-
∞
+
∞
‍
d
ω
′
1
+
f
B
ω
′
ω
-
ω
′
+
i
δ
-
f
B
ω
′
ω
-
ω
′
-
i
δ
×
I
b
R
k
,
ω
′
which leads to
(44)
H
R
k
,
ω
=
1
π
∫
-
∞
+
∞
‍
d
ω
′
1
ω
-
ω
′
+
i
δ
a
d
2
π
d
+
1
×
∫
Ω
B
‍
d
d
p
d
Ω
I
f
R
p
,
Ω
×
I
b
R
k
-
p
,
ω
′
-
Ω
×
tanh
β
Ω
2
+
coth
⁡
β
ω
′
-
Ω
2
.It is worth noting that, up to this point, the system of equations for the Green’s function and theanomalousself-energy is similar to the one derived in the two-particle self-consistent approach (TPSC) [86, 87, 160], the DMFT
+
Σ approach [163–167], and a Mori-like approach by Plakida and coworkers [6, 7, 9]. It would be the way to compute the dynamical spin and charge susceptibilities to be completely different as, instead of relying on a phenomenological model and neglecting the charge susceptibility as these approches do, we will use a self-consistent two-pole approximation. Obviously, a proper description of the spin and charge dynamics would definitely require the inclusion of a proper self-energy term in the charge and spin propagators too in order to go beyond both any phenomenological approch and the two-pole approximation (in preparation). On the other hand, the description of the electronic anomalous features could (actually will, see in the following) not need this further, and definitely not trivial, complication.


### 2.5. Dynamical Susceptibility

In this section, we will present a calculation of the charge-charge and spin-spin propagators (39) within the two-pole approximation. This approximation has shown to be capable of catching correctly some of the physical features of Hubbard model dynamics (for all details see [40]).Let us define the composite bosonic field:(45)
N
μ
i
=
n
μ
i
ρ
μ
i
,
n
μ
i
=
c
†
i
σ
μ
c
i
,
ρ
μ
i
=
c
†
i
σ
μ
c
α
i
-
c
†
α
i
σ
μ
c
i
.This field satisfies the Heisenberg equation:(46)
i
∂
∂
t
N
μ
i
=
J
μ
i
=
J
1
μ
i
J
2
μ
i
J
1
μ
i
=
-
4
t
ρ
μ
i
,
J
2
μ
i
=
U
κ
μ
i
-
4
t
l
μ
i
,
where the higher-order composite fields κ
μ
(
i
) and l
μ
(
i
) are defined as
(47)
κ
μ
i
=
c
†
i
σ
μ
η
α
i
-
η
†
i
σ
μ
c
α
i
+
η
†
α
i
σ
μ
c
i
-
c
†
α
i
σ
μ
η
i
,
l
μ
i
=
c
†
i
σ
μ
c
α
2
i
+
c
†
α
2
i
σ
μ
c
i
-
2
c
†
α
i
σ
μ
c
α
i
,
and we are using the notation
(48)
c
α
2
i
,
t
=
∑
j
‍
α
i
j
2
c
j
,
t
=
∑
j
l
‍
α
i
l
α
l
j
c
j
,
t
.We linearize the equation of motion (46) for the composite field N
(
μ
)
(
i
) by using the same criterion as in Section 2.2 (i.e., the neglected residual current δ
J
(
μ
)
(
i
) is orthogonal to the chosen basis (45))
(49)
i
∂
∂
t
N
μ
i
,
t
=
∑
j
‍
ɛ
μ
i
,
j
N
μ
j
,
t
,
where the energy matrix is given by
(50)
m
μ
i
,
j
=
∑
l
‍
ɛ
μ
i
,
l
I
μ
l
,
j
,
and the normalization matrix I
(
μ
) and the m
(
μ
) matrix have the following definitions:
(51)
I
μ
i
,
j
=
N
μ
i
,
t
,
N
μ
†
j
,
t
,
m
μ
i
,
j
=
J
μ
i
,
t
,
N
μ
†
j
,
t
.
As it can be easily verified, in the paramagnetic phase, the normalization matrix I
(
μ
) does not depend on the index μ: charge and spin operators have the same weight. The two matrices I
(
μ
) and m
(
μ
) have the following form in momentum space:
(52)
I
μ
k
=
0
I
12
μ
k
I
12
μ
k
0
,
m
μ
k
=
m
11
μ
k
0
0
m
22
μ
k
,
where
(53)
I
12
μ
k
=
4
1
-
α
k
C
c
c
α
,
m
11
μ
k
=
-
4
t
I
12
μ
k
,
m
22
μ
k
=
-
4
t
I
l
μ
ρ
μ
k
+
U
I
κ
μ
ρ
μ
k
.The parameterC
α is the electronic correlation function C
α
=
c
α
i
c
†
i. The quantities I
l
μ
ρ
μ
(
k
) and I
κ
μ
ρ
μ
(
k
) are defined as
(54)
I
l
μ
ρ
μ
k
=
F
l
μ
i
,
t
,
ρ
μ
†
j
,
t
,
I
κ
μ
ρ
μ
k
=
F
κ
μ
i
,
t
,
ρ
μ
†
j
,
t
.
Let us define the causal Green’s function (for bosonic-like fields we have to compute the casual Green’s function and deduce from this latter the retarded one according to the prescriptions in [33–35]),
(55)
G
μ
i
,
j
=
T
N
μ
i
N
μ
†
j
=
i
a
2
2
π
3
∫
‍
d
2
k
d
ω
e
i
k
·
R
i
-
R
j
-
i
ω
t
i
-
t
j
G
μ
k
,
ω
.
By means of the equation of motion (49), the Fourier transform of G
(
μ
)
(
i
,
j
) satisfies the following equation:
(56)
ω
-
ɛ
μ
k
G
μ
k
,
ω
=
I
μ
k
,
where the energy matrix has the explicit form
(57)
ɛ
μ
k
=
0
ɛ
12
μ
k
ɛ
21
μ
k
0
,
ɛ
12
μ
k
=
-
4
t
,
ɛ
21
μ
(
k
)
=
m
22
μ
k
I
12
μ
k
.
The solution of (56) is
(58)
G
μ
k
,
ω
=
Γ
μ
k
1
ω
+
i
δ
-
1
ω
-
i
δ
+
∑
n
=
1
2
‍
σ
n
,
μ
k
1
+
f
B
ω
ω
-
ω
n
μ
k
+
i
δ
-
f
B
ω
ω
-
ω
n
μ
k
-
i
δ
,
where Γ
(
μ
)
(
k
) is the zero frequency function (2
×
2 matrix) [33–35] and f
B
(
ω
)
=
[
e
β
ω
-
1
]
-
1 is the Bose distribution function. Correspondingly, the correlation function C
μ
(
k
,
ω
)
=
N
μ
i
N
μ
†
j has the expression
(59)
C
μ
k
,
ω
=
2
π
Γ
μ
k
δ
ω
+
2
π
∑
n
=
1
2
‍
δ
ω
-
ω
n
μ
k
1
+
f
B
ω
σ
n
,
μ
k
.
The energy spectra ω
n
(
μ
)
(
k
) are given by
(60)
ω
n
μ
k
=
-
n
ω
μ
k
,
ω
μ
k
=
ɛ
12
μ
k
ɛ
21
μ
k
,
and the spectral functions σ
(
n
,
μ
)
(
k
) have the following expression:
(61)
σ
n
,
μ
k
=
I
12
μ
k
2
ɛ
12
μ
k
ω
n
μ
k
1
1
ɛ
21
μ
k
ω
n
μ
k
.
Straightforward but lengthy calculations (see [40]) give for the 2D system the following expressions for the commutators in (54):
(62)
I
l
μ
ρ
μ
k
=
3
4
1
-
α
k
12
C
α
+
C
λ
+
6
C
μ
-
3
1
-
β
k
C
α
+
C
μ
-
3
4
1
-
η
k
C
α
+
C
λ
+
2
C
μ
+
1
4
1
-
λ
k
C
λ
+
3
2
1
-
μ
k
C
μ
,
I
κ
μ
ρ
μ
k
=
-
2
1
-
α
k
D
+
1
-
2
α
k
2
E
β
+
E
η
+
2
β
k
E
β
+
η
k
E
η
+
1
-
2
α
k
a
μ
+
1
4
b
μ
+
2
β
k
c
μ
+
η
k
d
μ
,
where α
(
k
), β
(
k
), η
(
k
), μ
(
k
), and λ
(
k
) are the Fourier transforms of the projectors on the first, second, third, fourth, and fifth nearest-neighbor sites. The parameters appearing in (62) are defined by
(63)
E
=
c
i
η
†
i
,
C
α
=
c
α
i
c
†
i
,
E
β
=
c
β
i
η
†
i
,
C
λ
=
c
λ
i
c
†
i
,
E
η
=
c
η
i
η
†
i
,
C
μ
=
c
μ
i
c
†
i
,
(64)
a
μ
=
2
c
†
i
σ
μ
c
α
i
c
†
i
σ
μ
c
α
i
-
c
†
α
i
σ
μ
σ
λ
σ
μ
c
α
i
n
λ
i
,
b
μ
=
2
c
†
i
σ
μ
c
†
i
σ
μ
c
i
c
i
α
-
c
†
i
σ
μ
σ
λ
σ
μ
c
i
n
λ
α
i
,
c
μ
=
2
c
†
i
σ
μ
c
†
i
η
σ
μ
c
i
α
c
i
α
-
c
†
i
σ
μ
σ
λ
σ
μ
c
i
η
n
λ
i
α
,
d
μ
=
2
c
†
i
σ
μ
c
†
i
β
σ
μ
c
i
α
c
i
α
-
c
†
i
σ
μ
σ
λ
σ
μ
c
i
β
n
λ
i
α
,
where we used the notation
(65)
i
=
i
x
,
i
y
,
t
,
i
β
=
i
x
+
a
,
i
y
+
a
,
t
,
i
α
=
i
x
+
a
,
i
y
,
t
,
i
η
=
i
x
+
2
a
,
i
y
,
t
.We see that the bosonic Green’s functionG
μ
(
i
,
j
)
=
T
N
μ
i
N
μ
†
j depends on the following set of parameters. Fermionic correlators: C
α, C
λ, C
μ, E
β, E
η, and D; bosonic correlators: a
μ, b
μ, c
μ, and d
μ; zero frequency matrix: Γ
(
μ
)
(
k
). The fermionic parameters are calculated through the Fermionic correlation function C
(
i
,
j
)
=
ψ
i
ψ
†
j. The bosonic parameters are determined through symmetry requirements. In particular, the requirement that the continuity equation should be satisfied and that the susceptibility should be a single-value function at k
=
0 leads to the following equations:
(66)
b
μ
=
a
μ
+
3
D
+
2
E
β
+
E
η
-
6
t
U
C
α
+
C
λ
-
2
C
μ
,
c
μ
=
a
μ
-
D
-
2
E
β
+
E
η
+
6
t
U
C
α
+
C
λ
-
2
C
μ
,
d
μ
=
a
μ
-
D
+
2
E
β
-
3
E
η
-
6
t
U
C
α
+
C
λ
-
2
C
μ
.The remaining parametersa
μ and Γ
11
(
μ
)
(
k
) are fixed by means of the Pauli principle
(67)
n
μ
i
n
μ
i
=
n
+
2
D
for
μ
=
0
n
-
2
D
for
μ
=
1,2
,
3
,
where D
=
n
↑
i
n
↓
i is the double occupancy, and by the ergodic value
(68)
Γ
11
μ
k
=
δ
μ
,
0
2
π
2
a
2
δ
2
k
n
2
.
By putting (68) into (58) and (59) we obtain
(69)
δ
n
μ
i
δ
n
μ
j
=
a
2
2
2
π
2
∑
n
=
1
2
‍
∫
‍
d
2
k
e
i
k
·
(
R
i
-
R
j
)
-
i
ω
n
(
μ
)
(
k
)
(
t
i
-
t
j
)
×
1
+
coth
⁡
ω
n
μ
k
2
k
B
T
σ
11
n
,
μ
k
(70)
R
δ
n
μ
i
δ
n
μ
j
=
i
a
2
2
π
3
∑
n
=
1
2
‍
∫
‍
d
2
k
d
ω
e
i
k
·
(
R
i
-
R
j
)
-
i
ω
(
t
i
-
t
j
)
(71)
×
σ
11
n
,
μ
k
ω
-
ω
n
μ
k
+
i
δ
.In conclusion, the dynamical susceptibilityχ
μ
(
k
,
ω
), which is independent of Γ
(
μ
) by construction, reads as
(72)
χ
μ
k
,
ω
=
-
F
R
δ
n
μ
i
δ
n
μ
j
=
16
t
1
-
α
k
C
α
ω
2
-
ω
μ
k
2
,
where no summation is implied on the μ index. It is really worth noticing the very good agreement between the results we obtained within this framework for the charge and spin dynamics of the Hubbard model and the related numerical ones present in the literature (see [40]).


### 2.6. Self-Consistency

In this section, we will give a sketch of the procedure used to calculate the Green’s functionG
(
k
,
ω
). The starting point is (19), where the two matrices I
(
k
) and ɛ
(
k
) are computed by means of the expressions (23). The energy matrix ɛ
(
k
) depends on three parameters: μ, Δ, and p. To determine these parameters we use the following set of algebra constraints:
(73)
n
=
2
1
-
C
11
-
C
22
,
Δ
=
C
11
α
-
C
22
α
,
C
12
=
〈
ξ
(
i
)
η
†
(
i
)
〉
=
0
,
where C
n
m and C
n
m
α are the time-independent correlation functions C
n
m
=
ψ
n
i
ψ
m
†
i and C
n
m
α
=
ψ
n
α
i
ψ
m
†
i. To calculate Σ
(
k
,
ω
) we use the NCA; the results given in Section 2.4 show that within this approximation Σ
(
k
,
ω
) is expressed in terms of the fermionic G
c
c
(
k
,
ω
) [cf. (38)] and of the bosonic χ
μ
(
k
,
ω
) [cf. (39)] propagators. The bosonic propagator is calculated within the two-pole approximation, using the expression (70). As shown in Section 2.5, χ
μ
(
k
,
ω
) depends on both electronic correlation functions [see (63)], which can be straightforwardly computed from G
(
k
,
ω
), and bosonic correlation functions, one per each channel (charge and spin), a
0 and a
3. The latter are determined by means of the local algebra constraints (67), where n is the filling and D is the double occupancy, determined in terms of the electronic correlation function as D
=
n
/
2
-
C
22.According to this, the electronic Green’s functionG
(
k
,
ω
) is computed through the self-consistency scheme depicted in Figure 1: 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
,
ω
). Usually, to get 6 digits precision for fermionic parameters, we need 8 full cycles to reach self-consistency on a 3D grid of 128
×
128 points in momentum space and 4096 Matsubara frequencies. Actually, many more cycles (almost twice) are needed at low doping and low temperatures.Figure 1
(Color online) Self-consistency scheme to compute the propagatorG in terms of the charge-charge and spin-spin propagator B and the residual self-energy Σ.Summarizing, within the NCA and the two-pole approximation for the computation ofχ
μ
(
k
,
ω
) we have constructed an analytical, completely self-consistent scheme of calculation of the electronic propagator G
c
c
(
k
,
ω
)
=
F
T
c
i
c
†
j, where dynamical contributions of the self-energy Σ
(
k
,
ω
) are included. All the internal parameters are self-consistently calculated by means of algebra constraints [cf. (67)-(68) and (73)]. No adjustable parameters or phenomenological expressions are introduced.


## 2.1. Hamiltonian

The Hamiltonian of the two-dimensional Hubbard model reads as(1)
H
=
∑
i
j
‍
-
μ
δ
i
j
-
4
t
α
i
j
c
†
i
c
j
+
U
∑
i
‍
n
↑
i
n
↓
i
,
where
(2)
c
i
=
c
↑
i
c
↓
i
is the electron field operator in spinorial notation and Heisenberg picture (i
=
(
i
,
t
i
)), i
=
R
i is a vector of the Bravais lattice, n
σ
(
i
)
=
c
σ
†
(
i
)
c
σ
(
i
) is the particle density operator for spin σ, n
(
i
)
=
∑
σ
n
σ
(
i
) is the total particle density operator, μ is the chemical potential, t is the hopping integral and the energy unit, U is the Coulomb on-site repulsion, and α
i
j is the projector on the nearest-neighbor sites:
(3)
α
i
j
=
1
N
∑
k
‍
e
i
k
·
R
i
-
R
j
α
k
,
α
k
=
1
2
cos
⁡
k
x
a
+
cos
⁡
k
y
a
,
where k runs over the first Brillouin zone, N is the number of sites, and a is the lattice constant.


## 2.2. Green’s Functions and Dyson Equation

Following COM prescriptions [33–35], we chose a basic field; in particular, we select the composite doublet field operator:
(4)
ψ
i
=
ξ
i
η
i
,
where η
(
i
)
=
n
(
i
)
c
(
i
) and ξ
(
i
)
=
c
(
i
)
-
η
(
i
) are the Hubbard operators describing the main subbands. This choice is guided by the hierarchy of the equations of motion and by the fact that ξ
(
i
) and η
(
i
) are eigenoperators of the interacting term in the Hamiltonian (1). The field ψ
(
i
) satisfies the Heisenberg equation:
(5)
i
∂
∂
t
ψ
i
=
J
i
=
-
μ
ξ
i
-
4
t
c
α
i
-
4
t
π
i
U
-
μ
η
i
+
4
t
π
i
,
where the higher-order composite field π
(
i
) is defined by
(6)
π
i
=
1
2
σ
μ
n
μ
i
c
α
i
+
c
i
c
†
α
i
c
i
with the following notation: n
μ
(
i
)
=
c
†
(
i
)
σ
μ
c
(
i
) is the particle (μ
=
0) and spin (μ
=
1
,
2
,
3) density operator and σ
μ
=
1
,
σ
→, σ
μ
=
-
1
,
σ
→, and σ
k  k
=
1
,
2
,
3 are the Pauli matrices. Hereafter, for any operator Φ
(
i
), we use the notation Φ
α
(
i
,
t
)
=
∑
j
‍
α
i
j
Φ
(
j
,
t
).It is always possible to decompose the sourceJ
(
i
) under the form:
(7)
J
i
=
ɛ
-
i
∇
ψ
i
+
δ
J
i
,
where the linear term represents the projection of the source on the basis ψ
(
i
) and is calculated by means of the equation:
(8)
δ
J
i
,
t
,
ψ
†
j
,
t
=
0
,
where ⋯ stands for the thermal average taken in the grand-canonical ensemble.This constraint assures that the residual currentδ
J
(
i
) contains all and only the physics orthogonal to the chosen basis ψ
(
i
). The action of the derivative operator ɛ
(
-
i
∇
) on ψ
(
i
) is defined in momentum space
(9)
ɛ
-
i
∇
ψ
i
=
ɛ
-
i
∇
1
N
∑
k
‍
e
i
k
·
R
i
ψ
k
,
t
=
1
N
∑
k
‍
e
i
k
·
R
i
ɛ
k
ψ
k
,
t
,
where ɛ
(
k
) is named energy matrix.The constraint (8) gives
(10)
m
k
=
ɛ
k
I
k
after defining the normalization matrix
(11)
I
i
,
j
=
ψ
i
,
t
,
ψ
†
j
,
t
=
1
N
∑
k
‍
e
i
k
·
R
i
-
R
j
I
k
,
and the m-matrix
(12)
m
i
,
j
=
J
i
,
t
,
ψ
†
j
,
t
=
1
N
∑
k
‍
e
i
k
·
R
i
-
R
j
m
k
.Since the components ofψ
(
i
) contain composite operators, the normalization matrix I
(
k
) is not the identity matrix and defines the spectral content of the excitations. In fact, the composite operator method has the advantage of describing crossover phenomena as the phenomena in which the weight of some operator is shifted to another one.By considering the two-time thermodynamic Green’s functions [186–188], let us define the retarded function
(13)
G
i
,
j
=
R
ψ
i
ψ
†
j
=
θ
t
i
-
t
j
ψ
i
,
ψ
†
j
.By means of the Heisenberg equation (5) and using the decomposition (7), the Green’s function G
(
i
,
j
) satisfies the equation
(14)
Λ
∂
i
G
i
,
j
Λ
†
∂
←
j
=
Λ
∂
i
G
0
i
,
j
Λ
†
∂
←
j
+
R
δ
J
i
δ
J
†
j
,
where the derivative operator Λ
(
∂
i
) is defined as
(15)
Λ
∂
i
=
i
∂
∂
t
i
-
ɛ
-
i
∇
i
,
and the propagator G
0
(
i
,
j
) is defined by the equation
(16)
Λ
∂
i
G
0
i
,
j
=
i
δ
t
i
-
t
j
I
i
,
j
.By introducing the Fourier transform(17)
G
i
,
j
=
1
N
∑
k
‍
i
2
π
∫
‍
d
ω
e
i
k
·
R
i
-
R
j
-
i
ω
t
i
-
t
j
G
k
,
ω
equation (14) in momentum space can be written as
(18)
G
k
,
ω
=
G
0
k
,
ω
+
G
0
k
,
ω
I
-
1
k
Σ
k
,
ω
G
k
,
ω
,
and can be formally solved as
(19)
G
k
,
ω
=
1
ω
-
ɛ
k
-
Σ
k
,
ω
I
k
,
where the self-energy Σ
(
k
,
ω
) has the expression
(20)
Σ
k
,
ω
=
B
irr
k
,
ω
I
-
1
k
with
(21)
B
k
,
ω
=
F
R
δ
J
i
δ
J
†
j
.The notationF denotes the Fourier transform and the subscript irr indicates that the irreducible part of the propagator B
(
k
,
ω
) is taken. Equation (18) is nothing else than the Dyson equation for composite fields and represents the starting point for a perturbative calculation in terms of the propagator G
0
(
k
,
ω
). This quantity will be calculated in the next section. Then, the attention will be given to the calculation of the self-energy Σ
(
k
,
ω
). It should be noted that the computation of the two quantities G
0
(
k
,
ω
) and Σ
(
k
,
ω
) is intimately related. The total weight of the self-energy corrections is bounded by the weight of the residual source operator δ
J
(
i
). According to this, it can be made smaller and smaller by increasing the components of the basis ψ
(
i
) (e.g., by including higher-order composite operators appearing in δ
J
(
i
)). The result of such a procedure will be the inclusion in the energy matrix of part of the self-energy as an expansion in terms of coupling constants multiplied by the weights of the newly included basis operators. In general, the enlargement of the basis leads to a new self-energy with a smaller total weight. However, it is necessary pointing out that this process can be quite cumbersome and the inclusion of fully momentum and frequency dependent self-energy corrections can be necessary to effectively take into account low-energy and virtual processes. According to this, one can choose a reasonable number of components for the basic set and then use another approximation method to evaluate the residual dynamical corrections.


## 2.3. Two-Pole Approximation

According to (16), the free propagator G
0
(
k
,
ω
) is determined by the following expression:
(22)
G
0
k
,
ω
=
1
ω
-
ɛ
k
I
k
.For a paramagnetic state, straightforward calculations give the following expressions for the normalizationI
(
k
) and energy ɛ
(
k
) matrices:
(23)
I
k
=
1
-
n
2
0
0
n
2
=
I
11
0
0
I
22
,
ɛ
11
k
=
-
μ
-
4
t
I
11
-
1
Δ
+
1
-
n
+
p
α
k
,
ɛ
12
k
=
4
t
I
22
-
1
Δ
+
p
-
I
22
α
k
,
ɛ
21
k
=
4
t
I
11
-
1
Δ
+
p
-
I
22
α
k
,
ɛ
22
k
=
U
-
μ
-
4
t
I
22
-
1
Δ
+
p
α
k
,
where n
=
n
i is the filling and
(24)
Δ
=
ξ
α
i
ξ
†
i
-
η
α
i
η
†
i
,
p
=
1
4
n
μ
α
i
n
μ
i
-
c
↑
i
c
↓
i
α
c
↓
†
i
c
↑
†
i
.
Then, (22) can be written in spectral form as
(25)
G
0
k
,
ω
=
∑
n
=
1
2
‍
σ
n
k
ω
-
E
n
k
+
i
δ
.The energy spectraE
n
k and the spectral functions σ
(
n
)
(
k
) are given by
(26)
E
1
k
=
R
k
+
Q
k
,
E
2
k
=
R
k
-
Q
k
,
σ
11
1
k
=
I
11
2
1
+
g
k
2
Q
k
,
σ
11
2
k
=
I
11
2
1
-
g
k
2
Q
k
,
σ
12
1
k
=
m
12
k
2
Q
k
,
σ
12
2
k
=
-
m
12
k
2
Q
k
,
σ
22
1
k
=
I
22
2
1
-
g
k
2
Q
k
,
σ
22
2
k
=
I
22
2
1
+
g
k
2
Q
k
,
where
(27)
R
k
=
-
μ
-
4
t
α
k
+
1
2
U
-
ɛ
12
k
2
I
11
,
Q
k
=
1
2
g
2
k
+
4
ɛ
12
2
k
I
22
I
11
,
g
k
=
-
U
+
1
-
n
I
11
ɛ
12
k
.The energy matrixɛ
(
k
) contains three parameters: μ, the chemical potential, Δ, the difference between upper and lower intrasubband contributions to kinetic energy, and p, a combination of the nearest-neighbor charge-charge, spin-spin, and pair-pair correlation functions. These parameters will be determined in a self-consistent way by means of algebra constraints in terms of the external parameters n, U, and T.


## 2.4. Noncrossing Approximation

The calculation of the self-energyΣ
(
k
,
ω
) requires the calculation of the higher-order propagator B
(
k
,
ω
) [cf. (21)]. We will compute this quantity by using the noncrossing approximation (NCA). By neglecting the pair term c
(
i
)
c
†
α
(
i
)
c
(
i
), the source J
(
i
) can be written as
(28)
J
i
,
t
=
∑
j
‍
a
i
,
j
,
t
ψ
j
,
t
,
where
(29)
a
11
i
,
j
,
t
=
-
μ
δ
i
j
-
4
t
α
i
j
-
2
t
σ
μ
n
μ
i
α
i
j
,
a
12
i
,
j
,
t
=
-
4
t
α
i
j
-
2
t
σ
μ
n
μ
i
α
i
j
,
a
21
i
,
j
,
t
=
2
t
σ
μ
n
μ
i
α
i
j
,
a
22
i
,
j
,
t
=
U
-
μ
δ
i
j
+
2
t
σ
μ
n
μ
i
α
i
j
.Then, for the calculation ofB
irr
(
i
,
j
)
=
R
δ
J
i
δ
J
†
j
irr, we approximate
(30)
δ
J
i
,
t
≈
∑
j
‍
a
i
,
j
,
t
-
a
i
,
j
,
t
ψ
j
,
t
.Therefore(31)
B
irr
i
,
j
=
4
t
2
F
i
,
j
1
-
σ
1
,
where we defined
(32)
F
i
,
j
=
R
σ
μ
δ
n
μ
i
c
α
i
c
†
α
j
δ
n
λ
j
σ
λ
with δ
n
μ
(
i
)
=
n
μ
(
i
)
-
n
μ
i. The self-energy (20) is written as
(33)
Σ
k
,
ω
=
4
t
2
F
k
,
ω
I
11
-
2
-
I
11
-
1
I
22
-
1
-
I
11
-
1
I
22
-
1
I
22
-
2
.In order to calculate the retarded functionF
(
i
,
j
), first we use the spectral theorem to express
(34)
F
i
,
j
=
i
2
π
∫
-
∞
+
∞
‍
d
ω
e
-
i
ω
(
t
i
-
t
j
)
×
1
2
π
∫
-
∞
+
∞
‍
d
ω
′
1
+
e
-
β
ω
′
ω
-
ω
′
+
i
ɛ
C
i
-
j
,
ω
′
,
where C
i
-
j
,
ω
′ is the correlation function
(35)
C
i
,
j
=
σ
μ
δ
n
μ
i
c
α
i
c
†
α
j
δ
n
λ
j
σ
λ
=
1
2
π
∫
‍
d
ω
e
-
i
ω
t
i
-
t
j
C
i
-
j
,
ω
.Next, we use the noncrossing approximation (NCA) and approximate(36)
σ
μ
δ
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(37)
F
k
,
ω
=
1
π
∫
-
∞
+
∞
‍
d
ω
′
1
ω
-
ω
′
+
i
δ
a
2
2
π
3
∫
‍
d
2
p
d
Ω
α
2
(
p
)
×
tanh
β
Ω
2
+
coth
⁡
β
ω
′
-
Ω
2
×
I
G
c
c
p
,
Ω
I
χ
k
-
p
,
ω
′
-
Ω
,
where G
c
c
(
k
,
ω
) is the retarded electronic Green’s function [cf. (13)]
(38)
G
c
c
k
,
ω
=
∑
a
,
b
=
1
2
‍
G
a
b
k
,
ω
,
(39)
χ
k
,
ω
=
∑
μ
‍
F
R
δ
n
μ
i
δ
n
μ
j
is the total charge and spin dynamical susceptibility. The result (37) shows that the calculation of the self-energy requires the knowledge of the bosonic propagator (39). This problem will be considered in the following section.It is worth noting that the NCA can also be applied to the casual propagators giving the same result. In general, the knowledge of the self-energy requires the calculation of the higher-order propagatorB
Q
(
i
,
j
)
=
Q
δ
J
i
δ
J
†
j, where Q can be R (retarded propagator) o
T (causal propagator). Then, typically we have to calculate propagator of the form
(40)
H
R
i
,
j
=
R
B
i
F
i
F
†
j
B
†
j
,
where F
(
i
) and B
(
i
) are fermionic and bosonic field operators, respectively. By means of the spectral representation we can write
(41)
H
R
k
,
ω
=
-
1
π
∫
-
∞
+
∞
‍
d
ω
′
1
ω
-
ω
′
+
i
δ
coth
⁡
β
ω
′
2
I
H
C
k
,
ω
′
,
where H
C
(
i
,
j
)
=
T
B
i
F
i
F
†
j
B
†
j is the causal propagator. In the NCA, we approximate
(42)
k
H
C
i
,
j
≈
f
C
i
,
j
b
C
i
,
j
,
f
C
i
,
j
=
T
F
i
F
†
j
,
b
C
i
,
j
=
T
B
i
B
†
j
.
Then, we can use the spectral representation to obtain
(43)
f
C
k
,
ω
=
-
1
π
∫
-
∞
+
∞
‍
d
ω
′
1
-
f
F
ω
′
ω
-
ω
′
+
i
δ
+
f
F
ω
′
ω
-
ω
′
-
i
δ
×
I
f
R
k
,
ω
′
,
b
C
k
,
ω
=
-
1
π
∫
-
∞
+
∞
‍
d
ω
′
1
+
f
B
ω
′
ω
-
ω
′
+
i
δ
-
f
B
ω
′
ω
-
ω
′
-
i
δ
×
I
b
R
k
,
ω
′
which leads to
(44)
H
R
k
,
ω
=
1
π
∫
-
∞
+
∞
‍
d
ω
′
1
ω
-
ω
′
+
i
δ
a
d
2
π
d
+
1
×
∫
Ω
B
‍
d
d
p
d
Ω
I
f
R
p
,
Ω
×
I
b
R
k
-
p
,
ω
′
-
Ω
×
tanh
β
Ω
2
+
coth
⁡
β
ω
′
-
Ω
2
.It is worth noting that, up to this point, the system of equations for the Green’s function and theanomalousself-energy is similar to the one derived in the two-particle self-consistent approach (TPSC) [86, 87, 160], the DMFT
+
Σ approach [163–167], and a Mori-like approach by Plakida and coworkers [6, 7, 9]. It would be the way to compute the dynamical spin and charge susceptibilities to be completely different as, instead of relying on a phenomenological model and neglecting the charge susceptibility as these approches do, we will use a self-consistent two-pole approximation. Obviously, a proper description of the spin and charge dynamics would definitely require the inclusion of a proper self-energy term in the charge and spin propagators too in order to go beyond both any phenomenological approch and the two-pole approximation (in preparation). On the other hand, the description of the electronic anomalous features could (actually will, see in the following) not need this further, and definitely not trivial, complication.


## 2.5. Dynamical Susceptibility

In this section, we will present a calculation of the charge-charge and spin-spin propagators (39) within the two-pole approximation. This approximation has shown to be capable of catching correctly some of the physical features of Hubbard model dynamics (for all details see [40]).Let us define the composite bosonic field:(45)
N
μ
i
=
n
μ
i
ρ
μ
i
,
n
μ
i
=
c
†
i
σ
μ
c
i
,
ρ
μ
i
=
c
†
i
σ
μ
c
α
i
-
c
†
α
i
σ
μ
c
i
.This field satisfies the Heisenberg equation:(46)
i
∂
∂
t
N
μ
i
=
J
μ
i
=
J
1
μ
i
J
2
μ
i
J
1
μ
i
=
-
4
t
ρ
μ
i
,
J
2
μ
i
=
U
κ
μ
i
-
4
t
l
μ
i
,
where the higher-order composite fields κ
μ
(
i
) and l
μ
(
i
) are defined as
(47)
κ
μ
i
=
c
†
i
σ
μ
η
α
i
-
η
†
i
σ
μ
c
α
i
+
η
†
α
i
σ
μ
c
i
-
c
†
α
i
σ
μ
η
i
,
l
μ
i
=
c
†
i
σ
μ
c
α
2
i
+
c
†
α
2
i
σ
μ
c
i
-
2
c
†
α
i
σ
μ
c
α
i
,
and we are using the notation
(48)
c
α
2
i
,
t
=
∑
j
‍
α
i
j
2
c
j
,
t
=
∑
j
l
‍
α
i
l
α
l
j
c
j
,
t
.We linearize the equation of motion (46) for the composite field N
(
μ
)
(
i
) by using the same criterion as in Section 2.2 (i.e., the neglected residual current δ
J
(
μ
)
(
i
) is orthogonal to the chosen basis (45))
(49)
i
∂
∂
t
N
μ
i
,
t
=
∑
j
‍
ɛ
μ
i
,
j
N
μ
j
,
t
,
where the energy matrix is given by
(50)
m
μ
i
,
j
=
∑
l
‍
ɛ
μ
i
,
l
I
μ
l
,
j
,
and the normalization matrix I
(
μ
) and the m
(
μ
) matrix have the following definitions:
(51)
I
μ
i
,
j
=
N
μ
i
,
t
,
N
μ
†
j
,
t
,
m
μ
i
,
j
=
J
μ
i
,
t
,
N
μ
†
j
,
t
.
As it can be easily verified, in the paramagnetic phase, the normalization matrix I
(
μ
) does not depend on the index μ: charge and spin operators have the same weight. The two matrices I
(
μ
) and m
(
μ
) have the following form in momentum space:
(52)
I
μ
k
=
0
I
12
μ
k
I
12
μ
k
0
,
m
μ
k
=
m
11
μ
k
0
0
m
22
μ
k
,
where
(53)
I
12
μ
k
=
4
1
-
α
k
C
c
c
α
,
m
11
μ
k
=
-
4
t
I
12
μ
k
,
m
22
μ
k
=
-
4
t
I
l
μ
ρ
μ
k
+
U
I
κ
μ
ρ
μ
k
.The parameterC
α is the electronic correlation function C
α
=
c
α
i
c
†
i. The quantities I
l
μ
ρ
μ
(
k
) and I
κ
μ
ρ
μ
(
k
) are defined as
(54)
I
l
μ
ρ
μ
k
=
F
l
μ
i
,
t
,
ρ
μ
†
j
,
t
,
I
κ
μ
ρ
μ
k
=
F
κ
μ
i
,
t
,
ρ
μ
†
j
,
t
.
Let us define the causal Green’s function (for bosonic-like fields we have to compute the casual Green’s function and deduce from this latter the retarded one according to the prescriptions in [33–35]),
(55)
G
μ
i
,
j
=
T
N
μ
i
N
μ
†
j
=
i
a
2
2
π
3
∫
‍
d
2
k
d
ω
e
i
k
·
R
i
-
R
j
-
i
ω
t
i
-
t
j
G
μ
k
,
ω
.
By means of the equation of motion (49), the Fourier transform of G
(
μ
)
(
i
,
j
) satisfies the following equation:
(56)
ω
-
ɛ
μ
k
G
μ
k
,
ω
=
I
μ
k
,
where the energy matrix has the explicit form
(57)
ɛ
μ
k
=
0
ɛ
12
μ
k
ɛ
21
μ
k
0
,
ɛ
12
μ
k
=
-
4
t
,
ɛ
21
μ
(
k
)
=
m
22
μ
k
I
12
μ
k
.
The solution of (56) is
(58)
G
μ
k
,
ω
=
Γ
μ
k
1
ω
+
i
δ
-
1
ω
-
i
δ
+
∑
n
=
1
2
‍
σ
n
,
μ
k
1
+
f
B
ω
ω
-
ω
n
μ
k
+
i
δ
-
f
B
ω
ω
-
ω
n
μ
k
-
i
δ
,
where Γ
(
μ
)
(
k
) is the zero frequency function (2
×
2 matrix) [33–35] and f
B
(
ω
)
=
[
e
β
ω
-
1
]
-
1 is the Bose distribution function. Correspondingly, the correlation function C
μ
(
k
,
ω
)
=
N
μ
i
N
μ
†
j has the expression
(59)
C
μ
k
,
ω
=
2
π
Γ
μ
k
δ
ω
+
2
π
∑
n
=
1
2
‍
δ
ω
-
ω
n
μ
k
1
+
f
B
ω
σ
n
,
μ
k
.
The energy spectra ω
n
(
μ
)
(
k
) are given by
(60)
ω
n
μ
k
=
-
n
ω
μ
k
,
ω
μ
k
=
ɛ
12
μ
k
ɛ
21
μ
k
,
and the spectral functions σ
(
n
,
μ
)
(
k
) have the following expression:
(61)
σ
n
,
μ
k
=
I
12
μ
k
2
ɛ
12
μ
k
ω
n
μ
k
1
1
ɛ
21
μ
k
ω
n
μ
k
.
Straightforward but lengthy calculations (see [40]) give for the 2D system the following expressions for the commutators in (54):
(62)
I
l
μ
ρ
μ
k
=
3
4
1
-
α
k
12
C
α
+
C
λ
+
6
C
μ
-
3
1
-
β
k
C
α
+
C
μ
-
3
4
1
-
η
k
C
α
+
C
λ
+
2
C
μ
+
1
4
1
-
λ
k
C
λ
+
3
2
1
-
μ
k
C
μ
,
I
κ
μ
ρ
μ
k
=
-
2
1
-
α
k
D
+
1
-
2
α
k
2
E
β
+
E
η
+
2
β
k
E
β
+
η
k
E
η
+
1
-
2
α
k
a
μ
+
1
4
b
μ
+
2
β
k
c
μ
+
η
k
d
μ
,
where α
(
k
), β
(
k
), η
(
k
), μ
(
k
), and λ
(
k
) are the Fourier transforms of the projectors on the first, second, third, fourth, and fifth nearest-neighbor sites. The parameters appearing in (62) are defined by
(63)
E
=
c
i
η
†
i
,
C
α
=
c
α
i
c
†
i
,
E
β
=
c
β
i
η
†
i
,
C
λ
=
c
λ
i
c
†
i
,
E
η
=
c
η
i
η
†
i
,
C
μ
=
c
μ
i
c
†
i
,
(64)
a
μ
=
2
c
†
i
σ
μ
c
α
i
c
†
i
σ
μ
c
α
i
-
c
†
α
i
σ
μ
σ
λ
σ
μ
c
α
i
n
λ
i
,
b
μ
=
2
c
†
i
σ
μ
c
†
i
σ
μ
c
i
c
i
α
-
c
†
i
σ
μ
σ
λ
σ
μ
c
i
n
λ
α
i
,
c
μ
=
2
c
†
i
σ
μ
c
†
i
η
σ
μ
c
i
α
c
i
α
-
c
†
i
σ
μ
σ
λ
σ
μ
c
i
η
n
λ
i
α
,
d
μ
=
2
c
†
i
σ
μ
c
†
i
β
σ
μ
c
i
α
c
i
α
-
c
†
i
σ
μ
σ
λ
σ
μ
c
i
β
n
λ
i
α
,
where we used the notation
(65)
i
=
i
x
,
i
y
,
t
,
i
β
=
i
x
+
a
,
i
y
+
a
,
t
,
i
α
=
i
x
+
a
,
i
y
,
t
,
i
η
=
i
x
+
2
a
,
i
y
,
t
.We see that the bosonic Green’s functionG
μ
(
i
,
j
)
=
T
N
μ
i
N
μ
†
j depends on the following set of parameters. Fermionic correlators: C
α, C
λ, C
μ, E
β, E
η, and D; bosonic correlators: a
μ, b
μ, c
μ, and d
μ; zero frequency matrix: Γ
(
μ
)
(
k
). The fermionic parameters are calculated through the Fermionic correlation function C
(
i
,
j
)
=
ψ
i
ψ
†
j. The bosonic parameters are determined through symmetry requirements. In particular, the requirement that the continuity equation should be satisfied and that the susceptibility should be a single-value function at k
=
0 leads to the following equations:
(66)
b
μ
=
a
μ
+
3
D
+
2
E
β
+
E
η
-
6
t
U
C
α
+
C
λ
-
2
C
μ
,
c
μ
=
a
μ
-
D
-
2
E
β
+
E
η
+
6
t
U
C
α
+
C
λ
-
2
C
μ
,
d
μ
=
a
μ
-
D
+
2
E
β
-
3
E
η
-
6
t
U
C
α
+
C
λ
-
2
C
μ
.The remaining parametersa
μ and Γ
11
(
μ
)
(
k
) are fixed by means of the Pauli principle
(67)
n
μ
i
n
μ
i
=
n
+
2
D
for
μ
=
0
n
-
2
D
for
μ
=
1,2
,
3
,
where D
=
n
↑
i
n
↓
i is the double occupancy, and by the ergodic value
(68)
Γ
11
μ
k
=
δ
μ
,
0
2
π
2
a
2
δ
2
k
n
2
.
By putting (68) into (58) and (59) we obtain
(69)
δ
n
μ
i
δ
n
μ
j
=
a
2
2
2
π
2
∑
n
=
1
2
‍
∫
‍
d
2
k
e
i
k
·
(
R
i
-
R
j
)
-
i
ω
n
(
μ
)
(
k
)
(
t
i
-
t
j
)
×
1
+
coth
⁡
ω
n
μ
k
2
k
B
T
σ
11
n
,
μ
k
(70)
R
δ
n
μ
i
δ
n
μ
j
=
i
a
2
2
π
3
∑
n
=
1
2
‍
∫
‍
d
2
k
d
ω
e
i
k
·
(
R
i
-
R
j
)
-
i
ω
(
t
i
-
t
j
)
(71)
×
σ
11
n
,
μ
k
ω
-
ω
n
μ
k
+
i
δ
.In conclusion, the dynamical susceptibilityχ
μ
(
k
,
ω
), which is independent of Γ
(
μ
) by construction, reads as
(72)
χ
μ
k
,
ω
=
-
F
R
δ
n
μ
i
δ
n
μ
j
=
16
t
1
-
α
k
C
α
ω
2
-
ω
μ
k
2
,
where no summation is implied on the μ index. It is really worth noticing the very good agreement between the results we obtained within this framework for the charge and spin dynamics of the Hubbard model and the related numerical ones present in the literature (see [40]).


## 2.6. Self-Consistency

In this section, we will give a sketch of the procedure used to calculate the Green’s functionG
(
k
,
ω
). The starting point is (19), where the two matrices I
(
k
) and ɛ
(
k
) are computed by means of the expressions (23). The energy matrix ɛ
(
k
) depends on three parameters: μ, Δ, and p. To determine these parameters we use the following set of algebra constraints:
(73)
n
=
2
1
-
C
11
-
C
22
,
Δ
=
C
11
α
-
C
22
α
,
C
12
=
〈
ξ
(
i
)
η
†
(
i
)
〉
=
0
,
where C
n
m and C
n
m
α are the time-independent correlation functions C
n
m
=
ψ
n
i
ψ
m
†
i and C
n
m
α
=
ψ
n
α
i
ψ
m
†
i. To calculate Σ
(
k
,
ω
) we use the NCA; the results given in Section 2.4 show that within this approximation Σ
(
k
,
ω
) is expressed in terms of the fermionic G
c
c
(
k
,
ω
) [cf. (38)] and of the bosonic χ
μ
(
k
,
ω
) [cf. (39)] propagators. The bosonic propagator is calculated within the two-pole approximation, using the expression (70). As shown in Section 2.5, χ
μ
(
k
,
ω
) depends on both electronic correlation functions [see (63)], which can be straightforwardly computed from G
(
k
,
ω
), and bosonic correlation functions, one per each channel (charge and spin), a
0 and a
3. The latter are determined by means of the local algebra constraints (67), where n is the filling and D is the double occupancy, determined in terms of the electronic correlation function as D
=
n
/
2
-
C
22.According to this, the electronic Green’s functionG
(
k
,
ω
) is computed through the self-consistency scheme depicted in Figure 1: 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
,
ω
). Usually, to get 6 digits precision for fermionic parameters, we need 8 full cycles to reach self-consistency on a 3D grid of 128
×
128 points in momentum space and 4096 Matsubara frequencies. Actually, many more cycles (almost twice) are needed at low doping and low temperatures.Figure 1
(Color online) Self-consistency scheme to compute the propagatorG in terms of the charge-charge and spin-spin propagator B and the residual self-energy Σ.Summarizing, within the NCA and the two-pole approximation for the computation ofχ
μ
(
k
,
ω
) we have constructed an analytical, completely self-consistent scheme of calculation of the electronic propagator G
c
c
(
k
,
ω
)
=
F
T
c
i
c
†
j, where dynamical contributions of the self-energy Σ
(
k
,
ω
) are included. All the internal parameters are self-consistently calculated by means of algebra constraints [cf. (67)-(68) and (73)]. No adjustable parameters or phenomenological expressions are introduced.


## 3. Results

In the following, we analyze some electronic properties by computing the spectral function(74)
A
k
,
ω
=
-
1
π
I
G
c
c
k
,
ω
,
the momentum distribution function per spin
(75)
n
k
=
∫
‍
d
ω
f
F
ω
A
k
,
ω
,
and the density of states per spin
(76)
N
ω
=
1
2
π
2
∫
‍
d
2
k
A
k
,
ω
,
where G
c
c
(
k
,
ω
)
=
G
11
(
k
,
ω
)
+
G
12
(
k
,
ω
)
+
G
21
(
k
,
ω
)
+
G
22
(
k
,
ω
) is the electronic propagator and f
F
(
ω
) is the Fermi function. We also study the electronic self-energy Σ
c
c
(
k
,
ω
), which is defined through the equation
(77)
G
c
c
k
,
ω
=
1
ω
-
ϵ
0
k
-
Σ
c
c
k
,
ω
,
where ϵ
0
(
k
)
=
-
μ
-
4
t
α
(
k
) is the noninteracting dispersion. Moreover, we define the quantity r
(
k
)
=
ϵ
0
(
k
)
+
Σ
c
c
′
(
k
,
ω
=
0
) that determines the Fermi surface locus in momentum space, r
(
k
)
=
0, in a Fermi liquid, that is, when lim
⁡
ω
→
0
Σ
c
c
′
′
(
k
,
ω
,
T
=
0
)
∝
ω
2 and lim
⁡
T
→
0
Σ
c
c
′
′
(
k
,
ω
=
0
,
T
)
∝
T
2. 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, at the same time and on equal footing, both the real and the imaginary parts of the self-energy and is directly related, within the sudden approximation and forgetting any selection rules, to what ARPES effectively measures.Finally, the spin-spin correlation functionn
z
n
z
α, the pole ω
(
3
)
(
k
=
Q
=
(
π
,
π
)
) of the spin-spin propagator, and the antiferromagnetic correlation length ξ are discussed. Usually, the latter is defined by supposing the following asymptotic expression for the static susceptibility:
(78)
lim
⁡
k
→
Q
χ
(
3
)
k
,
0
=
χ
(
3
)
Q
,
0
1
+
ξ
2
k
-
Q
2
,
where χ
(
3
)
k
,
0
=
-
G
(
3
)
k
,
0. It is worth noting that in our case (78) is not assumed, but it exactly holds [189].


### 3.1. Spectral Function and Dispersion

According to its overall relevance in the whole analysis performed hereinafter, we first discuss the electronic dispersion of the model under analysis or, better, its relic in a strongly correlated system. In general, the dispersion and its more or less anomalous features can be inferred by looking at the maxima of the spectral functionA
(
k
,
ω
). In Figures 2 and 3, the spectral function is shown, in scale of grays (increasing from white to black; red is for above-scale values), along the principal directions (Γ
=
(
0
,
0
)
→
M
=
(
π
,
π
), M
→
X
=
(
π
,
0
), X
→
Y
=
(
0
,
π
) and Y
→
Γ) for U
=
8, T
=
0.01 and (a) n
=
0.70, (b) n
=
0.78, (c) n
=
0.85, and (d) n
=
0.92 (T
=
0.02). In Figure 2, the whole range of frequencies with finite values of A
(
k
,
ω
) is reported, while in Figure 3 a zoom in the proximity of the chemical potential is shown. The light gray lines and uniform areas are labeled with the values of the imaginary part of the self-energy Σ
′
′
(
k
,
ω
) and give one of the most relevant keys to interpret the characteristics of the dispersion. The dark green lines in Figure 3 are just guides to the eye and indicate the direction of the dispersion justbefore the visible kink separating the black and the red areas of the dispersion.(Color online) Spectral functionA
(
k
,
ω
) along the principal directions (Γ
=
(
0
,
0
)
→
M
=
(
π
,
π
), M
→
X
=
(
π
,
0
), X
→
Y
=
(
0
,
π
) and Y
→
Γ) for U
=
8, T
=
0.01 and (a) n
=
0.70, (b) n
=
0.78, (c) n
=
0.85, and (d) n
=
0.92 (T
=
0.02).
(a)
(b)
(c)
(d)(Color online) Spectral functionA
(
k
,
ω
) close to the chemical potential (ω
=
0) along the principal directions (Γ
=
(
0
,
0
)
→
M
=
(
π
,
π
), M
→
X
=
(
π
,
0
), X
→
Y
=
(
0
,
π
) and Y
→
Γ) for U
=
8, T
=
0.01 and (a) n
=
0.70, (b) n
=
0.78, (c) n
=
0.85, and (d) n
=
0.92 (T
=
0.02).
(a)
(b)
(c)
(d)The red areas, as they mark the relative maxima ofA
(
k
,
ω
), can be considered as the best possible estimates for the dispersion. In a non- (or weakly-) interacting system, the dispersion would be a single, continuos, and (quite-) sharp line representing some function ɛ
k being the simple pole of G
c
c
(
k
,
ω
). In this case (for a strongly correlated system), instead, we can clearly see that the dispersion is well-defined (red areas) only in some of the regions it crosses in the (
k
,
ω
) plane: the regions where Σ
′
′
(
k
,
ω
) is zero or almost negligible. In the crossed regions where Σ
′
′
(
k
,
ω
) is instead finite, A
(
k
,
ω
) obviously assumes very low values, which would be extremely difficult (actually almost impossible) to detect by ARPES. Accordingly, ARPES would report only the red areas in the picture. This fact is fundamental to understand and describe the experimental findings regarding the Fermi surface in the underdoped regime, as they will be discussed in the next section, and to reconcile ARPES findings with those of quantum oscillations measurements.The two Hubbard subbands, separated by a gap of the orderU and with a reduced bandwidth of order 4
t, are clearly visible: the lower one (LHB) is partly occupied as it is crossed by the chemical potential and the upper one (UHB) is empty and very far from the chemical potential. The lower subband systematically (for each value of doping) loses significance close to Γ and to M, although this effect is more and more pronounced on reducing doping. In both cases (close to Γ and to M), A
(
k
,
ω
) loses weight as Σ
′
′
(
k
,
ω
) increases its own: this can be easily understood if we recall that both χ
0
k
,
ω and χ
3
k
,
ω have a vanishing pole at Γ (due to hydrodynamics) and that χ
3
k
,
ω is strongly peaked at M due to the strong antiferromagnetic correlations present in the system (see Section 3.6). The upper subband, according to the complementary effect induced by the evident shadow-bands appearance (signaled by the relative maxima of Σ
′
′
(
k
,
ω
) marked by 100 and the dark gray area inside the gap) due again to the strong antiferromagnetic correlations present in the system, displays a well-defined dispersion at M, at least for high enough doping, and practically no dispersion at all close to Γ. For low doping, the great majority of the upper subband weight is simply transferred to the lower subband. The growth, on reducing doping, of the undefined-dispersion regions close to Γ in the lower subband cuts down its already reduced bandwidth of order 4
t to values of the order J
=
4
t
2
/
U
=
0.5
t, as one would expect for the dispersion of few holes in a strong antiferromagnetic background. The shape of the dispersion is too compatible with this scenario: the sequence of minima and maxima is compatible, actually driven, by the doubling of the Brillouin zone induced by the strong antiferromagnetic correlations, as well as the dynamical generation of a t
′ diagonal hopping (absent in the Hamiltonian currently under study) clearly signaled by the more and more, on reducing doping, pronounced warping of the dispersion along the X
→
Y direction (more evident in Figure 3), which would be perfectly flat otherwise (for t
′
=
0).Moving to the zooms (Figure3), they show much more clearly the systematic reduction of the bandwidth on reducing doping, the doubling of the zone, the systematic increase of the warping along X
→
Y on reducing doping, and the extreme flatness of the dispersion at X coming from both Γ and M. This latter feature is in very good agreement with quantum Monte Carlo calculations (see [159] and references therein) as well as with ARPES experiments [190], which report a similar behavior in the overdoped region. Moreover, they show that, in contrast with the scenario for a non- (or weakly-) interacting system, where the doubling of the zone happens with the X
→
Y direction aspivot; here, the doubling is confined to the region close to M. The lower subband is completely filled at half-filling, while for an ordinary Slater antiferromagnet the gap opens at half-filling just on top of the van-Hove singularity along X
→
Y. In addition, the effective, finite value of t
′ makes the dispersion maximum close to S higher than the one present along the M
→
X direction, opening the possibility for the appearance of hole pockets close to S. Finally, it is now evident that the warping of the dispersion along X
→
Y will also induce the presence of two maxima in the density of states: one due to the van-Hove singularities at X and Y and one due to the maximum in the dispersion close to S. How deep is the dip between these two maxima just depends on the number of available well-defined (red) states in momentum present between these two values of frequency. This will determine the appearance of a more or less pronounced pseudogap in the density of states, but let us come back to this after having analyzed the region close to M.As a matter of fact, it is just the absence of spectral weight in the region close toM and, in particular and more surprisingly, at the chemical potential (i.e., on the Fermi surface, in contradiction with the Fermi-liquid picture), the main and more relevant result of this analysis, it will determine almost all interesting and anomalous/unconventional features of the single-particle properties of the model. The scenario emerging from this analysis can be relevant not only for the understanding of the physics of the Hubbard model and for the microscopical description of the cuprate high-T
c superconductors, but also for the drafting of a general microscopic theory of strongly correlated systems. The strong antiferromagnetic correlations (through Σ
′
′
(
k
,
ω
), which mainly follows χ
3
k
,
ω) cause a significative and anomalous/unconventional loss of spectral weight around M, which induces in turn the deconstruction of the Fermi surface (Section 3.2), the emergence of momentum selective non-Fermi-liquid features (Sections 3.4 and 3.5), and the opening of a well-developed (deep) pseudogap in the density of states (Section 3.3).Last, but not least, it is also remarkable 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 [82]. Such a phenomenon clearly signals the coupling of the electrons to a bosonic mode. In this scenario, the mode is clearlymagnetic in nature. The frequency of the kink, with respect to the chemical potential, reduces systematically and quite drastically on reducing doping, following the behavior of the pole of χ
3
k
,
ω (see Section 3.6). Combining the presence of kinks and the strong reduction of spectral weight below them, we see the appearance of waterfalls, in particular along the antinodal (X
→
Γ) direction, as found in some ARPES experiments [82]. Finally, the extension of the flat region in the dispersion around the antinodal points (X and Y), that is, at the van-Hove points, increases systematically on decreasing doping. This clearly signals the transfer of spectral weight from the Fermi surface, which is depleted by the strong antiferromagnetic fluctuations, which are also responsible for the remarkable flatness of the band edge.Before moving to the next section, it is worth noticing that similar results for the single-particle excitation spectrum (flat bands close toX, weight transfer from the LHB to the UHB at M, and splitting of the band close to X
,
…) were obtained within the self-consistent projection operator method [28, 29], the operator projection method [30–32], and within a Mori-like approach by Plakida and coworkers [7, 9].


### 3.2. Spectral Function and Fermi Surface

Focusing on the value of the spectral function at the chemical potential,A
(
k
,
ω
=
0
), we can discuss the closest concept to Fermi surface available in a strongly correlated system. In Figure 4, A
(
k
,
ω
=
0
) is plotted as a function of the momentum k in a quarter of the Brillouin zone for U
=
8 and four different couples of values of temperature and filling: (a) n
=
0.7 and T
=
0.01, (b) n
=
0.78 and T
=
0.01, (c) n
=
0.85 and T
=
0.01, and (d) n
=
0.92 and T
=
0.02. The Fermi surface, in agreement with the interpretation of the ARPES measurements within the sudden approximation, can be defined as the locus in momentum space of the relative maxima of A
(
k
,
ω
=
0
). Such a definition opens up the possibility of not only explaining ARPES measurements, but also going beyond them and their finite instrumental resolution and sensitivity with the aim at filling the gap with other kinds of measurements, in particular quantum oscillations ones, which seems to report results in disagreement, up to dichotomy in some cases, with the scenario depicted by ARPES.(Color online) Spectral function at the chemical potentialA
(
k
,
ω
=
0
) as a function of momentum k for U
=
8, T
=
0.01 and (a) n
=
0.70, (b) n
=
0.78, (c) n
=
0.85, and (d) n
=
0.92 (T
=
0.02).
(a)
(b)
(c)
(d)For each value of the filling reported, we can easily distinguish two walls/arcs; forn
=
0.92, they somewhat join. First, let us focus on the arc with the larger (by far) intensities; given the current sensitivities, this is the only one, among the two, possibly visible to ARPES. At n
=
0.7 (see Figure 4(a)), the 3D perspective allows to better appreciate the difference in the intensities of the relative maxima of A
(
k
,
ω
=
0
) between the region close to the main diagonal (M
→
X), where the signal is weaker, and the regions close to the main axes (Γ
→
X and Γ
→
Y), where the signal is stronger; this behavior has been also reported by ARPES experiments [82, 190] as well as the electron-like nature of the Fermi surface [190]. On decreasing doping, this trend reverses, passing through n
=
0.78, where the intensities almost match, and up to n
=
0.92, where the region in proximity of S is the only one with an appreciable signal. The less-intense (by far) arc, reported in [7] too, is the relic of a shadow band, as can be clearly seen in Figure 3, and, consequently, never changes its curvature, in contrast to what happens to the other arc, which is subject to the crossing of the van Hove singularity (n
≅
0.82) instead. Although the ratio between the maximum values of the intensities at the two arcs never goes below two (see Figure 6(b)), there is an evident decrease of the maximum value of the intensity at the larger-intensity arc on decreasing doping, which clearly signals an overall increase of the intensity and/or of the effectiveness (in terms of capability of affecting the relevant quasi-particles, which are those at the Fermi surface) of the correlations.Moving to a 2D perspective (see Figure5), we can add three ingredients to our discussion that can help us better understand the evolution with doping of the Fermi surface: (i) the n
(
k
)
=
0.5 locus (solid line), that is, the Fermi surface if the system would be noninteracting; (ii) the r
(
k
)
=
0 locus (dashed line), that is, the Fermi surface if the system would be a Fermi liquid or somewhat close to it conceptually; (iii) the values (grey lines and labels) of the imaginary part of the self-energy at the chemical potential Σ
c
c
′
′
(
k
,
ω
=
0
) (notice that T
≠
0). Combining these three ingredients with the positions and intensities of the the relative maxima of A
(
k
,
ω
=
0
), we can try to better understand what these latter signify and to classify the behavior of the system on changing doping. At n
=
0.7 (see Figure 5(a)), the positions of the two arcs are exactly matching r
(
k
)
=
0 lines; this will stay valid at each value of the filling reported, with a fine, but very relevant, distinction at n
=
0.92. This occurrence makes our definition of Fermi surface robust, but also versatile as it permits to go beyond Fermi liquid picture without contradicting this latter. The almost perfect coincidence, for the higher value of doping reported, n
=
0.7, of the n
(
k
)
=
0.5 line with the larger-intensity arc clearly asserts that we are dealing with a very-weakly-interacting Fermi metal. Σ
c
c
′
′
(
k
,
ω
=
0
) is quite large close to M (see Figure 3(a)) and eats up the weight of the second arc, which becomes aghost band more than ashadow one. The antiferromagnetic correlations are definitely finite (see Section 3.6) and, consequently, lead to the doubling of the zone, but not strong enough to affect the behavior of the ordinary quasiparticles safely living at the ordinary Fermi surface. Decreasing the doping, we can witness a first topological transition from a Fermi surface closed around Γ (electron-like, hole-like in cuprates language) to a Fermi surface closed around M (hole-like, electron-like in cuprates language) at n
≅
0.82, where the chemical potential crosses the van Hove singularity (see Figure 3). The chemical potential presents an inflection point at this doping (not shown), which allowed us to determine its value with great accuracy. In proximity of the antinodal points (X and Y), a net discrepancy between the position of the relative maxima of A
(
k
,
ω
=
0
) and the locus n
(
k
)
=
0.5 becomes more and more evident on decreasing doping (see Figures 5(b) and 5(c)). This occurrence does not only allow the topological transition, which is absent for the n
(
k
)
=
0.5 locus that reaches the antidiagonal (X
→
Y) at half-filling in agreement with the Luttinger theorem, but also accounts for the apparent broadening of the relative maxima of A
(
k
,
ω
=
0
) close to the antinodal points (X and Y). The broadening is due to the small, but finite, value of Σ
c
c
′
′
(
k
,
ω
=
0
) in those momentum regions (see Figures 3(b) and 3(c)), signaling the net increase of the correlation strength and, accordingly, the impossibility of describing the system in this regime as a conventional non- (or weakly-) interacting system within a Fermi-liquid scenario or its ordinary extensions for ordered phases. What is really interesting and goes beyond the actual problem under analysis (cuprates: 2D Hubbard model) is the emergence of such features only in well-defined regions in momentum space. This selectiveness in momentum is quite a new feature in condensed matter physics and its understanding and description require quite new theoretical approaches. In this system, almost independently from their effective strength, the correlations play a so fundamental role to come to shape and determinequalitatively the response of the system. Accordingly, any attempt to treat correlations without taking into account the level of entanglement between all degrees of freedom present in the system is bound to fail or at least to miss the most relevant features.(Color online) Spectral function at the chemical potentialA
(
k
,
ω
=
0
) as a function of momentum k for U
=
8, T
=
0.01 and (a) n
=
0.70, (b) n
=
0.78, (c) n
=
0.85, and (d) n
=
0.92 (T
=
0.02). The solid line marks the locus n
(
k
)
=
0.5, the dashed line marks the locus r
(
k
)
=
0, the gray lines are labeled with the values of Σ
c
c
′
′
(
k
,
ω
=
0
), and the dotted line is a guide to the eye and marks the reduced (antiferromagnetic) Brillouin zone.
(a)
(b)
(c)
(d)(Color online) (a) Density of statesN
(
ω
) as a function of frequency ω for U
=
8, (black squares) n
=
0.7 and T
=
0.01, (red circles) n
=
0.78 and T
=
0.01, (blue up triangles) n
=
0.85 and T
=
0.01, and (green down triangles) n
=
0.92 and T
=
0.02. (b) Spectral function in the proximity of the chemical potential A
(
k
_
,
ω
~
0
) at (black squares) k
_
=
S
=
(
π
/
2
,
π
/
2
), (red circles) S
_ (in the text), and (blue triangles) X
=
(
π
,
0
) for U
=
8, n
=
0.92, and T
=
0.02.
(a)
(b)Let us come to the most interesting result. Atn
=
0.92, the relative maxima of A
(
k
,
ω
=
0
) detach also from the r
(
k
)
=
0 locus, at least partially, opening a completely new scenario. r
(
k
)
=
0 defines a pocket (close, but absolutely not identical, read below, to that of an antiferromagnet), while the relative maxima of A
(
k
,
ω
=
0
) feature the very same pocket together with quite broad, but still well-defined, wings closing with one half of the pocket (the most intense one) what can be safely considered the relic of a large Fermi surface. This is the second and most surprising topological transition occurring to the actual Fermi surface: the two arcs, clearly visible for all other values of the filling, join and instead of closing just a pocket, as one would expect on the basis of the conventional theory for an antiferromagnet—here mimed by r
(
k
)
=
0 locus, develop (or keep) a completely independent branch. The actual Fermi surface is neither a pocket nor a large Fermi surface; for a more expressive representation, see Figure 4(d). This very unexpected result can be connected to the dichotomy between those experiments (e.g., ARPES) pointing to a small and those ones (e.g., quantum oscillations) pointing to a large Fermi surface. This result can be understood by looking once more at the dispersion for this value of the filling (see Figure 3(d)): the difference in height, induced by the effective finite value of t
′, which increases with decreasing doping, between the two highest maxima in the dispersion—one close to S and the other along the M
→
X direction—makes the latter to cross the chemical potential for larger values of the doping than the first, but given the significative broadening of the dispersion, even when the center mass of the second leaves the Fermi surface (disappearing from r
(
k
)
=
0 locus), its shoulders are still active and well-identifiable at the chemical potential—that is on the actual Fermi surface.The pocket is too far from beingconventional. It is clearly evident in Figure 5(d) that there are two distinct halves of the pocket: one with very high intensity pinned at S (again the only possibly visible to ARPES) and another with very low intensity (visible only totheoreticians and some quantum oscillations experiments). This is our interpretation for the Fermi arcs reported by many ARPES experiments [82, 93] and unaccountable for any ordinary theory relaying on the Fermi liquid picture, though being modified by the presence of an incipient spin or charge ordering. Obviously, looking only at the Fermi arc (as ARPES is forced to do), the Fermi surface looks ill defined as it does not enclose a definite region of momentum space, but having access also to the other half of the pocket, such a problem is greatly alleviated. The point is that the antiferromagnetic fluctuations are so strong to destroy the coherence of the quasiparticles in that region of momentum space as similarly reported within the DMFT
+
Σ approach [163–167] and a Mori-like approach by Plakida and coworkers [7, 9]. Moreover, we will see that (see Section 3.5) thephantom half is not simply lower in intensity because it belongs to a shadow band depleted by a finite imaginary part of the self-energy (see Figure 5(d)) but that it lives in a region of momentum where the imaginary part of the self-energy shows clear signs of non-Fermi-liquid behavior. The lack of next-nearest hopping terms (i.e., t
′ and t
′
′) in the chosen Hamiltonian (1), although they are evidently generated dynamically, does not allow us to perform a quantitative comparison between our results and the experimental ones, which refers to a specific material characterized by a specific set of hoppings. This also explains why our Fermi arc is pinned at S and occupies the outer reduced Brillouin zone, while many experimental results that report a Fermi arc occupying the inner reduced Brillouin zone. Actually, the pinning (with respect to doping within the underdoped region) of thecenter of mass of the ARPES-visible Fermi arc has been reported also from ARPES experiments [191].


### 3.3. Density of States and Pseudogap

The other main issue in the underdoped regime of cuprates superconductors is the presence of a quite strong depletion in the electronic density of states, known as pseudogap. In Figure6(a), we report the density of states 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 in proximity of the chemical potential. As a reference, we also report, in Figure 6(b), the spectral function in proximity of the chemical potential A
(
k
_
,
ω
~
0
) at k
_
=
S
=
(
π
/
2
,
π
/
2
), k
_
=
S
_ which lies where thephantom half of the pocket touches the main diagonal Γ
→
M (i.e., where the dispersion cuts the main diagonal Γ
→
M closer to M), and k
_
=
X
=
(
π
,
0
) for U
=
8, n
=
0.92, and T
=
0.02. As it can be clearly seen in Figure 6(a), the density of states presents two maxima separated by a dip, which plays the role of pseudogap in this scenario. Its presence is due to the warping in the dispersion along the X
→
Y direction (see Figure 3), which induces the presence of the two maxima (one due to the van-Hove singularity at X and one due to the maximum in the dispersion close to S—see Figure 6(b)), and to the loss of states, within this window in frequency, in the region in momentum close to M, as discussed in detail in the previous sections. As a measure of how much weight is lost because of the finite value of the imaginary part of the self-energy in the region in momentum close to M, one can look at the striking difference between the value of A
(
S
_
,
ω
=
0
) in comparison to that of A
(
S
,
ω
=
0
); see Figure 6(b). On reducing the doping, there is an evident transfer of spectral weight between the two maxima; in particular, the weight is transferred from the top of the dispersion close to S to the antinodal point X, where the van Hove singularity resides. At the lowest doping (n
=
0.92), a well-developed pseudogap is visible below the chemical potential and will clearly affect all measurable properties of the system. For this doping, we do not observe any divergence of Σ
c
c
′
(
k
,
ω
=
0
) in contrast to what is reported in [192] where this feature is presented as the ultimate reason for the opening of the pseudogap. In our scenario, the pseudogap is just the result of the transfer of weight from the single-particle density of states to the two-particle one related to the (antiferro) magnetic excitations developing in the system on decreasing doping at low temperatures (see Section 3.6).It is worth mentioning that an analogous doping behavior of the pseudogap has been found by theDMFT
+
Σ approach [163–167], a Mori-like approach by Plakida and coworkers [7, 9], and the cluster perturbation theory [86, 87, 175, 176].


### 3.4. Momentum Distribution Function

To analyze a possible crossover from a Fermi liquid to a non-Fermi-liquid behavior in certain regions of momentum space, at small dopings and low temperatures, and to better characterize the pocket forming on the Fermi surface at the lowest reported doping, we study the electronic momentum distribution functionn
(
k
) per spin along the principal directions (Γ
→
M, M
→
X, X
→
Y and Y
→
Γ) and report it in Figure 7, for U
=
8 and n
=
0.7, 0.75, 0.8, and 0.85 (T
=
0.01), and for n
=
0.92 (T
=
0.02). Figure 7(b) reports a zoom along the main diagonal (Γ
→
M). At the highest studied doping (n
=
0.7), n
(
k
) shows the usual features of a quasinoninteracting system, except for one single, but very important feature: the dip along the main diagonal (Γ
→
M) signaling the presence of a shadow band due to the weak, but anyway finite, antiferromagnetic correlations (see Section 3.6), as already discussed many times hereinbefore. The quite small height of thesecondary jumps along S
→
M and M
→
X directions (with respect to the height of themain jumps along Γ
→
S and Y
→
Γ directions) gives a measure of the relevance of this feature in the overall picture: not really much relevant, except at n
=
0.92 where it changes qualitatively. On increasing filling (reducing doping), the features related to the shadow band do not change much their positions and intensity, while the features related to the ordinary band change their positions as expected in order toaccommodate (i.e., to activate states in momentum for) the increasing number of particles. Finally, at n
=
0.92, the two sets of features close a pocket. At this final stage, what is very relevant, as it is very unconventional, is the evident and remarkable difference in behavior between themain jump along the Γ
→
S direction and thesecondary jump along the S
→
M direction (see Figure 7(b)): the former stays quite sharp (just a bit skewed by the slightly higher temperature) as the Fermi liquid theory requires, the latter instead loses completely its sharpness, much more than what would be reasonable because of the finite value of the temperature as it can be deduced by the comparison with the behavior of themain jump. Such a strong qualitative modification is the evidence of a non-Fermi-liquid-like kind of behavior, but only combining this occurrence with a detailed study of the frequency and temperature dependence of the imaginary part of the self-energy in the very same region of momentum space (see Section 3.5), we will be able to make a definitive statement about this.(Color online) Momentum distribution functionn
(
k
) for U
=
8 and n
=
0.7, 0.75, 0.8, n
=
0.85 (T
=
0.01), and n
=
0.92 (T
=
0.02): (a) along the principal directions (Γ
=
(
0
,
0
)
→
M
=
(
π
,
π
), M
→
X
=
(
π
,
0
), X
→
Y
=
(
0
,
π
), and Y
→
Γ); (b) along the principal diagonal (Γ
→
M).
(a)
(b)In Figure8, we try to summarize the scenario in the extreme case (n
=
0.92, T
=
0.01, and U
=
8), where all the anomalous features are present and well formed, by reporting the full 2D scan of the momentum distribution function n
(
k
) in a quarter of the Brillouin zone. We can clearly see now the pocket with its center along the main diagonal and the lower border touching the border of the magnetic zone at exactly S
=
(
π
/
2
,
π
/
2
). Actually, the 2D prospective makes more evident that there is a second underlying Fermi surface that corresponds to the ordinary paramagnetic one for this filling n
=
0.92 (large and hole-like) and touching the border of the zone between M
=
(
π
,
π
) and X
=
(
π
,
0
) (Y
=
(
0
,
π
)). This corresponds to the very small jump in Figure 9(a) along the same direction. It is worth mentioning that a similar behavior of the momentum distribution function has been found by means of a Mori-like approach by Plakida and coworkers [7, 9].Figure 8
(Color online) The momentum distribution functionn
(
k
) for n
=
0.92, T
=
0.01, and U
=
8.(Color online) The momentum distribution functionn
(
k
) along the main directions (Γ
=
(
0,0
)
→
M
=
(
π
,
π
)
→
X
=
(
π
,
0
)
→
Y
=
(
0
,
π
)
→
Γ
=
(
0,0
)) for different values of temperature T (a) and on-site Coulomb repulsion U (b) at n
=
0.92.
(a)
(b)In Figure9, we study the dependence of the momentum distribution function n
(
k
) on the temperature T (a) and the on-site Coulomb repulsion U (b) by keeping the filling n fixed at the most interesting value: 0.92. At high temperatures (in particular, down to T
=
0.4), the behavior of n
(
k
) is that of a weakly correlated paramagnet (no pocket along the S
→
M direction, no warping along the X
→
Y direction). For lower temperatures, the pocket develops along the S
→
M direction and the signal along the X
→
Y direction is no more constant, signaling the dynamical generation of a diagonal hopping term t
′, connecting same-spin sites in a newly developed antiferromagnetic background unwilling to be disturbed. In fact, such a behavior is what one expects when quite strong magnetic fluctuations develop in the system and corresponds to a well-defined tendency towards an antiferromagnetic phase (see Section 3.6). The M point becomes another minimum in the dispersion in competition with Γ and the dispersion should feature a maximum between them in correspondence to the center of the pocket. The whole bending of the dispersion confines the van Hove singularity below the Fermi surface in anopen pocket (it closes out of the actually chosen Brillouin zone) visible in the momentum distribution as a new quite broad maximum at X and Y. The dependence on U shows that, for n
=
0.92 and T
=
0.01, our solution presents quite strong antiferromagnetic fluctuations for every finite value of the Coulomb repulsion, although the two kinds of pockets discussed just above are not well formed for values of U less than U
=
3
÷
4.


### 3.5. Self-Energy

To obtain clear-cut pieces of information about the lifetime of the quasiparticles generated by the very strong interactions within this scenario, but even more to understand if it is still reasonable or not to discuss in terms of quasiparticles at all (i.e., if this Fermi-liquid-like concept still holds for each region of the momentum and frequency space), we analyze the imaginary part of the self-energyΣ
′
′
k
,
ω as a function of both frequency and temperature. In the top panels of Figure 10, we plot the imaginary part of the self-energy Σ
′
′
k
_
,
ω, together with its real part Σ
′
k
_
,
ω, the spectral function A
(
k
_
,
ω
), and the noninteracting dispersion ɛ
0
k
_
,
ω as functions of the frequency at the nodal point k
_
=
S (a) and at itscompanion position on thephantom half of the pocket k
_
=
S
_ (b) along the main diagonal Γ
→
M. In both cases, although k
_
=
S
_ is not visible in the picture but is very clear for k
_
=
S, the position of the relative/local maximum of A
(
k
_
,
ω
) coincides with the chemical potential (i.e., on the Fermi surface) and it is determined by the sum of ɛ
0
k
_
,
ω
=
0 and Σ
′
k
_
,
ω
=
0 as expected (i.e., both points belong to the r
(
k
)
=
0 locus). What is very interesting and somewhat unexpected and peculiar is that at the nodal point aparabolic-like (i.e., a Fermi-liquid-like) behavior of Σ
′
′
k
_
,
ω is clearly apparent, whereas, at k
_
=
S
_, the dependence of Σ
′
′
k
_
,
ω on frequency shows a predominance of a linear term giving a definite proof that this region in momentum space is interested to a kind of physics very different from what can be considered even by far Fermi-liquid-like.(Color online) Spectral density functionA
k
_
,
ω, real (
′
) and imaginary (
′
′
) parts of the self-energy Σ
k
_
,
ω, noninteracting dispersion ɛ
0
k
_
,
ω as functions of frequency at (a) k
_
=
S and (b) k
_
=
S
_ for n
=
0.92, T
=
0.02, and U
=
8. ((c) and (d)) Imaginary part of the self-energy Σ
′′
k
_
,
ω
=
0 as a function of temperature at (squares) k
=
S and (circles) k
=
S
_ for n
=
0.92 and U
=
8. The blue line is just a guide to the eye.
(a)
(b)
(c)
(d)In order to remove any possible doubt regarding the non-Fermi-liquid-like nature of the physics going on at thephantom half of the pocket, in the bottom panels of Figure 10, the imaginary part of the self-energy at the Fermi surface Σ
′
′
k
_
,
ω
=
0 is reported as a function of the temperature at the nodal point k
_
=
S and at itscompanion
k
_
=
S
_. The blue straight line in the right panel is just a guide to the eye. We clearly see that, in this case too, although it requires to move from a T to a T
2 representation (from (a) to (b)), the behavior of Σ
′
′
k
_
,
ω
=
0 shows rather different behaviors at the selected points. In particular, the temperature dependence of Σ
′
′
k
_
,
ω
=
0 is exactly parabolic (i.e., exactly Fermi-liquid) at the nodal point, while it exhibits a predominance of linear and logarithmic contributions at S
_. This is one of the most relevant results of this analysis and characterizes this scenario with respect to the others present in the literature.


### 3.6. Spin Dynamics

Finally, to analyze the way the system approaches the antiferromagnetic phase on decreasing doping and temperature and increasing correlation strength, we report the behavior of the nearest neighbor spin-spin correlation functionn
z
n
z
α, of the antiferromagnetic correlation length ξ and of the pole ω
(
3
)
(
Q
) (as defined in Section 3) as functions of filling n, temperature T, and on-site Coulomb repulsion U. In Figures 11, 12, and 13, respectively, such quantities are presented for values in the ranges 0.7
<
n
<
0.92, 0.01
<
T
<
1 and 0.1
<
U
<
8. The choice for the extremal values of low doping n
=
0.92, low temperature T
=
0.01, and strong on-site Coulomb repulsion U
=
8 has been made as for these values we find that all investigated single-particle properties (spectral density function, Fermi surface, dispersion relation, density of states, momentum distribution function, self-energy,
…) present anomalous behaviors. The nearest-neighbor spin-spin correlation function n
z
n
z
α is always antiferromagnetic in character (i.e., negative) and increases its absolute value on decreasing doping and temperature T and on increasing U as expected. The signature of the exchange scale of energy J
≈
4
t
2
/
U
≈
0.5 in the temperature dependence as a significative enhancement in the slope is rather evident. The analysis of the filling dependence unveils quite strong correlations at the higher value of doping too: ξ is always larger than one for all values of fillings showing that, at T
=
0.01 and U
=
8, we should expect antiferromagnetic fluctuations in the overdoped regime too as also claimed by recent experiments [193–195]. In the overdoped region, the antiferromagnetic fluctuations are quite less well defined, in terms of magnon/paramagnon width, than in the underdoped region [196–198] and than what is found in the current two-pole approximation. In fact, a proper description of the paramagnons dynamics would definitely require the inclusion of a proper self-energy term in the charge and spin propagators too (in preparation). Coming back to results (Figure 12), ξ overcomes one lattice constant at temperatures below J and tends to diverge for low enough temperatures. On the other hand, ξ seems to saturate for low enough values of doping. ξ equals one between U
=
3 and U
=
4 and again rapidly increases for large enough values of U. The pole ω
(
3
)
(
Q
) decreases on decreasing doping n and temperature T and on increasing U. In particular, it is very sensitive to the variations in temperature T and in on-site Coulomb repulsion U, which make the mode softer and softer clearly showing the definite tendency towards an antiferromagnetic instability.Figure 11
(Color online) The spin-spin correlation functionn
z
n
z
α as a function of filling n, temperature T, and on-site Coulomb repulsion U in the ranges 0.7
<
n
<
0.92, 0.01
<
T
<
1, and 0.1
<
U
<
8.Figure 12
(Color online) The antiferromagnetic correlation lengthξ as a function of filling n, temperature T, and on-site Coulomb repulsion U in the ranges 0.7
<
n
<
0.92, 0.01
<
T
<
1, and 0.1
<
U
<
8.Figure 13
(Color online) The poleω
(
3
)
(
Q
) of the spin-spin propagator as a function of filling n, temperature T, and on-site Coulomb repulsion U in the ranges 0.7
<
n
<
0.92, 0.01
<
T
<
1, and 0.1
<
U
<
8.


## 3.1. Spectral Function and Dispersion

According to its overall relevance in the whole analysis performed hereinafter, we first discuss the electronic dispersion of the model under analysis or, better, its relic in a strongly correlated system. In general, the dispersion and its more or less anomalous features can be inferred by looking at the maxima of the spectral functionA
(
k
,
ω
). In Figures 2 and 3, the spectral function is shown, in scale of grays (increasing from white to black; red is for above-scale values), along the principal directions (Γ
=
(
0
,
0
)
→
M
=
(
π
,
π
), M
→
X
=
(
π
,
0
), X
→
Y
=
(
0
,
π
) and Y
→
Γ) for U
=
8, T
=
0.01 and (a) n
=
0.70, (b) n
=
0.78, (c) n
=
0.85, and (d) n
=
0.92 (T
=
0.02). In Figure 2, the whole range of frequencies with finite values of A
(
k
,
ω
) is reported, while in Figure 3 a zoom in the proximity of the chemical potential is shown. The light gray lines and uniform areas are labeled with the values of the imaginary part of the self-energy Σ
′
′
(
k
,
ω
) and give one of the most relevant keys to interpret the characteristics of the dispersion. The dark green lines in Figure 3 are just guides to the eye and indicate the direction of the dispersion justbefore the visible kink separating the black and the red areas of the dispersion.(Color online) Spectral functionA
(
k
,
ω
) along the principal directions (Γ
=
(
0
,
0
)
→
M
=
(
π
,
π
), M
→
X
=
(
π
,
0
), X
→
Y
=
(
0
,
π
) and Y
→
Γ) for U
=
8, T
=
0.01 and (a) n
=
0.70, (b) n
=
0.78, (c) n
=
0.85, and (d) n
=
0.92 (T
=
0.02).
(a)
(b)
(c)
(d)(Color online) Spectral functionA
(
k
,
ω
) close to the chemical potential (ω
=
0) along the principal directions (Γ
=
(
0
,
0
)
→
M
=
(
π
,
π
), M
→
X
=
(
π
,
0
), X
→
Y
=
(
0
,
π
) and Y
→
Γ) for U
=
8, T
=
0.01 and (a) n
=
0.70, (b) n
=
0.78, (c) n
=
0.85, and (d) n
=
0.92 (T
=
0.02).
(a)
(b)
(c)
(d)The red areas, as they mark the relative maxima ofA
(
k
,
ω
), can be considered as the best possible estimates for the dispersion. In a non- (or weakly-) interacting system, the dispersion would be a single, continuos, and (quite-) sharp line representing some function ɛ
k being the simple pole of G
c
c
(
k
,
ω
). In this case (for a strongly correlated system), instead, we can clearly see that the dispersion is well-defined (red areas) only in some of the regions it crosses in the (
k
,
ω
) plane: the regions where Σ
′
′
(
k
,
ω
) is zero or almost negligible. In the crossed regions where Σ
′
′
(
k
,
ω
) is instead finite, A
(
k
,
ω
) obviously assumes very low values, which would be extremely difficult (actually almost impossible) to detect by ARPES. Accordingly, ARPES would report only the red areas in the picture. This fact is fundamental to understand and describe the experimental findings regarding the Fermi surface in the underdoped regime, as they will be discussed in the next section, and to reconcile ARPES findings with those of quantum oscillations measurements.The two Hubbard subbands, separated by a gap of the orderU and with a reduced bandwidth of order 4
t, are clearly visible: the lower one (LHB) is partly occupied as it is crossed by the chemical potential and the upper one (UHB) is empty and very far from the chemical potential. The lower subband systematically (for each value of doping) loses significance close to Γ and to M, although this effect is more and more pronounced on reducing doping. In both cases (close to Γ and to M), A
(
k
,
ω
) loses weight as Σ
′
′
(
k
,
ω
) increases its own: this can be easily understood if we recall that both χ
0
k
,
ω and χ
3
k
,
ω have a vanishing pole at Γ (due to hydrodynamics) and that χ
3
k
,
ω is strongly peaked at M due to the strong antiferromagnetic correlations present in the system (see Section 3.6). The upper subband, according to the complementary effect induced by the evident shadow-bands appearance (signaled by the relative maxima of Σ
′
′
(
k
,
ω
) marked by 100 and the dark gray area inside the gap) due again to the strong antiferromagnetic correlations present in the system, displays a well-defined dispersion at M, at least for high enough doping, and practically no dispersion at all close to Γ. For low doping, the great majority of the upper subband weight is simply transferred to the lower subband. The growth, on reducing doping, of the undefined-dispersion regions close to Γ in the lower subband cuts down its already reduced bandwidth of order 4
t to values of the order J
=
4
t
2
/
U
=
0.5
t, as one would expect for the dispersion of few holes in a strong antiferromagnetic background. The shape of the dispersion is too compatible with this scenario: the sequence of minima and maxima is compatible, actually driven, by the doubling of the Brillouin zone induced by the strong antiferromagnetic correlations, as well as the dynamical generation of a t
′ diagonal hopping (absent in the Hamiltonian currently under study) clearly signaled by the more and more, on reducing doping, pronounced warping of the dispersion along the X
→
Y direction (more evident in Figure 3), which would be perfectly flat otherwise (for t
′
=
0).Moving to the zooms (Figure3), they show much more clearly the systematic reduction of the bandwidth on reducing doping, the doubling of the zone, the systematic increase of the warping along X
→
Y on reducing doping, and the extreme flatness of the dispersion at X coming from both Γ and M. This latter feature is in very good agreement with quantum Monte Carlo calculations (see [159] and references therein) as well as with ARPES experiments [190], which report a similar behavior in the overdoped region. Moreover, they show that, in contrast with the scenario for a non- (or weakly-) interacting system, where the doubling of the zone happens with the X
→
Y direction aspivot; here, the doubling is confined to the region close to M. The lower subband is completely filled at half-filling, while for an ordinary Slater antiferromagnet the gap opens at half-filling just on top of the van-Hove singularity along X
→
Y. In addition, the effective, finite value of t
′ makes the dispersion maximum close to S higher than the one present along the M
→
X direction, opening the possibility for the appearance of hole pockets close to S. Finally, it is now evident that the warping of the dispersion along X
→
Y will also induce the presence of two maxima in the density of states: one due to the van-Hove singularities at X and Y and one due to the maximum in the dispersion close to S. How deep is the dip between these two maxima just depends on the number of available well-defined (red) states in momentum present between these two values of frequency. This will determine the appearance of a more or less pronounced pseudogap in the density of states, but let us come back to this after having analyzed the region close to M.As a matter of fact, it is just the absence of spectral weight in the region close toM and, in particular and more surprisingly, at the chemical potential (i.e., on the Fermi surface, in contradiction with the Fermi-liquid picture), the main and more relevant result of this analysis, it will determine almost all interesting and anomalous/unconventional features of the single-particle properties of the model. The scenario emerging from this analysis can be relevant not only for the understanding of the physics of the Hubbard model and for the microscopical description of the cuprate high-T
c superconductors, but also for the drafting of a general microscopic theory of strongly correlated systems. The strong antiferromagnetic correlations (through Σ
′
′
(
k
,
ω
), which mainly follows χ
3
k
,
ω) cause a significative and anomalous/unconventional loss of spectral weight around M, which induces in turn the deconstruction of the Fermi surface (Section 3.2), the emergence of momentum selective non-Fermi-liquid features (Sections 3.4 and 3.5), and the opening of a well-developed (deep) pseudogap in the density of states (Section 3.3).Last, but not least, it is also remarkable 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 [82]. Such a phenomenon clearly signals the coupling of the electrons to a bosonic mode. In this scenario, the mode is clearlymagnetic in nature. The frequency of the kink, with respect to the chemical potential, reduces systematically and quite drastically on reducing doping, following the behavior of the pole of χ
3
k
,
ω (see Section 3.6). Combining the presence of kinks and the strong reduction of spectral weight below them, we see the appearance of waterfalls, in particular along the antinodal (X
→
Γ) direction, as found in some ARPES experiments [82]. Finally, the extension of the flat region in the dispersion around the antinodal points (X and Y), that is, at the van-Hove points, increases systematically on decreasing doping. This clearly signals the transfer of spectral weight from the Fermi surface, which is depleted by the strong antiferromagnetic fluctuations, which are also responsible for the remarkable flatness of the band edge.Before moving to the next section, it is worth noticing that similar results for the single-particle excitation spectrum (flat bands close toX, weight transfer from the LHB to the UHB at M, and splitting of the band close to X
,
…) were obtained within the self-consistent projection operator method [28, 29], the operator projection method [30–32], and within a Mori-like approach by Plakida and coworkers [7, 9].


## 3.2. Spectral Function and Fermi Surface

Focusing on the value of the spectral function at the chemical potential,A
(
k
,
ω
=
0
), we can discuss the closest concept to Fermi surface available in a strongly correlated system. In Figure 4, A
(
k
,
ω
=
0
) is plotted as a function of the momentum k in a quarter of the Brillouin zone for U
=
8 and four different couples of values of temperature and filling: (a) n
=
0.7 and T
=
0.01, (b) n
=
0.78 and T
=
0.01, (c) n
=
0.85 and T
=
0.01, and (d) n
=
0.92 and T
=
0.02. The Fermi surface, in agreement with the interpretation of the ARPES measurements within the sudden approximation, can be defined as the locus in momentum space of the relative maxima of A
(
k
,
ω
=
0
). Such a definition opens up the possibility of not only explaining ARPES measurements, but also going beyond them and their finite instrumental resolution and sensitivity with the aim at filling the gap with other kinds of measurements, in particular quantum oscillations ones, which seems to report results in disagreement, up to dichotomy in some cases, with the scenario depicted by ARPES.(Color online) Spectral function at the chemical potentialA
(
k
,
ω
=
0
) as a function of momentum k for U
=
8, T
=
0.01 and (a) n
=
0.70, (b) n
=
0.78, (c) n
=
0.85, and (d) n
=
0.92 (T
=
0.02).
(a)
(b)
(c)
(d)For each value of the filling reported, we can easily distinguish two walls/arcs; forn
=
0.92, they somewhat join. First, let us focus on the arc with the larger (by far) intensities; given the current sensitivities, this is the only one, among the two, possibly visible to ARPES. At n
=
0.7 (see Figure 4(a)), the 3D perspective allows to better appreciate the difference in the intensities of the relative maxima of A
(
k
,
ω
=
0
) between the region close to the main diagonal (M
→
X), where the signal is weaker, and the regions close to the main axes (Γ
→
X and Γ
→
Y), where the signal is stronger; this behavior has been also reported by ARPES experiments [82, 190] as well as the electron-like nature of the Fermi surface [190]. On decreasing doping, this trend reverses, passing through n
=
0.78, where the intensities almost match, and up to n
=
0.92, where the region in proximity of S is the only one with an appreciable signal. The less-intense (by far) arc, reported in [7] too, is the relic of a shadow band, as can be clearly seen in Figure 3, and, consequently, never changes its curvature, in contrast to what happens to the other arc, which is subject to the crossing of the van Hove singularity (n
≅
0.82) instead. Although the ratio between the maximum values of the intensities at the two arcs never goes below two (see Figure 6(b)), there is an evident decrease of the maximum value of the intensity at the larger-intensity arc on decreasing doping, which clearly signals an overall increase of the intensity and/or of the effectiveness (in terms of capability of affecting the relevant quasi-particles, which are those at the Fermi surface) of the correlations.Moving to a 2D perspective (see Figure5), we can add three ingredients to our discussion that can help us better understand the evolution with doping of the Fermi surface: (i) the n
(
k
)
=
0.5 locus (solid line), that is, the Fermi surface if the system would be noninteracting; (ii) the r
(
k
)
=
0 locus (dashed line), that is, the Fermi surface if the system would be a Fermi liquid or somewhat close to it conceptually; (iii) the values (grey lines and labels) of the imaginary part of the self-energy at the chemical potential Σ
c
c
′
′
(
k
,
ω
=
0
) (notice that T
≠
0). Combining these three ingredients with the positions and intensities of the the relative maxima of A
(
k
,
ω
=
0
), we can try to better understand what these latter signify and to classify the behavior of the system on changing doping. At n
=
0.7 (see Figure 5(a)), the positions of the two arcs are exactly matching r
(
k
)
=
0 lines; this will stay valid at each value of the filling reported, with a fine, but very relevant, distinction at n
=
0.92. This occurrence makes our definition of Fermi surface robust, but also versatile as it permits to go beyond Fermi liquid picture without contradicting this latter. The almost perfect coincidence, for the higher value of doping reported, n
=
0.7, of the n
(
k
)
=
0.5 line with the larger-intensity arc clearly asserts that we are dealing with a very-weakly-interacting Fermi metal. Σ
c
c
′
′
(
k
,
ω
=
0
) is quite large close to M (see Figure 3(a)) and eats up the weight of the second arc, which becomes aghost band more than ashadow one. The antiferromagnetic correlations are definitely finite (see Section 3.6) and, consequently, lead to the doubling of the zone, but not strong enough to affect the behavior of the ordinary quasiparticles safely living at the ordinary Fermi surface. Decreasing the doping, we can witness a first topological transition from a Fermi surface closed around Γ (electron-like, hole-like in cuprates language) to a Fermi surface closed around M (hole-like, electron-like in cuprates language) at n
≅
0.82, where the chemical potential crosses the van Hove singularity (see Figure 3). The chemical potential presents an inflection point at this doping (not shown), which allowed us to determine its value with great accuracy. In proximity of the antinodal points (X and Y), a net discrepancy between the position of the relative maxima of A
(
k
,
ω
=
0
) and the locus n
(
k
)
=
0.5 becomes more and more evident on decreasing doping (see Figures 5(b) and 5(c)). This occurrence does not only allow the topological transition, which is absent for the n
(
k
)
=
0.5 locus that reaches the antidiagonal (X
→
Y) at half-filling in agreement with the Luttinger theorem, but also accounts for the apparent broadening of the relative maxima of A
(
k
,
ω
=
0
) close to the antinodal points (X and Y). The broadening is due to the small, but finite, value of Σ
c
c
′
′
(
k
,
ω
=
0
) in those momentum regions (see Figures 3(b) and 3(c)), signaling the net increase of the correlation strength and, accordingly, the impossibility of describing the system in this regime as a conventional non- (or weakly-) interacting system within a Fermi-liquid scenario or its ordinary extensions for ordered phases. What is really interesting and goes beyond the actual problem under analysis (cuprates: 2D Hubbard model) is the emergence of such features only in well-defined regions in momentum space. This selectiveness in momentum is quite a new feature in condensed matter physics and its understanding and description require quite new theoretical approaches. In this system, almost independently from their effective strength, the correlations play a so fundamental role to come to shape and determinequalitatively the response of the system. Accordingly, any attempt to treat correlations without taking into account the level of entanglement between all degrees of freedom present in the system is bound to fail or at least to miss the most relevant features.(Color online) Spectral function at the chemical potentialA
(
k
,
ω
=
0
) as a function of momentum k for U
=
8, T
=
0.01 and (a) n
=
0.70, (b) n
=
0.78, (c) n
=
0.85, and (d) n
=
0.92 (T
=
0.02). The solid line marks the locus n
(
k
)
=
0.5, the dashed line marks the locus r
(
k
)
=
0, the gray lines are labeled with the values of Σ
c
c
′
′
(
k
,
ω
=
0
), and the dotted line is a guide to the eye and marks the reduced (antiferromagnetic) Brillouin zone.
(a)
(b)
(c)
(d)(Color online) (a) Density of statesN
(
ω
) as a function of frequency ω for U
=
8, (black squares) n
=
0.7 and T
=
0.01, (red circles) n
=
0.78 and T
=
0.01, (blue up triangles) n
=
0.85 and T
=
0.01, and (green down triangles) n
=
0.92 and T
=
0.02. (b) Spectral function in the proximity of the chemical potential A
(
k
_
,
ω
~
0
) at (black squares) k
_
=
S
=
(
π
/
2
,
π
/
2
), (red circles) S
_ (in the text), and (blue triangles) X
=
(
π
,
0
) for U
=
8, n
=
0.92, and T
=
0.02.
(a)
(b)Let us come to the most interesting result. Atn
=
0.92, the relative maxima of A
(
k
,
ω
=
0
) detach also from the r
(
k
)
=
0 locus, at least partially, opening a completely new scenario. r
(
k
)
=
0 defines a pocket (close, but absolutely not identical, read below, to that of an antiferromagnet), while the relative maxima of A
(
k
,
ω
=
0
) feature the very same pocket together with quite broad, but still well-defined, wings closing with one half of the pocket (the most intense one) what can be safely considered the relic of a large Fermi surface. This is the second and most surprising topological transition occurring to the actual Fermi surface: the two arcs, clearly visible for all other values of the filling, join and instead of closing just a pocket, as one would expect on the basis of the conventional theory for an antiferromagnet—here mimed by r
(
k
)
=
0 locus, develop (or keep) a completely independent branch. The actual Fermi surface is neither a pocket nor a large Fermi surface; for a more expressive representation, see Figure 4(d). This very unexpected result can be connected to the dichotomy between those experiments (e.g., ARPES) pointing to a small and those ones (e.g., quantum oscillations) pointing to a large Fermi surface. This result can be understood by looking once more at the dispersion for this value of the filling (see Figure 3(d)): the difference in height, induced by the effective finite value of t
′, which increases with decreasing doping, between the two highest maxima in the dispersion—one close to S and the other along the M
→
X direction—makes the latter to cross the chemical potential for larger values of the doping than the first, but given the significative broadening of the dispersion, even when the center mass of the second leaves the Fermi surface (disappearing from r
(
k
)
=
0 locus), its shoulders are still active and well-identifiable at the chemical potential—that is on the actual Fermi surface.The pocket is too far from beingconventional. It is clearly evident in Figure 5(d) that there are two distinct halves of the pocket: one with very high intensity pinned at S (again the only possibly visible to ARPES) and another with very low intensity (visible only totheoreticians and some quantum oscillations experiments). This is our interpretation for the Fermi arcs reported by many ARPES experiments [82, 93] and unaccountable for any ordinary theory relaying on the Fermi liquid picture, though being modified by the presence of an incipient spin or charge ordering. Obviously, looking only at the Fermi arc (as ARPES is forced to do), the Fermi surface looks ill defined as it does not enclose a definite region of momentum space, but having access also to the other half of the pocket, such a problem is greatly alleviated. The point is that the antiferromagnetic fluctuations are so strong to destroy the coherence of the quasiparticles in that region of momentum space as similarly reported within the DMFT
+
Σ approach [163–167] and a Mori-like approach by Plakida and coworkers [7, 9]. Moreover, we will see that (see Section 3.5) thephantom half is not simply lower in intensity because it belongs to a shadow band depleted by a finite imaginary part of the self-energy (see Figure 5(d)) but that it lives in a region of momentum where the imaginary part of the self-energy shows clear signs of non-Fermi-liquid behavior. The lack of next-nearest hopping terms (i.e., t
′ and t
′
′) in the chosen Hamiltonian (1), although they are evidently generated dynamically, does not allow us to perform a quantitative comparison between our results and the experimental ones, which refers to a specific material characterized by a specific set of hoppings. This also explains why our Fermi arc is pinned at S and occupies the outer reduced Brillouin zone, while many experimental results that report a Fermi arc occupying the inner reduced Brillouin zone. Actually, the pinning (with respect to doping within the underdoped region) of thecenter of mass of the ARPES-visible Fermi arc has been reported also from ARPES experiments [191].


## 3.3. Density of States and Pseudogap

The other main issue in the underdoped regime of cuprates superconductors is the presence of a quite strong depletion in the electronic density of states, known as pseudogap. In Figure6(a), we report the density of states 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 in proximity of the chemical potential. As a reference, we also report, in Figure 6(b), the spectral function in proximity of the chemical potential A
(
k
_
,
ω
~
0
) at k
_
=
S
=
(
π
/
2
,
π
/
2
), k
_
=
S
_ which lies where thephantom half of the pocket touches the main diagonal Γ
→
M (i.e., where the dispersion cuts the main diagonal Γ
→
M closer to M), and k
_
=
X
=
(
π
,
0
) for U
=
8, n
=
0.92, and T
=
0.02. As it can be clearly seen in Figure 6(a), the density of states presents two maxima separated by a dip, which plays the role of pseudogap in this scenario. Its presence is due to the warping in the dispersion along the X
→
Y direction (see Figure 3), which induces the presence of the two maxima (one due to the van-Hove singularity at X and one due to the maximum in the dispersion close to S—see Figure 6(b)), and to the loss of states, within this window in frequency, in the region in momentum close to M, as discussed in detail in the previous sections. As a measure of how much weight is lost because of the finite value of the imaginary part of the self-energy in the region in momentum close to M, one can look at the striking difference between the value of A
(
S
_
,
ω
=
0
) in comparison to that of A
(
S
,
ω
=
0
); see Figure 6(b). On reducing the doping, there is an evident transfer of spectral weight between the two maxima; in particular, the weight is transferred from the top of the dispersion close to S to the antinodal point X, where the van Hove singularity resides. At the lowest doping (n
=
0.92), a well-developed pseudogap is visible below the chemical potential and will clearly affect all measurable properties of the system. For this doping, we do not observe any divergence of Σ
c
c
′
(
k
,
ω
=
0
) in contrast to what is reported in [192] where this feature is presented as the ultimate reason for the opening of the pseudogap. In our scenario, the pseudogap is just the result of the transfer of weight from the single-particle density of states to the two-particle one related to the (antiferro) magnetic excitations developing in the system on decreasing doping at low temperatures (see Section 3.6).It is worth mentioning that an analogous doping behavior of the pseudogap has been found by theDMFT
+
Σ approach [163–167], a Mori-like approach by Plakida and coworkers [7, 9], and the cluster perturbation theory [86, 87, 175, 176].


## 3.4. Momentum Distribution Function

To analyze a possible crossover from a Fermi liquid to a non-Fermi-liquid behavior in certain regions of momentum space, at small dopings and low temperatures, and to better characterize the pocket forming on the Fermi surface at the lowest reported doping, we study the electronic momentum distribution functionn
(
k
) per spin along the principal directions (Γ
→
M, M
→
X, X
→
Y and Y
→
Γ) and report it in Figure 7, for U
=
8 and n
=
0.7, 0.75, 0.8, and 0.85 (T
=
0.01), and for n
=
0.92 (T
=
0.02). Figure 7(b) reports a zoom along the main diagonal (Γ
→
M). At the highest studied doping (n
=
0.7), n
(
k
) shows the usual features of a quasinoninteracting system, except for one single, but very important feature: the dip along the main diagonal (Γ
→
M) signaling the presence of a shadow band due to the weak, but anyway finite, antiferromagnetic correlations (see Section 3.6), as already discussed many times hereinbefore. The quite small height of thesecondary jumps along S
→
M and M
→
X directions (with respect to the height of themain jumps along Γ
→
S and Y
→
Γ directions) gives a measure of the relevance of this feature in the overall picture: not really much relevant, except at n
=
0.92 where it changes qualitatively. On increasing filling (reducing doping), the features related to the shadow band do not change much their positions and intensity, while the features related to the ordinary band change their positions as expected in order toaccommodate (i.e., to activate states in momentum for) the increasing number of particles. Finally, at n
=
0.92, the two sets of features close a pocket. At this final stage, what is very relevant, as it is very unconventional, is the evident and remarkable difference in behavior between themain jump along the Γ
→
S direction and thesecondary jump along the S
→
M direction (see Figure 7(b)): the former stays quite sharp (just a bit skewed by the slightly higher temperature) as the Fermi liquid theory requires, the latter instead loses completely its sharpness, much more than what would be reasonable because of the finite value of the temperature as it can be deduced by the comparison with the behavior of themain jump. Such a strong qualitative modification is the evidence of a non-Fermi-liquid-like kind of behavior, but only combining this occurrence with a detailed study of the frequency and temperature dependence of the imaginary part of the self-energy in the very same region of momentum space (see Section 3.5), we will be able to make a definitive statement about this.(Color online) Momentum distribution functionn
(
k
) for U
=
8 and n
=
0.7, 0.75, 0.8, n
=
0.85 (T
=
0.01), and n
=
0.92 (T
=
0.02): (a) along the principal directions (Γ
=
(
0
,
0
)
→
M
=
(
π
,
π
), M
→
X
=
(
π
,
0
), X
→
Y
=
(
0
,
π
), and Y
→
Γ); (b) along the principal diagonal (Γ
→
M).
(a)
(b)In Figure8, we try to summarize the scenario in the extreme case (n
=
0.92, T
=
0.01, and U
=
8), where all the anomalous features are present and well formed, by reporting the full 2D scan of the momentum distribution function n
(
k
) in a quarter of the Brillouin zone. We can clearly see now the pocket with its center along the main diagonal and the lower border touching the border of the magnetic zone at exactly S
=
(
π
/
2
,
π
/
2
). Actually, the 2D prospective makes more evident that there is a second underlying Fermi surface that corresponds to the ordinary paramagnetic one for this filling n
=
0.92 (large and hole-like) and touching the border of the zone between M
=
(
π
,
π
) and X
=
(
π
,
0
) (Y
=
(
0
,
π
)). This corresponds to the very small jump in Figure 9(a) along the same direction. It is worth mentioning that a similar behavior of the momentum distribution function has been found by means of a Mori-like approach by Plakida and coworkers [7, 9].Figure 8
(Color online) The momentum distribution functionn
(
k
) for n
=
0.92, T
=
0.01, and U
=
8.(Color online) The momentum distribution functionn
(
k
) along the main directions (Γ
=
(
0,0
)
→
M
=
(
π
,
π
)
→
X
=
(
π
,
0
)
→
Y
=
(
0
,
π
)
→
Γ
=
(
0,0
)) for different values of temperature T (a) and on-site Coulomb repulsion U (b) at n
=
0.92.
(a)
(b)In Figure9, we study the dependence of the momentum distribution function n
(
k
) on the temperature T (a) and the on-site Coulomb repulsion U (b) by keeping the filling n fixed at the most interesting value: 0.92. At high temperatures (in particular, down to T
=
0.4), the behavior of n
(
k
) is that of a weakly correlated paramagnet (no pocket along the S
→
M direction, no warping along the X
→
Y direction). For lower temperatures, the pocket develops along the S
→
M direction and the signal along the X
→
Y direction is no more constant, signaling the dynamical generation of a diagonal hopping term t
′, connecting same-spin sites in a newly developed antiferromagnetic background unwilling to be disturbed. In fact, such a behavior is what one expects when quite strong magnetic fluctuations develop in the system and corresponds to a well-defined tendency towards an antiferromagnetic phase (see Section 3.6). The M point becomes another minimum in the dispersion in competition with Γ and the dispersion should feature a maximum between them in correspondence to the center of the pocket. The whole bending of the dispersion confines the van Hove singularity below the Fermi surface in anopen pocket (it closes out of the actually chosen Brillouin zone) visible in the momentum distribution as a new quite broad maximum at X and Y. The dependence on U shows that, for n
=
0.92 and T
=
0.01, our solution presents quite strong antiferromagnetic fluctuations for every finite value of the Coulomb repulsion, although the two kinds of pockets discussed just above are not well formed for values of U less than U
=
3
÷
4.


## 3.5. Self-Energy

To obtain clear-cut pieces of information about the lifetime of the quasiparticles generated by the very strong interactions within this scenario, but even more to understand if it is still reasonable or not to discuss in terms of quasiparticles at all (i.e., if this Fermi-liquid-like concept still holds for each region of the momentum and frequency space), we analyze the imaginary part of the self-energyΣ
′
′
k
,
ω as a function of both frequency and temperature. In the top panels of Figure 10, we plot the imaginary part of the self-energy Σ
′
′
k
_
,
ω, together with its real part Σ
′
k
_
,
ω, the spectral function A
(
k
_
,
ω
), and the noninteracting dispersion ɛ
0
k
_
,
ω as functions of the frequency at the nodal point k
_
=
S (a) and at itscompanion position on thephantom half of the pocket k
_
=
S
_ (b) along the main diagonal Γ
→
M. In both cases, although k
_
=
S
_ is not visible in the picture but is very clear for k
_
=
S, the position of the relative/local maximum of A
(
k
_
,
ω
) coincides with the chemical potential (i.e., on the Fermi surface) and it is determined by the sum of ɛ
0
k
_
,
ω
=
0 and Σ
′
k
_
,
ω
=
0 as expected (i.e., both points belong to the r
(
k
)
=
0 locus). What is very interesting and somewhat unexpected and peculiar is that at the nodal point aparabolic-like (i.e., a Fermi-liquid-like) behavior of Σ
′
′
k
_
,
ω is clearly apparent, whereas, at k
_
=
S
_, the dependence of Σ
′
′
k
_
,
ω on frequency shows a predominance of a linear term giving a definite proof that this region in momentum space is interested to a kind of physics very different from what can be considered even by far Fermi-liquid-like.(Color online) Spectral density functionA
k
_
,
ω, real (
′
) and imaginary (
′
′
) parts of the self-energy Σ
k
_
,
ω, noninteracting dispersion ɛ
0
k
_
,
ω as functions of frequency at (a) k
_
=
S and (b) k
_
=
S
_ for n
=
0.92, T
=
0.02, and U
=
8. ((c) and (d)) Imaginary part of the self-energy Σ
′′
k
_
,
ω
=
0 as a function of temperature at (squares) k
=
S and (circles) k
=
S
_ for n
=
0.92 and U
=
8. The blue line is just a guide to the eye.
(a)
(b)
(c)
(d)In order to remove any possible doubt regarding the non-Fermi-liquid-like nature of the physics going on at thephantom half of the pocket, in the bottom panels of Figure 10, the imaginary part of the self-energy at the Fermi surface Σ
′
′
k
_
,
ω
=
0 is reported as a function of the temperature at the nodal point k
_
=
S and at itscompanion
k
_
=
S
_. The blue straight line in the right panel is just a guide to the eye. We clearly see that, in this case too, although it requires to move from a T to a T
2 representation (from (a) to (b)), the behavior of Σ
′
′
k
_
,
ω
=
0 shows rather different behaviors at the selected points. In particular, the temperature dependence of Σ
′
′
k
_
,
ω
=
0 is exactly parabolic (i.e., exactly Fermi-liquid) at the nodal point, while it exhibits a predominance of linear and logarithmic contributions at S
_. This is one of the most relevant results of this analysis and characterizes this scenario with respect to the others present in the literature.


## 3.6. Spin Dynamics

Finally, to analyze the way the system approaches the antiferromagnetic phase on decreasing doping and temperature and increasing correlation strength, we report the behavior of the nearest neighbor spin-spin correlation functionn
z
n
z
α, of the antiferromagnetic correlation length ξ and of the pole ω
(
3
)
(
Q
) (as defined in Section 3) as functions of filling n, temperature T, and on-site Coulomb repulsion U. In Figures 11, 12, and 13, respectively, such quantities are presented for values in the ranges 0.7
<
n
<
0.92, 0.01
<
T
<
1 and 0.1
<
U
<
8. The choice for the extremal values of low doping n
=
0.92, low temperature T
=
0.01, and strong on-site Coulomb repulsion U
=
8 has been made as for these values we find that all investigated single-particle properties (spectral density function, Fermi surface, dispersion relation, density of states, momentum distribution function, self-energy,
…) present anomalous behaviors. The nearest-neighbor spin-spin correlation function n
z
n
z
α is always antiferromagnetic in character (i.e., negative) and increases its absolute value on decreasing doping and temperature T and on increasing U as expected. The signature of the exchange scale of energy J
≈
4
t
2
/
U
≈
0.5 in the temperature dependence as a significative enhancement in the slope is rather evident. The analysis of the filling dependence unveils quite strong correlations at the higher value of doping too: ξ is always larger than one for all values of fillings showing that, at T
=
0.01 and U
=
8, we should expect antiferromagnetic fluctuations in the overdoped regime too as also claimed by recent experiments [193–195]. In the overdoped region, the antiferromagnetic fluctuations are quite less well defined, in terms of magnon/paramagnon width, than in the underdoped region [196–198] and than what is found in the current two-pole approximation. In fact, a proper description of the paramagnons dynamics would definitely require the inclusion of a proper self-energy term in the charge and spin propagators too (in preparation). Coming back to results (Figure 12), ξ overcomes one lattice constant at temperatures below J and tends to diverge for low enough temperatures. On the other hand, ξ seems to saturate for low enough values of doping. ξ equals one between U
=
3 and U
=
4 and again rapidly increases for large enough values of U. The pole ω
(
3
)
(
Q
) decreases on decreasing doping n and temperature T and on increasing U. In particular, it is very sensitive to the variations in temperature T and in on-site Coulomb repulsion U, which make the mode softer and softer clearly showing the definite tendency towards an antiferromagnetic instability.Figure 11
(Color online) The spin-spin correlation functionn
z
n
z
α as a function of filling n, temperature T, and on-site Coulomb repulsion U in the ranges 0.7
<
n
<
0.92, 0.01
<
T
<
1, and 0.1
<
U
<
8.Figure 12
(Color online) The antiferromagnetic correlation lengthξ as a function of filling n, temperature T, and on-site Coulomb repulsion U in the ranges 0.7
<
n
<
0.92, 0.01
<
T
<
1, and 0.1
<
U
<
8.Figure 13
(Color online) The poleω
(
3
)
(
Q
) of the spin-spin propagator as a function of filling n, temperature T, and on-site Coulomb repulsion U in the ranges 0.7
<
n
<
0.92, 0.01
<
T
<
1, and 0.1
<
U
<
8.


## 4. Conclusions and Perspectives

We have reviewed and systematized the theory and the results for the single-particle and magnetic-response properties microscopically derived for the 2D Hubbard model, as minimal model for high-T
c cuprate superconductors, within the Composite Operator Method with the residual self-energy computed in the noncrossing approximation (NCA).Among the several scenarios proposed for the pseudogap origin [119], COM definitely falls into the AF scenario (the pseudogap is a precursor of the AF long-range order) [199–201] as well as the two-particle self-consistent approach (TPSC) [86, 87, 160], the DMFT
+
Σ approach [163–167], and a Mori-like approach by Plakida and coworkers [6, 7, 9].In the limit of strong on-site Coulomb repulsion and low doping, such results show the emergence of a pseudogap scenario, the deconstruction of the Fermi surface in ill-defined open arcs, and the clear signatures of non-Fermi-liquid features similarly to what has been found by ARPES experiments [82] and not only [84]. In particular, we have shown that a very low-intensity signal develops around M point and moves towards S nodal point on decreasing doping coming to close, together with the ordinary Fermi surface boundary, a pocket in the underdoped region. Whenever the pocket develops, it is just the remarkable difference in the intensity of the signal between the two halves of the pocket to make a Fermi arc apparent. As the doping decreases further, the arc shrinks into a point at S exactly at half-filling making it possible to reconcile the large-small Fermi surface dichotomy once the Fermi surface is defined as the relative maxima of the spectral function and the relic of the ordinary paramagnetic Fermi surface are also taken into account. The pseudogap develops since a region in momentum (and frequency) with a very low-intensity signal (it corresponds to thephantom half of the pocket at the chemical potential and to the shadow band out of it) is present between the van Hove singularity and the quite flat band edge (quite flat after the doubling of the Brillouin zone due to the very strong antiferromagnetic fluctuations). On changing doping, a spectral weight transfer takes place in the density of states between the two maxima corresponding to the van Hove singularity and the band edge, respectively. A crossover between a Fermi liquid and a non-Fermi-liquid can be clearly observed in the momentum distribution function (definitely not featuring a sharp jump on thephantom half of the pocket) and in the imaginary part of the self-energy (featuring linear and logarithmic terms in the frequency and temperature dependence instead of the ordinary parabolic term) on decreasing doping at low temperatures and large interaction strength. This crossover exactly corresponds to the process of deconstruction of the Fermi surface. We also report kinks in the dispersion along nodal and antinodal directions. In order to properly interpret the behavior of the spectral density function and of the momentum distribution function, we have also analyzed the characteristic features in the spin-spin correlation function, the antiferromagnetic correlation length, and the pole of the spin-spin propagator. As expected, on reducing doping or temperature and on increasing U, the correlations become stronger and stronger. The exchange scale of energy J is clearly visible in the temperature dependence of the spin-spin correlation function and drives the overall behavior of the magnetic response. These results also demonstrate that a properly microscopically derived susceptibility can give results practically identical or, at least, very similar to those attainable by means of phenomenological susceptibilities especially tailored to describe experiments. This is even more remarkable since COM brings the benefice of a microscopical determination of the temperature and filling dependencies of the correlation length.Many other issues should be addressed in the next future: to improve the Hamiltonian description by adding and fine-tuning longer-range hopping terms in order to quantitatively and not only qualitatively describe specific materials, to verify the stability and the modifications of this scenario with respect to the inclusion of a residual self-energy in the calculation of the charge and spin propagators closing a fully self-consistent cycle, to investigate the charge and spin responses in the full range of momentum and frequency relevant to these systems searching for hourglasses, thresholds, and all other peculiar features experimentally observed, to investigate the superconducting phase and establish that it is nature and relationship with the anomalous features of thenormal phases, to analyze in detail the transition/crossover between the quasiordinary antiferromagnetic phase at half-filling and the underdoped regime.



---

*Source: 515698-2014-11-10.xml*