# Coexistence of Nematic Order and Superconductivity in the Hubbard Model

###### Abstract

We study the interplay of nematic and superconducting order in the two-dimensional Hubbard model and show that they can coexist, especially when superconductivity is not the energetically dominant phase. Due to a breaking of the symmetry, the coexisting phase inherently contains admixture of the -wave pairing components. As a result, the superconducting gap exhibits very non-standard features including changed nodal directions. Our results also show that in the optimally doped regime the superconducting phase is typically unstable towards developing nematicity (breaking of the symmetry). This has implications for the cuprate high- superconductors, for which in this regime the so-called intertwined orders have recently been observed. Namely, the coexisting phase may be viewed as a precursor to such more involved patterns of symmetry breaking.

###### pacs:

74.20.-z, 74.20.Rp, 74.72.GhElectronic nematic instabilities are observed across several families of correlated electron systems Fradkin et al. (2010) including the cuprates Keimer et al. (2015); Peli et al. (); Fradkin et al. (2015); Vojta (2009), pnictides Chu et al. (2010); Fernandes et al. (2014), Ruthanate SrRuO Borzi et al. (2007), heavy fermionic URuSi Okazaki et al. (2011), as well as for dipolar gases in optical lattices Aikawa et al. (2014). In the nematic (N) phase the discrete lattice rotational symmetry is lowered while the system retains its translational symmetry. Such phases often appear close to superconductivity (SC) in the phase diagrams of these systems Keimer et al. (2015). A natural question therefore arises about the interplay of these two types of symmetry breakings, especially since it is widely assumed that both can have a common cause, a large on-site Coulomb interaction Yokoyama et al. (2013); Bünemann et al. (2012). This question is especially relevant for the high- cuprates, as in the optimal-doping regime both orderings were reported Keimer et al. (2015) (in different parameter ranges), as were also other microscopic orders. The coexistence of a number of them has led to the notion of ‘intertwined orders’ Keimer et al. (2015); Fradkin et al. (2015). Here we demonstrate the coexistence of the SC and N orderings in the minimal model for the description of strongly correlated/high- systems, the single-band Hubbard model.

Superconductivity and nematic order (in the form of an electronic ‘Pomeranchuk instability’ Pomeranchuk (1958); Halboth and Metzner (2000); Yamase and Kohno (2000a, b); Jȩdrak and Spałek (2010)) have frequently been considered as competing phenomena Yamase and Kohno (2000b); Jȩdrak and Spałek (2010); Miyanaga and Yamase (2006); Edegger et al. (2006). Consequently, the coexisting (N+SC) phase was studied usually for anisotropic models, where its description is much easier, as the rotational symmetry is broken already in the starting Hamiltonian Yamase and Kohno (2000b); Miyanaga and Yamase (2006); Yamase and Metzner (2006); Edegger et al. (2006); Okamoto et al. (2010); Su and Maier (2011); Yamase (2015). On the other hand, for isotropic models a N+SC phase was studied only in rather specific situations using (i) the perturbation expansion method Neumayr and Metzner (2003), applicable only in the weak-coupling limit, and (ii) a phenomenological model where interactions leading to both orderings were postulated separately Yamase and Metzner (2007).

The current state-of-the-art methods for strongly correlated electron systems are generally not well-suited to capture the subtle effects of the correlation-induced Fermi surface deformations related to a N phase, let alone those of a N+SC phase. The limitations stem mostly from finite system sizes, which translate to a low momentum space (-space) resolution ^{1}^{1}1cf. e.g. Fig. 2(b) of Ref. Edegger et al. (2006). For example, in Dynamical-Mean-Field Theory Okamoto et al. (2010) and Dynamic Cluster Approximation Okamoto et al. (2010); Su and Maier (2011); Gull et al. (2009) the system size is up to
^{2}^{2}2These methods suffer also from the sign problem, which can significantly restrict the temperature and parameter range and, as a result, can render the analysis of stabilities of considered phases difficult or impossible (cf. the discussion in Refs. Su and Maier (2011); Gull et al. (2009))
and in variational Monte Carlo (VMC) Edegger et al. (2006); Zheng et al. (2015); Xiao-Jun et al. (2015) up to . Consequently, finite-size errors can be significant even for anisotropic models ^{3}^{3}3cf. e.g. Fig. 7 of Ref. Edegger et al. (2006) or Fig. 1 of Ref. Zheng et al. (2015) and may even lead to qualitatively different results depending on the system size Su and Maier (2011).
Hence, it is an entirely open question in correlated-electron theory whether the SC and the N phase are generically competing or tend to stabilize each other Moon and Sachdev (2012); Wang and Liu (2013); Lederer et al. (2015); Yamase and Zeyher (2013).

In this work, we overcome the difficulties in describing a N+SC phase and show that N and SC orderings can coexist in the Hubbard model. To this end, we use a variational method based on Gutzwiller wave functions (GWF) Gutzwiller (1963) combined with an efficient diagrammatic expansion (DE) technique Gebhard (1990); Bünemann et al. (2012), which enables one to evaluate expectation values for GWF without any additional uncontrolled approximations. This DE-GWF method has been applied successfully to study Fermi-surface deformations, -wave superconductivity, and quasiparticle band structures in the Hubbard Bünemann et al. (2012); Kaczmarczyk et al. (2013); Kaczmarczyk (2015); Kaczmarczyk et al. (2015), - Kaczmarczyk et al. (2014), Anderson lattice Wysokiński et al. (2015, ), and multiband Münster and Bünemann () models. The method works in the thermodynamic limit, i.e., with no finite size limitations, which enables us to properly investigate the stability of a N+SC phase. In the single-band Hubbard model the dominant nematic order has a -wave form Bünemann et al. (2012); Halboth and Metzner (2000), meaning that the Fermi surfaces are stretched along one lattice axis and compressed along the other. Combined with a -wave superconducting pairing, symmetry requires that there is an additional induced -wave component of the SC gap with both on-site and long-range contributions (the on-site contribution is often neglected Yamase and Kohno (2000b); Edegger et al. (2006)). This leads to a non-trivial gap structure with features such as shifted nodal points and modified zero-gap regions, see below.

Our starting point is the Hubbard model on a two-dimensional, infinite square lattice, as given by the Hamiltonian

(1) |

where is the two-dimensional site-index, and are the hopping integrals for the nearest and next-nearest neighbors, respectively, is the Coulomb interaction, and is the spin quantum number.

To account for electronic correlations, the strength of which is determined by the ratio of , a ‘Jastrow correlator’ is used , where is a single-particle-product wave function (Slater determinant) to be defined later. We work with the Gutzwiller correlator , in which the local correlators can be expressed as . The parameters control occupancies of the four local states , whereas is related to the on-site pairing component Bünemann et al. (2005). The principal task is the evaluation of the expectation value of the Hamiltonian with respect to the Gutzwiller wave function . This evaluation remains a difficult many-particle problem. It has been shown in Gebhard (1990); Bünemann et al. (2012) that an efficient diagrammatic expansion scheme can be formulated for this purpose if the local correlator is chosen such that it fulfils the condition

(2) |

where is a variational parameter and the Hartree-Fock (HF) operators are defined by

(3) |

Here and we already allow for a breaking of the symmetry which, as mentioned above, leads to a finite on-site pairing ^{4}^{4}4It also leads to non-zero double occupancies, , and makes the analysis of a N+SC phase within the - model Yamase and Kohno (2000b); Edegger et al. (2006) inconsistent with the zero double occupancy condition.. In the following, we use the notation for expectation values with respect to and . With our choice of a correlator which satisfies (2) we eliminate on-site terms (the so-called ‘Hartree bubbles’) from the resulting diagrammatic expansion of expectation values. As a consequence, the results of the DE-GWF method converge rapidly
with an increasing order of the expansion parameter , as was demonstrated for one-dimensional systems Bünemann et al. (2012).
Note that with the condition (2), the parameters , in our Gutzwiller correlator are (for a given )
all determined as a function of , which serves as our only remaining variational parameter.
The DE-GWF method is systematic in the sense that in the zeroth order of the expansion it reproduces Kaczmarczyk et al. (2014); Wysokiński et al. (2015) the non-trivial results of the Gutzwiller approximation whereas, with an increasing order, the exact GWF solution is approached. For two-dimensional systems DE-GWF
gives results in agreement with VMC but with better accuracy Kaczmarczyk et al. (2013, 2014).

Within the DE-GWF method we obtain all expectation values, for example , as a power series in ,

(4) |

The explicit form of for states with a finite onsite pairing is given in Kaczmarczyk et al. (2015). The coefficients depend on the wave function or, more precisely, on the expectation values

(5) |

The intersite expectation values serve as lines in our diagrammatic expansion. The number of lines in the diagrams grows with the order . Instead of terminating the expansion (4) with some finite value of , it turns out to be more accurate to include all diagrams up to a certain maximum number of lines . We further need to introduce a real-space cutoff, i.e., we only include lines up to the maximum distance, here (measured in lattice constants).

In the presence of superconductivity we minimize the functional instead of Kaczmarczyk et al. (2013, 2014), where is the chemical potential. The minimization with respect to leads to the effective single-particle equation (cf. Appendix A of Ref. Schickling et al. (2014) and Ref. Seibold et al. (2008))

(6) |

where the effective Hamiltonian is given as

(7) |

Here we introduced the effective hopping and pairing parameters

(8) |

Let us underline that these parameters contain long-range components, with the same cutoff as for the lines (i.e. up to ). Such long-range components are usually neglected in other methods, but they turn out to be important for a proper description of the nematic phases. The remaining task is the self-consistent solution of Eqs. (6)-(8) in -space, together with the minimization condition (see Refs. Kaczmarczyk (2015); Kaczmarczyk et al. (2015, 2014) for details on the numerical procedures). From the final self-consistent solution we can calculate the effective dispersion and gap as Fourier transforms of and , as well as the ground-state energy labelled , , and for the three considered phases.

We select the value of the nearest-neighbor hopping, , as our unit of energy and choose the typical values of other parameters reflecting the cuprate high- superconductors, and (unless stated otherwise). Physical energies (in Kelvins) are obtained by assuming . The diagrammatic expansion was carried out up to terms with lines in the diagrams.

In Fig. 1 we compare the condensation energies (=energy gain relative to the phase without broken symmetry) of the three studied phases. We find the best conditions for a N+SC phase when as, e.g., for in Fig. 1(a) (marked with dot-dashed lines in this and some of the following figures). When is significantly higher than the coexisting phase becomes unstable, as for . In the regime of larger doping, where , an additional SC ordering on top of the N phase is stable even when the pure SC phase has a significantly lower condensation energy, as it is the case for . In fact, for the optimally-doped case the SC phase is higher in energy than the N (or N+SC) phase by independent of the expansion cutoff (as verified for ). The energy gain from developing additional order on top of the ‘optimal’ one-order phase (SC or N) is similar for and maximally equal to K in the vicinity of the crossing of the SC and N phase energies, as visualized in Fig. 1(b). Figure 1(c) shows how the boundary between the SC and N+SC phases evolves with , which is closely related to the crossing of the N and SC phases [cf. Fig. 1(a)]. This picture would be changed if the spin-exchange term was introduced into the starting Hamiltonian (i.e. in the -- model) to make up for the underestimation of spin-exchange effects by GWF. Such term favors the SC phase and we find that a N+SC phase is stable up to , above which the pure SC phase dominates. On the other hand, the N phase induces a distortion of the underlying lattice of ions Fradkin et al. (2010), which should favor the N and N+SC phases. To account for such effects the inclusion of electron-phonon coupling is required Ohgoe and Imada (2014); Watanabe et al. (2015), which is beyond the scope of the present paper.

In Fig. 2 we elucidate the non-standard gap structure in a N+SC phase. Breaking of the symmetry induces an additional -wave component of the SC gap. In Fig. 2(a) we show the integral of the magnitude of the - and -wave gap components, defined as along the Fermi surface. The ‘-wave contribution’ curve shows the ratio of the -wave integral to the sum of the two integrals to quantify the -wave input to the pairing. It is equal to () for a pure ()-wave state. Strikingly, although the energy gain from developing a SC gap on top of the N phase is rather small (maximally K), the value of the effective gap can be of the same order of magnitude as that in the pure SC phase. The deviation of the gap from the standard behavior along the Fermi surface is shown in Fig. 2(b). Such deviation is present even for the SC phase Kaczmarczyk et al. (2013, 2014), as also observed experimentally Mesot et al. (1999); Kondo et al. (2009); Hashimoto et al. (2014). For a N+SC phase additionally the nodal point is slightly shifted away from the diagonal direction (as marked with an arrow), in contrast to previous results Su and Maier (2011) for an anisotropic model, and the effect increases with doping. The inset shows the polar plots of the gap for N+SC and SC phases. In the former case the Fermi surface is open in the vertical and closed in the horizontal direction. To quantify the nematicity of the system we plot in Fig. 2(c) the dispersion relation along high-symmetry lines for . For the case without nematicity, the dispersions at the and points would be equal. The difference in the values of dispersion at these two points divided by the bandwidth is shown in Fig. 2(d) as a measure of the dispersion asymmetry for N+SC and N phases. The SC order does not modify the dispersion asymmetry significantly (nor the Fermi surface), unless the SC phase has lower energy than the N phase.

In Fig. 3 we show the effective gap in the Brillouin zone obtained from the DE-GWF method (with long-range contributions up to ) for the SC phase in Fig. 3(a) and N+SC phase in Fig. 3(b), as well as the corresponding gap structures with contributions only up to nearest neighbors for both phases in Fig. 3(c)-(d), as usually assumed in other methods (e.g. in VMC). It can be seen from Fig. 3 that such an assumption does not reflect all principal features of the optimal variational solutions. For example, the longer-range components of the gap are mostly opposite to the dominant component and this leads to circles with zero gap around the and points of the Brillouin zone in Fig. 3(a). For a N+SC phase the gap structure is significantly modified with respect to the pure SC phase: (i) the magnitudes of the gap values at and are different; (ii) the zero-gap direction is no longer a straight line along the diagonal but an irregular line along one of the axes (coinciding with the direction, in which the Fermi surface is open); (iii) a larger part of the Brillouin zone contributes to pairing.

The coexistence of the nematic and superconducting orders has been demonstrated in the Hubbard model by using the full Gutzwiller wave function (GWF). Application of the diagrammatic-expansion (DE) technique has enabled us to investigate the properties of the system without finite-size limitations, a condition crucial for the description of the nematic phases. We have shown that the superconducting and nematic orders coexist in the Hubbard model unless the pure superconducting phase is significantly lower in energy than the pure nematic phase. We have obtained the phase diagram, the energies and other properties of the investigated pure and coexisting phases. The gap structure in the coexisting phase is unconventional due to the breaking of the symmetry: the induced -wave gap component shifts the nodal point away from the diagonal direction and modifies the zero-gap region to form in the direction of open Fermi surface.

In the optimal doping regime pure superconductivity turns out not to be the dominating phenomenon as it is unstable against a -wave nematic instability with (as well as without) an additional superconducting order. This observation may be related to the fact that the cuprate high- superconductors develop additional orders in this regime. Namely, the investigated phase can be viewed as a precursor to more complicated orders including stripes and phases with charge density wave order involving more complex patterns of symmetry breaking.

The developed formalism can also be applied to other situations including dipolar Fermi gases in optical lattices where the anisotropy of dipolar interactions leads to appearance of the N phase Aikawa et al. (2014) and superconductivity can be induced by an attractive on-site interaction .

Acknowledgements. The authors are grateful to Florian Gebhard and Mikhail Lemeshko for discussions and critical reading of the manuscript. The work was supported by the Ministry of Science and Higher Education in Poland through the Iuventus Plus grant No. IP2012 017172 for the years 2013-2015, as well as by the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007-2013) under REA grant agreement n [291734]. JK acknowledges hospitality of the Leibniz Universität in Hannover where a large part of the work was performed.

## References

- Fradkin et al. (2010) E. Fradkin, S. A. Kivelson, M. J. Lawler, J. P. Eisenstein, and A. P. Mackenzie, Annu. Rev. Condens. Matter Phys. 1, 153 (2010).
- Keimer et al. (2015) B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, Nature 518, 179 (2015).
- (3) S. Peli et al., arXiv: 1508.03075 .
- Fradkin et al. (2015) E. Fradkin, S. A. Kivelson, and J. M. Tranquada, Rev. Mod. Phys. 87, 457 (2015).
- Vojta (2009) M. Vojta, Adv. Phys. 58, 699 (2009).
- Chu et al. (2010) J.-H. Chu, J. G. Analytis, K. De Greve, P. L. McMahon, Z. Islam, Y. Yamamoto, and I. R. Fisher, Science 329, 824 (2010).
- Fernandes et al. (2014) R. M. Fernandes, A. V. Chubukov, and J. Schmalian, Nat. Phys. 10, 97 (2014).
- Borzi et al. (2007) R. A. Borzi, S. A. Grigera, J. Farrell, R. S. Perry, S. J. S. Lister, S. L. Lee, D. A. Tennant, Y. Maeno, and A. P. Mackenzie, Science 315, 214 (2007).
- Okazaki et al. (2011) R. Okazaki, T. Shibauchi, H. J. Shi, Y. Haga, T. D. Matsuda, E. Yamamoto, Y. Onuki, H. Ikeda, and Y. Matsuda, Science 331, 439 (2011).
- Aikawa et al. (2014) K. Aikawa, S. Baier, A. Frisch, M. Mark, C. Ravensbergen, and F. Ferlaino, Science 345, 1484 (2014).
- Yokoyama et al. (2013) H. Yokoyama, M. Ogata, Y. Tanaka, K. Kobayashi, and H. Tsuchiura, J. Phys. Soc. Jpn. 82, 014707 (2013).
- Bünemann et al. (2012) J. Bünemann, T. Schickling, and F. Gebhard, Europhys. Lett. 98, 27006 (2012).
- Pomeranchuk (1958) I. Pomeranchuk, Sov. Phys. JETP 35, 524 (1958).
- Halboth and Metzner (2000) C. J. Halboth and W. Metzner, Phys. Rev. Lett. 85, 5162 (2000).
- Yamase and Kohno (2000a) H. Yamase and H. Kohno, J. Phys. Soc. Jpn. 69, 332 (2000a).
- Yamase and Kohno (2000b) H. Yamase and H. Kohno, J. Phys. Soc. Jpn. 69, 2151 (2000b).
- Jȩdrak and Spałek (2010) J. Jȩdrak and J. Spałek, Phys. Rev. B 81, 073108 (2010).
- Miyanaga and Yamase (2006) A. Miyanaga and H. Yamase, Phys. Rev. B 73, 174513 (2006).
- Edegger et al. (2006) B. Edegger, V. N. Muthukumar, and C. Gros, Phys. Rev. B 74, 165109 (2006).
- Yamase and Metzner (2006) H. Yamase and W. Metzner, Phys. Rev. B 73, 214517 (2006).
- Okamoto et al. (2010) S. Okamoto, D. Sénéchal, M. Civelli, and A.-M. S. Tremblay, Phys. Rev. B 82, 180511 (2010).
- Su and Maier (2011) S.-Q. Su and T. A. Maier, Phys. Rev. B 84, 220506 (2011).
- Yamase (2015) H. Yamase, Phys. Rev. B 91, 195121 (2015).
- Neumayr and Metzner (2003) A. Neumayr and W. Metzner, Phys. Rev. B 67, 035112 (2003).
- Yamase and Metzner (2007) H. Yamase and W. Metzner, Phys. Rev. B 75, 155117 (2007).
- (26) Cf. e.g. Fig. 2(b) of Ref. Edegger et al. (2006).
- Gull et al. (2009) E. Gull, O. Parcollet, P. Werner, and A. J. Millis, Phys. Rev. B 80, 245102 (2009).
- (28) These methods suffer also from the sign problem, which can significantly restrict the temperature and parameter range and, as a result, can render the analysis of stabilities of considered phases difficult or impossible (cf. the discussion in Refs. Su and Maier (2011); Gull et al. (2009)).
- Zheng et al. (2015) X.-J. Zheng, Z.-B. Huang, D.-Y. Liu, and L.-J. Zou, Phys. Rev. B 92, 085109 (2015).
- Xiao-Jun et al. (2015) Z. Xiao-Jun, H. Zhong-Bing, and Z. Liang-Jian, Chin. Phys. B 24, 017404 (2015).
- (31) Cf. e.g. Fig. 7 of Ref. Edegger et al. (2006) or Fig. 1 of Ref. Zheng et al. (2015).
- Moon and Sachdev (2012) E.-G. Moon and S. Sachdev, Phys. Rev. B 85, 184511 (2012).
- Wang and Liu (2013) J. Wang and G.-Z. Liu, New J. Phys. 15, 073039 (2013).
- Lederer et al. (2015) S. Lederer, Y. Schattner, E. Berg, and S. A. Kivelson, Phys. Rev. Lett. 114, 097001 (2015).
- Yamase and Zeyher (2013) H. Yamase and R. Zeyher, Phys. Rev. B 88, 180502 (2013).
- Gutzwiller (1963) M. C. Gutzwiller, Phys. Rev. Lett. 10, 159 (1963).
- Gebhard (1990) F. Gebhard, Phys. Rev. B 41, 9452 (1990).
- Kaczmarczyk et al. (2013) J. Kaczmarczyk, J. Spałek, T. Schickling, and J. Bünemann, Phys. Rev. B 88, 115127 (2013).
- Kaczmarczyk (2015) J. Kaczmarczyk, Philos. Mag. 95, 563 (2015).
- Kaczmarczyk et al. (2015) J. Kaczmarczyk, T. Schickling, and J. Bünemann, Phys. Status Solidi (b) 252, 2059 (2015).
- Kaczmarczyk et al. (2014) J. Kaczmarczyk, J. Bünemann, and J. Spałek, New J. Phys. 16, 073018 (2014).
- Wysokiński et al. (2015) M. M. Wysokiński, J. Kaczmarczyk, and J. Spałek, Phys. Rev. B 92, 125135 (2015).
- (43) M. M. Wysokiński, J. Kaczmarczyk, and J. Spałek, arXiv: 1510.00224 .
- (44) K. z. Münster and J. Bünemann, arXiv: 1510.07896 .
- Bünemann et al. (2005) J. Bünemann, F. Gebhard, K. RadnóÃ³czi, and P. Fazekas, J. Phys.: Condens. Matter 17, 3807 (2005).
- (46) It also leads to non-zero double occupancies, , and makes the analysis of a N+SC phase within the - model Yamase and Kohno (2000b); Edegger et al. (2006) inconsistent with the zero double occupancy condition.
- Schickling et al. (2014) T. Schickling, J. BÃ¼nemann, F. Gebhard, and W. Weber, New J. Phys. 16, 093034 (2014).
- Seibold et al. (2008) G. Seibold, F. Becca, and J. Lorenzana, Phys. Rev. Lett. 100, 016405 (2008).
- Ohgoe and Imada (2014) T. Ohgoe and M. Imada, Phys. Rev. B 89, 195139 (2014).
- Watanabe et al. (2015) H. Watanabe, K. Seki, and S. Yunoki, Phys. Rev. B 91, 205135 (2015).
- Mesot et al. (1999) J. Mesot et al., Phys. Rev. Lett. 83, 840 (1999).
- Kondo et al. (2009) T. Kondo, R. Khasanov, T. Takeuchi, J. Schmalian, and A. Kaminski, Nature 457, 296 (2009).
- Hashimoto et al. (2014) M. Hashimoto, I. M. Vishik, R.-H. He, T. P. Devereaux, and Z.-X. Shen, Nat. Phys. 10, 483 (2014).