High temperature superconductivity in the two-dimensional $t$-$J$ model: Gutzwiller wave function solution

A systematic diagrammatic expansion for Gutzwiller-wave functions (DE-GWF) proposed very recently is used for the description of superconducting (SC) ground state in the two-dimensional square-lattice $t$-$J$ model with the hopping electron amplitudes $t$ (and $t'$) between nearest (and next-nearest) neighbors. On the example of the SC state analysis we provide a detailed comparison of the method results with other approaches. Namely: (i) the truncated DE-GWF method reproduces the variational Monte Carlo (VMC) results; (ii) in the lowest (zeroth) order of the expansion the method can reproduce the analytical results of the standard Gutzwiller approximation (GA), as well as of the recently proposed"grand-canonical Gutzwiller approximation"(GCGA). We obtain important features of the SC state. First, the SC gap at the Fermi surface resembles a $d_{x^2-y^2}$-wave only for optimally- and overdoped system, being diminished in the antinodal regions for the underdoped case in a qualitative agreement with experiment. Corrections to the gap structure are shown to arise from the longer range of the real-space pairing. Second, the nodal Fermi velocity is almost constant as a function of doping and agrees semi-quantitatively with experimental results. Third, we compare the doping dependence of the gap magnitude with experimental data. Fourth, we analyze the $\mathbf{k}$-space properties of the model: Fermi surface topology and effective dispersion. The DE-GWF method opens up new perspectives for studying strongly-correlated systems, as: (i) it works in the thermodynamic limit, (ii) is comparable in accuracy to VMC, and (iii) has numerical complexity comparable to GA (i.e., it provides the results much faster than the VMC approach).


I. INTRODUCTION
The Hubbard and the t-J models of strongly correlated fermions play an eminent role in rationalizing the principal properties of high temperature superconductors (for recent reviews see [1][2][3][4][5] ). The relative role of the particles' correlated motion and the binding provided by the kinetic exchange interaction can be clearly visualized in the effective t-J model, where the effective hopping energy ∼ |t|δ ∼ 0.35 eV (δ ≡ 1 − 2n is the hole doping) is comparable or even lower than the kinetic exchange integral J ≈ 0.12 eV. Simply put, the hopping electron drags behind its exchange-coupled nearest neighbor (n.n.) via empty sites and thus preserves the locally bound configuration in such correlated motion throughout the lattice 6 . In effect, this real-space pairing picture is complementary to the more standard virtual boson (paramagnon) exchange mechanism which involves, explicitly or implicitly, a quasiparticle picture and concomitant with it reciprocal-space language [7][8][9][10] . Unfortunately, no single unifying approach, if possible at all, exists in the literature which would unify the Eliashberg-type and the real-space approaches, out of which a Cooper-pair condensate would emerge as a universal state for arbitrary ratio of the band energy ∼ W to the Coulomb repulsion U. The reason for this exclusive character of the approaches is ascribed to the presence of the Mott-Hubbard phase transition that takes place for W/U ≈ 1 (appearing for the half-filled band case) which also delineates the strong-correlation limit for a doped-Mott metallic state, for W substantially smaller than U. This is the regime, where the t-J model is assumed to be valid, even in the presence of partially-filled oxygen 2p states [10][11][12][13] . The validity of this type of physics is assumed throughout the present paper and a quantitative analysis of selected experimental properties, as well as a comparison with variational Monte-Carlo (VMC) results, is undertaken.
One of the approaches designed to interpolate between the W/U ≫ 1 and W/U 1 limits is the Gutzwiller wave function (GWF) approach 14 . Unfortunately, the method does not allow for an extrapolation to the W/U ≪ 1 limit, at least in the simpler Gutzwiller approximation (GA) 15 . Therefore, different forms of the GA-like approaches, appropriate for the t-J model, have been invented under the name of the Renormalized Mean Field Theory (RMFT) 4,[15][16][17][18][19][20] . The last approach based on the t-J model provides a rationalization of the principal characteristics of high temperature superconductors, including selected properties in a semiquantitative manner, particularly when the so-called statistically consistent Gutz-willer approach (SGA) [20][21][22][23][24][25] is incorporated into RMFT. However, one should also mention that neither GA nor SGA provide a stable superconducting state in the Hubbard model.
Under these circumstances, we have undertaken a project involving a full GWF solution via a systematic Diagrammatic Expansion of the GWF (DE-GWF), which becomes applicable to two-and higher-dimensional systems, for both normal 26 and superconducting 27 states. Previously, this solution has been achieved in one-spatial dimension in an iterative manner 28,29 . Obviously, the DE-GWF should reduce to SGA in the limit of infinite dimensions. In our preceding paper 27 we have presented the first results for the Hubbard model.
Here, a detailed analysis is provided for the t-J model, together with a comparison to experiment, as well as to the VMC and GA results. The limitations of the present approach are also discussed, particularly the inability to describe the pseudogap appearance.
The structure of the paper is as follows. In Sec. II we present the DE-GWF method (cf. also Appendices A and B). In Secs. III and IV (cf. also Appendices C, D, and E) we provide details of the numerical analysis and discuss physical results, respectively. In the latter section, we also compare our results with experiment and relate them to VMC and GA results. Finally, in Sec. V we draw conclusions and overview our approach.

A. t-J Model
We start with the t-J model Hamiltonian 30 on a two-dimensional, infinite square latticê whereν i =ν i↑ +ν i↓ withν iσ ≡n iσ (1 −n iσ ), the first term is the kinetic energy part and the second expresses the kinetic exchange. The spin operator is defined asŜ i = {Ŝ z i ,Ŝ + i ,Ŝ − i } and i,j denotes summation over pairs of n.n. sites (bonds). The parameter c is used to switch on (c = 1) or off (c = 0) the density-density interaction term reproducing the two forms and the translational invariance are not assumed and the analytical results presented in this section are valid for phases with broken symmetries. We study system properties in the thermodynamic limit, i.e., the system size L is infinite. We also neglect the three-site terms 30 .

B. Trial wave function
The principal task within a Gutzwiller-type 14 of approach is the calculation of the expectation value of the starting Hamiltonian with respect to the trial wave function, which is defined as where |Ψ 0 is a single-particle product state (Slater determinant) to be specified later. We define the local Gutzwiller correlator in the atomic basis of the form with variational parameters λ i,Γ ∈ {λ i,∅ , λ i,1↑ , λ i,1↓ , λ i,d }, which describe the occupation probabilities of the four possible local states A particularly useful choice of the parameters λ i,Γ is the one which obeyŝ where the Hartree-Fock operators are defined byd HF i ≡n HF i↑n HF i↓ andn HF iσ ≡n i,σ − n iσ with n iσ = Ψ 0 |n iσ |Ψ 0 . This form ofP 2 i decisively simplifies the calculations by eliminating the so-called 'Hartree bubbles' 26,31 .
For the t-J model, we work with zero double-occupancy, which sets λ i,d = 0 and eliminates x as a variational parameter from the solution procedure. Explicitly, from the conditions in Eqs. (5) and (6) we find λ 2 i,d = 1 + x(1 − n i↑ )(1 − n i↓ ) = 0. Calculating x and inserting to the expressions for λ i,1σ and λ i,∅ gives where n i = n i↑ + n i↓ .

C. Diagrammatic sums
Here we discuss the analytical procedure of calculating the expectation value in detail for the kinetic-energy term and we summarize the results for other terms. We start with expressions for the relevant expectation values of interest via the power series in x, i.e., i,σP i , and we have definedd HF l 1 ,...,l k ≡d HF l 1 · · ·d HF l k witĥ d HF ∅ ≡ 1, whereas the primed sums have the summation restrictions l p = l p ′ , l p = i, j for all p, p ′ .
Expectation values can now be evaluated by means of the Wick's theorem 32 and are carried out in real space. Then, in the resulting diagrammatic expansion, the k-th order terms of Eqs. (10)- (11) correspond to diagrams with one (or two) external vertices on sites i (or i and j) and k internal vertices. These vertices are connected with lines (corresponding to contractions from Wick's theorem), which in the case of the superconducting state with intersite pairing are given by where↑ =↓,↓ =↑. At this point, the application of the linked-cluster theorem 32 yields 26 the analytical result for the kinetic energy term where The diagrammatic sums appearing in Eq. (13) are defined by where and the k-th order sum contributions have the following forms where . . . c 0 indicates that only the connected diagrams are to be kept (see Appendix A for exemplary diagrams and their contributions to diagrammatic sums in the two lowest orders).
The notation (1)[(3)] means that for the index (3) also the term in square brackets needs to be taken into account, e.g. T HF j,σĉj,σd HF l 1 ,...,l k c 0 . In the following expressions we will drop the brackets in the upper indices of diagrammatic sums for the sake of brevity.
The exchange term can be rewritten in the form where the spin-component operators are given by The expectation values of the exchange term components can be expressed as For the expressions of the other components see Appendix B.
The diagrammatic sums appearing in the above expressions are defined by Eq.
and the k-th order sum contributions of the following forms In what follows, we evaluate these diagrammatic sums in particular situations.

D. Spin-isotropic case
The above expressions simplify significantly when a system with translational invariance and spin isotropy is considered. Explicitly, they become where n = n iσ = n jσ , g s = 1 In general, this situation is applicable when no Néel-type antiferromagnetism occurs, as for the spin-singlet paired state the spin isotropy is preserved.
In this situation, and if we additionally disregard the T 33 and I 44 terms, relations valid for isotropic system are obtained reproducing analytically the results of GA 15 . It is interesting to see how big is the difference between the exact expressions, Eqs. This difference is analyzed in Appendix C.
If we consider general phases, and we keep the T 33 term, then the expressions for the expectation values of the hopping and the exchange term become When the 4-line contribution from the diagrammatic sum I 44 ij (in Eq. (42)) is neglected, our method reproduces the GCGA results 33   The order k to which we carry out our expansion, is not the only parameter affecting convergence. Another one is the number of |Ψ 0 lines (defined in Eq. (12)) included when calculating the diagrammatic sums. Its effect on results for the spin-spin correlation function is analyzed in Appendix D.

III. VARIATIONAL PROBLEM
In the previous section we have provided analytical results for the expectation values of all terms appearing in the Hamiltonian (1) with respect to the assumed wave function (6).
These results enable us to calculate the ground state energy W ≡ Ĥ G for a fixed |Ψ 0 . The remaining task is the minimization of this energy (or of the functional F ≡ W −2µ G n G , with n G ≡ n iσ G ) with respect to the wave function |Ψ 0 . This wave function enters into the variational problem via n ≡ n iσ 0 and the lines P l,l ′ and S l,l ′ . In the following we consider only translationally invariant wave functions. Since we study superconducting states, the correlated and non-correlated numbers of particles (n G and n) may differ, and hence it is technically easier to minimize the functional F at a constant chemical potential µ G , and not the ground state energy at a constant number of particles n G .
The remaining variational problem leads (cf. e.g. Refs. 17, 35, and 36) to the effective single-particle Schrödinger-like equation with the self-consistently defined effective single-particle Hamiltonian The effective dispersion relation, the effective gap, and eigenenergies ofĤ eff 0 are defined as The resulting self-consistency loop is shown in Fig. 2. The convergence is achieved when the new |Ψ 0 lines differ from the previous ones by less than the assumed precision value, typically 10 −7 .

IV. RESULTS
The self-consistency loop in Fig. 2   in k space (this corresponds to an infinite system size, L → ∞). The typical accuracy of our solution procedure is equal to 10 −7 . We set |t| = −t as our unit of energy and, unless stated otherwise, and present the results for t ′ = 0.25 and J = 0.3. We consider the two cases with c = 0 and c = 1, but, as their results are very close, we show the c = 0 data only in Figs. 4 and 5a. In several figures we provide also the results of the GCGA (and GA) methods, which were obtained by the simplified zeroth order DE-GWF method (equivalent to GCGA or GA, as discussed in Sec. II E).
We carry out the expansion to the fifth order, which in most cases provides quite accurate results. The lower-order results are also exhibited in selected figures to visualize our method's convergence. To calculate the diagrammatic sums we need to neglect long-range |Ψ 0 -lines in real space. Namely, we take as nonzero only the lines , with 14 neighbors). The same condition applies for S i,j , t eff i,j , and ∆ eff i,j . We also define an additional convergence parameter. Namely, we take into account only those contributions to the diagrammatic sums, in which the total Manhattan distance (i.e., |X|+|Y |) of all lines is smaller than R tot typically set to R tot = 26.
In total, we have the three convergence parameters: (i) order k, (ii) |Ψ 0 cutoff radius R, In Fig. 5 we plot the two gaps: the effective gap at the antinodal point ∆ eff k=(π,0) and the correlated gap ∆ G . The effective gap agrees with the experimental values only after rescaling Therefore, a quantitative agreement with the experimental points in Fig. 5a should not be the goal in describing high-temperature superconductors, as including the pseudogap may change the picture essentially. One should also keep in mind that ∆ eff k=(π,0) depends on J. For lower J values we obtain much better agreement with the experiment (but at the same time, the agreement of the nodal Fermi velocity is then worse). Similarly as in VMC calculations 59 , we observe an exponential decay of the gap with the doping reaching the upper critical concentration δ c ∼ 1/2. We term as SC the phase with ∆ G > 10 −4 , which corresponds to gap values of the order of 0.4 K, below which other effects can destabilize the superconductivity. In our model situation however, we still have a stable superconducting solution even if we increase doping above such defined δ c by 8%. One must note that if the experimentally measured gap is usually determined for temperature T 1K, then the tail of ∆ G (δ) beyond δ c will not be detected as then effectively T > T c . In the inset of Fig. 5a and in the upper panel of Fig. 5b we show also the order-of-expansion dependence of the results. It can be seen that, for most of the doping values, the k-th order results are between the results obtained for order k − 1 and k − 2. Moreover, the difference between the orders diminishes with the increasing order. In the panel composing Fig. 6 we exhibit the doping dependence of the Fermi-surface topology, starting from the effective Hamiltonian (45). We also show results for the state with a spontaneously broken rotational symmetry, i.e., the appearance of the so-called Pomeranchuk phase 21,60,61 . This phase has also been investigated by VMC 39,62 . The drawback of using VMC in such calculations is that the finite-size effects become much more important than for the description of the SC phase (typically 12 × 12 points 62 or 8 × 8 points 39 are included within the quarter of the Brillouin zone, cf. also the discussion in Ref. 39).
Our method does not suffer from those finite-size limitations and therefore, it seems more appropriate for analyzing the Fermi-surface properties. It can be seen from Fig. 6b that the correlated Fermi surface differs essentially from the non-interacting one near half filling.
Namely, if we approach the half-filled case the Fermi surface becomes a line as in a bare Hamiltonian with the n.n. hopping only. This is caused by diminishing of certain effective hopping parameters in the vicinity of the half filling (as shown explicitly in Fig. 8b below).
The doping dependence of the Fermi surface in the Pomeranchuk phase is similar to that obtained in the Hubbard model 26 . The role of the Pomeranchuk instability will not be studied in detail here.  In Fig. 9 we present the diagram types for the kinetic energy (T 11 , T 13 , T 33 ), the potential energy (I 2 , I 4 ), and the "correlated delta" (S 11 , S 13 , S 33 ) diagrammatic sums.
We consider the first two orders (i.e., the diagrams with zero and one internal vertex).
For the paramagnetic phase we would have only the diagrams without dashed lines (and obviously, no correlated delta diagrams). The number of diagrams grows exponentially with the increasing order k, and therefore we determine these diagrams by means of a numerical procedure.
The general form of the resulting diagrammatic sums is obtained as e.g.
In order to perform the summations of diagrams over a lattice, we need to assume as nonzero the |Ψ 0 lines up to some finite distance. In the main text we have taken as nonzero the lines (S i,j ≡ S X,Y with X = (i 1 − j 1 ), Y = (i 2 − j 2 ), P X,Y -analogously) fulfilling If the cutoff distance is defined by X 2 + Y 2 ≤ 2, then they are as follows T 11 (0) = P 10 , T 33 (0) = −P 3 10 − P 10 S 2 10 , T 13 (1) = 2P 3 10 P 11 + 2P 10 P 3 11 + 2P 10 P 11 S 2 10 .
In our numerical procedure, when calculating the diagrammatic sums, we start from the general form (as in Eqs. (A1) -(A3)) and sum over the internal vertices positions (here over l 1 ) making sure that the condition X 2 + Y 2 ≤ 25 is fulfilled for all contributing lines P X,Y , S X,Y . vertices are marked with green (black) circles. The numbers in brackets below diagrams represent their multiplicity (a combinatorial factor).

Appendix B: Exchange term evaluation
The expressions for the components of the exchange term are as follows (with m i ≡ n i↑ − n i↓ ) .  (40)). Explicitly, we plot the following quantities where γ ≡ (1−2n) 2 2(n−1) 2 , by e.g. S 22 (0) we understand the zeroth-order diagrammatic sum, and by (...) we denote other diagrammatic sum terms, (see Eq. (32)). According to the above expressions, a situation in which GA approximates the average accurately corresponds to q = 1. If an average is overestimated (underestimated) by GA, this yields q < 1 (q > 1). It can be seen from Fig. 10 that for the exchange term averages q ≈ 1, and therefore GA works quite well for them. However, for the kinetic energy term averages GA largely overestimates the n.n. average (as also reported in Ref. 18) and underestimates the next n.n. average, especially for an underdoped system. This is the reason behind the large discrepancy of the GA and VMC results in this regime. The ratios q are quite similar in VMC-like and full DE-GWF methods. They are also similar in the PM phase (however, for the next n.n. hopping the ratio q 11 is substantially larger).  Nearly linear behavior of the differences in Fig. 11 suggests that the convergence is exponential (a logarithmic scale is used in Fig. 11). Note also that the higher-order results converge more slowly than the lower-order results, what indicates that to obtain the same accuracy (with respect to the complete |Ψ 0 results with all lines included) in a higher order we need to take into account more lines than in a lower order. Therefore, not only the inclusion of higher-order terms is important to improve accuracy, but also the inclusion of longer range lines.

Appendix E: Details of the VMC-like DE-GWF calculations
We set all parameters of the effective Hamiltonian to zero, except for n.n. pairing ∆ eff 10 and hoppings t eff 10 , t eff 11 , as well as t eff 00 playing the role of effective chemical potential. The n.n. hopping is kept fixed, whereas the other parameters are optimized variationally. In the resulting scheme the effective Hamiltonian contains the same variational parameters as that used in VMC 39 .
We have taken as nonzero the |Ψ 0 lines (S i,j ≡ S 0,(i−j) ≡ S XY with X = (i 1 − j 1 ), Y = (i 2 −j 2 ), P XY -analogously) fulfilling X 2 +Y 2 ≤ 25. In the situation when the number of |Ψ 0 lines does not match the number of effective parameters (t eff i,j and ∆ eff i,j ), the self-consistency loop would not find the true minimum of the energy and a more standard minimization of the energy with respect to ∆ eff 10 , t eff 00 , and t eff 11 is necessary. Namely, we numerically search for a minimum of the system grand canonical potential F by calculating its value for fixed ∆ eff 10 , t eff 00 , and t eff 11 . The flowchart of such calculations is presented in Fig. 12. Explicitly, having fixed effective parameters (step 1 in Fig. 12) we may construct the effective Hamiltonian (step 2), calculate the |Ψ 0 lines (step 3), and having them we can obtain the diagrammatic sums and the potential F (step 4). Finally, we choose the solution with ∆ eff 10 , t eff 00 , and t eff 11 corresponding to the lowest potential F .