Electron-phonon coupling and its implication for the superconducting topological insulators.

The recent observation of superconductivity in doped topological insulators has sparked a flurry of interest due to the prospect of realizing the long-sought topological superconductors. Yet the understanding of underlying pairing mechanism in these systems is far from complete. Here we investigate this problem by providing robust first-principles calculations of the role of electron-phonon coupling for the superconducting pairing in the prime candidate CuxBi2Se3. Our results show that electron-phonon scattering process in this system is dominated by zone center and boundary optical modes, with coexistence of phonon stiffening and softening. While the calculated electron-phonon coupling constant λ suggests that Tc from electron-phonon coupling is 2 orders smaller than the ones reported on bulk inhomogeneous samples, suggesting that superconductivity may not come from pure electron-phonon coupling. We discuss the possible enhancement of superconducting transition temperature by local inhomogeneity introduced by doping.

The successful fabrication of topological insulators (TIs)1234 has stimulated the search for other exotic topological matters. Among them the topological superconductors (TSCs)2 are of particular interest owing to its topologically protected gapless surface states consisting of massless Majorana fermions, which may find potential application in future topological quantum computation. The close relationship between TIs and TSCs256 indicates that topological superconductivity may born from TI. Indeed, there have been several attempts to induce topological superconductivity in TI789101112. In particular, the recent observation of superconductivity in doped topological insulators13141516 has received great attention and has been considered as prime candidates for TSCs.

The superconductivity in these doped topological insulators is unique since it occurs at relative low carrier concentration, and particularly the Dirac surface states remain intact with the onset of superconductivity in the case of CuBi2Se314. Shortly after the experimental discovery, Fu et al.17 argued that CuBi2Se3 favors spin-triplet paring with odd-parity owing to its strong spin-orbit coupled band structure, and may realize a time-reversal-invariant (TRI) TSC. Although there is controversy on whether the Fermi surface of superconducting CuBi2Se3 encloses odd number of TRI momenta1819, which is key to a true TSC17, the prospect for a topologically nontrivial paring in this system is quite appealing. Indeed, subsequent point-contact spectroscopy measurements15 seem observe the zero bias conductivity peak which signifies the presence of Majorana surface states. But the result was challenged by other tunneling spectra measurements2021, leaving the topological aspect of superconductivity in this system has yet to be confirmed.

Despite much work concerning these systems to date22232425262728, the microscopic paring mechanism of the superconducting topological insulators remains elusive. On the one hand, unconventional superconductivity is generally associated with systems with strong electron correlations which seems unlikely for the typical topological insulators with sp electrons like Bi2Se3. On the other hand, whether electron-phonon (e-ph) coupling, which plays a central role in conventional BCS superconductivity, would realize the required unconventional paring is still an open question172425. To address this problem, a comprehensive understanding of phonons and e-ph coupling in doped topological insulators is thus strongly called for.

In this Report, we present a systematic study of phonons and role of e-ph coupling for the representative system CuBi2Se3 for a range of doping with accurate first-principles calculations. The e-ph properties are evaluated robustly by sampling the Brillouin zone (BZ) with several tens of thousands of inequivalent electron and phonon momenta. Our results reveal interesting renormalization of phonon dynamics due to electron doping, with coexistence of optical phonon stiffening and softening at zone center and boundary, signifying that the whole e-ph scattering process is dominated by these modes. Nevertheless, the obtained e-ph coupling constant λ is 0.28 for the optimally doped (x = 0.12) Bi2Se3, resulting in a maximum T of only ~0.03 K, which is significantly smaller than the ones reported on inhomogeneous bulk13141516 but in line with very recent experiment where disorder and inhomogeneity are strongly suppressed29. Our results thus indicate that the superconductivity in CuBi2Se3 may not be explained by only resorting to e-ph mechanism. We discuss the possibility of enhancement of superconducting transition temperature by local inhomogeneity introduced by copper doping.

Results

Doping induced strong renormalization of phonons dynamics

The Bi2Se3 crystallize in a rhombohedral structure with stacked Se-Bi-Se-Bi-Se quintuple layer (QL) as building blocks. Experimental data have consistently shown that copper atoms are mainly interpolated into the van der Waals gap between quintuple layers of Bi2Se3 and act as electron donor131430. And the bulk electronic structures are virtually intact upon copper doping apart from a rigid shift of whole energy bands3031 (see Fig. 1(b)). Hence we simulate doped topological insulator Bi2Se3 with the rigid-band approximation, where the excess of electron is compensated by uniform background of positive charge.

Figure 1

The crystal structure of CuBi2Se3 and energy band evolution with electron doping.

The crystal structure (a) and Brillouin zone (c) of CuBi2Se3. (b) the electronic energy bands for 3 representative Cu doping (x = 0.05e for an ellipsoidal Fermi surface, x = 0.12e for a nearly cylindrical Fermi surface and x = 0.2e corresponding to even more complex Fermi surface topology). To highlight the rigid shift of whole energy bands with electron doping, the Fermi level of each case has been set to 0.

The vibrational properties of nominally undoped bulk Bi2Se3 have been extensively studied using Raman spectroscopy3233 and all 4 Raman active optical phonon have been observed33. On the other hand, most available phonon data from fully-relativistic first-principles calculations are, however, unsatisfactory due to the presence of dynamical instability3435. Our phonon dispersion does not suffer from this problem and the calculated zone center optical modes, which can be classified according to the irreducible representation of the point group of bulk Bi2Se3 as Γ = 2(E + A1 + E + A2), agree well with experimental one33. The phonon spectrum extends up to 21.8 meV, and the vibrations of Se atoms mainly account for the top 9 optical branches and well separated in energy from Bi due to large mass difference.

When doped with copper atoms, the bulk conduction band of parent compound Bi2Se3 is partially filled due to charge transfer, and eventually develops superconductivity with copper concentration exceeds 0.10 when cooled down. Fig. 2 shows the phonon dispersions of bulk CuBi2Se3 corresponding to several doping levels (0 ≤ x ≤ 0.14) from a fully-relativistic calculation. As can be seen, while leaving the electronic structure almost unchanged30, electron doping strongly alters the vibrational properties of doped system, especially as reflected in the pronounced softening of the highest A2 zone boundary modes (Z point) and the coexistence of modes stiffening (E) and softening (A1 and A2) around zone center (Γ point). This doping-dependent phonon renormalization is most effect when doping starts (x = 0.05) and quickly saturates as the doping continues (especially at Z point as can be seen from Fig. 2), reflecting different phonon dynamics of bulk insulating Bi2Se3 and electron-doped one due to electron screening. We note that in most experiments, the bulk sample of Bi2Se3 is usually slightly electron-doped owing to the presence of Se vacancies. Hence the abrupt change of phonon frequencies (Z point) may not be easily observed. Indeed, previous Raman spectroscopy measurement36 showed no sign of significant change of phonon frequencies around zone center with Cu doping. Nevertheless, the special phonon stiffening and softening with electron doping here signifies these modes would dominate the e-ph scattering process. Actually, according to previous first-principles calculations3, the energy bands around Fermi level are dominated by the hybridized p orbitals from outermost Se atoms in a unit cell. And the A2 mode at Z and Γ involves out-of-phase motions of Bi and Se atoms along the z-direction. Hence with electron doping, the partially occupied anti-bonding states would strongly interact with this phonon pattern, resulting in relatively large phonon linewidths as will be shown below.

Figure 2

The electron doping dependence of phonon dispersion.

Phonon modes of interest are labeled according to the irreducible representation of point group of bulk Bi2Se3. The black circles highlight the major modifications of phonon modes (A2 at Z; A1, A2 and E at Γ) with copper doping.

Phonon linewidths from electron-phonon couping

The electron doping induced stiffening and softening of specific phonon modes is a clear sign of moderate e-ph interaction in this system. This is manifested in the relatively large phonon linewidths around zone center and zone boundary for the optimally doped (x = 0.12) CuBi2Se3 (see Fig. 3). In calculating phonon linewidths, we have sampled the BZ with 50 × 50 × 50 inequivalent electron wave vectors and the δ in Eq. (4) is replaced with 0.001 Ryd. From Fig. 3 we can see the significant phonon linewidths mainly lie in the top 9 optical branches which involve in-plane Se phonons and restricted to certain modes around the zone center (Γ) and zone boundary (Z) where phonon stiffening and softening occur. The phonon linewidths at other regions where phonon momentum q > 2k (k is the electronic Fermi momentum) are negligible. This localized distribution of phonon linewidths resembles that in MgB23738, where the E2 in-plane phonons near zone center strongly couples with partially occupied σ-bonding states and results in the largest e-ph coupling. This may suggest that CuBi2Se3 is a good e-ph superconductor, but this is unfortunately not the case as will be discussed later.

Figure 3

The phonon linewidths for optimally doped (x = 0.12e) CuBi2Se3 and doping dependence of electron-phonon properties.

Upper left: The radius of red circle is proportional to the magnitude of phonon linewidths. Upper right: Comparison of Eliashberg spectral function and phonon density of states for optimally doped case (x = 0.12e). Lower: The doping dependence of electron-phonon coupling constant λ (red point line) and electronic density of states at Fermi level N(E) (blue point line).

There is another question need to be answered, namely how phonon linewidths evolve with electron doping, since the topology of Fermi surface (FS) of CuBi2Se3 undergoes significant change from being an 3-D ellipsoid to 2-D cylindrical one1939. It turns out that the doping dependence of phonon linewidths is negligible although the scattering phase space is increasing monotonically.

Role of e-ph couping for superconductivity

Now let us turn to the e-ph coupling contribution to the superconductivity in CuBi2Se3. It has been shown the superconductivity of CuBi2Se3 occurs for a wide range of doping level (0.10 < x < 0.60)40, and the highest superconducting transition temperature of 3.8 K is achieved at x = 0.12. Within the framework of Migdal-Eliashberg theory4142, the superconductivity arising from e-ph coupling is characterized by several quantities4344, the most important ones are Eliashberg spectral functionand the dimensionless e-ph coupling constant . As mentioned above, a reliable λ needs dense sampling of both phonon and electron momentum. This is achieved here with the aid of Wannier interpolation, where the electronic and phonon states as well as e-ph coupling matrix elements are calculated on a fine mesh with 36 × 36 × 36 and 30 × 30 × 30 inequivalent electron and phonon momenta, respectively. This dense grid ensures the convergence of λ to within 0.005.

Fig. 3 summarizes our e-ph results for CuBi2Se3. From the comparison of phonon DOS and Eliashberg spectral function we find that the α2F(ω) generally follow the trend of phonon DOS for the low-lying Bi modes (ω < 10 meV) while deviate for the high-lying Se modes, and the discrepancy is significant around 13 meV and 20 meV. Such behavior has been seen in MgB2, where the most significant contribution to the remarkable high T of 39 K comes from phonon modes around 60 meV38. Despite this resemblance, the e-ph coupling constant λ for the optimally doped (x = 0.12) CuBi2Se3 is 0.28. The corresponding transition temperature T can be obtained from Allen-Dynes formula43:whereis the logarithmically averaged phonon frequency. From Fig. 4 we can conclude that for typical retarded Coulomb repulsion μ* of 0.10, the estimated T is only 0.03 K, which is 2 orders smaller than experimental value. If we assume a constant ω (~11 meV), to match the experimental value of T the required λ would be 0.6–0.7 (with reasonable Coloumb μ* of 0.06–0.12, see Fig. 4), which is ~2 times larger than present one. Since the electronic DOS at Fermi level N(E) is 4.46 states/Ryd/spin/(unit cell) which is almost identical to that in MgB2, it is clear that the major obstacle prevents CuBi2Se3 from being a good e-ph superconductor as MgB2 (with λ ~ 1) is the lower characteristic phonon frequency (~60 meV37 for MgB2 but ~11 meV here) or Debye temperature (~900 K45 for MgB2 but ~180 K46 here).

Figure 4

The dependence of T on λ and μ*.

T has been plotted against λ and μ*, from which we can estimate the lower bound of λ to be ~0.6 to match the experimental T of 4 K. The inset indicates the calculated T is far less than the experimental one for optimally doped case.

In Fig. 3 we also plot the electron doping dependence of λ and electronic DOS at Fermi level N(E). At low doping (x = 0.05) where FS is a closed ellipsoid, the carrier concentration (n) is 3.5 × 1020 cm−3, slightly larger than the experimental one (~2.0 × 1020 cm−3)13. The corresponding N(E) is about 2.64 states/Ryd/spin/(unit cell), and the calculated λ for this doping is even smaller (being 0.16). When increasing electron doping to x = 0.10, the N(E) increased to 3.86 states/Ryd/spin/(unit cell), and the λ jumps to 0.26. Further increasing of doping will generally elevate N(E) and λ, but the λ is still unable to exceed 0.30 at x = 0.14, although n at this doping level is already one order larger than the experimental one, which is probably unrealistic since the ambipolar doping nature of Cu tends to make electron concentration saturate at x = 0.130. Even higher doping would in addition result in an complex Fermi surface structure which is probably inconsistent with experiments1939.

Thus, given the experimental facts (with n ~ 2.0 × 1020 cm−3 and an ellipsoidal or cylindrical FS), the value of 0.28 corresponding to optimal doping has been the upper bound for λ and actually has overestimated it since the n is ~4 times larger than experimental one.

Discussion

Our results are actually in line with recent experiments on epitaxial CuBi2Se3 with thickness between 6 quintuple layers (QL) to 13 QL29, where superconducting transition is never observed down to 0.8 K even though the n is already comparable with bulk one (~1020 cm−3), suggesting that electron doping alone could not afford for the observed T of ~4 K and in particular, a weak e-ph contribution to the observed superconductivity. Moreover, previous angle-resolved photoemission spectroscopy measurements of e-ph coupling in Bi2Se3 also suggest a relatively small λ of 0.2547 and 0.1748, which are very close to our first-principles results.

Therefore, we have to resort to other mechanisms besides e-ph coupling to recover the observed T of ~4 K. Since electronic states of CuBi2Se3 around Fermi level are dominated by p character3, we would expect minor contribution from spin-fluctuations either. On the other side, we note that the Cu-intercalated structure is formed both in bulk and films of CuBi2Se32931, and the only difference is that disorder and inhomogeneity are strongly suppressed in latter case29. Indeed, there has been consistent observation of inhomogeneity of superconductivity in CuBi2Se313141516. And the superconducting shielding fraction reported also varies from group to group, in particular, Kriener et al.40 have shown that shielding fraction strongly depends on doping level x and T exhibits an unusual monotonic decrease with x, raising the possibility of phase segregation.

Hence we could infer that the local inhomogeneity introduced by copper may play nontrivial role in the superconductivity of CuBi2Se32940. Indeed, the intimate relationship between local inhomogeneity and superconductivity has been extensively studied in the context of high T superconductors, and the possible enhancement of superconductivity by local inhomogeneities has been discussed recently by Martin I. et al.49 in the weak coupling BCS regime. Given many unusual properties of current system (relatively high T compared to its low carrier density and the unexpected drop of shielding fraction for doping level x > 0.5)40, it is likely that the local inhomogeneity is indispensable to the superconductivity and may enhance T of CuBi2Se3.

Regarding the symmetry of the paring of superconducting CuBi2Se3, no consensus has been reached since its first observation. Experimentally, point-contact spectroscopy measurements seem support the spin-triplet paring with odd-parity15, but results from scanning tunneling microscope and Andreev reflection spectroscopy show no evidence of characteristic zero-energy surface bound states2021, suggesting the picture of s-wave paring. On the other hand, based on a 2 band model, Fu et al.17 argued that CuBi2Se3 favors spin-triplet paring owing to its strong spin-orbit coupled band structure, and the possibility of phonon-mediated odd-parity paring has also been discussed theoretically2425. Our first-principles results here suggest the paring of this system may be unconventional, and this unconventionality may come from the local inhomogeneity introduced by copper doping. This is because the superconducting transition hasn't observed in epitaxial CuBi2Se329, which would rule out the possibility that superconductivity is mediated by phonon modes derived from pure Bi2Se3.

In summary, we have presented a systematic study of phonons and role of e-ph coupling in the Cu doped Bi2Se3. Our results show that strongly renormalized zone center and zone boundary modes with electron doping would dominate the whole e-ph coupling process. Despite moderate e-ph coupling in this system, our robust first-principles calculations of e-ph properties for wide range of copper doping suggest that CuBi2Se3 is not a conventional e-ph superconductor and e-ph coupling plays minor role in superconductivity of this system. We have also discussed the possible enhancement of T by the local inhomogeneity introduced by copper doping. Above all, our results rule out a conventional phonon-mediated superconductivity in CuBi2Se3 and point to a delicate interplay between e-ph coupling of parent compound Bi2Se3 and nontrivial role played by inhomogeneity.

Methods

The e-ph properties have been obtained with the isotropic approximation to Migdal-Eliashberg theory4142. In this framework, the phonon self-energy (Π) arising from e-ph coupling is expressed as44where is the electronic energy with band index n and crystal momentum k, and ω is the vibrational frequency with branch index v and crystal momentum q. The f is the Fermi-Dirac distributions, and δ is a positive infinitesimal. The is the e-ph coupling vertex, where ΔV is the variation of the self-consistent potential induced by a collective ionic displacement. Note the spin degree has been incorporated into the band index.

The imaginary part of Π corresponds directly to phonon half-width at half-maximum γ, and the phonon mode-resolved e-ph coupling constant is given bywith N(E) being the electronic density of states (DOS) at the Fermi level.

An accurate determination of phonon linewidths and hence e-ph coupling constant requires fine energy and momentum resolutions of electronic and phonon states as well as e-ph coupling matrix elements. We achieved this by the recently developed interpolation method through Maximally Localized Wannier Functions (MLWFs)5051, which has been demonstrated to be extremely successful in addressing e-ph properties52535455 and other electronic properties where ultra high density of momentum is needed5657. As first step of this method, we employs standard density functional theory (DFT)58 and density functional perturbation theory (DFPT)59 to obtain converged ground state electronic density and dynamical matrix. In this study, we use experimental crystal structure with a = 4.138 Å and c = 28.64 Å. Fully relativistic norm-conserving pseudopo-tentials is used for all the calculation here as relativistic corrections are necessary for a satisfactory quantitative description of topological insulator Bi2Se3. A kinetic energy cutoff of 35 Ryd with methfessel-paxton smearing widths of 0.01 Ryd and Monkhorst-Pack grids of 12 × 12 × 12 for k point sampling are used to ensure the convergence of total energy. The dynamical matrix is obtained on a relatively coarse grid of 4 × 4 × 4 phonon wave vectors (convergence has been achieved by comparing with phonon dispersion from 6 × 6 × 6 mesh).

Subsequently, the electronic and vibrational quantities obtained in first step are interpolated to dense grids which contains several tens of thousands of inequivalent phonon (electron) wave vectors. In this interpolation step, we first construct 30 spinor Wannier functions (using p-like atomic orbitals of Bi and Se) to span a subset of full Hilbert space around Fermi level, on which the operators in momentum space like Hamiltonian and electron-phonon coupling matrix are projected. After obtaining operators in Wannier representation, we then interpolate back to the momentum space (Bloch representation) to obtain the converged results5051. By doing this way, the e-ph coupling constant λ is calculated on fine mesh grid containing 36 × 36 × 36 and 30 × 30 × 30 inequivalent electron and phonon momenta respectively. This dense grid ensures the convergence of λ to within 0.005.

Author Contributions

X.L.Z. performed coding and calculations. X.L.Z., W.M.L. analyzed numerical results and contributed in completing the paper.