Superconductivity of repulsive spinless fermions with sublattice potentials

,

There have been substantial efforts [1][2][3][4] to understand superconductivity mechanisms beyond the conventional phonon-mediated [5] electron-electron attraction.In one category of mechanisms, bare repulsive electron-electron interaction becomes effectively attractive due to virtual processes after projections to the sublattice or bands [6][7][8].Recently, exact results for an effective attraction have been obtained for fermionic honeycomb lattice models with a large staggered sublattice potential [9][10][11].This mechanism can be essentially captured by a minimal model of spinless fermions [9], of which the low-energy physics projected to one sublattice shows effective attraction.Such a mechanism has been argued to be relevant for triplet pairing in materials [10,[12][13][14].
In this Letter, we study the pairing of spinless fermions on the square lattice in addition to the honeycomb lattice model studied in Ref. [9].Studying a different lattice can shed light on the relevance of the proposed pairing mechanism to layered materials, in which different lattice structures can be realized [15].Considering a different lattice contributes to further understanding the ingredients of the sublattice projection mechanism for superconductivity-and, as we show, reveals qualitatively different possibilities.The effective theory from a sublattice projection depends on the coordination number of the lattices; lattice symmetry is crucial for the realization of different types of unconventional superconductivity [16][17][18][19][20].
The overall result is summarized in Fig. 1.The quantum phases are inferred through infinite density matrix renormalization group (DMRG) [21,22] data for strong coupling combined with functional renormalization group (FRG) [23] data at weak coupling.Superconducting phases are found in a wide range of interaction parameters in the honeycomb model while its regime is limited to smaller interactions for the square model.Compared to a previous study [9], a significant difference is that there are two superconducting phases on the square lattice, the staggered d-wave and the p+ip topological phases, in contrast to the sole f -wave pairing on the honeycomb lattice.The d-wave pairing on the square lattice shares the same origin as the f -wave pairing on the honeycomb lattice in the sense of inter-valley pairing.The Cooper pair arises from an inter-valley attraction revealed by sublattice projection.This requires a next-nearest-neighbor hopping t to realize a two-valley band structure for the square lattice.Upon increasing doping, we observe a transition from staggered d to a topological p + ip [24] superconductor.With zero momentum, p + ip no longer results from the inter-valley attraction.It does not require the next-nearest-neighbor hopping.Moreover, at stronger interactions, we find evidence for a transition from superconductivity to inhomogeneous states.
Model and low-energy description.-Weuse the square lattice [Fig.1(ai)] as an example while the honeycomb model [Fig.1(aii)] can be found in Ref. [9] with the same form of Hamiltonian.The Hamiltonian is taken as where c i , (c † i ) is the fermionic annihilation (creation) operator on site i, and n i = c † i c i .The symbols i, j and i, j denote nearest neighbors and next-nearest neighbors, respectively.We limit our attention to repulsive arXiv:2208.02830v2[cond-mat.supr-con]26 Jan 2023 interaction V > 0 and sublattice potential D |t| > 0 on the sublattice B. At half filling and large D, the ground state is expected to have the A sublattice fully filled and the B sublattice unfilled.When t = 0, the Hamiltonian exhibits an explicit symmetry of particlehole transformation c † A → c A and c † B → −c B combined with spatial inversion that interchanges the sublattices.When t = 0, the combined particle-hole transformation equivalently changes the sign of t .In this work, we will only introduce t = 0 on the square lattice while t = 0 on the honeycomb lattice, which is motivated by the discussion below.
We focus on electron doping the system near half filling, where the low-energy physics is controlled by those extra electrons on the B lattice.The effective model is derived by a Schrieffer-Wolff transformation [25,26] (for details see the Supplemental Material [27]).Up to the second order of t, this effective Hamiltonian contains terms of hopping, correlated hopping, and interactions: Different parts of the Hamiltonian are introduced as follows (for the details of the coefficients see [27]): H hopping contains nearest neighbors ij and next-nearest neighbors ij terms for the sublattice B: where t B = 2λ 0 − t and t B = λ 0 with λ 0 = t 2 /(D + 2V ).
For most of our calculations, we will either fix t = 0 or t = λ 0 .The correlated hopping also includes two terms The combinations ijk and ijkl are summed over all possible ordered vertices of plaquettes in sublattice B, e.g., 1,2,3 and 4 in Fig. 1(ai).Finally, there are two-, three-, and four-body density interactions, The combinations [ijk] and [ijkl] are summed over all possible unordered vertices of plaquettes in sublattice B. The four-body interaction U 4 remains repulsive in the full parameter region, while other interaction terms turn from repulsion to attraction when increasing across The dispersion of the kinetic part H hopping depends on the next-nearest-neighbor hopping t .At t = 0 [shown in Fig. 2 (a)], the band minimum is located along the boundary of the Brillouin zone.The Fermi surface is connected and has an approximate rotation symmetry.By tuning t such that |t B /t B | > 0.5, two band minima appear at (0, ±π) and (±π, 0), respectively, where the unit of the wave vectors is 1/a.The low-energy physics is then controlled by these two valleys which are interchanged under a π/2 rotation.When tuning to higher doping, the Fermi surface includes the Van Hove singularities.They are located at (q, ±q) with q = ± arccos(−t B /(2t B )).The two-valley low-energy physics is replaced by the one exhibiting new instabilities driven by the larger density of states.We remark that introducing t on the honeycomb lattice only brings an overall factor to the band dispersion.Two-valley continuum theory of the square lattice model.-To construct the continuum theory in the case with two valleys, the degrees of freedom for doped electrons can be decomposed into two valleys: c j = σ a exp[iK σ • r j ]ψ σ (r j ) with K + = (0, π) and K − = (π, 0), where the fields ψ σ (r) vary slowly at the scale of a, the minimal distance between two B-sublattice sites.
At low doping, we ignore the three-and fourbody interactions in H U .The continuum Hamiltonian includes a kinetic part with anisotropic masses y ψ σ /(2m yy σ ) at two valleys and a two-body interaction term.There are two contributions to the two-body interaction in the continuum limit, the correlated hopping terms in Eq. ( 3) and the two-body repulsion terms Eq. ( 4).In the long-wavelength limit, the interaction can be written as where , where ∆ is odd under a π/2 rotation.While the pair has non-zero total momentum, the above reasoning for the pairing is the same as that for f -wave SC of low-doping honeycomb model [11].For finite doping, realizing pairing with (π, π) center-of-mass momentum is frustrated by the shape of the Fermi surfaces.This could lead to a transition to incommensurate (not observed) or other SC phases.Inferring the possible SC at finite doping from the bare Hamiltonian of the projected model is no longer simple.
The complication comes from the interactions projected on the Fermi surface.Nevertheless, we can show that for the intra-valley interaction, the correlated hopping in the projected model can induce bare attractive interaction term between pairs of fermion modes on the Fermi surface with zero net momentum, for details see the the Supplemental Material [27].Thereby the possibility of intra-valley pairing, likely p-wave pairing, is suggested.We will later discuss the role of Van Hove singularity for SC, which is independent of the role of projected interactions.
Binding energies for few-particle doping.-Next, we show our numerical results of pair and larger bound states formation in the dilute doping limit.Binding energies can be deduced from the difference between oneparticle doping energy and energy per particle of nparticle doping; the data for the effective model (D/t = ∞) are plotted in Fig. 3. (Our data for D/t = 5, 10 can be found in Ref. [27].)From the data, we can infer that at D/t = ∞, there can be a stable dilute pairing phase for the honeycomb lattice with V /D 1.The pairing phase is not favored for the square lattice with t = 0, but it can exist with t > 0. For t = λ 0 , the condition for pairing phase is V /D 0.6.
We also determine the momenta of the few-particle ground states.The momentum of a pair for the square lattice with t = λ 0 and the honeycomb lattice, are respectively (π, π) and (0, 0).Recall that two valleys of the honeycomb lattice are located at ±K (standard notation [28]), and those of the square lattice model are located at (π, 0) or (0, π).This along with finite pair binding energy results indicates an inter-valley pairing mechanism and explains the absence of it in the case of t = 0 with the absence of valleys.The two-valley structure allows stable pairing, for which a sufficient attraction between fermions in different valleys exists but no attractions sufficient for larger bound states.The latter condition can be usually met with weak coupling as the intra-valley coupling is less relevant in the dilute doping limit.
Numerical study of the quantum phases at finite doping.-Theabove few-body and continuum theory results provide an indication of superconductivity at low doping and its instability for large interactions.In the following, we apply DMRG and FRG to infer the quantum phases of the full models with D = 10 and t = λ 0 (square) and t = 0 (honeycomb) at finite fermion doping from weak to strong coupling (details see the Supplemental Material [27]); the results are summarized in Fig. 1.The d x 2 −y 2 -and f -wave superconductivity of square and honeycomb lattice expected at dilute doping are observed by both methods.Upon increasing doping of the square lattice by approximately 0.1, our FRG calculation indicates a p + ip superconducting phase.In the honeycomb lattice model, the f -wave superconductivity persists for higher doping, corroborating the main claim of Refs.[9,11]; but our DMRG data suggest the absence of superconductivity at the Van Hove singularity ν = 1  4 , in The ξ1 (ξ2) is the single-particle (pair) correlation length.For a larger bond dimension χ used in the iDMRG algorithm, a tighter lower bound of ξ is obtained.Here, a smaller ξ2(χ) for a larger nanotube is an artifact of the underestimation becomes more severe for larger systems for fixed χ.The other parameters are 1  8 doping and D = V = 5.
contrast to Refs.[9,11].Recall that near the Van Hove singularity the two-valley picture breaks down.
In the weak coupling regime (V /D 0.3), we perform FRG calculations [23,29] at the one-loop level.We only include the static self-energy and the static two-particle interaction.More particle processes are only included as virtual processes in the two-particle vertex.The inclusion of the static self-energy has been shown to cover already the relevant physics in one-dimensional (1D) systems [30] and can be argued to cover the relevant physics more generally by power counting arguments [23].The static selfenergy incorporates possible further increases of the band gap and deformations of the Fermi surface.For our simulations, we use the unified truncated unity (dubbed TU 2 ) approach [31] merging real and momentum space, which has been demonstrated to fulfill in the FRG equivalence class [32].We distinguish different phases in our FRG simulations by inspecting the eigenvectors corresponding to the largest eigenvalue of each diagrammatic channel.Each of these channels corresponds to a different type of instability and the symmetry of the eigenvector gives the symmetry of the ordering.The Fourier transformation of these eigenvectors at the B sublattice are visualized as insets in Fig. 1.In the strong-coupling regime (V /D 0.5), we use DMRG [21,22] to obtain ground states on infinite cylinder geometries [33].We consider cylinders with up to eight sites along the circumference.The counterpart of 2D superconductivity on the cylinder cannot retain long-range order because of the Mermin-Wagner theorem.However, the pair correlation is expected to be dominant over single-particle correlation.In most common cases, the single-particle excitation of a quasi-1D system is fully gapped (see, e.g., Ref. [34]); the single particle correlation length ξ 1 is finite, while the pair correlation length ξ 2 can diverge.Thus, observation of estimated ξ 2 ξ 1 serves as evidence for such pairing.
The DMRG estimation [35][36][37] usually sets lower bounds for correlation lengths, which become tighter with an increasing number of variational parameters characterized by bond dimension χ [37].
For the square lattice, we find the predicted d x 2 −y 2wave superconductor at low doping within our FRG simulations with t = λ 0 .However, the critical energy scale drops rapidly upon increasing doping and at higher doping (ν ≈ 0.1) we observe a transition to a p x + ip y topological superconductor with Chern number [38] C = 2.The low-doping phase is expected from the above bare interaction analysis.Besides the bare interaction terms revealed above, we speculate a density-fluctuation mediated mechanism [39] can be crucial for stabilizing the higher-doping SC.One observation is a clear increase of SC energy scale at higher doping closer to the density state maxima [Fig.3(b)].The transition between the phases seems to be driven by a change of weight within the particle-particle loop, whereupon doping the d x 2 −y 2 eigenvector will be increasingly suppressed while the p x /p y pair will increase in strength.At stronger interactions (V /D ≈ 0.25) our FRG breaks down, manifested as a linear ramp-up of the density-density interaction.This ramp-up marks the breakdown of the perturbative regime and hence FRG cannot be used to examine the phases.We additionally study the case t = 0 for which the two valleys are absent, the ramp-up problem exists at low doping even for weak interaction, and a C = 1 p + ip is observed for some higher doping [27].From the DMRG data (V /D 0.5 and t = λ 0 ), a finite single-particle correlation length ξ 1 is only consistently found in low doping approximately 1  64 and intermediate interaction; in this case, the pair correlation shows a dominant oscillatory part, supporting the staggered d x 2 −y 2 pairing [Fig.3(d)].We study the geometry with the axial direction along the shortest lattice unit vector (e.g., that connecting nodes 1 and 2 in Fig. 1(a1) right).We consider two-point correlations between ∆ x (i x , i y ) = c ix,iy c ix+1,iy , ∆ y (i x , i y ) = c ix,iy c ix,iy+1 and their Hermitian conjugation.Only the sites on B sublattice are considered.We observe that the signs of ∆ x (0, 0)∆ † x (l, 0) and ∆ y (0, 0)∆ † y (l, 0) oscillate in l; we also observe that the sign of ∆ x (0, 0)∆ † y (l, 0) is opposite to the previous two for a given l.For higher doping, no evidence of convergent ξ 1 is found and no evidence for time-reversal symmetry breaking is found for the implemented larger bond dimensions.While these can be features of a quasi-1-D analog of a Fermi liquid (FL), topological p + ip pairing cannot be excluded.The particle-number-conserved 1D analog of the topological p-wave state has been suggested to be adiabatically connected to an FL [40]; a deeper understanding of the quasi-1-D analog of p+ip is needed to better interpret the data for the p + ip SC or FL region of Fig. 1(b1).The region for large V /D denoted by DW in Fig. 1, is characterized by inhomogeneous densities within the implemented bond dimensions.The 2D phases are speculated to be charge density waves at sufficient commensurate doping; other doping could be Fermi liquids or phases separated by Maxwell construction.The density-wave patterns are difficult to determine as they may only fit on larger cylinders than those studied.
For the honeycomb lattice, we observe f -wave superconductivity in FRG for a broad range of doping, which exceeds Van Hove doping 1  4 .The range is slightly smaller than the random-phase approximation result Ref. [11].Similar to the square lattice, there is also a ramp-up refraining FRG prediction at stronger coupling.Our DMRG for stronger coupling shows a broad range for pairing with a single-particle gap.This is observed for all geometries we studied (e.g., Fig. 3(e)), including the zigzag and armchair nanotube geometries, denoted by (n, 0) and (n, n) [standard notation [33]] respectively where n characterizes the circumference.However, right at the Van Hove doping ν = 1  4 , most cylinder setups including the largest, point to insulating states [27].This feature indicates a possible mechanism of densityfluctuation-induced SC which can accompany a densitywave phase at a commensurate filling [18].This independent mechanism for f -wave SC provides an explanation why the phase extends to higher doping compared to the previous estimation of bare interaction [11].However, the SC energy scale is not largely enhanced closer to the density of states maxima [Fig.3(b)], in contrast to the square lattice.This point is further supported by the high doping state still being an f -wave superconductor [18], such that no competition between the mechanisms is realized.The CL indicates the collapse of fermions leaving part of the system with vanishing occupancy on the B sublattice; collapses are usually observed for models with strong attractive interactions [19,41,42].The observation is that the fermions on lattice B always concentrate on part of the unit cell [see Fig. 3(c)] when increasing iMPS unit cell size.
Discussion.-We examined fermion pairing driven by repulsive interaction and a strong sublattice potential for square lattices and honeycomb lattices.The honeycomb lattice is confirmed to show f -wave pairing, which can be interpreted as inter-valley pairing.The square lattice's counterpart of inter-valley pairing is found to give a low-doping d-wave superconductivity with (π, π) total momentum.Upon increasing doping, a p + ip topological superconductivity is found.Because of the role and existence condition of valleys, the square lattice model with next-nearest-neighbor hopping can exhibit an asymmetry for electron and hole doping.As an outlook, one may also include spin degrees of freedom [10,43] and more types of interactions and hoppings, which serve as extensions of ionic Hubbard models [44][45][46][47][48].This may have implications for real materials and provide the possibility of the sought-after p+ip superconductivity with topological order.
fully acknowledge the computing time granted by the Max Planck Computing and Data Facility.The authors gratefully acknowledge the computing time granted by the JARA Vergabegremium and provided on the JARA Partition part of the supercomputer JURECA [51]  To obtain S (1) , it is more convenient to decompose H k into different ladder operators of the onsite potential and the nearest-neighbor repulsion.We do this first for the hopping process from an A site to a nearby B site.It increases the onsite potential by D. The nearest-neighbor repulsion brought by this process depends on the number of neighbors of A and B, δE = (n B − n A )V .With these observations, the kinetic Hamiltonian is decomposed as where the next-nearest-neighbor hopping includes terms preserving H U and those connecting states with different distributions of A particles.The operators P A/B m,j are projection into the state where the particle A/B at unit cell j has m neighbors occupied.They can be written in terms of the sum over the product operators r n r r (1 − n r ), where r, r are the occupied/unoccupied neighbors.Notice that m, n cannot be larger than N − 1, where N is the number of nearest neighbors in this lattice, as the hopping operators in the middle of Eq. (S3) always eliminate one neighbor.One can verify the following commutation relations With these relations, we can define the leading order transformation for the AB hopping to be iS Such similar expressions S AA , S BB can also be obtained for the AA, BB hopping with different numbers of neighbors occupied.So S (1) is comprised of operators that create excitations of H U .The observation is that those S (1)

AA , S
(1) BB terms annihilate the state with all A sites occupied.So their contribution to the second order expansion [iS (1) , Hk ]/2 vanishes after the projection.What are left are terms diagonal in H U .The resulting expression [iS (1) , H (0) 00 ] t t/U is off-diagonal in the ground state of H U .It will be further eliminated to the order 1/U 2 by the second-order SW transformation.So we can neglect it at the order of 1/U .The Hamiltonian is simplified to a quadratic form of H ± m,n .As the low energy physics is obtained by projecting H to the state with all A sites occupied, this requires the total excitation should have equal numbers of + and − and the sum of m − n should vanish.As H − m,n annihilates the low-energy manifold, it ends up with the following equation As the occupation on site A must be conserved, there are two situations in the above summation.When the bond operators in H − m ,m +N −n and H + N,n are taking the same one, we obtain density interaction terms.When they differ, we have hopping terms for B fermions.It is more convenient to write out the Hamiltonian for the four neighbors around one A site.Choosing (i, j, k, l) to be the four neighbors of a particle A, we have the following terms for the density interaction part ) where PM means distinct combinations obtained by permuting i, j, k, l.For the hopping processes from i → j, we need that the A particle has neighbor i occupied before the hopping and j occupied after the hopping.These terms are given by Similarly, PM means distinct combinations obtained by permuting i, j, k, l Now we collect the contributions together.We have two-body, three-body and four-body interactions: The summation of i, j, k, l is defined by counting the different two-, three-and four-combinations of the neighbors around every A site.So the two-body interaction along the diagonal of the square lattice should be counted twice.Their coefficients are given by The effective hopping Hamiltonian can be assisted by the other two neighbors around each A site where the hopping parameters are All the parameters as a function of V /D are plotted in Fig. S1.

B. Continuum theory for low doping
The effective inverse mass tensor at the two valleys are where the subscript + (−) denotes the valley located at (0, π) ((π, 0)).Recall that the two valleys being minimum is given by the condition |t B /t B | > 0.5.We see that the mass tensor is diagonal in the coordinate we choose, and there is mass anisotropy for each valley.We introduce the center-of-mass coordinates: δr = r + − r − , R = m+ r + + m− r − , where m± = m ± / Tr(m ± ).The (first-quantized) kinetic Hamiltonian of the two par-ticles can be written as (∇ T δr µ ) and I is the 2 × 2 identity matrix.The continuum approximation of the two-particle problem is similar to that of the honeycomb model [9].The two-particle binding energies can be calculated in the center-of-mass reference frame [9]: where µ 0 = (8a 2 t B ) −1 is the relative mass.The result is independent of the specific value of t B (t in the full model) as long as the ground states of single fermion [two fermions] are in the (0, π) or (π, 0) [(π, π)] momentum sector, consistent with the microscopic Hamiltonian being independent of t B in the two-particle (π, π) sector.We fit e bp using exact data in the small V /D region and plot the continuum effective result together with the exact result in Fig. S2.Different from the honeycomb model, the binding energies will drop for large V because |g| will decay to zero for large V instead of converging to a constant.
To extract coefficients of interaction in the continuum theory, we write the correlated hopping as where a i is taken from the four vectors connecting A to its nearest B neighbours.The variables k, q, k are taken to be much smaller than |K + − K − |.Similarly, the repulsion term is rewritten as The continuum interactions between different valleys and inside the same valley are given by taking appropriate combinations of σ i and their anti-symmetrized partners.
The result of inter-valley interaction has been given in the main text.Now we consider intra-valley interaction with finite doping.We consider the weak interacting limit and discuss the interaction between modes on the Fermi surface of a valley, for example, the + valley.Here the momenta are defined as the deviation to (0, π).In general, we have.We focus on two-particle scattering with net zero deviation to K + and low doping (momenta is small enough to perform Taylor expansion). q1,q2 Set q 1 = q 2 and q 1 = −q 2 respectively, we can obtain the density-density interaction, ∝ [2λ 1 (q 2 1,x − q 2 1,y ) + U 2 (q 2  1,x + q 2 1,y )]n + (q 1 )n + (−q 1 ).With the Fermi surface shape close to an eclipse with the long axis along the y direction, and λ 1 ≈ U 2 for V /D < 1 (Fig. S1), such interaction in most momenta is attractive.
C. Few-particle-doping binding energy with finite D/t Here, we discuss few-particle-doping binding energies of the full model with finite D/t = 10, 5.The binding energies, in the unit of t 2 /D, are in general smaller for smaller D/t .However, even for D/t = 5, we find no substantial difference for inferred stable pairing region, compared to the effective model with D/t = ∞.
As an alternate of binding energy per particle, we represent the results as binding energies for forming bound states with composites.For two-particle and threeparticle bound states, the existence of bound states can be seen from a positive E b,1,1 = 2E 1 − E 0 − E 2 and E b,1,2 = E 1 +E 2 −E 0 −E 3 respectively.These quantities are plotted in Fig. S3.Binding energy per particle can be deduced from them.The existence of three-particle bound states does not mean that three-particle bound states are more favored than pairs for dilute doping.Favored bound states have the largest binding energy per particle.Some details of the numerical implementation are as follows.For the full (effective) lattice models, we obtain the ground-state energies E n of finite systems using DMRG (exact) diagonalization.To accurately compute binding energies, we need system sizes larger than the sizes of the bound states.We estimate the finitesize errors by doing 1/N extrapolation for data of the two largest systems we obtain.The extrapolated data in Fig. S3 have errors smaller than the size of markers.The periodic boundary condition is implemented in the exact diagonalization, which enables reading the momentum quantum numbers.DMRG is less efficient to deal with periodic conditions along two directions and we thus only implement a periodic boundary condition in one direction while implementing an open boundary condition on the other.In this case, to correctly calculate the bulk binding energy, we find it essential to eliminate the low-energy edge modes.Such modes can be understood by considering the potential and interaction part of the Hamiltonian Eq. 1: i,j V n i n j + i∈B Dn i .As A sublattice is almost fully-filled, if a fermion on B sublattice is located at the open boundary rather than in the bulk, it feels less repulsion from the fermions on B sublattice.Consequently, these configurations have lower energy.We find that introducing additional potential terms V n i on the boundary of B sublattice can eliminate low-energy edge modes.

D. Weak-coupling results
In the weak coupling regime, we apply a truncated unity functional renormalization group approach.The FRG flavor we employ is a method to construct unbi- ased single and two-particle interactions from a set of flow equations.While three-and four-body interactions are not counted as objects themselves, they are partially included via virtual processes in the two-particle interaction.We use a sharp energy cutoff [23], thus the critical energy scale can be interpreted as a critical temperature modulo an unknown scaling factor.We calculate the vertex on a 24×24 momentum mesh for both lattices, with a refinement for the bubble integration mesh of 45×45.On the square lattice, we include the 25/29 nearest neighbors in the truncated unity per site within the unit cell for the honeycomb/square lattice.We use a Bogacki-Shampine adaptive integrator for the integration of the flow equations, allowing for a maximal absolute error of 10 −2 per integration step.The results of the FRG simulations are visualized in Fig. S5.To distinguish different phases, we inspect the behavior of the maximal eigenvalues of each channel during the flow in combination with an inspection of the dominant eigenvectors at the end of the flow.These eigenvectors encode symmetries and specific types of instability.In the case of p y /p x we find the two eigenvectors to be exactly degenerate.In the real space representation we define p x/y = sign( v x/y • d)δ v x/y , d with v x = (1, 0), v y = (0, 1) and d is the vectorial distance between two sites.To distinguish all possible linear combinations cos(θ)p x + e iφ sin(θ)p y we perform a single-step mean-field calculation and compare the free energy of each starting configuration, as can be seen in Fig. S6.
To calculate the Chern number in the gapped phase, we employ the method described in Ref. [38].

E. Details of infinite DMRG calculations
We obtain approximate ground states of the Hamiltonian Eq. 1 defined on infinite cylinders.We do this by optimizing infinite matrix product states via two-site iDMRG algorithm [22].We implement the conservation of particle numbers, thus the phase diagram Fig. 1 is constructed in terms of doping densities ν.The accuracy of infinite matrix product states can be improved by increasing its bond dimensions (χ, size of the matrices); With efficient optimization, to reach a given accuracy, the required computational resource (e.g.χ) is exponentially large in cylinder circumference.Infinite matrix product states are constructed to be exactly translationally invariant by M lattice unit vector along the axial direction.To implement exact particle number conservation for doping density ν = p/q (irreducible fraction), M must be integer multiples of q/L y .(L y is the number of lattice unit vectors around the cylinder.)To be compatible with the possible spontaneous breaking of translational symmetry (e.g., charge-density-wave state), M has to be compatible with the enlarged unit cell.As mentioned in the main text, we estimate of correlation lengths of single-particle and pair to infer superconductivity.See Ref. [35,36] for the definition and extraction methods for one-dimensional correlation lengths.Here, we define the correlation length from the correlation along the axial (infinite) direction.We denote them estimated using bond dimension χ as ξ 1 (χ) and ξ 2 (χ) respectively; the larger the bond dimension, the more accurate the estimation.

FIG. 1 .
FIG. 1. Lattice structure and schematic phase diagrams.(a) Lattice structure for (ai) square and (aii) honeycomb lattices; sublattices are marked in orange and cyan.(b) Phase diagram for (bi) the square lattice with t = λ0 and (bii) the honeycomb lattice with t = 0, inferred from FRG (weak coupling) and DMRG (strong coupling).For superconducting (SC) phases, we plot the momentum dependence of the susceptibility from FRG.For the p + ip SC phase, a degenerate pair of dominant eigenvectors is found in FRG, and a mean-field analysis indicates the linear combination p + ip is favored.Here DW denotes the density-wave phase, CL denotes a phase separation via the collapse of electrons on the B sublattice, and FL denotes Fermi liquids.For the uncertainty in the DMRG data interpretation, see the discussion in the main text.

FIG. 3 .
FIG.3.Results from exact diagonalization (ED) of the effective model for few-particle doping, with FRG and DMRG of the full model for finite doping.(a) Energy per particle of the n-particle ground states; the unit is t 2 /D.Calculations are performed for the effective models using ED with finite-size extrapolation.From left to right are the square lattice model with t = λ0 and the honeycomb lattice t = 0. Finite E1 − E2/2 indicates the existence of two-particle bound states.Some paired phase near the dilute limit is indicated by the energy per particle E1 − En/n being a constant for every positive even n and larger than the values for odd n.On the other hand, if there is some n, with E1 − En/n greater than E1 − E2/2, larger bound states are favored.(b) FRG predicted phases and energy scales λc in units of t 2 /D.Shown on the left is the square lattice model with t = λ0 and on the right is the honeycomb lattice.We choose V = 2t as the interaction.(c) Cylinder geometry [(6, 0) honeycomb nanotube] and a density plot of an inhomogeneous density profile, indicating phase separation.The densities on B sublattice are represented by colors with green for small density and red for larger density.The parameters for the density plot are 1 8 doping, D = 10, V = 15, and unit cell size 32; only part of the unit cell is plotted.(d) Pair correlation functions for square lattice at 1 64 doping, V = 6, D = 10 and a tangential direction size eight unit cells.For the definitions of pair operators see the text; the subscripts m, n can be x, y.(e) Correlation length evidence for pairing on the infinite cylinder geometry of a honeycomb lattice.The indices for nanotubes are standard notations for their sizes and shapes.The ξ1 (ξ2) is the single-particle (pair) correlation length.For a larger bond dimension χ used in the iDMRG algorithm, a tighter lower bound of ξ is obtained.Here, a smaller ξ2(χ) for a larger nanotube is an artifact of the underestimation becomes more severe for larger systems for fixed χ.The other parameters are1  8 doping and D = V = 5.

4 FIG
FIG. S1.Parameters of the large-D effective theory of the square lattice model.

2 E
FIG.S2.Binding energies and bound states sizes of two fermions for the effective square lattice model with two valleys.The dashed line denotes the continuum effective theory with one fitting parameter given by matching small V /D data.
FIG. S3.Binding energies.The binding energies of two fermions are shown in the first row from left to right for the honeycomb lattice model and the square lattice model without t term and for t = λ0.Correspondingly, the binding energies for a fermion pair and a fermion to form a threefermion bound state are shown in the second row.The plot scale of the E b,1,1 of square lattice models is smaller than others'.
FIG. S5. Results of the FRG simulations for the three different setups and visualization of the f -wave superconductivity.We abbreviate a flow to strong coupling without divergent susceptibility as SCT, and a divergence of the f -wave component of the pairing susceptibility as f -SC.The y-axis displays the critical energy scale, which is linearly dependent on the critical temperature.The x-axis displays the doping.(a) shows the results for a square lattice with t = λ0 and (b) the square lattice with t = 0 results.(c) shows the results for the honeycomb lattice.In all simulations, we kept D = 10 and varied V in the given range.(d) visualizes the different superconducting order parameter symmetries encountered, as visualized by the eigenvector of the effective two-particle interaction at the orbital at the Fermi level.The upper left shows an f -wave on the honeycomb lattice, the upper right shows a d x 2 −y 2 -wave on the square lattice.The lower two plots correspond to the degenerate pair px and py with weak admixture of other dependencies.
FIG.S6.Free energy minimization in the square lattice with t = λ0, using cos(θ)px + e iφ sin(θ)py as starting values.The combination with the smallest free energy is1  2 (px + ipy) at Forschungszentrum Jülich.Y.H., J.B.H., and D.M.K. were supported by the Deutsche Forschungsgemeinschaft (German Research Foundation) under Grant No. RTG 1995, within the Priority Program SPP 2244 "2DMP" and under Germany's Excellence Strategy-Cluster of Excellence Matter and Light for Quantum Computing (ML4Q) Grant No. EXC 2004/1 -390534769.K.Y. and E.J.B. were supported by the Swedish Research Council (VR, 2018-00313), and the Wallenberg Academy Fellows program of the Knut and Alice Wallenberg Foundation (Grant No. 2018.0460).