Strongly correlated superconductivity with long-range spatial fluctuations

We review recent studies for superconductivity using diagrammatic extensions of dynamical mean field theory. These approaches take into account simultaneously both, the local correlation effect and spatial long-range fluctuations, which are essential to describe unconventional superconductivity in a quasi-two-dimensional plane. The results reproduce and predict the experimental phase diagrams of strongly correlated system such as cuprates and nickelates. Further studies reveal that the dynamical screening effect of the pairing interaction vertex has dramatic consequences for the transition temperature and may even support exotic mechanisms like odd-frequency pairing. We also discuss the dimensionality of layered materials and how to interpret the numerical results in two dimensions.


Introduction
Layered and thin-film superconductivity has been a central topic in condensed matter physics due to its fascinating properties. Especially after the discovery of the cuprate [1] and iron-based [2] superconductors, such quasi-two-dimensional systems can be considered central for unconventional superconductivity [3,4]. In 2019 superconductivity of thin nickelate films was discovered [5]. For understanding these systems, it is essential to describe appropriately correlated electrons of 3dtransition metals, including spatial fluctuation effects that are amplified by the quasi-low dimensionality.
For the theoretical understanding of unconventional superconductivity, Feynman diagrammatic perturbation techniques have often been used [6,7]. Such schemes can describe spatial fluctuation effects by collecting relevant diagrams, and they have succeeded in describing spin-fluctuation mediated superconductivity [8,9,10,11]. On the other hand, these methods have difficulties to capture non-perturbative phenomena such as the correlation-induced metal-insulator (Mott) transition [12,13], one of the cornerstones of correlated electron systems. In this regard, the dynamical mean field theory (DMFT) is the method of choice for treating Mott transitions [14,15]. DMFT becomes exact in the limit of infinite spatial dimensions [16]. In contrast, spatial fluctuation effects ignored by DMFT become more relevant for low dimensional systems, e.g., layered materials or thin films that are further away from this limit.
Against this background there have been many attempts to capture both the strong correlation effects and (long-ranged) spatial fluctuation effects simultaneously. In particular, two main types of extensions of DMFT for including such spatial fluctuations have been established: one class are the cluster extensions which use the numerically exact self-energy a self-consistently determined cluster model [17,18,19,20]. While it can describe the overall structure of the phase diagram of unconventional superconductors like cuprates [21,22], it has been gradually understood that the finite size effect of the cluster is severe, e.g. for capturing the gap structure induced by strong spin fluctuation in the intermediate coupling regime [23,24,25]. The other class of extensions is the diagrammatic route [26]. As ordinary diagrammatic methods, we need to choose how to collect diagrams beyond the DMFT starting point and a systematic improvement is cumbersome [27]. Considering additional physical requirements such as causality [28] is thus helpful. On the other hand, one can address more exotic situations including lower temperatures, longer ranged correlations and multi-orbital [29] physics. Thus these approaches can, e.g., capture accurately the (pseudog-)gap opening due to the strong long-ranged magnetic fluctuations [23,25,30,31]. One can analyze even extremely fluctuating systems like (quantum) critical regimes [32,33,34,35] and estimate the degree of anisotropy of electronic correlations in layered materials [36].
In this article, we review recent studies of unconventional superconductivity with diagrammatic extensions of DMFT. This article is organized as follows. In Sec. 2 we overview methods for normal state calculations and the analysis of the superconductivity. In Sec. 3 we argue some relevant (but rarely discussed) points for studying the layered and thin-film superconductivity within these frameworks. In Sec. 4 we show recent results on the superconductivity. We finally note the summary and outlooks in Sec. 5.

Normal state calculations
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 Γ 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 twoparticle Green functions G, χ ph through the inverse Bethe-Salpeter equation: where χ 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 2 RG [47], Γ in ladder DΓA [42,45] and the fully irreducible vertex in parquet DΓ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].

Superconductivity analysis
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 Γ pp More general particle-hole diagrams, including vertex corrections in strongly correlated systems (Γ: the irreducible vertex for the particle-hole channel). (c) Such spin fluctuations act as a pairing glue for superconductivity in the particle-particle channel.
Here, an exemplary diagram for the particle-particle irreducible vertex Γ pp (blue dotted box) is shown where Γ pp is made up from the particle-hole fluctuations from (b). Taken from Ref. [68].
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 (∆) is the normal (anomalous) Green function, and Γ pp is the irreducible vertex in the particle-particle (cooperon) channel, depending on the four-vectors k, k , q (momentum and Matsubara frequency), β is the inverse temperature, N k is the number of k-points, and λ is the eigenvalue of the kernel of the superconducting channel which is the key quantity. When λ reaches unity the gap function has a finite solution, indicating the superconducting phase transition. In other words, λ → 1 means the divergence of the Bethe-Salpeter equation in the particle-particle channel mediated by the pairing vertex Γ pp . A typical diagrammatic structure of the spin-fluctuation mediated Γ pp is shown in Fig. 1(c). The vertex Γ pp describes the scattering of Cooper pairs which corresponds to the pairing glue for superconductivity. As indicated in Fig. 1(c) in ladder DΓ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 Γ 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 c (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.

Remarks on calculating superconductivity
3.1. Eigenvalues with divergence of Γ 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 Γ 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 λ of the Eliashberg equation when such divergence occurs. The leading eigenvalue in the superconductivity channel λ is often referred to as a "measure of T c ", since λ usually monotonically increases as we decrease the temperature and λ = 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 loc ) depends on the Matsubara frequency ω n , and the local particle-particle irreducible vertex for the singlet channel (Γ loc,s pp ) can be calculated within DMFT.
In Fig. 2(b), we show the (three maximal and four minimal) eigenvalues λ 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 Γ 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 swave superconductivity. However, as we cross the orange line at U 1 /D ≈ 1.6 in Fig. 2(a) indicating the divergence in Γ pp , one large positive eigenvalue appears in the superconducting channel. When crossing the second divergence line at U 2 /D ≈ 1.9, we can see that another large eigenvalue appears and a large negative eigenvalue disappears. Concomitantly, the eigenvalue λ 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 λ. Similarly, for the negative eigenvalues, the second lowest λ within U 1 < U < U 2 continues as the lowest λ 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 λ remains = 1. This result can be naturally understood in terms of the divergence in Γ pp . Actually, the original article [70] analyzed the eigenvalues of the generalized susceptibility χ 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 1+Γppχ 0 . In the following we show how that the superconducting phase transition in the strongly correlated regime is actually governed by the eigenvalues λ which are close to unity, rather than by the leading (i.e. largest) eigenvalue. Since the pairing interaction matrix Γ singlet,triplet pp (k, k ) 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 {∆ i } are normalized as 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 ∆ i are orthogonal).
Eqs. (4,5) clearly indicate that the divergence of λ i only connects to the divergence of Γ [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 λ → 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].

Estimation of the superconducting transition temperature from the purely two-dimensional system
For purely two-dimensional systems the Mermin-Wagner theorem [79,80,81,82]  First, let us explain how the Mermin-Wagner theorem can be understood (not proven) within diagrammatic calculations from the paramagnetic region. The self-energy Σ(k) is connected to the fully reducible vertex function F (k, k , 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 χ is the susceptibility of the relevant channel and γ 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 ≈ Q where the relevant susceptibility can be approximated as where ξ is the correlation length. Please note that for ξ → ∞ we can replace q by Q in all terms but χ(q), perform the k 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 ( q≈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. (9) indicates that, for a two-dimensional system, the selfenergy 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 selfenergy 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 c for the superconductivity long-range order [78]. In this regard, we can understand that we now obtain a finite superconducting T 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 ⊥ /J → 0, satisfying the Mermin-Wagner theorem (J and J ⊥ are the in-plane and inter-plane Heisenberg coupling constant, respectively). On the other hand, they show quite different behaviors: probably reflecting the difference divergence strength in Eqs. (9).
Indeed, the authors of Ref. [89] found a universal behavior of the transition temperature T AF c in the weak inter-layer coupling region, indicating the relation between T AF c in quasi-1D/2D system and the singular behavior in the purely 1D/2D system. In this region, the T AF c is quite stable in 2D, which changes only 10 − 20% if we change the coupling strength orders of magnitude 0.001 < J ⊥ /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 twodimensional 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 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 c remains as a quasi-long-range order even in purely two dimensions, so that the quasitwo-dimensional T c is rather stable against the inter-layer coupling.
First, we start with the discussion of the phase diagram. For cuprates, the parent materials are antiferromagnetic charge-transfer insulators (often modelled as Mott insulators in a single-band model). As we dope carriers, superconductivity appears in some doping regime and then, upon further doping, the system eventually becomes a normal Fermi liquid. Many diagrammatic extensions of DMFT have succeeded in describing this behavior since they contain the spatial fluctuation beyond the local approximation. One crucial point is that these schemes can capture the critical divergent behavior of the spin susceptibility: χ ∝ exp(−∆/T ) [45,69], like two-dimensioal antiferromagnets [101]. For capturing such a critical behavior, the effect of longrange fluctuation is necessary. Properly evaluating such a strong spin fluctuation is essential for describing the strongly correlated superconductivity. It can be a glue for the superconductivity on one side, while at the same time, it can reduce the electron spectrum (and then T c ) dramatically on the other side. Another critical point is the dome structure of T c . In Fig. 4(a), we show the DΓA result for the leading eigenvalue in the two-dimensional square lattice Hubbard model without frustration, which exhibits the dome structure [68]. While it is widely observed in cuprate materials, this dome structure was difficult to be captured by the FLEX-like diagrammatic approaches. Since we now consider the diagrammatic expansion starting from the DMFT, which takes the local strongly correlated effect, many methods have succeeded in capturing the dome structure of T c [49,95,96,68,98,99,100]. Some FLEX+DMFT studies also discuss the relation between the Pomeranchuk (nematic) instability and superconductivity [95,100]. The TRILEX study showed that the charge fluctuation becomes more relevant to the superconductivity for the triangular lattice models [102].
One remarkable feature of the vertex correction is the dynamical screening effect on Dynamical vertex structure (Matsubara frequency dependence) of the particle-particle reducible (pairing) vertex, Γ pp,Q=(π,π) (ν n , ν n , ω = 0) which includes spatial (in particular spin) fluctuations by means of DΓA. We also show a typical structure of Γ pp of mean-field-like approaches in the inset of (b). Taken from Ref. [68] the pairing interaction [68,99]. In Fig. 4(b), we show the pairing interaction mediated by spin fluctuations at Q = (π, π), which depends on two fermionic Matsubara frequencies (with fixed bosonic frequency: ω = 0). We first see the strong intensity of the diagonal structure which corresponds the spin fluctuation mediated pairing: χ sp,Q=(π,π) (ν − ν ). Additionally, we can see the strong screening effect in the low-frequency regime. This structure cannot be captured by the ordinary mean-field like (paramagnon mediated) interaction (as shown in the inset of figure 4(b)) since paramagnons depend on one bosonic frequency only. This structure is general, and we observed it for all parameters of Fig. 4(a); see also Ref. [67]. One crucial point is that vertex corrections become less relevant for high frequencies, because complicated diagrams disappear due to the decay of the Green function G ∼ 1/ω n in high frequencies [64]. Eventually the asymptotic value of the diagonal structure becomes the same as the mean-field-like value. In this sense, with fixed spin fluctuation instability, the vertex corrections suppress the intensity of the pairing vertex at low frequencies. Such a screening decreases the transition temperature by orders of magnitude since the low-frequency value is the most important for Fermi surface instabilities, such as superconductivity. These results suggest the importance of the vertex structure of the pairing interaction for the quantitative estimation of the superconducting transition temperatures.
Considering the dynamical screening effect, we can now compare computational results with experiments quantitatively. Single-band DΓA Experiment value of 15K at 20%-doping (red diamond) [5]. Furthermore, the T c dome structure around 20%-doping is also compatible with subsequent experiments [103,104,105] (please see the review article [106] for details of nickelate calculations). We note that the ignored d-wave pairing fluctuation effect and the interlayer coupling will somewhat decrease T c while we expect these are secondary (minor) effects as mentioned in Sec. 3.2.
The DΓA result was also applied to analyze the recently found quintuple-layer nickelate superconductivity [107,108]. The dynamical vertex structure [68] may also change the qualitative nature of the pairing. Usually, the instability odd-frequency pairing is weak because of the node at ω = 0. On the other hand, there is a peak at ω = 0 for the even-frequency pairing, which is favorable for a Fermi surface instability. Since the suppression of T c by the dynamical vertex structure occurs in low frequencies, it would more strongly affect the even-frequency pairing, reducing the mentioned advantage of the even-frequency pairing and relatively supporting the odd-frequency pairing. This point is further confirmed in Ref. [67], where the author analyzed the dynamical screening effect on both even/oddfrequency-pairing explicitly. Indeed, the odd-frequency superconductivity is observed with the dual-fermion approach for the Kondo-lattice model [94]. There, a similar mechanism may enter while the model and relevant fluctuation are pretty different. Regarding the dynamical structure of the vertex, it would also be an interesting question about the relation between this dynamical screening effect and studies of non (simple-)bosonic pairing glue indicated from the real frequency structures [109].

Summary and outlook
In this article, we reviewed research for superconductivity obtained by diagrammatic extensions of the dynamical mean-field theory. These approaches can simultaneously capture the effects of strong correlations and spatial fluctuations, which are both relevant for describing layered or thin-film superconductivity. As a result, the dome structure of T c observed in experiments can be described. It is also suggested that we are now able to capture the phase diagram of unconventional superconductivity, not only qualitatively but even quantitatively. An in-depth analysis shows that the dynamical screening of the pairing vertex reduce the superconducting transition temperature by one order of magnitude. We can also apply these approaches to explore the possibility of qualitatively exotic states, e.g., odd-frequency pairing.
In the future, further analyses of the interrelationship among different non-local fluctuations are desirable. Indeed, it has been gradually recognized that the charge (stripe) instability is also relevant for cuprate superconductors both experimentally and theoretically. Stripe instabilities may even dominate the superconductivity in some cases [110,111,112]. Diagrammatically, spin-fluctuation-assisted charge instability has been proposed and studied recently but rarely starting from DMFT [113,69,114]. In such systems where many fluctuations compete, taking all channels and their mutual coupling into account will become more important and is possible through the parquet equations. Another direction is the extension to multi-orbital systems to study realistic systems and explore new physics, e.g., orbital fluctuation/Hund physics. Both of these directions are performed already in high-temperature regimes [115,114,29,116]. However, they have not been applied to low-temperature phenomena like unconventional superconductivity. For extending such studies to the low-temperature regime, numerical advance is crucial (c.f. Appendix).

Appendix A. Detail formulations of high-frequencies contributions
One big problem for treating the vertex function is convergence against the number of Matsubara frequencies. For example, we usually use a few thousand points when calculating the unconventional superconductivity by using FLEX. On the other hand, we can only handle a few hundred points for the vertex. Here, we will first explain the details of the schemes used in Ref. [68] (which we just briefly explained in the Supplemental Section S.II of that paper). This is a simpler version of the pioneering work [117], but instead, can be easily extended toward beyond DMFT and seems to be enough even for describing unconventional superconductivity. Afterwards, we will discuss our results and the relations to other recent progresses.
Similar to [117], we separate the frequency range into the directly calculated range and the region supplementing asymptotic structures to accelerate the convergence against frequency box size (ν, ν , ω)-range. On top of the exactly treated region (n core ), we consider the contribution from the bare U of the irreducible vertex Γ for large size (n outer ) box and further consider the bare bubble contribution for the susceptibility within an even larger n asympt range. The advantage of this kind of procedure is that we only need to solve n outer size of the local inverse Bethe-Salpeter equation once, and other parts require only n core size calculations. Therefore we can increase n outer up to more than 1000 points, and there, the Green's function is smoothly connected to the asymptotic behavior 1/(iω n ) (i.e., we can increase n asympt → ∞).
Later, we will focus on the Bethe-Salpeter equation in one specific channel, which is enough for treating ladder extensions. Then, the equations are independent for each bosonic frequency ω, and we can assume the vertex function as the matrix. Here, we first make some rules for changing the size of the matrix. For a large matrix A, we define A x , (x = core, outer, asympt) as a smaller size (n X ) matrix as, and when we sum up matrices of different size, we supplement 0 for mismatching regions. One remark is that "the size change and matrix inverse do not commute with each other A −1 core ≡ (A core ) −1 = (A −1 ) core , and we can change these order if the matrix is block diagonalized". This is an essential point regarding the frequency box size convergence, as we will see in the following.
Appendix A.1. Local part First, we need to extract the local Γ loc core from the generalized susceptibility of the impurity solver χ loc solver after convergence of the DMFT calculation. They are related by the local Bethe-Salpeter equation as Now, we want to consider the outside bare U contribution for Γ loc . Instead of simply solving the n core -size Bethe-Salpeter equation, we need to solve the n outer -size Bethe-Salpeter equation as where U is the matrix representation of the bare interaction and δΓ is the irreducible vertex after subtracting the bare-U contribution. χ outer must be the same as χ loc solver for the core region, while we do not know outsides as shown in Figure A1. Here, we can ignore the bare bubble contribution outside n outer , since χ is diagonalized and doesn't affect Γ core (i.e., (χ loc asympt ) −1 outer = (χ loc outer ) −1 ). At a glance, this inverse BSE is not well defined due to unknown regions. Since we fix the outer region of the irreducible vertex as the bare interaction, we can obtain Γ loc core by one-shot calculation as follows. First, we defineχ loc outer as irreducible vertex fully reducible vertex ncore nouter Figure A1. Schematic picture of how to extract the irreducible vertex. With fixed bosonic frequency, the vertex depends on two Matsubara frequencies, and we map these dependencies to a plane. Here, we need to solve the inverse Bethe-Salpeter equation from the fully reducible vertex F loc (which is obtained by the impurity solver). Solid lines indicate the frequency cutoff, which is necessary to perform practical calculations. The standard way takes the orange region (n core -Matsubara points) only and ignores outsides. Our approach takes some buffer region where we only consider the bare interaction as irreducible vertex and determine Γ loc from the Bethe-Salpeter equation with n outer -Matsubara points. Nevertheless, we can extract Γ deterministically (without a self-consistent calculation) as written in the text.
We considerχ loc for the full n outer , because we do not use χ loc solver here. From Eqs. From here on, we need to construct the equation only using n core size. Focusing on the n core range for (i, j ∈ core), we obtain  .7), we use the fact that δΓ loc core have finite value only for the n core range. Then, we finally obtain the n core -size equation as (χ loc outer ) core should be equal to χ loc solver which is only defined on n core δΓ loc core = (χ loc solver ) −1 − ((χ loc outer ) core ) −1 . used for the inverse calculation). In this regards, naive expectation of Eq. (A.8) is that δΓ would be more stable against the frequency-size, sinceχ has somewhat similar box-size dependence as χ and they will be cancelled out (instead subtracting box-size independent χ −1 0 for obtaining Γ in the usual protocol). In Figure A2, we checked the current procedure for the DMFT result on the square lattice Hubbard model at U = 1.5D, n = 0.85, βD = 100. We compare particular cuts of the irreducible vertex in the magnetic channel: Γ m , for (a,d) simply solving the inverse Bethe-Salpeter equation (without supplementing bare U ) and (b,e) with supplementing bare-U contribution as Figure A1. When simply solving the inverse Bethe-Salpeter equation, the result is not well converged even for n core = 160. On the other hand, the results in (b,e) are well converged already around n core = 40. In Figure (c,f), we further check the frequency range dependence of n = 0 (lowest Matsubara) values for (a-d). We can see that convergence is pretty good when we supplement the bare-U contribution. The converged value is indeed almost the same as the extrapolated value from the result of simply solving the inverse Bethe-Salpeter equation.

Appendix A.2. Non-local part
Let us move on to the formulation of the nonlocal diagrammatic extension. For evaluating the superconductivity eigenvalues, we need to calculate the nonlocal selfenergy Σ and the ladder expanded vertex F (for obtaining the pairing interaction Γ pp from Γ pp ≡ F − Φ DMFT pp where Φ DMFT pp is the reducible part of the two-particle vertex) given as [26,68] where χ r (q), (r = d, m) is the physical susceptibility defined as χ phys Since we want to consider high-frequency bare-U contributions for Γ loc , we write this effect explicitly like χ X , γ X , χ * X and then, we explain how to calculate these quantities with considering asymptotic behavior with small (n core ) size vertex calculations. Hereafter, we omit some subscripts (ω, d, m), since they are independent in the particle-hole channel and the following calculations apply equally for all degree of freedom.
Eq. (A.17) does not obey the well known asymptotic behavior γ ν → 1 as ν → ∞ [45]. This is because we consider the infinite range of bare-U interaction in Eq. (A.12) while we take n outer range for the whole formalism. Nevertheless, there is no ambiguity and this just comes from how to express the same quantity. Please note that γ appears as γ(1 − U r χ) in the self-energy calculation and γ asympt (1 − U r χ asympt ) = γ outer (1 − U r χ outer ), γ outer → 1. Then the current expression is fully consistent even if considering up to n outer range, and the exact asymptotic behavior holds: γ outer → 1.
On the other hand, we expect that χ asympt converges slightly faster, and the current expression is better for considering the physical susceptibility itself (for example, when applying sum-rule for them [45]). Finally, let's look at how the current formalism works for evaluating unconventional superconductivity. We show the numerical result employing the current scheme in Figure  A3 for βD = 400 [68]. Please note that it is impossible to get well converged results with simple implementation for such low temperatures. We can see that we can increase n outer around several thousand, and, thus, we only need a few hundred Matsubara grids for the direct treatment of the vertex.

Appendix A.3. Discussion
Figures A2 and A3 indicate that such the simple framework schematically vizualied in Fig. A1 efficiently accelerates the convergence. Indeed, figure A3 shows that the size of bare-U itself is much more important than the rest of the detailed asymptotics for Figure A3. Matsubara frequency cutoff (n core , n outer ) dependence of leading superconductivity eigenvalue in the square lattice Hubbard model for (a) n = 0.825, t /t = t /t = 0 and (b) n = 0.90, t /t = −0.20, t /t = 0.16. The purple lines show n outer -dependence with fixed n core = 80 and the blue lines show n core -dependence with fixed n outer = 1024 for (a) and n outer = 2048 for (b). Horizontal lines are guides to the eyes for the almost converged result at n core = 159. Taken from Ref. [68] capturing the high-frequency contributions. This fact can be understood as follows: For high frequencies, the fermion-boson three points vertex is reduced to the bare interaction, and in such regions, simply the summation of the vertex is an important quantity, for example, a part of the ladder: γχ 0 Γχ 0 γ (c.f., [118,119]) becomes Uχ 0 Γχ 0 U = U 2 ν,ν χ ν 0 Γ ν,ν χ ν 0 as γ → U. For irreducible vertices, many asymptotic structures persist but at a specific frequency. In the most case, the bare-interaction contribution dominates if we sum up against Matsubara frequencies except that the fluctuation in another channel diverges. Therefore, we expect supplementing the bare interaction for the irreducible vertex works efficiently at least for the most relevant channel.
There are many other works for treating the complicated vertex more efficiently. In addition to the pioneering work [117], Ref. [66] employed the vertex at large frequency to extrapolate it to its asymptotics in parquet calculations, Ref. [120] studied the effect of supplementing the vertex asymptotics [64] for all channels at the local level. Ref. [121] tries to use the analytical expression to take the large size limit of the outer-box. There has been some progress recently [122] solving the Bethe-Salpeter equation in the intermediate representation [123,124]. Also a compactification of the vertex in momentum space has been suggested using a form factor expansion [125]. Other than these, many kinds of the vertex decomposition have been proposed [126,127,128]. Ref. [129,130] extended the current (supplementing bare-U ) scheme to the parquet equations by using single-boson-exchange (SBE) decomposition [131] based on a concept of the irreducibility against the bare interaction [45]. Please note that SBE graphically has a similar structure as Eq. (5), and clarifying the relation between them would be important. Indeed, for the matrix, eigenvalue decomposition (more generally the singular value decomposition) is straightforward and almost equivalent to the solving the Bethe-Salpeter equation (c.f., Eqs. (4,5)), while it is not trivial for the relation between solving parquet equations and higher rank tensor decompositions of the irreducible vertices. A big step is still an efficient implementation for the parquet equations where memory constraints become even more relevant than for the ladder calculations and which likely, requires some further physical insight.