Strongly correlated superconductivity with long-range spatial fluctuations

Since unconventional superconductivity is regularly found in the vicinity of a magnetically ordered state [3, 4], it is widely believed that spin fluctuation effects should be an (if not the) important key for understanding unconventional superconductivity. Often these systems are modeled by a single-band Hubbard model [55, 56, 57, 58, 59]. Also for this model, fluctuation diagnostic analyses could show the importance of spin fluctuations in the pseudogap regime [60, 61, 62]. Diagrammatically a particle-hole ladder diagram is relevant for spin fluctuations. We show the simplest example in Fig. 1(a), where the ladder is mediated by the bare interaction (the on-site repulsion $U$ for the simple Hubbard model). Such diagrams are used, e.g., in the fluctuation exchange (FLEX) approximation [10]. When we consider the relevant diagrams from a DMFT starting point (for instance in the dynamical vertex approximation [42, 45]) similar particle-hole diagrams enter, but now mediated by the irreducible vertex $\Gamma$ instead of $U$. Here, an irreducible vertex is a diagram that can not be split into two pieces by cutting two internal Green function lines in the corresponding channel: particle-hole (ph) or particle-particle (pp). This quantity can be derived from the one- and two-particle Green functions $G,\chi_{\rm ph}$ through the inverse Bethe-Salpeter equation:

where $\chi_{0}$ is the bare (bubble) susceptibility. Please see Refs. [63, 64] for detailed explanations for the vertex function itself. Then, the diagram for the fully reducible vertex $F$ becomes like Fig, 1(b). From it we can calculate the self-energy by using the equation of motion (EOM, also known as Schwinger-Dyson equation in this context) which is an exact relation connecting the vertex and the self-energy. The basic concept is common to all diagrammatic extensions while the building block differs: the bubble (and ladder) diagrams in GW+DMFT (FLEX+DMFT) [37, 38, 48, 54, 49], the fully reducible local vertex in dual fermion [43] or DMF²RG [47], $\Gamma$ in ladder D$\Gamma$A [42, 45] and the fully irreducible vertex in parquet D$\Gamma$A [50], the three point fermion-boson vertex in TRILEX [51] etc.. In ladder schemes which we regularly use nowadays, we need to choose one relevant channel. Since the Coulomb interaction is repulsive and we often study the repulsive Hubbard model, we usually pick up particle-hole ladder diagrams. If we study the attractive model, the relevant channel will change and we should first pick up particle-particle diagrams [65]. Furthermore, solving the parquet extension (treating the contributions of all channels on equal footing) from DMFT result have been developed [50, 66] and recently applied to study superconductivity instability, however, for high temperatures [67].

The general way for approaching the symmetry broken superconductivity state is the extension to Nambu space, which is quite hard to treat. Instead, we often linearize anomalous components focusing on the infinitesimal gap function using the normal state quantities only, which is enough for calculating the transition temperature. In this case, we need to solve the eigenvalue problem (linearized gap equation) for the superconducting (particle-particle) channel [69, 49],

where $G$ ($\Delta$) is the normal (anomalous) Green function, and $\Gamma_{\rm pp}$ is the irreducible vertex in the particle-particle (cooperon) channel, depending on the four-vectors $k,k^{\prime},q$ (momentum and Matsubara frequency), $\beta$ is the inverse temperature, $N_{k}$ is the number of k-points, and $\lambda$ is the eigenvalue of the kernel of the superconducting channel which is the key quantity. When $\lambda$ reaches unity the gap function has a finite solution, indicating the superconducting phase transition. In other words, $\lambda\rightarrow 1$ means the divergence of the Bethe-Salpeter equation in the particle-particle channel mediated by the pairing vertex $\Gamma_{\rm pp}$. A typical diagrammatic structure of the spin-fluctuation mediated $\Gamma_{\rm pp}$ is shown in Fig. 1(c). The vertex $\Gamma_{\rm pp}$ describes the scattering of Cooper pairs which corresponds to the pairing glue for superconductivity. As indicated in Fig. 1(c) in ladder D$\Gamma$A calculations we use the particle-hole ladder diagrams [blue dashed box or Fig. 1(b)] that correspond to spin (and charge) fluctuations to calculate $\Gamma_{\rm pp}$. In a more complete but also numerically much more expensive parquet calculation [67], also the feedback of the particle-particle fluctuations ofn the particle-hole fluctuations is taken into account.

One crucial issue is how to perform low-temperature calculations. This is always a problem for studying unconventional superconductivity because of its typically low transition temperature: $T_{\rm c}\!\lessapprox\!({\rm bandwidth})/500$. Since vertices depend on three frequencies, we can store and calculate only a limited number of Matsubara frequencies which is usually insufficient for studying unconventional superconductivity. We explain how to technically overcome this problem and summarize recent progress on this topic in the Appendix.

As mentioned in the previous section, the fully reducible vertex function is calculated through the Bethe-Salpeter equation. The divergence of the Bethe-Salpeter series indicates a physical phase transition since it should connect to the divergence of the susceptibility in the corresponding fluctuation channel. However, it was found that the irreducible vertex itself may diverge without phase transitions [70, 71, 72, 73]. In Fig. 2(a) we show red/orange lines where $\Gamma_{\rm d/pp}$ diverges without any physical instabilities in the DMFT calculation for the two-dimensional square lattice Hubbard model (taken from Ref. [71]). Please note that this divergence is not an artifact of the DMFT approximation and are already present in the (exactly solvable) Hubbard atom [73]. Also, some recent studies show that these lines change if we include spatial fluctuations [74]. It has been demonstrated recently that these divergences are responsible for a reduction in the local charge response [75] and, at the same time, a strong enhancement of the charge compressibility in the vicinity of the Mott transition [76].

Here, we additionally analyze the eigenvalue $\lambda$ of the Eliashberg equation when such divergence occurs. The leading eigenvalue in the superconductivity channel $\lambda$ is often referred to as a “measure of $T_{\rm c}$”, since $\lambda$ usually monotonically increases as we decrease the temperature and $\lambda\!=\!1$ corresponds to the divergence of the superconducting susceptibility. However, this is not the case if divergences in the irreducible vertex occur, as shown in the following.

Let us consider eigenvalues in the local system as a simple example where the DMFT becomes exact for a corresponding Anderson impurity model. We then calculate the local Eliashberg equation

where the local Green function ($G^{\rm loc}$) depends on the Matsubara frequency $\omega_{n}$, and the local particle-particle irreducible vertex for the singlet channel ($\Gamma^{\rm loc,s}_{\rm pp}$) can be calculated within DMFT.

In Fig. 2(b), we show the (three maximal and four minimal) eigenvalues $\lambda$ of the local Eliashberg equation with the same model parameters as Fig. 2(a) at $T/D=0.1$, where $D=4t=1$ is the half-bandwidth. At first glance, the $U$-dependence of the eigenvalues is not systematic (dark red/blue points). For the region without any divergence in $\Gamma_{\rm pp}$ ($0<U/D<1.6$), all eigenvalues are smaller than zero, reflecting the fact that the repulsive interaction does not induce a momentum-independent $s$-wave superconductivity. However, as we cross the orange line at $U_{1}/D\!\approx\!1.6$ in Fig. 2(a) indicating the divergence in $\Gamma_{\rm pp}$, one large positive eigenvalue appears in the superconducting channel. When crossing the second divergence line at $U_{2}/D\!\approx\!1.9$, we can see that another large eigenvalue appears and a large negative eigenvalue disappears. Concomitantly, the eigenvalue $\lambda$ that was largest in $U_{1}<U<U_{2}$ smoothly connects to the second largest within $U_{2}<U<U_{3}$ ($U_{3}$: third divergence point), with the new large eigenvalue becoming the largest $\lambda$. Similarly, for the negative eigenvalues, the second lowest $\lambda$ within $U_{1}<U<U_{2}$ continues as the lowest $\lambda$ within $U_{2}<U<U_{3}$. The same happens again when crossing the third divergence line. This indicates that the large positive eigenvalue can appear through the infinitely large eigenvalue instead of crossing unity. The number of eigenvalues greater than unity corresponds to the number of vertex divergence lines crossed, but there is no phase transition as $\lambda$ remains $\neq 1$. This result can be naturally understood in terms of the divergence in $\Gamma_{\rm pp}$. Actually, the original article [70] analyzed the eigenvalues of the generalized susceptibility $\chi$ and found that its eigenvalues can be negative. This corresponds to a divergence in the irreducible vertex and, hence, is represented by the orange lines in Fig. 2(b) if we think about $\chi\sim\frac{\chi_{0}}{1+\Gamma_{\rm pp}\chi_{0}}$.

In the following we show how that the superconducting phase transition in the strongly correlated regime is actually governed by the eigenvalues $\lambda$ which are close to unity, rather than by the leading (i.e. largest) eigenvalue. Since the pairing interaction matrix $\Gamma^{\rm singlet,triplet}_{\rm pp}(k,k^{\prime})$ is Hermitian, we can perform an eigenvalue decomposition of the vertex. Taking care of the fact that the right and left eigenvectors of the Eliashberg equation (Eqs.(2,3)) are different, we can reconstruct the (irreducible and full) vertex by using the eigenvalues and (right) eigenvectors of the Eliashberg equation as

Here, the right eigenvectors of the Eliashberg equation $\{\Delta_{i}\}$ are normalized as $\sum_{k}\Delta^{*}_{i}(k)G^{*}(k)G(k)\Delta(k)_{j}=\delta_{i,j}$ where $\delta_{i,j}$ is the Kronecker delta (which means the right and left eigenvectors make a bi-orthogonal basis set). Such a low rank decomposition of the vertex function is useful for approximating the complicated vertex structure, and indeed had been applied in the FLEX+T-matrix approach [77, 78] (while they treated as if $\Delta_{i}$ are orthogonal).

Eqs. (4,5) clearly indicate that the divergence of $\lambda_{i}$ only connects to the divergence of $\Gamma$ [Eq. (4)] and not to the divergence of fully reducible vertex $F$ [Eq. (5)] which is immediately relevant for physical phase transitions (since it contains all the vertex corrections of the physical susceptibility). Therefore, for the purpose of analyzing phase transitions, we should not study the leading eigenvalue but eigenvalues close to unity since $\lambda\rightarrow 1$ still means a diverging $F$ and, thus, a phase transition. One important open question is how the eigenvalues above unity behave, which would, according to Eq. (5), be related to a negative contribution to the local susceptibility [75].

For purely two-dimensional systems the Mermin-Wagner theorem [79, 80, 81, 82] forbids continuous symmetry breaking at finite temperatures. Nevertheless, we sometimes refer to finite superconducting transition temperatures $T_{\rm c}$ [spontaneous continuous U(1)-gauge symmetry breaking]. To what extent is this justified?

First, let us explain how the Mermin-Wagner theorem can be understood (not proven) within diagrammatic calculations from the paramagnetic region. The self-energy $\Sigma(k)$ is connected to the fully reducible vertex function $F(k,k^{\prime},q)$ through the equation of motion, and, for the extremely fluctuating region, we can assume a specific channel dominates in the vertex contributions $F$ (Fig. 3) like,

where $\chi$ is the susceptibility of the relevant channel and $\gamma$ is the three-point vertex in that channel (c.f., [83, 84, 85]). For the transformation from Eq. (7) to Eq. (8), we assume an extremely fluctuating situation around $q\approx Q$ where the relevant susceptibility can be approximated as $\chi(q)\propto\frac{1}{(q-Q)^{2}+\xi^{-2}}$, where $\xi$ is the correlation length. Please note that for $\xi\rightarrow\infty$ we can replace $q$ by $Q$ in all terms but $\chi(q)$, perform the $k^{\prime}$ summation which yields a non-diverging term $f(k)$. Now the crucial difference regarding dimension is that the self-energy diverges in a two-dimensional system while it converges to a finite value in three dimensions (for detail differences of the spin-fluctuation mediated self-energy between in 2D and 3D, see, e.g., [84, 85]). Please note that this integral ($\int_{q\approx 0}dq/q^{2}$) is the key quantity for the rigorous proof [79] using Bogoliubov’s inequality. Also, such a difference behavior for the long-range fluctuation is essential for determining the critical behavior [86, 87, 88, 35].

The difference in Eq. (12) indicates that, for a two-dimensional system, the self-energy must diverge as the susceptibility approaches infinity while it has not to be the case for a three-dimensional system. As the relevant susceptibility increases, the self-energy damping effect will win at some point, and then the susceptibility is not enhanced anymore (For the specific method, e.g., fluctuation exchange approximation, there is a more detailed analytical explanation of why FLEX satisfies the Mermin-Wagner theorem for antiferromagnetic long-range orders [83]). Also, for the $d$-wave superconductivity, it has been discussed that the similar singular structure in 2D suppresses $T_{\rm c}$ for the superconductivity long-range order [78]. In this regard, we can understand that we now obtain a finite superconducting $T_{\rm c}$ if we do not include the $d$-wave superconducting fluctuation effect to the self-energy.

For performing more appropriate calculations, we need to include all fluctuation effects and study the (quasi-two-dimensional) three-dimensional model, which is too complicated for analyzing electronic lattice models. This question has been addressed more carefully for spin systems (c.f., a study even for the electron lattice model if limited to spin fluctuation effects [88]). In Ref. [89], the authors performed the (quantum) Monte-Carlo simulations for (quasi-one/two-dimensional) three-dimensional systems. For both cases the Néel temperature should go to zero as $J_{\perp}/J\rightarrow 0$, satisfying the Mermin-Wagner theorem ($J$ and $J_{\perp}$ are the in-plane and inter-plane Heisenberg coupling constant, respectively). On the other hand, they show quite different behaviors: $T_{\rm N}\propto-1/\ln{(\alpha J_{\perp}/J)}$ ($\alpha$: constant) in 2D while $T_{\rm N}\propto J_{\perp}/J$ in 1D probably reflecting the difference divergence strength in Eqs. (12).

Indeed, the authors of Ref. [89] found a universal behavior of the transition temperature $T^{\rm AF}_{\rm c}$ in the weak inter-layer coupling region, indicating the relation between $T^{\rm AF}_{\rm c}$ in quasi-1D/2D system and the singular behavior in the purely 1D/2D system. In this region, the $T^{\rm AF}_{\rm c}$ is quite stable in 2D, which changes only $10-20\%$ if we change the coupling strength orders of magnitude $0.001<J_{\perp}/J<0.1$, which we can assume as the typical quasi-two-dimensional materials. In this sense, even without the inter-layer coupling, we can roughly define the transition temperature as the stable value of weakly coupled systems. We should also mention that such a very weak inter-layer coupling is not relevant at all if we do not include extremely fluctuating critical $d$-wave pairing fluctuations, exemplified by the fluctuation exchange calculations [90].

The question remains of how much is the effect of introducing the $d$-wave superconducting fluctuation effect to the self-energy in quasi-two-dimensional systems (i.e., without diverging effect on the self-energy described above). While this point has not been fully understood yet, a dynamical cluster approximation (DCA, a cluster extension of DMFT) study [60] showed that the $d$-wave superconducting fluctuation has only a tiny effect on the self-energy even close to the superconducting transition temperature (please note that due to the finite size effect of the clusters, this argument can not be applied to the ideal two-dimensional physics described above, while the calculation was done in pure two dimension).

Summarizing the above-mentioned points, it is naively expected that without including the inter-layer coupling explicitly, we can roughly define the transition temperature in quasi-2D systems, which can be estimated by a purely two-dimensional system without $d$-wave pairing fluctuation effect on the self-energy (since the inter-layer coupling is too small to change the physics without critical $d$-wave pairing fluctuation [90], and for $d$-wave pairing fluctuation, the feedback effect for the self-energy is found to be minor [60] except for singular behavior of long-range fluctuation limit in a pure two-dimensional system). This point would be in sharp contrast to quasi-one-dimensional systems where it will be quite difficult to perform the quantitative argument without explicitly defining the inter-layer coupling.

Finally, $d$-wave pairing fluctuations would dominate if we approach $T_{\rm c}$ due to its divergent contribution in the purely two-dimensional model. In this sense, we can regard the current result as the extrapolation from a bit higher temperature calculations, where such critical fluctuations (introducing a singularity) are not relevant. The details of such a dimensional crossover on the superconductivity remain open. Further studies similar to Ref. [89] for the electronic lattice system are desired, i.e., including the $d$-wave superconducting fluctuation and comparing between purely two-dimensional systems and weakly coupled quasi-two-dimensional layered systems. We also note a possibility of the Berezinskii-Kosterlitz-Thouless (BKT) transition [91] regarding the superconductivity in the two-dimensional Hubbard model [92, 93]. In this case, $T_{\rm c}$ remains as a quasi-long-range order even in purely two dimensions, so that the quasi-two-dimensional $T_{\rm c}$ is rather stable against the inter-layer coupling.