Pairing in the two-dimensional Hubbard model from weak to strong coupling

The Hubbard model is the simplest model that is believed to exhibit superconductivity arising from purely repulsive interactions and has been extensively applied to explore a variety of unconventional superconducting systems. Here we study the evolution of the leading superconducting instabilities of the single-orbital Hubbard model on a two-dimensional square lattice as a function of onsite Coulomb repulsion U and band ﬁlling by calculating the irreducible particle-particle scattering vertex obtained from dynamical cluster approximation (DCA) calculations, and compare the results to both perturbative Kohn-Luttinger (KL) theory as well as the widely used random phase approximation (RPA) spin-ﬂuctuation pairing scheme. Near half-ﬁlling, we ﬁnd remarkable agreement of the hierarchy of the leading pairing states among these three methods, implying adiabatic continuity between weak- and strong-coupling pairing solutions of the Hubbard model. The d x 2 − y 2 wave instability is robust to increasing U near half-ﬁlling as expected. Away from half-ﬁlling, the predictions of KL and RPA at small U for transitions to other pair states agree with DCA at intermediate U as well as recent diagrammatic Monte Carlo calculations. RPA results fail only in the very dilute limit, where it yields a d xy ground state instead of a p -wave state established by diagrammatic Monte Carlo and low-order perturbative methods, as well as our DCA calculations. We discuss the origins of this discrepancy, highlighting the crucial role of the vertex corrections neglected in the RPA approach. Overall, a comparison of the various methods over the entire phase diagram strongly suggests a smooth crossover of the superconducting interaction generated by local Hubbard interactions between weak and strong coupling. results are interpolated to the Fermi surface of the noninteracting system to allow a direct comparison although the DCA calculation is only done for 8 × 8 discrete values.


I. INTRODUCTION
Since the theoretical proposal by Kohn and Luttinger (KL) [1,2] that superconductivity can arise from purely repulsive electron interactions and the subsequent discovery of superconductivity in materials like heavy fermions, cuprates, organic Bechgaard salts, and iron-based superconductors, superconducting instabilities in models of interacting fermions have been extensively studied. The Hubbard model [3] has played an exceptional paradigmatic role in this discussion. It is the simplest model of fermions with local interactions, and was argued furthermore to be the appropriate effective model to describe unconventional superconductivity in correlated electron systems, notably in cuprates [4]. The model has also been popular because the physics of pairing by spin fluctuations, originally suggested by Berk and Schrieffer [5,6] in continuum models and extended by Scalapino and others to lattice Hubbard-type models [7][8][9], is rather simple to capture within the straightforward random phase approximation (RPA) [10]. At present, theoretical studies of the Hubbard model constitute a growing research area as seen, for example, by several recent extensive comparisons between various state-of-the-art numerical methods providing updated benchmarks on the ground-state energy, the self-energy, and competing order in the Hubbard model [11,12]. An important next step is to compare and quantify the superconducting pairing instabilities within this model [13].
Close to half band filling, the 2D Hubbard model on a square lattice is known to exhibit strong d-wave pair correlations [14]. While a rigorous proof that d-wave superconductivity exists in the model at T = 0 is lacking, the preponderance of evidence from numerical calculations [15][16][17][18][19][20][21][22][23][24][25][26][27], as well as weak and intermediate coupling renormalization group studies [28][29][30][31][32][33][34][35], strongly support this conclusion. Rather less is known with high confidence at larger interaction U , further away from the half-filled state, or with regard to subleading pair channels throughout the phase diagram. It is convenient to study the latter two questions using controlled perturbative methods or via RPA due to physical transparency and ease of implementation. Several authors, including the current ones, have applied various weak-coupling schemes to map out the leading superconducting instabilities as a function of, e.g., doping and band parameters, displaying a rich mosaic of different pairing states [36][37][38][39][40][41][42]. Predictions of these studies for leading pairing instabilities throughout the phase diagram appear to agree rather well with recent diagrammatic Monte Carlo calculations that should be wellcontrolled and able to reach somewhat higher U [24]. However, the general question of how the pairing in the 2D Hubbard model changes as correlations are increased is still open.
In the special case close to half-filling, the Hubbard model reduces to the t-J model as U → ∞, and t-J studies have also found d-wave pairing in this doping regime [43][44][45][46][47][48][49][50]. Numerical results in Ref. [48] showed that the pairing interaction in the Hubbard and t-J models contains both a large retarded contribution arising from spin fluctuations and a smaller nonretarded contribution. For the Hubbard model, the "nonretarded" contribution arises from high-energy excitations involving the upper Hubbard band, while for the t-J model, there is a corresponding instantaneous exchange contribution. Moreover, this study found that the retarded spinfluctuation interaction provides the dominant contribution to the pairing interaction. This supports the intriguing conclusion that the physics of pairing at strong coupling is similar to, or at least evolves continuously from, that at weak coupling. However, away from half-filling little is known about this crossover.
In this study, we investigate if the evolution of the pairing interaction is continuous from weak to strong coupling. This is reflected in the continuity of the leading and subleading pairing instabilities upon increasing Coulomb interaction U . It is not our intention here to make detailed comparison of the relative stability of the various competing orders that are observed in cuprates and may be present in the Hubbard model or its generalizations, nor to explain the formation of the pseudogap state out of the Mott insulator at half-filling. Instead, we simply propose to understand the evolution of the pairing interaction in the Hubbard model as a function of doping and correlation strength.
With this question in mind, we calculate the superconducting pairing vertex of the Hubbard model via numerical solutions of the Bethe-Salpeter equation obtained in the dynamical cluster approximation (DCA) with quantum Monte Carlo (QMC) impurity solver. This approximation [51,52] is known to provide an accurate estimate of the pairing vertex over a range of intermediate strength U values appropriate for the cuprates [53]. We then compare these results with both perturbative Kohn-Luttinger (KL) theory and the RPA scheme. The latter breaks down at higher U due to the wellknown magnetic instability inherent to the approximation, but is thought to work well at smaller interaction strengths. We find that the hierarchy of pairing eigenvalues of the linearized gap equation match up well between the methods in the region where they can be compared, providing further evidence that the pairing evolution is smooth from weak to strong coupling. Results from the simple RPA agree spectacularly well with diagrammatic Monte Carlo [24] over the entire doping range except at the smallest doping, where p-wave spin triplet pairing is stable over a much narrower region in the RPA than obtained in asymptotically exact results [54,55]. We discuss the reasons for this discrepancy.

II. MODEL AND METHODS
We study the single-orbital Hubbard model defined on a 2D square lattice where c † iσ /c iσ creates/annihilates an electron at lattice site i with spin σ , and n iσ = c † iσ c iσ is the number operator of electrons with spin σ at site i. The nearest-neighbor hopping sets the energy unit, t = 1, and we include also next-nearestneighbor hopping t . The superconducting pairing originates from the repulsive Coulomb interaction and is treated numerically by three different approaches. First, we calculate the full energy-resolved pairing kernel with inclusion of selfenergy corrections, and solve the Bethe-Salpeter equation with Green's functions and irreducible particle-particle vertex obtained from quantum Monte Carlo simulations using the dynamic cluster approximation (DCA). The details are explained in Sec. II A below. Second, we apply perturbative KL theory, and lastly compare with the RPA spin-fluctuation method for pairing, both detailed in Sec. II B. In KL-theory, only the second-order diagrams enter the pairing theory, whereas in the RPA approach, the effective pairing interaction is evaluated diagrammatically by a selected class of diagrams that highlights the physics of nesting and pronounced spin fluctuations. Whereas KL theory is a controlled weak-coupling approach valid at small interactions U (compared to the bandwidth), the solution of the Bethe-Salpeter equation by DCA is not restricted to a certain regime of Hubbard U . However, this method is significantly heavier computationally and suffers from the sign problem [51]. This imposes constraints on the smallness of U as well as the size of t , doping, and cluster size [51,52].

A. Pairing within the dynamical cluster approximation
For the quantum Monte Carlo calculations, we use a dynamic cluster approximation [51] with a continuous-time auxiliary field (CT-AUX) QMC solver [52]. The DCA represents the bulk lattice by a finite size cluster and uses coarsegraining of the reciprocal space to retain information about the remaining bulk degrees of freedom. Within this cluster approach, the first Brillouin zone is divided into N c patches P K , each of which is represented by a cluster momentum K, and within which the self-energy (k, iω n ) is assumed to be constant and given by the cluster self-energy c (K, iω n ). One then averages the Green's function over the patches P K to determine the coarse-grained Green's function Here the sum is restricted to the N c /N momenta within the patch about the cluster momentum K. The corresponding bare propagator G 0 (K, iω n ) = [Ḡ −1 (K, iω n ) + c (K, iω n )] −1 is then used together with the interaction U to set up the effective cluster problem, in which the self-energy c (K, iω n ) = F[G 0 (K, iω n ), U ] is calculated with the CT-AUX QMC solver. This calculation is repeated iteratively until the selfenergy has converged. For further details, the reader is referred to Ref. [51]. After convergence, the two-particle Green's function in the particle-particle channel with zero center of mass momentum and energy, G c,2 (K, K ) = G ↑↓↓↑ c,2 (K, −K, −K , K ) is calculated for the cluster problem [53]. Here, K = (K, iω n ) and K = (K , iω n ). The irreducible particle-particle vertex pp (K, K ) is then extracted from the Bethe-Salpether equation and used in the DCA gap equation for the bulk lattice where the pairing kernel G(k, iω n )G(−k, −iω n ) has been coarse-grained over the momenta of the DCA patches P K to giveχ pp The solution of this eigenvalue equation gives the DCA results for the leading eigenvalues λ α and corresponding eigenvectors φ α (K ). We use a cluster size of N = 64 for U = 2 and N = 32 for U = 4, 6, and 8 and the temperature is set to T = 0.025, 0.05, 0.15, and 0.2 for U = 2, 4, 6, and 8, respectively. We stress that the DCA calculations in the intermediate to strong coupling regime (U = 6 and 8) were carried out at temperatures T 0.15. While there is mounting numerical evidence for static charge and spin stripe order in the ground state of the doped Hubbard model [56][57][58][59], particularly at 1/8 doping, fluctuating stripe correlations rather than static order was found at the elevated temperatures we have used in recent determinant QMC calculations [60]. We therefore ignore the possibility of competing static stripe order in the DCA calculations. Moreover, we only consider the on-site U repulsion. An additional nearest neighbor repulsion V term in the extended Hubbard model was considered in recent studies using DCA calculations [61] as well as cellular dynamical mean-field theory [62,63]. These calculations found that dwave superconductivity is only weakly suppressed by small to moderate V because of the strongly retarded nature of the d-wave pairing interaction.

B. Pairing within Kohn-Luttinger and RPA spin-fluctuation theory
In both weak-coupling KL theory as well as in RPA spin-fluctuation mediated superconductivity, one derives an effective Cooper pair term of the form with V (k, k ) denoting the effective pairing vertex. In the KL approach, the vertex is obtained to second order in U . Since Hubbard interactions connect propagators of opposite spin only, this amounts to an evaluation of the diagrams depicted in Fig. 1. The effective interaction obtained from the second order diagrams is given by where the upper (lower) sign corresponds to the singlet (triplet) channel. In the singlet channel, the effective interaction includes the first order repulsive term U . The bare spin susceptibility is given by the Lindhard function evaluated at zero energy Within RPA, the screening (bubble) and exchange (ladder) diagrams depicted in Fig. 1 are summed to infinite order in U . The bubble diagrams correspond to effective interactions through longitudinal fluctuations, while the ladder diagrams are due to exchange interactions mediated by transverse fluctuations. The final interaction between opposite spin electrons is As shown in Eq. (8), the contribution from screening can be written in terms of a spin-fluctuation term displaying the Stoner enhancement and a charge contribution, which generally remains small. In this study, we restrict ourselves to the paramagnetic phase, where the longitudinal and transverse spin susceptibilities are the same and all triplet channels are degenerate.
We calculate the spin-singlet (s) (spin-triplet (t )) gaps by symmetrizing (antisymmetrizing) the effective interactions, The superconducting gap equation with E k = √ ξ 2 k + | k | 2 is linearized by setting E k = |ξ k |. This gives the leading and subleading superconducting instabilities at T c and amounts to a calculation of the eigenvalues λ i and corresponding eigenvectors g i (k) of the matrix

III. RESULTS
We focus first on a case near half-filling with electron density n = 0.90 and nearest-neighbor hopping t = −0.15 and explore the role of increasing U from weak to strong coupling. The associated noninteracting Fermi surface and bare susceptibility χ 0 (q) featuring pronounced (π, π )-centered fluctuations are shown in Fig. 2 Table I and plotted in Fig. 3. Here, we limit the discussion to even-frequency solutions, but note that subleading odd-frequency solutions also exist. In contrast to DCA, both the KL and RPA schemes are limited to small values of the Coulomb interaction of U = O(t ), and RPA is additionally sensitive to the inherent magnetic instability that occurs upon increasing U . As seen from Table I and Fig. 3, in all the DCA cases the leading solution is the lowest order d x 2 −y 2 solution with four nodes along the Brillouin zone diagonals. Upon increasing U , the subleading DCA instabilities approach the leading d (4) instability as inferred from Fig. 3, and additional nodal singlet solutions appear in-between the d x 2 −y 2 state and the highest triplet state denoted p in Table I and Fig. 3. The number of gap nodes resolved at the Fermi level is sensitive to the cluster size; at U = 2 which is calculated for a cluster size of N = 64, twenty nodes are resolved at the Fermi surface for the third subleading d x 2 −y 2 DCA solution. In comparison, the third subleading solution at U = 4 exhibits only four nodes. While this could be a real effect due to the increased interaction strength, it may also simply be due to the smaller cluster size of N = 32. For U = 6, a larger number of subleading solutions appear, many of the same nodal structure, e.g., the d x 2 −y 2 solutions with twelve nodes at the Fermi surface, denoted d (12), which are distinguished by a change of spectral gap weight at different   (8)], A 2g [sin k x sin k y (cos(k x ) − cos(k y )) denoted by g (8)], B 1g [cos(k x ) − cos(k y ) denoted by d (4) and higher order (cos(2k x ) + cos(2k y ))(cos(k x ) − cos(k y )) denoted by d (12)] and the E u [(cos(k x ) − cos(k y )) sin(k x ) denoted by p (6), but with nodes displaced slightly away from the zone diagonal]. For simplicity, we only state instabilities that appear before the leading triplet solution (with the exception of U 0.1). In the last s solutions of the U = 8 column, the number of nodes is undecided due to system size limitations. (12) g (8) d (12) g (8) g (8) g (8) g(8) g (8) g (8) p ( parts of the Fermi surface. As mentioned above, the DCA result for U = 8 was obtained for t = 0 and T = 0.2 (as opposed to t = −0.15 and T = 0.15) in order to avoid the QMC sign problem. Thus, comparing the U = 8 result to the lower U cases, there are additional quantitative effects in the eigenvalues from this change of parameters. We do not, however, expect such changes to have a large effect on the  pairing instability for large U , and indeed the hierarchy of solutions remains unaltered as seen from Fig. 3.
Turning next to the KL and RPA results for the same band, the hierarchy of the leading pairing solutions are displayed also in Table I and Fig. 3. As expected for this filling, both methods predict a leading d x 2 −y 2 state and agree on the hierarchy of the subleading solutions at the lowest U . In contrast to DCA, however, in the low-U limit, the triplet solution denoted p becomes the second leading instability. This is due to the proximity to the van Hove singularity, which is known to enhance the effective triplet pairing interaction [39]. For t = −0.15, the critical density for which the van Hove saddle points reside at the Fermi surface is n van Hove = 0.875 and at n = 0.90 we are thus not far from this regime.
Upon increasing U (but still within the RPA regime), the p solution rapidly drops and, as seen from Table I and Fig. 3, there is excellent agreement with the DCA pairing hierarchy near the regime of U = 1. At larger U , still within the RPA calculation, the magnetic instability is approached, and the two nearly degenerate instabilities d (12) and g (8) are interchanged. This is an artifact of the RPA approach which can be understood in the following way. The g (8) solution increases more steeply because it takes full advantage of the strongly enhanced susceptibility. Unlike d (12), the g(8) solution does not have any nodes at the Fermi surface segments in the large gap regions close to (π, 0) and symmetry-related points. The symmetry-imposed nodes of the g (8) solution along the zone axes do not inhibit gap formation, since the Fermi surface segments do not close at (π, 0) and symmetryrelated points. Nevertheless, except from such caveats as (over)sensitivity to band details or magnetic instabilities, the overall evolution of the leading superconducting solutions as discussed here highlights the agreement of the methods, and points to an adiabatic continuity between weak-and strongcoupling pairing solutions of the Hubbard model near the half-filled regime.
Next, we compare the detailed properties of the gap solutions obtained by DCA to the results of RPA (at U = 1). In Fig. 4, the three leading instabilities of both approaches are displayed. As seen, there is remarkable agreement between the two methods, giving in both cases a leading d x 2 −y 2 solution with four diagonal nodes, i.e., d (k) = 2 [cos(k x ) − cos(k y )], but with strong gap enhancements around (0, ±π ) and (±π, 0). In RPA, this enhancement is caused by the large density of states present in those regions of k space due to the proximity of the van Hove singularity. A similar effect at different doping levels was discussed in Ref. [39]. In the linearized gap equation, this effect enters via the inverse of the Fermi speed in the matrix elements of Eq. (10) as well as in the amplification of the bare spin susceptibility. However, the enhancement is also found in the DCA approach where we interpolate the solution to the Fermi surface of the noninteracting system, see from Fig. 4(d). For U = 2, we expect the Fermi surface to be very similar to that of the noninteracting system. From Fig. 4(d), we see that the gap enhancements are robust towards the inclusion of self-energy effects.
The two sub-leading solutions shown in Fig. 4 consist of a d x 2 −y 2 -wave solution with twelve nodes d (12), and a lowestorder g-wave state with eight nodes g (8). These solutions are close in energy and are both strongly suppressed compared to the leading d x 2 −y 2 gap with four nodes. For the second and third leading solutions, the RPA approach produces strong gap enhancements at the Fermi surface points closest to (0, ±π ) and (±π, 0) whereas this effect is less pronounced in the DCA  [24] shown by open symbols and lines. The blue line displays the phase boundary between a simple harmonic p wave [sin(k x )/ sin(k y )] and d xy , and the dark green line indicates the phase boundary between d xy and d x 2 −y 2 . The yellow region indicates an s phase [40]. A region of higher-order p wave (sketch as inset [24]) is visible in both approaches for small values of U around a filling of n = 0.55. The triplet p state features six or ten nodes depending on temperature. In RPA, the region of leading triplet instability at very low filling is sensitive to resolution and temperature; it essentially vanishes as T is decreased to T = 0.0002 (small blue arrows and blue dashed line almost coinciding with the y axis). (b) Zoom in of the phase diagram close to n = 0.55 and comparison to Ref. [24]. The insets show solutions g(k) as obtained from RPA at the parameters marked by white points. The dashed red line (and the small red arrows) shows how the boundary between the p and d xy changes within RPA when T is lowered from T = 0.015 to 0.0002. calculations, especially for the third leading g-wave solution. This may arise from the fact that the g-wave solution has nodes along the zone axes and DCA, with fewer k points to sample the Brillouin zone as shown in the inset of Fig. 4(f), therefore does not capture the enhancement effect for this solution.
Encouraged by the overall agreement between KL/RPA and DCA, we next compare the RPA results obtained here with previous reports in the literature. In Fig. 5, we show a direct comparison of the ground-state phase diagram of Ref. [24] obtained from diagrammatic Monte Carlo simulations and our RPA calculations. As seen, there is qualitative and in most cases nearly quantitative agreement between the different methods. For example, upon hole doping the instability from the d x 2 −y 2 state to the d xy state occurs almost simultaneously in the two methods. Furthermore, both KL/RPA and diagrammatic Monte Carlo simulations find a p -wave triplet state to become leading for small values of U around a filling of n = 0.55.
Most often, triplet solutions become favorable when the system approaches a van Hove singularity regime [26,39,64]. To second order in U , the pairing vertices are given by In the absence of a q = 0 peak structure in the susceptibility, the triplet pairing cannot take advantage of the attractive contribution to the pairing kernel − U 2 2 χ 0 (k − k ). Usually, such a peak is what renders the triplet solution favorable in the vicinity of a van Hove instability. However, the region of triplet superconductivity evident from Fig. 5(a) around n = 0.55 at the smallest U has a different origin (at t = 0 the van Hove singularity occurs at the Fermi level for filling n = 1).
Around n = 0.55, the system is in an interesting crossover regime, where the spin susceptibility shows prominent features at nesting vectors Q (π, ± π 2 )/(± π 2 , π ), which lie right in-between nesting vectors along the zone diagonal and zone axes, driving the singlet d x 2 −y 2 and d xy solutions, respectively. The odd parity p solution most optimally accommodates this nesting structure, but since it is not supported by a q = 0 peak, it is rather fragile and becomes rapidly suppressed as U increases. This can be understood from the fact that the spin susceptibility exhibits extended ridgelike structures which are most pronounced around (±π, 0)/(0, ±π ). Upon increasing U , these ridges will dominate the effective pairing and drive the system into the singlet d xy solution.
In Fig. 5(b), we show a zoom-in of the phase diagram relevant to the triplet p phase. As seen, the phase boundary of the p phase to the d-wave states agrees remarkably well with the diagrammatic Monte Carlo calculations by Deng et al. [24] at the lowest temperatures. The detailed gap structure of the superconducting p state features ten nodes as shown by the inset in Fig. 5(b). This is slightly different from the illustration shown in Fig. 5(a) from Ref. [24], but similar to the gap structures discussed in Ref. [26].
We end with a brief discussion of the pairing instabilities in the low-density regime n 0.3. As seen explicitly in Fig. 5(a), the low-density limit hosts a triplet p-wave superconducting phase. We stress that this is a standard two-node sin(k x )/ sin(k y ) p-wave state distinct from the p triplet state discussed above. The possibility of a transition from d x 2 −y 2 to d xy or p-wave superconductivity in the low-density regime of the weakly repulsive 2D Hubbard model was discussed early on by Baranov and Kagan [65] and by Chubukov and Lu [54] who analyzed the behavior of the pairing vertex in symmetrydistinct pairing channels as a function of band parameters. The more recent diagrammatic Monte Carlo calculations by Deng et al. [24], mapped out the phase boundaries between the p-, d xy -, and d x 2 −y 2 -wave pairing solutions in the low-density and low-U limits, reproduced in Fig. 5. As seen from Fig. 5, even though there is substantial overall agreement with the RPA results, the low-density regime stands out as exceptional in this comparison between the methods. At the lowest T and in the limit n, U → 0 the preferred state from the RPA study is d wave with near degeneracy between d xy -, and d x 2 −y 2 -wave pairing. At larger U , however, as seen from Fig. 5, RPA does not capture the preference for p-wave pairing in the dilute limit. This result, however, is not surprising since the lowest-order diagrams included in the RPA procedure are known to not capture the tendency for p-wave pairing in the dilute limit. Only by including higher-order vertex renormalizations does p-wave pairing get supported. This was shown initially by Chubukov who analyzed the thirdorder diagrams for renormalization of the fermionic scattering amplitude in 2D, and found that the vertex renormalization in the particle-particle channel is crucial for realizing the pwave state at low densities [55]. Subsequent studies confirmed the importance of O(U 3 ) vertex corrections for stabilization of p-wave pairing at low density [66,67]. As a consistency check we applied the DCA machinery to calculate the leading instability at U = 4, t = 0, and T = 0.0125 at fillings n = 0.15 and n = 0.20 (due to resolution we cannot address lower n by this method). In the first case (n = 0.15), we obtained indeed a leading p-wave solution even for U = 4, while d xy is the preferred state at n = 0.20. This points to a rough agreement with the phase boundary obtained by diagrammatic Monte Carlo simulations in Ref. [24], and again suggests a smooth crossover of superconductivity from weak to strong interactions.

IV. CONCLUSIONS
While there is general agreement that the leading Cooper pairing instability of the Hubbard model close to half-filling is the d x 2 −y 2 -wave state, and work on the t-J model valid in this regime corresponding to very large U suggests the same, rather less is known consensually about the rest of the Hubbard model pairing phase diagram, including fillings far from n = 1 and intermediate to strong U . These regimes are not simply of academic interest, but may well represent reasonable descriptions of a variety of unconventional superconductors, including cuprates, organic Bechgaard salts, heavy fermion materials, iron-based superconductors, and ultracold fermionic gases. In this work, we have compared different approximate methods, expected to be valid in different correlation regimes, to predict the leading and subleading superconducting instabilities in these less-studied situations. We find that spin and charge fluctuation exchange pairing calculated from both KL (small U ) and RPA methods, and a DCA Quantum Monte Carlo approach (intermediate to strong U ) compare rather favorably to each other, suggesting a smooth crossover in pairing states within the Hubbard model from weak to strong coupling at all fillings. Our results compare well to recent diagrammatic Monte Carlo calculations, with the exception of the regime of very small electron density where weak-coupling approaches need to be cured by vertex corrections. The agreement with RPA allows for a transparent explanation of the physics of several of these less wellknown pairing phases. Clearly the hypothesis of adiabatic connectivity of pair states from weak to strong coupling needs further scrutiny and investigations by other methods capable of handling electron pairing in the strongly correlated regime.