Adiabatic and non-adiabatic phonon dispersion in a Wannier function approach
aa r X i v : . [ c ond - m a t . m t r l - s c i ] J u l Adiabatic and non-adiabatic phonon dispersion in a Wannier function approach
Matteo Calandra , Gianni Profeta , and Francesco Mauri CNRS and Institut de Min´eralogie et de Physique des Milieux condens´es,case 115, 4 place Jussieu, 75252, Paris cedex 05, France and Consiglio Nazionale delle Ricerche - Superconducting and Innovative Devices (CNR-SPIN), 67100 L’Aquila, Italy (Dated: November 3, 2018)We develop a first-principles scheme to calculate adiabatic and non-adiabatic phonon frequenciesin the full Brillouin zone. The method relies on the variational properties of a force-constantsfunctional with respect to the first-order perturbation of the electronic charge density and on thelocalization of the deformation potential in the Wannier function basis. This allows for calculationof phonon dispersion curves free from convergence issues related to Brillouin zone sampling. Inaddition our approach justify the use of the static screened potential in the calculation of thephonon linewidth due to decay in electron-hole pairs. We apply the method to the calculation ofthe phonon dispersion and electron-phonon coupling in MgB and CaC . In both compounds wedemonstrate the occurrence of several Kohn anomalies, absent in previous calculations, that aremanifest only after careful electron and phonon momentum integration. In MgB , the presenceof Kohn anomalies on the E g branches improves the agreement with measured phonon spectraand affects the position of the main peak in the Eliashberg function. In CaC we show that thenon-adiabatic effects on in-plane carbon vibrations are not localized at zone center but are sizablethroughout the full Brillouin zone. Our method opens new perspectives in large-scale first-principlescalculations of dynamical properties and electron-phonon interaction. PACS numbers: 63.20.dk, 71.15.-m, 63.20.kd
I. INTRODUCTION
Electron-phonon (EP) interaction is responsible ofmany important phenomena in solids. As an example,the temperature behavior of the electron relaxation timein metals is to a great extent due to the scattering be-tween carriers and atomic vibrations such that finitetemperature transport is largely ruled by the EP inter-action. Similarly, the high temperature heat capacity inmetals is enhanced by the increased electronic mass dueto the interaction with lattice vibrations . In metals, atlow temperatures, EP coupling can generate a supercon-ducting state in which electrons move with no electricalresistance . It also can increase the effective mass of thecarriers so much that the system is driven from a metallicto an insulating state, as it happens in case of polaronic orPeierls instabilities . Finally, the electron-phonon scat-tering is often the largest source of phonon damping inphonon-mediated superconductors .First-principles theoretical determination of theelectron-phonon coupling strength in solids requiresthe calculation of the electronic structure, the vibra-tional properties, and the electron-phonon couplingmatrix elements. In state-of-the-art electronic struc-ture calculations these quantities are obtained usingthe adiabatic Born-Oppenheimer approximation anddensity-functional theory (DFT) in the linear-responseapproach .More specifically, within the Born-Oppenheimer ap-proximation, the determination of phonon frequencies,that are related to the real part of the phonon self-energy, requires the calculation, in a self-consistent man-ner , of the variation of the Kohn-Sham potential V SCF with respect to a static ionic displacement u , namely δV SCF ( r ) δ u . As the displacement of the ions is static,the obtained δV SCF ( r ) δ u is real. Conversely, the phononlinewidth, related to the imaginary part of the phononself-energy, is obtained in a non self-consistent proce-dure using the electron charge density and the previouslydetermined δV SCF ( r ) δ u . The advantage of using a non self-consistent procedure to study phonon linewidth is mainlyrelated to the less expensive computational load with re-spect to a self-consistent one. In addition, as recentlydemonstrated, interpolation schemes of the electron-phonon matrix elements can be used to calculate theimaginary part of the phonon self-energy on ultra-densek-point grid.It is however unclear to what extent this procedure ofcalculating self-consistently the real part of the phononself-energy and non self-consistently the imaginary partis actually correct. Indeed, the proper way of treatingphonons should be to consider a monochromatic time-dependent displacement (at the phonon frequency ω ) andperform a time-dependent self-consistent linear-responsescheme. The resulting variation of the self-consistentpotential δV SCF ( r ,ω ) δ u would then be a complex quantity.The real part of the resulting phonon self-energy wouldthen determine the phonon frequencies while its imagi-nary part would lead to the phonon linewidth. In thisway, non-adiabatic (in the sense of Ref. 12) dynamical phonon frequencies could be accessed. Although feasiblein principle this procedure would require a full rewritingof the linear response code including time dependenceand it would also be more expensive then a standardstatic linear-response calculation. Moreover, in the pres-ence of Kohn anomalies and long-range force-constants,where extremely accurate k-point sampling of the Fermisurface is needed to converge the phonon self-energy, thecalculation would be unfeasible.Thus, it would be desirable to have a non -self consis-tent linear-response formulation to obtain both the realand imaginary phonon self-energy, both within the adia-batic/static and non-adiabatic/dynamic approximation.In this work we develop a scheme to calculate non self-consistently both the real and the imaginary part of thephonon self-energy by using a functional that is varia-tional with respect to the variation of the self-consistentcharge density. Our method opens the way to calculateadiabatic/static and non-adiabatic/dynamic phonon fre-quencies using an ultra-dense sampling of the electronand phonon wave-vector in a non self-consistent way,starting from the static self-consistent variation of thetime independent Kohn-Sham potential ( δV SCF ( r ) δ u ), ob-tained with a coarse sampling. In this way, the main com-putational load related to phonon frequencies calculationis drastically reduced. The efficiency of the method is fur-ther enhanced by using the interpolation of the electron-phonon matrix elements based on Wannier functions .The method is applied to study the dynamical andsuperconducting properties of MgB and CaC , two ofthe most studied superconducting materials in the lastfew years. These applications are meaningful and com-putationally challenging. In fact, although many exper-imental and theoretical studies appeared, there are stillimportant open issues and debated results.For example, conflicting results for the calculate MgB electron-phonon coupling ( λ ) are present in the litera-ture mainly due to difficulties in Brillouin Zone sampling.This, apart its fundamental importance, prevents a fullunderstanding of the normal and superconducting prop-erties of this material.In the case of CaC , very large non-adiabatic effectswere predicted and measured at zone-center. It isunclear to what extend non-adiabatic effects are sizablefar from Γ-point. In this last case, the adiabatic effectscan be relevant for thermodynamic properties resultingfrom an average of the phonon frequencies over the Bril-louin zone.The paper is organized as follows. In Section II wederive the variational formulation of the force-constantmatrix. This will be done starting from the formal def-inition of force constants in the case of a monochro-matic time-dependent ionic displacement (Sec.II A). InSec. II B we discuss when non-adiabatic/dynamic effectsshould be expected in the phonon spectra. The linearresponse equations for the dynamical matrix are intro-duced (Sec.II C) in the general case and then specializedin the density functional formulation (Sec.II D).In Section II E, we develop the variational formulationof the force constants in the density functional linear-response theory and outline the computational frame-work of phonon and electron-phonon coupling (Sec.II Fand II G). In Sec. II H we show how our approach justifythe use of the static screened potential in the calculationof the phonon linewidth due to decay in electron-hole pairs. Section III will be devoted to the description ofthe implementation of the theory using the Wannier in-terpolation scheme.In Section IV we reports results on dynamical and su-perconducting properties of MgB (Sec.IV A) and CaC (Sec.IV B).Section V summarizes our conclusions. II. THEORYA. Time dependent dynamical matrix and phonondamping.
We consider a crystal with N atoms in the unit cellsubmitted to a time dependent perturbation, in whichthe position of an atom is identified by the vector R I ≡ R L + τττ s + u I ( t ) , (1)where R L is the position of the L -th unit cell in theBravais lattice, τττ s is the equilibrium position of the s -thatom in the unit cell, u I ( t ) indicates the deviation fromequilibrium of the nuclear position and I = { L, s } . Theforce at time t acting on the J -th nucleus ( J = { M, r } )due to the displacement u I ( t ′ ) of the atom I -th at time t ′ is labeled F J ( t ). The force constants matrix is definedas: C IJ ( R L − R M ; t − t ′ ) = − δ F J ( t ) δ u I ( t ′ ) (2)where we used the translational invariance of the crystaland make evident the dependence of C IJ on the latticevector R L − R M (to lighten the notation we omit it inthe following equations where no confusion may arise).The ω -transform of the force-constants matrix is thus: C IJ ( ω ) = Z dte iωt C IJ ( t ) (3)While the force-constants matrix C IJ ( t ) is a real quan-tity, its ω -transform C IJ ( ω ) is not real and has both areal and imaginary part. The Fourier transform of theforce-constant matrix is C sr ( q , ω ) = X L e − i qR L C Ls,Mr ( ω ) (4)where, without loss of generality, we have chosen R M = . The Hermitian and anti-Hermitian combination of theforce-constant matrix in momentum space are: D sr ( q , ω ) = 12 √ M s M r [ C sr ( q , ω ) + C rs ( q , ω ) ∗ ] (5) A sr ( q , ω ) = 12 i √ M s M r [ C sr ( q , ω ) − C rs ( q , ω ) ∗ ] (6)where M s is the mass of the s-th atom in the unit cell.These quantities are associated to the real and imaginarypart of the dynamical matrix in coordinate space. If theimaginary part of the dynamical matrix is small withrespect to its real part, namely | A sr ( q , ω ) | << | D sr ( q , ω ) | (7)then the self-consistent conditiondet (cid:12)(cid:12) D sr ( q , ω q ν ) − ω q ν (cid:12)(cid:12) = 0 (8)determines non-adiabatic/dynamic phonon frequencies ω q ν and phonon eigenvectors (cid:8) e s q ν (cid:9) s =1 ,N and ν =1 , N indicates the phonon branches. The adia-batic/static phonon frequencies and eigenvectors are ob-tained considering a static perturbation, thus diagonal-izing D rs ( q , ω q ν = 0).On the other hand, the imaginary part of the force-constants matrix determines the phonon-damping γ q ν = 2 ω q ν X s,r e s q ν A sr ( q , ω q ν ) e r q ν (9)and the phonon linewidth. B. Adiabatic and non-adiabatic phonons
In the previous section we have defined the non-adiabatic/dynamic phonon frequencies as the eigenval-ues of the Hermitian part of the time dependent dynam-ical matrix, and the adiabatic/static phonon frequen-cies as the eigenvalues of the static time-independentdynamical matrix. Solid-state text-books and first-principles calculations of the phonon dispersion ,usually treat only the adiabatic/static case, since it iscommonly assumed that the adiabatic phonon frequen-cies coincide with the non-adiabatic ones.In insulators, where the fundamental energy gap be-tween the electronic ground state and the first availableexcited state is much larger than the phonon energy, theadiabatic/static approximation is well justified. In met-als the situation is more complex.
The crucial parameter in metals is the electron relax-ation time τ . In absence of electron-defect, electron-electron and electron-phonon scattering, τ is infinite. Ina real metallic system, the presence of these scatteringprocesses results in a finite relaxation time τ . We candefine three cases: (i) a clean limit, when the electronrelaxation time τ is much larger than the phonon pe-riod divided by 2 π (ii) a dirty limit, when the electronrelaxation time τ is much smaller than the phonon pe-riod divided by 2 π (iii) an intermediate regime, when theelectron relaxation time τ is comparable to the phononperiod divided by 2 π .It has been shown that in the dirty limit thenon-adiabatic phonon frequencies coincide with the adi-abatic ones also in metals. Instead, in the clean limitand in the intermediate regime, the adiabatic and non-adiabatic frequencies are, in general, different. The non-adiabatic calculations based on time-dependent DFT (described in the following sections) areperformed in the perfect clean limit. Indeed, in our time-dependent calculations, the electron relaxation time τ isinfinite, since we use an instantaneous (real in the fre-quency space) exchange-correlation Kernel, and we donot consider a broadening of the electronic levels due todefects and electron-phonon scattering.Thus, our non-adiabatic/dynamic DFT frequenciesshould be used to reproduce the phonon frequencies mea-sured in metal in the clean-limit. Instead, the adia-batic/static DFT frequencies (those generally computedwith DFT linear response codes) should be used to repro-duce the phonon frequencies measured in a metal in thedirty limit. Finally, to reproduce the phonon frequen-cies measured in a metal in the intermediate regime, oneshould explicitly include electron-lifetime effects in thelinear response calculation, as it has been done for zonecenter phonon in Refs. 12,22.By considering a single band and by linearizing theelectronic-band dispersion near the Fermi energy, it hasbeen shown that the differences between the adiabaticand non-adiabatic phonons are largest at the center ofthe Brillouin zone (BZ) and vanish for q ≫ ω q ν /v F ,where v F is the Fermi velocity. Such differences, at theBZ center and in the clean limit, have been computedwithin DFT , and can be very sizable (up to 30% of thephonon frequency). However a detailed DFT study ofnon-adiabatic effect away from the BZ center (beyond alinearized one band approximation ) is still missing. C. Time dependent linear response theory
The force-constant matrix can be evaluated in the lin-ear response theory, considering that the atomic displace-ment induces a perturbation in the external potential act-ing on the electrons.The Hellmann-Feynmann theorem states that theforce on atom J , F J ( t ) can be evaluated in terms of thevariation of the external potential: F J ( t ) = − Z d r n ( r , t ) δV ext ( r ) δ R J (10)where n ( r , t ) is the electronic charge density and V ext ( r )is the external potential, namely: V ext ( r ) = − X I Z | r − R I | + 12 X I = J Z I Z J | R J − R I | (11)It is worthwhile to recall that V ext ( r ) does not dependexplicitly on time but only through the dependence ontime of the phonon displacement u I ( t ). In linear responsethe force-constant matrix C IJ ( ω ) is written as : C IJ ( ω ) = Z δn ( r , ω ) δ u I δV ext ( r ) δ u J d r + Z n ( r ) δ V ext ( r ) δ u J δ u I d r . (12)where n ( r ) is the unperturbed charge density, n ( r , ω )is the ω -transform of the time-dependent charge den-sity n ( r , t ) and δV ext ( r ) δ u J and δ V ext ( r ) δ u I δ u J are real and time-independent quantities evaluated at the equilibrium po-sition of the nuclei. On the contrary δn ( r ,ω ) δ u I is complex. D. Time-dependent linear response in densityfunctional theory
In this section, we examine the linear-response calcu-lation of the real-space force-constants in the frameworkof density functional theory. The derivations are keptin real space as we want to keep the presence of imag-inary terms (damping) in the force-constants matrix asmanifest.In order to compute n ( r , ω ), we assume a monochro-matic perturbation of the form u I ( t ) = u I ( ω )( e iωt + e − iωt ) (13)where u I ( ω ) is real. In density functional theory theresulting monochromatic perturbing potential is the ex-ternal potential V ext ( r ) in Eq. 11 . The derivative of the ω -transform of the charge density with respect to a ionicdisplacement is written as n I ( r , ω, T ) = δn ( r , ω ) δ u I = 2 N k ( T ) X k i, k ′ j ( f k i ( T ) − f k ′ j ( T )) × h ψ k ′ j | δV SCF ( r , ω ) /δ u I | ψ k i i ǫ k i − ǫ k ′ j + ω + iη ψ ∗ k i ( r ) ψ k ′ j ( r ) (14)where k , k ′ label the crystal momentum, i , j are band in-dexes, the T is the electronic temperature (or broaden-ing) in the Fermi function f k i ( T ), η is an arbitrarily smallpositive real number, N k ( T ) is the number of k-pointsneeded to converge the sum in Eq. 14 at the electronictemperature T and the factor 2 (here and in the follow-ing) accounts for the spin-degeneracy. Finally, V SCF isthe Kohn-Sham self-consistent potential. The quantity δV SCF ( r , ω ) /δ u I is complex and given by δV SCF ( r , ω ) δ u I = δV ext ( r ) δ u I + Z K ( r , r ′ ) n I ( r , ω, T ) d r ′ . (15)where K ( r , r ′ ) = δE HXC [ n ] δn ( r ) δn ( r ′ ) is the kernel of the Hartreeand Exchange and correlation functional, E HXC [ n ]. Asusual we have assumed the Hartree and exchange-correlation Kernel to be instantaneous (real in ω space).Substitution of Eq. 14 in Eq. 12 leads to: C IJ ( ω, T ) = 2 N k ( T ) X k i, k ′ j ( f k i ( T ) − f k ′ j ( T )) × h ψ k ′ j | δV SCF ( r , ω ) /δ u I | ψ k i ih ψ k i | δV ext ( r ) /δ u J | ψ k ′ j i ǫ k i − ǫ k ′ j + ω + iη + Z n ( r ) δ V ext ( r ) δ u J δ u I d r . (16) where, from now on, we explicitly indicate the depen-dence on the electronic temperature T . This expressionof the force-constants matrix is normally used in stan-dard implementations of linear-response theory .Further substitution of Eq. 15 in Eq. 16 gives thefollowing alternative, but equivalent formulation for theforce-constants matrix in linear response theory: C IJ ( ω, T ) = 2 N k ( T ) X k i, k ′ j f k i ( T ) − f k ′ j ( T ) ǫ k i − ǫ k ′ j + ω + iη × h ψ k ′ j | δV SCF ( r , ω ) δ u I | ψ k i ih ψ k i | δV SCF ( r , ω ) δ u J | ψ k ′ j i + Z d r n ( r ) δ V ext ( r ) δ u I δ u J − Z Z n J ( r , ω, T ) K ( r , r ′ ) n I ( r ′ , ω, T ) d r d r ′ . (17)In Eq. 17 the term including the second derivative ofthe external potential is real whereas all the other termsare complex. The advantage of Eq. 17 is that it allowsus to introduce a variational formulation of the force-constants matrix as it will be shown in the following. E. Variational formulation of the force-constants.
We introduce the following force-constants functional F IJ : F IJ [ ρ ( r ) , ρ ′ ( r ) , ω, T ] == 2 N k ( T ) X k i, k ′ j f k i ( T ) − f k ′ j ( T ) ǫ k i − ǫ k ′ j − ω + iη × h ψ k ′ j | δV ext ( r ) δ u I + Z K ( r , r ′ ) ρ ( r ′ ) d r ′ | ψ k i i× h ψ k i | δV ext ( r ) δ u J + Z K ( r , r ′ ) ρ ′ ( r ′ ) d r ′ | ψ k ′ j i + Z d r n ( r ) δ V ext ( r ) δ u I δ u J − Z Z ρ ( r ) K ( r , r ′ ) ρ ′ ( r ′ ) d r d r ′ . (18)With this definition the force-constants matrix reads as: C IJ ( ω, T ) = F IJ [ n I ( r , ω, T ) , n J ( r , ω, T ) , ω, T ] (19)It is straightforward to note that the force-constantfunctional is quadratic with respect to ρ , namely: δF IJ [ ρ ( r ) , ρ ′ ( r ) , ω, T ] δρ ( r ) (cid:12)(cid:12)(cid:12)(cid:12) ρ ( r )= n I ( r ,ω,T ) ,ρ ′ ( r )= n J ( r ,ω,T ) = 0 . (20)The same relation holds upon derivation with respect to ρ ′ .The consequence of Eq. 20 is that C IJ ( ω, T ) is anextremal point and that a given error on n I ( r , ω, T ) oron n J ( r , ω, T ) affects the functional and the phonon fre-quencies only at second order.Our functional formulation is related to that used fordielectric tensors by Gonze et al. in Ref. 26. The maindifference is that, while the functional formulation of Ref.26 is variational with respect to the first-order perturba-tion of the wave-function, our formulation is variationalwith respect to the first-order perturbation of the elec-tronic charge density. F. Approximated force-constants functional.
At this point we can exploit the results of the pre-ceding section in order to develop a method to accu-rately calculate phonon and electron-phonon properties.Within density functional theory the most precise force-constant matrix C IJ ( ω, T ) is obtained when T coincideswith the physical temperature T ( e.g. room tempera-ture) of the system. However, in a metal, this tempera-ture is prohibitively small and the self-consistent calcu-lation of C IJ ( ω, T ) would be computationally unfeasi-ble as it requires a large number of k-points in order toproperly converge the summation in Eq. 17. In addition,the self-consistent calculation of the imaginary part of n I ( r , ω, T ) and of n J ( r , ω, T ) at finite ω = 0 requires toincrease further the number of k-points with respect tostandard static linear-response calculations of n I ( r , , T ).A remedy to these problems is to define an approxi-mated force constants matrix, ˜ C IJ , as:˜ C IJ ( ω, T ) = F IJ [ n I ( r , , T ph ) , n J ( r , , T ph ) , ω, T ] (21)that is, ( i ) the frequency dependence on n I ( r , ω, T ph ) isneglected and its static limit ( ω = 0) is considered and(ii) n I and n J are calculated at the temperature T ph in-stead of T . The quantity T ph is the electronic temper-ature commonly used in a self-consistent linear responsecalculation. Usually T ph is much larger then the phys-ical temperature T , ( e.g. T ph ≈ . The mainconsequence of the variational property of the functional F IJ is that it allows us to compute phonon frequenciesfrom an approximate force-constants matrix ˜ C IJ ( ω, T )in a non self-consistent way and performing an error thatis quadratic in | n I ( r , ω, T ) − n I ( r , , T ph ) | .Thus, a time consuming self-consistent calculation ofthe non-adiabatic force-constants C IJ ( ω, T ) has beenreplaced by an approximate non self-consistent one,˜ C IJ ( ω, T ), with an error that is negligible for the phononfrequencies, as it will be demonstrated in the applicationsconsidered in this work.We now first add and subtract C IJ (0 , T ph ) (the stan-dard linear-response adiabatic self-consistent force con-stants) in the left member of Eq. 21, and then perform a Fourier transform to obtain:˜ C sr ( q , ω, T ) = Π sr ( q , ω, T ) ++ C sr ( q , , T ph ) (22)where Π sr ( q , ω, T ) == 2 N k ( T ) N k ( T ) X k ij f k i ( T ) − f k + q j ( T ) ǫ k i − ǫ k + q j + ω + iη × d sij ( k , k + q ) d rji ( k + q , k ) − N k ( T ph ) N k ( T ph ) X k ij f k i ( T ph ) − f k + q j ( T ph ) ǫ k i − ǫ k + q j × d sij ( k , k + q ) d rji ( k + q , k ) (23)with the following definition of the deformation potentialmatrix element d smn ( k + q , k ) = h k + q m | δv SCF δ u q s | k n i (24)and where | k n i is the periodic part of the Bloch wave-function, i.e. | ψ k n i = e i k · r | k n i / √ N k , and u q s is theFourier transform of the phonon displacement u Ls . Theintegration in Eq. 24 is understood to be on the unitcell. The quantity v SCF is the periodic part of the staticself-consistent potential, namely V SCF ( r ) = v SCF ( r ) e i q · r .This static self-consistent potential is calculated usingstandard linear response at temperature T ph .The calculation of Π sr ( q , ω, T ) requires the knowledgeof band energies and eigenfunctions on a denser N k ( T )k-point grid using an electronic temperature T only inan energy window of width ∼ max( T ph , ω ) around theFermi level. Moreover it does not require to recalculatethe derivative of the self-consistent potential. As suchit is much less time consuming then an linear-responseself-consistent calculation to obtain C rs ( q , , T ). Thenon-adiabatic phonon frequencies are calculated (in theclean limit) non self-consistently from the Hermitian partof ˜ C sr ( q , ω, T ) (see Eq. 5).The expression of the adiabatic force-constants at thephysical temperature T is obtained setting ω = 0 in theEq. 22: ˜ C sr ( q , ω = 0 , T ) = Π sr ( q , , T ) ++ C sr ( q , , T ph ) (25)In this case the force constants are Hermitian and theerror in the approximated force constants is quadratic in | n I ( r , , T ) − n I ( r , , T ph ) | .The advantage of the present procedure is that thelinear-response self-consistent calculation is performedwith a small number of k-points N k ( T ph ) whereas thelow temperature force constants are obtained with a nonself-consistent calculation over much denser N k ( T ) k-points mesh. This approach requires, as in a conventionalelectron-phonon coupling calculation, the knowledge ofthe wavefunctions in a dense N k ( T ) k-points grid andin an energy window of the order of the maximum be-tween T ph and ω around the Fermi level, but does notrequire the self-consistent linear-response calculation ofthe derivative of the Kohn-Sham potential with respectto an atomic displacement. As such the procedure issubstantially less time consuming. G. Phonon frequencies interpolation over k-pointsat fixed phonon-momentum: practicalimplementation
The practical implementation of the above theoreticalformulation for the phonon frequencies proceeds as fol-lows:1. Perform a standard linear response calculation fora given phonon momentum q using a N k ( T ph ) k-points mesh and smearing T ph for the electronicintegration to obtain C sr ( q , , T ph ).2. In order to perform the second summation inΠ sr ( q , ω, T ) (Eq. 23), calculate deformation po-tential matrix element of Eq. 24 on the electron-momentum grid composed of N k ( T ph ) k-points us-ing a smearing T ph .3. Generate wavefunctions | k n i and energies ǫ k n onthe denser N k ( T ) electron-momentum k-pointsgrid.4. Perform a second non self-consistent calculation ofthe deformation-potential matrix-element on themore dense N k ( T ) k-points grid using smearing T and v SCF obtained on the N k ( T ph ) grid to obtainthe first summation in Π sr ( q , ω, T ).5. Calculate ˜ C sr ( q , ω, T ) (or ˜ C sr ( q , ω = 0 , T )) usingEq. 22 (or Eq. 25).6. Diagonalize the dynamical matrix to obtain phononfrequencies and phonon eigenvectors.The procedure illustrated in this section is used toobtain at a fixed phonon momentum q well convergedphonon frequencies with respect to electronic k-pointsand smearing. Then standard Fourier interpolation canbe performed to obtain dynamical matrices throughoutthe BZ. However, this procedure, has still two main com-putational shortcomings. First, the calculations of thewavefunctions and matrix elements (points 3 and 4) be-come cumbersome when very dense N k ( T ) k-points gridfor the electron-momentum are necessary to converge.Second, in presence of phonon anomalies and not-too-smooth phonon dispersion a dense k-sampling of thephonon BZ is required.An optimal solution to overcome these problems isrepresented by the use of maximally localized Wannierfunctions as an alternative electron basis function. The method and the implementation of the Wannierfunctions approach for the electron-phonon interpolationis presented in Sec. III. H. Phonon-linewidth and of the electron-phononcoupling
The imaginary part of the force-constants matrix is re-lated to phonon-damping, namely the energy-conservingdecay of a phonon in particle-hole pairs.The imaginary part of Eq. 16 is:Im( C IJ ( ω, T )) = 2 N k ( T ) X k i, k ′ j ( f k i ( T ) − f k ′ j ( T )) × (cid:18) π h ψ k ′ j | Re (cid:20) δV SCF ( r , ω ) δ u I (cid:21) | ψ k i ih ψ k i | δV ext ( r ) δ u J | ψ k ′ j i× δ ( ǫ k i − ǫ k ′ j + ω )+ h ψ k ′ j | Im (cid:20) δV SCF ( r , ω ) δ u I (cid:21) | ψ k i ih ψ k i | δV ext ( r ) δ u J | ψ k ′ j i× P (cid:20) ǫ k i − ǫ k ′ j + ω (cid:21)(cid:19) (26)where P indicates the principal part. To obtain Eq. 26we have assume that the unperturbed Hamiltonian istime-reversal symmetric so that a real set of eigenfunc-tions exists.This expression is exact but it has two disadvantages.First it requires the expensive calculation of the deriva-tive of the self-consistent potential at finite frequencyboth in its real and imaginary parts. Second the phonon-decay in electron-hole pairs is not manifest in the termincluding the principal part. These problems persist ifthe imaginary part of Eq. 17 is considered, since, inthis functional, both δV SCF ( r , ω ) /δ u I and n I ( r , ω, T ) arecomplex quantities. These two shortcomings are absentif the approximated force-constants matrix in Eq. 21 isused. Indeed the imaginary part of Eq. 21 is:Im( ˜ C IJ ( ω, T )) = 2 π N k ( T ) X k i, k ′ j ( f k i ( T ) − f k ′ j ( T )) ×h ψ k ′ j | δV SCF ( r , δ u I | ψ k i ih ψ k i | δV SCF ( r , δ u J | ψ k ′ j i× δ ( ǫ k i − ǫ k ′ j + ω ) (27)with an error quadratic in | n I ( r , ω ) − n I ( r , | . Fouriertransforming and replacing Eq. 27 in Eq. 9 gives the wellknown equation for the phonon linewidth γ q ν full-widthhalf-maximum (FWHM), that is: γ Fermi q ν = 4 πN k ( T ) N k ( T ) X k ,m,n | g νnm ( k , k + q ) | × ( f k n − f k + q m ) δ ( ǫ k + q m − ǫ k n − ω q ν ) , (28)where g νnm ( k , k + q ) = P s e s q ν · d smn ( k + q , k ) / p M s ω q ν .In absence of the electron-electron interaction a similarexpression involving only the derivative of the bare exter-nal potential is easily obtained using Fermi golden rule.In the presence of electron-electron interaction our for-malism justify the replacement of the derivative of thebare external potential by the derivative of the staticscreened potential V SCF ( r ,
0) in the Fermi golden rule,with a well controlled error.Under certain conditions illustrated in Refs. 29,30, thephonon linewidth can be reduced to:˜ γ q ν = 4 πω q ν N k ( T ) N k ( T ) X k ,m,n | g νnm ( k , k + q ) | × δ ( ǫ k n ) δ ( ǫ k + q m − ǫ k n − ω q ν ) (29)Finally, neglecting ω q ν in Eq. 29, we obtain Allenformula namely, γ q ν = 4 πω q ν N k ( T ) N k ( T ) X k ,n,m | g νnm ( k , k + q ) | δ ( ε k n ) δ ( ε k + q m )(30)The Allen-formula is widely used, and represents agood estimation of the phonon-linewidth due to electron-phonon effects, in the absence of anharmonic effects andother scattering processes. In addition, the Allen-formularelates the phonon-linewidth, as measured by inelasticNeutron or X-ray measurements, to the electron-phononcoupling as: λ q ν = γ q ν πω q ν N s , (31)where N s is the electronic density of states per spin atthe Fermi energy.With this definition, the isotropic Eliashberg-functionis α F ( ω ) = 12 N q X q ν λ q ν ω q ν δ ( ω − ω q ν ) (32)where N q is the number of phonon wavevectors. We alsodefine the integrated electron-phonon coupling as λ ( ω ) = 2 Z ω α F ( ω ′ ) ω ′ dω ′ . (33) III. WANNIER FUNCTIONSA. Wannier interpolation of the electron-phononmatrix element
In this section we explain how to interpolate theelectron-phonon matrix element throughout the BZ usingWannier functions .Using standard first-principles methods a set of Blochfunctions ψ k n are generated with k belonging to a uni-form N wk k-points grid centered in Γ. A set of Wannier functions centered on site R are defined by the relation | R m i = 1 p N wk X k n e − i k · R U nm ( k ) | ψ k n i (34)A suitable transformation matrix, U mn ( k ), must be de-termined in the so-called Wannierization procedure. Inthis work we choose to work with Maximally LocalizedWannier functions (MLWF) , although other Wan-nierization schemes are possible. In the MLWF case, thematrices U nm ( k ) are obtained following the prescriptionof Refs. 14,28, which guarantees the maximal space lo-calization of the final Wannier functions. This last re-quirement will be essential in the spirit of interpolationof the electron-phonon matrix elements.Specifically, if the band-structure is formed by a com-posite group of bands and the number of Wannier func-tions is identical to the number of bands of the compositegroup of bands, then the matrix U mn ( k ) is a square ma-trix. On the contrary if the desired bands are entangledin a larger manifold, a preliminary disentanglement pro-cedure from the manifold must be performed and thenthe square matrix U mn ( k ) is obtained .The deformation potential matrix-elements d smn ( k + q , k ) (see Eq. 24) are calculated using linear responsescheme. In this calculation particular care is neededfor the deformation potential at zone center, namely d smn ( k , k ), as explained in sec. III B.It is crucial to note that the periodic parts | k n i and | k + q m i entering in d smn ( k + q , k ) have to be exactly thesame wavefunctions used for the Wannierization proce-dure . If this is not the case, spurious (unphysical) phasesappear in the | k n i ’s ( i.e. a phase added by the diago-nalization routine or other computational reasons) andthe localization properties of the Wannier functions inreal space are completely lost. A way to enforce thiscondition is to use a uniform grid centered at Γ so that k + q = k ′ + G with k ′ still belonging to the originalgrid. In this way | k + q n i can be obtained from | k ′ n i sim-ply multiplying by a phase determined by the G -vectortranslation and no additional phases occur.Exploiting translational invariance, the deformation-potential matrix-element in the Wannier function basisis obtained by Fourier transform as, d smn ( R , R L ) = h m | δV SCF δ u s,L | R n i = 1 N wk N wk X k , q X m ′ ,n ′ e − i k · R + i q · R L × U ∗ mm ′ ( k + q ) d sm ′ n ′ ( k + q , k ) U n ′ n ( k ) (35)where R and R L belong to a N wk real-space supercell.At this point it is important to underline that the gridof N wk k − points on which the Wannierization processhas been carried out has to be exactly the same grid forthe phonon-momentum on which the dynamical matrices have been computed . One linear response calculation foreach k-point in the phonon irreducible BZ (IBZ) of the N wk electron k-points grid has to be carried out. If thisconstraint is not respected, the localization properties ofthe electron-phonon matrix element in real space is notguaranteed.Inverting Eq. 34, we obtain | ψ k n i = 1 p N wk X R X m e i kR U ∗ nm ( k ) | R m i (36)Noting that d smn ( k + q , k ) = N wk h ψ k + q m | δV SCF δ u q s | ψ k n i , andusing Eq. 36 one obtains: d smn ( k + q , k ) = 1( N wk ) X L X R X m ′ n ′ e i k · R + i q · R L U m ′ m ( k + q ) d sm ′ n ′ ( R , R L ) U ∗ nn ′ ( k ) (37)where R and R L belong to a N wk real-space supercell.Now, if d smn ( R , R L ) is localized inside the N wk real-spacesupercell then the interaction between different N wk real-space supercells can be neglected and d smn ( k + q , k ) canbe obtained from Eq. 37 via a slow Fourier transform. Inpractice this means that in Eq. 37 now k and q are any k-points in the BZ. As a result of the interpolation, d smn ( k + q , k ) can be used to calculate any physical property asin a simple tight-binding scheme. B. Deformation potential matrix elements ofoptical zone-center phonons
In this section we discuss the peculiarities related tothe calculation of the electron-phonon matrix elementsfor optical zone center phonons ( q = ). The peri-odic part of the self-consistent potential, induced bythe phonon displacement at wavelength q of the atom s can be decomposed in the Coulomb (Cl) potential (thesum of the bare and Hartree potentials) and exchange-correlation (XC) contributions: δv SCF ( r ) δ u q s = δv XC ( r ) δ u q s + δv Cl ( r ) δ u q s (38)The integral of the Coulomb contribution over the unitcell volume Ω is:∆ q s = 1Ω Z d r δv Cl ( r ) δ u q s . (39)In DFT linear response codes, such as QUANTUM-ESPRESSO , the phonon calculation with q = and q = are treated with two different approaches. At q = , the calculation is performed within the gran-canonical ensemble, with a constant electron chemical-potential (Fermi level) ǫ F . In a metal, the limitlim q → ∆ q s = ∆ + s (40) is well defined and independent on the direction of q . Ingeneral one has ∆ + s = 0.At q = , the calculation is performed in the canonicalensemble, with a constant number of electrons. In thiscalculation the variation of the average Coulomb poten-tial is conventionally set to zero, namely:∆ s = 0 (41)To keep constant the number of electrons the deriva-tive of the Fermi energy with respect to the phonondisplacement δǫ F /δ u0 s can be different from zero. Inparticular it must be: δǫ F δ u = − ∆ + s (42)In the QUANTUM-ESPRESSO implementation the δv SCF ( r ) /δ u q s stored-on-disk is discontinuous at q = ,since, in this case, it does not include (as it should) thecontribution from the average Coulomb potential. As aconsequence the electron phonon matrix elements com-puted with δv SCF ( r ) /δ u s are incorrect. Moreover suchdiscontinuity deteriorates the localization properties ofthe electron-phonon matrix elements in real space. Theproblem is easily solved by redefining the self-consistentpotential to be used in the calculation of the deformation-potential matrix-elements for q = as: δ ˜ v SCF ( r ) δ u s = δv SCF ( r ) δ u s − δǫ F δ u s . (43)Note that since δǫ F /δ u s transforms under symmetry op-eration as a force, it is different from zero only for theatoms for which the internal coordinates are not fixed bysymmetry. For this reason, δǫ F /δ u s is finite in CaC ,while it vanishes in MgB .We remark that the same problem occurs when afrozen-phonon calculation of δv SCF ( r ) δ u s is carried out as inall the electronic-structure codes the average Coulombpotential is conventionally set to zero. C. Wannier interpolation of the dynamical matrixat fixed phonon-momentum.
Steps 4 , , q the quantity Π sr ( q , ω, T ) to obtain Eq. 23. The ad-vantage of the Wannier-interpolation-based strategy withrespect to the procedure explained in sec. II G is that thetime needed to calculate the matrix element is now neg-ligible as there is no need to obtain, from first-principles,wavefunctions and energies on the dense N k ( T ) k-pointgrid. Still it is necessary to perform Fourier interpolationto obtain dynamical matrices throughout the full BZ. D. Wannier interpolation of the dynamical matrix.
In what follows we outline a strategy to obtain dynam-ical matrices throughout the BZ using both Fourier inter-polation and Wannier functions. The method solves themain problem related to the presence of Kohn-anomalies,originating from long range interactions in the dynamicalmatrix. The idea is to separate in the force-constant ma-trix obtained in Eq. 22 in the short and long range com-ponents. The long range force constants are associatedto Kohn anomalies driven by Fermi surface nesting, and,as such, they cannot be easily Fourier interpolated andrequire an accurate sampling of the Fermi surface. Thislast part will be treated using the Wannier interpolationscheme explained in section III A. On the contrary shortrange force constants produce smooth phonon dispersionsand can be safely Fourier interpolated and do not requireaccurate sampling of the Fermi surface. Treating sepa-rately the two length scales permits to obtain convergedphonon frequencies with respect to k-point sampling andelectronic temperature, everywhere in the Brillouin zone.This procedure is similar to the treatment of long-range dipolar forces due to the Born effective chargesin polar semiconductors . In this case, these long rangeforces generate a non-analytical behavior of the dynam-ical matrix at zone center. While Fourier interpolationof the dynamical matrices including the non-analyticalterms is impossible, it becomes possible if the long-range forces are subtracted to obtain short-range force-constants. The non-analytical part subtracted to per-form Fourier interpolation is added again after Fourierinterpolation.In this spirit, we rewrite Eq. 22 as,˜ C sr ( q , ω, T ) = Π sr ( q , ω, T ) − Π sr ( q , , T ∞ ) + ˜ C sr ( q , , T ∞ ) (44)where ˜ C sr ( q , , T ∞ ) = C sr ( q , , T ph ) ++ Π sr ( q , , T ∞ ) (45)and T ∞ is an electronic temperature large enough inorder to have only short range force constants named˜ C sr ( q , , T ∞ ), and no Kohn anomalies in the correspond-ing phonon branches. As a consequence, Fourier interpo-lation can be applied to ˜ C sr ( q , , T ∞ ).In practical implementation the procedure is the fol-lowing:1. Obtain C sr ( q , , T ph ) from self-consistent linear-response phonon calculations. Such force constantsare evaluated using a N wk k-point grid for thephonon momentum and an N k ( T ph ) k-point gridfor the electron momentum.2. Calculate Π sr ( q , , T ∞ ) using the Wannier func-tions as described in section III C. 3. Use Eq. 45 to compute ˜ C sr ( q , , T ∞ ) on the N wk phonon momentum grid.4. Fourier interpolate ˜ C sr ( q , , T ∞ ) to obtain dynam-ical matrices at any desired phonon momentum attemperature T ∞ .5. Obtain ˜ C sr ( q , ω, T ) in Eq.44 by Wannier interpo-lation of Π sr ( q , ω, T ) − Π sr ( q , , T ∞ ) at any desiredphonon momentum.6. Diagonalize the dynamical matrix associated to˜ C sr ( q , ω, T , T ph ) to obtain adiabatic, ω = 0, ornon-adiabatic (in the clean limit) ω = ω q ν , phononfrequencies. These phonon frequencies includeKohn-anomalies as long range terms are properlytreated by the Wannier-functions approach. IV. APPLICATIONS
We applied the theoretical and computational frame-work developed in the preceding sections to two wellknown superconductors, namely magnesium diboride(MgB ) and Calcium intercalated graphite (CaC ). Inthese systems phonon, electron-phonon and supercon-ducting properties are extensively investigated from boththeoretical and experimental point of view and a detailedcomparison can me made. At the same time, they rep-resent non-trivial examples, where convergence problemsin first-principles approaches are relevant for the calcu-lation of phonon and superconducting properties. A. MgB Magnesium diboride (MgB ) is, without doubts, oneof the most studied and investigated materials, since2001. The discovery of superconductivity at interme-diate temperatures (40 K) , in this compound, repre-sented an unexpected surprise in the scientific commu-nity. Easily available, a huge amount of experimentalcharacterizations were immediately possible . On theother hand, MgB being composed of light elements inan hexagonal unit cell with only three atoms, alloweda first-principles computational approach for the theo-retical description and understanding of its normal andsuperconducting properties. Density Functional Theoryand linear response techniques, revealed an extraordinaryhigh electron-phonon coupling between sp bonding andanti-bonding boron orbitals, with a particular in-planemode ( E g ) of the hexagonal boron sheet . Basedon model calculations of the electron-phonon couplingparameter ( λ ≃ λ ,were too low to explain this high critical temperature,unless unphysical low Coulomb pseudo-potentials were0used . The solution of this inconsistency, has come withthe discovery of two-band-superconductivity in MgB .Two different gaps develop on the σ and π bands of theboron sheets . This allows to gain new scattering chan-nels to enhance T c still containing the average value of λ . At the moment, MgB is classified as a two-bandelectron-phonon superconductor .This continuous improvement of the theoretical expla-nation of the normal and superconducting properties ofMgB was in part due to some difficulties in the propercalculations of the electron-phonon coupling and of thephonon frequencies. Soon after the discovery of super-conductivity in MgB , the first theoretical papers (basedon standard DFT codes) evidenced a relevant and anoma-lous difficulty to properly ”converge” the E g phononfrequency and the electron-phonon coupling. The rea-son is now clear: due to the two dimensional nature(almost cylindrical) of the σ ∗ Fermi surfaces electron-phonon coupling calculations require a careful integrationof the Fermi surface is needed. This translates, in prac-tical calculations, with the necessity to use an extremelydense k- space integration for electronic and phononicdegree of freedom . Although feasible in principle, itwould require a huge amount of resources and computa-tional time.So, still after about ten years of calculations, the funda-mental question about the value of the electron-phononcoupling in MgB , remains unanswered. This informa-tion is fundamental in order to validate the electron-phonon origin of the superconductivity in MgB and atthe same time, include other, still unexplored pairing in-teractions.The Wannnier interpolation scheme, developed inthe present paper, of phonon frequencies and electron-phonon matrix elements, represents the main tool tosolve this problem. It retains the first-principles char-acter, but at the same time it is affordable from a com-putational point of view.In MgB the BZ-center phonon-frequency, measuredby Raman spectroscopy, lies between the adiabatic andclean-limit-non-adiabatic phonon-frequency computedby DFT. Thus MgB is in the intermediate regime (seeSec. II B and the relaxation time estimated in Ref. 12).Since the experimental frequency at BZ center is closer tothe adiabatic than the clean-limit-non-adiabatic value, in MgB we will not consider non-adiabatic effects, andwe will present only the adiabatic dispersion.
1. Technical details
Wannier-interpolation of the MgB band-structure iscarried out using a N wk = 6 × × Γ K M G A L H A-2-1.5-1-0.500.511.52 E n e r gy ( e V ) FIG. 1: (color online) Wannier-interpolated (red-dashed) andfirst-principles-calculated bands (black-continuous) for MgB .The band-energies are plotted with respect to the Fermi level. trix element are calculated using linear-response onthe irreducible part of the N wk k-points grid centered atΓ. The δv SCF /δ u q s on the full grid are obtained ap-plying symmetry operations of the crystal. In the linearresponse calculation we perform electronic integration us-ing a N k = 16 × × k − point grid in the Brillouin zoneand an Hermite-Gaussian smearing T ph = 0 . T ∞ =0 .
11 Ryd to obtain smoothed “high temperature” phononfrequencies.
2. Adiabatic phonon dispersion
The Wannier-interpolated adiabatic phonon dis-persions on different k-point grids and smearingsare compared to state-of-the-art electronic structurecalculations and experimental data in Fig. 2.In the top panel of Fig. 2 we compare the Wannier in-terpolated phonon dispersion with T = 0 .
02 Ryd and N k ( T ) = 30 with linear response calculations per-formed on selected points in the Brillouin zone using thesame mesh and the same electronic temperature. Theerror in the interpolation scheme is always smaller then0 . k ( T ) = 45 at T = 0 .
01 Ryd with Fourier interpolated linear responsecalculations from a 6 × × A usingour new method but are completely missed by standardFourier interpolation. The main differences with respectto previous calculations are :1. a prominent softening (Kohn anomaly) of the twoE g modes at q ≈ . Γ K M Γ A LR (N k =30 , T ph =0.02)Wannier (N k =30 , T =0.02 Ryd) E n e r gy ( m e V ) Γ K M Γ A LR (N k =16 x12, T ph =0.025)ExperimentWannier (N k =45 , T ph =0.01) E n e r gy ( m e V ) FIG. 2: (color online) Top panel: Wannier-interpolated adi-abatic phonon dispersion (red line) compared to standardlinear-response calculations performed on selected points inthe Brillouin zone (black dots). Bottom-panel: Comparisonbetween state-of-the-art calculated phonon dispersion usingFourier interpolation (black line) and Wannier interpolation(red line). Both interpolation schemes start from the same linear response calculation. Experimental data in the bot-tom panel are from Ref. 5 (blue circles) and from Ref. 44(black circles). Electronic temperatures are in Rydberg. Theacoustic sum rule is not applied. Interpolated branches areobtained by connecting the points at which the interpolationhas been carried out in the order of increasing energy. the zone center. Its origin is due to a in-planenesting between the two σ cylinders. This Kohnanomaly could be probably inferred from previ-ous calculations on smaller k-points grids al-though the softening of the E g modes in theseworks is weaker.2. A second Kohn anomaly, absent in previouscalculations , is present on the E g and B g modes along the k z direction at q ≈ . . The data at q ≈ . even if this disagreement wasoverlooked in all previous publications.3. Overall the Wannier-interpolated phonon disper-sion is in better agreement against experimentswith respect to the linear-response calculated one.These differences point out the need of having an ultra-dense k-point sampling of the Brillouin zone not onlyfor the electron-phonon coupling but also for the phonondispersion.
3. Electron-phonon coupling
Having interpolated the electron-phonon matrix ele-ments and the phonon frequencies it is possible to calcu-late the phonon-linewidth and the electron-phonon cou-pling with a high degree of precision, eliminating sourceof errors coming from insufficient k − point sampling. InFig. 3 the Wannier-interpolated Eliashberg-function iscompared with previous calculations done with differentapproaches. There are two main improvements due tothe better k-point sampling.First, several previous calculations on smallerk-point grids generally overestimate the position of themain peak in α F ( ω ) (centered around ω ≈
63 meVin our work) due to the coupling of σ -bands with theE g phonon modes. This happens because in previouscalculations the Kohn anomalies of the E g branchesat q ≈ . K was underestimated and because theanomaly along Γ A is missing. This crucial dependenceof the Eliashberg function on the E g phonon frequencydemonstrates the need of having both an accurate deter-mination of λ and of the phonon frequencies to describesuperconducting properties.Second, in our work there is substantial more weightthen in previous works in the 90-100 meV region.This energy region is dominated by the coupling of π electrons to the E g vibrations.We now compare in more details our work with previ-ous calculations. In the work of Bohnen et al. all thepeaks are up-shifted not only by the limited k-points sam-pling but also by the choice of using the DFT-minimizedcrystal structure and not the experimental one. Kongand coworkers have calculated λ q ν and ω q ν on a 6 × × q of this grid be-longing to the ΓA direction the corresponding λ q ν val-ues were replaced with some of the 12 × × σ electronic states andthe E g mode. Indeed the average electron-phonon cou-pling is much larger than the majority of the calcula-2 ω (meV) α F ( ω ) α F ( ω ) λ ( ω) FIG. 3: (color online) Comparison of our Wannier-interpolated MgB Eliashberg function with previous calcu-lations. Kong et al. used first-principles calculations anda particular treatment for λ q E g for q along ΓA. Bohnen etal. used linear response calculations while Choi et al. per-formed frozen-phonon calculations. Finally Eiguren et al. used a different implementation of Wannier interpolation onthe electron-phonon matrix elements with 40 k-points gridfor electron-momentum and 40 k-points grid for phonon mo-mentum, but did not interpolate phonon frequencies. In thiswork we Wannier-interpolate both electron-phonon coupling(using 80 k-points grid for electron-momentum and 40 k-points grid for phonon momentum) and phonon frequencies(using 30 k-points grid for electron-momentum and 40 k-point grid for phonon momentum). The integrated λ ( ω ) fromour work is also shown in the bottom panel. tions present in literature (see Table I). Choi et al. used the LDA-minimized lattice structure but surpris-ingly obtained softer E g phonon frequencies at Γ thenall other works using experimental lattice parameters. InRef. 39 the electron-phonon coupling is calculated usingan interpolation scheme that does not assure the local-ization properties of the electron-phonon matrix elementin real space. Nevertheless the position of the main peakin α F ( ω ) is similar to what we obtain. On the con-trary the secondary peak related to coupling between π states and the E g mode is at too low energy. Finally,Eiguren et al. used a Wannier interpolation schemethat requires an explicit calculation of the wavefunctionsfor any phonon momentum. The interpolation schemewas used to calculate the average electron-phonon cou-pling but not to interpolate phonon frequencies as theywere obtained by Fourier interpolation. As a consequence Reference λ ω log (meV) N k N q Kong et al. [42] 0 .
87 62 . × Bohnen et al. [43] 0 .
73 60 . Choi et al. [39] 0 .
73 62 . ×
12 18 × et al. [11] 0 .
776 40 This work 0 .
741 58 . TABLE I: Comparison between different electron-phonon cou-pling calculations in MgB using different first principlesmethods. The label N k (N q ) indicates the number of k-pointused for electron (phonon) momentum integration. the phonon dispersion does not present the two observedKohn anomalies and the main peak of the α F ( ω ) is athigher energies respect to our work. This difference re-flects the need of having an accurate determination of thephonon frequency in order to have converged Eliashbergfunctions. This appears in the lower phonon frequencylogarithmic average obtained in our work with respect toprevious calculations as shown in table I. B. CaC Graphite intercalated compounds represent a wideclass of compounds. Due to the particular bonding prop-erties of carbon atoms, with strong and stable sp bondsand weak Van der Walls π z bonds, graphite allows incor-poration of many different atomic species between thehexagonal C sheets. Alkali metals and alkaline-earthmetals are the most common, but even H, Hg, Tl, Bi,As, O, S and halogens are reported The possibility to change the number of electrons inthe covalent C-C bonds, always stimulated the search ofsuperconductivity in intercalated graphite, but only in2005 calcium intercalated graphite (CaC ) was found tosuperconduct at temperatures as high as 11.5 K .First-principles theoretical investigations on CaC re-vealed that the superconducting critical temperaturecan be explained considering only the electron-phononmechanism , even if other pairing mechanisms wereproposed (excitonic and plasmonic).The pairing originates from the electron-phonon in-teraction between intercalant Fermi Surface and the Cout-of-plane modes and Ca in-plane modes . These lastmodes, with very low frequencies ( ≃
15 meV), accountfor about one-half of the total λ .In this section, we re-examine the phonon and electron-phonon properties of CaC , by means of the computa-tional framework based on the Wannier interpolation.In fact, from a computational point of view, the cal-culation of the energy position and q -dispersion of thelow energy Ca modes requires very high precision, notobtainable with the state-of-the-art phonon and electron-phonon calculations.3 L Γ X Γ T E n e r gy ( e V ) FIG. 4: (color online) Wannierized (red-dashed) and first-principles-calculated bands (black-continuous) for CaC . Thehorizontal dashed line marks the Fermi level. CaC is well described by the clean limit regime. In-deed, the CaC zone-center phonon frequency measuredby Raman, is significantly larger than the DFTadiabatic one but agrees very well with the DFT clean-limit-non-adiabatic frequency . For this reasonwe will discuss both the adiabatic and the non-adiabaticdispersion.
1. Technical details
Wannier-interpolation of the CaC band-structure iscarried out using a N wk = 6 × × z state on each Carbon. The Wannierized band structureas compared to first-principles calculations is shown inFig. 4.We perform one linear response calculation foreach point of the irreducible N wk k − point grid using N k ( T ph ) = 8 × × T ph = 0 .
05 Ryd. We then obtain δv SCF /δ u q s on the fullgrid applying the symmetry operations of the crystal. Atzone center, we added a Fermi level shift to the q = 0component of δv SCF /δ u q s to obtain an analytical behav-ior of at q = 0 (see section III B for more details).The interpolation of phonon frequencies is performedalong the strategy outlined in sec. III D and using T ∞ = 0 .
075 Ryd to obtain smoothed “high temperature”phonon frequencies. For the calculation of non-adiabaticphonon frequencies we chose ω in Eq. 21 as the E g phonon frequency at Γ and we also use N k ( T ) = 50 ,T = 0 .
01 Ryd and η = 0 .
01 Ryd in Eq. 23.
2. Adiabatic phonon dispersion
The phonon dispersions of CaC as obtained fromFourier-interpolation and from Wannier-interpolation areillustrated in Fig. 5. The Fourier interpolate phonondispersion is in agreement with previous linear re-sponse calculations . When compared with themore converged calculations performed using Wannier-interpolation several differences are found. The largestdeviation concerns the high-energy modes at ≈ − L . These modes involve car-bon vibrations parallel to the Graphite layers (C xy vibra-tions). There is a substantial hardening of the two high-est optical phonon-modes approaching zone-boundaryand a somewhat weaker softening close to zone center.The shape of the phonon dispersion of the two high-est optical modes along Γ L recalls the typical behav-ior of the real-part of the phonon self-energy due toelectron-phonon coupling (see for example for the caseof MgB Ref. 30, fig 3a.). The electron-phonon couplingto these optical modes comes essentially from π ∗ statesand mostly from intraband-coupling (small Fermi sur-face) so that a quite accurate BZ sampling is necessary todisplay this behavior. This explains the strong k − pointdependence of the result.Moreover the relevance of having an accurate samplingof the Fermi surface is evident from the large number ofKohn-anomalies (arrows in Fig. 5 mark the occurrenceof a Kohn-anomaly) present in the Wannier-interpolatedbranches and almost always absent in the Fourier inter-polated ones. At high energy (150-180 meV) we counta large number of anomalies, no one of which is presentin Fourier interpolated branches. These anomalies aremostly located close to the X and χ points. At low en-ergy we count 8 more softenings, 6 of which involve pointsclose to the X-point (although on different branches) andtwo occur at the χ point. In Fourier interpolated calcu-lations only a very broad anomaly is visible in the lowenergy modes, located exactly at the X-point, nothingin other energy modes. The facts that close to the Xpoint (i) several anomalies at the same vector occur ondifferent modes and that (ii) the phonon frequencies as-sociated to the anomaly are strongly dependent on theelectronic temperature suggest an origin dependent fromthe Fermi surface geometry. Indeed this is the case andthe anomaly is illustrated in the nesting vector in Fig. 6.
3. Non-adiabatic phonon dispersion
Non-adiabatic effects occur in CaC close to zonecenter . However it is unknown how strong are theseeffects far from zone center, namely what is the degree oflocalization of non-adiabatic effects around the Γ point.The question is meaningful as if non-adiabatic effects ex-tend substantially in the Brillouin zone then they canbe relevant for superconductivity and transport as theaverage electron-phonon coupling and the critical tem-4 L Γ χ X Γ T E n e r gy ( m e V ) WANIER INTERPOLATEDLINEAR RESPONSE L Γ χ X Γ T E n e r gy ( m e V ) WANIER INTERPOLATEDLINEAR RESPONSE L Γ χ X Γ T E n e r gy ( m e V ) WANIER INTERPOLATEDLINEAR RESPONSE
FIG. 5: (color online) Wannier-interpolated (red-continuous) and Fourier-interpolated (black-dashed) adiabatic phonon dis-persion in CaC . The center and right panels are zooms to the high and low energy parts of the phonon dispersion. TheFourier-interpolated phonon branches are obtained from the dynamical matrices calculated on a N wk = 6 phonon grid, using N k = 8 and Hermite-Gaussian smearing T ph = 0 .
05 Ryd. The Wannier interpolated bands are obtained on a N ˜ k = 40 andHermitian-Gaussian smearing T ph = 0 .
01 Ryd. The Kohn-anomalies are indicated by arrows.Reference λ ω log (meV) N k N q Calandra et al. [47] 0 .
83 24 . × × × × et al. [52] 0 .
84 26 .
28 8 × × × × et al. [48] 0 .
85 28.0 6 × × × × .
829 27 . TABLE II: Comparison between different electron-phononcoupling calculations for CaC using different first principlesmethods.The label N k (N q ) indicates the number of k-pointused for electron (phonon) momentum integration. perature could be affected.In Fig. 7 we plot the Wannier-interpolated non-adiabatic phonon dispersion. At Γ we obtain ω NA = 190meV ( ≈ − ), in good agreement with what foundin Ref. 12. Most interestingly the full E g branch isshifted to higher energies (particularly along Γ L , the k z direction) by the non-adiabatic effects. Even at half-wayto the zone border the correction is sizable and it can beeasily measured with inelastic X-ray or Neutron scatter-ing experiments. In the case of CaC as the high-energymodes are weakly-coupled to electrons the non-adiabaticcorrection, even if large, probably does not affect thecritical temperature. The situation can be substantiallydifferent in other superconductors so that non-adiabaticeffects have to be taken into account to describe theelectron-phonon interaction in layered superconductors.
4. Electron-phonon coupling
The Eliashberg function calculated with the adiabaticphonon frequencies and the integrated electron-phononcoupling for CaC are illustrated in Fig. 8. The situation is very similar to previous works in what concern theaverage value of the electron-phonon coupling, as shownin table II.Previous calculations performed using standard linear-response methods report similar values of λ and a 13%variation is present on ω log (see table II). Our calculationconfirms the largest value reported in literature. V. CONCLUSIONS
In this work we developed a first-principles variationalscheme to calculate adiabatic and non-adiabatic phononfrequencies in the full Brillouin zone. We demonstratedhow the method permits the calculation of phonon dis-persion curves free from convergence issues related toBrillouin zone sampling. Our approach also justify theuse of the static screened potential in the calculation ofthe phonon linewidth due to decay in electron-hole pairs.We apply the method to the calculation of the phonondispersion and electron-phonon coupling in MgB andCaC . In both compounds we demonstrate the oc-currence of several Kohn anomalies, absent in previousworks, that are manifest only after careful electron andphonon momentum integration. In MgB , the presence ofKohn anomalies on the E g branches improves the agree-ment with measured phonon spectra and affects the po-sition of the main peak in the Eliashberg function. InCaC we show that the non-adiabatic effects on in-planecarbon vibrations are not localized at zone center butare sizable throughout the full Brillouin zone. Thus NAeffects can be relevant for the determination of the su-perconducting properties of layered superconductors. Ingeneral our work underlines the important of obtainingaccurate (well converged) phonon frequencies to deter-mine the electron-phonon coupling and superconducting5 FIG. 6: (Color online) Fermi surface of CaC calculated using Wannier-functions on a 100 k-points grid. The green arrowsmark the nesting vectors responsible for the Kohn-anomaly close to X. properties in a reliable way.The large gain in the computational time neces-sary to obtain phonon dispersion provided by our the-oretical scheme opens new perspectives in large-scalefirst-principles calculations of dynamical properties andelectron-phonon interaction. Furthermore, the varia-tional property of the force constants functional withrespect to the first-order perturbation of the electroniccharge density is a property that can be generalized tothe calculation of other response functions such as op-tical properties and the calculation of Knight shifts asmeasured in nuclear magnetic resonance . VI. ACKNOWLEDGEMENTS
We acknowledge fruitful discussion with J. Yates andsupport from ANR PNANO-ACCATTONE. GP ac-knowledge support from CNRS. Calculations were per-formed at the IDRIS supercomputing center and atCineca (Bologna, Italy) through an INFM-CNR super-computing. J. M. Ziman, Electron and Phonons (Clarendon Press, Ox-ford, 2001). G. Grimvall,
The electron-phonon interaction in metals ,(North Holland, Amsterdam, 1981) J. R. Schrieffer, Theory of Superconductivity (AddisonWesley, 1988). R. E. Peierls, Quantum theory of solids (Oxford UniversityPress, 2001) Abhay Shukla, Matteo Calandra, Matteo d’Astuto,Michele Lazzeri, Francesco Mauri, Christophe Bellin,Michael Krisch, J. Karpinski, S. M. Kazakov, J. Jun, D.Daghero, and K. Parlinski, Phys. Rev. Lett. , 095506(2003) M. Born and J. R. Oppenheimer, Ann. Physik , 457(1927) S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, L Γ χ X Γ T E n e r gy ( m e V ) NON ADIABATICADIABATIC
FIG. 7: (color online) Wannier-interpolated non-adiabatic(red-continuous) and adiabatic (black-dashed) phonon disper-sion in CaC . α F( ω ) λ ( ω ) ω (meV) FIG. 8: (color online) Wannier-interpolated Eliashberg func-tion and integrated Eliashberg functions for CaC . Note thatthe wiggles in the 80-150 meV region are not due to poor con-vergence but to the presence of many well separated dispersivebranches in this energy region.Rev. Mod. Phys. , 515 (2001) S. Y. Savrasov and D. Y. Savrasov, Phys. Rev. B , 16487(1996) X. Gonze, C. Lee, Phys. Rev. B 55, 10355 (1997) F. Giustino, M. L. Cohen, and S. G. Louie, Phys. Rev. B , 165108 (2007) A. Eiguren and C. Ambrosch-Draxl, Phys. Rev. B ,045124 (2008). A. M. Saitta, M. Lazzeri, M. Calandra, and F. Mauri,Phys. Rev. Lett. , 226401 (2008) F. Giustino, J. R. Yates, I. Souza, M. L. Cohen, and S. G.Louie, Phys. Rev. Lett. , 047005 (2007) N. Marzari and D. Vanderbilt, Phys. Rev. B , 12847-12865 (1997) M. P. M. Dean, C. A. Howard, S. S. Saxena, and M.Ellerby, Phys. Rev. B , 045405 (2010) see e.g. N. W. Ashcroft and N. D. Mermin, Solid StatePhysics (Holt, Rhinehard and Winston, New York, 1976),Chap. 26, page 529, footnote 31. de Gironcoli, S., 1995, Phys. Rev. B , 6773. A. A. Quong and B. M. Klein, Phys. Rev. B 46, 10734(1992) S. Y. Savrasov, Phys. Rev. Lett. 69, 2819 (1992) E. Maksimov and S. V. Shulga, Solid State Communication , 553 (1996) M. Calandra, M. Lazzeri and F. Mauri, Physica C , 38(2007) E. Cappelluti, Phys. Rev. B 73, 140505 (2006). Hellmann, H., 1937,
Einf¨uhrung in die Quantenchemie (Deuticke, Leipzig). Feynman, R. P., 1939, Phys. Rev. , 340. M. Di Ventra and S. T. Pantelides, Phys. Rev. B , 16207(2000) X. Gonze, D. C. Allan and M. P. Teter, Phys. Rev. Lett. , 3603 (1992). In principle, the force constants matrix depends implicitlyon the temperature used in the self-consistent calculationof n ( r ), ǫ k ,i and ψ k ,i ( r ) 17. Moreover, the dependence offorce constants on this last temperature is weak and thusit is usually taken (as in the present derivation) equal to T ph . However, it should be set equal to the physical tem-perature, T as a single standard self-consistent calculationdoes not poses relevant problems in k-point convergence. N. Marzari and D. Vanderbilt, Phys. Rev. B , 035109(2001) P. B. Allen, Phys. Rev. B , 2577 (1972), P. B. Allen andR. Silberglitt, Phys. Rev. B , 4733 (1974). Matteo Calandra and Francesco Mauri, Phys. Rev. B ,064501 (2005) P. Giannozzi et al. J. Phys.: Condens. Matter 21, 395502(2009) J. Nagamatsu, N. Nakagawa, T. Muranaka, Y. Zenitani,and J. Akimitsu, Nature , 63 (2001). Recent Advances in MgB2 Research, Physica C (2007). J. Kortus, I. I. Mazin, K. D. Belashchenko, V. P. Antropov,and L. L. Boyer, Phys. Rev. Lett. , 4656 (2001). J. M. An and W. E. Pickett , Phys. Rev. Lett. , 4366(2001). Y. Kong, O.V. Dolgov, O. Jepsen, and O.K. Andersen,Phys. Rev. B , 020501 (2001). A. Y. Liu, I. I. Mazin, and J. Kortus, Phys. Rev. Lett. ,087005 (2001). F. Giubileo et al. , Phys. Rev. Lett. , 177008 (2001). H. J. Choi et al. , , D. Roundy, H. S., M. L. Cohen, and S.G. Louie, Nature (London) , 758 (2002), H. J. Choi,D. Roundy, H. S., M. L. Cohen, and S. G. Louie , Phys.Rev. B , 020513 (2002). A. Floris et al. , Phys. Rev. Lett. , 037004 (2005). M. Calandra and F. Mauri, Phys. Rev. B , 064501(2005). Y. Kong, O. V. Dolgov, O. Jepsen, and O. K. Andersen,Phys. Rev. B, , 020501(R) (2001). K. P. Bohnen, R. Heid and B. Renker, Phys. Rev. Lett. Matteo d’Astuto, Matteo Calandra, Stephanie Reich, Ab-hay Shukla, Michele Lazzeri, Francesco Mauri, JanuszKarpinski, N. D. Zhigadlo, Alexei Bossak, and Michael Krisch, Phys. Rev. B , 174508 (2007) . N. Emery et al. , Sci. Technol. Adv. Mater. T. E. Weller, M. Ellerby, S. S. Saxena, R. P. Smith and N.T. Skipper, Nat. Phys. , 39 (2005). M. Calandra and F. Mauri, Phys. Rev. Lett. , 237002(2005) A. Sanna, G. Profeta, A. Floris, A. Marini, E. K. U. Grossand S. Massidda, Phys. Rev. B , 020511 (2007) J. Hlinka, I. Gregora, J. Pokorny, C. H´erold, N. Emery,J. F. Marˆech´e, and P. Lagrange Phys. Rev. B 76, 144512(2007) A. Mialitsin, J. S. Kim, R. K. Kremer, and G. BlumbergPhys. Rev. B 79, 064503 (2009) Matteo Calandra and Francesco Mauri, Phys. Rev. B ,094507 (2006) J. S. Kim, L. Boeri, R. K. Kremer, and F. S. Razavi Phys.Rev. B , 214513 (2006) L. Boeri, G. B. Bachelet, M. Giantomassi, and O. K. An-dersen, Phys. Rev. B , 064510 (2007). M. d’Avezac, N. Marzari and F. Mauri, Phys. Rev. B76