Degenerate plaquette physics as key ingredient of high-temperature superconductivity in cuprates

We study the physics of high-temperature cuprate superconductors starting from the highly degenerate four-site plaquette of the t−t′−U\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t-t^{\prime} -U$$\end{document} Hubbard model as a reference system. The degeneracy causes strong fluctuations when a lattice of plaquettes is constructed. We show that there is a large binding energy between holes when a set of four plaquettes is considered. The next-nearest-neighbour hopping t′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t^{\prime}$$\end{document} plays a crucial role in the formation of these strongly bound electronic bipolarons whose coherence at lower temperature could be the explanation for superconductivity. A complementary approach is cluster dual fermion starting from a single degenerate plaquette, which contains the relevant short-ranged fluctuations from the beginning. It gives d-wave superconductivity as the leading instability under a reasonably broad range of parameters. The origin of the pseudogap is also discussed in terms of the coupling of degenerate plaquettes. Thus, some of the essential elements of cuprate superconductivity appear from the local plaquette physics.


INTRODUCTION
After 35 years since the discovery of the high-temperature superconductivity 1 , there is still no consensus on the nature of the mechanism of d-wave pairing in cuprates [2][3][4][5][6][7] . Nevertheless, new experimental findings clearly point to the existence of a quantum critical point around a hole doping of δ ≈ 0.24 [8][9][10] . This concentration separates the exotic bad-metal state for smaller doping from Fermi-liquid behaviour for larger hole concentration with "normal" Fermi-surface described, at least qualitatively, by conventional density-functional theory 11 . Moreover, the carrier density obtained from Hall effect measurements in the largedoping regime is equal to its nominal value n H ≈ 1 − δ while for smaller doping the bad-metal behaviour appears with Fermi arcs, an "enigmatic pseudogap phase" and n H ≈ δ at high temperature 8,12 . Recent investigations of highly overdoped cuprates show that this "strange metal phase" is located around the δ c ≈ 0.24 point 9 . For hole concentrations less than δ c superconducting pairs come entirely from the region of incoherent electrons at the antinode region (X-point) of the Brillouin zone (Planckian dissipators) 9 . Also at δ c ≈ 0.24, electronic specific heat measurements for many different cuprate superconductors have revealed a huge peak in the electron Density of States (DOS) at the Fermi energy in the normal phase, with strong evidence of a Quantum Critical Point (QCP) at this hole doping 8 . Thus, this critical concentration appears to be the first essential element for any theory of these superconductors.
Furthermore, a reliable theory of the high-T c cuprates needs to account for a mechanism for superconducting coupling, i.e., what provides the pair binding energy? The theory should also explain other key experimental observations such as nodal-antinodal dichotomy, pseudogap formation in the underdoped regime and strange metal behaviour.
First-principle electronic structure calculations 11 suggest that a single-band tight-binding model with next nearest neighbour (NNN) hopping and on-site Coulomb interaction, the so-called t À t 0 À U Hubbard model, has the ingredients to describe high-T c phenomena. Moreover, the case of t 0 =t ¼ À0:15 corresponds to the LSCO-cuprate family while one expects t 0 =t ¼ À0:3 to describe cuprate families with higher T c such as e.g., YBCO and Tl2201 13 . The precise role of t 0 is a second important question for the theory.
Finding the exact solution of the t À t 0 À U Hubbard model in the thermodynamic limit at arbitrary interaction strength, doping and temperature is tremendously difficult. Modern computational approaches, ranging from DMRG to Quantum Monte Carlo methods, are making great progress in this direction 6,14 , but a full solution of the problem still seems to be far away. Facing such a complex system that is hard to solve exactly, it is frequently helpful to use simpler, solvable reference systems instead, to get an understanding for the mechanisms responsible for the physics in question. Of course, the crux of the matter is to use the freedom to choose the reference system wisely, since it should contain the essential physics of the full system to be of any use. For example, for the Mott insulating phase, the minimal building block is a single interacting atom coupled to a (dynamical) bath, and extending this to a single bond explains how antiferromagnetic exchange interactions between local moments emerge 15 .
In the case of the doped t À t 0 À U Hubbard model, d-wave superconducting fluctuations are known to be important, and a four-site plaquette is the minimal reference system that contains their spatial structure 16 . This has led to a series of plaquette-based investigations over the years [17][18][19][20] . The plaquette is not just large enough to contain d-wave fluctuations, it also has an important degenerate point 17,21 reminiscent of the experimental picture at 1 δ ≈ 0.24. In analogy with the Kondo model 22 , where the degeneracy of the two spin states of a magnetic impurity plays a crucial role in the anomalous low-energy properties with a correspondingly divergent perturbation series, the degeneracy of the plaquette starting point gives rise to strong fluctuations in plaquette-based methods, which can reveal the nature of the anomalous behaviour of the interacting Hubbard model on a twodimensional lattice. In this manuscript, we discuss how several important aspects of the cuprate phenomenology can be seen to appear when spatial correlations are added to the plaquette. For this purpose, we use two complementary approaches. First, exact diagonalization of a 4 × 4 cluster, i.e., four coupled plaquettes, provides a way to add further short-ranged correlations to the plaquette starting point. It shows a large hole pair-binding energy at suitable values of t 0 and U. Secondly, the dual fermion 23 approach provides a recipe to start from an arbitrary local reference system 24 , in this case, the 2 × 2 degenerate plaquette, and to incorporate nonlocal corrections in a systematic fashion. Dual fermion perturbation theory 25,26 and the dual Bethe-Salpeter equation make it possible to study the momentum structure emerging from longer-ranged fluctuations. Here, it is important to state that the plaquette degenerate point leaves a clear lowtemperature signature in the two-particle correlation function, which is the basic building block of the dual fermion perturbation theory. Thus, large nonlocal corrections are expected to appear as the temperature is lowered.
The first attempt to discuss the plaquette physics as the main ingredient of the high-T c theory was done with the cluster dynamical mean-field theory (DMFT) scheme 16 , and later Altman and Auerbach analytically explained the importance of plaquette two-hole states with d x 2 Ày 2 symmetry 27 . Nevertheless, they did not consider the possibility of a degenerate ground state of the plaquette 17 , with N = 2, 3, 4 electrons per plaquette, at suitable values of t 0 , μ. and U.
We should point out that there is a curve of degenerate plaquettes in the t 0 , μ, U space. Here, we fix t 0 =t ¼ À0:15, the μ and U that correspond to a six-fold degenerate ground state of the plaquette are signified by the star in the Fig. 1. Since we use here periodic boundary conditions the critical Coulomb interaction for plaquette degenerate point becomes U/t = 5.56 in contrast with the case of an isolated plaquette 17 . This is in a very good agreement with the value of the Coulomb interaction U/t = 5.6 that was found in the diagrammatic Monte Carlo calculations 28 in a search of pseudogap formation, and the value of U/t ≈ 6 pointed out in the recent review 8 as the most reasonable value of the effective Hubbard interaction for cuprates. Note also that periodic boundary conditions effectively double t 0 compared to t, which explains the chosen value of the NNN hopping twice smaller than in ref. 17 . At a special value of the chemical potential 17 μ ≈ 0.48 the ground state for the half-filled N = 4 antiferromagnetic singlet is degenerate with the singlet for N = 2 electrons and with two doublets from N = 3 sector. For these values of the parameters the plaquette state corresponds to the hole doping of δ c = 0.25.
In the dual perturbation theory, starting from a degenerate plaquette point leads to divergences in the perturbation series. For the Kondo problem, the dual perturbation starting from the atomic limit 29 has a divergent local four-point vertex at low temperarure, while the Green's function is finite. In the case of the degenerate plaquette both the single-particle and two-particle Green's functions of the reference system are divergent.
We will also consider reference systems differing from the degenerate point in the value of the chemical potential. For smaller μ ≈ 0 (marked with the circle in Fig. 1) the lattice would tend to a metallic behaviour, for larger μ ≈ 0.8 (marked with the square) the perturbation for the lattice results in a superconducting d x 2 Ày 2 instability.

Short-ranged correlations: 4 × 4 cluster
To understand why superconductivity occurs, it is necessary to find a pairing mechanism, i.e., an attractive interaction between pairs of fermions. Although a phase transition can only occur in the thermodynamic limit, finite-size simulations can already point towards the energetic mechanism. We calculated the pairing energy of two holes on the 4 × 4 periodic cluster-which consists of four 2 × 2 plaquettes-through the exact diagonalization (ED) of ground state energies in the different occupation sectors 27,30 , where the energies are measured relative to the half-filled ground state E 0 with no holes,Ẽ Nh ¼ E Nh À E 0 . Note, that Δ 2h < 0 signals pairing. By construction, Δ 2h = 0 for U = 0 and U ≫ t, so it measures genuine correlation effects. Calculated energies for t 0 ¼ 0 are in perfect agreement with the standard ED results 30 . Figure 2 shows the pair binding energy Δ 2h between two holes for a 4 × 4 cluster t À t 0 À U Hubbard model with periodic boundary conditions as a function of interactions strength U for different next-nearest neighbours hopping t 0 . A striking observation is that switching on non-zero negative t 0 =t leads to a distinct minimum in the Δ 2h dependence on U. It can be attributed to the change of the ground state for the sector (7↑, 7↓) (see the  Supplementary Note 3, for a detailed consideration of the case U = 6 where it occurs at t 0 =t % À0:12). The binding of the two holes becomes extremely strong around U = 6 and t 0 =t ¼ À0:3, which is consistent with the estimate for NNN hopping in the cuprates 13 . The pairing energy of the 4 × 4 cluster is of the order of Δ 2h /t ≈ −0.7 which corresponds to ≈ 3000 K for t ≈ 0.4 eV from a generic cuprate model 11,13 , i.e., much larger than the superconducting temperature. This tells us that bound pairs exist for temperatures far above the superconducting region. The superconducting transition should then be seen as the condensation of these pairs. Thus, the cluster binding energy of two holes turns out to be much higher than the superconducting critical temperature which means that the pairs ("bipolarons") should be well-defined also in non-superconducting phase, a situation dramatically different form the conventional BCS superconductivity. The distinction is like the difference between purely itinerant weak ferromagnets and ferromagnets with local magnetic moments which exist until very high temperatures and only order, rather than appear, at the Curie temperature 31 . An important observation is that the pair binding energy on a single 2 × 2 plaquette is an order of magnitude smaller (see Supplementary Note 3), which is an indication that the bound pair of holes is spatially separated. Thus, while individual holes live in a single plaquette, hole-binding requires spatial correlation between neighbouring plaquettes. Plaquette-CDMFT scheme would be unlikely to capture this spatial correlation accurately, instead requiring a 4 × 4 cluster. The negative sign of t 0 =t is considered to be destructive for superconductivity in the t À t 0 À J model [32][33][34] , but it undoubtedly does not hinder hole binding in the t À t 0 À U model considered here. Figure 2 shows that the optimal U increases with t 0 , but also that the hole pair binding on the 4 × 4 cluster disappears for U ≫ W = 8t, i.e., in the strong coupling, t − J limit of the model. The positive sign of t 0 =t does not produce such a strong hole pair binding on the 4 × 4 cluster (see Supplementalary Note 3) opposite to the recent finding in the t − J model 33,34 . Altogether, these results explicitly show the combined importance of t 0 , which greatly increases the pairing energy gain, and U, since in a noninteracting systems Δ 2h = 0 by definition. Similar results of holebinding in a 4 × 4 Hubbard cluster were found recently 35 for a different model of inhomogeneous hopping 36 .
The drastic change of the magnetic correlations with negative t 0 =t is also observed. We analyzed the spin-spin correlation function in the sector (7↑, 7↓) with different NNN hoppings t 0 (Fig. 3) and clearly see a sharp change from antiferromagnertic correlations for t 0 ¼ 0 with a distinct "checkerboard" structure to stripe-like correlations for t 0 =t ¼ À0:3. These stripes point along the x or y directions, while magnetic correlations along other directions are small. A similar reduction of antiferromagnetic and appearance of ferromagnetic correlations with t 0 was found in a lattice QMC study 37 .

Dual perturbation theory-Bethe-Salpeter equation
The exact diagonalization results have shown the tendency towards inter-plaquette pair-binding. Complementary to exact diagonalization, it is possible to calculate the lattice two-particle quantities (vertex, susceptibility) using the dual theory. Recent studies have illustrated the value of the information encoded in vertices and susceptibilities 20,38-47 even in the case of a singleorbital model, which motivates our exploitation of these techniques in a cluster model.
As a starting point, we choose the plaquette parameters where the ground state is six-fold degenerate 17 with U = 5.56, t = −1, t 0 0 ¼ 0:15 μ 0 = 0.48 (the index "0" means the plaquette reference system is allowed to be different from the lattice model) and t 0 =t ¼ À0:15 or −0.3 with μ = 0.7 or 1.5 correspondingly, to keep the optimal doping δ ≈ 0.15 in the square lattice. Additional details of the calculations can be found in the Supplementary Note 2.
The dual Bethe-Salpeter equation shows the strength of interplaquette fluctuations in a lattice of plaquettes. The two ingredients are the single-particle Green's function and the twoparticle correlation function of the plaquette. For a degenerate plaquette, the latter is highly divergent in the limit T → 0, scaling as T −3 instead of T −1 , so strong fluctuations are expected at low and intermediate temperatures. From these plaquette correlation functions, we move to the lattice instabilities of interest via a Bethe-Salpeter equation. It has been shown that divergences of the dual two-particle quantities correspond to divergences of the lattice two-particle quantities 48 , so it is sufficient to consider the kernel of the dual Bethe-Salpeter equation. In the cluster dual fermion theory, a lattice instability (with wave-vector and frequency q, ω) manifests itself by the maximal eigenvalue of the following Bethe-Salpeter matrix Λ αβ reaching unity: λ max = 1 in the case of the particle-particle singlet channel (see the diagram in Fig. 4): With short notations for plaquette sites and fermionic Matsubara frequencies: α = (ij, ν), β ¼ ðkl; ν 0 Þ and taking the wave-vector and bosonic frequency where the transition occurs (q = 0, ω = 0), one needs to diagonalize the matrix Λ P αβ and find the maximum eigenvalue. The definitions ofG and γ are given in the Methods section. Here, we have used that the Bethe-Salpeter equation has an intertwined spin, site and frequency structure which can be simplified by looking at the different channels. Since our main interest is superconductivity, the formula given above corresponds to the singlet particle-particle channel. For comparison, we have also considered the particle-hole density and magnetic channels. Regarding the frequencies, we restrict ourselves to the lowest 10 Matsubara frequencies, since the vertex function γ P decays strongly with ðν; ν 0 Þ (see the Supplementary Note 2).
In this case, the matrix γ P is Hermitian (real for ω = 0), while the matrix Λ is not Hermitian, but the leading eigenvalues are still found to be real for all channels (we also calculate eigenvalues of corresponding Bethe-Salpeter equations in the density and magnetic particle-hole channels).
Results for the maximum eigenvalues of the Bethe-Salpeter matrix Λ at the critical point are presented in Fig. 4. The density and particle-particle eigenvalues cross unity at β ≈ 5 and 15 < β < 20.
The eigenvector corresponding to λ max for the particle-particle singlet case has d x 2 Ày 2 symmetry in the plaquette space. Exactly at the plaquette degenerate point the instability (signaled by λ crossing 1) in the density channel is very large because the N = 2, 3, 4 states are degenerate. We found that this density instability is not robust against changes of μ 0 and as soon as we shift it towards low hole doping μ 0 = 0.8 there is no density instability. On the other hand, the singlet superconducting instability is very robust and becomes the leading one for doping lower than δ = 0.25. The magnetic instability does not play any role for the doped case and becomes the leading one only in the half-filled case.

Spectral information
Due to the degeneracy of states with different particle numbers, the density of states of the plaquette is large close to the Fermi level. The availability of low-energy states is the driving force behind the instabilities that occur once a lattice of plaquettes is considered. At the same time, coupling of the electrons to these fluctuations moves spectral weight away from the Fermi level.
In Fig. 5 we compare the density of states (DOS) for plaquette DF second-order perturbation (DF2) with the so-called cluster perturbation theory (CPT) which corresponds to zero dual-self energy in Eq. (7) for a relatively high temperature (β = 3). To obtain the spectra, we used Padé analytical continuation from Matsubara to the real energy axis 15 . One can see that the DOS for the dual fermion theory is more sharply peaked near the Fermi level compared to the CPT-result. For comparison, we also show the ED result for the plaquette with a sharp peak exactly at Fermi level due to six-fold degenerate ground state. In this case, there is still no signature for a pseudogap and the lattice self-energy is "well-behaved".
In Fig. 6 we compare the DOS for the plaquette DF perturbation theory for a lower temperature (β = 5) with the ED results for the 4 × 4 cluster in the sector (7↑, 7↓), which corresponds to a 2 × 2 lattice of plaquettes. These two methods are complementary: the DF approach is perturbative in the interplaquette coupling and can handle large lattices, whereas the ED is exact but limited by the cluster size. From the comparison of the two curves, we conclude that the dual fermion theory shows a tendency towards pseudogap formation which is clearly seen in the ED results. It is possible to assume that the pseudogap in the 4 × 4 cluster is related to the coherent interactions of the large peak on the DOS in the reference plaquette or Fano-like effect of interactions with the "soft fermion mode" of the lowlying excitations which are encoded in the local vertex functions of the DF-approach. In this sense the pseudogap physics for the optimal doping may be related more to the "hidden fermion" physics 49,50 or "destructive interference phenomena" 51 than to (short-ranged) antiferromagnetic fluctuations, since one would expect those to already be present in the plaquette result. Furthermore, the analysis of the dual Bethe-Salpeter equation showed that the magnetic channel is not dominant at the parameters considered here. Leading instabilities. The maximum eigenvalues of the Bethe-Salpeter equation (BSE) for the particle-particle singlet (PPs), density (Den) and magnetic (Mag) channel for doped plaquette with U = 5.56 and t 0 = −0.15t, μ = 1.55. Insert: diagrammatic representation for BSE-matrix in the particle-particle channel.

DISCUSSION
There appears to be a close relation between the physics of cuprate superconductors, with the clear existence of a quantum critical point at δ c ≈ 0.24, and the degeneracy of the plaquette in the strong-coupling regime. In this sense, the plaquette and not the single site can be considered the minimal building block for cuprate physics, with pair binding arising in a lattice of plaquettes.
Exact diagonalization shows that the cluster pair binding energy is dramatically enhanced when four plaquettes are considered together, compared to the pair binding energy in a single plaquette. Given their large binding energy, these pairs should probably exist also at much higher temperatures than the superconducting critical temperature, remaining noncoherent. The exact diagonalization also shows the important role played by the next-nearest hopping t 0 , with a large pair-binding energy at t 0 =t % À0:3.
Dual fermion expansion starting from the plaquette reference system provides a complementary way to investigate interplaquette correlations. For the doping δ ≤ 0.25 the dual Bethe-Salpeter equation clearly shows the presence of a low temperature d x 2 Ày 2 instability, which has an eigenvalue substantially larger than the magnetic channel. Starting from the degenerate plaquette, fluctuations in the density channel are also very strong, but these seem to be less robust against changes in the filling.
The exact diagonalization of the 4 × 4 cluster as well as renormalized dual fermion perturbation starting from the plaquette reference system with δ = 0.25 also uncovers spectral consequences of this degeneracy. The formation of the pseudogap can be seen as the destructive interference or a Fano-like effect originating from the sharply peaked DOS in the isolated plaquette embedded into the band of surrounding fermions, as was hypothesised in ref. 17 . These observations about the mechanisms of superconductivity can all be made by starting the perturbation theory from an isolated plaquette. For more quantitative predictions of the theoretical phase diagram, the optimal dynamical embedding of the plaquette and the implications for the resulting perturbation theory need to be studied further.

Cluster Dual Fermion approach
We used the standard exact diagonalization method 30 for small systems as well as the special version of the cluster dual fermion scheme 23,26 for t À t 0 À U square lattice Hubbard model. The general strategy of the dual fermion approach is related to formally exact separation of the local-plaquette and non-local hybridization (Fig. 7). The details of the path-integral formulation of this approach can be found in the Supplementary Note 1.
We start from the following general lattice action and rewrite it as a sum of non-connected plaquette reference systems and the remaining coupling term: where ν = (2k + 1)π/β, with k 2 Z, are the fermionic Matsubara frequencies, β is the inverse temperature, τ is the imaginary time in the interval 0; β ½ Þ, μ is the chemical potential,t k is the hopping matrix downfolded onto the site-orbital space of the plaquette (see Eq. (11) below), and the Grassmann fields c, c * are vectors in the same space. The index i labels the lattice sites, σ is the spin projection and the k-vectors are supercell plaquette quasimomenta. In order to keep the notation simple, it is useful to introduce the combined index 1 j i i; n; σ; τ j i(n being the plaquette site index suppressed above) while assuming summation over repeated indices. The summation over Matsubara frequencies ν includes a normalization factor 1/β and the k integration is normalized by the volume of the reduced Brillouin zone. The general reference system is defined by a plaquette matrixΔ ν , which is also allowed to be instantaneous 24 (ν-independent). It can contain hopping inside the cluster as well as possible frequency-dependent connections to an auxiliary fermionic bath. The reference plaquette has the same local plaquette interaction matrixÛ, as illustrated in Fig. 7, and the corresponding action is: In this work, we restrict ourselves to instantaneousΔ. The main motivation for using the simple staticΔ is that such a reference system can be solved numerically using Exact Diagonalization (ED), without the introduction of "bath sites" and fitting parameters, and without the numerical costs and noise of continuous-time Quantum Monte Carlo (CT-QMC) 23,52 , which is able to treat general, frequency-dependent hybridizationΔ ν . In this work, we use an isolated plaquette cluster with periodic boundary conditions as the reference model, see Eq. (12).
Having solved the reference system exactly, including the calculation of all relevant correlation functions, we can derive an efficient perturbation series in the "coupling term"t kν t k ÀΔ ν À Á which is equivalent to solving of the effective dual fermion (d * , d) action and describes non-local correlation effects beyond the reference plaquette 23,23 : where the bare dual Green function has the form withĝ ν being the local Green's function matrix for the plaquette. The vertex γ P is given by the connected part of the local two-particle  Supplementary Notes 1 and 4). The action Eq. (5) allows us to calculate the dual self-energy,Σ kν with a level of approximation of our choice, e.g., using diagrammatic perturbation theory. Once this is done, the results are transformed back using an exact relation between the dual and the lattice Green's functions (see Supplementary Note 1): Perturbation in dual space The previous discussion left open the question of how to determineΣ. In this work, we use cluster dual fermion perturbation theory (Fig. 7) for the action (5). We start with the definition of interaction between dual fermions, using the particle-hole notation for the local vertex and writing explicit spin indices and Matsubara frequencies of the connected two particle Green's function for the original fermions (c * , c) as follows 23,53 : In the Matsubara space, the vertex depends on two fermionic frequencies, ν; ν 0 , as well as one bosonic frequency ω. For the sake of completeness and the reader's convenience, we mention the connection between the particle-particle and the particle-hole notation, γ 1234 ðν; ν 0 ; ωÞ ¼ γ P 1342 ðν; ν 0 ; ν þ ν 0 þ ωÞ. The bare vertex of the dual fermion perturbation theory is the full connected correlation function of the reference system. The present vertex differs from the usual dual fermion expression due to the different rescaling factor of the Hubbard-Stratonovich field. Here, we avoid amputation of the legs of the vertex, which requires division by Green's functions at all external points.
It is useful to symmetrize the vertex into charge density (d) and magnetic (m) channels: γ d=m 1234 ðν; ν 0 ; ωÞ ¼ γ "" 1234 ðν; ν 0 ; ωÞ ± γ "# 1234 ðν; ν 0 ; ωÞ Now we can write the first-order dual fermion self-energy which is local in plaquette space (Fig. 8): The second order Feynman diagram for the dual fermion perturbation (Fig. 8) in real space (R ij ) has density-and magnetic-channel contributions, with corresponding constants (c d ¼ À 1 4 and c m ¼ À 3 4 ): In principle, one can go beyond the second order perturbation expansion and include dual ladder diagrams 53,54 , dual parquet diagrams 55 or a stochastic sum of all dual diagrams with the two-particle vertex γ 1234 , using diagrammatic Monte Carlo in dual space [56][57][58] . In addition, the diagrammatic series can be made self-consistent, using dual skeleton diagrams and "bold" lines. Finally, one can also update the reference system (and obtain a frequency-dependent Δ). These approaches are all numerically more involved. As the main goal of the present work is not to present quantitatively reliable results but rather to highlight the connection between the degenerate reference system and the superconducting fluctuations, we will mostly stick to the second-order consideration.
The calculations shown here were performed using a Fortran implementation of dual fermions that uses the equivalence of the four sites in the plaquette to speed up the vertex calculation. The results were checked against an open-source implementation of the second-order dual fermion perturbation 59,60 , based on TRIQS 52 and with pomerol 61 as an impurity solver as well as cross-checked with the momentum-space cluster dual-fermion scheme 62 .
Since length scales are an important part of this manuscript, it is useful to discuss how these enter the dual (perturbation) theory. In particular, what are the spatial implications to cutting off the perturbation series at second order? First of all, correlations contained entirely within the size of the cluster reference model are handled by the impurity solver, so we should mainly worry about longer-ranged correlations. Away from U ≈ 0, electronic correlations lead to localization and the single-particle dual Green's function has a finite, short-range. For the second-order diagram, ΣðRÞ $G 3 ðRÞ, so in the second-order approximation, a short-ranged dual Green's function leads to a short-ranged dual self-energy. However, higherorder diagrams are able to cover the distance R by repeated short-ranged propagation, and therefore dominate at longer lengths. Length scale effects in dual fermion are discussed in more detail in refs. 63,64 . Perhaps the best example is the Heisenberg limit at half-filling and U ≫ t, where the spin susceptibility is very long-ranged even though the individual electrons are almost perfectly localized.

Plaquette tight-binding scheme
We study the optimally doped square lattice Hubbard model, with nearest neighbour hopping t and NNN hopping t 0 . As illustrated in Fig. 7, the original lattice can be reconsidered as a lattice of 2 × 2 plaquettes. Every unit cell of the plaquette lattice contains 4 atoms of the original lattice, as shown on the left-hand side of Fig. 7. The plaquette lattice has the following 4 × 4 hopping matrix (see Fig. 7), where the functions K mn k and L mn k , with m, n ∈ {−1, 0, +1}, are defined as K mn k ¼ 1 þ e iðmkx þnky Þ L mn k ¼ 1 þ e iðmkx þnky Þ þ e imkx þ e inky We will use a single plaquette as the reference system. Compared to the single-site dual fermion formalism, this plaquette reference system already encompasses the short-ranged correlations that are essential in this system.
In the dual fermion approach, there is a general freedom of choosing the most appropriate reference system. One way to construct a plaquette reference system would be to simply remove all black links in Fig. 7 (and attach the remaining sites to a bath). This is equivalent to the selfconsistent cluster-DMFT scheme 16 and corresponds to averaging over the supercell Brillouin zone. This scheme, however, eliminates exactly half of the nearest-neighbour hoppings and three quarters of the next-nearestneighbour hoppings.
Here we choose another path and consider plaquettes with periodic boundary conditions as a static reference system. In terms of the a b Fig. 8 Dual self-energy. Feynman diagram for the first order (a) and the second order (b) dual fermion perturbation for the self-energy e Σ 12 ðνÞ: a line represents the non-local dual Green's function e G 43 ðν 0 Þ and a box is the local plaquette vertex γ 1234 , ðσ; σ 0 Þ are spin-indices.
supercell Brillouin zone, this corresponds to achieving self-consistency for k = 0 only, instead of the momentum average. The intra-plaquette hopping reads Note that we include the possibility of using a different chemical potential μ 0 = −ε 0 in the reference system, compared to that of the lattice model μ =−ε to adjust the hole dopping to about δ = 0.15. We fix the nearest neighbour hopping t but retain the freedom of adjusting the next nearest neighbour hopping t 0 in the dual fermion transformation. For example this may be used to reduce the factor 4 for the t 0 hoppings for the periodic boundary conditions for 2 × 2 plaquette if we chose t 0 0 ¼ t 0 =2. With the plaquette as the reference system, one can use the exact diagonalization approach to calculate the dual Green's function and the plaquette vertex function 65 .

DATA AVAILABILITY
The data supporting the present work are available from the corresponding author upon request.