Competition of electron-phonon mediated superconductivity and Stoner magnetism on a flat band

The effective attractive interaction between electrons, mediated by electron-phonon coupling, is a well-established mechanism of conventional superconductivity. In metals exhibiting a Fermi surface, the critical temperature of superconductivity is exponentially smaller than the characteristic phonon energy. Therefore such superconductors are found only at temperatures below a few Kelvin. Systems with flat energy bands have been suggested to cure the problem and provide a route to room-temperature superconductivity, but previous studies are limited to only BCS models with an effective attractive interaction. Here we generalize Eliashberg's theory of strong-coupling superconductivity to systems with flat bands and relate the mean-field critical temperature to the microscopic parameters describing electron-phonon and electron-electron interaction. We also analyze the strong-coupling corrections to the BCS results, and construct the phase diagram exhibiting superconductivity and magnetic phases on an equal footing. Our results are especially relevant for novel quantum materials where electronic dispersion and interaction strength are controllable.


I. INTRODUCTION
The overarching idea in quantum materials is to design the electronic (or optical, magnetic, etc.) properties of materials to perform the desired functionality 1 . This goal is aided by generic models and concepts, such as specific lattice models that lead to certain topological phases. Often the studied models and the resulting topological phases for electronic systems are noninteracting, and do not include the possibility of spontaneous symmetry breaking. However, such noninteracting models are platforms for exotic electron dispersions that provide a basis for studying symmetry broken interacting phases. In particular, certain models support approximate flat bands 2-10 , and here we consider microscopic mechanisms for symmetry breaking phases in such systems.
We analyze the interplay of electron-phonon 11 and (screened) electron-electron interaction in providing means for a symmetry-broken phase transition, thereby coupling together works on flat band superconductivity 2,7,10,12 with those on flat band (Stoner) magnetism 9,[13][14][15][16][17] . In both cases the resulting mean-field critical temperature is linearly proportional to the coupling constant 18 , thus allowing for a very high critical temperature. The two types of interaction mechanisms work in opposite directions, and in the case of weak interactions in a symmetric way. However, upon increasing the coupling strength the retarded nature of the electron-phonon interaction shows upas opposed to the instantaneous electron-electron interaction -breaking the symmetry between the two. In particular, we generalize the Eliashberg's strong-coupling theory of superconductivity 19 , usually formulated for systems with a Fermi surface, for flat bands. As a result, we describe the dimensionless BCS attractive interaction 20 in terms of the electron-phonon coupling and the char- . T CE C is the temperature at which the TC 's of magnetic and superconducting order coincide. In the striped region these phases can form metastable domains inside the sample. This diagram is for N → ∞. For finite N the overlap region between the phases is smaller. acteristic phonon frequency [Eq. (8)]. In addition, we provide the generalization of the well-known McMillan formula of strong-coupling superconductivity (for Fermi surface systems) 21 to the case with flat bands in Eq. (14).
In addition to superconductivity, we consider flat-band Stoner magnetism. Because of the retarded nature of the electron-phonon interaction, the combined interaction can simultaneously have attractive and repulsive components, and thus the system can be unstable with respect to both singlet superconductivity and magnetism arXiv:1801.01794v5 [cond-mat.supr-con] 30 Nov 2018 (see a generic strong-coupling phase diagram in Fig. 1). Often one of the phases still dominates and suppresses the other, but we find that when the critical temperatures of the phases are similar, both phases are local minima of the free energy at low temperatures. We find that their bulk coexistence and the resulting odd-frequency triplet superconducting order 22,23 is only realized as an unstable solution. On the other hand, these phases can form metastable domains inside the sample, and therefore odd-frequency triplet order parameter can appear at the domain walls.
The structure of this paper is as follows. In Sec. II we introduce the model of surface bands with electronphonon and Coulomb interactions. In Sec. III we formulate the Eliashberg model extension for the surface bands, describe all possible ordered states that can appear within this model and calculate the critical temperatures of the superconducting and antiferromagnetic states. We study the competition and possible coexistence of these two types of ordering in Sec. IV. Conclusions are given in Sec. V.

II. MODEL
As a low-energy model for the flat band, we assume two sublattices coupled through an electronic Hamiltonian 3 where an integer N parametrizes the flatness of the dispersion, and ε 0 is the energy at p = p F B . The model is electron-hole symmetric and the two energy bands have the dispersions ±ε p . For large N , the states with low momenta, |p| < p F B , are almost at zero energy and the density of states is very high. The states with momenta larger than p F B do not contribute much to the momentum integrals due to their low density of states. Therefore, the results for large N do not depend much on the momentum cutoff, as long as it is larger than p F B . In our model we take the cutoff to infinity and consider only the cases N > 2. This is in contrast to models with isolated flat bands extending throughout the Brillouin zone. The effects discussed below are mostly applicable also to such models (provided they have the type of sublattice degree of freedom discussed below) for large N , as long as p F B is taken as the size of the Brillouin zone. Equation (1) is approximately realized for the surface states of N -layer rhombohedrally stacked graphite. In that system the surface states delocalize into the bulk at the edges of the flat band and this gives a momentum-dependent correction in the low-energy Hamiltonian 12,24 . In the case of N → ∞ the delocalization of the surface states to the bulk leads to strong amplitude mode fluctuations invalidating the mean-field theory 24 . Therefore, the theory considered in this paper is applicable to rhombohedral graphite only in the case where N is not too large.
We model the electron-electron interaction as a repulsive on-site Hubbard interaction 25 with energy U . The magnitude of U depends on the microscopic details of the system and its environment. The coupling between electrons and phonons, with strength g, creates an effective attraction between the electrons and makes the system susceptible to superconductivity 19 . We mostly consider Einstein phonons with constant energy ω q = ω E and discuss generalizations in the Supplementary Information.
The total Hamiltonian incorporating these effects is where N is the number of lattice points in the system and Ψ † pσ = (ψ † pσA , ψ † pσB ) is a pseudospinor in sublattice space. We assume that the low-energy states on the two sublattices ρ = A/B are spatially separated (e.g. localized on the two surfaces in rhombohedral graphite), so that neither the electron-electron interactions nor the phonons couple them. The only coupling between the sublattices comes from the off-diagonal dispersion relation. In the Supplementary Information we also show that the flat band phenomenology applies to linear, graphene-like dispersion with an electronic Hamiltonian with an energy cutoff ε c and Fermi velocity v F , provided the interaction energy scales are large compared to ε c . Hence, our results may also apply as an effective model for twisted bilayer graphene close to its "magic" angle. 26 In the theory of electron-phonon superconductivity of metals, the neglect of higher order diagrams in the perturbation theory is typically justified with the help of the Migdal theorem 27 . In that case, the expansion parameter gets an additional factor of ω E /E F , where E F is the Fermi energy. Because of the Migdal theorem, the theory of superconductivity for metals is not strictly limited to weak coupling with respect to the interaction parameter.
In the flat band, however, the chemical potential is located at the bottom of the band and there is no Fermi energy with which to compare the Debye energy. Migdal's theorem cannot be used in this case. In the intermediate case of narrow electronic bands, corrections in the higher orders of the adiabatic parameter ω E /E F have been studied in Refs. 28-31 and the Eliashberg theory has been found also to be in agreement with Monte Carlo results in the weak coupling regime when ω E /E F = 1 in Ref. 32. We find that the diagrams beyond the mean-field approximation do not influence the self-energies significantly if the effective pairing constant introduced below in Eq. (8) is small λ 1 and ω E , u ε 0 . Moreover, although the mean-field theory is applied beyond its formal limits of validity in the strong-coupling regime, this theory captures the interesting possibility that the retarded nature of the electron-phonon interaction can lead to the presence of attractive and repulsive components at the same time. As a result, the system can be simultaneously unstable with respect to the appearance of both singlet superconductivity and magnetism as discussed in Sec. IV.

III. ORDERED STATES
The Hamiltonian (2) allows for a number of spontaneous symmetry breaking phases. We restrict our study to spatially homogeneous phases. Therefore, the order parameter can appear in the spin, sublattice (pseudospin) and electron-hole (Nambu) spaces. The general selfenergy is where τ i , σ j and ρ k are the Pauli matrices in electronhole, spin and sublattice spaces, respectively. We characterize the different components Σ ijk and determine their values within the self-consistent Hartree-Fock model. This reduces to solving a set of non-linear integral equations, known as Eliashberg equations in the context of conventional superconductors.
To explore the possible phases of the system, we first assume that the U (1)-gauge symmetry is broken, but the SU (2)-spin rotation symmetry is not. After fixing the overall phase of the superconducting order parameter, we are left with the self-energy Σ 000 (iω n ) and three degrees of freedom for the superconducting singlet order parameter: the magnitudes of the order parameter on the sublattices ∆ A and ∆ B and the relative phase θ. Choosing θ = 0 leads to a gapped quasiparticle dispersion (Fig. 2c), whereas θ = π would imply a gapless dispersion (Fig. 2b). Thus, in the case of an instantaneous interaction the total energy is minimized when θ = 0 and ∆ A = ∆ B . Generalizing the above to the frequency dependent interactions, we choose the singlet to be proportional to the τ 2 σ 2 ρ 0component, whose magnitude and the functional form are obtained from the self-consistency equation. The selfenergy for the fermionic Matsubara frequency ω n is where Σ ω n = (1 − Z n )ω n is the frequency renormalization by the retarded interaction 19 . To simplify the equations, we define renormalized frequenciesω n = Z n ω n . We use the symbol φ n for the 'bare' singlet order parameter and ∆ for the maximum value of the renormalized singlet order parameter ∆ n ≡ φ n /Z n related to the energy gap.
When SU (2)-spin rotation symmetry is broken but U (1)-gauge symmetry is not, the self-energies describe the frequency renormalization and the magnetization. After fixing the direction of the magnetization on one sublattice, the relevant degrees of freedom are reduced to three similarly as in the superconducting case. These can be chosen as the magnitudes of the magnetizations in the two sublattices h A , h B and the relative angle ϕ between their directions. The quasiparticle dispersion in the magnetic case is the same as in the superconducting case if we identify ∆ A,B = h A,B and θ = π−ϕ (see Fig. 2). In this case, the relative angle ϕ = 0 leads to a gapless quasiparticle dispersion (Fig. 2b), and ϕ = π to a gapped dispersion (Fig. 2c). Thus, the energy minimum is obtained with h A = h B and ϕ = π. The stable magnetization is hence antiferromagnetic, with opposite magnetizations on the two sublattices, so that the selfenergy is where h n is the frequency dependent exchange field. This result agrees with DFT studies on rhombohedral graphite 33 and similar magnetization structure has been predicted also in the case of flat bands appearing at the zig-zag edges of graphene nanoribbons [34][35][36] . We also note that the AFM state is insulating (see Fig. 2c). If the noninteracting dispersion is completely flat at zero energy, the sublattices are uncoupled and the antiferromagnetic state is degenerate with the ferromagnetic ϕ = 0 state.
By calculating the Hartree-Fock self-energies, we find the self-consistency equations, from which we can determine the values of the self-energy terms. For the superconducting (SC) self-energy (4), they are where the interaction kernel is The functional form of the interaction kernel is determined by the phonon propagator from which it is derived. The width in frequency space is determined by the characteristic phonon frequency, which in this case is the Einstein frequency ω E . The effective interaction constants in the flat band are where Ω FB and Ω BZ are the momentum-space areas of the flat band and of the first Brillouin zone, respectively.
For an antiferromagnet (AFM) with self-energy (5), Quasiparticle dispersions E(p) for different kinds of symmetry breakings with N = 5. a) In the non-interacting case, the spin-bands are degenerate with E(p) = ±ε(p). b) For the ferromagnetic (FM) or the superconducting (SC) phase with a θ = π phase shift between the sublattices, one quasiparticle band is shifted up, and the other down in energy.
In this case, no energy gap is opened. c) For the antiferromagnetic (AFM) or the SC phase with θ = 0 an energy gap is opened and quasiparticle bands are doubly degenerate.
the self-consistency equations are Superconductivity and magnetism are thus symmetric with each other also on the level of the self-consistency equations, but with the roles of u and λ nm switched. Tovmasyan et al. have shown that this duality is also broken by taking into account higher order terms in the perturbation theory 37 .
To solve the self-consistency equations (6-10), we truncate the Matsubara sums with a cutoff ω C ∼ 10ω E . This causes no numerical error if we use the pseudopotential trick and simultaneously replace u with an effective value u * , which depends on the cutoff 38 . For superconductivity (magnetism) cutting off high-energy scatterings is compensated by a reduction (increase) in the low-energy effective interaction.
After the pseudopotential trick, the solutions are found by a fixed-point iteration. The iteration is continued until all of the components have converged. The fixed-point method only finds the stable solutions, to find the unsta- ble solutions we used a solver based on Newton's method.
For weak coupling λ 1 the frequency dependence of λ nm can be disregarded and we can approximate Z ≈ 1 and ∆ ≈ φ. Assuming λω E > u, the superconducting gap at T = 0 and the critical temperature are T sc These results are valid for N > 2 as the momentum integrals diverge without a cutoff for N ≤ 2. Note that the T = 0 limit can thus be taken before the flat-band limit of large N . Analogous results have been obtained before within the BCS model in Ref. 12. For large N , ∆ 0 is linear in the coupling and its magnitude is proportional to the phonon energy scale. Hence the associated critical temperature can be very large. Relabeling ∆ 0 → h 0 and λ ↔ũ, we find similar equations for magnetism. Here h 0 is the magnetic order parameter at T = 0. At strong coupling, the retardation matters and the results for magnetism and superconductivity diverge from each other. For superconductivity, we can improve on the weak coupling result by including some of the corrections from the Eliashberg theory when N →∞. We still neglect the full frequency dependence, but include the electron mass renormalization as a static factor Z 0 = 1+λ. The order parameter at zero temperature becomes In metals with a Fermi surface 39 , the electron-phonon interaction renormalizes the pairing potential with the factor of 1+λ instead of 1+2λ as in Eq. (13). Thus, for weak coupling, the electron-phonon renormalization is more effective in the flat band than in the usual metals. This difference is more pronounced at strong coupling, as we see next.
By linearizing Eqs. (6) and (7) with respect to φ, we can solve for the critical temperature, see Fig. 3a. We find that when N → ∞, the critical temperature scales as T sc C ∝ λ 0.2 ω E for large λ. In metals 39 the asymptotic scaling goes as T sc C ∝ λ 1/2 ω E . When u = 0, there is a critical point λ C such that for λ < λ C there is no superconducting transition at any temperature. For small u/ω E , λ C is linearly proportional to the Coulomb interaction. For large u, λ C increases sublinearly, see Fig. 3b.
An approximate numerical equation for T sc C is This is a flat band analogue of the McMillan equation 21 which for the the conventional superconductors incorporates the Eliashberg and Coulomb corrections to T sc C . The u 2 -term in the numerator accounts for the retardation correction to λ C as in Fig. 3b. The form of the denominator is chosen to show the λ 0.2 power law behaviour for large λ. The factor 2.6 is obtained by a fit in the region λ < 1 for u = 0. The fit is shown as the dashed line in Fig. 3a.
The ratio ∆ 0 /T sc C is not constant, but depends on both N and λ. For N → ∞, the ratio has the value 2 for weak coupling and increases as λ increases. For λ=1 the ratio is 2.56. For the critical temperature at finite N , see Fig. 4.
The phenomenology of the magnetism can be understood as follows. According to the Stoner criterion, the magnetization is related to the competition between the exchange energy gain and the kinetic energy penalty from moving electrons from one spin band to another. For a flat band with N → ∞, there is no kinetic energy penalty, and at zero temperature with λ = 0 even a small exchange interaction leads to a complete magnetization of the flat band. In the presence of the electron-phonon interaction the competition is between the exchange energy gain and the electron-phonon energy penalty which coincide at u = u C . If we can neglect the retardation, the total interaction in Eq. (9) is u − λω E . The flat band is completely magnetized when u > u C ≈ λω E . Due to retardation, for large λ the critical point is reduced from the linear estimate, see Fig. 3d.
Above, we have discussed the superconducting order parameter φ. The other important property of the superconducting state is the existence of a supercurrent. In the flat band the electronic group velocity vanishes and it is not immediately clear that there can be a finite supercurrent. However, the flat band surface states of superconducting rhombohedral graphite do support a finite supercurrent 40 and similarly it is known that quantum hall pseudospin ferromagnets can support a finite pseudospin supercurrent 16 . More generally, Peotta and Törmä 7 have shown that for a topological flat band there is an additional geometric contribution to the superfluid weight so that the critical current is finite. As we have not fixed the underlying topology in our model, it can be applied to topologically non-trivial flat bands.
As one can see the Eliashberg model describes the nucleation of both the magnetic and superconducting phases which can have rather close critical temperatures as shown in Fig. 3. In the next section we consider the non-linear problem by calculating the entire phase diagram of the ordered states to study the competition and the possible coexistence between the superconductivity and antiferromagnetism.

IV. COMPETITION BETWEEN THE PHASES.
If the electron-phonon interaction is approximated as instantaneous, we can sum the two interactions together and have either a total interaction which makes the normal state unstable to the superconducting transition (λω E −u > 0) or to the magnetic transition (λω E −u < 0), but not to both at once. On the other hand, if the electron-phonon interaction is retarded, the situation is different as the total interaction can be attractive for low frequencies, but repulsive for high frequencies. There is then a parameter range in which both phases are local minima of the free energy. This occurs when λ is large enough to overcome the suppressing effect of u in the case of superconductivity (λ > λ C (u) in Fig. 3b), but at the same time u is large enough to overcome the suppressing effect of λ and create a magnetic instability (u > u C (λ) in Fig. 3d).
We determine the phase diagram of the system by calculating the phase with a higher critical temperature as a function of u and λ (Fig. 5). The phase diagram is almost symmetric with respect SC and AFM phases except that the lack of retardation in electron-electron repulsion favors the AFM phase for strong coupling.
Even if there is a parameter region in T , u and λ where both SC and AFM self-consistency equations have a finite solution, it does not mean that both phases are necessarily simultaneously present. To determine the stability, we construct the coupled self-consistency equations in the case when both order parameters are non-zero and interact with each other 41 . By linearizing the coupled self-consistency equation with respect to SC, and solving the AFM part fully, the stability of the AFM phase with respect to SC transition can be determined, and vice versa. Fig. 1 shows the region in λ-T space with fixed u, where the two phases are stable. The figure shows that in the region where SC is dominant, the AFM phase is unstable near the expected second order transition (the solid line between the magnetic and paramagnetic phases), but becomes a local minimum of free energy at lower temperatures. The same happens for superconductivity when the AFM phase dominates. The transition between SC and AFM phases is of the first order.
When discussing superconductivity in the presence of an exchange field (either induced or spontaneous), we have an additional ingredient in the self-energy, namely the superconducting triplet order parameter 22,42 , which has been discussed in the context of the Eliashberg model in Ref. 43. The triplet is spatially isotropic, and in order to satisfy the fermionic antisymmetry, it has to be odd in frequency. It is generated in the self-energy only when there is an odd-frequency component in the interaction. In the retarded interaction, this is always satisfied. When calculating the stability of AFM with respect to SC, the triplet appears in the linear order. It hence modifies the boundaries of the region where both AFM and SC are stable. We have taken this effect into account in Fig. 1.
Besides the competition between AFM and SC phases, we need to consider the possibility of a coexistence phase in the dashed region of Fig. 1, where both phases can show up alone. We indeed have numerically found such a coexistence solution, but tests based on fixed-point iteration revealed it to be unstable at every temperature that we checked. This finding is in accordance with a simplified model where both interaction channels are instantaneous and independent of each other 41 .
However, the fact that the two phases are simultaneously local minima of the free energy suggests that this system could have domains of antiferromagnetic order coexisting with superconducting domains. Such domains would be separated by a domain wall mixing the two kind of phases and inducing odd-frequency triplet pairing, as schematically illustrated in Fig. 6. In addition to providing a mechanism for the appearance of odd-frequency triplet pairing, the domain walls can support interesting excitations. In particular, it is known that flat band ferromagnets can support interesting topological and domain wall excitations in the form of different kinds of spin textures 16,44 , and various combinations of spin textures and superconductivity may lead to the appearance of Majorana zero modes [45][46][47][48][49] . Also, alternatively to the intrinsic domain structure generation, the ferromagnetic superconductors can support different types of nonuniform magnetic order and spontaneous vortex states [50][51][52] . A detailed analysis of different possibilities goes beyond the scope of this paper.

V. CONCLUSIONS
We have proposed a simplified model of a flat band system with a retarded electron-phonon interaction and a repulsive Hubbard interaction. For this model, we have determined the self-consistency equations in the Hartree-Fock approximation and all the possible homogeneous phases. Antiferromagnetism and superconductivity are essentially symmetric in this system, with the only difference coming from the retardation of the electron-phonon interaction. For large λ, the retardation suppresses the increase in ∆ more effectively in a flat band than in metals with a Fermi surface. We find that the retardation also creates a situation in which both phases are separately local minima of the free energy suggesting a possibility of coexisting antiferromagnetic and superconducting domains inside the sample.
Our results indicate how flat-band superconductivity can be generated from electron-phonon interaction, and provides means to estimate the mean-field critical temperature when the details of the electron-phonon coupling and the screened interaction are known. The superfluid transition in low-dimensional systems occurs in the form of a Berezinskii-Kosterlitz-Thouless (BKT) transition at a temperature that is lower than the mean-field transition temperature. That the latter is non-zero is ensured by the possibility of having a non-vanishing supercurrent (see, for example, Refs. 7, 10, and 40) in a flat-band superconductor. Our results are of relevance in designing novel types of quantum materials for the interplay of superconducting and magnetic order, and the search of systems exhibiting exotic superconductivity with a very high critical temperature, up to room temperature. They may also shed light on recent evidence of high-temperature superconductivity in graphite interfaces 53 .
Our results could also explain some of the phenomena associated with the recent experiments on bilayer graphene 26,54 . (For a more microscopic description of that case within the BCS model, see Refs. 55 and 56.) In the experiment, the twist angle between two superimposed graphene layers is chosen to a certain magic angle, so that the two Dirac cones in the graphene layers hybridize, forming a pair of flat bands. Our model can be adjusted to describe this situation with small changes, see the Supplementary Information for details. When the chemical potential was tuned to lower of these bands, the system became an insulator. From our point of view, this could be the insulating AFM state we describe. When the chemical potential is tuned slightly off from the flat band, a superconducting dome in T − µ phase diagram was observed on both sides. These domes can be the s-wave SC phases we describe here. The competition between the particle-hole (AFM) and the particle-particle (SC) channels in the presence of the chemical potential has been considered by Löthman and Black-Schaffer in Ref. 8, and for a range of parameters, they reproduce a similar phase diagram near the flat band, with the AFM state at the level of the flat band and two superconducting domes with doping away from the flat band (see Fig  2b in the paper). In the experiments, SC domes are only observed on the hole-doped side. The electron-doped side exhibits only insulating behavior near the flat band. One possible explanation is the difference in screening, which changes the relative magnitude of repulsive and the attractive interactions, so that the AFM state covers the SC domes completely. However, we leave the detailed treatment of the effects of doping and screening (both intrinsic and that provided by the environment) for further work.
The Hamiltonian studied in the main text is derived from a tight binding model, in which the size of the system is naturally characterized by the number of lattice sites N . However, in approximating the sum over all momenta, it is more natural to use the area A of the system. The ratio A/N is the area of the real space unit cell A c , which in turn is inversely proportional to the area of the first Brillouin zone Ω BZ . For an infinite system, the momentum sum can be written in the following form.
where Ω FB ≡ πp 2 FB is the area of the flat band and p c is the momentum cutoff. We define a shorthand for the sum/integral over momenta and Matsubara frequencies (S2) We find that the effective interactions on the flat band are characterized by the constants is the phonon propagator. In other words, interactions are proportional to the ratio between the area of the flat band and that of the first Brillouin zone.

Appendix B: Self-energy components
In total, there are 4 3 = 64 combinations of Pauli matrices in spin, Nambu and sublattice spaces. The ones off-diagonal in sublattice space are not possible (see below), as the interactions in the model do not couple the two sublattices. This reduces the number by a factor of 2. We are left with 16 components symmetric and 16 components antisymmetric in sublattice index ρ. Of these 32 components, the 16 components diagonal in Nambu space are associated with non-superconducting properties. The τ 0 σ 0 ρ 0 -component renormalizes the frequencies in the propagator 39 . It vanishes if the interaction is instantaneous and is always present if the interaction has a nontrivial frequency structure. The τ 3 σ 0 ρ 0 -component on the other hand renormalizes the chemical potential, and is usually induced by finite temperature or interaction effects. In this model it vanishes because of the electron-hole symmetry of the model at half-filling. There could in principle also be a term proportional to τ 0 σ 0 ρ 3 . Its effect would be to renormalize frequencies antisymmetrically in the sublattices. However, it is not induced by any of the other terms, so the only way to have it would be by spontaneous symmetry breaking. It should be odd in frequency, and this makes it vanish in the BCS limit. Even in Eliashberg theory, it is unlikely, as it is supported only by the odd-frequency part of the interaction. In this text, we do not consider this and other antisymmetric frequency renormalization components any further.
The remaining 12 components diagonal in Nambu space are due to magnetism. The magnetization direction on one sublattice is described with the three components σ i τ 3 . We can parametrize the six degrees of freedom associated with magnetization with the overall magnitude and the relative magnitude of the order parameter, two angles for the overall magnetization direction, and two angles for the relative direction. The six other terms of the form σ i τ 0 ρ j are spin-antisymmetric frequency renormalization components.
The 16 components off-diagonal in Nambu space are associated with superconductivity. Four of these are associated with the singlet and its phase and the sublattice: overall phase, relative phase between the sublattices, overall magnitude of the order parameter and the relative magnitude between the sublattices. The remaining 12 off-diagonal Nambu components describe the three components of the triplet and its phase and sublattice degrees of freedom. Because both the spatial and the spin parts of the triplet are symmetric, it must be odd in frequency to preserve the fermionic antisymmetry 22 . Such components are supported by the odd-frequency part of the electron-phonon interaction 23 .
Appendix C: Hartree-Fock self-energies Starting from the Hamiltonian and using the above definitions for the interactions, we can write the self-energies in the Hartree-Fock approximation as where P ρ is the projection operator to sublattice ρ, ρ ∈ A, B. It is immediately clear from the above expressions that the self-energy cannot have terms with A−B mixing, as the projection operators force it to be diagonal in the sublattice space. We assume that the radius of the flat band is much smaller than the maximum phonon momentum, so that the phonon cutoff does not need to be enforced in the momentum sum. With a contact interaction, the only differences between the Hartree and Fock terms come from the sign change (from the fermionic loop in the Hartree term) and from the summation over spins. The total Coulomb self-energy is Note that the self-energy for up spin is determined from the propagator for the down spin and vice versa. For electron-phonon interaction we only include the Fock term. The Hartree term vanishes because there is no p=0 -phonon mediating the Hartree interaction.
(S10) This is written without particle-hole (Nambu) basis and must be extended to include superconductivity. For now, we consider the normal state to find what would be the stable phase if superconductivity would not be present. The total self-energy for spin σ and surface ρ is Below we use this to study the possibility of a pseudospin analogue of the magnetic state.
where we isolate a correction to the chemical potential as ∆µ. We assume a fixed particle number, and therefore this correction is counteracted by a shift in the chemical potential to the opposite direction and as a result, it vanishes. In the normal state, we can write G as a 4 × 4 matrix. Its inverse is where the basis is chosen as . The matrix can be inverted in 2 × 2 blocks labeled with spin: where The self-consistency equations for pseudo-spin magnetism are We see from these equations that a pure pseudo-spin magnetism can be ruled out, as the interactions do not support it: both interaction terms in Eq. (S20) are negative; compare with Eqs. (S23,S35). We would need to have a repulsive electron-phonon interaction or an attractive Coulomb interaction in order to obtain a nonzero solution. The unequal form of the interactions in spin-magnetism versus pseudospin-magnetism originates from the lack of an exchange term in the electron-phonon interaction.

Appendix E: Antiferromagnetism
For antiferromagnetism (AFM), the self-energy has the form where h n is the frequency-dependent exchange field. As for the pseudospin magnetism above, the self-consistency equations for AFM are In contrast to Eq. (S20), the sign of u is now positive, and this makes the normal state unstable to the AFM phase.
Here we take the momentum cutoff to infinity. For large N the results do not depend on the cutoff provided it is larger than p FB . For N ≤ 2 the momentum integration diverges without a cutoff, so the results below do not apply for those cases.
To numerically solve the above equations, we need to impose a cutoff ω max in the Matsubara summation. To do that, we need to replace the Coulomb interaction with a modified term using the pseudopotential trick 38 . The form of the pseudopotential depends on the details of the self-consistency equation, so it has to be formulated separately for different equations. The differences are in the details, and the basic idea stays the same: in the Coulomb part of the self-energy we divide the Matsubara summation to low-energy and high-energy parts. Then we solve for the self-energy term and define a new interaction constant, which takes into account the high-energy contribution 38 . The exact form of the contribution from the high energy sum is the part which differs between different equations. We formulate equations in a form which is easy to solve numerically.
For AFM equations, we can calculate the high energy part to be where the cutoff x max is chosen so that ε(p FB x max ) 2T . An argument larger than 4 is already large enough in order to approximate the hyperbolic tangent by unity, so we can choose x c = (8T /ε 0 ) 1/N . In total, the pseudopotential is For a better accuracy at the strong coupling and low-temperature regime, we also include the Coulomb part of the exchange field h c in Eq. (S24), so that α depends self-consistently on h c . The non-self-consistent α given above is sufficient for the calculation of T C . With cutoff and a pseudopotential, Eqs. (S22-S23) become which can be solved numerically.

Linearized equations for solving the critical temperature
To determine the critical temperature, the self-consistency equations can be linearized with respect to h. In this case the momentum integrals can also be done analytically. We obtain where α N = N sin(π/N )/π. If λ = 0, then Z = 1 and also the Matsubara summation can be done analytically. We find T m C in terms of the Riemann ζ-function, which approaches the value T m C = u/4 when N → ∞.
interaction vertex gets an additional Nambu structure; P ρ is replaced by P ρ τ 3 . The Hartree-Fock self-energy terms areΣ We note that the Hartree term only affects the normal-state self-energy components, and not the off-diagonal ones associated with superconductivity. Superconductivity is determined only from the Fock terms. The order parameter is the same on both sublattices, φ A = φ B . In principle, the order parameter could also have a different phase and magnitude on the two surfaces, but this choice is the one with the lowest energy 24 . The self-consistency equations for the superconducting phase are In Eq. (S35) we can impose a Matsubara cutoff ω max if we simultaneously replace u with a pseudopotential u + , like in Eqs. (S24, S27). Compared to the AFM case, there is a sign change in the pseudopotential, The interpretation of the pseudopotential is that for superconductivity, the scattering at high energies reduces the effect of interaction at lower energies, whereas for magnetism, the effect is reversed. If we remove the high-energy scattering from the theory, the interaction has to be replaced by an effective pseudopotential to account for their effect.

Correction to the critical temperature due to retardation
To obtain an approximation for the effect of the electron-phonon renormalization on superconductivity, we first solve the renormalization function Z at zero temperature by approximating the self-consistency equation as where = 0 + is used to regularize the integral. Inside the integral we approximate Z by its peak value Z 0 ≡ Z(iω = 0). The integral yields For ω = 0, we have Z 0 = 1 + λ/Z 0 , whose solution in the first order in λ is Z 0 = 1 + λ. Now the approximate self-consistency equation for φ is At zero temperature, for ∆ = φ/Z, we have

Effect of Debye dispersion
In a real material, the phonon dispersion is obviously not described by the Einstein model adopted in the text. To get an idea how much the exact phonon dispersion affects the results, we consider also the Debye model. For the Einstein model phonons there is no momentum dependence in the interactions, and for this reason the self-energy is also independent of momentum. For Debye phonons the interaction obtains a momentum dependence through the phonon dispersion. Without calculating the full momentum dependent theory, we can estimate the effect of the dispersion by estimating the typical energy of the exchanged phonon and the average interaction constant.
The maximum phonon energy exchanged within the flat band is limited by the flat band diameter 2p FB to be ω 0 = p FB ω D /q M , where ω D is the Debye energy and q M is the maximum phonon momentum. As typically q M is of the order of the size of the Brillouin zone, ω 0 ω D and the energy scale is reduced. On the other hand, we find that the dimensionless interaction constant is enhanced: λ ∝ 1/ω 2 0 . However, because ∆ 0 ∝ λ 0.2 ω 0 , the magnitude of ∆ 0 is restricted by ω 0 and the total effect is a smaller critical temperature than with the Einstein phonons with energy ω E equal to the Debye energy.
The above concerns interactions within one flat band. We can also consider Debye phonons in the context of two flat bands separated by the distance p d p F B in momentum space. The momentum range of Debye phonons connecting the two parts of the Brillouin zone is limited to p d , and they can be treated as if they had a constant energy ω 0 = ω D p d /q M . In this case they act essentially as Einstein phonons with an effective energy scale ω 0 .
The self-consistency equations are There are two local minima, which correspond to the magnetic and superconducting phases, and a saddle point which is the unstable coexistence solution. (b) Phase diagram with fixed u. When λ < u and T < T m C , the thermodynamically stable state is the magnetic state and when λ > u and and T < T sc C the stable state is the superconducting state. The dividing line λ = u is marked with the red dotted line. On the striped region, both phases are possible as metastable states.
Above, the product over the Matsubara frequencies is evaluated using a standard Matsubara trick. The overall constant cancels against the normal state partition function Z 0 . The free energy relative to the normal state is where the momentum sum gives the multiplicative factor C > 0. The self-consistency equations are given as derivatives of the free energy with respect to fields ∆ and h: As a check, we see that if there was no frequency dependence in the interactions, Eqs. (S55) and (S56) would give similar equations after the Matsubara summation. The difference is that in the full model with instantaneous interactions, there is effectively only one interaction constant λ eff = λω E −u.
The self-consistency equations (S62) and (S63) for ∆ and h should be solved simultaneously. At low temperatures, T h, ∆, we can approximate the hyperbolic functions with exponentials, and obtain where ∆ 0 = λ/2 is the zero-temperature order parameter for h = 0. For h in terms of ∆, an analogous expression can be found, with h 0 = u/2. At zero temperature, coexistence is only possible if u = λ and ∆ = h = λ/4. This coexistence point is a saddle point of free energy. This is illustrated in Fig. 7a. Similar kind of solutions are also found at a finite temperature. Assuming a second order phase transition from the normal state to a superconducting state, we can linearize (S62) with h = 0 to find a critical temperature T sc C = λ/4. We will see below that when u > λ the phase transition is actually of the first order from a magnetic state to the superconducting state, but we still use the above definition for T sc C to set a temperature scale. Similarly, assuming a second order phase transition from normal state to a magnetic state, we linearize (S63) to find a critical temperature T m C = u/4, which we take as the definition of T m C .
1. Stability of the phase with lower TC From the free energy we see that when λ < u, the magnetic phase is stable and the superconducting phase is either metastable or unstable. When λ > u, the roles are reversed. The solution is (meta)stable if it is a local minimum of the free energy. The criterion is When h = 0 or ∆ = 0 the cross derivatives vanish, and the stability condition becomes ∂ 2 F /∂h 2 > 0 and ∂ 2 F /∂∆ 2 > 0. Let us consider the (meta)stability of the superconducting phase when 0 < λ < u and T < T sc C , where T sc C = λ/4. The stability condition becomes At low temperatures, the first term in the parentheses dominates and the superconducting phase is metastable. Near T sc C the second term dominates and the superconducting phase becomes unstable against a spontaneous magnetization. In this case, ∆ ≈ 0, and The transition temperature T * at which the system becomes unstable can be determined from the condition ∂ 2 F /∂h 2 = 0, which is equivalent to solving the magnetic critical temperature by linearizing Eq. (S62) with ∆ as solved from Eq. (S63) with h = 0. For λ > u > 0, the magnetic phase is the metastable one, and the equations apply after interchanging u ↔ λ, h ↔ ∆ and T sc C ↔ T m C . The numerical solution for the transition temperature of both phases is shown in Fig. 7b as the dashed line.
Similar stability analysis can be used with the full model, except in that case the frequency dependence of the selfenergy functions complicates the situation. We can however simplify the problem by projecting the order parameters on the linearized solution.
To do this, we linearize Eqs. (S53)-(S56) with respect to h and d. We assume the singlet order parameter φ to be finite. The equations for Σ ω and φ are (S70) The equations for h and d are coupled to one 2 × 2 matrix, With the high-frequency cutoff, the equation can be written as a 2M × 2M matrix, where M is the number of Matsubara frequencies below the cutoff. The critical temperature T * of this AFM/triplet phase is then determined numerically by finding the temperature at which the largest eigenvalue of the matrix becomes larger than unity 23 .

Stability of the magnetic phase
Conversely, let us now consider parameters λ, u such that 0 < T m C < T sc C and determine the temperature T * above which the magnetic phase is unstable against a superconducting instability. Now we linearize Eqs. (S53)-(S56) with respect to φ and d. We assume the field h to be finite. The equations for Σ ω and h are (S73) The equations for φ and d are coupled, and can be written as a matrix equation Again, with the Matsubara cutoff, this is a matrix equation and the critical temperature T * can be solved by searching for the temperature at which the largest eigenvalue crosses 1. The solution curve for both cases, T m C < T sc C and T m C > T sc C , is shown as a dashed boundary line in the phase diagram, Fig. 1 in the main text.