Fermi surface evolution of the 2D Hubbard model within a novel four-pole approximation

We present a novel solution of the 2D Hubbard model in the framework of the Composite Operator Method within a four-pole approximation. We adopt a basis of four fields: the two Hubbard operators plus two fields describing the Hubbard transitions dressed by nearest-neighbor spin fluctuations. We include these nonlocal operators because spin fluctuations play a dominant role in strongly correlated electronic systems with respect to other types of nonlocal charge, pair and double-occupancy fluctuations. The approximate solution performs very well once compared with advanced (semi-) numerical methods from the weak-to the strong-coupling regime, being by far less computational-resource demanding. We adopt this solution to study the single-particle properties of the model in the strong coupling regime, where the effects of strong short-range magnetic correlations are more relevant and could be responsible for anomalous features. In particular, we will focus on the characterization of the Fermi surface and of its evolution with doping.We present a novel solution of the 2D Hubbard model in the framework of the Composite Operator Method within a four-pole approximation. We adopt a basis of four fields: the two Hubbard operators plus two fields describing the Hubbard transitions dressed by nearest-neighbor spin fluctuations. We include these nonlocal operators because spin fluctuations play a dominant role in strongly correlated electronic systems with respect to other types of nonlocal charge, pair and double-occupancy fluctuations. The approximate solution performs very well once compared with advanced (semi-) numerical methods from the weak-to the strong-coupling regime, being by far less computational-resource demanding. We adopt this solution to study the single-particle properties of the model in the strong coupling regime, where the effects of strong short-range magnetic correlations are more relevant and could be responsible for anomalous features. In particular, we will focus on the characterization of the Fermi surface and of it...


I. INTRODUCTION
Strongly correlated electron systems (SCES) exhibit a very rich phenomenology due to the interplay between different types of interactions and to the closeness of the related energy scales. [1][2][3][4] For instance, we have metal-insulator transitions (MIT) driven by charge, spin and/or orbital ordering, 5 or solely by the competition of kinetic energy and Hubbard repulsion (Mott-Hubbard insulators). 6 The competition or the coexistence between quite different (quasi-) ordered states often leads to very unconventional phases such as the still-puzzling pseudogap phase of underdoped high-temperature superconductors. [7][8][9][10][11][12][13][14][15][16] In condensed matter theory, the SCES prototype is the Hubbard model 17 that features only the electronic kinetic energy and the local Coulomb repulsion and, nonetheless, manages to capture relevant aspects of the basic phenomenology mentioned above, like the Mott-Slater and Mott-Hubbard MIT, 6 and non-Fermi-liquid normal phases characterized and dominated by shortrange spin and/or charge correlations. 2,10,13 Solving this correlated model is an extremely challenging task and definitely requires non-perturbative methods. 18 In this paper, we use the operatorial approach that studies the dynamics of the operators describing the relevant electronic transitions of the system under analysis. 19,20 Several methods use the operatorial approach: the Hubbard approximations, 17 the Projection Operator Method, 21,22 the works of Mori, 23 Rowe, 24 Roth, 25 the Spectral Density Approach 26 and the Composite Operator Method (COM). 10,[27][28][29] In this work, we adopt the COM that is based on the equations of motion and the Green's function (GF) formalism in the Heisenberg picture reformulated to be used for composite operators (CO). The latter are products of electronic operators describing the new elementary excitations appearing in the system due to the strong correlations.
Another key ingredient of the approach is the Algebra Constraints (AC) that are relations among correlation functions dictated by the non-canonical algebra obeyed by the CO. The AC can be used to determine unknown correlation functions and parameters in a self-consistent way. 10,[27][28][29] We apply the COM to the 2D Hubbard model, designing a novel four-pole (4p) approximation, performing very well with respect to numerical methods from the weak-to the strong-coupling regime, thanks to a careful treatment of nearest-neighbor spin fluctuations. In this work, we focus on the behavior of the Fermi surface on varying the doping, and discuss its deviation from the Fermi-liquid scenario and the appearance of peculiar features due to strong correlations.
The paper is organized as follows: in Sec. II, we apply the COM to the Hubbard model and devise a novel 4p approximation; in Sec. III, we show the accuracy of the related 4p solution comparing its results with numerical data and explore the doping evolution of its Fermi surface for a sizable value of the Coulomb repulsion. Then, in Sec. IV, we draw our conclusions.

II. MODEL AND METHOD
We consider the single-band 2D Hubbard model on a square L-site lattice: where is the electronic field operator in spinorial notation [· stands for the inner (scalar) product in spin space] and Heisenberg picture [i = (i, t i ), being i = r i a Bravais lattice vector and t i the time]; n σ (i) = c † σ (i)c σ (i) is the particle density operator for spin σ at site i and n(i) = σ n σ (i) = c † (i)·c(i) the total particle density operator at site i; U is the on-site Coulomb repulsion and µ the chemical potential; t is the nearest-neighbor hopping integral (energy unit hereafter), where k runs over the first Brillouin zone (BZ) and the lattice parameter has been set to 1. For any operator θ(i), we use the notation θ κ (i) = j κ ij θ(j, t i ).
Following COM prescription, 10,27-29 we adopt the basis Ψ: where the Hubbard operators ξ(i) = (1 − n(i))c(i) and η(i) = n(i)c(i) describe the electronic transitions with filling variation per site 1 → 0 and 2 → 1, respectively; ξ s (i) = 1 2 n k (i)σ k · ξ α (i) and η s (i) = 1 2 n k (i)σ k · η α (i) describe the same transitions dressed by nearest-neighbor spin fluctuations. n k (i) = c † (i) · σ k · c(i) is the spin density operator and σ k (k = 1, 2, 3) the Pauli matrices. The nonlocal basis components ξ s (i) and η s (i) have been included because spin fluctuations are key ingredients to capture the effects of strong correlations and influence the dynamics more substantially than other types of fluctuations. The vectorial current of the basis, The first term represents the projection of J on Ψ and δJ is named residual current. ε(i, j), the energy matrix, is obtained enforcing the constraint δJ(i, t), Ψ † (j, t) = 0, 10,[27][28][29] where . . . stands for the thermal average in the grand-canonical ensemble. Such a constraint assures that δJ contains only the physics orthogonal to the one described by the basis Ψ and leads to the relation ε If we neglect the residual current δJ, we can obtain the thermal retarded GF G(i; j) = R Ψ(i)Ψ † (j) and, in particular, its Fourier transform G(k, ω) within a n-pole approximation 10,27-29 solving its Dyson's equation in the frequency-momentum space: The eigenvalues E (ν) (k) of the energy matrix ε are the poles of G and represent the excitation energy spectrum, i.e. the electronic transition energies of the system. Having chosen an operatorial basis containing four operators, we end up with a G with just four poles: this is the essence of the 4p approximation we present in this paper. σ (ν) mn (k) are the spectral density weights per band: where Ω contains the eigenvectors of ε as columns. The correlation functions of the fields of the basis C mn (i, j) = Ψ m (i)Ψ † n (j) = 1 L k e ik·(r i −r j) C mn (k) can be determined in terms of G using the spectral theorem: , being f F the Fermi function. In Sec. III, we characterize the 4p paramagnetic homogeneous solution. Its I(k) and m(k) matrices, and the AC employed in the self-consistency procedure are given in the supplementary material.

III. RESULTS
In order to assess the quality of the 4p solution, we compare its internal energy per site E with data from advanced numerical methods. In Fig. 1, we show the behavior of E as a function of n for T = 0.25 and for several values of U, comparing our 4p results with Dynamical Cluster Approximation (DCA) and Diagrammatic Monte Carlo (DMC) data. 30 We obtain E as E = 1 L H + µn = 8tC α cc + UD (D is the double occupancy). We find a very satisfactory agreement with numerical data and an improvement with respect to the already good three-pole approximation: 20,29,31 the 4p results agree quantitatively with the numerical data for weak and moderate interactions (U = 2 and 4) at all fillings, and for larger interactions for n < 1. The 4p energy at n = 1 shows a deviation from the DCA one for U ≥ 6, because of the development of long-range (longer than the cluster size) magnetic order in the numerical solution. We point out that our 4p homogenous paramagnetic solution is devised to capture strong short-range spin fluctuations in the paramagnetic state and not, obviously, well-established magnetic orderings. 32,33 The overall outcome of this comparison assesses the very good accuracy of the 4p solution.
Once we have checked the reliability of the 4p solution, we can adopt it to study the evolution of the Fermi surface (FS) upon changing n. will determine the presence of two arcs in the FS. In fact, strong short-range antiferromagnetic (AF) fluctuations induce a bending of the band leading to the formation of a second minimum at M, which is anyway not as deep as the one at Γ. We also notice that the first cut, along Γ → M, has always a much larger spectral density weight with respect to the second one, which tends to vanish for n → 1. In fact, close to half-filling, the AF correlations increase and the related effective electronic self-energy reduces the weight of the second cut. It is worth noticing that the overall behavior is very different from the one expected for an ordinary AF solution because of (i) the inequality of the depth of the two minima and of the weight distribution and (ii) their evolution with doping.
In order to investigate the FS topology, we consider the spectral density function A(k, ω) = − 1 π Im[G cc (k, ω)], where G cc (k, ω) = 2 m,n=1 G mn (k, ω) is the electronic GF. We report A(k, ω = 0) contour plots: the positions of its maxima identify the effective FS as in Angle-Resolved PhotoEmission Spectroscopy experiments (see Fig. 3). Solid lines indicate the locus r(k) = 0, where r(k) is the real part of the G cc (k, ω) denominator. We also compute the momentum distribution function per spin n(k) = ∫ dωf F (ω)A(k, ω): the locus n(k) = 0.5 (dashed lines) determines the FS for an ordinary Fermi liquid (FL). We explore the evolution of the FS for T = 0, U = 6, upon changing n between 0.8 and 0.95. The strong short-range AF correlations force the system to depart from the usual FL-like scenario and lead to the formation of a small FS (SFS) arc centered at M. For n = 0.8, the correlation effects are still not so strong: the large FS (LFS) branch is centered at Γ and is not far from the FL result; the spectral density peak is high and well defined. Increasing n from 0.8 to 0.85, a major qualitative change in the FS topology occurs, upon crossing the van Hove singularity (vHs) for n ≈ 0.81 (not shown), i.e. a Lifshitz topological transition: the LFS concavity is reversed, so that for smaller dopings the LFS and the SFS are concentric and centered at M. The discrepancy with the FL case has increased, and the main A(k, ω = 0) peak is lower in intensity. Such a result is really remarkable as the vHs crossing occurs not at half-filling as expected (i) leading to the emergence of a finite critical doping and (ii) signaling the appearance of dynamically generated values of next-nearest neighbor hoppings typical of an AF background. For larger n, the deviation from the FL behavior is quite dramatic, in particular for n ≥ 0.9. In fact, for small doping, stronger correlations cause a profound reorganization of the spectral density: the LFS does not resemble the FS in the Landau's quasiparticle picture anymore and approaches the SFS, which is reducing its size without changing substantially its shape on raising n. The deviation from the FL regime and the presence of the SFS are due to strong correlations and induced by relevant nearest-neighbor spin fluctuations.

IV. CONCLUSIONS
In this paper, we have introduced a novel four-pole solution for the 2D Hubbard model within the Composite Operator Method: we have chosen a basis given by two Hubbard operators and two operators describing the Hubbard transitions dressed by nearest-neighbor spin fluctuations. This fourpole solution is very accurate once compared with advanced numerical methods, being by far less computational-resource demanding. Therefore, we have adopted this new solution to study the Fermi surface evolution varying the filling in the strong coupling regime (low temperature, small doping and sizable values of U). In such a case, the effects of the spin fluctuations are preeminent and induce anomalous features that could be put in connection with some of those measured in angle-resolved photoemission experiments in underdoped high-temperature superconducting cuprates. Remarkably, we have obtained that the Fermi surface exhibits a clear deviation from the Fermi-liquid behavior and develops a small additional arc centered at (π, π), with a spectral weight strongly doping dependent, as a consequence of strong short-range antiferromagnetic correlations.

SUPPLEMENTARY MATERIAL
In the supplementary material, we provide the explicit expressions for all equations defining, within the framework of the COM, the 4p approximation to the 2D Hubbard model presented in this work.