Spin and Charge Fluctuation Induced Pairing in ABCB Tetralayer Graphene

Motivated by the recent experimental realization of ABCB stacked tetralayer graphene [Wirth et al., ACS Nano 16, 16617 (2022)], we study correlated phenomena in moir\'e-less graphene tetralayers for realistic interaction profiles using an orbital resolved random phase approximation approach. We demonstrate that magnetic fluctuations originating from local interactions are crucial close to the van Hove singularities on the electron- and hole-doped side promoting layer selective ferrimagnetic states. Spin fluctuations around these magnetic states enhance unconventional spin-triplet, valley-singlet superconductivity with $f$-wave symmetry due to intervalley scattering. Charge fluctuations arising from long range Coulomb interactions promote doubly degenerate p-wave superconductivity close to the van Hove singularities. At the conduction band edge of ABCB graphene, we find that both spin and charge fluctuations drive $f$-wave superconductivity. Our analysis suggests a strong competition between superconducting states emerging from long- and short-ranged Coulomb interactions and thus stresses the importance of microscopically derived interaction profiles to make reliable predictions for the origin of superconductivity in graphene based heterostructures.

Introduction.The experimental discovery of cascades of correlated phases and superconductivity in ultra-clean graphene multilayers without [1][2][3][4][5][6] and with [7][8][9][10][11][12][13][14][15][16][17][18][19] a stacking twist has led to tremendous research interest.Experimental studies of multilayer graphene suggest the existence of ferromagnetic phases in the form of half-and quarter metals [1,20], Wigner crystals [1] and in particular superconductivity that emerges in Bernal bilayer (AB) and rhombohedral trilayer (ABC) graphene for external displacement fields and either in the presence of external magnetic fields [2,4] or proximity induced spin-orbit coupling [3,5].From a fabrication point of view, untwisted graphene stacks are simpler to handle as twist angle variations and other stacking imperfections are easier to control.It is therefore that experiments along with atomistic theoretical studies can provide a microscopic understanding of correlation effects, in particular as hindering disorder effects can be disregarded and multilayer graphene stacks omit nm-scale unit cells as well as a more sophisticated theory of topological obstruction of the low-energy bands when compared to their twisted counterparts [21][22][23].
Nonetheless band flattening in untwisted graphene stacks is restricted to small Fermi surface patches around the valleys  and  ′ , which hampers theoretical descriptions without losing contact to the microscopics.Previous works [24][25][26][27][28][29][30][31][32][33][34][35][36][37] resort to effective continuum or adapted Slonzecewski-Weiss-McClure models with ultra-violet cutoff that are only valid near the valleys  (′) and hence severely complicate a description of realistic orbital-resolved interaction profiles that account for short and long-ranged interactions consistently.Instead, an approximate spin/valley  (4) symmetry arises from considering the long-ranged tail of the Coulomb interaction alone.Motivated by the recent experimental characterization of different tetralayer graphene stacks by Wirth et al. [38], we remedy the aforementioned shortcomings and study correlated phenomena in moiré-less graphene tetralayers for ab-initio motivated models and realistic interaction profiles using an orbital resolved random phase approximation approach (RPA).
In particular, we focus on ABCB stacked graphene that breaks inversion symmetry and thus features intrinsic electrical fields.The latter open a band gap at charge neutrality and hence play the same role as external displacement fields that are necessary to stabilize superconductivity in bi-and trilayer graphene.While previous works mainly focused on pairing mediated by electron-phonon coupling [39][40][41], fluctuations around flavor symmetry broken phases [24][25][26][27] or charge fluctuations arising from repulsive Coulomb interactions [28][29][30][31][32][33][34][35][36][37], we discuss the interplay of spin and charge fluctuations to the formation of weak coupling instabilities on general grounds based on a miscroscopic model of the carbon   -orbitals.In view of the various spin-polarized phases governing large areas of the phase diagram of bi-and trilayer graphene [1], we first characterize the role of spin-fluctuations towards the formation of magnetic order in ABCB graphene.To disentangle the influence of long-and short-ranged interactions as well as their contributions to spin-and charge enhanced superconductivity, we next study three different pairing mechanisms based on (i) spin-fluctuations from local interactions (RPA) (ii) screening of long-ranged Coulomb interactions (PRA) and (iii) a combined approach (RPA) that captures both long and short ranged Coulomb interactions.Microscopic Model.We model the electronic band structure of graphene tetralayers employing a modified Slonzecewski-Weiss-McClure Hamiltonian [42].The resulting kinetic Hamiltonian represents an atomistic model of tetralayer graphene consistent with first principle calculations [43], see Supplementary Material (SM) [44] for details on the parametrization.The upper panels of Fig. 1 give an overview of the non-interacting electronic structure of the three inequivalent thermally stable stacking configurations of tetralayer graphene: ABCB (a), ABCA (b), and ABAB (c).We plot the low-energy band structure along the two irreducible directions in vicinity of the valley .The color encodes the polarization   of the dominant orbitals.For ABCB graphene we find that most weight is concentrated on the a-site in the first layer and the b-site in the third layer.Due to broken layer inversion symmetry, ABCB graphene features a band gap of Δ  ≈ 20 meV at charge neutrality without the application of an external electric field.This is reflected in the density of states [DOS; see Fig. 1 (d)] that indicates van Hove singu- larities (VHS) close to charge neutrality on the electron-and hole-doped side.As the internal electric field in ABCB causes valley band flattening, the VHS is stronger in ABCB compared to ABCA and ABAB graphene.Therefore we expect the effects of electron-electron interactions to be most significant in ABCB graphene.By inspection of the ABCB Fermi surface for various fillings we observe a Lifshitz transition from a (annular) single-pocket to a three-pocket structure, which signals an intriguing interplay between nesting-and DOS-driven Fermi surface instabilities and resembles AB graphene under an electrical field.As the latter hosts magnetic states [1], we start with a discussion of spin fluctuation induced instabilities in the following.

Spin fluctuations & magnetic instabilities.
To study magnetic correlations in few-layer graphene we use a random phase approximation (RPA) approach.We calculate the free electronic susceptibility χ0 () = Tr  [ Ĝ0 ( − ) ⊙ Ĝ 0 ()] in the static limit  0 → 0, where Ĝ0 () is the free Matsubara Green's function as a matrix in orbital space,  = ( 0 , ) the electronic "four-momentum", and "⊙" an element-wise matrix product.Assuming a local Hubbard interaction  allows us to resum the infinite ladder series in the crossed particlehole channel to arrive at a Stoner-like criterion [45] for the orbital-resolved susceptibility χ0 (): The critical interaction strength required for the onset of magnetic order is given by  c = 1∕  0 (  ), with   0 (  ) the largest eigenvalue of χ0 (  ) and   the transfer momentum at which the maximum occurs.We adapt  = 5 ⋅ 10 −5 eV as temperature broadening for χ0 () throughout this work, see SM [44] for more details on the orbital-space RPA.
Figure 2 (a) demonstrates that magnetism is potentially realized in ABCB tetralayer graphene near the valence band VHS.The critical interaction strength drops to  c ≈ 4.5 eV, which is of similar order as estimates of the effective (screened) onsite  in graphene [46,47].In Fig. 2 (b) we observe that the valence band VHS drives ferrimagnetic fluctuations within the C-layer of ABCB graphene, while away from the VHS layer-agnostic antiferromagnetic fluctuations gain in importance.The momentum structure of the magnetic susceptibility is shown in panel (c), where we plot the leading susceptibility eigenvalue for all  close to Γ and around the  (′) pointsthese are the areas in momentum space where we expect   0 () to pick up structure by either the VHSs or nesting.While the detailed weight distribution in   0 () evolves upon doping across the valence band VHS, we observe that the leading contribution always stems from momenta close to Γ. Therefore intravalley scatterings and the divergent DOS mainly drive magnetic order.Other orders as, e.g., intervalley coherent order, may emerge as subsequent instabilities from the magnetic one in the presence of purely local interactions.For the VHS on the electron-doped side, a similar picture of magnetic instabilities emerges (see SM [44]), allowing us to focus on the hole-doped side.
Pairing from local interactions.As established in the previous paragraph, local interactions cause ferrimagnetic correlations to prevail in ABCB graphene on the hole-doped side.This is particularly striking as large areas in the phase diagram of bilayer graphene are governed by symmetry-broken phases that show clear indications of hysteresis [1].Close to the ferrimagnetic instability, spin fluctuations arising from local, repulsive interactions can provide the pairing glue for an unconventional superconducting state [48][49][50][51][52][53].The effective pairing vertex in static approximation for this Kohn-Luttinger like mechanism includes contributions from local interactions to the exchange-and direct particle-hole channels [49] where  is the local onsite interaction,  () denotes the dominant momentum transfer in the exchange (direct) particle-hole channel and the hat symbol denotes the orbital structure of interactions and polarization function, see SM [44] for details.Within the spin fluctuation exchange formalism (xRPA), we finally solve the linearized gap equation in the projected subspace containing the bands in an energy window  around the Fermi level where V PP , ′ and  PP  ′ denote the band-projected effective pairing vertex and particle-particle susceptibility, respectively.The largest coupling constant therefore gives an estimate for the critical temperature   =   −1∕  of the superconducting (SC) transition [54] and the corresponding eigenfunction yields the symmetry of the SC order parameter Δ  ().
Since local interactions are flat in momentum space, they induce both inter-and intravalley coupling.In conjunction with single layer ferromagnetic fluctuations, the intervalley exchange enhances order parameters that are (i) spin-triplet and (ii) change sign under a valley flip.We find that local interactions exclusively promote valley-singlet, spin-triplet wave order for all fillings around the VHS in ABCB stacked graphene.The coupling constant   is enhanced close to the Stoner transition for  = 4 eV as plotted in purple in Fig. 3 (a), while it decreases continuously when doping away from the VHS.As soon as the local interaction strength  exceeds the critical value of the Stoner transition   , superconductivity is suppressed in the Stoner regime and the peak of the coupling constant   splits into two as shown in Fig. 3 (b).We note that further screening of local interactions may therefore shift the position of superconducting regions observed in experiment [11].
Pairing from screened Coulomb interactions.Screened Coulomb repulsion was discussed to provide the dominant pairing glue in bi-and trilayer graphene [28, 30-32, 35, 55].In the continuum theories employed by a majority of works, the additional flavor degeneracy due to the valley degree of freedom enhances the polarization function  0 () by a factor of 4 instead of 2 compared to the usual spin  (2), which boosts the contribution of charge-fluctuations to electron mediated pairing.However, effective continuum theories hinder a direct description of the microscopic interactions containing long-and short-ranged terms.This manifests in an artificial spin/valley  (4) symmetry that must be lifted by an phenomenological inter-valley Hunds coupling [35,55].Here, we remedy this shortcoming by exploiting the orbitalresolved RPA approach presented in this manuscript and use the "Ohno" interaction profile [46,56] with realistic parameters  = 0.3 0 ,  = 200 0 and  0 = 2.46 Å the graphene lattice constant.The parameter  controls the external screening from, e.g., metallic gates or the substrate,  controls the atomic scale decay of the Coulomb interaction in real space, and  sets the on-site (Hubbard) interaction strength.Pairing due to interaction-induced screening is then captured by the effective vertex where V  () is the discrete Fourier transform of the Ohno interaction matrix obtained from Eq. ( 3) and   is the relevant momentum transfer in the direct particle-hole channel.Taking only screened Coulomb interactions into account (dRPA), the superconducting coupling constant shows two peaks that are located near the VHS and the edge of the valence bands, see Fig. 3 (a,b).Near the VHS, the superconducting order parameter with highest transition temperature is a doubly degenerate spintriplet with  , symmetry that changes sign within each valley as shown in panel (c).
We argue that order parameters with an intravalley sign change are generically favored by the screened Coulomb interaction for realistic interaction parameters: Microscopically, screening induced pairing Eq. ( 4) is most relevant for momenta   where the polarization function χ0 (  ) peaks. Figure 2 (c) demonstrates that these momenta are   ≈ Γ. Therefore the repulsive Coulomb interaction V  (  ) gets suppressed significantly close to Γ. Nevertheless the screened interaction V PP  (,  ′ ) is still mostly positive in momentum space, c.f. Eq. ( 4) and SM [44], such that possible order parameters must change sign on the respective Fermi surface sheets connected by the vector   .Scattering of Cooper pairs is therefore selectively enhanced within the same valley ( ≈  ′ ).On the other hand, intervalley scattering   ≈  (′) caused by shortranged terms in the Coulomb interaction is suppressed.We note that this competition of inter-vs.intravalley scattering can only be resolved when considering the full Coulomb interaction with respective ratios of on-site and long range interaction strengths.
Indeed the two spin-triplet  , -wave order parameters are accompanied by two spin-singlet   2 − 2 , order parameters.The latter come with a slightly lower coupling constant   at the VHS, see Fig. 4. Figure 3 (a,b) further demonstrates that within the dRPA approach,  -wave and  + -wave symmetric gap functions are suppressed around the VHS for realistic values of on-site  such that the  , -wave pairing prevails.
Varying  has no influence on the dRPA results, see SM [44].
The situation at the valence band edge is fundamentally different: Here, the increase of the superconducting coupling constant   can be traced back to the trigonal warping related three pocket structure of the Fermi surface in ABCB.As the Fermi surface shrinks,  , -wave order is suppressed because an intravalley sign change is increasingly disfavored.Instead almost degenerate  -and  + -wave order prevails, which is driven by local interactions as in the xRPA (spin) case.The critical temperature reaches values up to   ≈ 10 mK and the predominant  ∕ + -wave order stays unmodified for different (numerical) momentum resolutions, see SM [44].
Competition of spin & charge fluctuations.Finally, we study the mutual interference of the charge and spin mechanisms discussed above by carrying out independent resummations of (i) short range interactions in the exchange channel and (ii) long range interactions in the direct channel (xdRPA), see SM [44].We observe that values of the coupling constant   lie between those of the spin fluctuation (xRPA) and the screening mechanism (dRPA) with critical temperatures reaching   ≈ 25 mK, whereas the order parameter symmetry shows  , -wave order near the VHS and  -wave order at the valence band edge.Even though the overall density dependence and symmetry of the superconducting state resembles pure screening driven SC, important insights can be gained from this analysis.First, we note that the coupling constant   is enhanced compared to the pure screening mechanism as shown in Fig. 3. Therefore, we conclude that spin-fluctuations support charge-fluctuations in the formation of superconductivity in ABCB, especially close to the Stoner transition.Second, the approximate degeneracies between and  + -wave order close to the valence band edge are significantly lifted in the combined xdRPA approach, as well as between  , -and   2 − 2 , -wave order around the VHS. Figure 4 presents a refined analysis of the superconducting tendencies for  =  VHS for different values of the on-site interaction strength  .Since fluctuations around the ferrimagnetic phase suppress spin-singlet superconductivity, the nearly degenerate  ∕ + -and  , ∕  2 − 2 , -wave order parameters acquire significant splitting with increasing  and only the triplet order parameters ( -and  , -wave) prevail.Third, the pairing glue arising from local interactions is significantly reduced by the long range tail of the Coulomb repulsion such that spin fluctuation induced pairing becomes effectively suppressed far from the Stoner transition.This is underlined by the observation that for larger on-site interactions [ = 5 eV, cf.Fig. 3 (b)], spin fluctuation driven SC (purple line) is visibly enhanced due to the vicinity to the magnetic state, while the effective coupling constant   within the combined xdRPA approach (orange line) is merely affected.At the VHS the situation is more delicate: As  →  c , we see the mutual enhancement of  , -and  -wave order parameters which benefit from the spin-fluctuation enhancement, see Fig. 4. In the immediate vicinity to the magnetic instability the xRPA behavior is recovered such that for  ≲  c the leading instability is again of  -wave character.While our combined xdRPA approach performs resummations in both particle-hole interaction channels (crossed & direct), it cannot capture inter-channel feedback, which could account for screened local interactions at the VHS, i.e., suppress the magnetic instability, and therefore boost  -wave SC at the VHS.
Discussion.In this work we discuss the role of spin fluctuations and screening on the formation of correlated magnetic and superconducting states in tetralayer graphene with emphasis on ABCB stacking.The particular shape and orbital polarization of the bands in ABCB on the hole-doped side resembles the bandstructure of bilayer graphene in an ex-ternal electric field.Since ABCB graphene is not topologically obstructed as relative twisted graphene stacks, we expect our result to be qualitatively transferable to other stacking configurations of few-layer graphene, such as AB bilayer graphene.In contrast to ABCB, superconductivity in AB graphene emerges at relatively high (experimentally accessible) displacement fields [2][3][4][5], implying that ABCB graphene provides an interesting avenue to further stabilize superconductivity in graphene-based multilayers.Due to the metastable nature of ABCB graphene [38], the experimental realization of transport devices is likely to pose special challenges for the sample preparation, especially for the encapsulation of ABCB graphene in hexagonal boron nitride like it is known from ABC trilayer graphene [20,57].
Our results further raise important consequences originating from the orbital-resolved weak coupling treatment that captures long and short-ranged Coulomb interactions consistently.First, we argue that for realistic values of the on-site interaction [46,47]  , -wave SC order will always prevail within a pure screening mechanism.This is in contrast to previous works [28,32] that show emerging  ∕ + -wave superconductivity at the VHS in bi-and tetralayer graphene stacks resulting from screened Coulomb interactions alone.Second, we observe another peak in the superconducting coupling constant   for low hole-doping with emergent  ∕ + -wave symmetry.Previous works by Ghazaryan et.al [35,55] report a slim peak near the valence band edge with dominant -wave symmetry, which contradicts our results and can be traced back to the effect of local interactions that are captured within the orbital-resolved RPA.Our xdRPA approach further highlights the importance of studying electronic phases in multilayer graphene stacks with methods that (i) allow for inter-channel feedback, e.g. the functional renormalization group [58,59], and (ii) incorporate local Coulomb repulsion on the same footing as the long range tail.

S1. NON-INTERACTING HAMILTONIAN
Our tight binding Hamiltonian for graphene tetralayers is a modified Slonzecewski-Weiss-McClure (SWMC) Hamiltonian, with parameters chosen based on Ref. [42] summarized in Table SI.The Hamilton matrix can be constructed as where  and  are site indices within the unit cell,  marks the type of the lattice and is either ABAB, ABCA or ABCB and  , is one if  is an a-site and 0 otherwise. 1 and  2 iterate over all unit cells in the infinite lattice, ( 1 ,  2 ) gives the vectorial distance between the unit cells and gives the distance between site  and the image of site  shifted by  1 and  2 unit-cell vectors.

S2. RANDOM PHASE APPROXIMATION
To calculate the non-interacting susceptibility χ0 () in static approximation, we must carry out the summation where   () denotes the dispersion relation and   () the Bloch function at momentum  for band  and orbital (site) .In the actual numerical implementation, we make use of the fact that the only momenta  relevant for magnetic instabilities are those where  ≈   ±  ′  , with  (′)  on the Fermi surface.We thus choose an energy widow  ≫  (in our calculation  = 3 meV and  = 0.025 … 0.1 meV) and obtain all momenta  ∈ { ±  ′ with ,  ′ ∈    } for the extended Fermi surface    = { ∈ ℤ s.t.∃ with |  () − | <  }.The calculation of χ0 () is perfectly parallel in .Moreover, we parallelize the summation over  to obtain higher performance in our custom CUDA kernels.We additionally note that we make use of a cached dispersion relation and obrital to band transform on a fixed momentum mesh (we checked 4800 × 4800 and 7200 × 7200 points).
From the free electronic susceptibility χ0 () on the relevant , we calculate the critical on-site interaction strength required for the onset of magnetic order.The crossed particle-hole channel RPA resummation of the effective interaction yields where we directly observe that for an on-site Hubbard interaction Û =  , the maximum eigenvalue of χ0 () over all  determines (i) the critical interaction strength  c = 1∕  0 (  ) and (ii) the type of magnetic order via the corresponding eigenvector ⃗   0 (  ).

A. Flowing random phase approximation
An alternative way to look at the RPA from above is via a single-channel functional renormalization group (FRG) flow.From the diagrams contained in the two particle irreducible flow with neglected self-energies [61] it is evident that the effective FRG vertex in static approximation sums up all diagrams contained in the effective RPA vertex Eq. (S3) when restricted to the crossed particle hole channel.We therefore rewrite the RPA as a single channel FRG, i.e., ̇Ŵ Λ () = − ŴΛ () ̇L Λ () ŴΛ () . (S4) One quickly verifies that the solution to this differential equation is given by From the FRG perspective, it is therefore sufficient to calculate the particle-hole loop integral to get an expression for the free susceptibility at scale  : Note that we employ a sharp frequency cutoff here.χFRG  () resembles an approximation to the Matsubara summation in the explicit RPA, where the frequencies above the lowest one are treated as continuum.Such an approximation is similar to the explicit Matsubara summation used in large unit cell systems in Ref. [45].We observe that in the flowing RPA formulation, though, we can make use of fast Fourier transforms (FFTs) to calculate the convolution over  and .For large momentum meshes, the numerical scaling of this operation is (  log   ) instead of (    ) for the explicit sum over  for each .Using the flowing RPA trick with FFTs allows us to bump the momentum resolutions significantly higher than using the explicit RPA: We calculated the flowing RPA susceptibilities for a momentum resolution of   = 12288 2 = (2 12 ⋅ 3) 2 points and confirmed that indeed our calculations are converged.Note that numerically, the flowing RPA covering roughly 400 adaptively chosen points in the Λ integral were significantly faster on a single CPU node of the JURECA cluster [60] than the explicit RPA CUDA kernels on four NVIDIA A100 GPUs with approximately half the resolution (  = 7200 2 ).

B. Magnetic instabilities for other parameters
To confirm that the conduction band physics is qualitatively equivalent the effects in the valence band, we show the magnetic instabilities for both sides at  = 10 −4 eV in Fig. S1.
The behavior of χ0 () as a function of temperature is depicted in Fig. S2.Note that we include high resolution, low temperature calculations that were obtained with the flowing RPA trick presented in Section S2 A in panels (d,e).

S3. SUPERCONDUCTIVITY
To study unconventional superconductivity driven by electronic interactions, we study different mechanisms that are based on (i) spin-fluctuations by local interactions (RPA) (ii) charge-fluctuations due to screening of the long-ranged Coulomb interactions (RPA) and (iii) a combined approach (RPA) that captures a superset of diagrams contained in the mechanisms above.This section is denoted to show the corresponding diagrams describing the contributions of either mechanism to superconducting pairing and discuss their practical implementation.

A. Local Interactions
In the vicinity of magnetic instabilities, spin and charge fluctuations that are driven by local Hubbard- interactions can provide the pairing glue for an unconventional superconducting state.To capture the effect of spin-fluctuation induced pairing in ABCB, we resort to the well-known diagrammatic expansion [49,62] based on transversal and longitudinal spin-fluctuations as shown in The effective spin fluctuation vertex Eq. (S7) is proportional to (1−  0 ) −1 reminiscent of the magnetic instability that is reached when the Stoner condition at  →   is met.For  ≥   the theory hence breaks down as magnetic fluctuations diverge and the system orders magnetically.To avoid numerical instabilities when calculating superconducting gap functions using the effective pairing vertex  PP  , we perform the matrix inversions 1∕(1 −   0 ) and 1∕(1 −  2 ( 0 ) 2 ) in Eq. (S7) in the eigenspace of  0 and add a small imaginary broadening constant   FLEX = 50 meV to all the eigenvalues in the denominator.Note that the energy scale of this broadening constant is to be compared with (and has to be very small in regard to) the vertex energy scales, i.e.

− 5 eV and not the valley-flat band energy scales.
As explained in the main text, local interactions are flat in momentum space such that both inter-and intravalley coupling are present.In conjunction with single layer ferromagnetic fluctuations, the intervalley exchange enhances order parameters that are (i) spin-triplet and (ii) change sign under a valley flip.We find that local interactions exclusively promote valley-singlet, spintriplet  -wave order for all fillings around the VHS in ABCB stacked graphene.We can understand the prevailance of  -wave superconductivity from a miscroscopic picture when analyzing the effective vertex structure Eq. (S7) as shown in Figure S3.
As known from the one-band Hubbard model [62] the spin fluctuation vertex remains purely positive in momentum space such that the symmetry of the superconducting order parameter must exhibit a sign change at the respective Fermi surface.We note that even though the effective vertex peaks around Γ, see Figure S3 (a,b and subsequent columns), the effective vertex has substantial amplitude for momentum transfers that connect the two valleys  (′) .Therefore, the superconducting condensate can minimize its free energy by changing sign between the valleys leading to dominante  -wave SC order.At the same time, the sub-leading instability within the RPA is of  , -wave type, see Figure 4 (b) in the main text.The  , spin-triplet instability is an anti-symmetric basis function of the  irreducible representation of  3 in momentum space.Similar to the RPA results, the  , -wave is driven by fluctuations that are centered at Γ.We hence note that the delicate interplay of short and long-ranged interactions, which we resolve within our orbital-resolved RPA, drive either order parameter.

B. Long-Ranged Interactions
Pairing mediated by screened Coulomb interactions can be accounted for by the effective pairing vertex where V  () =    1  2 () is the Fourier transform of the full Coulomb interaction and the factor of 2 in the denominator arises due to the internal trace over spin indices in the direct particle-hole channel.The latter is a direct consequence of the  (2) symmetry in the normal-state Hamiltonian.As described in the main text, we use a realistic (orbital-resolved) "Ohno" interaction profile [46,56] with realistic parameters  = 0.3 0 and  = 200 0 (where  0 = 2.46 Å is the graphene lattice constant).The parameter  controls the external screening from, e.g., metallic gates or the substrate,  controls the atomic scale decay of the Coulomb interaction in real space, and  sets the on-site (Hubbard) interaction strength.In particular, the 'Ohno" interaction profile ensures consistent treatment of long-and short-ranged terms of the Coulomb interaction.We further checked that our results do not change when varying the screening parameter in the realistic range of  = 50 … 200 0 , see Section S3 E. As argued in the main text, order parameters with an intravalley sign change are generically favored by the screened Coulomb interaction for realistic interaction parameters as screening induced pairing Eq. ( 4) is most relevant for momenta   where the polarization function χ0 (  ) peaks.To underline this argument, we show the momentum structure of the effective pairing vertex mediated by screened Coulomb interactions in Figure S3.Indeed, we observe that opposite to the RPA approach no amplitude exists at momentum transfers   =  (′) that connects the two valley.

C. Combined mechanism for long-and short-ranged interactions
Due to the different orbital-and momentum structures in the exchange and direct particle-hole channel, it is in general not possible to resum long-ranged interactions in the former channel.This is because the Fourier transform of the Ohno interaction    1  2 (  ) carries a channel-specific transfer momentum   =  −  ′ which is native in the direct particle-hole channel (each interaction line in Figure S4 carries transfer momentum ), but non-native in the exchange channel.Eq. (S7) shows that the native momentum transfer in the exchange channel is   =  +  ′ such that a simple RPA-like resummation with a single momentum transfer can no longer be achieved.In fact, it becomes necessary to unravel the full momentum and orbital dependence of the initial vertex and the susceptibilites, i.e. compute and store    1  2  3  4 (,  ′ , ) and  0  1  2  3  4 (, ).To resolve the small Fermi surface patches in multilayer graphene systems, we sample the entire BZ with   = 7200 2 point and it hence becomes numerically unfeasible to obtain the effective pairing vertex.Therefore, we propose an alternative approach which only accounts for local interactions in the exchange channel and takes the full long-ranged interaction in the direct particle-hole channel, i.e. the best of the previous approaches that can be handled numerically.
As local interactions are accounted for in the exchange channel only, their contribution can be simply ported from the RPA.For the direct particle-hole channel we would like to take the full long-ranged Coulomb interaction as in the RPA.This can not be done trivially as we already know from the RPA appraoch that the resummation within RPA of local interactions must be constrained such that only an even number of loops occur between the in-and out-going legs of the opposite spin vertex Eq. (S7).We can solve this problem by writing the spin-dependence of the susceptibility and the initial vertex explicitly and constrain the effective pairing vertex to the configuration (↑↓↑↓) afterwards.To this end we define where  ∕  1  2 () is the Fourier transform of the Coulomb interaction, but without the local Hubbard- term.The renormalized interaction within RPA is hence given by In the last step, we have restricted the vertex to the opposite spin configuration (↑↓↑↓) such that we recover all terms of the direct particle-hole channel that were present in the RPA and the RPA.

D. Linearized gap equation
To obtain information about the pairing symmetry of the underlying Cooper pairs as well as the superconducting coupling constant, i.e. the critical temperature of the superconducting transition, we solve a linearized gap equation in the projected subspace containing the bands in an energy window  around the Fermi level where V PP  ′ and  PP  ′ denote the band-projected effective pairing vertex and particle-particle susceptibility, respectively.The largest eigenvalue  SC > 0 will lead to the highest transition temperature   and the corresponding eigenfunction Δ(, ) determines the symmetry of the gap, which can be classified according to the irreducible representations of the point group of the normal-state Hamiltonian.The band-projected quantities in Eq. (S12) are defined as In ABCB, only a single band contributes to the Fermi surface on the hole-doped side in the vicinity of the VHS.Therefore, the linearized gap equation Eq. (S12) can effectively be constrained to only capture points on the Fermi surface within the respective band, i.e.only intraband pairing of electrons becomes significant.To this end, we project the effective pairing vertex from orbital to band space Eq.(S13).The vertex V PP  ′ (,  ′ ) hence describes scattering of Cooper pairs with momenta (, −) on Fermi surface patch   to ( ′ , − ′ ) on Fermi surface patch   ′ .Moreover, we can restrict the momentum dependence of the vertex to points on the Fermi surface.In fact, the problem of setting up the pairing vertex reduces to determine proper points on the Fermi surface, i.e. momentum points that satisfy  min <   () <  max .To treat this cutoff problem consistently, we choose  min = −5 meV and  max = −2 meV to capture all fillings on the hole-doped side for which we analyze the superconducting instabilities.In particular, we ensure that Fermi surface broadening may not be smaller than thermal broadening in the system:  FS >  .The explicit contraction of the effective pairing vertex with the particle-particle loop Eq. (S12) will suppress all contributions that are not in the immediate vicinity of the Fermi surface ∝  −1 .
t e x i t s h a 1 _ b a s e 6 4 = " + 3 T o j D K X L 8 w t m D i o W f q N U S 3 D H 4 s= " > A A A C 7 X i c h V F N T 9 t A E H 2 4 Q I H y k Z Z j L y s i J C 5 E d h S V X i q h A h W X S i C R D 4 l A t H a W 1 L J j W / Ym a o i 4 9 w / 0 h r j 2 x r X 9 K / B b O P C 8 O J U K q r L W e m b f v H k 7 s + M m o Z 9 p 2 7 6 b s V 7 N z s 2 / X l h c e r O 8 s r p W e v u u k c W D 1 F N 1 L w 7 j t O X K T I V + p O r a 1 6 F q J a m S f T d U T T f Y y + P N o U o z P 4 5 O 9 C h R Z 3 3 Z i / w L 3 5 O a U K e 0 0 d 5 X o Z Y i E J 9 E T Q S d 7 0 J 2 b N H 2 u r E W j n 0 + 3 q 5 e d U p l u 2 K b J V 4 6 T u G U U a y j u H S P N r q I 4 W G A P h Q i a P o h J D J + p 3 B g I y F 2 h j G x l J 5 v 4 g p X W G L u g C x F h i Q a 8 N / j 6 b R A I 5 5 z z c x k e 7 w l 5 E 6 Z K b D J / c U o u m T n t y r 6 G e 0 D 9 6 X B e v + 9 Y W y U 8 w p H t C 4 V F 4 3 i V + I a 3 8 i Y l t k v m J N a p m f m X W l c 4 K P p x m d 9 i U H y P r 2 / O v u M p M Q C E x E 4 M M w e N V x z H v I F I t o 6 K 8 h f e a I g T M d d W m m s M i p R o S i p l 9 L m r 8 9 6 O G b

FIG. 1 .
FIG. 1. Band structure and density of states (DOS) for different tetralayer graphene structures without an applied electric field.Panels (a-c): Low-energy band structure near valley  with the relative orbital polarization   of the bands indicated by the diverging colormap.The different stacking configurations of ABCB, ABCA and ABAB tetralayer graphene are shown above while the two sites with dominant spectral weight are indicated.Panel (d): DOS within the range of the valley-flat bands.Intrinsic electrical fields caused by broken inversion symmetry in ABCB stacked graphene open a bandgap of Δ  ≈ 20 meV and flatten the low-energy bands around the valleys ,  ′ .Hence, ABCB has the highest DOS near the VHS on the electron-and hole-doped side.The inset highlights different Fermi surfaces around valley  as the chemical potential is tuned through the VHS on the hole-doped side.Starting from the edge of the valence band, the Fermi surface undergoes a Lifshitz transition as the the three pockets originating from trigonal warping continuously transform to an annular Fermi surface.

FIG. 2 .
FIG.2.Spin correlations in the presence of local interactions in ABCB tetralayer graphene.Panel (a) shows the critical on-site interaction  c required for the onset of magnetic order for  close to the van-Hove filling of the valence band (purple).We additionally include the density of states (DOS, gray).We expect an increased tendency towards magnetic order at the van-Hove filling VHS ≈ −3.9 meV.The leading magnetic correlations are displayed in panel (b): Each curve corresponds to the magnetization within one layer.The insets depict the three different magnetic states found: At  ≈  VHS , local interactions drive ferrimagnetism in the C-layer (3).At  <  VHS , antiferromagnetic fluctuations in layer 3 are leading, and for states close to charge neutrality ( >  VHS ) the antiferromagnetic fluctuations extend over all layers.In panel (c) we show the momentum structure of the leading eigenvalue of the susceptibility matrix χ0 () for  <  VHS ,  ≈  VHS , and  >  VHS [c.f.markers in (a)].The transfer momenta supported by the Fermi surface are  ≈ Γ and  ≈  (′) .The -valley (left subpanel) includes the Fermi contour.

FIG. 3 .
FIG. 3. Superconducting instabilities mediated by spin &charge fluctuations in ABCB graphene.Panels (a,b) show the leading SC coupling constant   as function of the chemical potential .Local Hubbard- interactions particularly support spin-fluctuations, which generically enhance spin-triplet  -wave superconductivity due to the allowed intervalley scattering (xRPA).Screening from long-ranged Coulomb interactions (dRPA) instead favors doubly-degenerate  ,wave superconducting order near the VHS and almost degenerate  ∕ + -wave SC order in regions of the three-pocket Fermi surface near the valence band edge.Combining the effect of spin &charge fluctuations (xdRPA), we find that local interactions are screened significantly by the long-ranged Ohno interaction and that SC order similar to the screening scenario (dRPA) prevails and is further amplified by spin-fluctuation near the Stoner instability.For  = 5 eV, (b) the Stoner instability suppresses superconductivity from spinfluctuations due to the presence of ferromagnetic ordering and only screening driven SC order survives.The corresponding instabilities in the valleys  and  ′ are displayed in panel (c).

FIG. 4 .
FIG. 4. Competition of superconducting order at the VHS of ABCB for realistic values of the on-site interaction  separated into contribution from spin-fluctuations and screening.(a) Screening of long-ranged Coulomb interactions (dRPA) generically favors  , (  2 − 2 , )-wave order parameters that change sign within the same valley, while  ∕ + order is suppressed for all values of  .(b) Spinfluctuation exchange (xRPA) is dominated by local interactions and favors spin-triplet  -wave order.(c) Combining screening and spinfluctuations (xdRPA) reveals a delicate interplay of above mechanisms.The degeneracy between ∕-wave is lifted as spin-singlet order is suppressed by spin-fluctuations.Instead,  -wave order is enhanced and eventually exceeds -wave order at the Stoner transition.
FIG.S1.Critical interaction strength (red) and layer magnetization curves for the hole-doped (a) and electron-doped (b) van Hove singularity in ABCB graphene.In contrast to the hole-doped side, the electron-doped carries most magnetization in layer A(1).Apart from that, the qualitative features are equivalent.The dominant ferrimagnetism in layer 1 [panel (a)] corresponds to the layer polarization of the Bloch functions being mostly in layer 1 (cf.Fig.1).
FIG. S2.Critical interaction strength as a function of filling for a range of temperatures and momentum resolutions.Panels (a-c) use the explicit RPA summations and panels (d,e) the flowing RPA.
FIG. S6.Temperature dependence for the different superconducting mechanisms RPA, RPA and RPA.Each panel shows the leading superconducting coupling constant   as function of chemical potential  on the hole-doped side of ABCB.Different columns display the results for different values of the local Hubbard- interactions, whereas different rows show the different mechanisms based on spin fluctuations (RPA), screened Coulomb interactions (RPA) and a combined mechanisms (RPA) that both long and short ranged Coulomb interactions.