Tunneling and electronic structure of the two-gap superconductor MgB

J. A. Silva-Guillén,1,2 Y. Noat,3 T. Cren,3 W. Sacks,4 E. Canadell,5 and P. Ordejón1,2 1ICN2–Institut Català de Nanociència i Nanotecnologia, Campus UAB, 08193 Bellaterra, Spain 2CSIC–Consejo Superior de Investigaciones Cientı́ficas, Campus UAB, 08193 Bellaterra, Spain 3Sorbonne Universités, UPMC Université Paris 6, CNRS,UMR 7588, Institut des Nanosciences de Paris, 4 Pl. Jussieu 75005 Paris, France 4Sorbonne Universités, UPMC Université Paris 6, CNRS, Institut de Minéralogie, de Physique des Matériaux et de Cosmochimie, 4 Pl. Jussieu 75005 Paris, France 5Institut de Ciència de Materials de Barcelona (CSIC), Campus UAB, 08193 Bellaterra, Spain (Received 28 May 2015; published 20 August 2015)


I. INTRODUCTION
MgB 2 has been known for a long time [1,2], but interest in this material was boosted by the discovery of superconductivity with a high critical temperature (T c = 39 K). [3] It exhibits the AlB 2 -type crystal structure [4], where hexagonal layers of graphenelike boron atoms alternate with hexagonal layers of magnesium atoms sitting on top of the center of the boron hexagons (see Fig. 1).The simplicity of the structure made possible very detailed theoretical studies.An and Pickett [5] proposed that superconductivity originates in the boron (p x ,p y ) bands and Liu et al. [6] suggested the possibility of two-gap superconductivity for this compound.
The two-band superconductivity scenario for MgB 2 soon received support from different experimental studies using scanning tunneling spectroscopy (STS) [7], point-contact spectroscopy [8][9][10], specific heat measurements [11], and Raman spectroscopy [12].Iavarone et al. [13] reported tunneling studies on both c-axis-oriented films and compact pellets and provided evidence for directional tunneling with respect to the crystallographic orientation of the grains in the pellets.It is now well established that there are two distinct superconducting (SC) energy gaps, S = 2.3 meV and L = 7.1 meV [13,14].As pointed out by Koshelev and Golubov [15], this peculiarity leads to two different length scales, revealed experimentally in the vortex structure [16].
Although there have been several first-principles theoretical studies [5,6,17,18] of the electronic structure of MgB 2 , to our knowledge there is only one study of the tunneling images by Li et al. [19].As noted by Schmidt et al. [8], the calculated contributions of the two different types of bands to the bulk density of states at the Fermi level are insufficient to account for the tunneling experiments.In addition, tunneling studies on other superconductors such as 2H -NbSe 2 [20] or Ba 8 Si 46 [21] report a single SC gap measured by STS, while two gaps are rather observed in specific heat measurements, strongly suggesting the existence of tunneling selectivity at the surface [22].Therefore, further work is needed in order to understand properly the tunneling spectra of multigap superconductors, which explicitly takes into account the shape and symmetry of the electron wave functions at the surface and allows to calculate the different contributions of each band as a function of the distance from the surface.
In this work, we report tunneling measurements of MgB 2 measured at low temperatures in the superconducting state.We fit the experimental spectra with the McMillan equations for the proximity effect in reciprocal space [23], which describes properly the SC state of two-gap superconductors [22,24].We also present a first-principles density functional theory (DFT) [25,26] study of the electronic structure in the normal state, from which we identify the character of the electronic states associated with the two-gap superconductivity.These calculations allow us to explain the ratio of the interband coupling parameters deduced from the fit of various STS experiments of the literature, including ours, with the McMillan equations.Finally, we calculate using DFT the tunneling selectivity arising from the wave function at the MgB 2 surface, allowing a full understanding of the tunneling spectra.This combined experimental-theoretical approach is also useful to understand other multigap superconductors such as 2H -NbSe 2 , Ba 8 Si 46 , the FeSe-based superconductors, and so on.

II. EXPERIMENTAL TUNNELING SPECTRA
As mentioned above, MgB 2 is a multigap superconductor with two gaps, S (small gap) and L (large gap).Several models have been used to describe two-gap superconductivity, which imply either a pair coupling between bands [27] or, alternatively, an interband quasiparticle coupling [28].In the Suhl-Matthias-Walker (SMW) model, the SC density of states (DOS) is simply given by a sum of BCS-like DOS arising from each band, which does not fit properly the experimental spectra.In the second case, the two-gap superconductor is described by the McMillan equations and superconductivity is induced from one band to the other by the proximity effect in reciprocal space.This proximity effect is mediated by quasiparticle scattering from one band to the other.As stressed by Schmidt et al. [24], only the latter can describe properly the peculiar shape of the excitation spectrum deduced from STM or point contact measurements.
In Fig. 2, we show two typical tunneling spectra in superconductor-insulator-superconductor (SIS) geometry with a SC MgB 2 tip and a MgB 2 surface [Fig.2(a)] or V 3 Si surface [Fig.2(b)] measured using our home-built STM/STS setup at low temperature (T = 5 K).The shape of the MgB 2 −MgB 2 tunneling spectrum [Fig.2(a)] is very similar to the one reported by Schmidt et al. [24].It deviates strongly from what is expected from a single gap s-wave superconductor.On the other hand, it can be described satisfactorily with the two-gap McMillan model (see the fit in Fig. 2).The parameters deduced from the McMillan fits are very similar in both cases (MgB 2 or V 3 Si surface).We include an additional broadening parameter attributed to either a small gap anisotropy and/or a slight oxidation of the V 3 Si surface.The latter are taken into account by using a small imaginary part to the gap − iσ , with σ = 0.5 meV [29].The normalized tunneling DOS probed by STS for a twogap superconductor can be written, according to the McMillan model, as: where the S and L indexes indicate the "small" and "large" gaps, respectively, and the normal-state bands from which they originate.N n (E F ) is the normal-state density of states at the Fermi level.T S,L = α S,L N S,L n (E F )/N n (E F ) account for the partial DOS at the Fermi energy for each of the bands N S,L n (E F ) as well as the averaged tunneling probability α S,L (which result from the band structure, the symmetry of the bands, and also the tip position and electronic structure).
S,L (E) are the energy-dependent gaps, which are obtained from the self-consistent equations [23]: (2) In addition, in this model the ratio of the interband coupling parameters is related to the ratio of the DOS at the Fermi level for the two bands: For MgB 2 , from the McMillan fit to the experimental STS spectra, we find a S / L ratio close to 1, which is in agreement with the data inferred from specific heat and penetration depth measurements (see Ref. [30] and references therein).It is also found that T S is close to 1, while T L is nearly zero, meaning that the contribution of the small gap states strongly dominates in the tunneling process.The latter is true for tunneling measurements along the c-axis for a typical tip-surface distance.A significant contribution of the large gap band can arise for a different crystal orientation [see Fig. 3(c)] or near the contact regime [31].In order to test the robustness of these results, we compare our values with those obtained by fitting data from the literature [13,32,33].The results are summarized in Table I.
In Fig. 3(a), we show a typical low-temperature spectrum obtained by Iavarone et al. [13] for a c-axis-oriented MgB 2 sample and the corresponding McMillan fit.One can note that the model describes well the experimental spectrum.As for the SIS measurement, the shape of the tunneling spectrum clearly deviates from what is expected from a single gap s-wave superconductor.For this orientation, the tunneling selectivity coefficient towards the small gap is nearly 1.
Some spectra [see Fig. 3(a)] exhibit a small contribution of the small gap band, probably as a result of a local defect breaking the planar symmetry.This gives rise to a hump around 7 meV and to a small change in the quasiparticle band coupling.On the contrary, there is no hump [Fig.3(a)] when the spectrum reveals the full contribution of the small gap (planar symmetry).On an a-oriented crystallite [Fig.3 there is a more significant contribution of the large gap band (T L ≈ 0.4 − 0.5), the tunneling selectivity coefficients being modified as a result of the crystal orientation.Although the spectral weight deviates slightly above 10 meV, the overall double-gap features are clearly matched.However, as a result of the quasiparticle interband coupling, the small gap S (E) depends on the large gap L (E).This is not the case in the Suhl's model.Variations in T S can be attributed to the random orientation of the MgB 2 grain attached to the PtIr tip.Similarly, as for the SIS junctions, the DOS contributions are the same for the small and large gaps.
As shown in Table I, the parameters are very similar for all experiments.Therefore, we obtain two important conclusions from these data: (i) the tunneling process is dominated by the small gap and (ii) the ratio between the interband coupling parameters is close to 1.In the following we will explain these observations on the basis of DFT calculations.

III. ELECTRONIC STRUCTURE OF MgB 2
The electronic structure of MgB 2 has been discussed a number of times in the literature [5,6,17,18,34].Here we will touch upon those aspects that are relevant to understand the signatures of multigap superconductivity in tunneling spectroscopy.For this purpose, we compute the electronic structure of MgB 2 using a DFT method that uses numerical atomic orbitals as basis sets, implemented in the SIESTA code [35,36].We have used the local density approximation (LDA) to DFT and, in particular, the functional of Ceperly-Alder [37].Only the valence electrons are considered in the calculation, with the core being replaced by norm-conserving scalar relativistic pseudopotentials [38] factorized in the Kleinman-Bylander form [39].The nonlinear core-valence exchange-correlation scheme [40] was used for all elements.We have used a split-valence double-ζ basis set, including polarization functions [41].The energy cutoff of the real-space integration mesh was 300 Ry.The experimental crystal structure was used for the bulk calculations.For the simulation of the STM images, we use a B-terminated symmetric slab, containing 9 and 10 layers of Mg and B, respectively.The Brillouin zone was sampled using 30 × 30 × 30 and 30 × 30 × 1 k-points in the Monkhorst-Pack scheme [42] for the bulk and slab calculations, respectively.
The calculated band structure is shown in Fig. 4.There are three partially filled bands.Two of them are associated with the σ bonds within the boron layers, and they have mostly weight on the boron p x and p y orbitals (bands shown with blue circles in Fig. 4).Therefore, they have a strong two-dimensional (2D) character because of the very small interaction with the orbitals of the Mg atoms intercalated between successive B planes (see the small dispersion along the to A line).The Fermi surface sheets corresponding to these bands are warped cylinders along the c * -axis, centered around the point (see Fig. 5).In contrast, the remaining partially filled band (shown with red circles in Fig. 4) is built from the boron p z orbitals and exhibits dispersion along both the plane of the boron layers (because of their π -type interactions along the boron layers) and the interlayer direction (because of the good overlap between the out-of-plane pointing B p z and Mg orbitals).Consequently, the corresponding Fermi surface sheets are a more complex 3D network.We will refer to these two sets of bands and their associated Fermi surface sheets as σ and π , respectively.Despite the very different shapes and sizes of the Fermi surfaces, the contribution from σ cylinders and the π 3D network to the DOS at the Fermi level are very similar: 45% and 55%, respectively (see Table II), because the Fermi velocity of the cylinders is very low along the c * -direction.The ratio of the partial DOS of the two sets of bands is, therefore, N σ (E F )/N π (E F ) ≈ 0.84.If we identify the σ and π states with the large and small superconducting gaps (see Sec. IV), we obtain a value which is in good agreement with the ratio extracted from the experimental tunneling spectra shown in Sec.II, which was close to 1.

IV. TUNNELING IMAGES AND SELECTIVITY
We now focus on the experimental fact that the tunneling process is dominated by the contribution of the small gap (see Sec. II).To understand the origin of this observation, which lies in the different tunneling sensitivity of the σ and π states, we have calculated the STM current in the normal state using a boron-terminated slab.Note that, because of the weak contribution of the Mg atoms to the DOS at the Fermi level [see Fig. 4(b)], even in the case of a surface partially or completely covered by Mg atoms, these are only visible at very low surface-tip distances.As shown by Li et al. [19], Mg-terminated surfaces should exhibit the same type of STM images as the B-terminated ones for any experimentally reasonable height in STS experiments.
Our calculations were done using the Tersoff-Hamann approximation [43], where the current at a given tip position is proportional to the local density of states at that point, integrated over the energy window defined by the tip-surface potential difference (which we take as ±0.1 eV).The images correspond to iso-DOS plots, showing the map of heights that produce a constant tip-surface current.In our images, shown in Fig. 6, we separate the contribution of the σ and π bands.As expected, in the STM images generated from the σ bands, the brightest positions are the center of the bonds between surface B atoms.On the other hand, for the π states, the brightest positions are on top of the surface B atoms.The full image (not shown) would correspond to the sum of the two contributions.Experimentally, STM images with atomic resolution are unfortunately difficult to achieve with this material.Images of the MgB 2 surface are presented in Ref. [44], showing at least a qualitative agreement with our calculation and those of Ref. [19].
To address the tunneling sensitivity, we now compare the values of the contributions to the current due to tunneling through the σ and π states.Figure 7 shows the ratio of these currents I π /I σ , as a function of the tip-surface distance, for several locations of the tip on the surface plane.It is clear that, for any reasonable tip-surface distance, the contribution to the current from the π states is about three orders of magnitude larger than that from the σ states, independently of the position of the tip along the surface.The tunneling selectivity therefore indicates that the tunneling images will be dominated by the electronic states associated with the π bands.In the superconducting regime, therefore, we expect that the STS spectra will only be sensitive to one of the gap components.As we described in Sec.II, this is actually what one observes in the experiment, where only the small gap component is observed in the tunneling spectra.This allows us to conclude that the small superconducting gap is associated with the Fermi surface sheets from the π states, whereas the large gap is associated with the σ states.
The stronger contribution of the π states to the STM current, compared to that of the σ states, is partly explained by the  intrinsic difference of the decay of the B p z and (p x ,p y ) orbitals along the c-axis (perpendicular to the surface), as shown in Fig. 8.The decay into vacuum of these atomic orbitals is governed by the product of the radial part (which is common for both sets of orbitals in the free atom) and the angular part (spherical harmonics).The radial part decays exponentially into the vacuum for large distances with the same rate, but the angular part introduces a further decay along the c axis, which is much stronger for the (p x ,p y ) orbitals than for p z .
The symmetry of the lattice and the wave functions at the Fermi surface further enhance this different decay.The wave functions of the Fermi surface cylinders have a pseudo etype symmetry.Because of their nodal properties, these e-type functions have a nil value along the threefold symmetry axis parallel to the c-direction and going through one boron atom, as well as the axis passing through the center of the boron hexagons.In contrast, no such symmetry restrictions apply to the wave functions based on the boron p z orbitals.The implications for the tunneling intensities can be understood with the following argument, based on the decay of the wave function into vacuum at large distances from the surface typical of tunneling experiments.Far from the surface, where the potential is roughly the vacuum potential, the wave functions can be expressed analytically as: where C k (G) is the Fourier component of the wave function at a reference plane (taken as z = 0), G = (G x ,G y ) the surface reciprocal lattice vector, and r = (x,y) the in-plane position.Each Fourier component decays into vacuum with a decay factor which depends on k and G, given by where κ is the standard inverse decay length determined by the work function φ (i.e., κ = √ 2mφ/ ).For the σ bands of MgB 2 , the expansion in Eq. ( 4) will contain a large weight of Fourier components with large G vectors due to their nodal structure.These components will decay faster into vacuum, according to Eq. ( 5).However, the π wave functions, as they are smoother and do not have such nodal structure, will have larger components for small G vectors and therefore a slower decay into vacuum.

V. CONCLUSIONS
A combined experimental (SIS tunneling spectra) and theoretical (DFT) study of the two-gap superconductor MgB 2 clearly shows that the small gap is associated with the band mostly based on the boron p z orbitals leading to the 3D component of the Fermi surface.This band strongly dominates the tunneling current along the c-direction, so it is the major contributor to the STM/STS experiments at this surface.Both the directional shape of the boron p z orbitals and the symmetry properties of the lattice and the wave functions determine this selectivity.The experimental observation of only the small gap states in STS spectra allows us to identify these with the π bands associated with the B p z orbitals.In addition, the relative value of the theoretical density of states at the Fermi level for the π and and σ bands are in good agreement with those inferred from the experimental spectra for the small and large gap components.The latter two gaps follow directly using the self-consistent McMillan equations.
Our approach is thus useful to unravel the nature of the different gaps in multigap superconductors and to interpret the tunneling spectra.The method can be applied to other materials such as transition metal dichalcogenides, borocarbides, or pnictides.

FIG. 2 .
FIG. 2. (Color online) (a) Dots: Experimental SIS tunneling spectrum measured at T = 5 K with an MgB 2 tip and a MgB 2 caxis-oriented film.Lines: Fit to a two-gap McMillan model.The parameters used for the fit are S = 2.6 meV; L = 2.6 meV; 0 S = 1.2 meV; 0 L = 6.5 meV; T S = 1; T L = 0. (b) Dots: Experimental SIS tunneling spectrum measured at T = 5.5 K with an MgB 2 tip and a V 3 Si c-axis-oriented surface.Lines: Fit to a two-gap McMillan model.The parameters used for the fit are S = 2.6 meV; L = 2.6 meV; 0 S = 1.3 meV; 0 L = 7.0; T S = 0.8; T L = 0.0.
FIG. 3. (Color online) Three different SIN tunneling spectra from Iavarone et al. [13] using a normal tip on a MgB 2 film with different cristal orientations and revealing the two-gap features.In (a) and (b) SIN spectra obtained along the c-axis orientation and the fit with the two-gap McMillan model.For the spectrum (a) there is a full contribution of the small gap band.The parameters used for the fit are S = 2 meV; L = 2 meV; 0 S = 0.7 meV; 0 L = 6.5 meV; T S = 1; T L = 0.For the spectrum (b) there is a small contribution of the large gap band.The parameters used for the fit are S = 2.9 meV; L = 2.41 meV; 0 S = 0.6 meV; 0 L = 7.5 meV; T S = 0.95; T L = 0.05.The hump around 7 meV is due to a small contribution of the large gap, probably due to a local defect.As a result, there is small change in the quasiparticle coupling compared to spectrum (a).(c) SIN tunneling spectrum obtained along the a-axis orientation.The parameters used for the fit are S = 2.5 meV; L = 2.08 meV; 0 S = 1.1 meV; 0 L = 7.5 meV; T S = 0.55; T L = 0.45.

FIG. 5 .
FIG. 5. (Color online) Calculated Fermi surface for MgB 2 .(a) Top view; (b) perspective view.The red and green 3D sheets come from the boron p z contribution, while the warped cylinders sheets in purple and blue come from the (p x ,p y ) bands.The solid black line shows the boundary of the Brillouin zone.

FIG. 6 .
FIG. 6. (Color online) Calculated constant-current STM images for the σ (a) and π (b) states, obtained for an iso-DOS value of 10 −4 e/eV/unit cell.The maxima in (a) correspond to the position on top of the center of the B-B bonds and in (b) to the position on top of the B atoms.

FIG. 7 .
FIG. 7. (Color online) Ratio of the calculated tunneling current for the π and σ states, as a function of the tip to surface distance, for different positions on the surface plane; z = 0 indicates the position of the first atomic layer.

FIG. 8 .
FIG. 8. (Color online) Shape φ(z) of the boron 2p x , 2p y , and 2p z orbitals along the c axis at different positions on the surface: (a) on top a (sub-surface) Mg atom; (b) at a quarter distance from two neighbor surface B atoms.z = 0 indicates the position of the first atomic layer.

TABLE I .
Parameters obtained from the fit of experimental spectra in the literature with the two-gap McMillan model.
a The film has crystallites oriented randomly.b Typical spectrum.c Unusual spectrum (large gap is sometimes observed).

TABLE II .
Bulk total and partial DOS at the Fermi level of MgB 2 for the σ and π bands [in e/eV/unit cell].