Stability of the coexistent superconducting-nematic phase under the presence of intersite interactions

We analyze the effect of intersite-interaction terms on the stability of the coexisting superconucting-nematic phase (SC+N) within the extended Hubbard and $t$-$J$-$U$ models on the square lattice. In order to take into account the correlation effects with a proper precision, we use the approach based on the \textit{diagrammatic expansion of the Gutzwiller wave function} (DE-GWF), which goes beyond the renormalized mean field theory (RMFT) in a systematic manner. As a starting point of our analysis we discuss the stability region of the SC+N phase on the intrasite Coulomb repulsion-hole doping plane for the case of the Hubbard model. Next, we show that the exchange interaction term enhances superconductivity while suppresses the nematicity, whereas the intersite Coulomb repulsion term acts in the opposite manner. The competing character of the SC and N phases interplay is clearly visible throughout the analysis. A universal conclusion is that the nematic phase does not survive within the $t$-$J$-$U$ model with the value of $J$ integral typical for the high-T$_C$ cuprates ($J\approx 0.1$eV). For the sake of completeness, the effect of the correlated hopping term is also analyzed. Thus the present discussion contains all relevant two-site interaction terms which appear in the parametrized one-band model within the second quantization scheme. At the end, the influence of the higher-order terms of the diagrammatic expansion on the rotational symmetry breaking is also shown by comparing the DE-GWF results with those corresponding to the RMFT.


I. INTRODUCTION
The nematic ordering is believed to appear in a number of strongly correlated compounds such as URu 2 Si 2 1 , ironpnictides 2,3 , cuprates 4-6 , Sr 3 Ru 2 O 7 7 , as well as quantum Hall systems 8 . Nematicity is characterized by a spontaneous rotational symmetry breaking of the electronic structure, with the preservation of the translational symmetry imposed by the crystal lattice. This condition excludes positional or magnetic orderings such as those appearing in the cases of spin-density-wave (SDW) or charge-density-wave (CDW) phases. However, it has been argued that in the cuprates the CDW phase may be formed through a precursor state which has a nematic character 6 . In some of the copper based compounds a small anisotropy of the lattice makes it difficult to validate the nematic behavior of the electronic wave function, as the C 4 symmetry of the Cu-O planes is already broken by the crystal structure. Nevertheless, in spite of such a small structural anisotropy, a large anisotropy of various physical properties has been observed 5,[9][10][11] . This fact, together with recent research on LESCO, LNSCO, and LBCO compounds 6,12 indicate, that the anisotropic character of electronic properties of Cu-O planes is not a trivial consequence of the lattice distortions. Instead, it may be caused by an intrinsic susceptibility towards the nematic order appearance and may be due to the inter-electronic interactions.
For the copper-based materials the appearance of superconducting phase can also be ascribed to the interelectronic correlation effects. Therefore, the question of the SC and N phases interplay/coexistence within typical models referring to strongly correlated systems is worth exploring. The mean-field analysis of SC+N appearance for the case of phenomenological model suggests that the two phases compete with each other 15 . Other investigations, going beyond the mean-field approach, included methods limited only to weak or intermediate interactions 13,14,16,17 . The SC+N phase induced solely by strong correlations has been analyzed recently 18 for the case of Hubbard model (with intrasite repulsion only), by using the diagrammatic expansion of the Guwtziller wave function (DE-GWF) approach. The same method has been applied by us to the t-J-U model what has lead to a very good quantitative agreement between theory and experiment for the selected principal properties of the superconducting phase in the cuprates 19,20 . Namely, it has been found that, the presence of both the J term and the possibility of having a small but non-zero number of double occupancies at the same time was indispensable in order to obtain the proper quantitative agreement. One should also note that additional interactions terms which are frequently omitted, may affect the stability of various correlationinduced phases 21,22 .
Here we use the DE-GWF method in order to carry out a detailed analysis of nematic and superconducting phases coexistence/competition in the presence of all significant two-site interaction terms, i.e., the antiferromagnetic exchange, the intersite Coulomb repulsion, and the correlated hopping. To show that the C 4 symmetry breaking presented here is due to interelectronic effects, we focus mainly on the square-lattice structure. However, the influence of the preexistent lattice distortion is also discussed. To show that the higher order terms of the diagrammatic expansion are essential to induce the ten-dency towards the spontaneous C 4 symmetry breaking, we compare the obtained results with those calculated within the RMFT method equivalent to the zeroth order expansion of the GWF 23 .
The structure of the paper is as follows. In the next Section we introduce the t-J-U -V model and the DE-GWF method of its solution. In Sec. III we discuss the resulting phase diagram and related physical properties comprising the regimes of pure-and coexisting-phases stability. Conclusions are contained in Sec. IV.

II. MODEL AND METHOD
The most general form of the Hamiltonian considered here is given beloŵ (1) The first two terms contain the single-particle and the correlated-hopping (∼ K) contributions, respectively, the third term represents the antiferromagnetic exchange interaction, and the last two terms refer to the intraand inter-site Coulomb repulsions, respectively. By ... and ... we denote the summations over the nearestneighbors and next-nearest-neighbors, respectively. For J = K = V ≡ 0 we obtain the Hubbard model which constitutes the reference point of our analysis of the particular interaction terms and their influence on the SC+N phase. With the increasing U → ∞ the model reduces to an extended t-J model. In order to take into account the inter-electronic correlations we use the description based on the Gutzwillertype wave function defined by where |Ψ 0 is the non-correlated wave function (to be defined later) and the correlation operatorP G is provided belowP where λ i,Γ ∈ {λ i∅ , λ i↑ , λ i↓ , λ id } are the variational parameters which correspond to four states of the local basis |∅ i , | ↑ i , | ↓ i , | ↑↓ i at site i, respectively. An important step of the DE-GWF method is the application of the condition 24P where x is yet another variational parameter andd HF i ≡ n HF i↑n HF i↓ ,n HF iσ ≡n iσ − n 0 , with n 0 ≡ Ψ 0 |n iσ |Ψ 0 . One should note that λ Γ parameters are all functions of x which results in only one variational parameter of the wave function. As it has been shown in Refs. 24 and 25, condition (4) leads to rapid convergence of the resulting diagrammatic expansion with the increasing order in the resultant variational parameter x.
Within this approach, the expectation value in the correlated state from any two local operators,ô i andô j , can be expressed in the following form The primmed summation has the restrictions l p = l p , l p = i, j for all p and p .
The averages in the non-correlated state on the righthand side of Eq. (5) can be decomposed by the use of the Wick's theorem applied directly in real space and expressed in terms of the correlation functions P ij ≡ ĉ † iσĉ jσ 0 and S ij ≡ ĉ † i↑ĉ † j↓ 0 . Such a procedure allows us to express the ground state energy Ĥ G ≡ Ψ G |Ĥ|Ψ G / Ψ G |Ψ G as a function of P ij , S ij , n 0 , and x. It has been shown that the desirable convergence can be achieved by taking the first 4-6 terms of the expansion in x appearing in Eq. (5). Here the first 5 terms of the diagrammatic expansion (5) have been taken into account when carrying out the calculations.
The effective Schrödinger equation can be derived from the minimization condition of the ground-state energy functional F ≡ Ĥ G − µ G N G , where µ G and N G are the chemical potential and the total number of particles determined in the state |Ψ G , respectively 26,27 . The explicit form the equation is given beloŵ where the effective single-particle Hamiltonian has the formĤ with the effective parameters It is necessary to introduce the real-space cutoff for the parameters P ij and S ij , which are going to be taken into account while executing explicitly the Wick's decomposition of expansion (5). Here, in order to carry out calculations in a reasonable time the maximum distance has been taken as R 2 max = 5a 2 , where a is the lattice constant.
The self consistent equations for all the parameters S ij and P ij are derived after transforming the effective Hamiltonian (7) to the reciprocal space. The solution of self-consistent equations is concomitant with the minimization over variational parameter x. After calculating P ij , S ij , x, µ G , and P ii = n 0 for a selected set of microscopic parameters (t , K, J, U , V ), we can determine the value of the so-called correlated SC gaps j↓ G , as well as the correlated-hopping averages P G,ij ≡ ĉ † iσĉ jσ G . The d-wave gap symmetry is most widely used for the description of high-T C superconductivity in the cuprates. Here, small corrections to the bare d-wave symmetry appear due to the fact that farther-distance than the nearest averages are included, i.e., those corresponding to atomic sites up to |R ij | 2 ≡ |R i − R j | 2 = 5a 2 . In spite of that, the dominant contribution to the pairing amplitude arises from the nearest-neighbor SC averages: However, in general, when the C 4 symmetry is broken, an s-wave admixture to the d-wave component appears. In such a situation it is convenient to introduce the d-wave and s-wave correlated gap parameters, respectively Also, since for the nematic phase the (1, 0) and (0, 1) directions are not equivalent, the corresponding hopping averages will also differ and the following parameter charactering the nematicity can be introduced in the form: In the pure SC phase ∆ G d = 0, ∆ G s ≡ 0, and δP G ≡ 0, whereas in the coexistent SC+N phase: ∆ G d = 0, ∆ G s = 0, and δP G = 0. For the case of pure nematic phase (without the SC order) one obtains ∆ G d = ∆ G s ≡ 0 and δP G = 0, while for the pure paramagnetic (normal) phase with neither SC nor N we have that ∆ G d = ∆ G s ≡ 0 and δP G ≡ 0. In what follows we study systematically the phase diagram involving all the mentioned phases.

III. RESULTS
In our analysis we have selected the hopping parameters as t = −0.35eV and t = 0.25|t| (unless stated otherwise) which are typical for the copper based compounds. All the energies in the presented results are in the units of nearest-neighbor hopping integral |t|. The calculations correspond to the case of square lattice. However, at the end we also discuss the influence of the lattice distortion towards the orthorhombic structure.
First, we analyze the SC+N phase coexistence in the Hubbard model defined by Hamiltonian (1), i.e., with J = K = V = 0 and for the case of square lattice. These results constitute the reference point for the subsequent analysis focused on the influence of particular two-site terms on the onset of nematicity in the extended models.
In Fig. 1 we plot the phase diagrams on the (U, δ) plane, in which we mark the stability region of the nematic phase coexistent with superconductivity (region labelled by SC+N, with ∆ G s = 0, δP G = 0, and ∆ G d = 0). As one can see, the appearance of the (1, 0) and (0, 1) directions inequivalence, which manifests itself by the nonzero values of δP G [shown in Fig. 1 (c)], is accompanied by comcomitant appearance of the s-wave component of the SC correlated gap [shown in Fig. 1 (a)]. However, the s-wave gap amplitude is two orders of magnitude smaller than that corresponding to the d-wave symmetry. For large values of Coulomb repulsion (U 10) superconductivity wins with the nematic phase in the underdoped regime and appears in the pure d-wave form [region labelled by SC in Fig. 1 (a) 1 (d) we show the correlated gap for the case when the nematic phase is not taken into account. In such a situation only d-wave component of the SC gap appears and its values are significantly larger as compared to that in the SC+N phase [c.f. Figs. 1 (b) and (d)]. This means that the adjustment of the SC phase to the C 4 symmetry breaking leads to the weakening of the d-wave SC, what in turn indicates the competing character of SC and N phases interplay.
In Fig. 2 we analyze the effect of the J-term on the C 4 symmetry breaking for two significantly different values of Hubbard U (U = 11.5 and U = 21.5). As one can see, with increasing J the d-wave superconductivity is enhanced while the nematicity gets reduced substantially. Above the value of J ≈ 0.15 the latter is completely destroyed leaving only the pure SC phase without any s-wave component of the gap. For larger U values [ Figs. (b), (d), (f)] the effect of N phase suppression is even stronger. As a result, the nematicity is already destroyed for the set of parameters for which the proper agreement between theory and experiment has been obtained in Ref. 19 with respect to high-T C superconductivity in the copper-based compounds (t = −0.35eV, t = 0.25|t|, U = 22, J = 0.25|t|). Hence, the latter results are not affected by the s-wave SC gap component appearance which would be destructive for the nodal direction presence observed in the cuprates.  Fig. 3 (e)]. Therefore, in the model with both J-and V -terms included, the competition between N and SC phases is determined by both these factors. As a consequence, the SC+N phase can be sustained for values of J typical for the cuprates (J ≈ 0.3) when sufficiently strong intersite Coulomb integral is considered. In Fig. 4 we show such a situation which represent the t-J-U -V model case. However, here the nematicity appears in the overdoped regime which would be against the experimental findings for the cuprates.
For the sake of completeness we also analyze the influence of electronic-structure details on the SC+N phase stability. Namely, in Figs. 3 (b), (d), and (f) we show the doping dependences of the correlated gap components and the nematicity factor, all as functions of hole doping for selected values of the next-nearest-neighbor hopping integral t . As one can see, with the decreasing t value the stability region of the SC+N phase is narrowed down. However, the lower critical concentration of the nematicity onset is not affected and is close to δ = 0.05 [see Fig. 3 (f)], which is similar to the upper critical concentration for the AF phase appearance observed in experiments on the cuprates. Such result differs from the one obtained recently in Ref. 14, where it was shown that the lower critical concentration for the N phase appearance is moving together with the filling value (tuned by t ) which corresponds to the van Hove singularity. This discrepancy can be caused by the differences of the details of the two approaches. Namely, in the above mentioned work the FLEX+DMFT method has been used in the intermediate correlations regime (U = 4) and at higher temperature βt = 20. The off-diagonal elements of the Coulomb interaction between the nearest neighboring lattice sites i, j , with the corresponding two-site integral K ij ≡ ii|V (r − r )|ij introduce the so-called correlated hopping term which also has been studied by us 21 . In Fig. 5 we show the parameters which characterize the SC+N phase as func-tions of both hole doping δ and the correlated hopping integral K. As one can see the influence is not significant up to the values of K ≈ 1, close to which the nematicity is destroyed and the SC order is being reduced. It should be noted that within the present approach the appearance of the nematic phase is not induced by any straightforward mechanism such as the lattice distortion. Instead, the C 4 symmetry of the electronic wave function is broken spontaneously for high enough values of the Hubbard U . Nevertheless, in Fig. 6 we also provide the results with inclusion of the orthorhombic lattice distortion since often such distortion appears in the cuprates. For simplicity, our analysis is carried our for the case when t = 0 and the lattice structure is changed by tuning the t 0,1 /t 1,0 ratio. One can see that when t 0,1 /t 1,0 = 1, the d-wave gap is decreased mainly in the region of SC+N phase stability and the s-wave gap component changes sign [see Figs. 6 (a) and (b)]. In Figs. 6 (c) and (d) we show how the anisotropy in the hopping integrals affect the anisotropy of the hopping averages in the correlated state P G 0,1 , P G 1,0 . As one can see, in the doping range close to δ ≈ 0.1 even for very small lattice anisotropy (t 0,1 /t 1,0 ≈ 0.95, t 0,1 /t 1,0 ≈ 0.97) we obtain a substantial anisotropy of the hopping averages (P G 0,1 /P G 1,0 ≈ 0.6). Moreover, as shown before even for the case of square lattice, one obtains P G 0,1 /P G 1,0 1 which signals a spontaneous C 4 symmetry breaking [red solid line in 6 (c)] and leads to the SC+N phase.
It is not clear what determines the optimal values of doping which lead to the tendency towards anisotropic character of the electronic properties. Nevertheless, significance of the electronic correlations taken into account by higher order terms of the diagrammatic expansion (5) is evident, since the analyzed result can be obtained only by going beyond the RMFT approach. We show this in Fig. 6 (d), where the comparison of the two methods is provided. Since within the RMFT approach no stability of the SC phase is obtained in the Hubbard model, we compare the two methods limiting to the pure nematic phase only. As one can see, for the case of square lattice (t 0,1 /t 1,0 = 1), no nematic behavior (P G 0,1 /P G 1,0 = 1) is obtained according to the RMFT method, whereas within the DE-GWF approach the anisotropic behavior of the electronic system is sustained. Also, as shown in the inset to Fig. 6 (d), in RMFT we obtain P G 0,1 /P G 1,0 ≈ t 0,1 /t 1,0 in the whole doping range, while the DE-GWF approach leads to a large enhancement of the electronic anisotropy for δ 0.15. , components of the correlated gap, as well as P G 0,1 /P G 1,0 (c) all as functions of hole doping for different values of the lattice distortion rate, t0,1/t0,1. For t0,1/t1,0 < 1 the lattice distortion is introduced which enhances nematicity and suppresses d-wave superconductivity. In (d) we show P G 0,1 /P G 1,0 for the case of pure nematic phase for δ = 0.1 as a function of the lattice distortion rate for the case of DE-GWF and RMFT calculations. The inset shows the doping dependence of P G 0,1 /P G 1,0 for the case of pure nematic phase for the selected value of t0,1/t1,0 = 0.98. The results are for the Hubbard model (J = V = K = 0) with U = 11.5.

IV. CONCLUSIONS
This paper is a continuation of our detailed studies of high-T C SC within an extended t-J (or extended Hubbard) model treated within the diagrammatic expansion of the Gutzwiller wave function (DE-GWF) in two dimensions, that goes beyond the renormalized mean-field theory in a systematic manner [19][20][21][22] . Explicitly, we have analyzed the effect of all the significant intersite interaction terms on the coexistence of superconducting (SC) and nematic (N) phases within that method. As a starting point of our analysis we have determined the stability range of the coexistent phase on the (δ, U ) plane for the case of Hubbard model. The coexistent SC+N phase appears for high enough values of the Coulomb repulsion (U 6) and in a wide doping range. Due to the C 4 symmetry breaking the d-wave pairing amplitude is suppressed and the s-wave component of the SC gap appears in the SC+N phase (cf. Fig. 1). This signals a competing character of the SC and N phases interplay. Moreover, the appearance of the s-wave SC component with the onset of nematicity in addition to the d-wave SC hampers the gapless character of the latter in the nodal direction. Nevertheless, the s-wave amplitude is about two orders of magnitude smaller than that of the d-wave.
For the case of the extended model the competition between SC and N is determined by both the exchange interaction and the intersite Coulomb repulsion terms. Namely, the J-term enhances SC and suppresses nematicity, whereas for the V -term the opposite is true (cf. Figs 2 and 3 (a), (c), (d)). According to our analysis of the t-J-U model, the nematicity survives up to J ≈ 0.15 what means that the SC+N phase is already destroyed for the parameter set, for which a good agreement between theory and experiment has been achieved for the copper-based superconductors 19 . Hence, in such as situation the s-wave gap component is absent (only pure d-wave SC survives) and the nodal direction is well defined. Nevertheless, by adding the V -term to the t-J-U model, one could still sustain the stability of the SC+N phase for the values of J ≈ 0.3, typical for the cuprates.
Our analysis of the effect of electronic structure details on the SC+N phase have shown that there is no influence of the van Hove singularity position on the lower critical doping for the coexistent-phase onset. For all the considered t values the lower critical doping remains almost constant and equal to δ c ≈ 0.05, the value close to that, below which the antiferromagnetic phase appears in the cuprates. This result differs with that presented in Ref. 14, where the FLEX+DMFT method has been used. However, as mentioned earlier, the results obtained within the latter method are limited to small Hubbardmodel U values.
The influence of the correlated hopping term on the SC+N phase is not significant up to the value K ≈ 1, where the nematicity is destroyed and the d-wave gap is suppressed (cf. Fig. 5).
As could be expected, the assumed from the start anisotropy of the lattice induces anisotropy of the electronic properties in the whole doping range. However, a substantial increase of the electronic anisotropy is obtained close to δ ≈ 0.1, both for the case of the coexistent SC+N phase [cf. Fig. 6 (c)] and for the pure nematic phase [cf. inset to Fig. 6 (d)]. Such a result brings into mind the experimental data for the cuprates, where a very small structural anisotropy leads to a large effect for selected physical properties 5,[9][10][11] . The latter result is not reproduced within RMFT method, where we obtain P G 0,1 /P G 1,0 ≈ t 0,1 /t 1,0 in the whole doping range. Moreover, within the RMFT, no spontaneous C 4 symmetry breaking appears for the case of square lattice [cf. Fig.  6 (d)]. This in turn demonstrates that the correlation effects taken, into account in the higher-order of the DE-GWF approach, are responsible for the nematic phase appearance for the square-lattice case.
It would be interesting to investigate if the susceptibility towards the C 4 -symmetry breaking of the electronic system can also induce the orthorhombic crystal distortion of the lattice within the present approach. In order to take into account the subtle interplay between the electronic system and the lattice structure, one would have to calculate the hopping integrals in an ab-initio fashion instead of treating them as model parameters as here. Such an analysis could be carried out by combining the DE-GWF method with the EDABI 28-30 approach. Moreover, such a method could also be used to analyze theoretically the interplay between the unconventional superconductivity and lattice distortion, which is observed in the copper based compounds 31 .
The question of connection between the nematicity and the CDW appearance is quite involved and requires a separate detailed analysis.