SO(4) FLEX+DMFT formalism with SU(2)$\otimes$SU(2)-symmetric impurity solver for superconductivity in the repulsive Hubbard model

TL;DR

This paper introduces a novel SO(4) symmetric FLEX+DMFT formalism with an SU(2)$_{ m spin} imes$SU(2)$_{ m eta spin}$ impurity solver, enabling simultaneous investigation of superconductivity and antiferromagnetism in the Hubbard model.

Contribution

The paper develops a full-SU(2) slave-boson formalism integrated into FLEX+DMFT, allowing for a unified treatment of competing quantum orders.

Findings

Incorporates SO(4) symmetry into FLEX+DMFT framework.

Uses an SU(2)$_{ m spin} imes$SU(2)$_{ m eta spin}$ impurity solver.

Enables study of coexistence and competition of quantum orders.

Abstract

Here we have developed a FLEX+DMFT formalism, where the symmetry properties of the system are incorporated by constructing a SO(4) generalization of the conventional fluctuation-exchange approximation (FLEX) coupled self-consistently to the dynamical mean-field theory (DMFT). Along with this line, we emphasize that the SO(4) symmetry is the lowest group-symmetry that enables us to investigate superconductivity and antiferromagnetism on an equal footing. We have imposed this by decomposing the electron operator into auxiliary fermionic and slave-boson constituents that respect SU(2)$_{\rm spin}\otimes$SU(2)$_{\eta{\rm spin}}$. This is used not in a mean-field treatment as in the usual slave-boson formalisms, but instead in the DMFT impurity solver with an SU(2)$_{\rm spin}\otimes$SU(2)$_{\eta{\rm spin}}$ hybridization function to incorporate the FLEX-generated bath information into DMFT iterations. While there have been attempts such as the doublon-less SU(2) slave-boson formalism, the present "full-SU(2)" slave-boson formalism is expected to provide a new platform for addressing the underlying physics for various quantum orders, which compete with each other and can coexist.

Figures

Equations

$$\mathrm{SO}(4)\otimes \mathrm{U}(1)$$

$$\mathrm{SO}(5)$$

$${\mathrm{SU}(2)}_{\mathrm{s}}\otimes {\mathrm{SU}(2)}_{\eta }$$

$$H=-$$

$$+$$

$${\mathbf{\eta }}_i=({\eta }_i^x,{\eta }_i^y,{\eta }_i^z)$$

$$=\frac{1}{2}\left({ϵ}_i{c}_{i\downarrow }^{†}{c}_{i\uparrow }^{†}+{ϵ}_i{c}_{i\uparrow }{c}_{i\downarrow },i{ϵ}_i{c}_{i\downarrow }^{†}{c}_{i\uparrow }^{†}-i{ϵ}_i{c}_{i\uparrow }{c}_{i\downarrow },1-n_{i\uparrow }-n_{i\downarrow }\right),$$

$${\eta }_i^{+}={\eta }_i^x+i{\eta }_i^y={ϵ}_i{c}_{i\uparrow }{c}_{i\downarrow },$$

$${\eta }_i^{-}={\eta }_i^x-i{\eta }_i^y={ϵ}_i{c}_{i\downarrow }^{†}{c}_{i\uparrow }^{†},$$

$${ϵ}_i=\{\begin{array}{l}1i\in \mathrm{A},\\ -1i\in \mathrm{B}.\end{array}$$

$$H_{\mathrm{U}}=$$

$$=$$

$$\displaystyle+\frac{U}{2}\big{(}c^{\dagger}_{i\uparrow}c^{\dagger}_{i\downarrow}c_{i\downarrow}c_{i\uparrow}+c_{i\downarrow}c_{i\uparrow}c^{\dagger}_{i\uparrow}c^{\dagger}_{i\downarrow}\big{)}\Big{]}$$

$$=\sum _i\frac{2U}{3}{\mathbf{\eta }}_i^2+\sum _{i}2\mu {\eta }_i^z,$$

$${\tilde{C}}_i=\left(\begin{array}{cc}{c}_{i\uparrow }& {ϵ}_i{c}_{i\downarrow }^{†}\\ c_{i\downarrow }& -{ϵ}_i{c}_{i\uparrow }^{†}\end{array}\right),$$

$$H_{SU2}=$$

$$\displaystyle-\left(\sum\limits_{i}^{\rm A}\sum\limits_{j(\neq i)}^{\rm A}+\sum\limits_{i}^{\rm B}\sum\limits_{j(\neq i)}^{\rm B}\right)v_{|i-j|}{\rm Tr}\big{(}\widetilde{C}^{\dagger}_{i}\widetilde{C}_{j}\big{)}+H_{\rm U},$$

$$\displaystyle{\bm{\eta}}_{i}=\frac{1}{4}{\rm Tr}\big{(}\widetilde{C}_{i}{\bm{\tau}}\widetilde{C}^{\dagger}_{i}\big{)}.$$

$$\displaystyle{\bm{S}}_{i}=\frac{1}{4}{\rm Tr}\big{(}\widetilde{C}^{\dagger}_{i}{\bm{\tau}}\widetilde{C}_{i}\big{)}$$

$$=\frac{1}{2}\left(c_{i\downarrow }^{†}{c}_{i\uparrow }+c_{i\uparrow }^{†}{c}_{i\downarrow },i{c}_{i\downarrow }^{†}{c}_{i\uparrow }-i{ϵ}_i{c}_{i\uparrow }^{†}{c}_{i\downarrow },n_{i\uparrow }-n_{i\downarrow }\right),$$

$$S_i^{+}=S_i^x+i{S}_i^y=c_{i\uparrow }^{†}{c}_{i\downarrow },$$

$$S_i^{-}=S_i^x-i{S}_i^y=c_{i\downarrow }^{†}{c}_{i\uparrow }.$$

$$\mathrm{SU}(2{)}_{\mathrm{s}}\otimes \mathrm{SU}(2{)}_{\eta }/{\mathbb{Z}}_2\approx \mathrm{SO}(4),$$

$$C_{\mathrm{N}i}^{†}=\left(\begin{array}{cccc}{c}_{i\uparrow }^{†},& -{ϵ}_i{c}_{i\uparrow }^{†},& c_{i\downarrow },& {ϵ}_i^{\ast }{c}_{i\downarrow }\end{array}\right),$$

$${\Psi }_{1i}=(0,-{ϵ}_i{c}_{i\uparrow }^{†},c_{i\downarrow },0),$$

$${\Psi }_{2i}=(c_{i\uparrow },0,0,{ϵ}_i{c}_{i\downarrow }^{†}),$$

$$C_{\mathrm{N}i}^{†}={\Psi }_{2i}^{†}+{\Psi }_{1i}.$$

$$C_{\mathrm{N}i}^{†}=\left(\begin{array}{cccc}{c}_{iA\uparrow }^{†},& -{\widetilde{ϵ}}_B{c}_{iB\uparrow }^{†},& c_{iA\downarrow },& {\widetilde{ϵ}}_B^{\ast }{c}_{iB\downarrow }\end{array}\right).$$

$$C_{\mathrm{N}\mathbf{k}}^{†}(\tau )=$$

$$\left(\begin{array}{cccc}{c}_{\mathbf{k}A\uparrow }^{†}(\tau ),& -c_{\mathbf{k}+\mathbf{Q}B\uparrow }^{†}(\tau ),& c_{-\mathbf{k}A\downarrow }(\tau ),& c_{-(\mathbf{k}+\mathbf{Q})B\downarrow }(\tau )\end{array}\right).$$

$$G_{\alpha \beta \sigma \sigma }(\mathbf{k},{\mathbf{k}}^{\prime };\overline{\tau })=-⟨{\mathcal{T}}_{\tau }{c}_{\mathbf{k}\alpha \sigma }(\tau )c_{{\mathbf{k}}^{\prime }\beta \sigma }^{†}({\tau }^{\prime })⟩,$$

$$F_{\alpha \beta \sigma \overline{\sigma }}(\mathbf{k},{\mathbf{k}}^{\prime };\overline{\tau })=-⟨{\mathcal{T}}_{\tau }{c}_{\mathbf{k}\alpha \sigma }(\tau )c_{{\mathbf{k}}^{\prime }\beta \overline{\sigma }}({\tau }^{\prime })⟩.$$

$${\mathcal{G}}_{4\times 4}(\mathbf{k},\overline{\tau })$$

$${\mathcal{G}}_{4\times 4}(\mathbf{k};\overline{\tau })=$$

$$\left(\begin{array}{cccc}{G}_{AA\uparrow \uparrow }(\mathbf{k},\mathbf{k};\overline{\tau })& -G_{AB\uparrow \uparrow }(\mathbf{k},\mathbf{k}+\mathbf{Q};\overline{\tau })& F_{AA\uparrow \downarrow }(\mathbf{k},-\mathbf{k};\overline{\tau })& F_{AB\uparrow \downarrow }(\mathbf{k},-(\mathbf{k}+\mathbf{Q});\overline{\tau })\\ -G_{BA\uparrow \uparrow }(\mathbf{k}+\mathbf{Q},\mathbf{k};\overline{\tau })& G_{BB\uparrow \uparrow }(\mathbf{k}+\mathbf{Q},\mathbf{k}+\mathbf{Q};\overline{\tau })& -F_{BA\uparrow \downarrow }(\mathbf{k}+\mathbf{Q},-\mathbf{k};\overline{\tau })& -F_{BB\uparrow \downarrow }(\mathbf{k}+\mathbf{Q},-(\mathbf{k}+\mathbf{Q});\overline{\tau })\\ F_{AA\downarrow \uparrow }^{\ast }(-\mathbf{k},\mathbf{k};\overline{\tau })& -F_{AB\downarrow \uparrow }^{\ast }(-\mathbf{k},\mathbf{k}+\mathbf{Q};\overline{\tau })& -G_{AA\downarrow \downarrow }(-\mathbf{k},-\mathbf{k};\overline{\tau })& -G_{AB\downarrow \downarrow }(-\mathbf{k},-(\mathbf{k}+\mathbf{Q});\overline{\tau })\\ F_{BA\downarrow \uparrow }^{\ast }(-(\mathbf{k}+\mathbf{Q}),\mathbf{k};\overline{\tau })& -F_{BB\downarrow \uparrow }^{\ast }(-(\mathbf{k}+\mathbf{Q}),\mathbf{k}+\mathbf{Q};\overline{\tau })& -G_{BA\downarrow \downarrow }(-(\mathbf{k}+\mathbf{Q}),-\mathbf{k};\overline{\tau })& -G_{BB\downarrow \downarrow }(-(\mathbf{k}+\mathbf{Q}),-(\mathbf{k}+\mathbf{Q});\overline{\tau })\end{array}\right).$$

$${\mathcal{G}}_{4\times 4}(\mathbf{k};i{\omega }_n)={\int }_0^{\beta }\mathrm{d}\tau e^{i{\omega }_n\tau }{\mathcal{G}}_{4\times 4}(\mathbf{k};\tau ),$$

$${\mathcal{G}}_{4\times 4}(k{)}^{-1}=i{\omega }_n{\mathbb{I}}_{4\times 4}-{\epsilon }_{4\times 4}(\mathbf{k})-{\Sigma }_{4\times 4}^{\mathrm{eff}}(k),$$

$${\epsilon }_{4\times 4}(\mathbf{k})=$$

Full text

I Introduction

It has been recognized that several high-temperature superconductor families share a universal property of the phase diagram with superconductivity in proximity to magnetic phases Dai (2015); Keimer et al. (2015); Stewart (2017); Capone et al. (2009). Hence a theoretical formalism that treats the spin-fluctuation-mediated pairing and magnetic or other phases arising from strong repulsive interactions on an equal footing is highly desired. In the hole-doped cuprates, the phase diagram features spin- and charge-density waves, the pseudogap (PG) region, and the strange metal on top of the $d$-wave superconductivity. Although there exist various materials with different crystal structures within the cuprate family, the crucial ingredient in the whole family is the two-dimensional $\rm CuO_{2}$ plane Tsukada et al. (2005); Wesche (2015). The three-band model for the copper oxide is usually simplified into a two dimensional square lattice, which is widely believed to be the building block from which high-temperature superconductivity originates Schrieffer (2007). Countless theoretical investigations have addressed various experimental observations, specifically the structure of gap-functions in both superconducting phases Mesot et al. (1999); Fettes, W. and Morgenstern, I. (1999); Kitatani et al. (2015), pseudogap regions Ping et al. (2002); Matsuura et al. (2017); Gull et al. (2013); Wu et al. (2018), and the antiferromagnetic phase Schrieffer (2007); Krahl et al. (2009); Kuroki et al. (1999). Many-body numerical algorithms including, among others, (i) extensions of the mean-field approximation Krull et al. (2016); Sentef et al. (2017); Potthoff et al. (2003), (ii) various generalizations of the dynamical mean-field theory (DMFT) Gull et al. (2013); Kitatani et al. (2015); Jarrell et al. (2001); Maier et al. (2002); Staar et al. (2013); Ayral et al. (2013); Capone and Kotliar (2006); Kotliar et al. (2001), (iii) diagrammatic extensions Rubtsov et al. (2008); Vučičević et al. (2017); Kitatani et al. (2019), and (iv) quantum Monte Carlo (QMC) methods Tahara and Imada (2008); Kuroki and Aoki (1996); Ying et al. (2014). These approaches mostly involve various approximations, but there is no general consensus about the capability of the present approaches to investigate the multiple phases on an equal footing. This is imperative, since the glue for the pair formation as well as the origin of pseudogap phase should emerge by treating superconducting (SC) and antiferromagnetic (AF) phases in a unified framework.

Along with the advances of the above-mentioned numerical toolboxes, the phenomenologies approaches in terms of competing and/or synergistic order parameters based on generalizations of the Ginzburg-Landau (GL) theory have been developed for providing insights into various phases. Such approaches introduce the order parameters of the system based on symmetry considerations Zacher et al. (2000); Sun et al. (2005, 2008); Demler et al. (2004); Zhang (1997); Podolsky et al. (2004). In the GL theory, the respected symmetries of the relevant order parameters indeed govern the diversity of the associated phase diagram. As a result, for systems whose phase diagram is largely known from experiments, a reverse strategy is to find a large symmetry group which will be reduced to one of its subgroups upon the emergence of distinct phases. It is also crucial to terminate these subgroups such that all the smallest subgroups obey the conservation laws, namely conservation of charge ($\rm U(1)$) and spin (SU(2)s). It has been discussed that for families of high-temperature superconductors, the generators of the supergroup should be determined such that p- and d-wave superconductivity, staggered magnetization, charge density waves, spin, charge, and number operators can be expressed in their terms Guidry et al. (2001). For this reason, it has been argued that evoked symmetries can be classified Wu et al. (2003); Guidry et al. (2001) in terms of subgroups of SU(4) as

[TABLE]

Here, ${\rm SO(4)}\otimes{\rm U(1)}$ represents the AF order ($\rm SO(4)$) generated by the total and relative spin operators on even and odd sites, along with ${\rm U(1)}$ the charge group. SU(2)η is the second SU(2) symmetry inherent in the Hubbard model involving what is called $\eta$ spin Yang (1989). Indeed, ${\rm SU(2)}_{\rm s}\otimes{\rm SU(2)}_{\eta}$ introduces an algebra for describing the d-wave SC order as well as the AF order Sachdev et al. (2009); Chatterjee et al. (2017); Sachdev et al. (2019); Wu et al. (2018). (We briefly recall the symmetries associate with the different subgroups. More details will be given in the following.) The well-known SO(5) theory Zhang (1997); Demler et al. (2004) is one of the prominent examples, where the two-dimensional superconducting gap function and the three-component antiferromagnetic order constitute the generators. Introducing the primary orders in any subgroups of the SU(4) theory, one can interpret and predict the emergence of different orders. Although this approach can elegantly explain various aspects of the phase diagram as far as the Ginzburg-Landau (GL) type phenomenologies are concerned, a microscopic theory is necessary to give a solid basis to the GL approach and to provide information about the origin of the order parameters.

Thus it is crucial to combine a microscopic approach with a symmetry-group theoretic approach. One class of methods which is suitable to achieve this goal is formed by the slave-particle methods Kim (2007); Hermele (2007); Lechermann et al. (2007); Lee et al. (2006); Lee and Salk (2001); Wen and Lee (1996); de’Medici et al. (2005); Rüegg et al. (2010); Yu and Si (2012); de’ Medici (2017). In this framework, an electron operator is regarded to be comprising various auxiliary particles, where each component conveys a specific symmetry property of the original electron. As this procedure enlarges the Hilbert space, the factorized particles are subject to a constraint for suppressing unphysical states. Under such constraints, the method is dubbed the slave-particle decomposition.

A caution in advance: usual slave-particle works primarily use mean-field approaches. In the present work, by contrast, we use the slave-particles in the impurity solver for the FLEX-DMFT framework as we shall elaborate. Having said that, the slave building blocks so far considered are typically SU(2) slave-rotor Kim (2007); Hermele (2007), or SU(2) doublon-less slave-boson Lee et al. (2006); Lee and Salk (2001); Wen and Lee (1996), can then shed light on the origin of different order parameters, where each of the order parameters involves one of the introduced auxiliary particles. One should note that in these approaches it is a challenge to characterize the nature of order parameters which are not necessarily associated with one particular auxiliary particle. While the computational cost will obviously blow up as one increases the number of involved slave-particles, an appropriate choice may control the added cost. Of particular interest is the doublon-less SU(2) slave-boson decomposition, which was initially introduced to treat the $t$-$J$ model under the Gutzwiller projection that eliminates doubly-occupied states. This decomposition respects the rotational symmetry of real spins (which we denote SU(2)s). When one performs the mean-field treatment for this model, in the doublon-less SU(2) slave-boson picture, the solution shows that the pseudogap and superconducting transition temperatures can be attributed to the auxiliary spin-less and charge-less particles, respectively Lee and Salk (2001); Wen and Lee (1996). Align with the slave-particle approaches, SU(2) gauge theory of fluatuating antiferromagnetism Sachdev et al. (2009); Chatterjee et al. (2017); Sachdev et al. (2019); Wu et al. (2018) attempts to address the underlying physics of by fractionalizing the physical operators into spin-less and charge-less particles with SU(2) symmetry which can be spontanously broken by the condensation of the Higgs field.

However, while the results, obtained using the SU(2) gauge theory or the slave-particle methods, give a satisfactory picture of the various competing phases, these have been obtained within the mean-field treatment, which has to be examined if we are seriously interested in the correlation physics. For instance, one conspicuous feature in the phase diagram for the SU(2) doublon-less slave-boson is the superconducting Tc dome that is entirely covered by the pseudogap regime Lee and Salk (2001), which agrees with some experimental pieces of evidence Taillefer (2010); Kordyuk (2015), but it contrasts with other experiments Ramshaw et al. (2015). Thus this makes us to question if the result would be an artifact of mean-field approximations, and whether a more appropriate treatment of the electron correlation should be required for the d-wave SC and other orders. This has precisely motivated us to propose in the present paper a new approach, in which we introduce two concepts: (i) We first improve the slave-particle decomposition itself to allow double occupancies to study moderate repulsive interactions away from the strong-coupling limit, and (ii) we then apply the formalism not to a mean-field treatment, but to the FLEX+DMFT algorithm, which combines the fluctuation-exchange approximation (FLEX) and the dynamical mean-field theory (DMFT). We opt here for the FLEX+DMFT formalism Gukelberger et al. (2015); Kitatani et al. (2015), where FLEX can treat the momentum-dependent pairing interaction for d-wave SC, while DMFT can treat the Mott transition, thereby allowing us to incorporate both spatial (FLEX) and dynamical (DMFT) fluctuations.

Specifically, inspired by the phenomenological SU(4) theory of superconductivity Sun et al. (2005, 2008); Guidry et al. (2001), here we extend the slave-boson decomposition to capture

• (i)

the bipartite nature of the AF and d-wave SC phase, along with

• (ii)

SU(2)${}_{\rm s}\bigotimes$ SU(2)η symmetry [section(II) below].

We call the formalism a “bipartite full-SU(2) slave-boson”, which enables us to investigate AF and d-wave SC in the Hubbard model on an equal footing. Our approach substantially improves a previous mean-field work Sentef et al. (2017) which indeed invoked the SO(4) symmetry. While the present bipartite treatment maintains the SU(2)s and SU(2)η on a bipartite lattice as in Ref. (45), the latter is a mean-field treatment. Another difference is, instead of the O(4) rotor as the slave-particle employed in Ref. (45), we introduce here two species of bosons to convey charge-related information.

The organization of this paper is as follows. In Sec. II we introduce the model Hamiltonian and relevant symmetry properties. Section III starts with constructing a new SO(4) slave-boson formalism for Green’s function and other quantities that respects full-SU(2) symmetry and the bipartite structure. We then present the SO(4) DMFT+FLEX formalism by explaining the SO(4) DMFT method, followed by introducing the SO(4) FLEX. We shall then explain how the DMFT self-consistency loop is performed in the slave-particle space for the DMFT impurity solver. We conclude the paper in Sec. IV.

II Model Hamiltonian

The one-band Hubbard Hamiltonian on the square lattice is given by

[TABLE]

where $c^{\dagger}_{i\sigma}$ creates an electron with spin $\sigma$ at site $i$, and $n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$. The interaction ($U$) describes the on-site Coulomb repulsion between electrons with opposite spins. The density of electrons controlled by the doping is specified by the chemical potential ($\mu$). Electrons hop from site $i$ to $j$ on the lattice with the amplitude $v_{|i-j|}$, here taken into account up to the third-neighbors. We take the nearest-neighbor $v_{1}$ as the unit of energy hereafter.

This model has not only the usual SU(2) spin-rotational symmetry, but another SU(2) symmetry. For the latter, we can define a pseudospin (often called $\eta$-spin),

[TABLE]

where $\epsilon_{i}$, the bipartite factor, is defined on two (A,B) sublattices as

[TABLE]

If we rewrite the Hubbard interaction as

[TABLE]

we can immediately see that, at half-filling ($\mu=0$), $\sum_{i}{\bm{\eta}}_{i}$ commutes with the Hamiltonian.

In order to treat both of antiferromagnetism and d-wave superconductivity, namely, to incorporate the two SU(2) symmetries simultaneously, we then introduce, following Ref. Hermele (2007), a representation of an electron operator at site $i$,

[TABLE]

where we express $2\times 2$ matrices with tilde hereafter. Then the Hubbard Hamiltonian can be cast into a manifestly SU(2)$\otimes$SU(2) representation as

[TABLE]

where we have decomposed the hopping terms into those across different sublattices (the nearest-neighbor hopping) as represented by the first term, and those within the same sublattice (second- and third-neighbor hoppings) represented by the second term. The factor $\epsilon_{i}\epsilon_{j}$, arising from $\widetilde{C}^{\dagger}_{i}\widetilde{C}_{j}$ in these terms, depends on whether the hopping is inter- or intra-sublattice, so that for the former we have inserted $\tau^{z}$ to compensate $\epsilon_{i}\epsilon_{j}=-1$ in terms of Pauli matrices $\bm{\tau}=(\tau^{x},\tau^{y},\tau^{z})$ on the $2\times 2$ space in Eq.5. In this representation, ${\bm{\eta}}_{i}$ is the right-generator of the SU(2)ηYang (1989); Kitamura and Aoki (2016) as

[TABLE]

Similar to the $SU(2)_{\rm s}$ symmetry which allows transforming the spin-up electrons into spin-down electrons, SU(2)η symmetry enables doubly occupied states to be converted into empty states. Thus, at finite doping, with more electrons than number of sites, the SU(2)η symmetry will be lowered to the $U(1)$ charge symmetry. However, to enable our formalism to treat all doping regimes, including the half-filling, we incorporate the SU(2)η symmetry in our formalism. Even when we express the Hamiltonian in terms of $\widetilde{C}$ operators, the formalism does not enforce the symmetry and can describe the case of broken SU(2)η symmetry away from half-filling.

The usual spin-rotational SU(2)s, on the other hand, is the left-generator as

[TABLE]

Thus Eq. 6, at half-filling, enjoys a symmetry Feldbacher et al. (2003),

[TABLE]

where $\mathbb{Z}_{2}$ is the sublattice symmetry.

III Numerical method

In order to introduce a numerical formalism that treats AF and SC phases on an equal footing, we now formulate the FLEX+DMFT method Gukelberger et al. (2015) with our novel bipartite full-SU(2) slave-boson impurity solver in the ${\rm SO(4)}(\approx{\rm SU(2)_{s}}\otimes{\rm SU(2)_{\eta}}/\mathbb{Z}_{2})$ representation. To preserve the SO(4) symmetry in the whole formalism, we first introduce a spinor operator in the four-component Nambu (N) representation as

[TABLE]

where we again incorporate the bipartite factor $\epsilon_{i}$. This representation puts the two Nambu doublets together, namely

[TABLE]

with which Eq. 10 is expressed as

[TABLE]

As the superconducting gap function in the high-Tc cuprates is known to be singlet Keimer et al. (2015), we represent Eq. 10 in such a way that the singlet pairing correlation is included in the two-point propagator (see Eq. 17 below).

Since antiferromagnetism and d-wave pairing both involve neighboring sites, which we call $i$ and $i+\hat{x}$ (and $i+\hat{y}$) where $\hat{x}(\hat{y})$ is the nearest-neighbor vector along $x(y)$ axis, we can cast, for treating them on an equal footing, the ${C}_{{\rm N}i}$ defined in Eq.(10) into a bipartite representation where, at cluster $i$, ${C}_{{\rm N}i}$ involves both of $(A,B)$ sublattice sites. We thus redefine Eq. 10 as

[TABLE]

where $c_{iA\sigma}$ denotes electron at sublattice $A$ on cluster site $i$ with spin $\sigma$. To incorporate the rotational symmetry of our square lattice, we take $c_{iB\sigma}\equiv(c_{i+\hat{x}A\sigma}+c_{i+\hat{y}A\sigma})/2$, where the factor $2$ imposes having only two sites per unit cell, and $\widetilde{\epsilon}_{B}\equiv(\epsilon_{i+\hat{x}}+\epsilon_{i+\hat{y}})/2$. Note that the coordinate of cluster $i$ ($\bm{r}_{\rm i}$) is the same as its $A$ sublattice ($\bm{r}_{\rm iA}$), see also Fig. 1. Applying the Fourier transform to the electron operators on each sublattice as $c^{\dagger}_{\bm{k}A(B)}=(1/N)\sum\limits_{l}^{\rm A(B)}e^{i\bm{k}\cdot\bm{r}_{l}}c^{\dagger}_{lA(B)}$ with momentum $\bm{k}\equiv(k_{x},k_{y})$ and $N$ being the total number of lattice sites, we obtain

[TABLE]

Here we have used $\widetilde{\epsilon}_{B}=e^{i\bm{Q}\cdot\bm{r}_{\rm iB}}$ with $\bm{Q}=(\pi,\pi)$ being the AF (Brillouin-zone corner) wave vector, and $\bm{r}_{\rm iB}$ being the coordinate of the sublattice B on bipartite site $i$.

To calculate Green’s function, we use the the Heisenberg representation of electrons, $c(\tau)=e^{\tau H}ce^{-\tau H}$, where $H$ is given by Eq. 1, and $\tau$ is the Matsubara time. In the present four-component Nambu representation, the Green’s function is a $4\times 4$ matrix, ${\cal G}_{4\times 4}$, which is given in terms of the normal component defined as

[TABLE]

along with the anomalous component for treating SC phases,

[TABLE]

Here $\overline{\tau}=\tau-\tau^{\prime}$, $\alpha$ and $\beta$ denote sublattice indices, ${\cal T}_{\tau}$ is the time-ordering operator, and $\overline{\sigma}=-\sigma$. In the momentum space the $4\times 4$ Green’s function is concisely expressed as

[TABLE]

with $\overline{\tau}\equiv\tau-\tau^{\prime}$, for which the matrix elements are explicitly given as

[TABLE]

Equivalently one can transform the imaginary-time Green’s function into Matsubara-frequency space as

[TABLE]

with $\beta=1/T$, where the fermionic Matsubara frequencies are given by $\omega_{n}=\beta(2n-1)/\pi$.

The lattice Green’s function in Eq. 18 satisfies the Nambu-Dyson equation,

[TABLE]

where $k\equiv(\bm{k},i\omega_{n})$, and $\mathbb{I}_{4\times 4}$ denotes the identity matrix. In the above, $\varepsilon_{4\times 4}(\bm{k})$ is the dispersion due to the hopping terms in Eq. 1,

[TABLE]

with

[TABLE]

where $v_{1},v_{2},v_{3}$ are first, second, and third neighbor hopping amplitudes, respectively.

Another term, ${\Sigma}^{\rm eff}_{4\times 4}(k)$, in Eq. 19 is the self-energy, and in FLEX+DMFT formalism which we adopt here, this comprises a combination of the self-energies in DMFT and FLEX,

[TABLE]

where we have defined the non-local part of the FLEX self-energy as

[TABLE]

where we subtract the local part of the FLEX contribution, $\Sigma_{4\times 4}^{\rm FLEX,local}$, to avoid double counting of local Feynman diagrams.Kitatani et al. (2015) This local contribution is obtained by following the FLEX self-energy prescription as introduced below [Eq.(48)] after substituting the lattice Green’s functions with the local Green’s function,

[TABLE]

The local Green’s function has a matrix representation,

[TABLE]

where $G_{AB\sigma\sigma^{\prime}}$, and $F_{AB\sigma\sigma^{\prime}}$ are defined as

[TABLE]

Here $c_{A(B)\sigma})$ denotes the annihilation operator of an electron with spin $\sigma$ residing on sublattice $A\leavevmode\nobreak\ (B)$ of the impurity site.

The FLEX+DMFT self-consistency iteration is described in detail in the following sections, where the outline, see Fig. 2, is as follows:

Initialize $\Sigma^{\rm nl}_{4\times 4}$, and obtain the lattice Green’s function (${\cal G}_{4\times 4}$) without the impurity self-energy. 2. 2.

Insert the Green’s function into the FLEX loop to compute a new $\Sigma^{\rm FLEX}$, which is then plugged in the ${\rm DMFT}$ iteration for computing the hybridization function. 3. 3.

With the hybridization function, solve the impurity problem to evaluate the impurity self-energy. 4. 4.

The effective lattice self-energy is then determined by summing the impurity and nonlocal FLEX self-energies. 5. 5.

With $\Sigma^{\rm eff}$, Eq. 22, we solve the Dyson equation, Eq. 19, to update the lattice Green’s function. 6. 6.

Repeat the above double (FLEX+DMFT) self-consistency loops until the convergence is attained.

A difference from the FLEX+DMFT in Ref.Kitatani et al. (2015) is that here we explicitly treat the hybridization function ($\cal{D}$ in Fig. 2).

III.1 DMFT with the bipartite full-SU(2) slave-boson solver

In the usual DMFT, the many-body problem is mapped onto an impurity problem which is embedded in a self-consistent medium (bath). In this procedure, we have thus a mean-field treatment in real space, while we do retain temporal (dynamical) quantum fluctuations, which enables the method to treat Mott transitions. Thus the approach is a computationally less demanding algorithm. We shall later improve the method by combining with the FLEX framework to incorporate spatial fluctuations. An implementation of DMFT requires to solve the impurity model and compute the Green’s function. Several impurity solvers have been implemented, each with its own advantages and disadvantages. Numerically exact approaches including continuous-time Quantum Monte Carlo, exact diagonalization and Numerical Renormalization Group provide numerically accurate results at least for single-site DMFT Georges et al. (1996), but the numerical cost increases rapidly with the number of orbitals and sites in the unit cell and/or in approaches where DMFT is combined with non-local effects. Approximate analytical methods can provide the desired information at a smaller cost.

To allow our DMFT framework to treat AF and SC on an equal footing, we should consider a two-site impurity cluster (inset of Fig.1). Thus we expound here the SO(4) DMFT steps to solve the two-site impurity problem using the bipartite full-SU(2) slave-boson solver.

The impurity action for the Hamiltonian (1) is given by

[TABLE]

where, in this section, the electron annihilation ($c$) and number ($n$) operators refer to the impurity cluster, and the mathfrac site index ${\mathfrak{i}}$ refers to A and B sublattice sites in the impurity cluster. The impurity Hamiltonian reads

[TABLE]

$\mathfrak{G}$ and $\mathfrak{F}$ appearing in Eq.(28) are normal and anomalous hybridization functions, respectively, which convey the bath information, and are expressed, in a $4\times 4$ matrix form in the present formalism, as

[TABLE]

where $\mathfrak{G}(\mathfrak{F})$ appear in the diagonal (off-diagonal) blocks. In the FLEX+DMFT iteration we determine this matrix through

[TABLE]

where “$*$” denotes the convolution integral on the Matsubara frequency, and we have employed Eq.(25) to obtain Eq.(30) (see Appendix for detailed derivation). If we define, as in ${\cal G}^{\rm loc}_{4\times 4}$ in Eq.(25), the impurity Green’s function in a $4\times 4$ matrix form with

[TABLE]

we can then impose that

[TABLE]

in the DMFT scheme. To solve the impurity problem, we shall introduce in the next section the bipartite full-SU(2) slave-boson impurity solver.

III.1.1 Bipartite full-SU(2) slave-boson impurity solver

Within the slave-particle formalism, we decompose the electron operator into fermionic and bosonic particles such that the associated matrix elements of physical states be equal to their counterparts in the auxiliary slave-particle Hilbert space. Within this prescription, the fermion statistics for the physical particle may or may not be satisfied depending on what kind of slave-particles we adopt. The doublon-less SU(2) slave-boson, in particular, represents the electron operator as a composite of charge-less fermions (spinons) with up and down spins, and two-flavored spin-less bosons (holons) with flavor indices ($1,2$). The two species of holons are required to maintain the SU(2)s symmetry of the system. The doublon-less SU(2) slave-boson representation is actually introduced to examine the limit of the large repulsive interaction, where we eliminate double occupancies by applying the Gutzwiller projection Lee and Salk (2001); Lee et al. (2006); Wen and Lee (1996). One consequence of this is that the fermion statistics is violated. To remedy this, here we introduce another set of spin-less bosons (doublons) that take care of double occupancies. This procedure will also incorporate the SU(2)η symmetry into the formalism. Then the original electron creation and annihilation operators are expressed as

[TABLE]

where we have introduced the spinon operator ($f$) that satisfies the fermionic commutation relation, along with the holon operators $(b_{1},b_{2})$ and doublon operators ($d_{1},d_{2}$) that respectively satisfy the bosonic commutation relation. In the above we have maintained the bipartite formalism, where the index in mathfrac (${\mathfrak{i}}\in\{A,B\}$) refers to impurity sites, and $\overline{\sigma}\equiv-\sigma$. Then the fermionic commutation relation for the original electron operator is preserved. This decomposition thus realizes our perception of the real-spin SU(2) symmetry described by spinons along with the charge ($\eta$) SU(2) symmetry taken care of by creation and annihilation of holon and doublon operators. On top of this, the $\mathbb{Z}_{2}$ symmetry is imposed by exploiting the bipartite factor ($\epsilon_{\mathfrak{i}}$) in the above representation. We can make its structure more transparent by adopting the $2\times 2$ matrix forms as in Eq.5 to have

[TABLE]

The bipartite full-SU(2) slave-boson representation has a Hilbert space that is larger than the physical one, so we need to eliminate the unphysical states by a constraint, which can be obtained as follows You et al. (2012). For the SU(2) gauge-invariance of the electron doublet ($\widetilde{C}$ above) we can construct the SU(2) generators, ${\widetilde{K}}_{x/y/z}$, in such a way that we have commutation relations,

[TABLE]

We can show that this equation has a solution in a vector form,

[TABLE]

where the three components have explicit forms of

[TABLE]

The noninteracting states ($|\rangle$’s) of the system on a bipartite unit cell can also be translated in the SU(2) slave-boson decomposition as

[TABLE]

where $\mathfrak{i}\in\{A,B\}$, and the vacuum of the bipartite impurity site is defined as $|{\rm vac}_{A},{\rm vac}_{B}\rangle\equiv|{\rm vac}\rangle_{A}\otimes|{\rm vac}\rangle_{B}$. Here the slave-boson vacuum is decomposed as

[TABLE]

where $|{\rm vac}\rangle_{f\mathfrak{i}(b\mathfrak{i},d\mathfrak{i})}$ is the associated vacuum of the spinons (bosons) at impurity site $\mathfrak{i}$. One should note that the bipartite SU(2) impurity solver, with the bipartite factor $\epsilon_{\mathfrak{i}}$ inserted, results in a sublattice dependence of the projected states in Eq.(42). In addition, one should note that the transformed states of the SU(2) slave-boson can be occupied by more than one slave-particles. Obviously, removing $d$-bosons from our slave-boson picture would exclude the doubly-occupied sites, which is precisely the distinction between the present formalism and the conventional doublon-less SU(2) slave-boson Wen and Lee (1996); Lee and Salk (2001).

The total density of holes at $\mathfrak{i}$ is equivalent to the total number of holons,

[TABLE]

Similarly, the double occupancy is expressed as

[TABLE]

The impurity Hamiltonian (29) can be expressed up to a constant as

[TABLE]

where $H_{\rm hop}$ is the hopping term in Eq. 29 translated into the slave-boson language, see Supplementary material, while ($\lambda^{z}_{A/B},\lambda^{+}_{A/B},\lambda^{-}_{A/B}$) are Lagrange multipliers introduced to impose the SU(2) constraints in Eqs. 41 on each sublattice. Now we can express the action in Eq. 28 in the slave-boson language as

[TABLE]

Employing this action we can derive the equations of motion for the spinon, holon, and doublon operators, respectively, for details see Supplementary material. If Lagrange multipliers are correctly determined, the obtained results should satisfy $\langle\bm{K}\rangle_{i}={\bf 0}$ in terms of expectation values, which approximately satisfy Eq. 40. Alternatively, the auxiliary equations of motion can be iteratively solved by a new set ($\lambda_{A/B}^{+/-/z}$) until constraints are satisfied. With the final slave-particle Green’s functions, we can update the electron’s Green’s functions (${\cal G}_{4\times 4}^{\rm loc}$), see again Supplementary material. Finally, the DMFT self-energy is obtained as

[TABLE]

III.2 Fluctuation-exchange approximation (FLEX)

The fluctuation-exchange approximation provides a self-energy of an interacting system by summing over bubble and ladder diagrams. The self-energy in this formalism can be obtained from close-linked diagrams known as the Luttinger-Ward functional $\Phi$, so that the scheme is a conserving approximation Wang et al. (2015). Studies of the normal states of the Hubbard model reveal that including only the particle-hole channels as dominating contributors enables us to study the properties of the incommensurate antiferromagnetic spin structure as well as the superconducting instabilities in the overdoped regime Kuroki et al. (1999); Yanase and Yamada (2001); Oka and Aoki (2010); Kitatani et al. (2015). It has been discussed that, due to insufficiently treated dynamical and pairing fluctuations in this formalism, the superconducting dome in the phase diagram can only be partially captured with the AF region overestimated Yanase and Yamada (2001).

In the following, we shall propose an extended SO(4) FLEX self-energy, which treats the superconducting pairing and spin fluctuations on an equal footing. Hence not only the fluctuations in the electron pairs having momenta $\bm{k}$ and $-\bm{k}$, as treated in Ref. (65), but also the pairing between $\bm{k}$ and $-(\bm{k}+\bm{Q})$, known as $\eta$-pairing Sentef et al. (2017), are incorporated into the present formalism.

Now we delve into the present FLEX formalism by introducing the FLEX self-energy in a $4\times 4$ matrix form as

[TABLE]

Here $k+\bm{Q}\equiv(\bm{k}+\bm{Q},i\omega_{n})$, the Hartree self-energy ($\Sigma^{\rm H}$) is a diagonal matrix with momentum-independent elements as

[TABLE]

which is just equal to $Un_{4\times 4}$ with $n_{4\times 4}$ being the diagonal elements of ${\cal G}_{4\times 4}$ in Eq. 17 at $\tau=\beta$.

The normal self-energy,

[TABLE]

and the anomalous self-energy,

[TABLE]

are obtained from the functional derivative of the Luttinger-Ward functional $\Phi$, for which the expansion up to the forth-order in $U$ is displayed in Fig. 3. Here, the normal component of the self-energy consisting of the particle-hole and transverse-spin contributions is given by

[TABLE]

while the anomalous component has

[TABLE]

Here, $N$ is the total number of lattice sites, and

[TABLE]

which comes from the structure of Eqs(17,48) and we have to make an appropriate choice depending on which matrix element in Eq(48) is considered, and $q\equiv(\bm{q},i\nu_{n})$ with $\nu_{n}=2\pi n/\beta$ being the Matsubara frequency for bosons.

Since the Luttinger-Ward functional here incorporates the anomalous part with the anomalous self-energy given in Eq. 51, we should consider the local correction to the anomalous self-energy $\Sigma^{F}$ as $\Sigma^{F}_{\rm FLEX+DMFT}=\Sigma^{F}_{\rm FLEX}+\Sigma^{F}_{\rm loc}$. Now, our interest here is the anisotropic, $d$-wave pairing instability in the repulsive model, for which we can ignore the local correction to the anomalous self-energy $\Sigma^{F}_{\rm loc}$ which does not depend on momentum. The remaining term, $\Sigma^{F}_{\rm FLEX}=\delta\Phi_{\rm FLEX}[G,F^{\dagger},F]/\delta F^{\dagger}$, is the same as the right-hand side of the linearized Eliashberg equation (III.2 below) if we linearize the anomalous part. Then our formalism treats the normal and anomalous self-energies consistently, as functional derivatives of the same Luttinger-Ward functional $\Phi_{\rm FLEX+DMFT}$ Kamenev (2011); Kitatani et al. (2015).

The particle-hole (ph), and superconducting pairing (scp) vertex functions are computed as Yan (2005); Yan and Ting (2006)

[TABLE]

[TABLE]

There, the spin ($s$) and charge ($c$) susceptibilities are given by

[TABLE]

with

[TABLE]

In order to obtain the transition temperature for the converged Green’s functions, we solve the linearized Eliashberg equation,

[TABLE]

where $\lambda_{\rm EL}$ denotes the largest eigenvalue of the linearized Eliashberg equation, $\Delta_{\sigma\overline{\sigma}}$ is the the gap function, and $q=\overline{k}-\overline{k}^{\prime}$ with choices of $\overline{k}$ and $\overline{k}^{\prime}$ depending on the involved normal ($G$) matrix elements of the Green’s function ${\cal G}_{4\times 4}$ in Eq. III.2. The effective singlet-pairing interaction $V^{\rm eff}$ is

[TABLE]

The superconducting transition occurs when the maximum eigenvalue of the diagonalized Eliashberg equation reaches $\lambda_{\rm EL}=1$.

IV Conclusion

In conclusion, we have proposed a novel formalism to explore correlated systems as exemplified by the one-band repulsive Hubbard model. We have presented the SO(4) generalization of the FLEX+DMFT method so that both antiferromagnetism and superconductivity can be treated on an equal footing. Namely, in the FLEX+DMFT formalism, we solve the bipartite impurity problem in the SO(4) DMFT by introducing a novel “full-SU(2)” slave-boson impurity solver. This impurity solver respects the group-symmetry properties of the Hubbard model, namely spin SU(2) and pseudospin SU(2) symmetries. We have introduced bosonic and fermionic auxiliary particles to convey all the charge- and spin-related information in our impurity solver. This approach is particularly suitable to study the interplay between AFM and d-wave SC treating the two phases on equal footing, while previous calculations based on cluster extension of DMFT Capone and Kotliar (2006) can suffer of some bias due to the choice of the cluster. Going over to FLEX+DMFT with the slave-boson framework with double self-consistent loops then incorporates improved k-dependent fluctuations. As our slave-particle decomposition is treated within the FLEX+DMFT approach, addressing correlated physics of antiferromagnetism and d-wave superconductivity is feasible. Furthermore, the decomposed nature of our impurity solver may shed light on the origin of, still puzzling, pseudogap physics. Numerical study will be desirable as a future work.

Extending our formalism to address multi-band physics Bunemann (2011) is another exciting direction to pursue. This will enable us not only explore the three-band model, more suitable for cuprates, but also enable us to examine multi-band superconductors, e.g., iron-based superconductors Dai (2015), where the superconducting phase also sits adjacent to the magnetic phase.

Acknowledgements.

We thank Martin Eckstein and Abolhassan Vaezi for a collaboration in early stages of this work. We also appreciate fruitful discussions with Motoharu Kitatani. H.A. wishes to thank hospitality of Department of Physics, ETH Zürich, where the present manuscript was started. Sh.S. would like to thank Patrick A. Lee for pointing out the possibility of adding one more boson to the formalism. She also appreciates Makiko Nio and Tatsumi Aoyama for helpful instructions on performing diagrammatic calculations with FORM Vermaseren (2000); Aoyama et al. (2012). Sh.S. and H.A. acknowledge a support from the ImPACT Program of the Council for Science, Technology and Innovation, Cabinet Office, Government of Japan (Grant No. 2015-PM12-05-01) from JST. H.A. is also supported by JSPS KAKENHI Grant No. JP26247057. M.C. acknowledges support from H2020 Frame-work Programme, under ERC Advanced Grant No.692670 “FIRSTORM” and MIUR PRIN 2015 (Prot.2015C5SEJJ001) and SISSA/CNR project ”Superconductivity, Ferroelectricity and Magnetism in bad metals” (Prot.232/2015).

I Bipartite full-SU(2) slave-boson impurity solver for the DMFT

Employing the notation of the bipartite full-SU(2) slave-boson, we can rewrite the impurity Hamiltonian in Eq. (29) as

[TABLE]

where ($\lambda^{z},\lambda^{+},\lambda^{-}$) are Lagrange multipliers for enforcing the SU(2) vector constraints in Eq. (40). The SU(2) nature of the doubly-occupied state imposes the sextic terms (the last term in the first line) in the impurity Hamiltonian. As this term is local in time for simplicity we would apply a mean-field contraction in the spinon sector such that

[TABLE]

where density operators read

[TABLE]

This impurity Hamiltonian is associated with the action in the slave-boson language,

[TABLE]

Here electron operators are defined in Eq.(36) in the main text, and the $4\times 4$ matrix representation of the hybridization function has the form of

[TABLE]

Using Eq. (S4), we derive the equations of motion for all of spinon, holon, and doublon flavors. Here we employ the generic notation for the auxiliary Green’s functions,

[TABLE]

where $a\text{ and }a^{\prime}\in\{f,b,d\}$ label the auxiliary particles, $\alpha$ is the flavor of the particular slave-particle which is $\{\uparrow,\downarrow\}$ for spinons and $\{1,2\}$ for bosonic operators, and $i,j$ denote sublattice indices. When Green’s functions consist only one species of the auxiliary particle, namely both $a$ and $a^{\prime}$ are spinon, holon, or doublon, we refrain from repeating these labels for the Green’s functions as $G^{aa^{\prime}}\leavevmode\nobreak\ (F^{aa^{\prime}})$, with a shorthand $G^{aa}\equiv G^{a},F^{aa}\equiv F^{a}$ when $a=a^{\prime}$. The self-energy diagrams for these auxiliary Green’s function, denoted generically as $\cal G^{\rm aux}$, are exemplified in Figs. S1 and S2. Using the auxiliary particles, we can construct the Green’s functions as presented in Sec. I.13 and plotted in Fig. S3.

I.1 Equation of motion for $f_{A\uparrow}$ particles

[TABLE]

I.2 Equation of motion for $f_{A\downarrow}$ particles

[TABLE]

I.3 Equation of motion for $f_{B\uparrow}$ particles

[TABLE]

I.4 Equation of motion for $f_{B\downarrow}$ particles

[TABLE]

I.5 Equation of motion for $b_{A1}$ particles

[TABLE]

I.6 Equation of motion for $b_{2A}$ particles

[TABLE]

I.7 Equation of motion for $b_{1B}$ particles

[TABLE]

I.8 Equation of motion for $b_{2B}$ particles

[TABLE]

I.9 Equation of motion for $d_{A1}$ particles

[TABLE]

I.10 Equation of motion for $d_{2A}$ particles

[TABLE]

I.11 Equation of motion for $d_{1B}$ particles

[TABLE]

I.12 Equation of motion for $d_{2B}$ particles

[TABLE]

I.13 Electron Green’s functions in the slave-boson language

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]