TRILEX and $GW$+EDMFT approach to $d$-wave superconductivity in the Hubbard model

TL;DR

This paper extends the TRILEX method to superconducting phases in the Hubbard model, successfully capturing $d$-wave superconductivity at strong coupling and exploring its dependence on doping, dispersion, and AF correlations.

Contribution

The paper introduces a generalized TRILEX approach for superconductivity, capable of modeling $d$-wave pairing with a single-site impurity, and compares it with $GW$+EDMFT results.

Findings

TRILEX produces a $d$-wave superconducting dome at strong coupling.

The method captures $d$-wave pairing using only a single-site impurity model.

Superconducting temperature depends on hopping amplitudes and AF correlations.

Abstract

We generalize the recently introduced TRILEX approach (TRiply Irreducible Local EXpansion) to superconducting phases. The method treats simultaneously Mott and spin-fluctuation physics using an Eliashberg theory supplemented by local vertex corrections determined by a self-consistent quantum impurity model. We show that, in the two-dimensional Hubbard model, at strong coupling, TRILEX yields a $d$-wave superconducting dome as a function of doping. Contrary to the standard cluster dynamical mean field theory (DMFT) approaches, TRILEX can capture $d$-wave pairing using only a single-site effective impurity model. We also systematically explore the dependence of the superconducting temperature on the bare dispersion at weak coupling, which shows a clear link between strong antiferromagnetic (AF) correlations and the onset of superconductivity. We identify a combination of hopping amplitudes particularly favorable to superconductivity at intermediate doping. Finally, we study within $GW$+EDMFT the low-temperature $d$-wave superconducting phase at strong coupling in a region of parameter space with reduced AF fluctuations.

Equations

$$H=\sum _{ij\sigma }{t}_{ij}{c}_{i\sigma }^{†}{c}_{j\sigma }-\mu \sum _{i\sigma }{n}_{i\sigma }+U\sum _i{n}_{i\uparrow }{n}_{i\downarrow }$$

$$t_{ij}=\left\{\begin{array}[]{cc}t,&\mathbf{r}_{i}=\mathbf{r}_{j}\pm\mathbf{e}_{x,y}\\ t^{\prime},&\mathbf{r}_{i}=\mathbf{r}_{j}\pm\mathbf{e}_{x}\pm\mathbf{e}_{y}\\ t^{\prime\prime},&\mathbf{r}_{i}=\mathbf{r}_{j}\pm 2\mathbf{e}_{x,y}\\ 0,&\mathrm{otherwise.}\end{array}\right.$$

$${\epsilon }_{\mathbf{k}}$$

$$+2t^{\prime \prime }(\cos 2k_x+\cos 2k_y)$$

$$S_{\mathrm{eb}}[c,c^{\ast },\varphi ]$$

$$+{\lambda }_{\mu \nu \alpha }{c}_{\mu }^{\ast }{c}_{\nu }{\varphi }_{\alpha }$$

$$\boldsymbol{\Psi}_{i}(\tau)\equiv\left[\begin{array}[]{c}c_{\uparrow i}^{*}(\tau)\\ c_{\downarrow i}(\tau)\\ c_{\downarrow i}^{*}(\tau)\\ c_{\uparrow i}(\tau)\end{array}\right]$$

$$S_{\mathrm{eb}}^{\mathrm{Nambu}}[\mathbf{\Psi },\varphi ]$$

$$+\frac{1}{2}{\varphi }_{\alpha }{\mathbf{\Psi }}_u{\mathbf{\lambda }}_{uv\alpha }{\mathbf{\Psi }}_v$$

$$\int \mathcal{D}[\mathbf{\Psi }]e^{\frac{1}{2}{\mathbf{\Psi }}_u{\mathbf{A}}_{uv}{\mathbf{\Psi }}_v}={\left(\det \mathbf{A}\right)}^{\frac{1}{2}}$$

$${\mathbf{G}}_{uv}$$

$$W_{\alpha \beta }$$

$${\mathbf{\chi }}_{uv\alpha }^{3,\mathrm{conn}}$$

$${\mathbf{\Lambda }}_{uv\alpha }\equiv {\left[{\mathbf{G}}^{-1}\right]}_{uw}{\left[{\mathbf{G}}^{-1}\right]}_{xv}{\left[W^{-1}\right]}_{\alpha \beta }{\mathbf{\chi }}_{wx\beta }^{3,\mathrm{conn}}$$

$$Z=\int \mathcal{D}[\mathbf{\Psi },\varphi ]e^{-S_{\mathrm{eb}}^{\mathrm{Nambu}}[\mathbf{\Psi },\varphi ]}=\int \mathcal{D}[c,c^{\ast },\varphi ]e^{-S_{\mathrm{eb}}[c,c^{\ast },\varphi ]}$$

$${\mathbf{\Sigma }}_{uv}$$

$$P_{\alpha \beta }$$

$$U{n}_{i\uparrow }{n}_{i\downarrow }=\frac{1}{2}\sum _I{U}^I{n}_i^I{n}_i^I$$

$$U^{\mathrm{ch}}-U^{\mathrm{sp}}$$

$$U^{\mathrm{ch}}-3U^{\mathrm{sp}}$$

$$U^{\mathrm{ch}}=U/2,U^{\mathrm{sp}}=-U/2$$

$$U^{\mathrm{ch}}=U/2,U^{\mathrm{sp}}=-U/6.$$

$$\boldsymbol{G}_{0,ij}(\tau)=\left[\begin{array}[]{cccc}0&0&0&-G_{0,ji}(-\tau)\\ 0&0&G_{0,ij}(\tau)&0\\ 0&-G_{0,ji}(-\tau)&0&0\\ G_{0,ij}(\tau)&0&0&0\\ \end{array}\right]$$

$$G_{0,ij}(\tau )$$

$$G_{0\mathbf{k}}(i\omega )$$

$${\mathbf{\lambda }}_{uv\alpha }$$

$$\displaystyle\boldsymbol{\delta}_{\tau_{u},\tau_{v}}=\left[\begin{array}[]{cccc}\delta_{\tau_{u},\tau_{v}^{+}}&&&\\ &\delta_{\tau_{u}^{+},\tau_{v}}&&\\ &&\delta_{\tau_{u},\tau_{v}^{+}}&\\ &&&\delta_{\tau_{u}^{+},\tau_{v}}\end{array}\right]$$

$$\boldsymbol{\lambda}^{I}=\left[\begin{array}[]{cccc}&\sigma_{\uparrow\downarrow}^{I}&&\sigma_{\uparrow\uparrow}^{I}\\ -\sigma_{\uparrow\downarrow}^{I}&&-\sigma_{\downarrow\downarrow}^{I}\\ &\sigma_{\downarrow\downarrow}^{I}&&\sigma_{\downarrow\uparrow}^{I}\\ -\sigma_{\uparrow\uparrow}^{I}&&-\sigma_{\downarrow\uparrow}^{I}\end{array}\right]$$

$$W_{0,ij}^I(\tau )={\delta }_{ij}{\delta }_{\tau }{U}^I$$

$$S_{\mathrm{eb}}^{\mathrm{Nambu}}[\mathbf{\Psi },\varphi ]=$$

$$\frac{1}{2}\sum _{i,j,a,b}\iint \mathrm{d}\tau \mathrm{d}{\tau }^{\prime }{\mathbf{\Psi }}_{ia}(\tau )[-{\mathbf{G}}_0^{-1}{]}_{ia,jb}(\tau -{\tau }^{\prime }){\mathbf{\Psi }}_{jb}({\tau }^{\prime })$$

$$+\frac{1}{2}\sum _i\sum _I\int \mathrm{d}\tau {\varphi }_i^I(\tau )[-(U^I{)}^{-1}]{\varphi }_i^I(\tau )$$

$$+\frac{1}{2}\sum _i\sum _I\int \mathrm{d}\tau {\varphi }_i^I(\tau ){\mathbf{\Psi }}_{ia}(\tau ){\mathbf{\lambda }}_{ab}^I{\mathbf{\Psi }}_{ib}(\tau )$$

$${\mathbf{\Sigma }}_{ab,\mathbf{k}}(i\omega )$$

$$\times W_{\mathbf{q}}^I(i\Omega ){\mathbf{\Lambda }}_{db,\mathbf{kq}}^I(i\omega ,i\Omega ),$$

$$P_{\mathbf{q}}^I(i\Omega )$$

$$\times {\mathbf{G}}_{cd,\mathbf{k}}(i\omega ){\mathbf{\Lambda }}_{db,\mathbf{kq}}^I(i\omega ,i\Omega ).$$

$${\mathbf{\Lambda }}_{\mathbf{kq},ab}^I(i\omega ,i\Omega )$$

$$⟨c_{\uparrow }^{\ast }(\tau )c_{\uparrow }^{\ast }(0)⟩=⟨c_{\downarrow }^{\ast }(\tau )c_{\downarrow }^{\ast }(0)⟩=0$$

Full text

Strongly-correlated electrons systems such as high-temperature superconductors pose a difficult challenge to condensed-matter theory.

One class of theoretical approaches for this problem focuses on the effect of long-range spin-fluctuationsChubukov et al. (2002); Efetov et al. (2013); Wang and Chubukov (2014); Metlitski and Sachdev (2010); Onufrieva and Pfeuty (2009, 2012). They neglect vertex corrections in an Eliashberg-like approximation for the electronic self-energy and predict a $d$-wave superconducting order.

Another class of approaches focuses, following the seminal work of AndersonAnderson (1987), on the fact that high-temperature superconductors are doped Mott insulators. In the recent years, progress has been made in this direction with cluster extensionsHettler et al. (1998, 2000); Lichtenstein and Katsnelson (2000); Kotliar et al. (2001); Maier et al. (2005a) of dynamical mean field theory (DMFT) Georges et al. (1996). These methods have been shown to capture the essential aspects of cuprate physics, such as Mott insulating, pseudogap and $d$-wave superconducting phasesKyung et al. (2009); Sordi et al. (2012a); Civelli et al. (2008); Ferrero et al. (2010); Gull et al. (2013a); Macridin et al. (2004); Maier et al. (2004, 2005b, 2006a); Gull et al. (2010); Yang et al. (2011); Macridin and Jarrell (2008); Macridin et al. (2006); Jarrell et al. (2001); Bergeron et al. (2011); Kyung et al. (2004, 2006); Okamoto et al. (2010); Sordi et al. (2010, 2012b); Civelli et al. (2005); Parcollet et al. (2004); Ferrero et al. (2008, 2009); Gull et al. (2009); Chen et al. (2015, 2017). Cluster DMFT methods can be converged with respect to the cluster size at relatively high temperatureLeBlanc et al. (2015); Koch et al. (2008), including in the pseudogap region Wu et al. , but not at lower temperatures and in particular in the superconducting phase.

Several approaches beyond cluster DMFT have been proposed recently Rubtsov et al. (2008, 2012); van Loon et al. (2014); Stepanov et al. (2016); Toschi et al. (2007); Katanin et al. (2009); Schäfer et al. (2015); Valli et al. (2015); Li et al. (2016); Rohringer and Toschi (2016); Ayral and Parcollet (2016a); Biermann et al. (2003); Sun and Kotliar (2002, 2004); Ayral et al. (2012, 2013); Biermann (2014); Ayral et al. ; Rohringer et al. . In Refs Ayral and Parcollet, 2015, 2016b, the Triply Irreducible Local Expansion (TRILEX) approach was introduced. It consists in a local approximation of the electron-boson vertex extracted from a quantum impurity model with a self-consistently determined bath and interaction, in the spirit of DMFT. TRILEX interpolates between DMFT at strong interaction and the weak-coupling Eliashberg-like spin-fluctuation approximation at weak interaction. It is able to simultaneously describe Mott physics and the effect of long-range bosonic fluctuations. Hence, it unifies the two theoretical approaches mentioned above in the same formalism.

The main purpose of this paper is to study $d$-wave superconductivity in the Hubbard model within the single-site TRILEX approach. Contrary to DMFT, where $d$-wave superconducting correlations can by construction be captured only within multi-site (cluster) impurity models, here we only need to solve a single-site impurity model. We also compare TRILEX to two simpler approaches, $GW$+EDMFT and $GW$, which can be viewed as further approximations of the electron-boson vertex in TRILEX. We show that TRILEX yields a $d$-wave superconducting dome at strong coupling.

We also study the dependence of the superconducting critical temperature $T_{c}$ on the choice of the tight-binding parameters at weak coupling using the $GW$ method. While $T_{c}$ is enhanced by strong antiferromagnetic fluctuations, we find a region of parameter space where the superconducting transition occurs at a higher temperature than the antiferromagnetic instability of the method. At this point, we stabilize and study a superconducting solution below $T_{c}$ within $GW$+EDMFT. We also identify a choice of dispersion where, at $16\%$ doping, we have a pronounced maximum of $T_{c}$ in the space of hopping parameters, which seems to persist even at strong coupling.

The paper is organized as follows: In Section I, we describe the Hubbard model studied in this paper. In Section II, we generalize the TRILEX equations to superconducting phases via the Nambu formalism, and discuss their simplifications $GW$ and $GW$+EDMFT. In Section III, we describe the numerical methods and details used to solve the equations. In Section IV, we turn to the results. We first describe the phase diagram (subsection IV.1) within TRILEX and $GW$+EDMFT, and then focus on the weak-coupling regime (subsection IV.2) where, using the $GW$ method, we scan the space of the nearest and next-nearest-neighbor hopping parameters in search of dispersions with a weak antiferromagnetic instability where it is possible to reach a paramagnetic superconducting phase. The two dispersions which we thus identify are investigated in more detail at strong coupling with $GW$+EDMFT in subsections IV.3 and IV.4.

I Model

We solve the Hubbard model on the square lattice with longer range hoppings, defined by the Hamiltonian:

[TABLE]

with $i,j$ indexing lattice sites. $c_{\sigma i}^{\dagger}$($c_{\sigma i}$) denote creation (annihilation) operators, $n_{\sigma i}=c_{\sigma i}^{\dagger}c_{\sigma i}$ the density operator, $\mu$ the chemical potential, and $U$ the on-site Hubbard interaction. The hopping amplitudes, depicted on Fig. 1, are given by

[TABLE]

where $\mathbf{e}_{x,y}$ are the lattice vectors in the $x$ and $y$ directions. The bare dispersion is therefore

[TABLE]

When $t^{\prime}=t^{\prime\prime}=0$, the half-bandwidth is $D=4|t|$, but non zero $t^{\prime},t^{\prime\prime}$ in general make the bandwidth larger. Hereinafter, we express all quantities in units of $D$, unless stated differently.

II Formalism

The main goal of this paper is to study the superconducting (SC) phase of the two-dimensional Hubbard model within the TRILEX approach introduced in Refs. Ayral and Parcollet, 2015 and Ayral and Parcollet, 2016b. TRILEX is based on a bosonic decoupling of the interaction and a self-consistent approximation of the electron-boson vertex $\Lambda$ with a quantum impurity model. The decoupling of the on-site interaction is done by an exact Hubbard-Stratonovich transformation, leading to a model of non-interacting electrons coupled to some auxiliary bosonic modes representing charge and spin fluctuations.

We also study two methods which can be regarded as simplifications of the TRILEX method, namely $GW$+EDMFTSun and Kotliar (2002); Biermann et al. (2003); Sun and Kotliar (2004); Ayral et al. (2012, 2013); Biermann (2014) and $GW$Hedin (1965, 1999). In $GW$+EDMFT, vertex corrections are neglected in the non-local part of the self-energy and polarization. As both decay to zero, this additional approximation is negligible at very long distances. Due to the full treatment of the local vertex corrections, $GW$+EDMFT can capture the Mott transition, and we use it to obtain superconducting results in the doped Mott insulator regime. In the $GW$ method, vertex corrections are neglected altogether, and the self-energy and polarization are entirely calculated from bold bubble diagrams. The $GW$ equations do not require the solution of an auxiliary quantum impurity model, and are therefore less costly to solve. This additional approximation is justified only at weak coupling (see e.g. Ref. Ayral et al., 2012 for an illustration of its failure at large $U$), and there we use it to explore a large region of $(t^{\prime},t^{\prime\prime},T,n_{\sigma})$ parameter space ($T$ denotes temperature, $n_{\sigma}$ occupancy per spin).

Finally, let us stress that, in this paper, we use only single-site impurity models. Cluster extensions of TRILEX are discussed in our different work, Ref.Ayral et al., 2017. They naturally incorporate the effect of short-range antiferromagnetic exchange $J$ and give a quantitative control on the accuracy of the solution.

II.1 Superconducting Hedin equations

In this section, we derive the Hedin equationsHedin (1965, 1999); Aryasetiawan and Biermann (2008) which give the self-energy and polarization as functions of the three-leg vertex function. The derivation holds in superconducting phases and is relevant for fluctuations not only in the charge channelLinscheid and Essenberger , but also in the longitudinal and transversal spin channels.

II.1.1 The electron boson action

The starting point of the TRILEX method, as described in Ref. Ayral and Parcollet, 2016b, is the following electron-boson action:

[TABLE]

where $c^{*}_{\mu}$ and $c_{\nu}$ are Grassmann fields describing fermionic degrees of freedom, while $\phi_{\alpha}$ is a real bosonic field describing bosonic degrees of freedom. Indices $\mu,\nu$ stand for space, time, spin, and possibly other (e.g. band) indices $\mu\equiv(\mathbf{r}_{\mu},\tau_{\mu},\sigma_{\mu},\dots)$, where $\mathbf{r}_{\mu}$ denotes a site of the Bravais lattice, $\tau_{\mu}$ denotes imaginary time and $\sigma_{\mu}$ is a spin index ($\sigma_{\mu}\in\{\uparrow,\downarrow\}$). Indices $\alpha,\beta$ denote $\alpha\equiv(\mathbf{r}_{\alpha},\tau_{\alpha},I_{\alpha},\dots)$, where $I_{\alpha}$ indexes the bosonic channels. Repeated indices are summed over. Summation $\sum_{\mu}$ is shorthand for $\sum_{\mathbf{r}\in\mathrm{BL}}\sum_{\sigma}\int_{0}^{\beta}\mathrm{d}\tau$. $G_{0,\mu\nu}$ (resp. $W_{0,\alpha\beta}$) is the non-interacting fermionic (resp. bosonic) propagator.

Action (4) can result from the exact Hubbard-Stratonovich decoupling of the Hubbard interaction of Eq. (1) with bosonic fields $\phi$, but it can also simply describe an electron-phonon coupling problem.

In the present work, we are interested in a generalization of TRILEX able to accommodate superconducting order. To this purpose, we rederive the TRILEX equations starting from a more general action, written in terms of Nambu four-component spinors. The departure from the usual two-component Nambu-spinor formalism is necessary to allow for spin-flip electron-boson coupling in the action. Such terms do appear in the Heisenberg decoupling of the Hubbard interaction (see Section II.1.2).

We define a four-component Nambu-Grassmann spinor field as a column-vector

[TABLE]

where $i$ stands for the lattice-site $\mathbf{r}_{i}$. In combined indices, analogously to (4), a general electron-boson action can be written as

[TABLE]

where $u,v$ is a combined index $u\equiv(\mathbf{r}_{u},\tau_{u},a_{u},\dots)$, with $a,b,c,...\in\{0,1,2,3\}$ a Nambu index comprising the spin degree of freedom. The sum is redefined to go over all Nambu indices $\sum_{u}\equiv\sum_{\mathbf{r}\in\mathrm{BL}}\sum_{a}\int_{0}^{\beta}\mathrm{d}\tau$. Bold symbols are used for Nambu-index-dependent quantities.

This action does not depend on the conjugate field of $\boldsymbol{\Psi}$, because $\boldsymbol{\Psi}_{i}$ already contains all the degrees of freedom of the action (4) at the site $i$. The partition function corresponding to the bare fermionic part of the action has the following formZinn-Justin (2002)

[TABLE]

which is valid for any anti-symmetric matrix $\boldsymbol{A}$. Due to the unusual form of the action (no conjugated fields), the right-hand side is not the determinant of $\boldsymbol{A}$, but its square root, i.e. the Pfaffian. We can redefine the propagators/correlation functions of interest as

[TABLE]

The “conn” superscript denotes the connected part of the correlation function. The renormalized vertex is defined by

[TABLE]

Actions (6) and (4) are physically equivalent, namely their partition functions coincide:

[TABLE]

for an appropriate choice of $\boldsymbol{G}_{0}$ and $\boldsymbol{\lambda}$. Yet, they are not formally identical to each other, i.e. one cannot reconstruct (6) from (4) by mere relabeling $c\rightarrow\boldsymbol{\Psi}$, $\mu\nu\rightarrow uv$ (note the absence of Grassmann conjugation and the additional prefactors in the Nambu action). Therefore, one must rederive the Hedin equations which connect the self-energy and polarization with the full propagators $\boldsymbol{G}$ and $W$ and the renormalized vertex $\boldsymbol{\Lambda}$. We present the full derivation using equations of motion in Appendix A.2; here we just present the final result:

[TABLE]

Compared to the expressions in the normal case, there are extra factors $\frac{1}{2}$ in the Hartree term (second line in (13)) and polarization (13b). These factors come from the fact that with four-spinors, the summation over spin is performed twice. Note that the Hartree term can in principle have a frequency dependence if the bare electron-boson vertex has a dynamic part. On the other hand, the term beyond Hartree may as well contribute to the static part of the self-energy, if the bosonic propagator and the bare electron-boson vertex contain a static part. In all the calculations in this paper, the Hartree term is static and is the sole contributor the static part of self-energy. We will thus henceforth omit the Hartree term, as it can be absorbed in the chemical potential.

II.1.2 Connection to the Hubbard model

In this section, we specify the bare propagators and vertices such that action (6) corresponds to the Hubbard model Eq.(1). We then rewrite the Hedin equations under the assumption of spatial and temporal translational symmetry.

The Hubbard-Stratonovich transformation leading from Eq.(1) to an action of the form Eq.(4) relies on decomposing the Hubbard interaction as follows

[TABLE]

with $n_{I}\equiv\sum_{\sigma\sigma^{\prime}}c_{\sigma}^{\dagger}\sigma_{\sigma\sigma^{\prime}}^{I}c_{\sigma^{\prime}}$, and $I$ running within $\{0,z\}$ (“Ising decoupling”) or $\{0,x,y,z\}$ (“Heisenberg decoupling”) ($\sigma^{0}$ is the $2\times 2$ identity matrix, $\sigma^{x/y/z}$ are the usual Pauli matrices). This identity is verified, up to a density term, whenever

[TABLE]

in the Heisenberg decoupling. We have defined $U^{\mathrm{ch}}\equiv U^{0}$ and $U^{\mathrm{sp}}\equiv U^{x}=U^{y}=U^{z}$. Eqs (15a-15b) leave a degree of freedom in the choice of $U^{\mathrm{ch}}$ and $U^{\mathrm{sp}}$. Here, the choice $U^{x}=U^{y}=U^{z}$ stems from the isotropy of the Heisenberg decoupling (contrary to the Ising decoupling); it can describe SU(2) symmetry-broken phases. In the rest of the paper, we denote all quantities diagonal in the channel index with the channel as a superscript.

To make contact with the results of Ref. Kitatani et al., 2015, for $GW$ we will use the Ising decoupling with

[TABLE]

because the AF instabilities discussed in Sec. III.4, which violate the Mermin-Wagner theorem, are weaker in this scheme.

The equivalence of the action (6) with the Hubbard model is accomplished by setting

[TABLE]

The $4\times 4$ matrices are written in Nambu indices. The bare vertex reads:

[TABLE]

Thus, this vertex is local and static. The bare bosonic propagators are also local and static, as well as diagonal in the channel index:

[TABLE]

Our Hubbard lattice Nambu action reads (in explicit indices)

[TABLE]

II.1.3 Translational invariance, singlet pairing and SU(2) symmetry

In this paper, we restrict ourselves to phases with no breaking of translational invariance. With translational invariance in time and space, the propagators depend on frequency and momentum, and are matrices only in the Nambu index. We rewrite the Hedin equations derived above in the special case of the Hubbard action

[TABLE]

Similarly (see Appendix A.4 for details),

[TABLE]

Furthermore, we restrict ourselves to SU(2) symmetric phases, and allow only for singlet pairing, therefore

[TABLE]

We allow no emergent mixing of spin

[TABLE]

These assumptions simplify the structure of the Green’s function in Nambu space

[TABLE]

where the normal and anomalous Green’s functions read

[TABLE]

Under the present assumptions $G_{\mathbf{k}}(\tau)$ is real, therefore $G_{\mathbf{k}}(-i\omega)=G_{\mathbf{k}}^{*}(i\omega)$. Here note that SU(2) symmetry and lattice inversion symmetry imply $F_{ij}(\tau)=F_{ij}(-\tau)=F_{ji}(\tau)$ (this can be proven by rotating $c_{\sigma}\rightarrow(-)^{\delta_{\uparrow,\sigma}}c_{\bar{\sigma}}$). Therefore, if $F_{ij}(\tau)$ is real, $F_{\mathbf{k}}(i\omega)$ is also purely real. In this paper we consider only purely real $F_{ij}(\tau)$.

Similarly, the block structure of the self-energy is given by:

[TABLE]

$\Sigma$ and $S$ are the normal and anomalous self-energies defined by the Nambu-Dyson equation

[TABLE]

where the inverse is assumed to be the matrix inverse in Nambu indices. Component-wise, under the present assumptions, the Nambu-Dyson equation reads

[TABLE]

Furthermore, due to SU(2) symmetry, the full bosonic propagator will be identical in the $x$, $y$ and $z$ channels, so we define

[TABLE]

and have $W^{x}=W^{y}=W^{z}=W^{\mathrm{sp}}$, and similarly for the renormalized vertex. This will simplify the calculation of the self-energy in the Heisenberg decoupling scheme, as the contribution coming from $x$ and $y$ bosons is the same as the one coming from the $z$ boson. The bosonic Dyson equation is then always solved in only two channels

[TABLE]

II.2 TRILEX, $GW$+EDMFT and $GW$ equations

II.2.1 Single-site TRILEX approximation for $d$-wave superconductivity

The single-site TRILEX method consists in approximating the renormalized vertex by a local quantity, obtained from an effective single-site impurity model

[TABLE]

Solving the TRILEX equations amounts to finding $\boldsymbol{\mathcal{G}}(i\omega)$ and $\mathcal{U}(i\Omega)$ such that the full propagators in the effective impurity problem (II.2.1) coincide with the local components of the ones obtained on the lattice, namely, we want to satisfy

[TABLE]

where the vertex of Eq. (21) is approximated by the impurity vertex:

[TABLE]

In this paper, we allow only strictly $d$-wave superconducting pairing. Thus

[TABLE]

which means that the anomalous components of the local Green’s function $\boldsymbol{G}_{\mathrm{loc}}$ will be zero. Therefore, at self-consistency (Eq. (42a)), the impurity’s Green’s function is normal and thus the anomalous components of the bare propagator on the impurity must vanish

[TABLE]

This means that the impurity problem will be identical to the one in the normal-phase calculations, which can be expressed in terms of the original Grassmann fields

[TABLE]

where the bare vertices (slim symbols denote the impurity quantities) are given by Pauli matrices $\lambda^{I}_{\sigma\sigma^{\prime}}=\sigma^{I}_{\sigma\sigma^{\prime}}$. After integrating out the bosonic degrees of freedom, one obtains an electron-electron action with retarded interactions

[TABLE]

This single-site impurity problem is solved using the numerically exact hybridization-expansion continuous-time quantum Monte Carlo (CTHYB or HYB-CTQMCWerner et al. (2006); Werner and Millis (2007)), employing the segment algorithm. The transverse spin-spin interaction term is dealt with in an interaction-expansion mannerOtsuki (2013). See Ref. Ayral and Parcollet, 2016b for details.

Under the present assumptions, the approximation for the renormalized vertex entering the Hedin equations Eq. (21) is

[TABLE]

where $\circ$ denotes the elementwise product $[A\circ B]_{ij}=A_{ij}B_{ij}$ (see Appendix A.5 for details).

We obtain $\Lambda_{\mathrm{imp}}^{\eta}$ from the three-point correlation function on the impurity using

[TABLE]

where

[TABLE]

and

[TABLE]

[TABLE]

We can now write the final expressions for the self-energy and polarization:

[TABLE]

with $p_{\mathrm{ch}}=1$, $p_{\mathrm{sp}}=0$, $m_{\mathrm{ch}}=1$. These equations hold in both the Heisenberg ($m_{\mathrm{sp}}=3$) and Ising ($m_{\mathrm{sp}}=1$) decoupling schemes. In the expression for the polarization (Eq. (61c)) we have used lattice inversion symmetry and the symmetries of $\Lambda$ and $\boldsymbol{G}$. Under the present assumptions, $P$ is purely real (see Appendix A.3 for details).

II.2.2 $GW$+EDMFT

The $GW$+EDMFT approximation can be regarded as a simplified version of TRILEX where, in the calculation of the non-local ($\mathbf{r}\neq 0$) part of self-energy and polarization (second line of Eqs. (64a),(64b) and (64c) below), an additional approximation is made:

[TABLE]

The efficiency is gained because one need not measure the three-point correlator $\tilde{\chi}^{3,\eta,\mathrm{conn}}$ in the impurity model. The local self-energy and polarization still have vertex-corrections, but are identical to $\Sigma$ and $P$ on the impurity, which can be computed from only two-point correlators. Furthermore, the calculation of the non-local parts of the self-energy and polarization can now be performed in imaginary time, as opposed to the explicit summation over frequency needed in Eqs. (64a), (64b) and (64c).

II.2.3 $GW$

If we approximate the renormalized vertex by unity even in the calculation of the local part of self-energies, we obtain an approximation similar to the $GW$ approximation, with the important difference that we are using a decoupling in both charge and spin channels, unlike the conventional $GW$ approaches which are limited to the charge channel. This additional approximation eliminates the need for solving an impurity problem, as now even the local self-energy and polarization are calculated by the bubble diagrams Eq. (61a), (61b) and (61c), simplified by Eq. (62).

To summarize, the exact expressions for the self-energy and boson polarization are compared to the approximate ones in $GW$, EDMFT, $GW$+EDMFT, and TRILEX in Fig. 2.

II.2.4 Normal phase calculation

In the normal phase, the further simplification is that $F_{\mathbf{k}}(i\omega)=0$. Therefore, $S_{\mathbf{k}}(i\omega)=0$ and the Dyson equation (38a) reduces to the familiar form

[TABLE]

III Methods

III.1 Numerical implementation of the Hedin equations

As shown in Ref. Ayral and Parcollet, 2016b, it is numerically advantageous to perform the computation in real space and to split the self-energy and polarization in the following way:

[TABLE]

where $\tilde{X}_{\mathbf{r}}(i\omega)\equiv(1-\delta_{\mathbf{r}})X_{\mathbf{r}}(i\omega)$. In the presence of lattice inversion symmetry, $X_{\mathbf{r}}=X_{-\mathbf{r}}$. The impurity’s self-energy and polarization are defined as

[TABLE]

III.2 Solution by forward recursion

In practice, the TRILEX, $GW$+EDMFT and $GW$ equations can be solved by forward recursion:

Start with a given $\Sigma_{\mathbf{k}}(i\omega)$ and $P^{\eta}_{\mathbf{q}}(i\Omega)$, and (for SC phase only) $S_{\mathbf{k}}(i\omega)$ and (for TRILEX and $GW$+EDMFT only) $\Sigma_{\mathrm{imp}}(i\omega)$ and $P^{\eta}_{\mathrm{imp}}(i\Omega)$ (for instance set them to zero, or use EDMFT results) 2. 2.

Compute the new $G_{\mathbf{k}}(i\omega)$ and $W_{\mathbf{q}}^{\eta}(i\Omega)$ and (for SC phase only) $F_{\mathbf{k}}(i\omega)$ from Eqs. (38a, 40, 38b). 3. 3.

(TRILEX/$GW$+EDMFT only) Impose the self-consistency condition Eq. (42a, 42b) by reversing the impurity Dyson equations (65a, 65b), such that

[TABLE] 4. 4.

(TRILEX/$GW$+EDMFT only) Solve the impurity model with the above bare fermionic and bosonic propagators: compute $G_{\mathrm{imp}}$, $\chi_{\mathrm{imp}}^{\eta}$, $\langle n^{\eta}\rangle_{\mathrm{imp}}$ and (for TRILEX only) $\tilde{\chi}^{3,\eta,\mathrm{conn}}$ and from them $\Sigma_{\mathrm{imp}}$ (Eq. 65a), $P_{\mathrm{imp}}^{\eta}$ (Eq. 65b) and (TRILEX only) $\Lambda_{\mathrm{imp}}^{\eta}$ (Eq. 54); 5. 5.

Compute $\Sigma_{\mathbf{k}}(i\omega)$ and $P_{\mathbf{q}}^{\eta}(i\Omega)$ and (for SC phase only) $S_{\mathbf{k}}(i\omega)$ with Eqs. (64a, 64c, 64b); 6. 6.

Go back to step 2 until convergence is reached.

III.3 Superconducting temperature $T_{c}$

In order to determine the superconducting transition temperature $T_{c}$, we solve a linearized gap equation (LGE). At $T=T_{c}$, the anomalous part of the self-energy $S$ vanishes. Linearizing Eq. (38b) with respect to $S$ and plugging it into Eq. (64b) leads to an implicit equation for $T_{c}$, featuring only the normal component of the Green’s function

[TABLE]

Using four-vector notation $k\equiv(\mathbf{k},i\omega)$, we obtain

[TABLE]

This is an eigenvalue problem for $S$. In practice, it is more convenient to consider the spectrum of the operator $A$,

[TABLE]

The eigenvalues $\lambda$ and the eigenvectors $S_{k}^{\lambda}$ depend on the temperature $T$. The critical temperature $T_{c}$ is therefore given by

[TABLE]

where $\lambda_{m}$ is the largest eigenvalue of $A$. In other words, $T=T_{c}$ when the first eigenvalue crosses 1. In addition, the symmetry of the superconducting instability is given by the $k$ dependence of $S$ for the corresponding eigenvector.

In practice, we first solve the normal-phase equations, and then solve the LGE Eq. (III.3) by forward substitution. Starting from an initial simple $d_{x^{2}-y^{2}}$-wave form

[TABLE]

we use the power method Mises and Pollaczek-Geiringer (1929) to compute the leading eigenvalue of the operator $A$. We do this in a select range of temperature for the given parameters $(U,n,t,t^{\prime},t^{\prime\prime})$ and monitor the leading eigenvalue $\lambda_{m}(T)$. If we observe a $T_{c}$ ($\lambda_{m}(T)>1)$), we can then use the eigenvector $S$ as an initial guess to stabilize the superconducting solution using the algorithm from section III.2. We have also examined other irreducible representations of the symmetry group and found that this $d$-wave representation is the one with highest $T_{c}$, in agreement with Refs. Otsuki et al., 2014; Rømer et al., 2015.

III.4 The AF instability

As documented in Refs. Ayral and Parcollet, 2015, 2016b, the TRILEX equations present an instability towards antiferromagnetism below some temperature $T_{\mathrm{AF}}$ (see also Refs Kitatani et al., 2015; Otsuki et al., 2014). The antiferromagnetic susceptibility $\chi^{\mathrm{sp}}$ is related to the propagator of the boson in the spin channel via

[TABLE]

They both diverge at $T=T_{\mathrm{AF}}$ because the polarization becomes too large (the denominator in (40) vanishes). This instability, which is an artifact of the approximation for the two-dimensional Hubbard model, violates the Mermin-Wagner theorem. For many values of $t^{\prime},t^{\prime\prime}$, this AF instability prevents us from reaching the superconducting temperature $T_{c}$.

This AF instability also exists in conventional cluster DMFT methods (cellular DMFT, DCA) Staar et al. (2014); Maier et al. (2005c); Capone and Kotliar (2006). Yet, in most works, it is simply ignored by enforcing a paramagnetic solution (by symmetrizing up and down spin components). In TRILEX, however, we do not have this possibility. Indeed, the antiferromagnetic susceptibility directly enters the equations (via $W$), and its divergence makes it impossible to stabilize a paramagnetic solution of the TRILEX equations at a temperature lower than $T_{\mathrm{AF}}$. For a precise definition of $T_{\mathrm{AF}}$ in the present context, see Appendix C.

In the following, we circumvent this issue in two ways: either by extrapolating the temperature dependence of the eigenvalue of the linearized gap equation to low temperatures, despite the AF instability (section IV.1, with tight-binding values $t^{\prime},t^{\prime\prime}$ relevant for cuprate physics), or, in section IV.2, by finding other values of $t^{\prime},t^{\prime\prime}$, where the Fermi surface shape is qualitatively similar to the cuprate case, but where the AF instability occurs at a temperature lower than $T_{c}$.

IV Results and discussion

IV.1 Phase diagram

First, using the linearized-gap equation (LGE) method described in Sect. III.3, we compute the SC phase boundary from high temperature, for $t^{\prime}=-0.2t,t^{\prime\prime}=0$, a physically relevant case for the physics of cuprates. We set $U/D=4$ in order to be above the Mott transition threshold at half filling (we recall that for the square lattice, $U_{c}/D\approx 2.4$ within single-site DMFTSchäfer et al. (2015)). The results are presented on Fig. 3.

The top panel presents the largest eigenvalue of the LGE as a function of temperature, for TRILEX and $GW$+EDMFT. The calculation becomes unstable due to AF instability before we can observe $\lambda_{m}>1$. The extrapolation of $\lambda_{m}$ towards low temperature is not straightforward. We use an empirical law

[TABLE]

to fit the data and extrapolate to lower temperature. This form can be shown (see Appendix C) to provide a very good fit to similar computations in the DCA and DCA*+* methods, from the data of Refs. Maier et al., 2006b; Staar et al., 2014. We perform the fit and extrapolation with $\gamma=0.3$ for $GW$+EDMFT and $\gamma=0.45$ for TRILEX, and get the result for $T_{c}$ reported with solid lines on the bottom panel. The error bars shown are obtained by fitting and extrapolating with $\gamma$ varied in the window 0.3-0.6. The error bars coming from the uncertainty of the fit for a fixed $\gamma$ and a detailed discussion of the fitting procedure can be found in Appendix C. The dashed lines denote the temperature of the antiferromagnetic instability, below which no stable paramagnetic calculation can be made.

For all values of $\gamma$, the raw data at high temperature for both methods indicate a similar dome shape for $T_{c}$ vs $\delta$, where $\delta$ is the percentage of hole-doping, $\delta[\%]=(1-2n_{\sigma})\times 100$ ($n_{\sigma}=\frac{1}{2}$ corresponds to half-filling). The fact that $T_{c}$ vanishes at zero $\delta$ can be checked directly, but we cannot exclude that it vanishes at a finite, small value of $\delta$. The optimal doping in both methods is found to be around 12%. At half-filling, both methods recover a Mott insulating state, and $\lambda_{m}(T)$ is found to be very small. We observe that TRILEX has a higher $T_{c}$ than $GW$+EDMFT, showing that the effects of the renormalization of the electron-boson vertex are non-negligible in this regime.

These results for $T_{c}(\delta)$ are qualitatively comparable to the results of cluster DMFT methods, e.g. the 4-site CDMFT + ED computation of Refs Capone and Kotliar, 2006; Sakai et al., 2016; Kancharla et al., 2008, or the 8-site DCA results of Ref. Gull and Millis, 2015. In particular, Ref. Sakai et al., 2016 reports a $T_{c}/D\approx 0.0125$ at doping $\delta=13\%$ in a doped Mott insulator, which falls half-way between the TRILEX and $GW$+EDMFT results. Furthermore, the optimal doping in Ref. Capone and Kotliar, 2006 seems to coincide with our result, while in Ref. Sakai et al., 2016 it is somewhat bigger (around $20\%$). We emphasize however that here we solve only a single-site quantum impurity problem, and obtain the $d$-wave order, which is not possible in single-site DMFT due to symmetry reasons.

Let us now turn to the weak-coupling regime ($U/D=1$). We present in Fig. 4 the SC temperature in the $GW$ and $GW$+EDMFT approximation within the Ising decoupling (for the $\lambda(T)$ plot, see Appendix C). Both methods give similar results, which justifies using the faster $GW$ at weak coupling. In contrast to the larger-$U$ case, one does not obtain the dome versus doping due to the absence of Mott insulator at $\delta=0$.

We compare our results with the order parameter at $T=0$ obtained from a $2\times 2$ CDMFT+ED calculationCapone and Kotliar (2006). The general trend observed is similar: optimal doping is zero, and there is a quick reduction of $T_{c}$ between 12 and 16% doping.

As for the value of $T_{c}$, we compare to the result presented in Ref. Staar et al., 2014. Here, a DCA*+* calculation with a 52-site cluster impurity, at $U/D=1$,$t^{\prime}=t^{\prime\prime}=0$, $\delta=10\%$, predicts $T_{c}/D\approx 0.06$. With the same parameters, $GW$ gives $T_{c}/D\approx 0.21$, $GW$+EDMFT gives $T_{c}/D\approx 0.27$, hence overestimating $T_{c}$.

IV.2 Weak coupling

As explained in Sec. III.4, in order to study the SC phase itself, we need to identify a dispersion for which $T_{c}$ is above $T_{\mathrm{AF}}$. To achieve this, we first scan a large set of parameters $t^{\prime},t^{\prime\prime}$ with the $GW$ approximation at weak coupling. Indeed, at weak coupling, we can approximate TRILEX by $GW$, which is faster to compute (there is no quantum impurity model to solve). We look for a $(t^{\prime},t^{\prime\prime})$ point for which not only $T_{\mathrm{AF}}>T_{c}$, but also the shape of the Fermi surface is qualitatively compatible with cuprates. We find a whole region of parameters where this is satisfied, and then use these parameters in a strong-coupling computation with $GW$+EDMFT and TRILEX. Whether a weak coupling computation is a reliable guide in the search for $t^{\prime},t^{\prime\prime}$ with maximal $T_{c}$ at strong coupling remains open and would require a systematic exploration with cluster methods. However, at least in one example (shown below), this assumption will provide us with an appropriate choice of hopping amplitudes that allows us to stabilize a superconducting solution in the doped Mott insulator regime.

Fig. 5 presents the computation of the AF instability ($T_{\mathrm{AF}}$) and the SC instability ($T_{c}$) in $GW$, for $U/D=1$ and various $t^{\prime},t^{\prime\prime}$ ($t=-1.0$ is held fixed) and various dopings. The temperature is taken from $0.2$ down to the lowest accessible temperature, but not below $0.01$ in cases where the extrapolation of $\lambda(T)$ yielded no finite $T_{c}$. The temperature step depends on $T$ (smaller step at lower $T$; see Appendix C for an example of raw data).

The first observation is that the region of high $T_{c}$ broadly coincides with the region of high $T_{\mathrm{AF}}$. This is expected as in $GW$ the attractive interaction comes from the spin-boson, and a high-valued and sharply-peaked $W^{\mathrm{sp}}$ is clearly necessary for satisfying the gap equation Eq. (III.3) with $\lambda=1$. However, the maximum of $T_{c}$ with respect to $(t^{\prime},t^{\prime\prime})$ at a fixed $n$ does not coincide with the maximum of $T_{\mathrm{AF}}$, thus indicating that there are factors other than sharpness (criticality) of the spin-boson which contribute to the height of $T_{c}$. While the maximum of $T_{\mathrm{AF}}$ is found rather close to $t^{\prime}=t^{\prime\prime}=0$ at all dopings, the maximum in $T_{c}$ starts from $(t^{\prime},t^{\prime\prime})=(-0.6,-0.4)$ at $n=0.38$ and gradually moves as $n$ is increased. It is only at half-filling that the two maxima are found to coincide. Furthermore, while at around $t^{\prime}=t^{\prime\prime}=0$ and $t^{\prime}\approx t^{\prime\prime}$ one sees $T_{\mathrm{AF}}>T_{c}$, this trend is gradually reversed as $t^{\prime\prime}$ is made more and more negative, such that around $t^{\prime}\approx t^{\prime\prime}+0.4$ one usually sees a finite $T_{c}$ in the absence of a finite $T_{\mathrm{AF}}$.

In Fig. 6, we plot $T_{\mathrm{AF}}$ and $T_{c}$ vs doping for different values of $t^{\prime},t^{\prime\prime}$. The corresponding dispersion (color map) and Fermi surfaces (gray contours; red for the maximal $T_{c}$) are presented in the insets.

Finally, in Fig.7, we summarize the observations from Fig. 5. The blue dot denotes the global maximum of $T_{c}$ and $T_{\mathrm{AF}}$. The dashed gray lines denote the directions of the slowest and quickest decay of antiferromagnetism. The red ellipses denote the regions of maximal $T_{c}$, at various dopings. The yellow region is where one finds little antiferromagnetism, but still a sizable $T_{c}$. The green region corresponds to dispersions relevant for cupratesPavarini et al. (2001). The points A,B, and C are the dispersions that we focus on and for which we perform TRILEX and $GW$+EDMFT computations. Pt. B is most relevant for the cuprates, and was analyzed in Fig.3. Pt. C has $T_{\mathrm{AF}}<T_{c}$ which allows us to converge a superconducting solution at both weak and strong coupling. We analyze it in the next subsection. Pt. A is where we observe an maximal $T_{c}$ at $16\%$ doping, and we focus on it in Section IV.4.

IV.3 The nature of the superconducting phase at strong coupling

In this section, we study the dispersion C $(t,t^{\prime},t^{\prime\prime})=(-1,-0.3,-0.6)$. In Figure 6, we have determined that at weak coupling ($U/D=1$), the superconducting temperature $T_{c}$ is larger than the AF temperature: we can therefore reach the superconducting phase numerically (see Appendix D). It turns out that at strong coupling, the AF instability is also absent. This allows us to stabilize superconducting solutions in the doped Mott insulator regime. We also perform a calculation restricted to the normal phase for all parameters in order to compare results to the ones in the SC phase. For simplicity, in this section we will present only $GW$+EDMFT results for $U/D=4$.

In Fig. 8, we show the superconducting temperatures at $U/D=1$ and $U/D=4$. Contrary to pt.B, in pt.C strong coupling seems to strongly enhance superconductivity. Also, the SC dome extends to higher dopings.

In Fig. 9 we show the results for the both the anomalous self-energy and Green’s function, as well as the imaginary part of the normal self-energy, in both the normal phase and superconducting solution, anti-nodal and nodal regions.

The imaginary part of the normal self-energy is larger at antinodes than at nodes and is growing when approaching the Mott insulator. When going from the normal phase to the SC phase, the imaginary part of the self-energy is strongly reduced at the antinode and weakly reduced at the node. The difference between the normal and SC solution (light blue area) is roughly proportional to the anomalous self-energy in the SC phase (blue line). Note that we observe a similar phenomenon even at weak coupling (see Appendix D).

In Fig. 10, we plot the spectral function at the antinodes at low temperature, in the normal and in the superconducting phase. At low doping, we observe at low energy a pseudo-gap in the normal phase and the superconducting gap in the SC phase. The result obtained here is qualitatively different to the one obtained using 8 sites DCA cluster by Gull et al. Gull et al. (2013b); Gull and Millis (2015). In the cluster computations, the superconducting gap is smaller than the pseudogap, i.e. the quasi-particle peak at the edge of the SC gap appears within the pseudogap. It is not the case here. Also, we do not see any “peak-dip-hump” structure. Note that we are however using different parameters (for the hoppings $t^{\prime},t^{\prime\prime}$, the interaction $U$ and the doping $\delta$). It is not clear at this stage whether these qualitative differences are due to this different parameter regime or to an artifact of the single-site TRILEX method, e.g. the lack of local singlet physics in a single-site impurity model. Further investigations with cluster-TRILEX methods are necessary in the SC phase.

In Fig. 11, we plot various quantities at the lowest Matsubara frequency, as a function of $\mathbf{k}$. In the first two rows we compare the anomalous self-energy and the pairing amplitude. Both are clearly of $d$-wave symmetry. The pairing amplitude has a different order of magnitude (see Appendix A.6 for an illustration of the dependence between $F$,$G$,$\Sigma$ and $S$). In the third and fourth row we show the imaginary part of the Green’s function in the SC and normal phase. Due to the absence of long-lived quasiparticles in this sector, the maximum of $F_{\mathbf{k}}$ is moved towards the nodes, and does not coincide with the maximum of $S_{\mathbf{k}}$. At small doping, the Fermi surface in both cases becomes less sharp and more featureless, due to proximity to the Mott insulator. In the next two rows we show the imaginary part of the normal self-energy. In the superconducting phase, $\mathrm{Im}\Sigma_{\mathbf{k}}$ is strongly reduced in only anti-nodal regions, and thus flattened (made more local). In the last row, we show the non-local part of the propagator for the spin boson. At large doping we observe a splitting of resonance at $(\pi,\pi)$ which corresponds to incommensurate AF correlations (see e.g. Ref. Onufrieva and Pfeuty, 2002 for a similar phenomenon). Having that the Green’s function at around $\mathbf{k}=(0,0)$ is quite featureless, and that the boson is sharply peaked at zero frequency, the shape of the spin-boson around $\mathbf{q}=(\pi,\pi)$ is similar to the self-energy at around $\mathbf{k}=(\pi,\pi)$. This pattern is observed at all three dopings.

IV.4 Strong-coupling $T_{c}$ at pt.A

At weak coupling, we have observed in section IV.2 that the dispersion pt.A ($(t,t^{\prime},t^{\prime\prime})=(-1,-0.5,-0.2)$) presents a pronounced maximum in $T_{c}(t^{\prime},t^{\prime\prime})$ at $16\%$ doping. Here, we investigate that point at strong coupling using $GW$+EDMFT and TRILEX and find that also at $U/D=4$, the $T_{c}$ is substantially higher than in pt.B and pt.C. Here $T_{c}$ is below $T_{\mathrm{AF}}$ and the result is again based on extrapolation of $\lambda$. The proposed fitting function in this case does not perform as well and the extrapolation is less reliable, but $GW$+EDMFT and TRILEX are in better agreement than in the case of pt.B. A further investigation using cluster methods is necessary since, apart from Refs. Kancharla et al., 2008; Chen et al., 2013; Rømer et al., 2015, little systematic exploration of $T_{c}(t^{\prime},t^{\prime\prime})$ has been performed.

V Conclusion

In this work, we have generalized the TRILEX equations and their simplifications $GW$+EDMFT and $GW$ to the case of paramagnetic superconducting phases, using the Nambu formalism. We also generalized the corresponding Hedin equations. We have then investigated within TRILEX, $GW$+EDMFT and $GW$ the doping-temperature phase diagram of the two-dimensional single-band Hubbard model with various choices of hopping parameters. In the case of a bare dispersion relevant for cuprates, in the doped Mott insulator regime, both TRILEX and $GW$+EDMFT yield a superconducting dome of $d_{x^{2}-y^{2}}$-wave symmetry, in qualitative agreement with earlier cluster DMFT calculations. Let us emphasize that this was obtained at the low cost of solving a single-site impurity model. At weak coupling, we have performed a systematic scan of tight-binding parameter space within the $GW$ approximation. We have identified the region of parameter space where superconductivity emerges at temperatures higher than antiferromagnetism. With one of those dispersions, we studied the properties of the superconducting phase at strong coupling with $GW$+EDMFT. We also addressed the question of the optimal dispersion for superconductivity in the Hubbard model at weak coupling. At $16\%$ doping we identify a candidate dispersion for the highest $d$-wave $T_{c}$, which remains to be investigated in detail at strong coupling (e.g. with cluster DMFT methods).

The next step will be to solve in the SC phase the recently developed cluster TRILEX methods Ayral et al. (2017). Indeed, the single-site TRILEX method contains essentially an Eliashberg-like equation with a decoupling boson, and a local vertex (computed from the self-consistent impurity model) which has no anomalous components. The importance of anomalous vertex components and the effect of local singlet physics (present in cluster methods) is an important open question. Note that the framework developed in this paper can also be used to study more general pairings and decoupling schemes in TRILEX, e.g. the effect of bosonic fluctuations in the particle-particle (i.e. superconducting) channel.

Finally, let us emphasize that the question of superconductivity in multi-orbital systems like iron-based superconductors is another natural application of the TRILEX method, in particular in view of the strong AF fluctuations in these compounds. In this multi-orbital case, being able to describe the SC phase without having to solve clusters (which are numerically very expensive within multi-orbital cluster DMFTNomura et al. (2015); Sémon et al. (2017)) could prove to be very valuable.

Acknowledgements.

We thank M. Kitatani for useful insights and discussion. This work is supported by the FP7/ERC, under Grant Agreement No. 278472-MottMetals. Part of this work was performed using HPC resources from GENCI-TGCC (Grant No. 2016-t2016056112). The CT-HYB algorithm has been implemented using the TRIQS toolboxParcollet et al. (2015).

Appendix A Details of derivations

A.1 Relation between $\boldsymbol{\chi}^{3}$ and $\tilde{\boldsymbol{\chi}}^{3}$

Let us define the following correlation functions:

[TABLE]

In this section, we derive useful relations between these quantities.

Let us introduce source fields in the electron-boson action, Eq. 6:

[TABLE]

We may now write

[TABLE]

Let us now integrate out the bosonic degrees of freedom in Eq. 74. We obtain:

[TABLE]

with

[TABLE]

We now perform the derivatives of Eqs. 75 and 76 using the new expression Eq. 77, yielding:

[TABLE]

Thus, we have, for the full correlator, as well as for the connected and disconnected parts:

[TABLE]

A.2 Derivation of Hedin equations from Equations of motion

In this section, we derive the Hedin equations of the main text using the Dyson-Schwinger equation of motion techniqueZinn-Justin (2002) already used in Ref. Ayral and Parcollet, 2016b.

A.2.1 E.O.M. for the self-energy

Since the functional integral of a total derivative vanishes:

[TABLE]

for any $f$ and $g$, we have

[TABLE]

which comes directly from the Leibniz derivation rule for Grassmann variables. $\deg f$ denotes the degree of the polynomial $f$ in the variable $\boldsymbol{\Psi}$. Let us now assume $f[\Psi]=e^{-S_{0}[\boldsymbol{\Psi}]}=e^{\frac{1}{2}\boldsymbol{\Psi}_{u}\boldsymbol{G}_{0,uv}^{-1}\boldsymbol{\Psi}_{v}}$ and $g[\boldsymbol{\Psi}]=h[\Psi]e^{-V[\boldsymbol{\Psi}]}$, with $h$ containing an odd number of Grassmann fields. $f$ has an infinite number of terms but all are products of an even number of $\boldsymbol{\Psi}$ fields. We obtain

[TABLE]

On the l.h.s. we have again used the Leibniz rule with $\deg h$ assumed to be odd, hence the extra minus sign. On the r.h.s similarly, $\deg\boldsymbol{\Psi}=1$, and $\boldsymbol{G}_{0,uv}^{-1}=-\boldsymbol{G}_{0,vu}^{-1}$, so the $\frac{1}{2}$ prefactor is canceled. Both integrals are now averages with respect to the action $S=S_{0}+V$, namely

[TABLE]

Let us now consider the case when $h\equiv\boldsymbol{\Psi}_{v}$, and $V$ is the interacting part of the electron-electron action (A.1), with the source field $H$ set to zero, i.e. $V\equiv\frac{1}{8}\left[W_{0}\right]_{\alpha\beta}\left(\boldsymbol{\Psi}_{u}\boldsymbol{\lambda}_{uw\alpha}\boldsymbol{\Psi}_{w}\right)\left(\boldsymbol{\Psi}_{y}\boldsymbol{\lambda}_{yz\beta}\boldsymbol{\Psi}_{z}\right)$. We get

[TABLE]

Multiplying both sides by $\boldsymbol{G}_{0}$ and using Eqs. (73a) and (81):

[TABLE]

Since the self-energy is defined as

[TABLE]

we obtain

[TABLE]

The second term is the Hartree term (note the 1/2 factor). The Fock term is included in the first term.

A.2.2 E.O.M. for the polarization

Real fields $\phi$ commute with the derivative, so the Leibniz rule is simpler. Analogously to Eq. (83)

[TABLE]

Similarly to Eq. (84), by taking $f[\phi,\boldsymbol{\Psi}]=e^{-S_{0}[\boldsymbol{\Psi},\phi]}$, where $S_{0}$ is the non-interacting part of the electron-boson action (4), and $V[\boldsymbol{\Psi},\phi]=\frac{1}{2}\boldsymbol{\Psi}_{u}\boldsymbol{\lambda}_{uv\delta}\boldsymbol{\Psi}_{v}\phi_{\delta}$, one has

[TABLE]

Again, note the minus sign in the left-hand side (to be compared with Eq. (84)) coming from the bosonic nature of the field $\phi$. For $h\equiv\phi_{\alpha}-\langle\phi_{\alpha}\rangle$,

[TABLE]

Multiplying by $W_{0}$ and using Eqs. 10 and 11, we obtain

[TABLE]

With the definition of $P$ as

[TABLE]

we identify

[TABLE]

Note the extra prefactor $\frac{1}{2}$ compared to the normal-case expression.

A.3 Proof that $P$ is real

In the derivation of Eq. (61c) we have used the symmetries of $G$,$F$ and $\Lambda$. It turns out that the imaginary part of $\Lambda$ does not play a role in the summation and that the polarization is strictly real.

The renormalized vertex has the following symmetriesAyral and Parcollet (2016b)

[TABLE]

Under the present assumptions, all components of the Green’s function ($G$ and $F$) have the property

[TABLE]

[TABLE]

Therefore

[TABLE]

which proves that the polarization is real. In the derivation of the first term in Eq. (61c), we have used the equality between Eq. (95) and Eq. (96). Then, for any real-valued $F$, we furthermore have $F_{\mathbf{k}}(i\omega)=F_{\mathbf{k}}(-i\omega)$, which gives us

[TABLE]

which is what we use in the derivation of the second term in Eq. (61c).

A.4 Fourier transforms: Hedin equations with translational symmetry

Here, we derive Eq. (22). A completely analogous derivation can be used for Eqs. (21).

For the sake of clarity, we omit the spatial indices, as the spatial Fourier transform (FT) is completely analogous to the temporal FT.

[TABLE]

Applying the (implicit) integration over times produces Kronecker delta functions at $\omega^{\prime\prime}=\omega$, $\omega^{\prime}=\omega+\Omega$ and $\Omega=\Omega^{\prime}$. Therefore

[TABLE]

We now reinstate the momentum indices, and obtain Eq. (22) by identifying the summands on both sides of the equation

[TABLE]

Here, summation over $c,d$ is implicit.

A.5 $\boldsymbol{\Lambda}_{\mathrm{imp}}$ from $\Lambda_{\mathrm{imp}}$

Here we prove Eq. (48). In the Hubbard model we’ll have

[TABLE]

On the impurity (II.2.1), where we have no anomalous components

[TABLE]

The $\frac{1}{2}$ prefactor in (81) cancels the prefactor $2$ in (101). If we define

[TABLE]

[TABLE]

we can rewrite

[TABLE]

More compactly

[TABLE]

where $\boldsymbol{\lambda}^{I}$ and $\circ$ have been defined in main text. For $I=0,z$, we have used

[TABLE]

and

[TABLE]

and similar considerations for $I=x,y$. Expressions completely analogous to (103e) and (103j) hold for the connected part of $\boldsymbol{\chi}^{3}$. Plugging these in Eq. (22) together with Eq. (81),

[TABLE]

and

[TABLE]

immediately yields Eq. (48). Eq. (111) holds in presence of SU(2) symmetry. It can be proven by applying a $\pi/2$ rotation around the $y$ axis ($n^{z}\rightarrow-n^{x}$,$n^{x}\rightarrow n^{z}$, $n^{y}\rightarrow n^{y}$), i.e $c_{\sigma}\rightarrow[\exp(-\frac{i}{2}\frac{\pi}{2}\sigma^{y})]_{\sigma,\sigma^{\prime}}c_{\sigma^{\prime}}=\frac{1}{\sqrt{2}}(c_{\sigma}+(-)^{\delta_{\sigma,\uparrow}}c_{\bar{\sigma}})$:

[TABLE]

and then rotating the operators of the first term on the r.h.s. by $\pi$ around the $y$ axis ($c_{\sigma}\rightarrow[\exp(-\frac{i}{2}\pi\sigma^{y})]_{\sigma,\sigma^{\prime}}c_{\sigma^{\prime}}=(-)^{\delta_{\uparrow,\sigma}}c_{\bar{\sigma}}$).

A.6 Relation between $S$, $F$, $\Sigma$ and $G$

Here we emphasize that the order of magnitude of the anomalous self-energy $S$ and that of the pairing amplitude $F$ are not the same, as illustrated on Fig. 13. The pairing amplitude has a strongly non-monotonous dependence on the anomalous self-energy. At a given normal self-energy, there is a “sweet spot” where a small anomalous self-energy produces a very strong superconducting pairing. As soon as the anomalous self-energy starts gapping out the Green’s function, this affects also the pairing amplitude as no pairing is possible in the absence of long-lived quasi-particles. In general, strong superconducting gap and normal self-energy diminish both the Green’s function and the pairing amplitude.

Appendix B Numerical details

The numerical parameters in our calculations include

• •

the number of $\mathbf{k}$-points in the irreducible Brillouin zone, $N_{k}$; we take it to be temperature dependent, growing as temperature is lowered, to be able to capture increasingly sharp Fermi surface, and gain extra precision when the spin boson is nearly critical.

[TABLE]

• •

the cutoff frequency $i\omega_{\mathrm{max}}$ for the Green’s functions, and the frequency above which the data is replaced by the high-frequency tail fit $i\omega_{\mathrm{fit}}$. Throughout the paper we use $i\omega_{\mathrm{fit}}=14.0$ and $i\omega_{\mathrm{max}}=30$. The actual number of Matsubara frequencies taken is therefore temperature dependent.

• •

the number of $\tau$-points is taken simply as the number of frequencies times 3.

• •

the mixing ratio for the polarization between iterations; in $GW$ we take $P^{\mathrm{old}}:P^{\mathrm{new}}=0.95:0.05$. In $GW$+EDMFT and TRILEX, we use $P^{\mathrm{old}}:P^{\mathrm{new}}=0.7:0.3$.

• •

number of iterations performed and the level of convergence reached; in $GW$ we start from the non-interacting solution, and perform up to 70 iterations. In the superconducting phase, we perform 150 iterations. In $GW$+EDMFT and TRILEX, we start from DMFT solution at the highest temperature, and then use the $GW$+EDMFT solution as the initial guess at lower temperature, and perform up to 30 iterations. In all cases, we reach convergence level $\mathrm{max}_{i\omega_{n}}|G^{\mathrm{loc,new}}(i\omega_{n})-G^{\mathrm{loc,old}}(i\omega_{n})|\lesssim 10^{-3}$.

• •

the parameter $\gamma$ used in the LEV extrapolation; in $GW$ for Fig.5 we use $\gamma=0.5$.

Appendix C Extrapolation of the lowest eigenvalue

Because of the AF instability in the methods used in the present paper, there is a need for extrapolating the results for the leading eigenvalue (LEV, $\lambda(T)$) in the linearized gap equation (LGE) to lower temperatures. In Fig. 14 we show some examples of this procedure. The $\lambda(T)$ results are contrasted with $\mathrm{max}_{\mathbf{q},i\nu_{m}}U^{\mathrm{sp}}P^{\mathrm{sp}}_{\mathbf{q}}(i\nu_{m})$ which is shown to approach $1$ at finite temperature. Below this temperature, a stable calculation is not possible. For the precise definition of $T_{\mathrm{AF}}$ shown in figures in IV.2 and IV.1, we follow Ref.Kitatani et al., 2015, and identify it with the condition $\mathrm{max}_{\mathbf{q},i\nu_{m}}U^{\mathrm{sp}}P^{\mathrm{sp}}_{\mathbf{q}}(i\nu_{m})=0.99$ (this value is denoted with a horizontal black line in the bottom two panels of Fig. 14).

The LEV $\lambda(T)$ is found to follow a simple law and we perform a parabola fit

[TABLE]

with $\hat{\theta}=a,b,c$, to extrapolate it to lower temperatures.

Interestingly, a similar $\lambda(T)$ behavior is observed in DCA and DCA+ calculations (see Fig. 15). The fact that the general temperature-dependent behavior of the LEV (as found in the LGE) is captured correctly with respect to DCA, indicates that the leading contribution to $\Gamma^{pp}_{\sigma\bar{\sigma}}$, and therefore the superconducting glue, is indeed bosonic-like, dominated by the RPA-like processes. Otherwise, one would expect a slower decay of $\lambda(T)$ with temperature in DCA than observed in $GW$, as here the decay is determined primarily by the gradual decondensation of the spin boson. This notion has been investigated thoroughly in Ref. Maier et al., where the authors have found both the spin-spin correlation and the $pp$-irreducible vertex from a full DCA calculation to be in excellent agreement with simple random-phase approximation estimates.

In the main text (section III.3), we have estimated the error bar on the extrapolation of the lowest eigenvalue by varying the parameter $\gamma$ (see Fig. 3). Here, we give a method to determine the prediction interval for the extrapolation at fixed $\gamma$. We choose the parameters corresponding to pt.B (Fig.7) to illustrate this method.

Following standard statistics (see e.g. Ref. Ryan, 1997, sec. 13.8.1), we proceed as follows:

(i) for a given doping $n$, we carry out a least-squares fit of the $N$ data points ($T_{i}$,$\lambda_{i}$) to Eq. (113): this yields optimal least-square parameters $\hat{\theta}=a^{*},b^{*},c^{*}$.

(ii) for a given temperature $T_{0}$ (not necessarily in the same range as the data points), the prediction interval at $100\cdot(1-\alpha)$% is given by the two extremal values

[TABLE]

where $\overline{\sigma}$ is the empirical variance

[TABLE]

$t_{\alpha,k}$ is defined as

[TABLE]

where $P_{N}(t)$ is the probability density function of the Student distribution function. $V$ is the $N\times 3$ matrix

[TABLE]

and $v_{0}$ the column vector:

[TABLE]

The corresponding prediction intervals (at 68%) are shown in the upper panel of Fig. 16. They are used to compute the error bars shown in the lower panel of the same figure.

Especially in $GW$+EDMFT, the fit is found to be of high quality and as the extrapolation is not carried far away from the range of data points, the prediction intervals are found to be small. In TRILEX, the fit is of poorer quality and the prediction intervals are comparable to the uncertainty due to free parameter $\gamma$.

Appendix D Superconducting phase at weak coupling

Here, we compare the results of the below-$T_{c}$ calculation: $GW$ at weak coupling (Fig.17) vs. $GW$+EDMFT at strong coupling (Fig.9), at the same dispersion, pt.C. We observe that in the weak coupling case, the normal self-energy remains constant with doping, while at strong coupling it grows by a factor of about 5 in a similar range of doping, as Mott insulating phase at half-filling is approached. In the normal phase and at weak-coupling, the self-energy becomes smaller as half-filling is approached, while the trend is the opposite at strong coupling. On the other hand, the onset of the anomalous self-energy in the anti-nodal regions also seems to reduce the normal self-energy in these regions, therefore making the normal-self energy more local. This seems to be a generic feature, not only associated with the doped-Mott insulator regime. It is particularly interesting that the reduction in $\mathrm{Im}\Sigma$ seems proportional to $S$ in both cases.