Gap Anisotropy in Multiband Superconductors Based on Multiple Scattering Theory.

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 MgB 2 . In these materials multiple distinct peaks at the gap in the density of states were observed, showing signiﬁcant 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
The structure of the superconducting gap in s-wave phonon mediated superconductors may show a surprisingly complex structure. Not only can it show a large degree of anisotropy on the Fermi surface of a given material but several known examples exhibit clear signatures of two superconducting gaps. The most notable of cases are Pb [1][2][3][4] and MgB 2 [5][6][7] and to a lesser extent Nb [8][9][10][11][12] with conflicting reports [13,14]. While it is understood that the gap anisotropy arises in these systems from multiple Fermi surface sheets, not many phonon mediated s-wave superconductors showing multiple gaps have been identified.
In the original BCS theory [15] only a single spherical band was considered and so it was impossible to obtain any gap anisotropy. In addition, it was limited to the weak coupling regime, which later was rectified by the Eliashberg theory [16,17] and more recently by full DFT approaches to electron phonon coupling driven superconductivity [18,19]. Although gap anisotropy analysis is present in such codes, much of the DFT based ab initio work focused on other aspects. These aspects being; expanding the description of s-wave superconductivity with emphasis on the correct description of the driving mechanism [20][21][22][23][24], or on the extension to unconventional pairing pushing the boundaries in our treatment of iron based [25], spin [26] and magnetic effects [27,28].
Beyond the ab initio work on superconductivity, there is also extensive work on using parametrised models to describe superconductivity [29][30][31][32][33][34][35][36][37][38]. These models are a powerful tool for describing unconventional superconductors and some of the basic principles of standard s-wave superconductors. However these models have the drawback that they require the normal state to be parametrised from either experimental data or a density * t.saunderson@bristol.ac.uk functional theory (DFT) calculation. This often results in over-simplification of the normal state band structure in order to construct an efficient model to describe the superconducting state. While this usually leads to deep understanding of aspects of the superconductivity, it might obscure the importance of the complexity and orbital hybridisation of the underlying electronic structure.
In this work we aim to follow a route distinct to those two main directions. We aim to describe the full complexity of the normal state electronic structure whilst treating the superconducting pairing interaction within a simplified model. This idea is very similar to the LDA+U method where the Hubbard U interaction is also treated as a tunable parameter [39]. This has proven to be very successful in modelling strongly correlated materials without compromising on the full electronic structure within a strongly correlated state.
In analogy to this, our method leads to a full electronic structure within the superconducting state giving access to the full gap anisotropy in multi sheeted s-wave superconductors. We build upon previous work from some of the authors on the implementation of the scalar relativistic [40][41][42] and fully relativistic [43,44] BdG equations in a layered KKR DFT code. The advantage of this approach is the access to the full normal state electronic structure from first principles without using simplified models. By just considering s-wave pairing and one effective pairing parameter, we show quantitatively the complexity of the superconducting gap including its full anisotropy even in simple elemental crystals.
The work is structured as follows. First, we will outline the basic theoretical background of the method including the differences in implementation to the earlier work [40][41][42]. This will be followed by numerical tests to determine the robustness and accuracy of the theoretical 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 some simple arguments on how the normal state proper-ties drive the observed anisotropies. In all cases we will relate our work to experimental observations.

II. METHOD
This section will describe the implementation of superconductivity into the existing KKR code [45,46]. The theory for the single site solver and multiple scattering terms have been introduced by G. Csire et al [40]. All equations are given in Rydberg atomic units.
The density functional theory for superconductors was initially presented by L. N. Oliveira, E. K. U. Gross and W. Kohn [47], who introduced the effective pairing potential ∆ ef f (r, r ), describing the superconducting state in addition to the conventional Kohn-Sham potential V ef f (r). These two potentials are defined as where χ(r, r ) is the anomalous density and E xc [ρ, χ] is the exchange correlation functional for a superconductor. The exchange correlation functional was later approximated [48] as where E 0 xc [ρ] is the normal state exchange-correlation functional, and Λ[ρ, χ](r 1 , r 1 , r 2 , r 2 ) is the pairing kernel. Similarly to ref. [48], in the KKR framework we consider the atomic sphere approximation (ASA). The kernel is then approximated to Λ[ρ, χ](r 1 , r 1 , r 2 , r 2 ) = Λδ(r 1 − r 1 )δ(r 1 − r 2 )δ(r 1 − r 2 ), (4) inside the ASA spheres and outside, where Λ is called the interaction parameter. A further conventional simplification is for ∆ ef f (r, r ) and χ(r, r ) to be local, resulting in Equation (2) therefore becomes ∆ ef f (r) = Λχ(r).
When the spherical symmetry (ASA) and local approximation for the pairing potential are combined, we restrict ourselves to s-wave superconductivity. Any kind of non-s wave superconductivity would need a non-spherical pairing potential, coupling between different orbital channels or non-locality in the pairing potential. In this report Λ is 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 this are discussed in section III. The KKR method exploits a local atomic basis that uses multiple scattering theory, which gives access to the full Green's function of the system. For details on the implementation we refer to G. Csire et al [40]. Here we restrict the discussion to the essential components to highlight differences in the implementation and to shed light onto the results. The Green's function for the system is defined aŝ whereĤ BdG (r) = r|Ĥ BdG |r and Here, µ is the chemical potential and z = + iδ. The densities ρ(r) and χ(r) can be calculated by taking the trace of different components of G BdG (z, r, r ), where the limit is taken such that δ → 0 + . We use equations (12) and (13) to find new expressions for V ef f (r) and ∆ ef f (r) using equations (1) and (8). From here a newĤ BdG (r) is constructed to calculate a new Green's functionĜ BdG ( , r, r ), and thus self consistency can be achieved.
The ASA approximation sets a boundary to each atomic site, i, called the ASA radius r ASA i . This approximation implies that the potentials V ef f (r) and ∆ ef f (r) can be written in sums and ensures that V i (r) = 0 and ∆ i (r) = 0 if |r| ≥ r ASA i . The analogue of equation (8) becomes, The resulting code is therefore able to perform calculations which relax both charge and anomalous densities 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 resulting gaps within our framework. We define the average gap [48] where V W 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 actual gap in the density of states. Fig. 1 shows the results of the convergence tests for non-relativistic (NR) Niobium. For low Λ, all calculations with 30 energy points are suppressed relative to calculations with a higher energy meshes. At high Λ the solutions with the same sized k-mesh group together. We decided that 50 energy points and 200 × 10 3 k-points is the best trade off between accuracy and computational cost.
FIG. 1. The convergence test performed on a non-relativistic calculation of Niobium. The parameters are the energy and k-mesh. The k-mesh is the actual number of k-points used to calculate the energy points closest to the superconducting gap. Plot symbols denote 2 × 10 4 (circle), 2 × 10 5 (cross) and 5 × 10 5 (square) k-points, colours denote 30 (red), 50 (green) and 60 (blue) energy points.
For the remainder of the report, we perform calculations using the scalar-relativistic (SR) BdG equations. A comparison between NR and SR results is shown in Fig. 2 for Nb, V and Cu. There is no clear trend in terms of the gap size going from NR to SR. For V and Nb the gap is reduced, for Cu it is increased. The influence becomes more dominant with larger atomic number as expected when including scalar relativistic corrections. Within this framework, Λ is the only free parameter. In general this parameter is tuned such that the resulting gap matches the experimental zero temperature gap size. The complication, especially for anisotropic gaps or multigap systems, is what definitions are used experimentally and theoretically to establish the gap size. In the following we discuss how in practice we constrain our free parameter Λ. The first definition is using the average gap from equation (17). This quantity, however, is just the integrated anomalous density of states χ(r) multiplied 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 T c 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 quantity is distinct to the average∆ r calculated from the real space integral over ∆ ef f (r) (see Table I, column 5). For Nb and Pb,∆ r is smaller than∆ k but for MgB 2 it is larger.
Finally, we would like to summarize the numerical parameters used throughout the papers. After convergence of the potentials is reached, the imaginary part of the energy, δ = 2µRy, 0.2µRy, 10µRy is used for Nb, Pb and MgB 2 respectively. These are the parameters for the DOS and Bloch spectral function calculations when focussing on the gap structure as well as for ∆(k). The number of k-points in the irreducible Brillouin zone for DOS calculations are 1 × 10 7 for Nb and Pb and 2.7 × 10 6 for MgB 2 .  (17). ∆ exp are average gaps from experiments [4,7,12,49],∆ k are average gaps from theoretical ∆(k) integrations [49][50][51].

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 superconducting 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, ∆ 2 = 1.69meV and ∆ 3 = 1.79meV, which is in reasonable quantitative agreement to the experiment. In comparison to other literature the gap anisotropy is also well matched [8][9][10][11]. The lattice constant of a = 6.23685a.u. was used for the bcc structure of Nb. The lower panel of Fig. 4 shows the Bloch spectral function for a large energy window but in the top panel we focus on the region around the Fermi level along the Γ to N direction with a smaller energy window (−50meV to 50meV). On that scale the opening of the superconducting gap is just about resolved. In order to investigate the gap and the associated anisotropy as highlighted in  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 contributing 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 points of the Fermi surface associated with panels 1-3 from Fig. 5 are labelled. Evidently panels 1 and 2 have similar velocities at the point of crossing, whereas panel 3 has a velocity which is noticeably larger. Since the DOS is inversely proportional to the Fermi velocity, and within the BCS theory [15,52] the gap scales with the DOS we can infer a gap anisotropy which is related to the distribution of Fermi velocities on the Fermi surface sheets. In Fig. 7 we extend this analysis to the full Fermi surface and in comparing it to Fig. 6 a strong correlation between v F (k) and ∆(k) is observed.  Fig. 5.

B. Lead
Lead often has disordered or amorphous crystal structures. Experimentally, this has been avoided in STM experiments by M. Ruby et al [4]. In this experiment the authors could identify two distinct peaks of the superconducting 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 future work will focus on describing the surface explicitly as measured in the experiment. For the bulk material we assume that the gap anisotropy in Pb should at least be of similar order as found in the STM experiments. For an interaction parameter of Λ = 0.351Ry the over 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 comparable to the 150µeV identified by Ruby et al [53]. In order to resolve such small separation a fine mesh of 1 × 10 7 k-points in the irreducible part of the Brillouin zone was required. Following the same process as for Nb the distinct gap sizes can be traced to different points in k-space as shown in Fig. 9. However, to get the full picture 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 associated 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 forward 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 obvious 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.
Interestingly, the Fermi velocities as shown in Fig. 11 together with the simple argument developed for Nb 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 connecting the Fermi velocity to the size of the gap largely holds the relationship is less convincing for the open sheet. It is interesting to note the large value of the interaction parameter of lead compared to niobium despite the fact that the superconducting gap of lead is smaller than the gap of niobium (see Table I.). This underlines the importance of dealing with the realistic band structure and also suggests that the superconducting state of  Fig. 9 and identify the points on the Fermi surface where the gaps in Fig. 9 appear. lead falls far from the weak coupling limit. Since BCS theory is valid in the weak coupling limit, therefore the relation between the Fermi velocity and the superconducting gap is expected to be less robust. This aspect deserves further investigation and in particular the influence of surfaces and spin-orbit coupling has to be addressed in future work. Other theoretical work on Pb by A. Floris at al. [50] established the anisotropic electron phonon coupling from fully first principles calculations for a bcc crystal structure. Since this is a different crystal structure, it complicates the comparison. The extracted degree of anisotropy, measuring the relative difference between the gaps on the two Fermi surface sheets, was ≈ 20%. For the fcc crystal considered here, the same anisotropy measure is ≈ 4%.

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 E 2g mode derived from the Boron atom [54][55][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 experimental zero temperature gap size [5][6][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. There is very good agreement between the experimental gap sizes and the smallest gap at ∆ 1 = 3.27meV and the largest at ∆ 3 = 7.00meV, which can be clearly identified in Fig. 12. In addition, there is a third peak associated with a third superconducting gap. Fig. 13 shows three band gaps in some high symmetry directions associated with each of the three peaks. As before Fig. 14 extends this analysis to the full Fermi surface. Each sheet can be associated with one of the three distinct gaps going from ∆ 1 = 3.27meV (blue) over ∆ 2 = 4.77meV (green) to ∆ 3 = 7.00meV (red).
Comparing this result to the Fermi velocities (Fig. 15) the simple relation, as established for Nb, holds to a certain extent. However for the top right panel the Fermi velocities vary quite remarkably across the Fermi surface sheets whereas the gap is relatively constant across each of the individual sheets (see Fig. 14).
The parameter-free calculation from A. Floris et al [51] accurately predicts both gaps at T = 0K and derives the correct transition temperature T c of MgB 2 . This theoretical study, among others [6,20,55,57], predicts two distinct gaps in the superconducting state. J. Bekaert et al [58] identifies a third gap for thin films of MgB 2 . This 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 experimental value rather using structural relaxation. Relaxed structures would show a smaller lattice constant possibly increasing out of plane hybridization and suppressing the third gap. On the other hand, in experiments impurity scattering will broaden any gap structures making it difficult to resolve a possible third peak.

V. DISCUSSIONS AND CONCLUSIONS
We have implemented a self-consistent solution of the BdG equations into the 3D bulk screened KKR formalism. In order to model realistic bulk systems, the BdG equation was solved self-consistently by choosing a simple exchange correlation functional (3) to model both the normal and superconducting order parameters ρ(r) and χ(r) respectively. The key parameter is the interaction Λ which was tuned such that the zero temperature gap size matched the predicted experimental gaps.
Calculating the gap self-consistently for Niobium we found two superconducting gaps from distinct Fermi surface sheets. We argued that the Fermi velocity of these bands at the Fermi level plays a key role in the indicated gap anisotropy. The full anisotropic gap on the Fermi surface supported the simple picture connecting the gap size to the inverse of the Fermi velocity. The difference between the ∆ 3 and ∆ 1 is 0.36meV, and is comparable to experimental results [8][9][10][11]. The most recent tunnelling experiment gives roughly the same gap sizes [12] which are ∆ 001 = 1.20meV, ∆ 110 = 1.79meV FIG. 14. The Fermi surfaces of MgB2 in the normal state with the gap size in the superconducting state superimposed as a colour scale on top. The labels refer to the panels in Fig. 13 and identify the points on the Fermi surface where the gaps in Fig. 13 appear. and ∆ 111 = 1.64meV for the different planes. Similarly our calculations gave 3 distinct gaps of ∆ 1 = 1.43meV, ∆ 2 = 1.69meV and ∆ 3 = 1.79meV. These gaps are very comparable to the tunneling experiment.
For fcc lead we identified three main gaps with a total separation of 140µeV. This is comparable to tunnelling experiments performed by Ruby et al [4] on single crystalline lead surfaces where they found a difference of 150µeV. However, our calculations do not support their Fermi-surface analysis. The future aim will be to investigate the real surfaces based on the method established here.
For MgB 2 we established gap sizes of similar order as found in the literature [5][6][7]. Within this framework we identified three superconducting gaps which is not unexpected 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].
Experiments by Ruby et al probed the influence of magnetic impurities placed on the surface of superconducting Pb [59,60]. The presented formalism here can be readily extended to such impurity systems, enabling the study of localised Yu-Shiba-Rusinov states [53,59,61]. Furthermore, implementing the fully relativistic BdG equations [43] will include spin-orbit interaction. This, coupled with magnetic impurities, will make it possible to induce a triplet order parameter in these s-wave superconductors.
In summary, we showed that using a fully ab initio model to describe the normal state and a simple approximation for the superconducting exchange correlation functional produces gap anisotropy in Nb, Pb and MgB 2 . This gap anisotropy is of the same order as experimental data for each of these systems. One of the key features of the formation of the gap anisotropy is the Fermi velocity at the Fermi surface which roughly correlates inversely with the magnitude of the gap. With just one free parameter, Λ, we have shown that it is possible to model anisotropy in superconductors to a high level of accuracy. In future it will be possible to extend this formalism to model superconductivity including impurities.

VI. ACKNOWLEDGEMENTS
The above work was supported by the Centre for Doctoral Training in Condensed Matter Physics, funded by EPSRC EP/L015544/1. B.Újfalussy was supported by the Hungarian National Research, Development and Innovation Office under contract OTKA K115632 and the BME Nanotechnology FIKP grant (BME FIKP-NAT).
G. Csire gratefully acknowledges support from the European Union's Horizon 2020 research and innovation programme under the Marie-Sklodowska-Curie grant agree-ment No. 754510. The authors would like to thank M. Czerner and Prof Heiliger. In addition, thanks to Ming-Hung Wu and Reena Gupta for many helpful discussions.