Superconductivity and magnetism in the surface states of ABC-stacked multilayer graphene

ABC-stacked multilayer graphene (ABC-MLG) exhibits topological surface flat bands with a divergent density of states, leading to many-body instabilities at charge neutrality. Here, we explore electronic ordering within a mean-field approach with full generic treatment of all spin-isotropic, two-site charge density and spin interactions up to next-nearest neighbor (NNN) sites. We find that surface superconductivity and magnetism are significantly enhanced over bulk values. We find spin-singlet $s$ wave and unconventional NNN bond spin-triplet $f$ wave to be the dominant superconducting pairing symmetries, both with a full energy gap. By establishing the existence of ferromagnetic intra-sublattice interaction, $(J_2<0)$ we conclude that the $f$-wave state is favored in ABC-MLG, in sharp contrast to bulk ABC-graphite where chiral $d$- or $p$-wave states, together with s-wave states, display stronger ordering tendencies albeit not achievable at charge neutrality. We trace this distinctive surface behavior to the strong sublattice polarization of the surface flat bands. We also find competing ferrimagnetic order, fully consistent with density functional theory (DFT) calculations. The magnetic order interpolates between sublattice ferromagnetism and antiferromagnetism, but only with the ratio of the sublattice magnetic moments ($R$) being insensitive to the DFT exchange correlation functional. We finally establish the full phase diagram by constraining the interactions to the $R$-value identified by DFT. We find $f$-wave superconductivity being favored for all weak to moderately strong couplings $J_2$ and as long as $J_2$ is a sufficiently large part of the full interaction mix. Gating ABC-MLG away from charge neutrality further enhances the $f$-wave state over the ferrimagnetic state, establishing ABC-MLG as a strong candidate for $f$-wave superconductivity.


I. INTRODUCTION
Graphene, the single layer version of graphite, has ever since its isolation [1] continued to provide surprises.Recently, electron interaction effects have been at the center of attention, as electronically ordered states, including superconductivity as well as various magnetic or charge ordered or correlated insulating states, have been found in multiple few-layer graphene systems.This includes twisted bilayer and trilayer graphene [2][3][4][5][6], where a small, magic-angle, twist between layers generates a moiré pattern hosting low-energy flat bands [7], which dramatically boosts the low-energy density of states (DOS).With the kinetic energy quenched within these flat bands, interactions effects are bound to become prominent, resulting in the already observed plethora of ordered states.Likewise, bilayer and trilayer graphene strongly biased such that a low-energy partially flat band emerges have recently also been found to host both superconductivity and magnetic or charge ordered states [8][9][10][11][12].While most properties of the ordered states in twisted or biased few-layer graphene systems are still intensively investigated , these results attest to the importance of electron interactions in graphene systems, whenever the low-energy DOS is dra-matically enhanced.
Just beyond the few-layer graphene systems, ABCstacked multilayer graphene (ABC-MLG) [54], a finite stacking sequence of rhombohedral graphite, has notably been known for some time to host topological flat band surface states at zero energy (or charge neutrality) [55,56], also more recently observed in experiments [57,58].The bulk topology generates the flat bands on surfaces perpendicular to the stacking direction, given a sufficient number of layers and thereby clearly appearing in thicker ABC-stack, beyond five to six layers.These surface flat bands generate a large density of states (DOS) at charge neutrality, thus opening the possibility for interaction effects and electronic ordering without any need for twists or bias.While ABC-stacking is not the most common or stable stacking, recent synthesis and characterization have identified both tall and wide area samples of ABC-MLG, even exceeding a dozen layers [57][58][59][60][61][62][63][64].This makes ABC-MLG a natural contender for considering interaction effects in an unperturbed and structurally clean graphene system, beyond the need for externally applied bias or moiré structures.
A straightforward testimony to the power of flat bands for ordering comes from the Stoner criterion for ferromagnetism, stating that a magnetic state forms as soon as the interaction strength exceeds the inverse Fermienergy DOS [65].Magnetic ordering is therefore always a possibility in flat band systems with large zeroenergy DOS, including ABC-MLG, as it only requires weak interactions.In fact, ABC-MLG surfaces have al-ready been found to ferrimagnetically order in ab initio density functional theory (DFT) calculations, with works focusing on systems with 3 to 12 layers [66][67][68][69].Experiments have also found evidence of magnetism in trilayer [70,71], tetralayer [72], and many-layer ABC-MLG [62,64].While this magnetic ordering clearly demonstrates the enhanced ordering tendency, there exist many other possible many-body instabilities, including superconductivity, as also clearly demonstrated by recent results on twisted and biased few-layer graphene [3,9,11].In fact, a flat band generically enhances all ordering channels [73], including superconductivity, implying a wealth of putative states also in ABC-MLG, at least at low temperatures.
In terms of superconductivity, a flat band dispersion has been found to give a linear dependence for the superconducting critical temperature T c on the pairing interaction strength [73][74][75][76], in contrast to the exponentially suppressed BCS result, which thus enables highertemperature superconductivity.However, boosting T c via flat bands has historically been dismissed, since the conventional contribution to the superfluid weight also vanishes when flattening the band dispersion [77].This barrier imposed by theory was recently removed by uncovering an additional geometric contribution to the superfluid weight, induced by non-trivial band topology [78][79][80][81], also present in ABC-MLG [82].
When it comes to the actual possibilities of superconductivity in ABC-MLG, the discoveries of superconductivity in twisted and biased few-layer graphene are very promising.Not only do they show that superconductivity is compatible with a flat band dispersion, but, most importantly, show that all carbon-based material, and graphene systems in particular, have at least one mechanism for producing superconducting pairing [3, 18-20, 83, 84].In fact, superconductivity found in biased trilayer graphene is in an ABC stack, where two distinct superconducting regions are reached by strong gating and displacement field bias [11].While there is as of yet no consensus on the superconducting mechanism, or even the superconducting order parameter symmetry, these results have already attracted considerable theoretical attention [36][37][38][39][40][41][42].However, with three layers being insufficient to develop the topologically protected surface flat bands of ABC-MLG, trilayer graphene requires a strong displacement field bias to produce its partial band flattening leading to superconductivity.
Considerations of superconductivity in ABC-MLG not requiring any twist or bias, but instead utilizing the natural topological flat band regime in thicker stacks, have so far been limited to modeling of conventional (onsite) s-wave superconductivity [73][74][75][76]85], from effective phonon-mediated pairing.While a phonon-mediated pairing mechanism remains plausible, electronic interactions are also known to be important, with estimates from ab initio calculations suggesting strong local and non-local repulsive Coulomb interaction in both single and multilayer graphene systems [86].Additionally, the rich phase diagrams of twisted and biased few-layer graphene, clearly show the importance of electron interactions, likely also for superconductivity.For example, intriguing similarities exists between twisted magic-angle bilayer graphene and the high-T c cuprate superconductors, including characteristic domes of superconductivity next to correlated insulating states, pseudogap physics, and strange metal behavior [2,3], all hinting at an electronic mechanism for superconductivity.With the large surface zero-energy DOS in ABC-MLG, we may expect similar prospects for electronic-driven superconductivity in the surface states of ABC-MLG.
In this work we address the possibility of electronic ordering in the surface flat bands of ABC-MLG systems arising from repulsive electronic interactions.We use a combination of approaches, including fully self-consistent tight-binding mean-field calculations and group theory analysis, as well as complementary DFT calculations.For a fully unbiased treatment, we start from completely general two-site spin-isotropic charge density and spin interactions.We then analyze within mean-field theory all superconducting, charge, and magnetic ordering channels, retaining all channels out to the next-nearestneighbor (NNN) range.
We generally find strongly enhanced surface superconducting and magnetic orders due to the topological flat bands in thick ABC-MLG systems.For superconductivity we find finite-ordering temperatures also for weak interactions for most superconducting symmetries.This is in sharp contrast to both few-layer graphene and bulk ABC-graphite where, due to their semimetallic dispersion, pairing is only possible beyond a critical interaction strength at charge neutrality.In particular, we establish isotropic spin-singlet s-wave pairing and unconventional NNN bond spin-triplet f x(x 2 −y 2 ) wave to be the dominating pairing symmetries, both having a full energy gap.We find that the spin-triplet f -wave state is favored by ferromagnetic NNN bond Heisenberg interactions (J 2 ), while the s wave requires NNN antiferromagnetic interactions or on-site attraction from electron-phonon interactions.With intra-sublattice interactions expected to be ferromagnetic in graphene, i.e.J 2 < 0, and the fwave state also surviving significant on-site repulsion U 0 expected to be present in graphene, while on-site repulsion instead hurts the s-wave state, we conclude that the spin-triplet f -wave pairing is the most likely superconducting symmetry in ABC-MLG.This is very different from bulk ABC-graphite and also few layer graphene, where instead the chiral (d x 2 −y 2 + id xy )-or (p x + ip y )wave states together with an s-wave state display the best ordering tendencies.We trace this different bulk versus surface behavior to the strong sublattice polarization of the surface states.
For putative charge and magnetic ordering, we find a dominating ferrimagnetic state on the surface, also strongly enhanced by the surface flat band and not present in the bulk or few-layer graphene.This is in agreement with earlier DFT calculations, which captures all mean-field charge and magnetic orders (but not superconductivity) [66][67][68][69].We further find that the ferrimagnetic ordering interpolates between pure sublattice ferromagnetism and complete antiferromagnetism, as any of the interaction strengths is tuned from weak to strong coupling.The ferrimagnetic state can therefore be quantified by its sublattice magnetic moment ratio R = |m B1 /m A1 |, with m A1(B1) being the magnetic moments on the surface A (B) sublattice.We also find that our complementary DFT calculations gives highly varying m A1(B1) depending on the choice of exchangecorrelation functional, but their ratio notably stays approximately constant at R ≈ 0.42.
We finally treat the competition between ferrimagnetism and f -wave superconductivity in ABC-MLG.A ferromagnetic J 2 drives both states, but for a large range of weak to moderate J 2 , we find the superconducting f -wave state to have a higher transition temperature.Thus, we find f -wave superconductivity as long as the proportion of J 2 in the full mix of interaction parameters is sufficiently large.We find that superconductivity is further favored over magnetic ordering by gating the system, such that the Fermi level moves away from the flat band.We finally establish the full phase diagram, including gating effects, as a function of relevant interaction parameters by utilizing the constant R-value extracted from the DFT calculations, and find that it contains a substantial portion of f -wave superconductivity.Taken together, our results establish that spin-triplet fwave superconductivity is not only feasible from purely repulsive interactions on the surface of ABC-MLG but even the dominating ordered state in a large portion of the phase diagram.Gating the system away from charge neutrality further enhance the f -wave state compared to a competing ferrimagnetic state.Due to the distinctive properties of the surface flat band in ABC-MLG these results are notably different from ordered states in biased or twisted few layer graphene.Beyond direct predictions for ABC-MLG, our results also clearly illustrate the strong dependence on the normal state behavior for electronic ordering in graphene systems, leading directly to a surprising range of variations of electronic ordering in graphene-based systems.
We organize the remaining of this work as follows.We start in Section II by recapitulating the properties of ABC-MLG, focusing especially on the surface flat bands.Then in Section III we introduce a fully general treatment of all electron-electron interaction parameters up to the NNN range and perform a complete mean-field theory decomposition, taking all possible orders into account.In Section IV we investigate all possible superconducting orders, starting with bulk ABC-graphite as a reference, then focusing on ABC-MLG and its surface superconducting state, elucidating the layer dependence, superconducting pairing symmetries, and superconducting phase diagram containing both spin-singlet s-and spin-triplet f -wave pairing.In Section V we provide estimation for the various interaction parameters and are thereby able to single out the spin-triplet f -wave state as the leading superconducting order.We then move onto charge and magnetic ordering in Section VI where we establish a surface ferrimagnetic state based on our meanfield framework, and show how it can be quantified by the ratio R. We study the resulting competition between the f -wave superconducting and ferrimagnetic state in Section VII.With the help of both mean-field and DFT results we establish the overall phase diagram, also taking gating effects into account, which overall contains a large part with dominating f -wave superconductivity.Finally, in Section VIII we summarize our results.

II. NORMAL STATE ABC-MLG
In ABC-MLG, the graphene layers alternate cyclically between three positions: A, B, and C, see Fig. 1(a).Each graphene layer is in turn bipartite and when viewed from the side, ABC-MLG can be visualized as a staircase arrangement of dimers, see Fig. 1(b), in which the sublattice A atom is located above (or below) the B site atom of the adjacent layer.A consequence of this repeated staircase structure is that infinite ABC-MLG, i.e. rhombohedral-graphite, has a simple unit cell with only two carbon atoms, with an in-plane unit cell length a.To model the normal state electronic structure of ABC-MLG, we use a tight-binding Hamiltonian.The energy scale is foremost set by the intralayer nearestneighbor (NN) hopping amplitude t = 3.0 eV within each graphene layer.Thereafter follows the interlayer hopping between different sublattices of neighboring layers t ⊥ = 0.1t = 0.3 eV, since each carbon atom has another atom directly above or below it.A Hamiltonian consisting only of these two hopping amplitudes is already a good approximation for the non-interacting normal state electronic structure of ABC-MLG, as it captures the bulk topology and the associated topological surface flat bands [56,74].For N L layers, this Hamiltonian reads [56,73,87]: where a † iLα creates an electron with spin α at site i in layer L = 1, . . ., N L and neighboring sites (layers) are denoted by ⟨• • • ⟩.The number of electrons is regulated by the chemical potential µ.
The low-energy electronic structure of normal state ABC-MLG originates from the single-layer graphene Dirac cones found at the K-and K ′ corner points of the graphene Brillouin zone.In ABC-MLG, the interlayer hybridization t ⊥ transforms the individual graphene Dirac points into Fermi spirals, twisting around the crystal momentum axis of the stacking direction k z , connecting the H − K − H and H ′ − K ′ − H ′ points [55].Asso- ciated with these Fermi spirals is a topological winding number with the consequence that a topological surface flat band is formed within the projection of the Fermi spirals on each stack surface [55,56,88].This mimics the physics of the one-dimensional Su-Schrieffer-Heeger (SSH)-model, giving ABC-MLG very different properties compared to other graphene stackings [55,56,67,89].A finite stack of ABC-MLG thus hosts two topological flat bands, one on each perpendicular surface, which are located around the K-and K ′ -points and tied to the charge neutral zero energy point (the graphene Dirac point).This is illustrated in Fig. 1(c), where we plot the lowest positive energy electronic band of eight layer ABC-MLG.The vanishing dispersion of the flat bands is further seen in Fig. 1(d), where we plot the band structure along a line cut passing through one of the K-points, ξ (0, k y ), as indicated by the dashed line in Fig. 1(c).
The flat dispersion results in a significantly enhanced DOS around zero energy, seen by plotting the local density of states (LDOS) in Fig. 1(e) for the same eight layer ABC-MLG.The large difference between the surface (cyan and purple) and third (red and green) layer LDOS clearly display the surface localization, as well as a strong sublattice dependence.Resolving the states further reveals that the flat bands exactly at the two K points are almost completely localized on of the outer surface sites A1 and B8, marked by X in Fig. 1(b).Moving slightly away from the K points, other sites also contribute, shown by using the band color in Fig. 1(d) to indicate the A1 site probability density.
We further quantify the sublattice layer dependence by integrating the site-dependent LDOS in a narrow energy window ξ ∈ [E F − 3 meV, E F + 3 meV] around zero energy, where E F denotes the Fermi energy.We plot the resulting integrated density, P, as a function of layer in Fig. 1(f) for the two sublattices A and B. Together Figs.1(d-f) demonstrate that the topological surface states are strongly localized to the two surfaces and that the two surfaces have a near complete sublattice polarization that is reversed between them.We note that this sublattice polarization clearly sets many layer ABC-MLG apart from ABC-graphite, twisted bilayer graphene, and also to a large extent from few layer ABC-MLG where the sublattice polarization is non-existent or much less pronounced.
Due to their topological origin, the surface flat bands in ABC-MLG depends on having a sufficient number of layers in the stack to establish the bulk topology.To track the development of the surface flat band, we plot the flat band radius against the number of layers in Fig. 1(g), where we define the radius to be the length of the wavevector around the K-point for which the band energy is below a fixed threshold, using the same energy window as in Figs.1(e) and 1(f).As seen, the flat band radius grows monotonically with N L , slowly saturating to the projected Fermi spiral area.Thus, there is a rapid increase in the DOS, which is proportional to the flat band area, with each added layer of ABC-MLG when going from only a few layers to thick samples.
Finally, we comment that in the NN hopping model of Eq. ( 1), the bulk topology of ABC-MLG depends on particle-hole symmetry, while longer-ranged intrasublattice hopping terms break this symmetry.However, such longer range hopping is relatively small and adds only a small curvature to the flat bands.In fact, the surface flat bands have been shown to clearly develop also in DFT calculations at both surfaces as well as interfaces of ABC-MLG [67], showing that the essential features are not dependent on the assumptions underpinning Eq. (1).

III. ELECTRON INTERACTIONS
The large surface DOS makes ABC-MLG very susceptible to develop a variety of ordered states from interactions, a select few of which have already been explored [55,56,73].Here we proceed within a more systematic approach.We first note that graphene and graphite systems [86,90], including ABC-MLG [91,92], have been found to have strong local and non-local Coulomb interactions.We therefore analyze all possible symmetry-breaking orders, including superconductivity, arising from local and non-local electronic interactions in ABC-MLG, by adding these interactions to the noninteracting normal state model of Eq. ( 1).This allows us to determine the leading superconducting and magnetic instabilities, as well as their competition.
We keep our analysis generic and therefore include all spin-isotropic two-site effective density-density and spinspin interactions within each graphene layer.We are able to limit our analysis to these intralayer interactions, since the individual graphene layers of ABC-MLG are only coupled with weak van der Waals forces.Since the inplane interaction potential decays with distance [86], we are also able to restrict ourselves to relatively short-range interactions.Still, we note that non-local interactions are important for ordering in ABC-MLG.In particular, we find that the strong surface sublattice polarization of the ABC-MLG surface states means that the next-nearestneighbor (NNN) range is in fact the shortest interaction range beyond on-site interactions that connects carbon sites with a significant low-energy DOS, thereby producing a strong ordering susceptibility.For this reason, we consider interaction ranges up to NNN to capture these important effects.Thus, we extend the non-interacting normal state model of Eq. ( 1) by adding an interaction Hamiltonian H I of the form: Here the matrices σ 0 and σ ν , ν = x, y, z, are the 2 × 2 identity matrix and Pauli matrices in the spin degrees of freedom (indexed by α, β, γ, δ), respectively, which we also collect in the four-vector σ η = {σ 0 , σ ν }.The total interaction potentials Γ ij between sites i and j can always be divided into two distinct parts: U ij being the effective density-density Coulomb interactions and J ij being the effective isotropic Heisenberg spin-spin (or exchange) interactions.This captures all two-site electron interactions and we will consider all such terms up to NNN distance.

A. Mean-field decomposition
To proceed we perform a complete Hartree-Fock-Bogoliubov mean-field decomposition of the total Hamiltonian H = H 0 + H I .This allows us to access all possible orders, superconducting as well as charge and magnetic orders, in the ABC-MLG.Doing so we arrive at the mean-field interacting Hamiltonian: where In order of appearance, the terms in the mean-field Hamiltonian (3) are: the particle-particle (superconducting or Bogoliubov) channel with the order parameters ∆ η jiL (here ση = iσ y σ η ), the direct particle-hole (Hartree) channel with order parameters m η ijL , and the exchange particle-hole (Fock) channel with the order parameters mη ijL .The order parameters in each channel are defined self-consistently through Eqs.(4).For each channel, the total coupling strengths, Γ η,± ij are linear combinations of the effective interaction strengths U ij and J ij , as tabulated in Table I.In the table we divide the channels further into spin-singlet and spin-triplet configurations, Table I.Total mean-field coupling strengths in Eq. ( 4) for both superconducting and particle-hole (PH) channels in terms of the effective density-density Coulomb interactions Uij and isotropic Heisenberg spin-spin interactions Jij in Eq. ( 2).
based on the sigma matrices.Note that in Eq. ( 4), nonlocal interactions connect order parameters at different sites.Additionally, the Bogoliubov and Fock channels, but not the Hartree channel, contain bond order parameters that also connect two different sites directly in the operator expectation value.The former is illustrated in Fig. 2(a).

B. Linearized gap equation
The complete mean-field decoupling of Eq. ( 3) reduces the interacting system to one of non-interacting quasiparticles moving through the collective mean-field potentials defined by the non-linear self-consistency equations of Eq. ( 4).These mean fields will often only renormalize the quasiparticle properties.There is, however, also the possibility that the interactions produce a qualitative different state through spontaneous symmetry breaking.Since we seek to analyze the possibility of such symmetrybreaking orders, we hereafter assume that all renormalizations have already been included in the normal state Hamiltonian (1) and its parameters, which therefore corresponds to what would be measured as the normal state above the ordering temperatures, as opposed to the bare parameters.
The full nonlinear self-consistency equations of Eq. ( 4) are valid at any temperature as they completely define each order parameter.However, the order parameter of any symmetry breaking order only become finite below a certain critical transition temperature, T c .Right at T c the order parameter is therefore vanishingly small.This allows the self-consistency of Eq. ( 4) to be linearized by a perturbation expansion in H MF (see Appendix A for details).After such linearization we arrive at the generic linearized gap equation (LGE) [73]: with the order parameters ∆ η collected in a column vector D η,+ , and similarly m η and mη collected in D η,− .Here, ī, j label the bands of the normal state, H 0 .In particular, because only the superconducting orders breaks the U(1) gauge symmetry, there is no coupling between the superconducting particle-particle (+ superscript) and particle-hole, i.e., the density and magnetic, (− superscript) orders and thus Eq. ( 5) is in fact two equations for these two separate channels.Central to the LGE is M η,± , the stability matrix, which is the temperature dependent response matrix.Its matrix elements are a combination of the interaction coupling strengths Γ η,± ij , stored in Γ η,± , the structure factors expressed by the symmetry matrix S η,± (derived in Appendix A), and the temperature-dependent particleparticle and particle-hole susceptibilities: where β = (k B T ) −1 is the inverse temperature, k B is the Boltzmann constant, and ξ ī are the band energies of the normal state H 0 .
For a given set of interactions and at an arbitrary temperature, the LGE Eq. ( 5) does not generically have a solution.Rather, Eq. ( 5) can only be satisfied when the stability matrix M has at least one eigenvalue that is unity.At high temperatures, all the eigenvalues are generally suppressed by the temperature factor (β) in the susceptibility, Eq. ( 6).With decreasing temperatures, the eigenvalues tend to grow.The temperature where an eigenvalue eventually reaches unity marks the critical temperature T c where Eq. ( 5) then has a solution and an ordering transition takes place.The type of order and its spatial pattern is completely specified by that solution's corresponding eigenvector D η,± .The LGE Eq. ( 5) is therefore an eigenvalue problem of the temperature dependent stability matrix, M η,± , and by looking at the eigenvalues of M it is possible to determine the order competition.The order with the largest eigenvalue becomes the dominant order, but can also be followed by sub-dominant orders, arising at lower temperatures.Two orders with similar eigenvalues can likewise be viewed as directly competing.

C. Order parameter symmetries
The LGE Eq. ( 5) has the same symmetry as the normal state Hamiltonian Eq. ( 1) and therefore the stability matrix M is block diagonal in a symmetry basis of the irreducible representations (irreps) of the normal state Hamiltonian.Only below T c is the symmetry reduced, when a symmetry-breaking order parameter becomes finite.
Since the normal state Hamiltonian Eq.( 1) is spinisotropic, the mean-field channels first split into two spinsymmetry sectors, spin-singlet (η = 0) and spin-triplet (η = ν).The spin-triplet sector is further completely degenerate, and we can thus here focus only on the ẑ channel, with no loss of generality.With only in-plane interactions, we can resolve the order parameters as functions of the in-plane bond ranges by writing any parameter as D L,ij = D L,i,i+r = D L,r , with r = 0, 1, 2 for on-site (ON), nearest-neighbor (NN), and next nearest neighbor (NNN) orders, respectively.Because each of these bond ranges is a close set under the symmetry transformations of the normal state Hamiltonian H 0 , we can construct a real space symmetry basis for all orders labeled by spin sector η, bond range r, and irreps of the lattice symmetry of H 0 .Consequently, any order parameter D η,± can be written as a linear combination in this basis: where χ labels all the in-plane irreps in layer L. Since any solution to Eq. ( 5) has a definite symmetry, the nonzero amplitudes A η L can only belong to one irrep, which is then the same in all layers of a ABC-MLG stack.We can therefore drop the L subscript.The real-space symmetry basis states in Eq. ( 7) also translate to a corresponding expression in reciprocal space: where k is the in-plane crystal momentum in the first Brillouin zone and Dη,χ,± is the on-site or bond amplitude on the in-plane on-site or bond vector g (r) n , as illustrated in Fig. 2(a) for sublattice A.
We concretize this general derivation by tabulating in Table II all possible basis states in Eq. ( 7) for the superconducting order possibilities, using the notation ⃗ ∆ η r for Dη,χ,+ r , where the rows represent the possible irreps χ for each range r and divided into the two different spin sectors.Note that in the spin-singlet sector only even irreps (gerade) are allowed, while spin-triplet only allows for odd representations (ungerade), due to the Fermi-Dirac statistics of the superconducting Cooper pairs.For each irrep and bond range, we extract the corresponding real-space structure for ⃗ ∆ η r (second column).In reciprocal space we tabulate the lowest order expansion around the Brillouin zone center Γ (third column).This choice of high-symmetry Brillouin zone point makes it possible to directly match up the resulting basis function to the irreps of the D 6h point group of the honeycomb lattice (fourth column) [83,[93][94][95][96][97].Other reciprocal space choices, such as the K point, do not change the overall decomposition but can mix up the gerade (ungerade) and spin singlet (spin triplet) as they are not centered at zero momentum.
Our goal next is to analyze which of the ordering channels have the highest T c and are thus dominant, starting with superconductivity.Before proceeding we note that below T c all multidimensional irreps generally split [97], with the dominating channel being the one with the lowest energy determined by higher order (non-linear) terms in the self-consistent equations, Eq. ( 4).Moreover, below T c it is in principle possible to also produce co-existing orders from different spin channels and irreps, which can only be accessed by solving the full non-linear selfconsistency equations (4).In practice this is, however, seldom seen and then only at notably lower temperatures, and we will therefore for now ignore that possibility and focus on the order with the highest critical temperature.

IV. SUPERCONDUCTIVITY
Having performed a full symmetry analysis of all possible order parameter symmetries, we next determine which superconducting orders can actually be realized in ABC-stacked graphene.We consider here both finite and infinite ABC-stacked graphene.The infinite-layer case, corresponding to bulk ABC-graphite, serves as a baseline to separate bulk from surface contributions.Since the spin-singlet and -triplet sectors never mix at T c , we consider each spin sector separately.To clearly identify each contribution and possible order, we at first also turn off the coupling between the different pairing range r (ON, NN, NNN), treating them also separately before a joint analysis in Sec.IV B 3. Note however that spin-triplet ON ordering is forbidden by Fermi-Dirac statistics, resulting in total five possible ordering channels, all tabulated in Table II.For each of these channels, we solve Eq. ( 5), using 320×320 in-plane k-points to find both the complete spatial structure of the leading superconducting order parameter ∆ η = D η,+ , i.e. both layer distribution and in-plane irrep, and its T c as a function of the coupling strength Γ η,+ r in that particular channel, given by Table I.

A. Bulk ABC-stacked graphite
We start by considering bulk ABC-stacked graphite.For the pristine situation with µ = 0 we plot T c for each channel in Fig. 2(b) (dashed lines) as a function of the absolute value of that channel's coupling strength Γ η,+ r .As seen, in each pairing channel there is no superconductivity below a channel dependent critical coupling strength, Γ η,+ r,c .Estimates of the critical coupling strengths can be extracted from Fig. 2(b) as: ,c |} = {1.8t,0.6t} for the spin-triplet channel.Assuming bulk ABC-stacked graphite to not be superconducting, these critical coupling strengths set a natural upper limit for the values we later investigate for the coupling strengths.It is the nodal electronic structure of ABC-graphite with its vanishing DOS at µ = 0 that produces these critical coupling strengths, thus with no superconducting states possible for any interactions lower than these critical couplings.This is in contrast to the standard BCS relationship, which produces a finite (but possibly very small) T c for any infinitesimal attraction.The situation in bulk ABC-stacked graphite is analogous to that of single layer graphene where the Dirac cones similarly form a nodal electronic structure.To emphasize this similarity, we find that the spin-singlet NN channel critical coupling strength |Γ 0,+ 1,c | ≈ 1.9t is almost identical to that of single layer graphene reported before [83,93].
Since we find that ABC-stacked graphite has large critical coupling strengths in all superconducting channels, it means that pristine ABC-stacked graphite is by itself an unlikely host of superconductivity.To later contrast with finite ABC-MLG, we however still report the superconducting symmetries above the critical coupling strengths in ABC-stacked graphite and summarize the result in the bulk prevalence column of Table II.First, the spin-singlet ON channel has necessary s-wave spatial symmetry.Next, we find that the leading order in the spin-singlet NN channel is the two-fold degenerate d-wave solution belonging to the E 2g representation.At µ = 0 this d-wave solution is fully degenerate with the rotationally symmetric s ext -wave solution (A 1g ), whereas the d-wave solution becomes dominant for any finite µ, just as in single-layer graphene [83,93].In contrast, in the spin-singlet NNN channel, the d-wave order is always the dominating solution over the symmetric s ext -wave solution, with further increasing dominance at finite µ.Both the ON and NNN pairing have equal order parameter magnitudes on both sublattices.For spin-triplet superconductivity we find that the two-fold degenerate p-wave solution of the E 1u representation to be dominant in both the NN and NNN channels.

B. ABC-MLG
In contrast to bulk ABC-graphite, the topological flat bands on the surfaces of sufficiently thick stacks of ABC-MLG drastically increase the possibilities for superconductivity.In Fig. 2(b) we also plot the critical temperatures T c of the superconducting channels for three (dotted) and eight (solid) layer ABC-MLG as a function of the coupling strengths Γ η,+ r .It is immediately evident that the onset of superconductivity, at any specified temperature, occurs at much lower coupling strengths compared to the bulk ABC-stacked graphite results (dashed) for almost all symmetries.In particular, for thick slabs, we find that the surface states initially generate a linear increase in T c for weak coupling strengths, which is consistent with the T c dependence earlier derived for flat band systems [73,75].For increasing coupling strengths, T c instead starts grows exponentially with the channel coupling strength, leading to a more traditional BCS-like increase in T c .This is expected, as for larger strengths, also higher energy states besides the flat band states begin to contribute to the susceptibility and hence the superconducting T c .The significant enhancement of the superconducting susceptibilities channels becomes particularly clear when we in Fig. 2(c) plot T c as a function of the normalized coupling strengths Γη,+ r = Γ η,+ r /Γ η,+ r,c for each channel.By doing this normalization with respect to the bulk ABC-stacked graphite critical couplings, we directly quantify the large effect of the surface states in ABC-MLG, including clearly illustrating how eight lay-ers are much more prone to superconductivity than three layers.This strong surface enhancement of superconductivity is the first important result of this work.
The notable exception to the strong enhancement of superconductivity in ABC-MLG is for NN range pairing, for which there is instead no change compared to the bulk, independent on spin channel (blue and green curves).This is in sharp contrast to both doped graphene and magic-angle twisted bilayer graphene, which both have been shown to host enhanced NN pairing compared to the pristine graphene baseline, resulting in either chiral (d + id) or nematic d-wave superconducting pairing [27,34,83].We attribute the lack of any possibility for d-wave pairing in ABC-MLG to its strong, near complete, sublattice polarization of the surface flat bands (see Fig. 1).This sublattice polarization makes it essentially impossible to have NN, or any inter-sublattice, superconducting pairing.

Layer dependence
We next quantify the layer and sublattice dependence of the superconducting states.In Figs.2(d We here illustrate the behavior by plotting the NNN spintriplet state, but find similar behavior also for the spinsinglet ON and NNN states.In the weak coupling regime, which we define as 0 < |Γ η,+ 0,2 | ≪ |Γ η,+ 0,2,c |, i.e., very weak coupling relative to the critical coupling in bulk ABCstacked graphite, we find superconductivity only localized to the outer surface layers.In fact, we find that the order parameter magnitude generally follows the surface state LDOS [see Fig. 1(f)].Additionally, as seen in Fig. 2(d), the order parameter has a large sublattice staggering within each layer that reflects the strong normal state sublattice polarization.This means that both the LDOS and superconducting state live essentially exclusively on one sublattice on each surface.Increasing the coupling strength and the system enters an intermediate regime, 0 ≪ |Γ η,+ 0,2 | ≲ 0.85|Γ η,+ 0,2c |, where the superconducting order parameters on the two outer surfaces begin to couple, even for an eight layer stack [see Fig. 2(e)].Simultaneously, the order parameter also begins to spread to both sublattices, even if it maintains a large sublattice staggering.Finally, in the strong coupling regime, the superconducting order parameter becomes finite throughout the entire ABC-MLG stack [see Fig. 2(f)].This happens because the coupling strength approaches that of the critical coupling of the bulk.As seen, for eight layers this occurs around |Γ η,+ 0,2 | ≈ 0.85|Γ η,+ 0,2c |, while thicker stacks naturally requires a coupling strength closer to the bulk critical coupling for the full stack stack to become superconducting.The results in Figs.2(d) -2(f) clearly illustrate the large effect of the surface states on superconductivity.ABC-MLG does however require a certain thickness before the topological surface flat band states are fully formed, see Fig. 1(g).Thus, T c is also expected to grow with additional layers in ABC-MLG.To illustrate this stack thickness dependence, we plot the T c of all three enhanced superconducting channels: spin-singlet ON (ONs), NNN (NNNs), and spin-triplet NNN (NNNt) pairing, all as a function of the number of layers in Fig. 3. To isolate the layers dependence from other factors, we fix the coupling strengths Γ η,+ r of each channel to produce a fixed T c of either 20K or 10K for a 20 layer stack.As seen, all superconducting channels display the same T c -layer dependence.Thus, the topological flat band is equally important for enhancing all superconducting tendencies in ABC-MLG, independent on both the superconducting pairing range and spatial symmetry.This complements our first important result of this work, stating that the surface states enhance all superconducting states equally.Moreover, to achieve a higher T c , thicker slabs are generally required to reach saturation towards the maximum T c for surface superconductivity.We also note in Fig. 3 that for ABC-MLG with three layers or less we find only a vanishingly small T c , consistent with essentially no flat band formed for such thin stacks.Thus, to find superconductivity in pristine (undoped and unbiased) ABC-MLG, we need to consider stacks substantially thicker than trilayer graphene.

Superconducting symmetries
Having shown the strongly enhanced critical temperature in ABC-MLG, we next look at the order symmetry in each of the achievable surface superconducting channels, i.e., ON and NNN pairing.In the spin-singlet channel, we find that the surface flat bands support either on-site s on -wave or NNN extended s ext -wave symmetry super-  conductivity.Both of these states are isotropic, preserving the full symmetry of the normal state Hamiltonian, Eq. ( 1), thus belonging to the A 1g irrep of the D 6h symmetry group.They also both have a fully gapped energy spectrum.For the s on -wave state the gap is constant.In contrast, for the s ext -wave state, we display the reciprocal space form factor ∆ 0 2 (k) in Fig. 4(a), illustrating a nodal line encircling the Γ point.Yet, close to charge neutrality, the Fermi surface is far away from this nodal line and the spectrum is thus fully gapped.
In the spin-triplet channel, here only for the NNN pairing range, we find that the dominant order is an unconventional f x(x 2 −3y 2 ) -wave symmetry superconducting state belonging to the B 2u irrep.When measured against the normalized coupling strength Γη,+ r in Fig. 2(c), we find that this f -wave state is the most enhanced superconducting channel, giving the highest T c for all normalized coupling strengths.One explanation for this significant enhancement of the f -wave state is that it is also a fully gapped superconducting state, for all filling factors around the charge neutral level.The full gap of the fwave state can be traced back to the fact that the three nodal lines of the pairing potential, shown in Fig. 4  around the K and K ′ -points.This full gap results in a large superconducting condensation energy, which in turn makes the f -wave energetically favorable.To summarize these results, we tabulate the symmetries of the leading superconducting states for ABC-MLG in each channel in Table II in the surface prevalence column.We note the dominant surface superconducting symmetries are essentially completely opposite to results in the bulk.While the bulk will host NN or NNN spin-singlet d-wave or spin-triplet p-wave states, the surface will only host NNN spin-singlet s ext or spin-triplet f -wave states.
In fact, the only dominant superconducting symmetries appearing for both bulk and surface is the on-site spinsinglet s-wave state, which is to be expected since it does not have any real-or reciprocal space structure.This distinct difference between surface and bulk superconductivity is surprising and a second major result of this work.
With the f -wave state being the most enhanced surface superconducting state, we next investigate it further by additionally solving self-consistently for spin-triplet NNN pairing at zero temperature T = 0 for an eight layer stack.For this we use Eq.(4a) within the Bogoliubovde-Gennes (BdG) formalism [98], which then gives an excellent complement to the LGE results in Fig. 2 representing the solution at T c .In Fig. 5(a) we show the T = 0 self-consistent result for the amplitude |∆ z 2 | on the surface (black) and middle (red) layers as a function of coupling strength Γ z,+ 2 with the full-layer profile in the inset, which fully corroborates the results at T c Moreover, tracing the order parameter as a function of Γ z,+ 2 , we find that it increases monotonically everywhere, but for 0 < |Γ z,+ 2 ≲ |0.4t the order parameter is only localized to the surfaces, illustrating the strongly enhanced surface superconductivity.(The slight difference between self-consistent upper coupling bound |Γ z,+ 2 | ≈ 0.42t and that of the LGE, |Γ z,+ 2 | ≈ 0.85|Γ z,+ 2,c | = 0.51t obtained in Section IV B 1, stems from the inclusion of all higher order terms in self-consistency calculation while LGE in contrast contains only first order terms.)The f -wave superconducting state also opens a superconducting gap on the surfaces, splitting the surface flat bands, as illustrated in Fig. 5(b), where we plot the BdG quasi-particle spectrum E (k).The color code shows the probability density of the bands on site A1 (same result for B8, thus summing to 1), which show that the gapped flat band states fully reside on the outer layers with essentially full sublattice polarization.Finally, in Fig. 5(c), we show the DOS as a function of energy, E, for the normal state (blue), and the f -wave superconducting state at |Γ z,+ 2 | = 0.2t, 0.4t (red, black).While the normal state DOS is gapless with large surface DOS at E = 0, the superconducting DOS has a finite spectral gap which increases with Γ z,+ 2 .

Phase diagram
In the previous subsections, we showed that the topological flat band surface states of ABC-MLG generate an enhanced susceptibility in the spin-singlet on-site s-wave and NNN extended s-wave channels, as well as in the unconventional spin-triplet NNN f -wave superconducting channel.So far, we have however only considered each pairing-range r and spin-channel separately and also regarded the channel coupling strengths Γ η,+ r as independent system parameters.While this is appropriate for a general analysis of the available possibilities, we note that the different ordering channels are generally not independent.Rather, the coupling strengths Γ η,+ r are linear combinations of the effective interactions, U r and J r , see Table I.Next, we therefore analyze all pairing channels jointly to show their possible couplings and dependence on the effective interactions J r and U r in order to extract the phase diagram.
When considering the pairing channels together, they are also allowed to couple in the LGE, especially mixing different ranges.We however find that the NN pairing ranges do not couple sufficiently to the other pairing ranges to affect the analysis for interaction weaker than their corresponding bulk critical coupling strength, which we use as an upper limit of parameters that we consider.We therefore exclude the NN pairing channels from the rest of our analysis.The NNN spin-triplet range is therefore the only spin-triplet channel left in our analysis, and with spin-channels not mixing with each other, we can treat the NNN spin-triplet channel as uncoupled.In contrast, the two spin-singlet s-wave pairing states, the ON s on -wave and NNN s ext -wave states, couple even to first order in the LGE in a joint analysis, since both states belong to the same A 1g symmetry, see Table II.
In Fig. 6 we analyze the phase diagram for the relevant combinations of U r and J r .We start in Fig. 6(a) by plotting the joint T c produced by the two spin-singlet s-wave states as a function of the two normalized coupling strengths Γ0,+ 0 and Γ0,+ 2 .We here use the normalized strength to be able to cut off the phase diagram at Γ0,+ 0,2 = 1 as that corresponds to the bulk critical coupling strength.We find that the spin-singlet s-wave state is favored by either an effective on-site attraction, U 0 < 0, a NNN attraction U 2 < 0, or an antiferromagnetic NNN interaction J 2 > 0. These values all produce the highest T c , appearing in the lower right corner of the phase diagram.Moreover, Fig. 6(a) also shows that a net attraction for one of the pairing ranges, i.e.U 0 or U 2 , can overcome a net repulsion in the other range, i.e.U 2 or U 0 , to still produce a finite T c .
To draw a complete phase diagram, the spin-singlet T c of Fig. 6(a) has to be compared with the T c of the spintriplet channel.In Figs.6(b-d), we overlay the T c contours of both the spin-singlet s-wave and the spin-triplet f -wave states for three different fixed values of the on-site interaction U 0 and as a function of the remaining interactions J 2 and U 2 .For U 0 = 0, shown in Fig. 6(c), the dominant superconducting order is directly determined by the sign of J 2 .A ferromagnetic NNN interaction, J 2 < 0, favors spin-triplet f -wave pairing, while an antiferromagnetic, J 2 > 0, favors the spin-singlet s-wave pairing.A repulsive U 2 > 0 reduces the total coupling constant of both pairing channels, which produces a wedge-shaped region with boundaries Γ z,+ 2 = 0 and Γ 0,+ 2 = 0 (hatchedmarked region) inside of which both the spin-singlet and -triplet channels have a net repulsion, leading to a phase space region with no pairing.
Adding a repulsive on-site term U 0 = t in Fig. 6(d), we find that the spin-singlet s-wave state is suppressed.The T c of the f -wave is, however, completely unaffected by U 0 , since by symmetry the spin-triplet f -wave pair amplitude vanishes on-site and thus completely avoids any on-site repulsion.The result is that the spin-singlet state is suppressed and the non-superconducting wedgeshaped region grows at the expense of the spin-singlet state region.We note that, technically, this wedge shape region actually splits in two.Inside the original wedge, both channels have a net repulsion, devoid of all pairing.In contrast, between the old and new wedge-boundary (light blue hatched region), one channel is still technically attractive, which in principle gives a finite T c , but due to the strong on-site repulsion, this T c is vanishingly small.To be able to draw a realistic phase diagram in Fig. 6(d), we have included a cut-off of T c ≤ 10 −2 K to define this second boundary.
Finally, if we instead add an attractive U 0 = −t, as in Fig. 6(b), then the T c of the spin-singlet channel increases.Such attractive U 0 can be viewed as the effective interaction arising from putative electron-phonon coupling in ABC-MLG graphite and is, as such, still relevant to consider.Since the ON spin-singlet channel is now attractive everywhere in the phase diagram, the entire prior non-superconducting wedge is now transformed into a region with finite but vanishingly small T c , indicated by the light blue region in Fig. 6(c).Moreover, the attractive on-site interaction further shifts the wedge shaped region which now also encroaches on the f -wave region which consequently shrinks.Together, Figs.6(a) -6(d) give a complete picture of the superconducting phase diagram for all interactions out to the NNN coupling range, also including the effect of electron-phonon attraction.We find that the spin-singlet s wave (jointly ON and NNN pairing) and the spin-triplet f wave (NNN pairing) are both present even for repulsive Coulomb interactions U 0 and U 2 and that the sign of the Heisenberg interaction J 2 determines which of these states is materialized.This constitutes a third important result in this work.

V. INTERACTION ESTIMATES
In the previous section, we established the generic superconducting phase diagram of ABC-MLG.We on purpose kept the discussion general, without assuming any specific interaction strengths in order to show the full range of possibilities, finding that electronic interactions in ABC-MLG can support both spin-singlet isotropic swave and unconventional spin-triplet f -wave superconductivity depending on the interactions.Here we aim to give an estimate of the actual interaction strengths in ABC-MLG, in order to finally discern the most likely superconducting state in ABC-MLG.
We note first that the effective density-density Coulomb interactions have already been calculated from ab initio using the constrained random phase approximation (cRPA) [86], which find a strong on-site interaction in both single-layer free-standing graphene (U 0 = 9.3 eV) and in graphite (U 0 = 8.1 eV).Similar estimates were later also found in a DFT+U + V framework, for both graphene and graphite (U 0 = 7.6 eV) [90].These estimates produce a dielectric constant in graphene close to the experimentally observed value [86,90] and are also comparable to values (U 0 ∼ 10) eV used in quantum chemistry models of hydrocarbon and the Pariser-Parr-Pople model [99,100].
The very strong on-site repulsive Coulomb interactions found in graphene and graphite significantly impact the range of likely ordering also in ABC-MLG.Foremost, it suggests that the effective Coulomb interactions is repulsive and therefore significantly detrimental to the isotropic spin-singlet pairing in ABC-MLG, as found in Fig. 6.Conventionally, the main pairing mechanism for the isotropic spin-singlet pairing channel would be an effective electron-phonon coupling, which here is effectively modeled as an attractive U 0 < 0. We note, however that in single-layer graphene the strong out-of-plane vibrations do not couple to the out-of-plane p z orbitals of the low-energy electronic structure to first order [101,102].There is likewise no coupling between the out-of-plane p z orbitals and in-plane vibrational modes due to symmetry [103].Consequently, the effective dimensionless electron-phonon pairing is small in graphene and any Coulomb repulsion has further been shown to dramatically suppress both the pairing and the critical temperature of phonon-driven superconductivity in single layer graphene [102].Despite this lack of electron-phonon coupling in graphene, we nonetheless include the possibility of an effective on-site attraction U 0 < 0 in Fig. 6(b), since there exist additional phonon modes in ABC-MLG due to the additional layers.For instance, in-plane acoustic longitudinal phonon modes have been suggested as a mechanism for superconductivity in both trilayer and tetralayer ABC graphene [37,48].This might thus open for phonon-driven s-wave superconductivity also in ABC-MLG, although any such electron-phonon coupling would have to first overcome the strong repulsive Coulomb interaction before generating superconducting pairing.
Beyond suppressing any putative s-wave state, the strong on-site repulsive Coulomb interactions also have large implications for the exchange interactions in ABC-MLG.The direct Coulomb exchange interaction produces a ferromagnetic exchange between orbital sites [92,104].On the other hand, the strong coupling limit of the onsite Coulomb repulsion U 0 is captured by the t-J model [105], featuring an antiferromagnetic super-exchange on NN bonds.Thus, while the NN direct exchange in graphene has been calculated to be ferromagnetic with J NN,FM ≈ 1.6 eV, using ab initio cRPA, this contribution is dwarfed by the antiferromagnetic NN superexchanged J NN,AF = 4t 2 /U ≈ 3.5 eV [106].We can thus here estimate the net effective Heisenberg interaction J 1 = J NN,AF − J NN,FM ≈ 1.9 eV to be antiferromagnetic on NN bonds [106].
In contrast, the NNN bonds super-exchange is suppressed by the smaller NNN hopping amplitude t 2 ∼ 0.1t [107,108], which can thus be expected to be 100 times smaller than the NN super-exchange.The NNN direct ferromagnetic Coulomb exchange remains finite but is also similarly small [104].The question here remains as to their balance.From a basic point of view, the bipartite lattice structure of graphene has been shown to dictate the sign of the bare spin susceptibility with an alternating sign between the two sublattices [109].Due to the role the bare spin susceptibility plays in generating effective interactions, e.g., in the random phase and fluctuation-exchange approximation [110][111][112], the effective exchange is expected to inherit this alternating sign structure of the bare susceptibility.Consequently, we expect the NN effective exchange J 1 to be antiferromagnetic J 1 > 0, just as we concluded in the previous paragraph, while we expect the NNN term to be ferromagnetic, J 2 < 0. Corroborating this basic expectation, the effective interaction generated in the functional renormalization group (fRG) flow has been shown to have an alternating sign structure in single layer, bilayer, and in both ABC-and ABA-stacked trilayer graphene [113][114][115], when starting either from a single repulsive on-site Hubbard term or a full set of longer range ab initio cRPA Coulomb repulsion terms [86].Additional results for trilayer ABC-MLG also point to the projected Coulomb interactions being overall ferromagnetic for NNN [92].Additionally, a ferromagnetic J 2 would also help explain the ferrimagnetic state found in ab initio spin-polarized DFT calculations of ABC-MLG [66][67][68][69], as we discuss further in Sec.VI.Taken together, we expect antiferromagnetic J 1 > 0, while a ferromagnetic J 2 < 0.
Notably, an antiferromagnetic J 1 has been proposed as a mechanism for generating chiral d + id superconductivity in doped graphene [83], as well as chiral or nematic d-wave superconductivity in magic-angle twisted bilayer graphene [27,34].However, our results in Fig. 2 show that a similar mechanism cannot exist for ABC-MLG, since no NN pairing is enhanced by the surface state due to its sublattice polarization.This effectively excludes NN or J 1 driven pairing in ABC-MLG.Left are then only ON and NNN pairing.With the expectation that J 2 < 0 and ferromagnetic, we conclude from the phase diagrams in Fig. 6, that the most likely superconducting state in ABC-MLG is the spin-triplet f -wave state.The only other realistic possibility is a phonon-driven spin-singlet s-wave, but due to weak electron-phonon coupling and strong repulsive Coulomb interactions in graphene, we expect that to be a less likely option.This argument for spin-triplet f -wave pairing in ABC-MLG is the fourth important result of this work.

VI. MAGNETISM
In the previous two sections we have investigated the superconducting channels in ABC-MLG and given an estimate of the interaction strengths, respectively, together strongly pointing to a spin-triplet f -wave state.However, the divergent DOS of the topological flat bands also enhance orders in the particle-hole channels, capturing charge and magnetic orders [73].Intense order competitions are therefore likely.In this section, we therefore start by addressing all possible ordering in the particlehole channel in ABC-MLG.
Suggestively, ab initio spin-polarized DFT calculations of ABC-MLG have already found a magnetically ordered electronic structure [58,[66][67][68][69], and a Hubbard model, producing an intrinsic magnetic gap, has been invoked to explain transport experiments of ABC-trilayer graphene [87].Further signatures of magnetic ordering in ABC-MLG have also been found in multiple experiments [62,64,70,72].These results all indicate a strong tendency for magnetic ordering, although we at the same time note that spin-polarized DFT calculations are unable to capture superconducting orders and hence these results cannot resolve any competition between superconductivity and magnetism, which is instead our goal of the next section.
To be able to compare superconductivity with magnetism, we need to treat them on an equal footing.We therefore here first need to solve our tight-binding self-consistent mean-field equations for the particle-hole channels.Ordering in the particle-hole channels originates from the Hartree, or direct, and Fock, or exchange, terms in the mean-field decomposition of Eq. ( 3).Accounting for all these direct and exchange terms out to the NNN range, tabulated for each channel in Table I, we solve the self-consistency equations (4b) and (4c).We find that the magnetic orders (spin-triplet in Table I) generally dominate over the charge orders (spinsinglet in Table I), in agreement with previous results [66][67][68][69].Additionally, we find that the Fock exchange bond magnetization orders mν ij are negligible, compared to the orbital magnetization from the direct terms, except the on-site term where we have mν ii ≡ m ν ii .Based on these observations, we are able to greatly simplify our analysis and hereafter focus exclusively on the orbital magnetization with only the channel coupling strengths Here we have used the fact that for r = 0, the direct and exchange terms contribute additively, but for simplicity we will ignore the factor of 1  2 henceforth.We have also already evoked the expected signs of the interactions in ABC-MLG, i.e., a negative sign has been included in U 0 and J 2 with regards to Table I, as established in the previous section: ON repulsion U 0 > 0, antiferromagnetic NN Heisenberg term J 1 > 0, and ferromagnetic NNN Heisenberg term J 2 < 0.
To proceed, we first analyze the magnetic ordering that arises by considering each individual interaction range separately, just as we did earlier for superconductivity.In Fig. 7(a) we show the magnetic critical temperature T c in both bulk (dashed) and eight layer ABC-MLG (solid), produced by either ON, NN, or NNN interactions.In the case of bulk ABC-stacked graphite, there exist a critical coupling for each interactions below which there is no ordering, very similar to the case of superconductivity in Fig. 2. Numerically, we find that the bulk critical coupling strengths are: Γ 1,c = 0.75t, and Γ z,− 2,c = −J − 2,c = 0.37t.Since bulk ABCstacked graphite is not magnetic, we will henceforth assume that Γ z,− r are all below these values.In contrast, for eight layer ABC-MLG, we find that the flat band surface states enhance the susceptibilities of all interaction ranges, which now all display ordering far below their corresponding bulk critical coupling strengths.We note that this is different from the case of superconductivity in Fig. 2, where NN pairing was not at all affected by the surface states.However, as we scale the couplings strengths by their corresponding bulk critical coupling, i.e.Γz,− r = Γ z,− r /Γ z,− r,c , and use this as the normalized interaction strength to plot T c the inset of Fig. 7(a), we find that NN range magnetism produces lower T c compared to both the ON and NNN range.Just as in the superconducting case, we attribute this relative suppression of NN magnetism to the large sublattice polarization on the surfaces of ABC-MLG.Additionally, we note that ON and NNN magnetism both map on to equivalent mean-field models for periodic systems, despite them originating from very different interactions, but with distinct couplings strengths.
We next consider the type of magnetic order generated by the ON, NN, and NNN interactions.In Fig. 7(b) we plot the sublattice resolved layer profile of the magnetic moments for the different interactions.Note that we only show the moments for the ON and NN ranges, since the ON and NNN produce the same moments for the fixed normalized coupling strength Γz,− r = 0.5.Overall, we find very similar magnetic state for the ON, NN, and NNN ranges: large magnetic moments are localized on the outermost surfaces and the magnetization profile generally follows the LDOS of the topological flat bands in Fig. 1(f) and therefore also the profile of the superconducting order displayed in Fig. 2(d).Moreover, we find that the magnetic moments change signs between the two sublattices, but the sublattices do not have the same moment magnitudes, thus creating an in-plane ferrimagnetic ordering.In-between layers we find an antiferromagnetic ordering since an A site sits next to B-sites of the neighboring layer.This magnetic order, emerging from our unconstrained self-consistent tight-binding model matches the ground state found in spin-polarized DFT [66][67][68][69], both in terms of symmetry and in the qualitative profile of magnetic moments, both of which remain qualitatively unchanged for all interaction strengths below the bulk critical coupling strengths, where instead the whole stack becomes magnetic.
In particularly, the close agreement in spatial magnetic moment structure between all mean-field results can be quantified by plotting the ratio between the surface ferrimagnetic moments, R = |m B1 /m A1 | as a function of the normalized coupling strengths Γz,− r in the inset of Fig. 7(b).Specifically, we find the ratio R to be an approximately linear function that is zero at no interaction Γz,− r = 0 and unity at the bulk critical coupling Γz,− r = 1.This means that there is a gradual, approximately linear, connection between a ferromagnetic state with no magnetic moment on the minority lattice to a completely balanced antiferromagnetic state close to the bulk critical coupling.Inbetween these two extremes, the ABC-MLG surface state is in a ferrimagnetic state.
So far we have treated each interaction range separately and discussed its magnetic order.Finally, we here consider all interactions up to the NNN range jointly.We find that the magnetizations of all three different interaction ranges add linearly and additively, to an excellent approximation.Consequently, we find that the surfaces of constant T c as a function of all three Γz,− r=0,1,2 are all simple simplexes when plotted in a three-dimensional space with the three different interactions ranges as orthogonal axes in Fig. 8(a).In Fig. 8(b) we also show that the same holds true for the surface ratio R. We thus conclude that the surface state of ABC-MLG easily harbors a magnetic surface state, with an in-plane ferrimagnetic structure.This state is driven by any of the interactions U 0 > 0, J 1 > 0, J 2 < 0 and if present together, the resulting magnetization combine in an additive fashion.Peculiarly, we find that the ratio R of magnetic moments on the two sublattice sites is independent of the interaction range and increases approximately linearly with the total interaction strength.Taken together, this conclusion about the magnetic state constitute the fifth significant result of this work.

VII. SUPERCONDUCTIVITY VERSUS MAGNETISM
Having established both superconductivity and magnetism in the previous sections, we in this section address the resulting order competition and how it changes with interaction strengths and also gating.The latter is particularly easy in stacked systems, where a dual gate setup enables simultaneous control of both band filling and band tuning by field gating and displacement field bias [62,91,116,117].Such control has already been used to map the phase diagram of twisted bilayer graphene which includes superconducting domes [3,[18][19][20]84], as well as used to uncover two distinct superconducting phase space regions in trilayer graphene [11], where a displacement field bias was used to produce a partial band flattening as the stack is not thick enough for zeroenergy surface states.We here first continue within the same mean-field framework as in the previous sections, and then also complement with ab initio DFT calculations.For the gating we primarily study a homogenous shift of the chemical potential, but note that a displacement field has similar effects [73].

A. Mean-field analysis
With both superconductivity and magnetism both generated by the J 2 interaction, we start by limiting our analysis to when only J 2 is present, as that will inform the rest of our analysis.In Fig. 9(a) we show the critical temperature of both the spin-triplet f -wave superconducting state, now called T c,sc , and the ferrimagnetic order, now called T c,mag , as a function of J 2 for eight layer ABC-MLG.At charge neutrality, i.e. µ = 0 (solid lines), T c,sc is higher for all |J 2 | ≲ 0.25t.This is in spite of the magnetic state having a lower bulk critical coupling strength of |J − 2,c | ≈ 0.4t compared to the critical coupling strength |J + 2,c | ≈ 0.6t for the superconducting state.In fact, the magnetic state only overtakes the superconducting state for J 2 approaching the bulk critical strength J − 2,c .Thus, we find that the f -wave superconducting state is the dominant order for a wide range of weak to moderate J 2 values.
Electronic gating further offers an experimentally accessible tool to tune ABC-MLG and the competition between magnetism and superconductivity.We here investigate the effects of applying a homogeneous chemical potential µ to the whole stack to illustrate the possibilities.In Fig. 9(a) we also display T c,sc and T c,mag for µ = 0.35 meV (dashed lines).As seen, both T c,sc and T c,mag are reduced for all J 2 .However, the reduction of T c,mag is much larger compared to the reduction in T c,sc .For example, even if finite doping requires a finite interactions strength for ordering, |J 2 | ∼ 0.15t is needed for magnetic ordering at µ = 0.35 meV, which is considerably larger than the |J 2 | ∼ 0.05t needed for superconductivity to first appear at the finite doping.In fact, only at |J 2 | ∼ 0.15t does the magnetic state appear, but now only via a first order transition line marked by 1st in Fig. 9(a) where the free energy of the magnetically ordered state crosses the free energy of the normal state [73].Consequently, a finite µ strongly increases the advantage of the superconducting state, although both T c s are reduced.This advantage of superconductivity over magnetism at finite µ can be attributed to a universal feature of ordering in flat band systems.Within mean-field theory it has been shown that all superconducting, or particleparticle, orders maintain a higher critical temperature following a shift of the Fermi energy away from the flat band or DOS peak, compared to any competing charge or magnetic, i.e. particle-hole, orders [73].The result is that a flat band supports superconducting orders over a much wider range of filling factors and superconducting domes are thus likely to form on the flanks of flat bands.To summarize, the sixth important result of this work is that f -wave superconductivity is the dominant order for weak to moderate J 2 interaction in ABC-MLG and that it becomes further favored over the ferrimagnetic state with gating or doping away from charge neutrality.
So far in this section we have only investigated the individual effect of J 2 .We next consider an arbitrary interaction mix of the expected interactions in ABC-MLG as established in Sec.V: ON repulsion U 0 > 0, antiferromagnetic NN Heisenberg term J 1 > 0, and ferromagnetic NNN Heisenberg term J 2 < 0. From Sec.VI we already know that each of the interactions additively increases T c,mag .In contrast, J 2 is the only interaction that generates the f -wave superconducting state.Because of this fact, we anticipate to find that a generic interaction mix needs a minimal proportion, here defined as γ min of J 2 to maintain the T c,sc > T c,mag .To fully establish this general statement for an arbitrary interaction mix however requires a careful parametrization of the interactions.
To find a suitable way to study any interaction mix beyond just J 2 , we start by noting that the analysis is greatly aided by the observations that the constant T c,mag surfaces in Fig. 8(a) are simple simplexes.To an excellent approximation, T c,mag is therefore given by the additive contributions of U 0 , J 1 , and J 2 .Thus, if we define Γ z,− r (T c,mag ) to be the value of Γ z,− r that generates the critical temperature T c,mag at µ = 0 when Γ z,− r is the only coupling present, then it is also true that the mix of interactions given by {Γ z,− 0 , Γ z,− 1 , Γ z,− 2 } = {αΓ z,− 0 (T c,mag ), βΓ z,− 1 (T c,mag ), γΓ z,− 2 (T c,mag )} also generates T c,mag over the entire simplex spanned α + β + γ = 1, here assuming α, β, γ ≥ 0. As a consequence, we first choose to parametrize any mix of interactions first by the relative proportion of Γ z,− defined by the coefficients α, β, and γ respectively, and study its influence.Second, we also study the total net interaction strength given by the critical temperature T c,mag that this mix generates at µ = 0.
Having established a useful interaction parametrization, we first look at the effect of tuning the relative proportions of the individual interactions, keeping the total net strength constant such that T c,mag = 3 K at µ = 0.In Fig. 9(b) we show the full phase diagram for eight layer ABC-MLG in the α, β, γ parameter space, or equivalently viewed as plotting the phase diagram on a constant T c,mag simplex surfaces such as those in Fig. 8(a).At the top vertex the interaction mix consists solely of  5), but to an unstable equilibrium where the magnetic state free energy is higher than the normal state free energy; instead magnetic state has a first order transition to the normal state (vertical dashed red line, marked by 1st).(b) Ternary plot of the phase diagram on the simplex α + β + γ = 1, (α, β, γ ≥ 0) corresponding to the interaction strengths {Γ z,− 0 , Γ z,− 1 , Γ z,− 2 } = {U0, J1, −J2} = {αΓ z,− 0 (Tc,mag), βΓ z,− 1 (Tc,mag), γΓ z,− 2 (Tc,mag)}, where Γ z,− r (Tc,mag) is the coupling strength that by itself gives Tc,mag = 3 K.Note that Tc,mag is here constant over the entire simplex, as established in Fig. 8. Region I (green) is superconducting, with Tc,sc > Tc,mag both at µ = 0 and finite µ.Region II (red) is magnetic, with Tc,sc < Tc,mag at all µ.Region III (light green) is magnetic or superconducting depending on the chemical potential µ, with Tc,sc < Tc,mag at µ = 0, but for a finite µ, Tc,sc > Tc,mag.(c) Region I boundary γmin (solid) and region III boundary γ * min (dashed), defined in (b) as a function of Tc,mag with regions I-III indicated.Note that solutions with γmin > 1 have no region I in the phase diagrams, since the region I boundary γmin lies outside the phase diagram; the system is therefore magnetic at µ = 0.
J 2 and the superconducting order has the highest critical temperature with T c,sc > T c,mag , as already known from Fig. 9(a).Away from the top vertex, the proportion γ of J 2 in the interaction mix decreases and, since J 2 is the only interaction supporting the superconducting state, there exists a minimal proportion γ min below which the magnetic state instead has the highest critical temperature.Region I spanning γ > γ min therefore represents the parameter space where ABC-MLG is superconducting at charge neutrality.Conversely, for smaller proportions J 2 , in regions II and III ABC-MLG is magnetically ordered at µ = 0.However, because a finite µ give a competitive advantage to superconductivity, as already established in Fig. 9(a), we find that in the intermediate region III, defined as γ min > γ > γ * min ), ABC-MLG is superconducting for at least at some finite µ.
In order to easily establish the position of γ * min in parameter space, we use Ref. [73], which derives that, assuming an isolated flat band at the Fermi level for µ = 0 and that the interactions are independent on µ, superconducting domes exists at finite µ flanking a magnetic, or more generically a particle-hole, order whenever T c,sc > T c,mag /2 at µ = 0. Using this result together with the data in Fig. 9(a) we can position γ * min in Fig. 9(b).We note, however, that the assumptions behind this derivation start to breakdown for large enough T c,mag due to the proximity to the bulk critical coupling strengths in combination with the flat band not being entirely isolated.Thus, γ * min is best viewed as a lower bound, which we estimate should be valid for T c,mag ≲ 10 K. Nonetheless, Fig. 9(b) clearly illustrates that there exists a large extended region III, where J 2 comprises only a modest proportion of the total interaction, yet ABC-MLG is still superconducting and not magnetic.
Next, we show in Fig. 9(c) how the phase diagram in Fig. 9(b) evolves as we tune the total interaction strength, here parameterized by T c,mag at µ = 0.In particular, we show how the phase boundaries γ min and γ * min evolve with T c,mag .We find that γ min and γ * min increases with T c,mag , a behavior that can be derived to the fact that the difference between T c,sc and T c,mag decreases as J 2 increases, as seen in Fig. 9(a).In particular, the region I of Fig. 9(c) remains finite until the two curves cross, T c,sc = T c,mag , in Fig. 9(a), beyond which the systems is magnetic at µ = 0. Still, even as T c,mag increases further there is a considerable Region III for which ABC-MLG is superconducting for some finite doping.Thus, there exists a wide range of coupling strengths for which only a modest proportion of J 2 generates superconducting phase space regions (green and light green).Finally, we note that we have in the above results completely treated the most relevant interactions of Eq. ( 2) out to the NNN range.In terms of superconductivity, we however also know from Fig. 6 that U 2 is an important parameter, although it has little to no effect on the ferrimagnetic state.In Fig. 6 it is evident that the effect of U 2 on the NNN spin-triplet f -wave state is a simple shift of the coupling strength.
Thus, the net effect of incorporating a finite This somewhat shrinks the superconducting regions of the phase diagram.
The final result of incorporating all relevant interaction strengths out to the NNN range is that the superconducting f -wave state surprisingly dominates over the ferrimagnetic state, in fact, as soon as the system exceeds a minimal proportion of J 2 in the total interaction strengths.The superconducting state is further favored by finite doping µ, such that a substantial part of the interaction parameter space results in a superconducting state, as displayed in Figs.9(b) and 9(c).This abundance and resilience of the superconducting state when incorporating the relevant interactions constitutes the seventh important result of this work.

B. Density functional theory analysis
In the previous subsection we established the general phase diagram and competition between the ferrimagnetic and the spin-triplet NNN f -wave states in ABC-MLG within unconstrained mean-field theory for all interaction ranges out to NNN.Here we complement this analysis with ab initio spin-polarized DFT calculations.While DFT cannot capture the superconducting state, it can provide additional information on the magnetic state.Multiple DFT results already exist [66][67][68][69], but we are here primarily concerned with being able to explicitly match with our mean-field theory results, as they offer the direct comparison with the superconducting state.To proceed, we model an eight layer ABC-MLG system with the intralayer (interlayer) carbon distance a CC = 1.42 Å c = 3.347 Å [118] in a slab geometry using QUANTUM ESPRESSO [119,120].We assume in-plane translational invariance and use a 204 × 204 inplane k-point sampling.Additionally, we surround the ABC-MLG stack with 30 Å of padded vacuum.Due to highly varying results in the previous literature [66][67][68][69], we use and compare the results of multiple different exchange-correlation (XC) functionals, including also van der Waals corrections to accurately capture the weak but important interlayer couplings [121].
In Fig. 10 we show the sublattice and layer-resolved magnetic moments calculated for six different choices of XC functionals (see legend in Fig. 10).For all choices, the spin-polarized DFT calculations give a magnetically ordered state at T = 0, with intralayer ferrimagnetism localized to the surfaces with a strong decay of the magnetic moments in to the center of the ABC-MLG slab.This is qualitatively the same magnetic state we find within the mean-field theory modeling in Sec.VI and also consistent with previous DFT results [66][67][68][69].Yet, we also find that the size of the magnetic moments depends strongly on the choice of XC functional.Overall, the sizes of the magnetic moments are small and of the order of ∼ 10 −3 µ B , µ B being the Bohr magneton, but Figure 10.Layer profile of sublattice resolved magnetic moments m (units of Bohr magneton µB) calculated within DFT for different XC functionals.Legend summarizes magnetic moments on A1 and B1 surface sites as well as the ratio R = |mB1/mA1|, sorted from largest to smallest mA1.Here, PBE-SSSP is a PBE XC functional based on the Standard Solid State Pseudopotentials (SSSP) [122,123], while all other XC functionals are built-in in QUANTUM ESPRESSO.Note that vdW-F2 and vdW-b86r are special XC functionals for van der Waals systems.We also included PBE vdW, which includes a non-local van der Waals interaction, in addition to the pure PBE.
there is a factor of 2 difference depending on the choice of XC functional.Notably, while van der Waals interactions are expected to be important in ABC-MLG, we find no notable improvement in the consistency of the results when including or not including van der Waals interactions.With this strong variations in magnetic moments, and the fact that all interaction ranges investigated within our mean-field treatment produces the ferrimagnetic state, we can unfortunately not use these DFT-calculated magnetic moments to get a better handle on the realistic magnitudes of the mean-field interaction parameters U 0 , J 1 , and J 2 .
Despite the large variation in the calculated magnetic surface moments, we can still actually extract useful data.For this, we turn to our mean-field results in Fig. 7(b), which established that the ratio of the surface magnetic moments R = |m B1 /m A1 | to an excellent approximation are independent on the interaction range and instead only depend on normalized channel coupling strengths.We therefore also compute R from our DFT results and we find that, despite the large magnitude variations, R remains remarkably stable across all different XC functionals, ranging from just R = 0.39 to 0.46.This both establishes consistency between the mean-field theory and DFT results and that R is an efficient tool to capture the magnetic structure of the ferrimagnetic state.
We next use the consistency of R to extract the full ternary phase diagram of magnetism and superconductivity.In Fig. 8(b) we established that constant R surfaces are to excellent approximation simple simplexes when plotted against all three ranges of interactions, just as the T c,mag surfaces in Fig. 8(a).In Fig. 9(b) we utilized the latter result and plotted the full ternary phase diagram with magnetism and superconductivity of ABC-MLG for a fixed T c,mag .Here we similarly extract the phase diagram using a constant R simplex, choosing R = 0.425 as an average extracted from the DFT calculations.Note that we could not easily have performed this analysis earlier only within the mean-field theory results as we then did not know the size of R. In Fig. 11, we plot the ternary phase diagram plot of ABC-MLG on the constant R = 0.425 surface.Each point in Fig. 11 thus correspond to an interaction mix, {Γ and the coupling strengths Γ z,− r (R) are defined to produce a magnetic order with R = 0.425 when it is the only coupling in the mix.We find that the constant R ternary phase diagram of Fig. 11 shares its basic features with the constant T c,mag ternary phase diagram of Fig. 9(b).Both phase diagrams contain three different regions.In region I, close to the top vertex, superconducting order dominates both at µ = 0 and for finite µ, since then the interaction mix consist primarily of J 2 which favors the superconducting state.
In contrast, close to the other two other vertices, U 0 or J 1 is the main ingredient of the interaction mix, and the magnetic state dominates in Region II, while in region III a superconducting state is still present at finite µ.The only real difference between the constant T c,mag and constant R phase diagram is that the phase boundaries in Fig. 11 are slanted instead of horizontal as in Fig. 9(b).
The origin of this slant can be found in the normal directions of the constant R and T c,mag simplexes in Fig. 8, which in turn can be traced back to the lower T c,mag produced by the NN ranged interaction compared to the ON and NNN interactions seen in Fig. 7(a), while R is less affected.The existence of a constant R value, largely independent on XC functional, and thus the possibility to extract the phase diagram in Fig. 11 is the eighth major result of this work.It thoroughly establishes that only a moderate portion of J 2 is needed in order to achieve superconductivity in ABC-MLG as guided by DFT calculations, possibly as domes flanking a ferrimagnetic state.

VIII. CONCLUSIONS
In this work we analyze the electron interaction driven ordering possibilities in ABC-stacked multilayer graphene (ABC-MLG), focusing on thick stacks that host topological surface flat bands at charge neutrality.We consider all spin-symmetric short-ranged twosite interactions, out to the next-nearest-neighbor (NNN) range within a fully general mean-field treatment, complemented with general group theory analysis, as well as spin-polarized DFT calculations.For these interactions, the infinite layered ABC-MLG, i.e., bulk ABC-stacked graphite, remains unordered in all ordering channels, due but where instead the simplex Tc,mag = 3 K was chosen.Note that R is here constant over the entire simplex, as established in Fig. 8. Region I (green) is superconducting, with Tc,sc > Tc,mag both at µ = 0 and finite µ.Region II (red) is magnetic, with Tc,sc < Tc,mag at all µ.Region III (light green) is magnetic or superconducting depending on the chemical potential µ, with Tc,sc < Tc,mag at µ = 0, but for a finite µ, Tc,sc > Tc,mag.
to its vanishing DOS at charge neutrality demanding a too large finite critical interaction strength for ordering.In contrast, the topological flat band surface states of finite stacks of ABC-MLG strongly enhance most order possibilities.
For superconducting ordering, we find that all pairing channels hosts a finite critical temperature T c,sc for surface superconductivity also for weak interactions, with the sole exception of pairing in the nearest-neighbor (NN) range, which we attribute to an almost complete sublattice polarization of the surface states preventing NN pairing.We further establish that the superconducting phase diagram is dominated by an isotropic spinsinglet s-wave pairing and unconventional NNN bond spin-triplet f x(x 2 −y 2 ) -wave pairing, which both have a full energy gap.We find that the spin-triplet f -wave state is favored by predominantly ferromagnetic NNN bond Heisenberg interactions (J 2 ), while the s wave is favored by antiferromagnetic such interactions.With intra-sublattice interactions having been identified as ferromagnetic in graphene, i.e.J 2 < 0, we conclude that the fully gapped unconventional spin-triplet f -wave state is the most likely superconducting state on the surface of ABC-MLG.We further note that this f -wave state survives significant on-site repulsion U 0 that hurts the spinsinglet s-wave state, which further promotes tenability of the f wave.Finally, electron-phonon coupling could in principle favor an s-wave state, but due to it being very weak in graphene and in face of strong on-site repulsion, we estimate that the f -wave superconductivity is the most favorable superconducting state in ABC-MLG.These results are in notable contrast with results in bulk ABC-stacked graphite as well as few layer graphene without the ABC-MLG surface state, where instead the chiral d x 2 −y 2 + id xy -or the (p x + ip y )-wave state displays the best ordering tendencies.We establish that this very different bulk versus surface behavior stems from the sublattice polarization of the surface states.
The topological surface flat bands also generically promote particle-hole ordering in ABC-MLG.By treating all charge and magnetic ordering possibilities on an equal footing, we find that ABC-MLG has a ferrimagnetic order on the surface for all relevant spin-symmetric short-ranged interactions with mean-field theory, with results fully agreeing with spin-polarized DFT calculations.Thus, while bulk ABC-graphite is not magnetically ordered, ABC-MLG exhibits a finite T c,mag also for weak interactions.We further find that the ferrimagnetic order interpolates between complete sublattice ferromagnetism for weak coupling to a complete antiferromagnetic state close to the critical bulk coupling strengths, quantified by the ratio between the surface sublattice magnetic moments R = |m B1 /m A1 | ranging from zero and one.Additionally, we find that differently ranged interactions contribute additively to the ferrimagnetic order and therefore the constant T c,mag surface is a simple simplex in the {U 0 , J 1 , −J 2 } interaction space, which encodes the relevant interactions: repulsive on-site Coulomb U 0 > 0, antiferromagnetic J 1 > 0, and ferromagnetic J 2 < 0.
Having both f -wave superconductivity and ferrimagnetic ordering as possibilities in ABC-MLG, we finally establish their competition and the resulting phase diagram as a function of interaction strengths.Both ferrimagnetism and f -wave superconductivity are driven by a ferromagnetic J 2 , but we find that the susceptibility for the superconducting f -wave state is larger, resulting in T c,sc > T c,mag over a large range of weak to moderate J 2 interaction strengths.Consequently, we find that the f -wave state is the dominant order in ABC-MLG, as long as the proportion of J 2 in the full mix of interaction parameters is sufficiently large.Moreover, we find that superconductivity is further favored by gating, resulting in a homogeneous doping or an electric field gradient perpendicular to the ABC-MLG stack, which suppresses the magnetic order much more than the superconducting order.Thus, f -wave superconductivity generally appears upon doping, even if the J 2 proportion is relatively small.By using spin-polarized DFT calculations we are able to further specify the phase diagram.We find that the sizes of the magnetic moments are highly sensitive to the choice of exchange-correlation functional used for the DFT calculations.Still, the surface sublattice moment ratio R is remarkably stable, which we utilize to construct an improved phase diagram, with a similarly notable superconducting portion.Overall, our results establish that spin-triplet f -wave superconductivity is possible from purely repulsive interactions on the surface of ABC-MLG and feasible over a large region of the phase diagram in competition with a ferrimagnetic state.The f -wave state further gains a competitive advantage when doping slightly away from from charge neutrality.We finally note that due to the special properties of the surface flat band in ABC-MLG, these results are notably different from doped or twisted few layer graphene, which clearly illustrates the fascinating possibilities and variations of electronic ordering in graphene systems.
But we also that Tr [ρ 0 a i a j ] = 0 since the order parameters vanish in the noninteracting phase.Thus, only the second term in the expansion in Eq. (A1) contributes to the first order perturbation and thus the LGE.
To proceed, we let Ŝ (k) be the unitary matrix that diagonalizes the noninteracting Hamiltonian, given in real space by Eq. ( 1) in the main text, and ξī (k) the corresponding energy bands, where ī labels the bands.Then, the real space operators a jLα used in the main text can be written in terms of band operators c k jα as a jLα = k j Ŝj L,j (k) e ik•rj c k jα , (A3) where r j is position vector of the site j.Using Eq. (A1) and Eq.(A3) in Eq. (A2), and putting the result into the self-consistency equations, Eq. ( 4) in the main text, we obtain the LGE as where the order parameter D η,+ ij = ∆ η ij D η,− ij = m η ij , mη ij encodes superconducting (particle-hole) order.The meanfield decomposed coupling strengths Γ η,± ij are listed in Table I, while the susceptibilities χ ± ī, j are given in Eq. ( 6), both in the main text.Finally, the factor S contains symmetry-related information and is given by S ī, j,η,+ L,ij;L,sr (k) = The exponential factors in S η,− differentiates between bond-independent direct (Hartree) particle-hole order parameter, m η ij , and bond-dependent exchange (Fock) particle-hole order parameter, mη ij .We also note that S η,− is explicitly spin independent, implying the sign of the coupling differentiates between particle-hole orders.
From Eq. (A4) we see that each order parameter on the left-hand side of the equation is a linear combination of the order parameters on the right-hand side and each term in the linear combination is weighted by Γ η,± S η,± χ η,± .Thus, the LGE in Eq. (A4) is a system of linear equations and by collating the order parameters on both sides as column vectors, D η,± , the LGE can be written in the compact form D η,± = M η,± D η,± , which is presented in Eq. ( 5) in the main text.

Figure 1 .
Figure 1.(a) Top view of the atomic positions for three layers of graphene (blue,black, green) with ABC stacking order.Red dashed region shows the unit cell of length a.(b) Side view of atomic positions for eight layers with bonds labeled by their tight-binding hopping amplitudes of Eq. (1).(c-f) Normal state electronic structure of eight layer ABC-MLG: (c) Contour plot of lowest positive energy band ξ( ⃗ k).Dashed brown hexagon marks the first Brillouin zone with high symmetry points Γ, K, K ′ , and M .(d) Energy band plot ξ = ξ (kx = 0, kya) along the dashed vertical line in (c), passing through the K ′ -point.Band line color shows probability density on surface site A1 (or equivalently B8), marked by X in (b).(e) Local density of states (LDOS) on surface sites A1 and B1 and third layer sites A3 and B3.(f) Integrated LDOS, P = ξ LDOS (ξ) resolved by sublattice and layer for narrow low-energy window ξ ∈ [EF − 3 meV, EF + 3 meV] around EF = 0. (g) Surface flat band radius, rK as a function of number of layers NL in a ABC-MLG stack.The radius rK , illustrated by a red line in (d), is estimated as the momentum where the flat band exits the narrow energy window used in (f).
) -2(f) we plot the layer and sublattice resolved magnitude of the normalized superconducting order parameter going from weak [Fig.2(d)] and intermediate [Fig.2(e)] to strong coupling [Fig.2(f)] strengths, where we use spatially normalized magnitudes as

2 T
ky pyTableII.Symmetry allowed superconducting states for pairing with range r = 0, 1, 2 encoding on-site (ON), nearest-neighbor (NN), and next-nearest-neighbor (NNN) pairing, respectively, divided into spin-singlet η = 0 and spin-triplet η = ν, with ν = 1, 2, 3. Results are based on D 6h point group for the normal state Hamiltonian (1).Columns represent range r, order parameter vector in real space ⃗ ∆r [see Fig.2(a) and with notation h ′ = −h], func.meaning the expansion in reciprocal space of the basis function around the Γ point up to cubic order, irrep.meaning the irreducible representation, char.meaning the nature of the orbital part of the superconducting order parameter, and their prevalence in bulk (graphite) and surf.(ABC-MLG).For the prevalence, ✓ implies dominant state at Tc, • implies degenerate with dominant state at µ = 0 but sub-dominant at µ ̸ = 0 at Tc, ⋆ implies sub-dominant at Tc, and − implies the state does not occur or its eigenvalue at Tc is less than 10% of that of the dominant state.Note that E1u and E2g irreps are two-fold degenerate.
(b) pass through the M -points of the Brillouin zone, which therefore do not intersect with the flat bands residing (a) (b)

Figure 9 .
Figure 9. Order competition between ferrimagnetism and f -wave superconductivity in eight layer ABC-MLG.(a) Critical temperatures Tc,mag of magnetism (Mag) and Tc,sc of triplet f -wave superconductivity (SC) as a function of the effective ferromagnetic NNN Heisenberg interaction J2 at charge neutrality µ = 0 (solid line) and µ = 0.35 meV (dashed line).Grey line are also solution to the LGE, Eq. (5), but to an unstable equilibrium where the magnetic state free energy is higher than the normal state free energy; instead magnetic state has a first order transition to the normal state (vertical dashed red line, marked by 1st).(b) Ternary plot of the phase diagram on the simplex α + β + γ = 1, (α, β, γ ≥ 0) corresponding to the interaction strengths {Γ z,− 0 , Γ z,− 1 , Γ z,− 2 } = {U0, J1, −J2} = {αΓ z,− 0 (Tc,mag), βΓ z,− 1 (Tc,mag), γΓ z,− 2 (Tc,mag)}, where Γ z,− r (Tc,mag) is the coupling strength that by itself gives Tc,mag = 3 K.Note that Tc,mag is here constant over the entire simplex, as established in Fig.8.Region I (green) is superconducting, with Tc,sc > Tc,mag both at µ = 0 and finite µ.Region II (red) is magnetic, with Tc,sc < Tc,mag at all µ.Region III (light green) is magnetic or superconducting depending on the chemical potential µ, with Tc,sc < Tc,mag at µ = 0, but for a finite µ, Tc,sc > Tc,mag.(c) Region I boundary γmin (solid) and region III boundary γ *