{'fraction_non_alphanumeric': 0.04373320302956267, 'fraction_numerical': 0.017444417297825556, 'mean_word_length': 4.077152071446291, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High-temperature surface superconductivity in rhombohedral graphite

Surface superconductivity in rhombohedral graphite is a robust phenomenon which can exist even when higher order hoppings between the layers lift the topological protection of the surface flat band and introduce a quadratic dispersion of electrons with a heavy effective mass. We show that for weak pairing interaction, the flat band character of the surface superconductivity transforms into a BCS-like relation with high critical temperature characterized by a higher coupling constant due to a much larger density of states than in the bulk. Our results offer an explanation for the recent findings of graphite superconductivity with an unusually high transition temperature.

Surface superconductivity in rhombohedral graphite is a robust phenomenon which can exist even when higher order hoppings between the layers lift the topological protection of the surface flat band and introduce a quadratic dispersion of electrons with a heavy effective mass. We show that for weak pairing interaction, the flat band character of the surface superconductivity transforms into a BCS-like relation with high critical temperature characterized by a higher coupling constant due to a much larger density of states than in the bulk. Our results offer an explanation for the recent findings of graphite superconductivity with an unusually high transition temperature. A low critical temperature of conventional superconductors results from a constant density of states (DOS) due to a linear-in-momentum electronic spectrum near the Fermi energy. With a higher-order dispersion, the relation between the critical temperature and the coupling constant becomes stronger, boosting the superconductivity. The extreme case would be a completely dispersionless energy spectrum, a flat band, which has been predicted in many condensed matter systems, see e.g. Refs. 1-4. In some cases the flat bands are protected by topology in momentum space; they emerge in gapless topological matter [5][6][7][8][9][10][11][12][13]. A singular DOS associated with the dispersionless spectrum was recently shown [14] to essentially enhance the transition temperature opening a new route to room-temperature superconductivity.
The problem is to find the metal with such a higherorder dispersion around the Fermi sea. Refs. 9 and 14 have shown that within the nearest-neighbour approximation, rhombohedral graphite (RHG) has topologically protected surface states with a flat band at the Fermi energy, and these surface states support high-temperature superconductivity where the superconducting order parameter is concentrated around the surfaces. A flat band forms out of a low dispersive band that appears on the surface of a multilayered rhombohedral graphene structure with a large number of layers. The corresponding critical temperature depends linearly on the pairing interaction strength and can be thus considerably higher than the usual exponentially small critical temperature in the bulk. Flat-band superconductors can carry quite high surface supercurrent with the critical value proportional to the critical temperature [15].
Experimental evidence of a high-temperature superconductivity in graphite in the form of a small Meissner effect and of a sharp drop in resistance appeared in the literature during past years [16,17]. Recently, these findings have been ratified by observations of zero resistance in graphitic samples up to 175 K [18] and indi-cations of even room-temperature superconductivity in specially prepared graphite samples [19]. The enhanced superconducting density has been also reported on twin boundaries in Ba(Fe 1−x Co x ) 2 As 2 [20]. In this Letter, we argue that the high-temperature superconductivity in graphite can be related to surface superconductivity that may form either on the outer surfaces of the sample or on twin boundaries of or on grain boundaries between inclusions of RHG. We demonstrate that the surface superconductivity is a robust phenomenon which survives even when the topological protection of the flat band itself is lifted. In particular, the next-nearest neighbour hoppings in RHG can break the exact topological protection and, therefore, the flat-band mechanism of superconductivity could be destroyed. We study these higher-order interactions [21] and show that, though breaking the flat-band scenario at sufficiently low values of the coupling energy, they provide another mechanism of surface superconductivity which is of the BCS type but still has a much larger coupling constant than the usual superconductivity in bulk, thus favoring high-temperature superconductivity. The enhanced coupling constant comes from a high DOS associated with a heavy effective mass of surface quasiparticles emerging on the background of the pre-existing flat band. Our results help to identify the regime of parameters where extremely high-temperature surface superconductivity may be found.
Electron dispersion in RHG. The RHG lattice and the tight-binding couplings are depicted in Fig. 1. We label the layers (starting from the bottom) by n, the atoms A in layer n are on top of atoms B in layer n − 1, and the vector between the layers is d. In our estimates, we use the tight-binding parameters denoted in Fig. 1. They satisfy γ 0 ≫ γ 1 ∼ γ 3 ≫ γ 4 [22]. In the tightbinding numerics below, we use the values γ 0 = 2.58 eV, γ 1 = 0.34 eV, γ 3 = 0.17 eV, and γ 4 = 0.04 eV, which give the best fit to the density functional theory (DFT) calculation of the surface state dispersion (see Appendix). The RHG is a multilayered graphene structure. The conical spectrum near the Dirac points K and K ′ of the Brillouin zone of a single-layer graphene (for details, see review [22] and references therein) is transformed into low-dispersion, low-energy bands (see Fig. 2) which determine the unique features of this system. Since we are interested in low energies we concentrate on the in-plane momenta p = (p x , p y ) close to one of these Dirac corners. A standard Fourier series expansion near K yields [21] whereĤ mn (K, p) = 4 l=0Ĥ (l) mn (K, p) and Hereγ 3 = γ 3 /γ 0 ,γ 4 = γ 4 /γ 0 , p ± = p x ± ip y = pe ±iφ , and v F = 3a 0 γ 0 /2 . The Pauli matricesσ and 2σ ± = σ x ± iσ y act on pseudo-spinorsψ n = (ψ 1 n , ψ 2 n ) T ,ψ † n = (ψ 1 * n , ψ 2 * n ), where ψ 1 n = ψ A n , ψ 2 n = e iπ/6 ψ B n . To construct the associated Bogoliubov-de Gennes (BdG) Hamiltonian for the superconducting state we need also the time-reversed "hole" Hamiltonian for the Dirac point K. It follows from the particle Hamiltonian in a vicinity of the opposite Dirac point −K which is equivalent to K ′ . The wave function ψ h K of a hole excitation near K is ψ h K =ψ * −K . One can check that the hole Hamiltonian is H h mn (K, p) = H * mn (−K, −p). Therefore, In what follows, we denote the electron wave function bŷ u n =ψ n and the hole wave function byv n =ψ h n . In the next section we study the normal state using Hamiltonian Eq. (1) and the relative magnitudes of the coupling constants listed above. The results of the numerical solution are displayed in Fig. 2, along with the corresponding analytical approximations and DFT calculations. Low-energy spectrum in the normal state. The Schrödinger equation takes the form mĤ nm (K, p)û m (p) = (ǫ + µ)û n (p) . ( The energy is measured from the chemical potential µ. The energy spectrum in bulk is obtained by ignoring the outermost layers n = 1 and N and using the ansatẑ u n ∝ e ipz dn , where p z is out-of-plane momentum. For zero doping µ = 0, the Fermi surface is determined by ǫ(p, p z , φ) = 0. If γ 3 = γ 4 = 0, the Fermi surface shrinks to a spiral line v F p = γ 1 , φ = p z d + π/6. Projection of this spiral onto the momentum plane q = 0 determines the area of a flat band for surface states [9] in the limit N → ∞. If only γ 4 = 0 while γ 3 = 0, equation ǫ = 0 for µ = 0 can still be shown to give a Fermi surface in the form of a (corrugated) spiral whose projection determines the surface flat band [21]. This is because the interaction comes with matricesσ x andσ y so that the full Hamiltonian obeys the same anti-commutation rule [σ z , (Ĥ (0) +Ĥ (1) +Ĥ (3) )] + = 0 as the initial Hamiltonian H (0) +Ĥ (1) . This preserves the same topological invariant and the same topology of the Fermi surface [9]. Since γ 3 does not affect the presence of a flat band, we restrict our analytical consideration to the case when only γ 4 is nonzero while γ 3 = 0 for simplicity. However, our numerical analysis is carried out using the full Hamiltonian.
Surface states have complex p z = p ′ z + ip ′′ z and decay into the bulk. For γ 3 = 0 Eq. (3) in the particle channel, for low energies has a solution in the form Here the out-of-plane momentum At the outermost layers, the terms withψ 0 andψ N +1 in Eq. (4) disappear. The components which do not have γ 1 couple the constants A + and A − in Eq. (5) and determine the energy of the surface states for ξ p , ǫ ≪ γ 1 . Here The interaction γ 4 breaks the symmetry between the conduction and valence bands in a way similar to a shift in µ due to doping. The spectrum ǫ p has a quadratic dispersion with the effective mass m * on a background of a much weaker high-order dispersion ξ p . The latter transforms into a flat band ξ p = 0 with a radius p < p FB for an infinite number of layers, N → ∞. The effective mass is much larger than the characteristic band mass m 3 in 3D graphite. Indeed, we have m * /m 3 ∼ γ 1 /γ 4 where we estimate 2 /(m 3 a 2 0 ) ∼ γ 0 as the conduction band width in graphite. We see that m * /m 3 ≫ 1. This dispersion is compared with the results of numerical diagonalization of H(K, p) in Fig. 2 usingγ 3 = 0.066, andγ 4 = 0.016.
BdG equations are constructed using the particle and hole Hamiltonians (1) and (2) coupled through the superconducting order-parameter field ∆. As distinct from the quasiparticle energy measured from the chemical potential upwards, E = µ + ǫ, the energy of holes is measured from µ downwards, E = µ − ǫ. We have Here we introduce objects in the Nambu spacě Each component of the Nambu vectorΨ n is a pseudospinor. For analytical consideration we assume that ∆ n = 0 for n = 1, N . This is justified by the numerical solution of the self-consistency equation using the full BdG equations [14,21]. In this case, Eq. (9) for n = 1, N does not contain ∆, so that one can use the normalstate solution, Eq. (5), where we have Nambu vectorš A ± = (A ± , B ± ) T instead of the corresponding scalars, and the Nambu matrixζ, Eq. (6), withτ 3 ǫ instead of ǫ.
At the outermost layers, the terms withû 0 ,v 0 and u N +1 ,v N +1 in Eq. (9) disappear. The components which do not contain γ 1 yielď (11) provide the surface-state spectrum and determine four independent surface states. If ∆ 1 = ∆ N , the spectrum isǫ 2 = (μ p ± ξ p ) 2 + |∆| 2 . If the number of layers N is large, ξ p → 0 for p < p FB , the two surface states decouplẽ In this case, Eq. (10) at layer n = 1 yields A + = U , Surface superconductivity. The surface states discussed above form a basis for the superconducting gap localized near outer surfaces. Within the mean-field approximation, the gap at layer n is determined by the self-consistency equation. As was shown in Ref. 14, the surface states dominate due to a much larger DOS. For a large number of layers when ξ p = 0, the self-consistency equation for the gap at the surface takes the form Here we used Eq. (13) to find the integrand, ǫ is one of the spectral branches in Eq. (12), the integration is carried out over momenta within the flat band, p < p FB , and W is the 3D coupling potential. This is the central result of our Letter. The superconducting coupling is described by the energy g = (W/d)p 2 F B / 2 . It can also be expressed in terms of the usual BCS coupling constant λ = ν 3 W where ν 3 = m 3 p 3F /2π 3 3 is the 3D density of states and p 3F is the Fermi momentum in 3D graphite. Assuming the conduction band width in 3D graphite of the order of γ 0 we have g/γ 1 ∼ λ(γ 1 /γ 0 ) if /a 0 p 3F ∼ 1.
The overall behavior of ∆ vs. the coupling energy is plotted in Fig. 3. Let us consider first the resulting ∆ for zero doping µ = 0. The quadratic dispersion due to γ 4 comes with an energy scale α = 2γ 4 γ 1 , which determines a crossover between exponentially suppressed  14) at µ = 0 and the red (top) line corresponds to µ = µopt that maximizes ∆ for given g. For large g, the gap tends towards the flat-band limit ∆ ∝ g. For g 4πα (upper inset), the gap is exponentially suppressed, ∆ ∝ exp(−4πα/g). The lower inset shows the (normalized) gap as a function of µ for a few values of g. and flat band superconductivity. For g ≫ 4πα and for zero doping, Eq. (14) yields the flat-band result [14], ∆ = g/8π for T = 0, and the critical temperature satisfying ∆ = 3k B T C . Due to its linear dependence on the interaction strength, the critical temperature is proportional to the area of the flat band and can be essentially higher than that in the bulk. Doping in the flat band regime destroys the surface superconductivity [14]. Both ∆ 0 and T c vanish at the critical doping level |µ| = 2k B T c .
For g ≪ 4πα the weak dispersion Eq. (7) with a heavy mass m * dominates. The integral in Eq. (14) is logarithmic which results in a BCS-like expression (for T = 0) where α = 2γ 1γ4 . The estimate for g gives λ 2 ∼ λ(γ 1 /γ 4 ). This is a much larger coupling constant than λ for bulk superconductivity. The gap disappears at µ = α. The crossover from the BCS-like to the flat-band regime occurs at g ∼ 4πα, and ∆ ∼ α. The coherence length ξ 0 = v g /∆ ∼ a 0 (γ 0 /γ 1 )e 1/λ2 is much longer than the interatomic distance a 0 . These analytical results are compared in Fig. 3 to the numerical solution of the selfconsistency equation using the full Hamiltonian, Eq. (9). In contrast to the flat-band and the BCS regimes, the gap in the intermediate region α ∼ g is enhanced by an optimum doping, i.e., the critical temperature is very sensitive to the presence of impurities (lower inset in Fig. 3). This complies with the reports of high-temperature superconductivity in doped graphite, Ref. 16.
Effect of fluctuations. The quality of the mean-field approximation used above is determined by the Ginzburg number Gi which is a measure of the relative magnitude of order-parameter fluctuations. For usual 3D superconductors Gi ≪ 1 due to a small ratio of the critical temperature to characteristic energy of electrons (i.e., the Fermi energy). Here we demonstrate that the mean field approach also works well when the quadratic dispersion dominates over the flat band. In this case the fluctuation free energy density for T not too close to T c is F 1 ∼ ν 2 ∆ 2 1 /2, where ν 2 = m * /2π 2 is the 2D DOS and the effective mass m * is determined by Eq. (8). The energy of an area πξ 2 0 with a radius of When the quadratic dispersion dominates, one has T c ≪ γ 1γ4 with the Ginzburg number Gi = e −1/λ2 ≪ 1, thus the average fluctuation of the order parameter is small compared to its mean-field value. However, at the crossover to the flat band regime, the fluctuation becomes of the same order as the mean-field ∆ 0 ∼ γ 1γ4 . Therefore, the mean-field approach is not exact for the flat-band regime. Nevertheless, it is used here as an initial step towards a full theory of high temperature surface superconductivity.
Summary. Rhombohedral graphite is a promising candidate for high-temperature surface superconductivity due to its (approximate) topologically protected flat band. Besides surfaces, similar type of superconductivity may arise around stacking faults and interfaces between differently stacked regions of graphite, as long as the system contains more than a few layers of RHG regions [21]. Recent observations of high-temperature superconductivity in graphite [16][17][18][19] are compatible with surface or interface superconductivity described by our theory if there are RHG regions embedded inside otherwise Bernally stacked graphite. Our predictions can be used for search or for an artificial fabrication of layered and/or twinned systems with high-and even room-temperature superconductivity. With the hopping parameters used above, the crossover between the flat-band and BCS-like regimes takes place around g c ∼ 4πα ≈ 0.39γ 1 ≈ 0.15 eV (see Fig. 3 and Appendix) corresponding to the mean-field T c (g c ) ≈ 20 K. For g > g c we thus find T c ∼ (g/γ 1 ) × 50 K. This is much greater than the expected gap in the bulk for the same magnitude of coupling.
We thank G. Volovik for helpful comments and the collaboration that initiated this project. We also acknowledge fruitful discussions with F. Mauri. This work is supported in part by the Academy of Finland through its Centre of Excellence Program (projects no.  The in-plane lattice parameter was optimized for monolayer graphene and it was found to be 2.466Å. The same lattice parameter was used for the rhombohedral slabs but the distance between adjacent layers was allowed to relax freely until forces acting on atoms were less than 0.001 eV/Å. This resulted in variations of the order 0.01Å in the interlayer distances, the outer layers being slightly expanded from the optimized interlayer distance 3.332Å of rhombohedral graphite bulk. The use of the lattice parameter of graphene leads to negligible errors, as the difference to the optimized bulk parameter, 2.462Å, is very small, and the multilayer structures are expected to interpolate between the bulk and monolayer limits. The band structure in the vicinity of the K-point was fitted to a tight-binding model, in which hoppings between nearest neighbors in-plane (γ 0 ) and out-of-plane (γ 1 ), as well as between next-nearest neighbors out-of-plane (γ 3 and γ 4 ) were included. The fit was concentrated to the region around K using a Gaussian function to weigh the squared error between the tight-binding and DFT bands, and the minimal energy of the parabolic region was shifted to match the corresponding DFT energy. The width of the Gaussian was chosen such that the weight outside the flat band region was practically zero. Fig. 4(a) demonstrates a comparison between the DFT band structure and the corresponding tight-binding fit, as well as the shape of the weighing function on a 20-layer slab.
Two parameter sets with comparable agreement around the K-point were found when the width of the weighing Gaussian was altered. As Fig. 4(b) and (c) show, one of them better captures the overall band structure. It was thus chosen for the further calculations. Both parameter sets are reported in Table I. The differences in the parameters is insignificant for the considerations of the surface superconductivity in the main paper, and they both give reasonable Fermi velocities (v F,1 =0.84·10 6 m/s and v F,2 =1.04·10 6 m/s, respectively) .

DETERMINING THE EFFECTIVE MASS DIRECTLY FROM A PARABOLIC FIT TO THE DFT BANDS AROUND K
Additionally, the DFT-calculated dispersion very close to the K-point was fitted directly using a parable [Eq. (8) in the main text], and a value for the effective mass m * and chemical potential µ were extracted separately for both the K − Γ and K − M directions, in addition to a mean value fitted simultaneously to both directions. These values were compared to those obtained from the In both plots, blue symbols refer to fit 1 and red symbols to fit 2 (Table I). tight-binding fit. Fig. 5(a) shows the m * and µ as a function of the slab thickness as well as those calculated from Eq. (9) based on the fitted tight-binding parameters γ i . The effective mass increases with an increasing number of layers, as the breadth of the flat band increases. In DFT, there is also a small linear component in the parabolic fit E − E F = β 0 (∆k) 2 + β 1 (∆k) − µ due to the longer-range couplings neglected in the tight-binding model. The effect of this component is more important further away from K-point, as it decreases when the fitting region is made narrower. In DFT, the minimum of the parabolic bands does not lie exactly at the Fermi energy. This finite doping is, however, small, µ/α < 0.1 in all DFT calculations.
The parameter α determines the crossover between the BCS-and flat-band regimes. It is worth noting that even though the values for the hopping amplitudes differ in the tight-binding parameters sets, both yield comparable values for α, as illustrated in Fig. 5(b). Moreover, the