# Beyond Eliashberg superconductivity in MgB: anharmonicity, two-phonon scattering, and multiple gaps

###### Abstract

Density-functional calculations of the phonon spectrum and electron-phonon coupling in MgB are presented. The E phonons, which involve in-plane B displacements, couple strongly to the electronic bands. The isotropic electron-phonon coupling constant is calculated to be about 0.8. Allowing for different order parameters in different bands, the superconducting in the clean limit is calculated to be significantly larger. The phonons are strongly anharmonic, and the non-linear contribution to the coupling between the modes and the p bands is significant.

###### pacs:

The recent discovery of superconductivity near 40 K in MgB has generated much interest in the properties of this simple intermetallic compound [1]. A significant B isotope effect strongly suggests phonon-mediated pairing [2]. To explain the large , an electron-phonon coupling (EPC) constant of is needed. Yet estimates of the coupling strength based on the latest measurements of the low-temperature specific heat [3], combined with the density of states (DOS) from density-functional calculations [4], yield 0.6-0.7. Further, the measured temperature dependence of the electrical resistivity [5] is consistent with . First-principles calculations of the EPC give [4, 6, 7]. Clearly there is a problem in reconciling all these numbers. Another puzzle involves tunneling measurements of the gap. Values of ranging from 1.2 to 4 have been reported. The values below the BCS weak coupling limit of 3.5 have been attributed to surface effects, but the best-quality spectra[8] show a very clean gap with Sharvin contact measurements[9] reveal a gap at 4.3 meV (, and additional structures at and 3, raising the possibility of multiple gaps. Careful analysis of the temperature and magnetic field dependence of the specific heat suggests anisotropic or multiple gap structure as well [3]. Thus, even if superconductivity in MgB is phonon-mediated, it is likely that an analysis beyond the simple isotropic Eliashberg model is needed.

The MgB lattice consists of two parallel systems of flat layers. One layer contains B atoms in a honeycomb lattice, the other Mg atoms in a triangular lattice halfway between the B layers. First-principles calculations[4] find that the electronic states near the Fermi level are primarily B in character and the Fermi surface (FS) comprises four sheets: two nearly cylindrical hole sheets about the -A line arising from quasi-2D B bands, and two tubular networks arising from 3D bonding and antibonding bands [4]. The difference in character between the sheets raises the possibility that each has a distinct gap that could be observed in the clean limit. Such interband anisotropy enhances the effective EPC constant relevant to superconductivity and decreases the coupling constant for transport, compared to the average values [10, 11, 12]. This could explain the discrepant values of deduced from different types of experiments.

In this paper, we report first-principles calculations of the EPC in MgB. The coupling constant is decomposed into contributions from the four different bands crossing the Fermi level, allowing for an analysis of the effects of interband anisotropy on the superconducting and gap structure. The strongest coupling arises from the phonon modes along the to A line, which strongly interact with the quasi-2D electronic states. This phonon mode is calculated to be highly anharmonic, and it also has significant nonlinear contributions to the EPC[13].

Harmonic phonon frequencies and linear EPC parameters were calculated using the linear-response method within the local density approximation [14]. Norm-conserving pseudopotentials [15] were used, with a plane-wave cutoff of 50 Ry. We used the experimental crystal structure, with Å and [16]. The electronic states were sampled on grids of up to k-points in the full Brillouin zone, and the dynamical matrix was calculated on a grid of phonon wavevectors q [17].

The calculated phonon density of states and Eliashberg function are plotted in Fig. 1. The results are similar to those reported in Refs. [6] and [7]. All of these calculations give a slightly softer phonon spectrum than what is observed in neutron experiments [18]. While and are similar in shape in many materials, they are strikingly different in MgB. In particular, has a pronounced peak in the range of 60 to 70 meV arising from dispersive optic modes that do not give rise to large structures in . Correspondingly, the average phonon frequency meV is less than the logarithmically averaged frequency meV, despite the fact that logarithmic averaging preferentially weights lower frequencies. The isotropic EPC constant, which determines in the dirty limit, is found to be 0.77, in reasonable agreement with other calculations [4, 6, 7].

The peak in between 60 and 70 meV arises from the phonon modes with q along the -A line. This Raman active phonon mode, doubly degenerate at , involves in-plane, hexagon-distorting displacements of the B atoms. In fact, by symmetry, this is the only mode at that has a linear EPC. Going away from the -A line the EPC drops sharply when the phonon wave vector q becomes larger than the diameter of the 2D Fermi surface; at the same time the frequency increases by roughly 30%. This indicates that the reason why this B-B bond-stretching mode is not the highest-frequency mode at is because of softening due to EPC. However, this softening should weaken in the superconducting state, since some of the screening electrons form Cooper pairs and are removed from the Fermi sea. Such hardening of optical phonons below was first discussed by Zeyher and Zwicknagl [19]. The overall scale of the relative hardening, is set by a specific EPC constant, where is the EPC matrix element. (The Fermi level is set to zero.) In the BCS limit, is a known analytical function[20] of . We calculate , for the mode. Taking meV we predict about a 12% hardening of this mode below . This shift should be easily observable in Raman or neutron experiments.

Since the 2D FSs are calculated to play an important role in the EPC, we have decomposed the relevant electronic characteristics in terms of the four sheets of the FS. We list in Table I the partial DOS , and plasma frequencies referring to the light(heavy)-hole 2D sheets of the Fermi surface, and to the bonding (antibonding) sheets. The EPC constant was also decomposed into contributions from scattering of an electron from band to band : , with

Here is the frequency of the corresponding phonon, and is the standard (Eliashberg) isotropic coupling constant. Allowing for interband anisotropy of the order parameter (clean limit), the effective coupling constant for superconductivity is given by the maximum eigenvalue of the matrix , which is always larger than . Assuming the same interaction parameters for transport properties, the lowest order variational approximation for the Boltzmann equation corresponds to the transport EPC constant . On the other hand, allowing variational freedom for the different sheets of the Fermi surface yields an effective transport coupling constant which is always smaller than . In effect, the different bands provide parallel channels for conduction, so that when “scattering-in” is neglected, [12].

The calculated interaction parameters are listed in Table II. Due to similarities between the two 2D sheets, and between the two 3D sheets, we have simplified the model to allow for two different order parameters for these two sets of bands. This gives Ry, Ry, and Ry, where and stand for the 2D and 3D bands, respectively. Then and suggesting de Haas-van Alphen mass renormalizations of and for the two sets of bands, and specific heat renormalization of 1.77[21]. The resulting anisotropic effective coupling constant for superconductivity is . Using the Allen-Dynes approximate formula for [22], we find that to have K, a Coulomb pseudopotential of is needed. This is more physically plausible than the required when is used. For transport, interband anisotropy reduces the in-plane coupling constant from 0.70 to 0.57, but has essentially no effect on the out-of-plane (Table III). This is because the anisotropic formula accounts for the fact that the transport is mostly due to the 3D bands in any direction, simply because they couple less with phonons. The measured resistivity [5] can be fit remarkably well with the Bloch-Grüneisen formula using the calculated isotropic , with the in-plane and out-of-plane contributions appropriately averaged for polycrystalline samples [23]. However, with band anisotropy, the resistivity is slightly underestimated.

The temperature dependence of the individual gaps in the weak-coupling multigap model is defined by As shown in Fig. 2, the larger 2D gap is calculated to be BCS-like, with a slightly enhanced , while the 3D gap is about three times smaller in magnitude. Thus, in the clean limit, MgB should have two very different order parameters, which in turn should affect thermodynamic properties in the superconducting state. Experiments indicate that the coherence length in MgB is close to 50 Å. The mean free path corresponding to the residual resistivity observed in Ref. [5] is more than 1000 Å, so that This is in the reasonably clean regime, and it is likely that the intrinsic resistivity is even smaller. However, stronger defect scattering should be detrimental to superconductivity: using the Allen-Dynes formula with the same we get an isotropic K. Indeed, irradiation has been found to drastically reduce [24]. Some of the experimental manifestations of multigap superconductivity would be a reduced and impurity-sensitive specific-heat jump at , a deviation of the critical-field temperature dependence from the Hohenberg-Werthamer formula, a reduction of the Hebel-Slichter peak in NMR, and a substantial difference between the in-plane and out-of-plane tunneling spectra. In particular, the latter should see only the smaller gap[25].

Note that is in the intermediate coupling regime. Furthermore, the multigap scenario suggests particular sensitivity to impurity scattering. This means one should really solve the anisotropic Eliashberg equations with impurity scattering, rather than the weak coupling BCS equations we used. Thus we do not make any quantitative thermodynamic and spectroscopic predictions here.

We focus now on the phonon modes, which carry the lion share of the EPC. We have examined this mode at in detail using frozen-phonon calculations. The calculations were carried out using a general potential LAPW code with the same setup as in Ref. [4]. This mode has substantial anharmonicity. A fit of the total energy for B displacements between a.u. to a fourth-order polynomial () gives Ry/a.u., Ry/a.u. and Ry/a.u.. Anharmonicity increases the frequency of this mode by about 15%, which should result in an overall reduction of by 10%, and an increase of by 6%.

More interestingly, the modes have a significant nonlinear coupling with electrons. The linear coupling vertex, , corresponding to scattering by a single phonon, is proportional to matrix elements of , where while the second-order coupling, involving exchange of two phonons, is proportional to matrix elements of . At , the Hellman-Feynman theorem allows the calculation of via deformation potentials. This is no longer the case for . One can use only as a qualitative estimate of . For the cylindrical sheets of the Fermi surface, and 20 mRy as compared to and 16 mRy. This suggests that nonlinear pairing via two-phonon exchange is comparable to or even larger than the linear coupling[26]. The reason for this anomalous behavior lies in the specifics of the band structure of the 2D bands. In the nearest-neighbor tight-binding approximation it can be described as

where and are the three smallest lattice vectors. At the point and there are two doubly degenerate states. The bonding pair forms the 2D Fermi surfaces, and it appears that the main effect of the phonons is to lift the degeneracy at , thereby changing the overall splitting between the two subbands. This effect does not depend on the sign of the ionic displacement (the degeneracy is lifted either way) and thus is nonlinear by definition. This is illustrated in Fig. 3 by the two Fermi surface plots for two opposite phonon patterns, both corresponding to . One can see that while for the 3D sheets the coupling is mostly linear (changes of the Fermi surface are opposite), for the 2D cylinders it is mostly quadratic (changes are the same, cf. the undistorted Fermi surface in Ref.[4]). Nonlinear EPC is also a likely source of anharmonicity: as discussed above, the contribution to the phonon self-energy from the EPC with the 2D FSs amounts to as much as 20 meV, as evidenced by the softening of this phonon at . A sizeable part of this softening probably comes from two-phonon processes. Quartic anharmonicity of the ion-ion interaction arises in the fourth order in the linear interaction constant, but in the second order in the nonlinear one.

In summary, we have presented a first-principles investigation of the electron-phonon coupling in MgB. Interband anisotropy enhances the coupling constant from its isotropic dirty-limit value of to an effective clean-limit value of for superconductivity. With meV, this is arguably consistent with the measured of nearly 40 K. In the clean limit, we predict two different superconducting order parameters: a larger one on the 2D FSs and a smaller one (by approximately one third) on the 3D FSs. Since current experiments suggest that MgB is indeed in the clean limit, multiple gaps should be observable. There are hints of this in both the tunneling and thermodynamic data.

The phonon mode involving in-plane B motion provides the strongest coupling. We predict a hardening of % of this mode at the zone center below . In addition, this mode is highly anharmonic and it may also have significant nonlinear electron-phonon coupling. The former likely reduces the linear EPC, while the latter increases the total EPC. Further work is needed to better elucidate the effect of the anharmonicity and nonlinear coupling on the superconducting properties of this material.

We thank Jim Freericks for helpful discussions. This work was partially supported by the National Science Foundation under Grant DMR-9973225, and by ONR. AYL acknowledges support from the ASEE-US Navy Faculty Sabbatical Program.

total | 1 | 2 | 3 | 4 | |
---|---|---|---|---|---|

4.83 | 0.66 | 1.38 | 1.26 | 1.52 | |

7.21 | 2.91 | 2.95 | 3.05 | 5.04 | |

6.87 | 0.44 | 0.52 | 4.62 | 5.06 |

11 | 12 | 13 | 14 | 22 | |
---|---|---|---|---|---|

(Ry) | 0.676 | 0.419 | 0.064 | 0.096 | 0.477 |

23 | 24 | 33 | 34 | 44 | |

(Ry) | 0.064 | 0.097 | 0.113 | 0.106 | 0.092 |

isotropic | 0.77 | 0.70 | 0.46 | 0.58 |
---|---|---|---|---|

multi-gap | 1.01 | 0.58 | 0.46 | 0.53 |

## References

- [1] J. Nagamatsu et al., Nature (London) 410, 63 (2001).
- [2] S. L. Bud’ko et al., Phys. Rev. Lett. 86, 1877 (2001).
- [3] Y. Wang, T. Plackowski and A. Junod cond-mat/0103181; N.E. Phillips, to be published.
- [4] J. Kortus et al., Phys. Rev. Lett. xx, xxx (2001).
- [5] P. Canfield et al., cond-mat/0102289.
- [6] Y. Kong et al., cond-mat/0102499 (2001).
- [7] K.-P. Bohnen, R. Heid, B. Renker, cond-mat/0103319
- [8] G. Rubio-Bollinger, H. Suderow and S. Vieira, cond-mat/0102242.
- [9] H. Schmidt et al., cond-mat/0102389.
- [10] H. Suhl, B. T. Matthias, and L. R. Walker, Phys. Rev. Lett. 3, 552 (1959).
- [11] W. H. Butler and P. B. Allen, in Superconductivity in d- and f-band Metals, edited by D. H. Douglass (Plenum, New York, 1976).
- [12] I.I. Mazin et al., Physica C 209, 125 (1993).
- [13] T. Yildirim et al, cond-mat/0103469.
- [14] A.A. Quong and B.M. Klein, Phys. Rev. B 46, 10734 (1992); A.Y. Liu and A.A. Quong, ibid. 53, 7575 (1996).
- [15] N. Troullier and J.L. Martins, Phys. Rev. B 43, 8861 (1991).
- [16] Muetterties, The Chemistry of Boron and Its Compounds (Wiley, New York, 1967).
- [17] Neglecting interband scattering, it is not possible for an optic mode at q=0 to couple to electrons since energy cannot be conserved. The contribution to the EPC from the region near is estimated by assuming a constant matrix element in this region and approximating the Fermi surface as two cylindrical sheets [20].
- [18] R. Osborn et al., cond-mat/0103064.
- [19] R. Zeyher and G. Zwicknagl, Solid State Commun. 66, 617 (1988); Z. Phys. B 78, 175 (1990).
- [20] C. O. Rodriguez et al., Phys. Rev. B 42, 2692 (1990).
- [21] Current experiments[3] indicate renormalization of 1.6, but the LDA may overestimate the DOS by 10-15%. The overestimation of NMR relaxation rates by the LDA provides indirect evidence of this (E. Pavarini, unpublished).
- [22] P.B. Allen and R.C. Dynes, Phys. Rev. B 12, 905 (1975).
- [23] D. Stroud, Phys. Rev. B 12, 3368 (1975).
- [24] A.E. Karkin et al., cond-mat/0103344
- [25] In Ref. [8], where high-quality spectra with a gap three times smaller than the BCS value were obtained, a unique procedure was used for sample preparation, which may have preferentially oriented individual microcrystals.
- [26] For linear coupling, the total isotope effect of all components must add up to 0.5; nonlinear coupling may change this total in either direction.