This is the submitted version of the article: Saunderson T.G., Annett J.F., Újfalussy B., Csire G., Gradhand M.. Gap anisotropy in multiband superconductors based on multiple scattering theory. Physical Review B, (2020). 101. 064510: - . 10.1103/PhysRevB.101.064510. Available at: https://dx.doi.org/10.1103/PhysRevB.101.064510 Gap Anisotropy in Multiband Superconductors Based on Multiple Scattering Theory Tom G. Saunderson,1, ∗ James F. Annett,1 Balázs Újfalussy,2 Gábor Csire,3, 2, 1 and Martin Gradhand1 1 HH Wills Physics Laboratory, University of Bristol, Tyndall Ave, BS8 1TL, United Kingdom 2 Institute for Solid State Physics and Optics, Wigner Research Centre for Physics, Hungarian Academy of Sciences, P.O. Box 49, H-1525 Budapest, Hungary 3 Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC, BIST, Campus UAB, Bellaterra, Barcelona, 08193, Spain (Dated: July 1, 2020) We implement the Bogoliubov-de Gennes (BdG) equation in a screened Korringa-Kohn-Rostoker (KKR) method for solving, self-consistently, the superconducting state for 3d crystals. This method combines the full complexity of the underlying electronic structure and Fermi surface geometry with a simple phenomenological parametrisation for the superconductivity. We apply this theoretical framework to the known s-wave superconductors Nb, Pb, and MgB2 . In these materials multiple distinct peaks at the gap in the density of states were observed, showing significant gap anisotropy which is in good agreement with experiment. Qualitatively, the results can be explained in terms of the k-dependent Fermi velocities on the Fermi surface sheets exploiting concepts from BCS theory. I. INTRODUCTION functional theory (DFT) calculation. This often results in over-simplification of the normal state band structure The structure of the superconducting gap in s-wave in order to construct an efficient model to describe the su- phonon mediated superconductors may show a surpris- perconducting state. While this usually leads to deep un- ingly complex structure. Not only can it show a large derstanding of aspects of the superconductivity, it might degree of anisotropy on the Fermi surface of a given ma- obscure the importance of the complexity and orbital hy- terial but several known examples exhibit clear signatures bridisation of the underlying electronic structure. of two superconducting gaps. The most notable of cases In this work we aim to follow a route distinct to those are Pb [1–4] and MgB2 [5–7] and to a lesser extent Nb two main directions. We aim to describe the full complex- [8–12] with conflicting reports [13, 14]. While it is un- ity of the normal state electronic structure whilst treat- derstood that the gap anisotropy arises in these systems ing the superconducting pairing interaction within a sim- from multiple Fermi surface sheets, not many phonon plified model. This idea is very similar to the LDA+U mediated s-wave superconductors showing multiple gaps method where the Hubbard U interaction is also treated have been identified. as a tunable parameter [39]. This has proven to be In the original BCS theory [15] only a single spherical very successful in modelling strongly correlated materi- band was considered and so it was impossible to obtain als without compromising on the full electronic structure any gap anisotropy. In addition, it was limited to the within a strongly correlated state. weak coupling regime, which later was rectified by the In analogy to this, our method leads to a full electronic Eliashberg theory [16, 17] and more recently by full DFT structure within the superconducting state giving access approaches to electron phonon coupling driven supercon- to the full gap anisotropy in multi sheeted s-wave su- ductivity [18, 19]. Although gap anisotropy analysis is perconductors. We build upon previous work from some present in such codes, much of the DFT based ab initio of the authors on the implementation of the scalar rela- work focused on other aspects. These aspects being; ex- tivistic [40–42] and fully relativistic [43, 44] BdG equa- panding the description of s-wave superconductivity with tions in a layered KKR DFT code. The advantage of emphasis on the correct description of the driving mecha- this approach is the access to the full normal state elec- nism [20–24], or on the extension to unconventional pair- tronic structure from first principles without using sim- ing pushing the boundaries in our treatment of iron based plified models. By just considering s-wave pairing and [25], spin [26] and magnetic effects [27, 28]. one effective pairing parameter, we show quantitatively Beyond the ab initio work on superconductivity, there the complexity of the superconducting gap including its is also extensive work on using parametrised models to full anisotropy even in simple elemental crystals. describe superconductivity [29–38]. These models are The work is structured as follows. First, we will outline a powerful tool for describing unconventional supercon- the basic theoretical background of the method includ- ductors and some of the basic principles of standard ing the differences in implementation to the earlier work s-wave superconductors. However these models have [40–42]. This will be followed by numerical tests to de- the drawback that they require the normal state to be termine the robustness and accuracy of the theoretical parametrised from either experimental data or a density framework, namely the solution of the scalar relativistic BdG equations. In the next part we will present several examples of simple s-wave superconductors, the resulting gap anisotropy on the various Fermi surface sheets, and ∗ t.saunderson@bristol.ac.uk some simple arguments on how the normal state proper- 2 ties drive the observed anisotropies. In all cases we will or non-locality in the pairing potential. In this report Λ is relate our work to experimental observations. tuned such that the gap in the density of states matches experimental results for the zero temperature gap size ∆(T = 0) of the material in question. Further details of II. METHOD this are discussed in section III. The KKR method exploits a local atomic basis that This section will describe the implementation of super- uses multiple scattering theory, which gives access to the conductivity into the existing KKR code [45, 46]. The full Green’s function of the system. For details on the im- theory for the single site solver and multiple scattering plementation we refer to G. Csire et al [40]. Here we re- terms have been introduced by G. Csire et al [40]. All strict the discussion to the essential components to high- equations are given in Rydberg atomic units. light differences in the implementation and to shed light The density functional theory for superconductors was onto the results. The Green’s function for the system is initially presented by L. N. Oliveira, E. K. U. Gross defined as and W. Kohn [47], who introduced the effective pair-  ee  ing potential ∆ef f (r, r0 ), describing the superconducting Ĝ (z) Ĝeh (z) −1 ĜBdG (z) = he hh = z Iˆ − ĤBdG , (9) state in addition to the conventional Kohn-Sham poten- Ĝ (z) Ĝ (z) tial Vef f (r). These two potentials are defined as Z where ĤBdG (r) = hr|ĤBdG |ri and ρ(r) δExc [ρ, χ] Vef f (r) = Vext (r) + d3 r 0 + , (1) |r − r | δρ(r)   Ĥ(r) ∆ef f (r) ĤBdG (r) = , (10) δExc [ρ, χ] ∆ef f (r)∗ −Ĥ(r)∗ ∆ef f (r, r0 ) = , (2) δχ(r, r0 ) Ĥ(r) = −∇2 + Vef f (r) − µ. (11) 0 where χ(r, r ) is the anomalous density and Exc [ρ, χ] is the exchange correlation functional for a superconductor. Here, µ is the chemical potential and z =  + iδ. The The exchange correlation functional was later approxi- densities ρ(r) and χ(r) can be calculated by taking the mated [48] as trace of different components of GBdG (z, r, r0 ), 1 ∞ Z Z 0 Exc [ρ, χ] =Exc [ρ] − d3 r1 d3 r2 d3 r10 d3 r20 (3) ρ(r) = − df ()ImTrGee (, r, r0 ) π −∞ χ∗ (r1 , r01 )Λ[ρ, χ](r1 , r01 , r2 , r02 )χ(r2 , r02 ), 1 ∞ Z − d[1 − f ()]ImTrGhh (, r, r0 ), (12) π −∞ 0 where Exc [ρ] is the normal state exchange-correlation Z ∞ 1 functional, and Λ[ρ, χ](r1 , r01 , r2 , r02 ) is the pairing kernel. χ(r) = − d[1 − 2f ()]ImTrGeh (, r, r0 ) Similarly to ref. [48], in the KKR framework we consider 4π −∞ Z ∞ the atomic sphere approximation (ASA). The kernel is 1 − d[1 − 2f ()]ImTrGhe (, r, r0 ), (13) then approximated to 4π −∞ Λ[ρ, χ](r1 , r01 , r2 , r02 ) = Λδ(r1 − r01 )δ(r1 − r2 )δ(r1 − r02 ), where the limit is taken such that δ → 0+ . We use equa- (4) tions (12) and (13) to find new expressions for Vef f (r) inside the ASA spheres and and ∆ef f (r) using equations (1) and (8). From here a Λ[ρ, χ](r1 , r01 , r2 , r02 ) = 0, (5) new ĤBdG (r) is constructed to calculate a new Green’s function ĜBdG (, r, r0 ), and thus self consistency can be outside, where Λ is called the interaction parameter. A achieved. further conventional simplification is for ∆ef f (r, r0 ) and The ASA approximation sets a boundary to each χ(r, r0 ) to be local, resulting in atomic site, i, called the ASA radius riASA . This approx- imation implies that the potentials Vef f (r) and ∆ef f (r) ∆ef f (r, r0 ) = ∆ef f (r)δ(r − r0 ), (6) can be written in sums χ(r, r0 ) = χ(r)δ(r − r0 ). (7) X Vef f (r) = Vi (r), (14) Equation (2) therefore becomes i X ∆ef f (r) = Λχ(r). (8) ∆ef f (r) = ∆i (r), (15) i When the spherical symmetry (ASA) and local approxi- mation for the pairing potential are combined, we restrict and ensures that Vi (r) = 0 and ∆i (r) = 0 if |r| ≥ riASA . ourselves to s-wave superconductivity. Any kind of non-s The analogue of equation (8) becomes, wave superconductivity would need a non-spherical pair- ing potential, coupling between different orbital channels ∆i (r) = Λi χi (r). (16) 3 The resulting code is therefore able to perform calcu- more dominant with larger atomic number as expected lations which relax both charge and anomalous densities when including scalar relativistic corrections. along with the Fermi energy. In practice, we first relax the densities and only in the last steps allow the Fermi energy to relax as well. We found that relaxing the Fermi energy does not significantly change the solution and is computationally expensive. The results we present in this report are therefore fixed Fermi level calculations. III. COMPUTATIONAL DETAILS First, we test the numerical parameters and the result- ing gaps within our framework. We define the average gap [48] Z ∆¯r = 1 ∆ef f (r)d3 r , (17) FIG. 2. The comparison of the average gap size of Nb, V and VW S W S Cu for non relativistic and scalar relativistic calculations as a function of Λ. where VW S represents the volume of the Wigner-Seitz cell. This is a practical expression, however, in a real bandstructure this average gap is not the same as the Within this framework, Λ is the only free parameter. actual gap in the density of states. In general this parameter is tuned such that the result- Fig. 1 shows the results of the convergence tests for ing gap matches the experimental zero temperature gap non-relativistic (NR) Niobium. For low Λ, all calcula- size. The complication, especially for anisotropic gaps tions with 30 energy points are suppressed relative to or multigap systems, is what definitions are used exper- calculations with a higher energy meshes. At high Λ the imentally and theoretically to establish the gap size. In solutions with the same sized k-mesh group together. We the following we discuss how in practice we constrain our decided that 50 energy points and 200 × 103 k-points is free parameter Λ. The first definition is using the aver- the best trade off between accuracy and computational age gap from equation (17). This quantity, however, is cost. just the integrated anomalous density of states χ(r) mul- tiplied by a constant, and does not necessarily relate to the gap seen in the density of states. It can been easily proven comparing the gaps in Fig. 3 and ∆ ¯ r for Nb in table I. Clearly these numbers are different from each other. In fact, by comparing the values for ∆ ¯ r for Pb and Nb it is possible to see that Nb has a smaller ∆ ¯r than Pb. As Nb has a larger Tc than Pb it is obvious that the real gap size must be larger. We believe that the fundamental reason for this discrepancy lies in the fact that the relationship between the effective pairing potential ∆ef f (r) and the k-dependent gap ∆(k) is not trivial. Similarly, the average gap found in the literature (∆¯ k ) and shown in table I, column 6, is conventionally obtained from a k-space integration of ∆(k). This quan- tity is distinct to the average ∆¯ r calculated from the real FIG. 1. The convergence test performed on a non-relativistic calculation of Niobium. The parameters are the energy and space integral over ∆ef f (r) (see Table I, column 5). For k-mesh. The k-mesh is the actual number of k-points used Nb and Pb, ∆ ¯ r is smaller than ∆ ¯ k but for MgB2 it is to calculate the energy points closest to the superconducting larger. gap. Plot symbols denote 2 × 104 (circle), 2 × 105 (cross) and Finally, we would like to summarize the numerical pa- 5 × 105 (square) k-points, colours denote 30 (red), 50 (green) rameters used throughout the papers. After convergence and 60 (blue) energy points. of the potentials is reached, the imaginary part of the energy, δ = 2µRy, 0.2µRy, 10µRy is used for Nb, Pb For the remainder of the report, we perform calcula- and MgB2 respectively. These are the parameters for the tions using the scalar-relativistic (SR) BdG equations. A DOS and Bloch spectral function calculations when fo- comparison between NR and SR results is shown in Fig. 2 cussing on the gap structure as well as for ∆(k). The for Nb, V and Cu. There is no clear trend in terms of number of k-points in the irreducible Brillouin zone for the gap size going from NR to SR. For V and Nb the gap DOS calculations are 1 × 107 for Nb and Pb and 2.7 × 106 is reduced, for Cu it is increased. The influence becomes for MgB2 . 4 ¯ r (meV) Element D(F ) (Ry−1 ) Λ (Ry) ∆ ∆exp (meV) ¯ k (meV) ∆ Nb 17.72 0.0861 1.05 1.56[49],(1.79, 1.64, 1.20)[12] 1.53[49] Pb 6.77 0.351 1.15 1.33[49],(1.40, 1.27)[4] 1.35[49],(1.0, 0.8)[50] MgB2 3.83 0.291 7.93 (7.0, 3.0)[7] (6.8, 2.45)[51] TABLE I. D(F ) is the density of states at the Fermi level in the normal state obtained from our calculations. Λ is the interaction parameter used in this investigation to match our calculations to the experimental zero temperature gap, more ¯ r , is calculated using (17). ∆exp are average gaps from experiments [4, 7, 12, 49], ∆ details in section IV. The average gap, ∆ ¯k are average gaps from theoretical ∆(k) integrations [49–51]. IV. RESULTS A. Niobium For calculations in the case of Nb, we tuned Λ such that the gap in the density of states around the Fermi level was matched to the experimental gap size via tunnelling experiments [12]. It predicted different sizes for the su- perconducting gap depending on the exposed surface of the single crystalline Nb. The different crystal planes investigated were [001], [110] and [111]. The values for the superconducting gap are given by ∆001 = 1.20meV, ∆110 = 1.79meV and ∆111 = 1.64meV. We chose to tune Λ such that the outer peak of our superconducting gap matched ∆110 . The result is shown in Fig. 3. In the inset we can identify two clear peaks and one weak shoulder corresponding to three distinct gaps at ∆1 = 1.43meV, FIG. 4. A Bloch spectral function of Niobium in the super- ∆2 = 1.69meV and ∆3 = 1.79meV, which is in reason- conducting state, showing the bandstructure in directions of able quantitative agreement to the experiment. In com- high symmetry. The units of the spectral function are arbi- parison to other literature the gap anisotropy is also well trary. In the top panel with a higher energy resolution the matched [8–11]. The lattice constant of a = 6.23685a.u. superconducting gap is just about visible. was used for the bcc structure of Nb. Fig. 3 the relevant energy resolution is −4.5meV to 0meV and for the bandstructure we focus on the high symmetry line Γ to N in Fig. 5. A double peaked superconducting gap is clearly resolved, with band gaps 1 and 2 contribut- ing to the outermost gap, and the 3rd band gap relating to the inner peak. In order to understand this effect we consider the orbital character associated with each band. While the inner peak is ‘p-d’ hybridised, the outer peak is almost entirely of a d-electron character. Typically, ‘p’ character bands will show a larger Fermi velocity at the Fermi level compared to the ‘d’ character bands. This is confirmed by showing the Fermi velocities as a colour map on the distinct Fermi surface sheets in Fig. 6, where FIG. 3. A scalar relativistic calculation of the electronic den- points of the Fermi surface associated with panels 1-3 sity of states of Niobium in the superconducting state, with from Fig. 5 are labelled. Evidently panels 1 and 2 have the inset showing the gap at the Fermi level F . similar velocities at the point of crossing, whereas panel 3 has a velocity which is noticeably larger. Since the DOS The lower panel of Fig. 4 shows the Bloch spectral is inversely proportional to the Fermi velocity, function for a large energy window but in the top panel we focus on the region around the Fermi level along the 1 Γ to N direction with a smaller energy window (−50meV D(F ) ∝ , (18) to 50meV). On that scale the opening of the supercon- vF ducting gap is just about resolved. In order to investigate the gap and the associated anisotropy as highlighted in and within the BCS theory [15, 52] the gap scales with 5 the DOS   1 ∆ ≈ 2kB ΘD exp − , (19) V D(F ) we can infer a gap anisotropy which is related to the dis- tribution of Fermi velocities on the Fermi surface sheets. In Fig. 7 we extend this analysis to the full Fermi sur- face and in comparing it to Fig. 6 a strong correlation between vF (k) and ∆(k) is observed. FIG. 7. The Fermi surface sheets of Nb in the normal state with the gap size of Nb in the superconducting state superim- posed as a colour scale on top. The labels refer to the panels in Fig. 5 and identify the points on the Fermi surface where the gaps in Fig. 5 appear. experiments by M. Ruby et al [4]. In this experiment the authors could identify two distinct peaks of the su- perconducting gap at an energy separation of 150µeV. Within our calculations the fcc crystal structure with a lattice constant of 4.95Å is used. Here we restrict the computation to three dimensional periodic crystals but FIG. 5. Left Panel: The total and orbital resolved density of future work will focus on describing the surface explicitly states of Niobium near the Fermi level. Panel 1-3: Band gaps in the Γ to N direction. Panels 1 and 2 are associated with as measured in the experiment. For the bulk material we the larger gap in the DOS, and panel 3 is associated with the assume that the gap anisotropy in Pb should at least be smaller gap. The units of the spectral function are arbitrary. of similar order as found in the STM experiments. FIG. 8. A scalar relativistic calculation of the electronic den- sity of states of Lead in the superconducting state with an inset figure showing the peaks of the superconducting gap. FIG. 6. The Fermi surface sheets with the color scale show- ing the Fermi velocity of Nb. The labels refer to the gaps For an interaction parameter of Λ = 0.351Ry the over identified in Fig. 5. all gap size is found to be comparable to Ruby et al [4]. The inset of Fig. 8 displays the gap structure at this interaction parameter where the separation between the distinct gaps is 140µeV. This energy separation is of B. Lead comparable to the 150µeV identified by Ruby et al [53]. In order to resolve such small separation a fine mesh of Lead often has disordered or amorphous crystal struc- 1 × 107 k-points in the irreducible part of the Brillouin tures. Experimentally, this has been avoided in STM zone was required. 6 FIG. 10. The Fermi surface sheets of Pb in the normal state with the gap size of Pb in the superconducting state superim- posed as a colour scale on top. The labels refer to the panels FIG. 9. Left Panel: The total and orbital resolved density of in Fig. 9 and identify the points on the Fermi surface where states of Pb near the Fermi level. Panels 1-3: Band gaps in the the gaps in Fig. 9 appear. Γ to K, X to U and Γ to X directions. Panel 1 is associated with the central coherence peak, panel 2 is associated with the small shoulder in the DOS and panel 3 is associated with the largest shoulder. The units of the spectral function are lead falls far from the weak coupling limit. Since BCS arbitrary. theory is valid in the weak coupling limit, therefore the relation between the Fermi velocity and the supercon- ducting gap is expected to be less robust. This aspect Following the same process as for Nb the distinct deserves further investigation and in particular the in- gap sizes can be traced to different points in k-space fluence of surfaces and spin-orbit coupling has to be ad- as shown in Fig. 9. However, to get the full picture dressed in future work. of the anisotropy Fig. 10 shows the size of the gap on both Fermi surface sheets, with the relevant directions and band crossings from Fig. 9 highlighted. It is clear that the Fermi sheet in the left panel is mainly associ- ated with the larger of the two gaps, and the sheet in the right panel contributes to the smaller shoulder. However, both sheets do contribute to a lesser extent to the other gaps as well. This is in contrast to the argument put for- ward by Ruby et al [4] who argued that the closed Fermi sheet (left panel of Fig. 10) contributes to the smaller gap and the open sheet (right panel of Fig. 10) relates to the outer peak. There are many factors which could be responsible for this discrepancy. Perhaps the most ob- vious is the fact that we performed a bulk calculation, whereas their experiment probes the surface states. An accurate calculation of the surface is therefore crucial to shed light on this aspect. FIG. 11. The Fermi surfaces with the colour scale showing the Interestingly, the Fermi velocities as shown in Fig. 11 Fermi velocity of Pb. The labels refer to the gaps identified together with the simple argument developed for Nb in Fig. 9. would suggest the open band to show the larger gap in agreement with the experimental observation. While for the closed Fermi surface sheet the simple argument con- Other theoretical work on Pb by A. Floris at al.[50] es- necting the Fermi velocity to the size of the gap largely tablished the anisotropic electron phonon coupling from holds the relationship is less convincing for the open fully first principles calculations for a bcc crystal struc- sheet. It is interesting to note the large value of the in- ture. Since this is a different crystal structure, it teraction parameter of lead compared to niobium despite complicates the comparison. The extracted degree of the fact that the superconducting gap of lead is smaller anisotropy, measuring the relative difference between the than the gap of niobium (see Table I.). This underlines gaps on the two Fermi surface sheets, was ≈ 20%. For the the importance of dealing with the realistic band struc- fcc crystal considered here, the same anisotropy measure ture and also suggests that the superconducting state of is ≈ 4%. 7 C. Magnesium Diboride The third example is Magnesium diboride, which is well established as a phonon-mediated high temperature superconductor. Its superconductivity is mainly driven via a large E2g mode derived from the Boron atom [54– 56]. For this reason, we model the superconducting state with the interaction on the Mg atoms ΛM g = 0 and the interaction parameter for the Boron atoms as ΛB = 0.288Ry. This parameter is tuned to fit the exper- imental zero temperature gap size [5–7]. In experimental studies, the smaller gap at zero temperature ranges from 1.8meV [6] to approximately 3meV [5, 7], whereas the larger gap is around 7meV. In Fig. 12 we present the density of states within the superconducting state with the lattice constants a = 5.8317 a.u. and c = 6.6594 a.u. in a hexagonal lattice. FIG. 13. Left Panel: The total and atom orbital resolved density of states of MgB2 near the Fermi level. Panels 1- 3: Band gaps in the Γ to K and A to L directions. Panel 1 contributes to the largest gap in the DOS, panel 2 contributes to the middle gap and panel 3 is associated with the smallest gap. The units of the spectral function are arbitrary. third peak vanishes going beyond a thickness of 3 MLs highlighting the importance of out of plane hybridization. This suggests that any result will subtly depend on the out of plane lattice constant which we fixed to the exper- imental value rather using structural relaxation. Relaxed structures would show a smaller lattice constant possibly increasing out of plane hybridization and suppressing the FIG. 12. The electronic density of states of MgB2 in the su- third gap. On the other hand, in experiments impurity perconducting state. Inset shows the density of states around scattering will broaden any gap structures making it dif- the superconducting gap. ficult to resolve a possible third peak. There is very good agreement between the experimen- tal gap sizes and the smallest gap at ∆1 = 3.27meV and V. DISCUSSIONS AND CONCLUSIONS the largest at ∆3 = 7.00meV, which can be clearly iden- tified in Fig. 12. In addition, there is a third peak asso- We have implemented a self-consistent solution of the ciated with a third superconducting gap. Fig. 13 shows BdG equations into the 3D bulk screened KKR formal- three band gaps in some high symmetry directions as- ism. In order to model realistic bulk systems, the BdG sociated with each of the three peaks. As before Fig. 14 equation was solved self-consistently by choosing a sim- extends this analysis to the full Fermi surface. Each sheet ple exchange correlation functional (3) to model both the can be associated with one of the three distinct gaps normal and superconducting order parameters ρ(r) and going from ∆1 = 3.27meV (blue) over ∆2 = 4.77meV χ(r) respectively. The key parameter is the interaction (green) to ∆3 = 7.00meV (red). Λ which was tuned such that the zero temperature gap Comparing this result to the Fermi velocities (Fig. 15) size matched the predicted experimental gaps. the simple relation, as established for Nb, holds to a cer- Calculating the gap self-consistently for Niobium we tain extent. However for the top right panel the Fermi found two superconducting gaps from distinct Fermi sur- velocities vary quite remarkably across the Fermi surface face sheets. We argued that the Fermi velocity of these sheets whereas the gap is relatively constant across each bands at the Fermi level plays a key role in the indi- of the individual sheets (see Fig. 14). cated gap anisotropy. The full anisotropic gap on the The parameter-free calculation from A. Floris et al [51] Fermi surface supported the simple picture connecting accurately predicts both gaps at T = 0K and derives the the gap size to the inverse of the Fermi velocity. The correct transition temperature Tc of MgB2 . This theo- difference between the ∆3 and ∆1 is 0.36meV, and is retical study, among others [6, 20, 55, 57], predicts two comparable to experimental results [8–11]. The most re- distinct gaps in the superconducting state. J. Bekaert et cent tunnelling experiment gives roughly the same gap al [58] identifies a third gap for thin films of MgB2 . This sizes [12] which are ∆001 = 1.20meV, ∆110 = 1.79meV 8 FIG. 14. The Fermi surfaces of MgB2 in the normal state with FIG. 15. The Fermi surfaces with the colour scale showing the the gap size in the superconducting state superimposed as a Fermi velocity of MgB2 . The labels refer to the gaps identified colour scale on top. The labels refer to the panels in Fig. 13 in Fig. 13. and identify the points on the Fermi surface where the gaps in Fig. 13 appear. coupled with magnetic impurities, will make it possible and ∆111 = 1.64meV for the different planes. Similarly to induce a triplet order parameter in these s-wave su- our calculations gave 3 distinct gaps of ∆1 = 1.43meV, perconductors. ∆2 = 1.69meV and ∆3 = 1.79meV. These gaps are very In summary, we showed that using a fully ab initio comparable to the tunneling experiment. model to describe the normal state and a simple ap- proximation for the superconducting exchange correla- For fcc lead we identified three main gaps with a to- tion functional produces gap anisotropy in Nb, Pb and tal separation of 140µeV. This is comparable to tun- MgB2 . This gap anisotropy is of the same order as ex- nelling experiments performed by Ruby et al [4] on single perimental data for each of these systems. One of the crystalline lead surfaces where they found a difference of key features of the formation of the gap anisotropy is the 150µeV. However, our calculations do not support their Fermi velocity at the Fermi surface which roughly corre- Fermi-surface analysis. The future aim will be to inves- lates inversely with the magnitude of the gap. With just tigate the real surfaces based on the method established one free parameter, Λ, we have shown that it is possible here. to model anisotropy in superconductors to a high level of For MgB2 we established gap sizes of similar order as accuracy. In future it will be possible to extend this for- found in the literature [5–7]. Within this framework we malism to model superconductivity including impurities. identified three superconducting gaps which is not unex- pected in a system with three bands crossing the Fermi level and rather varied Fermi velocities. Three gaps have been predicted theoretically before in thin films only [58]. VI. ACKNOWLEDGEMENTS Experiments by Ruby et al probed the influence of magnetic impurities placed on the surface of supercon- The above work was supported by the Centre for Doc- ducting Pb [59, 60]. The presented formalism here can be toral Training in Condensed Matter Physics, funded by readily extended to such impurity systems, enabling the EPSRC EP/L015544/1. B. Újfalussy was supported by study of localised Yu-Shiba-Rusinov states [53, 59, 61]. the Hungarian National Research, Development and In- Furthermore, implementing the fully relativistic BdG novation Office under contract OTKA K115632 and the equations [43] will include spin-orbit interaction. This, BME Nanotechnology FIKP grant (BME FIKP-NAT). 9 G. Csire gratefully acknowledges support from the Euro- ment No. 754510. The authors would like to thank M. pean Union’s Horizon 2020 research and innovation pro- Czerner and Prof Heiliger. In addition, thanks to Ming- gramme under the Marie-Sklodowska-Curie grant agree- Hung Wu and Reena Gupta for many helpful discussions. [1] P. G. Tomlinson and J. P. Carbotte, Phys. Rev. B 13, [28] A. Linscheid, A. Sanna, and E. K. Gross, Phys. Rev. B 4738 (1976). 92, 024506 (2015). [2] B. L. Blackford, Physica 55, 475 (1971). [29] M. Sigrist and K. Ueda, Rev. Mod. Phys. 63, 239 (1991). [3] J. D. Short and J. P. Wolfe, Phys. Rev. Lett. 85, 5198 [30] A. P. Mackenzie, Rev. Mod. Phys. 75, 657 (2003). (2000). [31] J. F. Annett, G. Litak, B. L. Györffy, and K. I. [4] M. Ruby, B. W. Heinrich, J. I. Pascual, and K. J. Franke, Wysokiński, Phys. Rev. B 73, 135401 (2006). Phys. Rev. Lett. 114, 157001 (2015). [32] A. D. Hillier, J. Quintanilla, B. Mazidian, J. F. Annett, [5] P. Szabó, P. Samuely, J. Kačmarčı́k, T. Klein, J. Marcus, and R. Cywinski, Phys. Rev. Lett. 109, 097001 (2012). D. Fruchart, S. Miraglia, C. Marcenat, and A. G. M. [33] B. Mazidian, J. Quintanilla, A. D. Hillier, and J. F. Jansen, Phys. Rev. Lett. 87, 137005 (2001). Annett, Phys. Rev. B 88, 224504 (2013). [6] H. J. Choi, D. Roundy, H. Sun, M. L. Cohen, and S. G. [34] R. P. Singh, A. D. Hillier, B. Mazidian, J. Quintanilla, Louie, Nature 418, 758 (2002). J. F. Annett, D. M. Paul, G. Balakrishnan, and M. R. [7] D. Mou, R. Jiang, V. Taufour, S. L. Bud‘Ko, P. C. Can- Lees, Phys. Rev. Lett. 112, 107002 (2014). field, and A. Kaminski, Phys. Rev. B 91, 214519 (2015). [35] Z. F. Weng, J. L. Zhang, M. Smidman, T. Shang, J. Quin- [8] E. Dobbs and J. Perz, Rev. Mod. Phys. 36, 257 (1964). tanilla, J. F. Annett, M. Nicklas, G. M. Pang, L. Jiao, [9] M. L. A. MacVicar and R. M. Rose, J. Appl. Phys. 39, W. B. Jiang, Y. Chen, F. Steglich, and H. Q. Yuan, 1721 (1968). Phys. Rev. Lett. 117, 027001 (2016). [10] H. Berndt and F. Sernetz, Phys. Lett. A 33, 427 (1970). [36] J. Robbins, J. F. Annett, and M. Gradhand, Phys. Rev. [11] S. Klaumünzer, G. Ischenko, and P. Müller, Zeitschrift B 96, 144503 (2017). für Phys. 268, 189 (1974). [37] A. Kreisel, B. M. Andersen, P. O. Sprau, A. Kostin, J. C. [12] A. Hahn, S. Hofmann, A. Krause, and P. Seidel, Phys. Davis, and P. J. Hirschfeld, Phys. Rev. B 95, 174504 C 296, 103 (1998). (2017). [13] A. Anderson, C. Satterthwaite, and S. Smith, Phys. Rev. [38] P. M. Brydon, D. F. Agterberg, H. Menke, and C. Timm, B. 3, 3762 (1971). Phys. Rev. B 98, 224509 (2018). [14] V. Novotny and P. P. Meincke, J. Low Temp. Phys. 18, [39] V. I. Anisimov, F. Aryasetiawan, and A. I. Lichtenstein, 147 (1975). J. Phys. Condens. Matter 9, 767 (1997). [15] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. [40] G. Csire, B. Újfalussy, J. Cserti, and B. Gyorffy, Phys. Rev. 108, 1175 (1957). Rev. B 91, 165142 (2015). [16] G. M. Eliashberg, Sov. Phys. JETP 11, 696 (1960). [41] G. Csire, J. Cserti, and B. Újfalussy, J. Phys. Condens. [17] J. P. Carbotte, Rev. Mod. Phys. 62, 1027 (1990). Matter 28, 495701 (2016). [18] M. Lüders, M. A. Marques, N. N. Lathiotakis, A. Floris, [42] G. Csire, J. Cserti, I. Tütto, and B. Újfalussy, Phys. G. Profeta, L. Fast, A. Continenza, S. Massidda, and Rev. B 94, 104511 (2016). E. K. Gross, Phys. Rev. B 72, 024545 (2005). [43] G. Csire, A. Deák, B. Nyári, H. Ebert, J. F. Annett, and [19] M. A. Marques, M. Lüders, N. N. Lathiotakis, G. Pro- B. Újfalussy, Phys. Rev. B 97, 024514 (2018). feta, A. Floris, L. Fast, A. Continenza, E. K. Gross, and [44] G. Csire, B. Ujfalussy, and J. F. Annett, Eur. Phys. J. S. Massidda, Phys. Rev. B 72, 024546 (2005). B 91, 217 (2018). [20] E. R. Margine and F. Giustino, Phys. Rev. B 87, 024505 [45] M. Gradhand, M. Czerner, D. V. Fedorov, P. Zahn, B. Y. (2013). Yavorsky, L. Szunyogh, and I. Mertig, Phys. Rev. B 80, [21] J. A. Flores-Livas, A. Sanna, and E. K. Gross, Eur. Phys. 224413 (2009). J. B 89, 63 (2016). [46] D. Sébilleau, K. Hatada, and H. Ebert, Multiple Scat- [22] A. Sanna, J. A. Flores-Livas, A. Davydov, G. Profeta, tering Theory for Spectroscopies (Springer Proceedings K. Dewhurst, S. Sharma, and E. K. Gross, J. Phys. Soc. in Physics, 2018) p. 381. Japan 87, 041012 (2018). [47] L. N. Oliveira, E. K. U. Gross, and W. Kohn, Phys. Rev. [23] F. Schrodi, A. Aperis, and P. M. Oppeneer, Phys. Rev. Lett. 60, 2430 (1988). B 98, 094509 (2018). [48] M. B. Suvasani, W. M. Temmerman, and B. Gyorffy, [24] M. N. Gastiasoro, A. V. Chubukov, and R. M. Fernan- Phys. Rev. B 48, 1202 (1993). des, Phys. Rev. B 99, 094524 (2019). [49] S. Savrasov and D. Savrasov, Phys. Rev. B 54, 16487 [25] F. Essenberger, A. Sanna, P. Buczek, A. Ernst, L. San- (1996). dratskii, and E. K. Gross, Phys. Rev. B 94, 014503 [50] A. Floris, A. Sanna, S. Massidda, and E. K. Gross, Phys. (2016). Rev. B 75, 054508 (2007). [26] F. Essenberger, A. Sanna, A. Linscheid, F. Tandetzky, [51] a. Floris, G. Profeta, N. N. Lathiotakis, M. Lüders, G. Profeta, P. Cudazzo, and E. K. Gross, Phys. Rev. B M. a. L. Marques, C. Franchini, E. K. U. Gross, A. Con- 90, 214504 (2014). tinenza, and S. Massidda, Phys. Rev. Lett. 94, 037004 [27] A. Linscheid, A. Sanna, F. Essenberger, and E. K. Gross, (2005). Phys. Rev. B 92, 024505 (2015). [52] J. Ketterson and S. N. Song, Superconductivity (Cam- bridge University Press, 1999) pp. 195–254. 10 [53] M. Ruby, F. Pientka, Y. Peng, F. von Oppen, B. W. [57] A. Y. Liu, I. I. Mazin, and J. Kortus, Phys. Rev. Lett. Heinrich, and K. J. Franke, Phys. Rev. Lett. 115, 197204 87, 87005 (2001). (2015). [58] J. Bekaert, A. Aperis, B. Partoens, P. M. Oppeneer, and [54] S. L. Bud’ko, G. Lapertot, C. Petrovic, C. E. Cunning- M. V. Milošević, Phys. Rev. B 96, 094510 (2017). ham, N. Anderson, and P. C. Canfield, Phys. Rev. Lett. [59] M. Ruby, Y. Peng, F. von Oppen, B. W. Heinrich, and 86, 1877 (2001). K. J. Franke, Phys. Rev. Lett. 117, 186801 (2016). [55] I. I. Mazin and V. P. Antropov, Phys. C 385, 49 (2003). [60] M. Ruby, B. W. Heinrich, Y. Peng, F. von Oppen, and [56] A. Shukla, M. Calandra, M. D’Astuto, M. Lazzeri, K. J. Franke, Phys. Rev. Lett. 120, 156803 (2018). F. Mauri, C. Bellin, M. Krisch, J. Karpinski, S. M. Kaza- [61] B. W. Heinrich, J. I. Pascual, and K. J. Franke, Prog. kov, J. Jun, D. Daghero, and K. Parlinski, Phys. Rev. Surf. Sci. 93, 1 (2018). Lett. 90, 095506 (2003).