Superconductivity and intra-unit-cell electronic nematic phase in the three-band model of cuprates

The intra-unit-cell nematic phase is studied within the three-band Emery model of the cuprates with the use of the approach based on the diagrammatic expansion of the Gutzwiller wave function (DE-GWF). According to our analysis the spontaneous $C_4$ symmetry breaking of the electronic wave function, leading to the nematic behavior, can appear due to electron correlations induced mainly by the onsite Coulomb repulsion, even in the absence of the corresponding intersite oxygen-oxygen repulsion term. The latter has been considered as the triggering factor of the nematic state formation in a number of previous studies. Also, we show that, at the transition to the nematic phase electron concentration transfer from $d$- to $p$- orbitals takes place, apart from the usually discussed $p_x/p_y$ polarization. The determined stability regime of the nematic phase appears in the doping range similar to that of the paired phase, showing that both phases have a common origin, even though they compete. Also, we show that in a significant doping range a coexistence region of superconductivity and nematicity appears. The results are discussed in the view of the experimental findings considering the relation between nematicity and pseudogap behavior.


I. INTRODUCTION
A number of broken symmetry states appear in the cuprate high temperature superconductors. One of the key issues is to identify the mechanism of their creation and how they are interrelated. The four-fold (C 4 ) rotational symmetry breaking with the preservation of the translational symmetry, which is observed in the cuprates 1 and titanium-oxypnictides 2 corresponds to the creation of the so-called nematic phase. Due to the structural LTT phase transition in La-based cuprates or the orthorombic distortion in YBCO, it is difficult to validate the intrinsic nematic behavior of the electronic wave function, as the C 4 symmetry is already broken by the crystal lattice distortion. Nevertheless, it has been argued that a significant contribution to the nematicity is distinct from the effects related to the lattice 3,4 . The STM measurements on Bi-2212 and NCCOC seem to show a more direct evidence of electronic nematicity, which is not related to the structure 1,5 . This suggests that one of the generic features of the copper-based compounds may be the intrinsic susceptibility towards the C 4 symmetry breaking of the electronic wave function in the CuO 2 planes.
It has been established that the nematic ordering in the cuprates arises from the differences in electron concentrations at the two oxygen sites within each unit cell of the copper-oxygen plane 1 . An analogus situation takes place in the titanium-based materials 2 . Such a charge shift between the p x /p y orbitals is also reported in the charge-density wave (CDW) phase of the cuprates 6 . Therefore, the connection between the two phases has been discussed 3,7,8 . In particular, it has been suggested, that the nematicity may be understood as a precursor state before the formation of charge ordering, in which both the C 4 and translational symmetries are broken 8 . Also, in some analysis the nematic phase has been related to the appearance of the so-called pseudogap phase 1,4,9 . In fact, a strong thermodynamic evidence for the nematic character of the pseudogap state has been reported recently 4 . However, the question if the C 4 symmetry breaking is the primary cause or a secondary effect of the pseudogap behavior still remains open. Nevertheless, since the pseudogap phenomenon is reported down to T ≈ 0 K and is connected with the C 4 symmetry breaking 1,10 , than both superconductivity and nematicity should appear simultaneously in a significant doping range. Again, it is not clear if the pairing appears inside the nematic domains leading to a coexistent superconducting-nematic phase or a phaseseparation scenario is realized.
The nematicity has been studied theoretically in the single band models, used for the effective description of the Cu-O planes of the cuprates [11][12][13][14][15][16][17] . Due to the intra-unit-cell character of the nematic phase formation, a more realistic description should include explicitly the oxygen degrees of freedom. Therefore, the threeband Emery model has also been applied with respect to the considered symmetry breaking within the mean field approach 18,19 or more sophisticated methods 20,21 . In these considerations the Coulomb repulsion between the oxygen orbitals played an important role leading to the nematic instability [18][19][20] or the spin-fluctuation-driven mechanism has been proposed in the strong coupling limit 21 .
Here, we analyze the C 4 symmetry breaking resulting from the p-orbital polarization, n px = n py , in the three band Emery model, with the values of the microscopic parameters appropriate for the cuprates. To focus purely on the susceptibility towards the nematic instability of the electronic wave function we consider the ideal square lattice situation. Within our approach the nematicity appears as a result of strong inter-electronic correlations, which are taken into account by the higher order terms of the diagrammatic expansion of the Gutzwiller wave function (DE-GWF method). The method has been recently applied to both the single-and three-band descriptions of the paired phase of the cuprates and leads to good agreement with the principal experimental observations [22][23][24] . In contrast to the previous results obtained for the Emery model [18][19][20][21] , we show that the nematic behavior of the electronic wave function can be induced by strong interelectronic correlations, with the dominant role of the onsite Coulomb repulsion at the copper d-orbitals, even without the corresponding intersite oxygen-oxygen term. Such a result is also supported by previous analysis carried out for the single-band Hubbard model [14][15][16] . We also study the interplay between the d-wave pairing and nematic phase. In particular, according to our interpretation the C 4 symmetry breaking and the paired state seem to have the same origin and in a significant doping range superconductivity and nematic phase coexist (SC+N), in spite of the circumstance that the two compete. This last result is discussed in view of the experimental findings considering the relation between the nematicity and the pseudogap behavior 4, 10 .
The paper is organized as follows. In the subsequent Section we present the theoretical model and provide some details of the DE-GWF calculation scheme. Next, the results corresponding to the pure nematic phase are discussed, before turning to the analysis of the nematicity-superconductivity interplay. The conclusions are deferred to the last Section.

III. MODEL AND METHODS
We start from the three-band Emery model in the electron representation of the form whereĉ † ilσ (ĉ ilσ ) are the creation (anihilation) operators of electrons with spin σ at the i-th atomic site and orbital denoted by l ∈ {d, p x , p y }. The notation il, jl means that the summation is carried out only for the interorbital nearest neighbor hoppings (cf. Fig. 1). The notation of the corresponding hopping energies is shown in Fig 1. The p orbitals are located at the oxygen atomic sites which reside in between every two nearest neighbor copper sites (containing the d orbital states) located at the nodes of the square lattice. In such a structure a single unit cell contains one copper and two oxygen atomic sites. The second term of the Hamiltonian defines to the d and p x /p y atomic levels ( px = py ≡ p , d − p ≡ dp ), together with the chemical potential contribution. The interaction parameters U d and U px = U py ≡ U p represent the intrasite Coulomb repulsions between two electrons with opposite spins located on the d and p x /p y orbitals, respectively.
The model represents an effective description of a single copper-oxygen plane of the cuprates. Here, we take the typical values of the hopping energies t dp = 1.13 eV, t pp = 0.49 eV, and the charge-transfer energy dp = 3.2 eV. The commonly used values of the interaction parameter U d (U p ) range between 8 − 10.5 eV (4 − 6 eV), depending on the particular approach [25][26][27] .
To take into account the inter-electronic correlations resulting from the significant onsite Coulomb repulsion at the copper atomic sites, we use the approach based on the so-called diagrammatic expansion of the Gutzwiller wave function (DE-GWF method). The method has been applied to both single-and multi-band models 22,23,28-31 as well as recently to the description of the superconducting phase of the cuprates within the three band Emery model 24 , which is also considered here.
The Gutzwiller-type projected many particle wave function is taken in the form where |Ψ 0 represents the wave function of uncorrelated state. The intra-site intra-orbital projection operator has the following form where λ Γ|il are the variational parameters determining relative weights corresponding to |Γ il , which in turn represent states of the local basis on the atomic sites with the three types of orbitals (l ∈ {d, p x , p y }) The consecutive states represent the empty, singly, and doubly occupied local configurations, respectively. As can be seen, the variational parameters, which tune the local electronic configurations in the resulting wave function, are orbital-dependent. By minimizing the energy of the system over the variational parameters one reduces the number of configurations which correspond to increased interaction energies. The details of the DE-GWF calculation scheme as applied to the d-wave superconducting phase in the three-band Emery model are provided in Ref. 24. It should be noted, that the C 4 symmetry breaking leads to much more involved calculations, since the number of the so-called hopping and pairing lines which determine the |Ψ 0 wave function (cf. Ref. 22 and 24 for the definition of the lines) is twice as large as that corresponding to the situation in which the C 4 symmetry is preserved. In the considered model the nematicity is realized by the p-orbital polarization which means that n px = n py within each unit cell. Thus, the i site index in the variational parameters λ Γ|il can be dropped and we end up with three sets of variational parameters λ Γ|d , λ Γ|px , and λ Γ|py , which correspond to different electronic configurations at the three orbitals appearing in the model. As we have shown in Ref. 24, due to the fact that the Coulomb repulsion at the copper orbitals is the dominant energy in the system, the projection at the oxygen orbitals can be omitted by taking λ Γ|px = λ Γ|py ≡ 1. This assumption is also applied here. However, since the oxygen degrees of freedom are particularly important for the creation of the nematic phase, in the Appendix A we show explicitly that the results are not altered significantly by including the p-orbital projection

IV. RESULTS AND DISCUSSION
A. Intra-unit-cell nematicity In this subsection, we analyze the spontaneous formation of the intra-unit-cell nematicity without the inclusion of superconducting pairing. The effect of the latter is studied in the next subsection. In all the figures hole doping is defined in the following manner: δ = 5 − n px − n py − n d , hence the parent compound corresponds to five electrons on each CuO 2 complex. In the nematic phase, the electronic concentration is shifted between the p x and p y orbitals, which induces the C 4 symmetry breaking. The corresponding nematic order parameter is thus defined in the following manner: η ≡ (n px − n py )/(n px + n py ), and represents the normalized p-orbital polarization. For nonzero η also the d-p hopping expectation values in the (1, 0) and (0, 1) directions differ. The parameter corresponding to the latter effect is defined in the analogous manner γ = (P dpx − P dpy )/(P dpx + P dpy ), The nematic order parameter η together with the hopping asymmetry parameter γ, defined by Eqs. (5) and (6) respectively, as a function of hole doping; (b-e) The electronic concentrations np x , np y , n d , the hopping expectation values P dpx , P dpx , and the double occupancies at the d-orbital d 2 d , all as functions of hole doping. Additionally, the values of all the physical quantities in (b-e) for the case of non-nematic (normal) phase are also provided (in blue). They are: np x = np y ≡ np, P dpx = P dpy ≡ P dp , n d (non-N), d 2 d (non-N); (f) The energy difference between the nematic phase and the non-nematic phase ∆E = EN − Enon−N as well as the corresponding kinetic energy contribution ∆E0, both vs. hole doping.
where P dpx and P dpy are the nearest-neighbor hopping expectation values in the correlated state |Ψ G . They are defined in the following manner: P dpx = ĉ † idσĉ ipxσ G and P dpy = ĉ † idσĉ ipyσ G , where ... G = Ψ G |...|Ψ G / Ψ G |Ψ G . We carry out our analysis for the typical values of the model parameters. If not stated otherwise they are set to: t dp = 1.13 eV, t pp = 0.49 eV, dp = d − p = 3.2 eV, U d = 8 eV, U p = 4.1 eV.
As we show in Fig. 2 the nematic phase appears in a significant hole doping range below δ 0.3, where both η = 0 and γ = 0. Relatively small normalized p-orbital polarization (η) induces significantly larger values of the normalized hopping asymmetry (γ). This can be also seen in (b) and (c) where we show explicitly the values of the electronic concentrations n px , n py and hopping expectation values P dpx ,P dpx . For the sake of comparison, in (b-e) we provide the corresponding results for the non-nematic state with the C 4 symmetry constraint (n px = n py ≡ n p , P dpx = P dpy ≡ P dp ). As seen in (d), apart from the usual p-orbital polarization at the transition to the nematic state, there is also a relatively small electron concentration transfer from d to p orbital. This results in a reduced number of double occupancies at the d orbitals (d 2 d ≡ n id↑nid↓ G ) in the nematic state with respect to the normal, non-nematic state (e). The latter effect decreases the interaction energy resulting from the Coulomb repulsion at the copper atomic sites. However, the interaction energy loss at the transition to the nematic phase is at the expense of the kinetic energy gain. Nevertheless, the overall energy balance is advantageous leading to the nematic behavior of the system. This is explicitly shown in (f) where the energy difference between the nematic and non-nematic states is plotted (∆E = E N − E non−N ) as a function of hole doping. Additionally, the kinetic energy gain is also provided in the Figure (∆E 0 ).
One should note that the only interaction terms included in the model have an onsite character, with the dominant contribution coming from the Coulomb repulsion on the copper atomic sites. Therefore, the appearance of nematicity shown in Fig. 2 suggests that the electron correlations resulting from the high U d value play the dominant role in the spontaneous C 4 symmetry breaking. Such a conclusion is distinct from the analysis presented in Refs. 18-21, where the role of the intersite oxygen-oxygen Coulomb repulsion in creating the nematic phase has been emphasized.
To analyze in detail the influence of the U d -term on the nematic behavior, in Fig. 3 we have plotted the order parameter η on the (U d , δ) plane. As one can see, the intrasite Coulomb repulsion integral has to be large enough to induce the onset of nematic phase, what again indicates the significant role of the onsite electronic correlations in the C 4 symmetry breaking. Furthermore, the region of the nematic phase stability is very similar to that corresponding to the superconducting phase stability determined by us very recently in Ref. 24 (Fig. 12 in that paper) within the same model. Furthermore, by reducing the energy difference between the copper and oxygen atomic levels ( dp = d − p ) one moves the nematic phase stability regime towards larger values of U d [cf. Figs. 3 (a) and (b)]. Again, the same effect is seen for the case of the paired phase 24 . As discussed previously 24 , the lowest-energy excitation for the parent compound (∆E = U d − U p + dp ) should be considered as the factor that determines the strength of the electronic correlations in the model. Therefore, by decreasing dp one has to provide higher values of U d so that to achieve large enough ∆E to induce the nematicity. The similarity between nematic phase and superconducting phase behaviors in this respect points to the common origin of both. In the considered scenario such an origin would be the inter-electronic correlations, with the dominant con- tribution comming from the onsite Coulomb repulsion at the copper sites. This interpretation is also consistent with the determined U p dependence of the nematic order parameter (cf. Fig. 4)

B. Coexistence of superconductivity and nematicity
Since both the nematic ordering and the d-wave superconductivity seem to have the same origin in the considered approach, and they reside at the similar area of the (U d ,δ)-phase diagram (cf. Fig. 3 here and Fig. 12 in Ref. 24), the question of interplay between the two is in place here. Therefore, we have carried out calculations in which both the superconducting pairing and the C 4 symmetry breaking may appear together. As shown in our recent paper 24 , within the three-band description various pairing amplitudes contribute to the resultant superconducting phase. They correspond to the intra-and interorbital pairing between subsequent nearest-neighboring atomic sites. Nevertheless, the dominant contribution results from the pairing between the nearest-neighbor dorbitals due to copper. Therefore, here we focus on the analysis of the latter pairing amplitude and do not show the remaining ones, which are much smaller. Since in the nematic phase the (1, 0) and (0, 1) directions within the Cu-O plane are not equivalent, a mixed d-and s-wave pairing may appear in the coexistent superconductingnematic phase. The correlated d-wave and s-wave gap parameters that are going to be analyzed have the following form where the summations run over the nearest-neighbor d-d orbitals only, ∆ with R ij = R i − R j . Since we are considering a homogeneous situation, the i indices in the expressions for the gap parameters can be dropped (∆ As shown in Fig. 5, the d-and s-wave pairing amplitudes, as well as the nematic order parameter η, all become non-zero in the doping range below ∼ 0.3, which indicates the appearance of the coexistent superconducting-nematic phase (SC+N). For comparison, in Fig. 5 (a) we show the d-d pairing amplitude for the case of pure superconducting d-wave state for which the C 4 symmetry is preserved (in blue). Above δ ≈ 0.3 the d-wave SC gap increases with decreasing doping, however, below δ ≈ 0.3 where the nematicity sets in the ∆ (d) dd becomes very weakly dependent on the doping and is significantly reduced with respect to the gap corresponding to the pure SC state (cf. green and blue lines in Fig. 5a). Additionally, in the region of the SC+N phase stability, a small s-wave contribution to the pairing appears (yellow line in Fig. 5a). It can be concluded from the experimental research that superconductivity and nematicity appear simultaneously in a wide doping range reaching above the optimal doping for T ≈ 0 K in the cuprates 1,4,10 . However, it is not clear if in fact the coexistent superconducting-nematic phase appears in the experiments or a phase separation scenario is realized with a purely nematic non-superconductig domains  residing inside an essentially d-wave superconducting environment.
As we have show, within the present approach the two broken-symmetry states can coexist (cf. Fig. 5). Nevertheless, the suppression of the d-wave pairing amplitude in the SC+N phase should be considered as a signature of competition between the d-wave superconductivity and nematic order, which is also seen in the single-band models of the Cu-O plane 16 . In the latter model, the exchange term between the nearest-neighbor atomic sites works in favor of the superconducting phase, reducing the nematic behavior. In order to analyze if a similar effect can be seen here we have carried out the calculations for the SC+N phase in the three-band model (1) with the inclusion of a similar exchange term between the copper atomic-sitesĤ where J is the exchange integral, the summation is carried out over the nearest-neighbor copper atomic sites, andŜ id are the spin operators at the d orbitals. The estimates for the J value via Raman scattering experiments for the undoped situation varies between 0.1 eV and 0.14 eV, depending on particular compound [32][33][34] , which is consistent with the theoretical estmates 35,36 . The model defined by Eq. (1) supplemented with the term given by (9) constitutes the three-band correspondent of the so-called t-J-U model 23 .
As one can see in Fig. 6 for non-zero values of J the d-wave pairing amplitude is enhanced in the wide range of the hole doping, in contrast to both the s-wave pairing amplitude and the nematic order parameter. Above J ≈ 0.16 eV the latter is completely suppressed and the stability of the pure SC phase is restored. Nevertheless, as shown previously 19 the inter-site oxygen-oxygen Coulomb repulsion ∼ V strengthens the nematic phase. Therefore, one can expect that the V -term can lead to the appearance of the nematic behavior even for J ≈ 0.16 eV. Hence, the final form of the ground state with respect to the C 4 symmetry breaking may result from a subtle interplay between the two factors.

V. CONCLUSIONS
We have shown that the nematic phase can appear in the three-band Emery model in the absence of the intersite oxygen-oxygen Coulomb repulsion, which has been considered as the triggering effect of nematicity in the previous study [18][19][20] . Also, as shown here, at the transition to the nematic phase electron concentration transfer between the d-and p-orbitals takes place in addition to the usually discussed p x /p y polarization (cf. Fig. 2e). Such an effect leads to a decrease of the number of double occupancies on the d orbitals, which is energetically favorable due to the strong onsite Coulomb repulsion at those orbitals.
According to our analysis, a spontaneous C 4 symmetry breaking appears due to inter-electronic correlations, strength of which is determined by the energy value corresponding to the electron transfer from the oxygen p-to the copper d-orbitals (for the parent compound ∆E = U d − U p + dp ). A significant value of ∆E has to be reached to induce the nematic phase -a condition that is met for the model parameters corresponding to the copper-oxides. Such an interpretation is consistent with the fact that by decreasing dp one moves the whole nematic phase stability regime towards higher U d values (cf. Fig. 3). Also, the effect of U d and U p parameters on the nematic phase is of opposite character, due to opposite signs of the two factors when entering the ∆E expression (cf. Fig. 4).
The results analyzed here and those presented in our previous report 24 point to a common origin of both the superconducting and nematic phase (cf. Fig. 3 here and Fig. 12 in Ref. 24). Also, we show that the superconducting and nematic phases may coexist in a significant doping range leading to a suppression of the d-wave pairing amplitude and the appearance of a small s-wave contribution to the pairing. Similarly as in the single band picture, the competition between the d-wave pairing and C 4 symmetry breaking may be tuned by the exchange interaction term and the intersite Coulomb repulsion with the former working in favor of the pairing and the later enhancing nematicity.
One should note that according to experimental research both superconductivity and the pseudogap phase appear in a wide doping range (cf . Fig 6d in Ref. 10) reaching above the optimal doping for T T C (with T C being the superconducting critical temperature). On the other hand, a strong evidence of nematicity in the pseudogap state has been provided quite recently 4 . Such experimental picture could be reconciled with the results presented here, where the superconducting-nematic coexistent phase appears also in a relatively wide doping range for T = 0 K (cf. Fig. 5). Furthermore, the weak doping dependence of the correlated d-wave pairing amplitude inside the SC+N phase shown in Fig. 5a can be related to the experimental result presented in Fig 6e of Ref. 10, where it is reported that the gap slope near the nodal direction, corresponding to the d-wave symmetry of the pairing, shows a similar behavior due to the presence of the pseudogap phase. On the other hand, the scenario with no coexistence region but with purely nematic domains residing inside a d-wave superconducting environment, would be in accord with the measured d-wave symmetry of the gap in the whole doping range where the pairing appears. Moreover, it should be noted that within our approach we do not analyze directly the pseudogap behavior. Therefore, the definite answer to the question of the relation between the nematic phase and the pseudogap behaviour is beyond the scope of this paper.
At the end, one should note that our results correspond to the ground state of the system (T = 0 K). A natu- ral question concerns the effect of finite temperatures on both the pairing and C 4 symmetry breaking considered by us here. Unfortunately, the Gutzwiller wave function is designed to describe the ground state properties and the application of such an approach to non-zero temperatures is not resolved as yet.