Density Functional Theory for Superconductors with Particle-hole Asymmetric Electronic Structure
aa r X i v : . [ c ond - m a t . s up r- c on ] M a y Density Functional Theory for Superconductors with Particle-hole AsymmetricElectronic Structure
Ryosuke Akashi and Ryotaro Arita , Department of Applied Physics, The University of Tokyo, Hongo, Bunkyo-ku, Tokyo 113-8656, Japan and JST-PRESTO, Kawaguchi, Saitama 332-0012, Japan (Dated: September 15, 2018)To extend the applicability of density functional theory for superconductors (SCDFT) to systemswith significant particle-hole asymmetry, we construct a new exchange-correlation kernel enteringthe gap equation. We show that the kernel is numerically stable and does not diverge even in the lowtemperature limit. Solving the gap equation for model systems with the present kernel analyticallyand numerically, we find that the asymmetric component of electronic density of states, which hasnot been considered with the previous kernel, systematically decreases transition temperature ( T c ).We present a case where the decrease of T c amounts to several tens of percent. PACS numbers: 74.20.-z, 74.62.Dh, 74.62.Yb
I. INTRODUCTION
Non-empirical calculation of superconducting transi-tion temperature ( T c ) has been a great challenge in theo-retical physics. Conventionally, T c s of phonon-mediatedsuperconductors are estimated by the McMillan-Allen-Dynes formula, T c = ω ln . (cid:18) − . λ ) λ − µ ∗ (1 + 0 . λ ) (cid:19) , where ω ln , λ and µ ∗ (Ref. 3) are the logarithmicallyaveraged phonon frequency, electron-phonon coupling,and Coulomb pseudo-potential, respectively. This for-mula clearly indicates the extreme sensitivity of T c tothe phonon and electron properties; to achieve a reliablecalculation of T c , we need to evaluate these factors withhigh accuracy. For the phonon properties, developmentsin methods based on density functional theory enablesone to obtain reliable values very efficiently. On the otherhand, for µ ∗ , while some calculation based on the staticrandom-phase approximation has been made, quantita-tive description of the retardation effect still remains aserious obstacle for the calculation of T c .Recently, the situation has changed by the progressin the density-functional theory for superconductors(SCDFT). There, a non-empirical scheme treatingboth the electron-phonon and electron-electron inter-actions was formulated referring the Migdal-Eliashberg(ME) theory of strong-coupling superconductivity.
Subsequent applications have shown that this schemecan reproduce experimental T c s of various conventionalphonon-mediated superconductors with error less than afew K. An important advantage of the current SCDFT-basedscheme is that the particle-hole symmetric component ofelectronic structures around the Fermi level is accuratelytreated on the basis of the density functional calculations.In the context of the ME theory, its effect on T c hasbeen intensively studied for A15 compounds. How-ever, the particle-hole antisymmetric component ignored in the current SCDFT has also been shown to have a sig-nificant impact on T c . Moreover, there has also beena great experimental progress of electric-field-inducedcarrier-doping technique, and it has stimulated currentinterest to superconductivity in doped insulators havingsignificant particle-hole asymmetry.
With this back-ground, it is important to extend the applicability ofSCDFT to systems with particle-hole asymmetric elec-tronic structure.In Ref. 9, the assumption of particle-hole symmetrywith respect to the Fermi energy ( E F ) was originally in-troduced to avoid numerical divergence in the exchange-correlation kernel. In this paper, we show that we canconstruct an exchange correlation functional for systemswith asymmetric electronic structure which does not haveany divergences. In Sec. II, we show how the kernel isconstructed, and discuss how it is numerically stable.In Sec. III, we perform calculations of T c for the caseswhere the present improvement becomes indeed impor-tant. Section IV is devoted to summary and conclusion. II. GAP-EQUATION KERNEL FORPARTICLE-HOLE ASYMMETRY
In this section, we describe the improvement of thegap equation based on the current SCDFT to include theparticle-hole asymmetry effect. In the current SCDFT, we solve the following gap equation:∆ n k = −Z n k ∆ n k − X n ′ k ′ K n k n ′ k ′ tanh[( β/ E n ′ k ′ ] E n ′ k ′ ∆ n ′ k ′ . (1)Here, n and k denote the band index and crystal mo-mentum, respectively, ∆ is the gap function, and β isthe inverse temperature. The energy E n k is defined as E n k = p ξ n k + ∆ n k and ξ n k = ǫ n k − µ is the one-electronenergy measured from the chemical potential µ , where ǫ n k is obtained by solving the normal Kohn-Sham equa-tion in density functional theory (DFT).Functions Z and K are the exchange-correlation ker-nels, which represent the effects of mass renormalizationand effective pairing interactions, respectively. In theSCDFT, Z is treated by only the phonon contribution as Z n k = Z ph n k (2)with Z ph n k representing the renormalization of electronicstates by the phonon exchange, and K consists of boththe electron-phonon K ph n k n ′ k ′ and electron-electron K el n k n ′ k ′ exchange contributions K n k n ′ k ′ = K ph n k n ′ k ′ + K el n k n ′ k ′ . (3)Formally, Z ph n k is derived from the functional derivativeof the free energy given by the Kohn-Sham perturbationtheory with respect to the anomalous density as Z ph , PT n k ≡ Z ph1 ,n k + Z ph2 ,n k (4)with Z ph1 ,n k = 1tanh (cid:2) ( β/ ξ n k (cid:3) X n ′ k ′ X ν | g n k ,n ′ k ′ ν k − k ′ | × (cid:26) ξ n k [ I ( ξ n k , ξ n ′ k ′ , ω ν k − k ′ ) − I ( ξ n k , − ξ n ′ k ′ , ω ν k − k ′ )] − I ′ ( ξ n k , ξ n ′ k ′ , ω ν k − k ′ ) (cid:27) (5)and Z ph2 ,n k = − P n ′ k ′ ( β/ / cosh (cid:2) ( β/ ξ n ′ k ′ (cid:3) × (cid:20) ξ n k − β/ (cid:2) ( β/ ξ n k (cid:3) cosh (cid:2) ( β/ ξ n k (cid:3) (cid:21) × X m l m l X ν | g m l ,m l ν l − l | I ′ ( ξ m l , ξ m l , ω ν l − l ) . (6)Here, g n k ,n ′ k ′ ν k − k ′ denotes the electron-phonon coupling, and I and I ′ are defined as I ( ξ, ξ ′ , ω )= f β ( ξ ) f β ( ξ ′ ) n β ( ω ) (cid:20) e βξ − e β ( ξ ′ + ω ) ξ − ξ ′ − ω − e βξ ′ − e β ( ξ + ω ) ξ − ξ ′ + ω (cid:21) , (7) I ′ ( ξ, ξ ′ , ω )= ∂∂ξ I ( ξ, ξ ′ , ω ) , (8)with f β and n β denoting the Fermi and Bose distributionfunctions, respectively. We note that Z ph1 originates fromexplicit dependence of the free energy on the anomalousdensity, whereas Z ph2 comes from the implicit dependencevia the chemical potential.L¨uders et al. reported that this form suffers from nu-merical divergence and that the divergent contribution isantisymmetric with respect to E F . Accordingly, in theirpaper they symmetrized Z ph1 in ξ n k , and dropped Z ph2 because of its antisymmetric dependence on ξ n k . Theresulting form is Z ph , sym n k = − (cid:2) ( β/ ξ n k (cid:3) X n ′ k ′ X ν | g n k ,n ′ k ′ ν k − k ′ | × (cid:2) I ′ ( ξ n k , ξ n ′ k ′ , ω ν k − k ′ )+ I ′ ( ξ n k , − ξ n ′ k ′ , ω ν k − k ′ ) (cid:3) . (9)Further, they modified an unphysically large componentof Z ph , sym n k around E F , and obtained the following form Z ph n k = − (cid:2) ( β/ ξ n k (cid:3) X n ′ k ′ X ν | g n k ,n ′ k ′ ν k − k ′ | × (cid:2) J ( ξ n k , ξ n ′ k ′ , ω ν k − k ′ )+ J ( ξ n k , − ξ n ′ k ′ , ω ν k − k ′ ) (cid:3) , (10)where J is defined by Eq. (80) in Ref. 9. So far, only thisform has been used in practice.In Eqs. (9) and (10), the antisymmetric parts of theelectronic density of states (DOS) and | g n k ,n ′ k ′ ν k − k ′ | with re-spect to E F within the phonon energy scale are ignored.Within the standard ME theory, the effect of their an-tisymmetric parts is assumed to be minor and oftenignored for simplicity. On the other hand, in systemswith rapidly varying DOS, this treatment is also knownto lead to a considerable error in the estimate of T c andthe isotope-effect coefficient. Below, we will show a nu-merically stable form of Z without the symmetrization,with which the asymmetry effect is properly treated. A. Cancellation of divergent terms
First, we analytically examine the divergence of Z ph , PT . For simplicity, we deal with the n k -averagedform of Z ph , PT , defined by Z ph , PT ( ξ )= Z ph1 ( ξ ) + Z ph2 ( ξ ) , (11) Z ph i ( ξ ) ≡ N ( ξ ) X n k δ ( ξ − ξ n k ) Z ph i,n k ( i = 1 , , (12)where N ( ξ ) is the electronic DOS. We also introduce anapproximation for the electron-phonon coupling X n k n ′ k ′ ν δ ( ξ − ξ n k ) δ ( ξ ′ − ξ n ′ k ′ ) δ ( ω − ω ν k − k ′ ) | g n k ,n ′ k ′ ν k − k ′ | ≃ N ( ξ ) N ( ξ ′ ) N (0) α F ( ω ) , (13)where the antisymmetric component of DOS is retained.The function α F ( ω ) denotes the Eliashberg function,with which the electron-phonon coupling coefficient λ and the characteristic frequency ω ln are defined as λ =2 Z dω α F ( ω ) ω , (14) ω ln =exp (cid:20) λ Z dω α F ( ω ) ω ln ω (cid:21) . (15)Starting from the above form for Z , in brief, we showthat the numerically divergent terms analytically cancelwith each other so that we can construct a numericallystable form. Readers who are not interested in the de-tail of this analysis can skip the remaining part of thissubsection.Let us first focus on Z ph1 . As we show below, the diver-gence in this part basically comes from the slow decay oforder O ( ξ ′− ) of the integrands. In order to extract this,we transform the following factor appearing in I and I ′ [Eqs. (7) and (8)] as1 ξ − ξ ′ ± ω = P ξ ( ξ − ξ ′ ± ω )( ξ ′ ∓ ω ) − P 1 ξ ′ ∓ ω , (16) The principal value integral P has been introduced todeal with the poles at ξ ′ = ± ω in the right hand side;note that the expression in the left hand side doesnot have these poles. Substituting this expression intoEq. (12), inserting identities 1 = R dξ ′ δ ( ξ ′ − ξ n ′ k ′ ) and1 = R dωδ ( ω − ω ν k − k ′ ), and using Eq. (13), we obtain Z ph1 ( ξ )= 1tanh[( β/ ξ ] Z dω Z dξ ′ N ( ξ ′ ) N (0) α F ( ω ) (cid:26) I ( ξ, ξ ′ , ω ) + I ( ξ, ξ ′ , ω ) − J ( ξ, ξ ′ , ω ) + J ( ξ, ξ ′ , ω ) + J ( ξ, ξ ′ , ω )] (cid:27) , (17) I ( ξ, ξ ′ , ω )= − ξ f β ( ξ ) f β ( ξ ′ ) n β ( ω ) (cid:2) P e βξ − e β ( ξ ′ + ω ) ξ ′ + ω − P e βξ ′ − e β ( ξ + ω ) ξ ′ − ω − P 1 − e β ( ξ + ξ ′ + ω ) ξ ′ + ω +P e β ( ξ + ξ ′ ) − e β ( ω ) ξ ′ − ω (cid:3) , (18) I ( ξ, ξ ′ , ω )= f β ( ξ ) f β ( ξ ′ ) n β ( ω ) (cid:20) P e βξ − e β ( ξ ′ + ω ) ( ξ − ξ ′ − ω )( ξ ′ + ω ) − P e βξ ′ − e β ( ξ + ω ) ( ξ − ξ ′ + ω )( ξ ′ − ω ) − P 1 − e β ( ξ + ξ ′ + ω ) ( ξ + ξ ′ + ω )( ξ ′ + ω ) +P e β ( ξ + ξ ′ ) − e βω ( ξ + ξ ′ − ω )( ξ ′ − ω ) (cid:21) , (19) J ( ξ, ξ ′ , ω )= f ′ β ( ξ ) f β ( ξ ′ ) n β ( ω ) (cid:2) P 1 + e β ( ξ ′ + ω ) ξ ′ + ω + P e βξ ′ + e βω ξ ′ − ω (cid:3) , (20) J ( ξ, ξ ′ , ω )= − f ′ β ( ξ ) f β ( ξ ′ ) n β ( ω ) ξ (cid:20) P 1 + e β ( ξ ′ + ω ) ( ξ − ξ ′ − ω )( ξ ′ + ω ) + P e βξ ′ + e βω ( ξ − ξ ′ + ω )( ξ ′ − ω ) (cid:21) , (21) J ( ξ, ξ ′ , ω )= − f β ( ξ ) f β ( ξ ′ ) n β ( ω ) (cid:20) e βξ − e β ( ξ ′ + ω ) ( ξ − ξ ′ − ω ) − e βξ ′ − e β ( ξ + ω ) ( ξ − ξ ′ + ω ) (cid:21) . (22)Here I i ( i =0 ,
1) and J i ( i =0 , ,
2) denote the terms origi-nating from I and I ′ , respectively. From this expression,the divergent contribution Z ph , div1 is extracted as Z ph , div1 ( ξ )= 1tanh[( β/ ξ ] Z dω Z dξ ′ N ( ξ ′ ) N (0) α F ( ω ) × (cid:26) I ( ξ, ξ ′ , ω ) − J ( ξ, ξ ′ , ω ) (cid:27) . (23)We do not explicitly write the variables ( ξ, ξ ′ , ω ) unlessnecessary in the following.Let us examine Z ph , div1 . For simplicity, we start fromthe case of constant electronic DOS N ( ξ )= N (0) withasymmetric energy cutoffs [ − L , L ] as depicted in Fig. 1,which is the simplest case of particle-hole asymmetry. Wethen get to Z ph , div1 ( ξ ) = 1tanh[( β/ ξ ] Z dωα F ( ω ) Z L − L dξ ′ × (cid:26) I − J (cid:27) . (24)Next we carry out the ξ ′ integral. Since the integrand isof order O ( ξ ′− ), integration results in cutoff-dependent DOS
L-L ξ N ( ξ ) FIG. 1: Constant DOS with asymmetric cutoff energies, L and − L . terms as Z L − L dξ ′ I ≃ ξ (cid:2) f β ( ξ ) − f β ( − ξ ) (cid:3) ln (cid:12)(cid:12)(cid:12)(cid:12) L + ωL + ω (cid:12)(cid:12)(cid:12)(cid:12) , (25) Z L − L dξ ′ J ≃ f ′ β ( ξ )ln (cid:12)(cid:12)(cid:12)(cid:12) L + ωL + ω (cid:12)(cid:12)(cid:12)(cid:12) . (26)Here we have neglected terms of order O ( e − βω ), O ( e − βL ) and O ( e − βL ), and this approximation is im-plicitly used in the following equations with “ ≃ ” in thissection. We consequently get to Z ph , div1 ( ξ ) ≃− (cid:26) ξ − β/ β/ ξ ]cosh[( β/ ξ ] (cid:27) × Z dωα F ( ω )ln (cid:12)(cid:12)(cid:12)(cid:12) L + ωL + ω (cid:12)(cid:12)(cid:12)(cid:12) . (27)The ξ -dependent factor here is reminiscent of Z ph2 [Eq. (6)]. In fact, we prove that this term exactly cancels Z ph2 ( ξ ) in the following.Next, we consider Z ph2 ( ξ ) using the approximations ofthe same level. The n k -averaged form is explicitly writ-ten as Z ph2 ( ξ )= − (cid:20) ξ − β/ (cid:2) ( β/ ξ (cid:3) cosh (cid:2) ( β/ ξ (cid:3) (cid:21) × Z dωα F ( ω ) Z L − L dξ ′ Z L − L dξI ′ ( ξ, ξ ′ , ω ) , (28)where we have used X n ′ k ′ ( β/ / cosh (cid:2) ( β/ ξ n ′ k ′ (cid:3) ≃ N (0) . (29)By carrying out the integration R dξ R dξ ′ and ignoringterms of order O [( T /L ) , ( T /L ) ], we get Z L − L dξ ′ Z L − L dξI ′ ( ξ, ξ ′ , ω ) ≃ ln (cid:12)(cid:12)(cid:12)(cid:12) L + ωL + ω (cid:12)(cid:12)(cid:12)(cid:12) , (30)and immediately obtain Z ph2 ( ξ ) ≃ (cid:20) ξ − β/ (cid:2) ( β/ ξ (cid:3) cosh (cid:2) ( β/ ξ (cid:3) (cid:21) × Z dωα F ( ω )ln (cid:12)(cid:12)(cid:12)(cid:12) L + ωL + ω (cid:12)(cid:12)(cid:12)(cid:12) . (31)This expression exactly cancels out Z ph , div1 ( ξ ) [Eq. (27)].One can observe the same cancellation for the noncon-stant DOS case. The cutoff-dependent parts Z ph , div1 and Z ph2 then read Z ph , div1 ( ξ ) ≃ β ξ Z dωa F ( ω ) Z dξ ′ N ( ξ ′ ) N (0) (cid:2) I − J (cid:3) , (32) Z ph2 ( ξ ) ≃− (cid:26) ξ − β/ β/ ξ ]sinh[( β/ ξ ] (cid:27)Z dωα F ( ω ) × Z L − L dξ N ( ξ ) N (0) Z L − L dξ ′ N ( ξ ′ ) N (0) I ′ ( ξ, ξ ′ , ω ) . (33) =+ ++ ++ + DOS
FIG. 2: Decomposition procedure of DOS. Typical DOS isrepresented as red-shaded area.
Here, we approximate the DOS as a step-like func-tion with a certain set of energy tics {± ǫ i } ( i =1 , , ..., N ǫ ; ǫ i < min[ L , L ]). The approximate functionis then decomposed into functions θ + i ( ξ ) ≡ θ ( ξ − ǫ i ) + θ ( − ξ − ǫ i ) and θ − i ( ξ ) ≡ θ ( ξ − ǫ i ) − θ ( − ξ − ǫ i ) as N ( ξ ) = N (0) (cid:26) N ǫ X i =1 [ N − i θ − i ( ξ ) + N + i θ + i ( ξ )] (cid:27) . (34)This procedure is schematically depicted in Fig. 2. For Z ph , div1 , straightforward calculations yield Z ph , div1 ( ξ ) ≃ − (cid:26) ξ − β/ β/ ξ ]sinh[( β/ ξ ] (cid:27)Z dωα F ( ω ) × (cid:26) ln (cid:12)(cid:12)(cid:12)(cid:12) L + ωL + ω (cid:12)(cid:12)(cid:12)(cid:12) + X i (cid:20) N − i ln (cid:12)(cid:12)(cid:12)(cid:12) ( L + ω )( L + ω )( ǫ i + ω ) (cid:12)(cid:12)(cid:12)(cid:12) + N + i ln (cid:12)(cid:12)(cid:12)(cid:12) L + ωL + ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:21)(cid:27) . (35)For Z ph2 , terms proportional to N ± i N ± j (any double sign)seemingly become nonzero, but careful calculations showthat these terms cancel out each other (See AppendixA). Finally, we see that the cutoff-dependent part in Z ph2 exactly cancels Eq. (35). By taking limit for the numberof tics N ǫ → ∞ , the present cancellation is proved evenfor general DOS. Now let us point out the important aspect of Eq. (11):With the present approximation, the cutoff dependent di-vergence cancels if analytically calculated, and the appar-ently divergent contribution is in fact nonsingular withrespect to ξ , since it is proportional to a nonsingularfunction ξ − β/ β/ ξ ]sinh[( β/ ξ ] . In order to verify this,we have calculated Z ph , PT ( ξ ) [Eq. (11)] with a certainasymmetric model DOS [depicted later in Sec. III]. Wehave generated the energy grid so that it becomes uni-form in logarithmic scale. The calculated result is shown [eV] -7 -8 -9 -10
10 10 10 10 -10 -6 -2 2 -10 -10 -10-10-10 -6-22 cutoff= Z ph , PT ( ξ ) ξ FIG. 3: Calculated values of Z ph , PT ( ξ ) [Eq. (11)] with the ap-proximation given in Eq. 13. Model DOS depicted in Fig. 5with E F = 0 and model Eliashberg function given by Eq. (56)with λ = 1, ω E = 400cm − , and σ = ω E /
10 were used. Thecalculations were done with T = 0 . − (redsolid line), 10 − (green dashed line), 10 − (blue dotted line),and 10 − (violet, dot-dashed line) were used. The DOS val-ues for | ξ | < T were fixed to N (0) so that it is consistentwith the condition ( ǫ i ≫ T ) assumed in Sec. II A. in Fig. 3. By lowering the minimum energy cutoff for thegrid, convergence within order of λ has been achieved,where the apparent divergence gradually vanishes. How-ever, we have also found that the convergence requiresformidably accurate integration within the energy scalesmaller than temperature, and so it is extremely difficultto achieve in practical calculations. B. Z kernel for particle-hole asymmetic systems The above analysis specifies the numerically diverging(or unstable ) but analytically irrelevant terms. In or-der to circumvent the numerical difficulty, we propose tosimply subtract Z ph , div1 and Z ph2 , obtaining Z ph , aux ( ξ )= 1tanh[( β/ ξ ] Z dωα F ( ω ) Z dξ ′ N ( ξ ′ ) N (0) × (cid:2) I − J − J (cid:3) . (36)This auxiliary expression shows stable convergence, in-cludes the effect of the antisymmetric part of electronicstates, and agrees with Eq. (11) with error of order O [( T /L ) , ( T /L ) ] in the limit where Eq. (13) becomesexact. The form given in Eq. (36) is a generalization ofthe n k -averaged form of Eq. (9) Z ph , sym ( ξ )= − (cid:2) ( β/ ξ (cid:3) Z dωα F ( ω ) Z dξ ′ N ( ξ ′ ) N (0) × (cid:2) I ′ ( ξ, ξ ′ , ω )+ I ′ ( ξ, − ξ ′ , ω ) (cid:3) : (37)Symmetrization of Eq. (36) in ξ yields Eq. (37), the lim-iting value lim ξ → [ Z ph , aux ( ξ ) − Z ph , sym ( ξ )] = 0 because of the fact that the integrand for ξ =0 is symmetric in ξ ′ , and the two forms show almost the same temperaturedependence.Following the procedure of Ref. 9, we next modify Z ph , aux ( ξ ) around E F and obtain Z ph , aux2 ( ξ )= 1tanh[( β/ ξ ] Z dωα F ( ω ) Z dξ ′ N ( ξ ′ ) N (0) × (cid:2) I − J ′ − J (cid:3) , (38) J ′ = f ′ β ( ξ ) ξ (cid:20)(cid:8) f β ( ξ − ω )+ n β ( − ω ) (cid:9) P 1( ξ − ξ ′ − ω )( ξ ′ + ω ) − (cid:8) f β ( ξ + ω )+ n β ( ω ) (cid:9) P 1( ξ − ξ ′ + ω )( ξ ′ − ω ) (cid:21) . (39)In practice, we also found that the principal-value in-tegral around the singularity at ξ ′ = ± ω is numericallyunstable and its convergence is difficult to achieve with apractical computational cost. We therefore introduce aneven smoothing function p ( ξ ′ ± ω ) whose value is unityfor | ξ ′ ± ω | & T and of order O [( ξ ′ ± ω ) ] for | ξ ′ ± ω | ≪ T .Substituting p ( ξ ′ ± ω ) / ( ξ ′ ± ω ) for P[1 / ( ξ ′ ± ω )] and re-arranging the terms yield Z ph , new ( ξ )= 1tanh[( β/ ξ ] Z dωα F ( ω ) Z dξ ′ N ( ξ ′ ) N (0) (cid:2) I − J (cid:3) , (40) I ( ξ, ξ ′ , ω )=˜ I ( ξ, ξ ′ , ω ) − ˜ I ( ξ, ξ ′ , − ω ) − ˜ I ( ξ, − ξ ′ , ω )+ ˜ I ( ξ, − ξ ′ , − ω ) , (41)˜ I ( ξ, ξ ′ , ω )=[ f β ( ξ )+ n β ( ω )] f β ( ξ ′ ) − f β ( ξ − ω ) ξ − ξ ′ − ω p ( ξ ′ + ω ) ξ ′ + ω , (42) J ( ξ, ξ ′ , ω )= ˜ J ( ξ, ξ ′ , ω ) − ˜ J ( ξ, ξ ′ , − ω ) , (43)˜ J ( ξ, ξ ′ , ω )= − f β ( ξ )+ n β ( ω ) ξ − ξ ′ − ω p ( ξ ′ + ω ) (cid:20) f β ( ξ ′ ) − f β ( ξ − ω ) ξ − ξ ′ − ω − βf β ( ξ − ω ) f β ( − ξ + ω ) ξξ ′ + ω (cid:21) . (44)We have used p ( x ) = [tanh(500 βx )] in this paper,though the calculated value is not sensitive to a specficchoice of p ( x ).With the above form, the asymmetry effect is properlytreated in a numerically stable manner; the smoothnessof the integrands at ξ ′ = ± ξ ± ω, ± ω , and ξ = 0 is re-tained. In addition, for approximately symmetric DOS,the calculated results are guaranteed to be quite close tothose with the n k -averaged form of Eq. (10) Z ph ( ξ )= − β/ ξ ] Z dωα F ( ω ) Z dξ ′ N ( ξ ′ ) N (0) × (cid:2) J ( ξ, ξ ′ , ω )+ J ( ξ, − ξ ′ , ω ) (cid:3) (45)because of the following properties: (i)lim ξ → [ Z ph , new ( ξ ) − Z ph ( ξ )] ≃
0, and particu-larly for constant DOS, the desirable behavior lim ξ → Z ph , new ( ξ ) ≃ lim ξ → Z ph ( ξ ) ≃ λ is satisfied,(ii) temperature dependence of the calculated value issimilar to that of Z ph ( ξ ), and (iii) the symmetric partof the integrand, that is, I ( ξ, ξ ′ , ω )+ I ( ξ, − ξ ′ , ω ) − J ( ξ, ξ ′ , ω ) − J ( ξ, − ξ ′ , ω )2 (46)agrees with J ( ξ, ξ ′ , ω )+ J ( ξ, − ξ ′ , ω ) in Eq. (45) when we im-pose a condition ω & | ξ |≫ T . These properties are demon-strated in the following calculations. Further, we havefound that the n k -resolved counterpart of Eq. (40) givenby Z ph , new n k = 1tanh[( β/ ξ ] X n ′ k ′ ν | g k − k ′ νn k ,n ′ k ′ | × (cid:2) I ( ξ n k , ξ n ′ k ′ , ω ν k − k ′ ) − J ( ξ n k , ξ n ′ k ′ , ω ν k − k ′ ) (cid:3) (47)is also numerically stable, which is a reasonable general-ization of Eq. (10). Equations (40) and (47) are those wepropose as the new forms of Z for particle-hole asymme-try. III. ASYMMETRY EFFECT ON T c In the rest of the present paper we get insights into theasymmetry effect included in the present improvement.In Sec. III A, we analytically construct a formula of T c to discuss in what situation the asymmetry significantlyaffects T c . In Sec. III B, a typical model case is presentedwhere our new kernel and the previous one give quitedifferent T c , and some remarks on the application to ele-mental metals are given. We here note that we focus onthe n k -averaged form [Eq. (40)], and thereby the presentscope includes only the asymmetry of DOS, not that ofthe electron-phonon matrix elements. A. Analytic calculation
We first analytically solve the energy-averaged gapequation ∆( ξ )= −Z ( ξ )∆( ξ ) − Z dξ ′ N ( ξ ′ ) K ( ξ, ξ ′ ) tanh[( β/ ξ ′ ] ξ ′ ∆( ξ ′ ) , (48)with the electron-phonon kernels only. We introduce asimple model for which DOS and kernels are given bythe constant part and the antisymmetric part, given as Z ( ξ ) = Z c ( ξ ) + Z a ( ξ ), N ( ξ ) = N c ( ξ ) + N a ( ξ ), K ( ξ, ξ ′ ) = K c ( ξ, ξ ′ ) + K a ( ξ, ξ ′ ), with each part defined within the K c ( ξ, ξ ) K a ( ξ, ξ ) Z c ( ξ ), N c ( ξ ) Z a ( ξ ), N a ( ξ ) ω D ωω D ω D - ω - ω D - ω - ω D ω ω D ξξ ξξ ξξ K a K a K a - K a - FIG. 4: (Color online) Schematic description of terms definedby Eqs. (49)–(51).
Debye frequency ω D by Z c ( ξ )= Z c , Z a ( ξ ) = − Z a , ( ξ ≤ − ω ′ )0 , ( | ξ | < ω ′ ) Z a , ( ξ ≥ ω ′ ) (49) N c ( ξ )= N c , N a ( ξ ) = − N a , ( ξ ≤ − ω ′ )0 , ( | ξ | < ω ′ ) N a , ( ξ ≥ ω ′ ) (50) K c ( ξ, ξ ′ )= K c , K a ( ξ, ξ ′ )= − K a , (cid:2) ( ξ ≤ − ω ′ and ξ ′ ≤ − ω ′ )or ( ξ ≥ ω ′ and ξ ′ ≥ ω ′ ) (cid:3) K a , (cid:2) ( ξ ≤ − ω ′ and ξ ′ ≥ ω ′ )or ( ξ ≥ ω ′ and ξ ′ ≤ − ω ′ ) (cid:3) , (cid:0) | ξ | < ω ′ or | ξ ′ | < ω ′ (cid:1) . (51)Here, in addition to ω D , we have also introduced ω ′ whichspecifies the energy scale where DOS and kernels sub-stantially deviate from the values at E F . These formsare depicted in Fig. 4.Similarly to McMillan’s procedure, we assume the fol-lowing rectangular form for the solution∆( ξ ) = ∆ − , ( − ω ′ ≥ ξ ≥ − ω D )∆ , ( | ξ | < ω ′ )∆ + , ( ω D ≥ ξ ≥ ω ′ ) . (52)The condition that a nonzero solution for ∆ exists gives T c . By retaining the lowest order term with respect tothe values with subscript “a”, we obtain T c ∝ exp (cid:26) Z c K c N c − N a Z a ln[ ω D /ω ′ ] N c (1 + Z c ) + Z ln[ ω D /ω ′ ](1 + Z c ) (cid:27) . (53)Since Z a is dependent on N a through Eq. (40), we cansubstitute an approximate form for Z a in terms of N a /N c ,and then we obtain T c ∝ exp ( − λλ − (cid:18) N a N c (cid:19) Λ( r, λ )ln r ) . (54)Here, r ≡ ω D /ω ′ , and we have used the following relation K c N c = − λ and Z c = λ (Ref. 9). The first term in thebracket corresponds to the zeroth order, which appearsin literatures. The function Λ( r, λ ) is a positive functionwhich monotonically increases by increasing r or λ andconverges to a finite value < r → ∞ and λ → ∞ . The detail of the derivation of Eqs. (53) and(54) is given in Appendix B.From Eq. (54), we see how the antisymmetric part af-fects T c : It reduces T c regardless of the sign of N a andits amount becomes substantial when the ratios N a /N c or r ≡ ω D /ω ′ are large. In addition, we find the anti-symmetric part of the nondiagonal kernel ( K a ) does notcontribute to T c within the lowest order. B. Numerical calculation
The above analysis tells us that the asymmetry effectbecomes pronounced for the cases where N a /N c and r are large. To quantify the effect more explicitly, we cal-culated T c by numerically solving Eq. (48) using a modelDOS. We included both the phonon and electron con-tribution to the kernels, K = K ph + K el and Z = Z ph .We employed the n k -averaged form for K ph [Eq. (23) inRef. 14], and we treated K el by an approximate form K el ( ξ, ξ ′ ) = µN (0) (55)with a certain energy cutoff. For Z ph , we employed theprevious form [Eq. (45)] and the new one [Eq. (40)].We introduced a step-like DOS as depicted in the up-per panel of Fig. 5 so that N a /N c and r can be large.This form is characterized by two parameters; the ratioof the high-energy and low-energy values N + /N − , andthe energy range where DOS changes d . For the Eliash-berg function, we used the following model function α F ( ω ) = λ σ √ π ω exp " − (cid:18) ω − ω D σ (cid:19) (56)with width σ = ω D /
10. We set ω D = 400[cm − ], d = ω D /
2, and N + /N − = 6. This setting is realistic inthat one often see a similar dependence in doped semi-conductors with quasi-2D Fermi surfaces. With these settings, we calculated T c for various set-tings of E F . The calculation was performed with a log-arithmically uniform energy grid, where the minimumenergy cutoff was set to | ξ | = 10 − eV and the numberof points per digit was 20. The energy cutoff for K el wasfixed to 20eV. The parameters λ and µ were fixed to 1 . . E F so that the resulting T c becomes & T c using Z ph , new ( ξ ) [Eq. (40)] to that using Z ph ( ξ ) [Eq. (45)] as a function of E F , where the absolutevalues of T c are given in the inset. The ratio is system-atically less than unity; the effect of asymmetry lowers T c for any E F , which is consistent with Eq. (54). As theFermi level approaches to the lower edge of the slope,the decrease of the ratio becomes significant and finallyamounts to more than 20%. When the Fermi level is lo-cated on the upper edge of the slope, the decrease is notappreciable. This E F dependence can also be consistentlyunderstood with Eq. (54); the characteristic energy scale( ω ′ ) decreases as the Fermi level approaches to the slope,and the ratio N a /N c is relatively small when the Fermilevel is on the upper edge. Thus, the change of T c by thepresent improvement becomes crucial when E F is closeto a point where DOS increases rapidly.We show the calculated value of Z ph for E F = − . ω D in Fig. 6 to examine the properties of the present form.The model DOS is also depicted in the same energy scale(green shaded area). For both T =0 .
01K (thick lines) and T ≃ T asymc (thin lines), Z ph , new ( ξ ) becomes asymmetric in ξ , but they are similar to Z ph ( ξ ). In particular, the lim-iting value for ξ → ξ → Z ph ( ξ ). Thesefeatures well represent the properties (i)–(iii) describedin the last paragraph of Sec. II B. We observe a differ-ence between Z ph , new ( ξ ) and Z ph ( ξ ) around | ξ | & ω D ,which is responsible for the difference of T asymc and T symc in Fig. 5. Clearly, Z ph , new ( ξ ) becomes larger for the en-ergy region with larger density of states, whereas Z ph ( ξ )does not. The dependence of the former is reasonablebecause larger DOS should result in stronger mass renor-malization due to stronger total electron-phonon cou-pling. Thus we were able to see that the particle-holeasymmetric electronic structure is properly treated with Z ph , new in a numerically stable manner.Finally, let us comment on the application to actual su- E [ ω ] F D E [ ω ] F D c [ K ] T asym c T sym c a s y m c s y m c / TT T DO S N - N + d FIG. 5: (Color online) (Bottom) Ratio of T c calculated us-ing Z ph , new [Eq. (40)] ( T asymc ) to that calculated using Z ph [Eq. (45)] ( T symc ) plotted as a function of E F . In the inset,absolute values of T asymc and T symc are shown. (Top) DOSused for the calculations. We set N ( ξ )= N ( E F ) for | ξ | > [ ω D ] 0 1 2 3 -4 -2 0 2 4 T c ( 0.5) T c ( 0.5) Z ph , new ( ξ ) Z ph ( ξ ) Z ×× ξ asymasym FIG. 6: (Color online) Calculated values of Z ph , new ( ξ ) (redsolid line) and Z ph ( ξ ) (blue dotted line) for E F = − . ω D . Thegreen shaded area depicts DOS (arb. unit.) in the sameenergy scale. Note that the values for T = T asymc (thin lines)are multiplied by 0.5. perconductors. For typical superconductors, aluminumand niobium, we have carried out stable calculations of T c (see appendix C), confirming again that the listed prop-erties (i)–(iii) of Z ph , new are valid. The asymmetry ef-fect on T c was, however, estimated to be less than 0.1 %for the both systems. This is mainly due to the small N a /N c in these systems. The present weak asymmetryeffect gives support to the validity of the particle-holesymmetrizing treatment for many cases of con-ventional superconductors. IV. SUMMARY AND CONCLUSIONS
We constructed a new exchange-correlation kernel Z ph , new in SCDFT to treat the particle-hole asymme-try of electronic structure. The obtained n k -averagedforms [Eq. (40)] and n k -resolved form [Eq. (47)] do notshow any divergences or instabilities, and, for the systemswith good symmetry, well agree with the previous sym-metrized forms. By analytically deriving T c formula, wefound that the asymmetry systematically decreases T c .The amount of the decrease becomes substantial whenthe following two ratios are large: The ratio of the anti-symmetic component of DOS to the constant component,and that of the Debye frequency to the characteristic en-ergy scale of the DOS variation. We also calculated T c with a model step-like DOS, showing that the amountof the reduction of T c can be more than 20%. With thepresent work, we successfully reinforced the theoreticalfoundation of SCDFT to extend its applicability to sys-tems with significant particle-hole asymmetric electronicstructure. Acknowledgments
The authors thank Kazuma Nakamura, Shiro Sakai,Hideyuki Miyahara, and Jianting Ye for enlighteningcomments. This work was supported by Grants-in-Aidfor Scientic Research from JSPS (No. 23340095), Fund-ing Program for World-Leading Innovative R & D on Sci-ence and Technology (FIRST program) on Quantum Sci-ence on Strong Correlation”, JST-PRESTO and the NextGeneration Super Computing Project and NanoscienceProgram from MEXT, Japan.
Appendix A: Evaluation of Z ph2 ( ξ ) for nonconstantDOS We here describe the detail of evaluation of Z ph2 ( ξ )[Eq. (33)] for the nonconstant DOS. Under the de-composition approximation of DOS introduced inSec. II A [Fig. 2, Eq. (34)], Z ph , ( ξ ) is given by Z ph , ( ξ )= − (cid:26) ξ − β/ β/ ξ ]sinh[( β/ ξ ] (cid:27) × Z dωα F ( ω ) Z L − L dξ (cid:26) X i [ N − i θ − i ( ξ )+ N + i θ + i ( ξ )] (cid:27) × Z L − L dξ ′ (cid:26) X j [ N − j θ − j ( ξ ′ )+ N + j θ + j ( ξ ′ )] (cid:27) I ′ ( ξ, ξ ′ , ω ) . (A1)With omitting the temperature-dependent terms of or-der O [( T /ǫ i ) ], O ( e − βǫ i ), O ( e − βL ) and O ( e − βL ), theintegral for each term is carried out as Z L − L dξ Z L − L dξ ′ θ ± i ( ξ ) θ ± j ( ξ ′ ) I ′ ( ξ, ξ ′ , ω ) ≃ ± ln (cid:12)(cid:12)(cid:12)(cid:12) ( L + ǫ i + ω )( L + ǫ j + ω )( L + ǫ j + ω )( L + ǫ i + ω ) (cid:12)(cid:12)(cid:12)(cid:12) , (A2) Z L − L dξ Z L − L dξ ′ θ ± i ( ξ ) θ ∓ j ( ξ ′ ) I ′ ( ξ, ξ ′ , ω ) , ≃ ± ln (cid:12)(cid:12)(cid:12)(cid:12) ( L + ǫ i + ω )( L + ǫ j + ω )( L + ǫ i + ω )( L + ǫ j + ω )( L + L + ω ) ( ǫ i + ǫ j + ω ) (cid:12)(cid:12)(cid:12)(cid:12) , (A3) Z L − L dξ Z L − L dξ ′ θ − i ( ξ ) I ′ ( ξ, ξ ′ , ω ) , ≃ ln (cid:12)(cid:12)(cid:12)(cid:12) ( L + L + ω ) ( ǫ i + ω ) ( L + ω )( L + ǫ i + ω )( L + ω )( L + ǫ i + ω ) (cid:12)(cid:12)(cid:12)(cid:12) , (A4) Z L − L dξ Z L − L dξ ′ θ − i ( ξ ′ ) I ′ ( ξ, ξ ′ , ω ) , ≃ ln (cid:12)(cid:12)(cid:12)(cid:12) ( L + ǫ i + ω )( L + ǫ i + ω )( L + L + ω ) (cid:12)(cid:12)(cid:12)(cid:12) , (A5) Z L − L dξ Z L − L dξ ′ θ + i ( ξ ) I ′ ( ξ, ξ ′ , ω ) , ≃ ln (cid:12)(cid:12)(cid:12)(cid:12) ( L + ǫ i + ω )( L + ω )( L + ǫ i + ω )( L + ω ) (cid:12)(cid:12)(cid:12)(cid:12) , (A6) Z L − L dξ Z L − L dξ ′ θ + i ( ξ ′ ) I ′ ( ξ, ξ ′ , ω ) , ≃ ln (cid:12)(cid:12)(cid:12)(cid:12) L + ǫ i + ωL + ǫ i + ω (cid:12)(cid:12)(cid:12)(cid:12) . (A7) Using these relations, Z ph2 ( ξ ) is transformed to the sameform as Eq. (35) with the opposite sign. Appendix B: Derivation of Eq. (53) and Eq. (54)
Let us first plug the kernels and the gap function givenin Eqs. (49)–(52) into the energy-averaged gap equation[Eq. (48)]. The condition to have a nonzero solution isdet M =0 ,M ≡ Z c − Z a + ( N c − N a )( K c − K a ) q N c K c p ( N c + N a )( K c + K a ) q ( N c − N a ) K c q Z c + N c K c p ( N c + N a ) K c q ( N c − N a )( K c + K a ) q N c K c p Z c + Z a + ( N c + N a )( K c − K a ) q , (B1)with p and q defined by p = Z ω ′ dξ ′ tanh[( β c / ξ ′ ] ξ ′ , q = Z ω D ω ′ dξ ′ tanh[( β c / ξ ′ ] ξ ′ . (B2)Here β c denotes the inverse of T c . In order to treat thelowest-order contribution, we keep the terms up to thesecond order with respect to the values with subscript“a” asdet M =0 ≃ (1+ Z c ) + [ K c N c C − K a N c q ](1+ Z c ) +[ − K c K a N qC − K c N a Z a q − Z ](1 + Z c ) − K c N c Z p, (B3) where C ≡ p + q has been introduced. Retaining thelowest-order terms, we obtain C ≃− Z c K c N c + N a Z a qN c (1 + Z c ) − Z q (1 + Z c ) . (B4)Using C = ln [2 e γ ω D / ( πT c )] with γ being the Euler con-stant, we get Eq. (53).Next we consider the relation between Z a and N a . Letus start from our newly developed form [Eq. (40)]. Theantisymmetric part of the integrand in ξ ′ can be written0as˜ I ( ξ, ξ ′ , ω ) − ˜ I ( ξ, ξ ′ , − ω ) − ˜ I ( ξ, − ξ ′ , ω ) + ˜ I ( ξ, − ξ ′ , − ω ) − ˜ J ( ξ, ξ ′ , ω )+ ˜ J ( ξ, ξ ′ , − ω )+ ˜ J ( ξ, − ξ ′ , ω ) − ˜ J ( ξ, − ξ ′ , − ω ) . (B5)Using the facts ˜ I ( − ξ, − ξ ′ , − ω ) = ˜ I ( ξ, ξ ′ , ω ) and˜ J ( − ξ, − ξ ′ , − ω ) = ˜ J ( ξ, ξ ′ , ω ), one can easily find theabove part divided by tanh[( β/ ξ ] is antisymmetric in ξ . The symmetric part [Eq. (46)] divided by tanh[( β/ ξ ]is, on the other hand, symmetric in ξ . Consequently, theantisymmetric part of Z ph , new ( ξ ) is yielded by only theantisymmetric part of N ( ξ ′ ).Substituting N ( ξ ′ ) = N c ( ξ ′ ) + N a ( ξ ′ ), for the antisym-metric contribution, we obtain Z dξ ′ N a ( ξ ′ ) N [ I ( ξ, ξ ′ , ω ) − J ( ξ, ξ ′ , ω )]= 2 N a N c Z ω D ω ′ dξ ′ [(antisym.)] , (B6)where (antisym.) represents the terms in Eq. (B5). Thisintegration is easily performed since we can assume ξ ≫ T , ξ ′ ≫ T , and ω ≫ T , so that we get2 N a N c Z ω D ω ′ dξ ′ [(antisym.)] ≃ N a N c Z ω D ω ′ dξ ′ ξ + ξ ′ + ω ) = 2 N a N c (cid:20) ξ + ω ′ + ω − ξ + ω D + ω (cid:21) . (B7)Finally, by using the Einstein spectrum α F ( ω ) = λ ω D δ ( ω − ω D ), we obtain12 (cid:2) Z ph , new ( ξ ) −Z ph , new ( − ξ ) (cid:3) = sgn( ξ ) λ N a N c ω D ( ω D − ω ′ )( ξ + ω ′ + ω D )( ξ +2 ω D ) , (B8)where the sign function comes from lim β →∞ / tanh[( β/ ξ ].Subsequently, we consider to approximate this formas the rectangular form [Eq. (49)]. Using the value at ξ = ω D yields Z a = λ N a N c r − r + 1) ≡ λ N a N c h ( r ) , (B9)where r = ω D /ω ′ >
1. By using this form in Eq. (53) andsubstituting K a N c = − λ , Z c = λ , we get to T c ∝ exp ( − λλ − (cid:18) N a N c (cid:19) " λh ( r )1+ λ − (cid:18) λh ( r )1+ λ (cid:19) ln r ) . (B10) This expression is equivalent to Eq. (54), where the func-tion Λ( r, λ ) is defined as the part in the square bracketin Eq. (B10). Its positiveness, monotonic dependencesand the convergence in the limit r → ∞ and λ → ∞ referred in Sec. III A are easily confirmed with this ex-pression. We also note that Λ( r, λ ) has these propertiesfor any values of ξ ( ω D >ξ>ω ′ ) in Eq. (B8) when we getan approximate form of Z a . Appendix C:
Ab initio calculation
We applied our formalism to typical weak- and strong-coupling superconductors, Al and Nb. We solved the gapequation [Eq. (1)] with the energy-averaged approxima-tion for the phonon kernels K n k ,n ′ k ′ = K ph ( ξ n k , ξ n ′ k ′ ) + K el n k ,n ′ k ′ , (C1) Z n k = Z ( ξ n k ) , (C2)where K ph ( ξ n k , ξ n ′ k ′ ) and K el n k ,n ′ k ′ are defined byEq. (23) in Ref. 14 and Eq. (13) in Ref. 38, respectively.For the diagonal part Z ( ξ n k ), we employed Z ph , new ( ξ n k )[Eq. (40)] and Z ph ( ξ n k ) [Eq. (45)]. We used an accuraterandom-sampling scheme given in Ref. 35, and the result-ing sampling error was not more than a few percent. Thedetailed condition of the calculations is given in Ref. 39.In Fig. 7, the calculated DOS for Al [(a)] and Nb[(b)] are shown. The calculated values of Z n k with the [eV] [eV] DO S -10 -5 0 5 10Energy[eV] DO S (a) (b)(c) (d) Z n k Z n k α Fα F ξ n k ξ n k FIG. 7: (Color online) Calculated DOS for Al [(a)] and Nb[(b)]. Z kernels for Al [(c)] and Nb [(d)] calculated fromEq. (40) (solid line) and Eq. (45) (dashed line): The in-put Eliashberg function (dotted line) and DOS (green shadedarea) are given in the same energy scale, with which λ and ω ln were estimated for Al (Nb) as 0.417 (1.267) and 314 K(176 K). α F . For the bothsystems, the DOS (green shaded area) has small butnonzero antisymmetric component in the correspondingenergy scale. However, the two forms of Z yield approx-imately the same value. For the calculated T c , regard-less of the form of Z n k , we obtained 0.7 and 8.7 K for Al and Nb, respectively. Although the gap values calcu-lated with Z ph , new ( ξ n k ) were approximately 0.1% smallerthan those calculated with Z ph ( ξ n k ) for each set of thesampling points, this difference was within the samplingerror. The present result indicates that the asymmetryeffect is not significant for T c in Al and Nb. W. L. McMillan, Phys. Rev. , 331 (1968). P. B. Allen and R. C. Dynes, Phys. Rev. B , 905 (1975). P. Morel and P. W. Anderson, Phys. Rev. , 1263(1962); N. N. Bogoliubov, V. V. Tolmachev, and D. V.Shirkov,
A New Method in the Theory of Superconductiv-ity (1958) (translated from Russian: Consultants Bureau,Inc., New York, 1959). S. Baroni, S. deGironcoli, A. Dal Corso, and P. Giannozzi,Rev. Mod. Phys. , 515 (2001). K. Kunc and R. M. Martin, in
Ab initio Calculation ofPhonon Spectra , edited by J. T. Devreese, V. E. van Doren,and P. E. van Camp (Plenum, New York, 1983), p. 65. K. H. Lee, K. J. Chang, and M. L. Cohen, Phys. Rev. B , 1425 (1995). L. N. Oliveira, E. K. U. Gross, and W. Kohn, Phys. Rev.Lett. , 2430 (1988). T. Kreibich and E. K. U. Gross, Phys. Rev. Lett. , 2984(2001). M. L¨uders, M. A. L. Marques, N. N. Lathiotakis, A. Floris,G. Profeta, L. Fast, A. Continenza, S. Massidda, and E.K. U. Gross, Phys. Rev. B , 024545 (2005). A. B. Migdal, Sov. Phys. JETP , 996 (1958). G. M. Eliashberg, Sov. Phys. JETP , 696 (1960). D. J. Scalapino, in
Superconductivity edited by R. D. Parks,(Marcel Dekker, New York, 1969) VOLUME 1. J. R. Schrieffer,
Theory of superconductivity; Revised Print-ing , (Westview Press, Colorado, 1971) M. A. L. Marques, M. L¨uders, N. N. Lathiotakis, G. Pro-feta, A. Floris, L. Fast, A. Continenza, E. K. U. Gross, andS. Massidda, Phys. Rev. B , 024546 (2005). A. Floris, G. Profeta, N. N. Lathiotakis, M. L¨uders, M. A.L. Marques, C. Franchini, E. K. U. Gross, A. Continenza,and S. Massidda, Phys. Rev. Lett. , 037004 (2005). A. Sanna, G. Profeta, A. Floris, A. Marini, E. K. U. Gross,and S. Massidda, Phys. Rev. B , 020511(R) (2007). C. Bersier, A. Floris, A. Sanna, G. Profeta, A. Continenza,E. K. U. Gross, and S. Massidda, Phys. Rev. B , 104503(2009). L. R. Testardi, Rev. Mod. Phys. , 637 (1975) K. M. Ho, M. L. Cohen, and W. E. Pickett, Phys. Rev.Lett. , 815 (1978). W. E. Pickett, Phys. Rev. B , 3897 (1980); , 1186(1982); W. E. Pickett and B. M. Klein, Solid State Com-mun. , 95 (1981). S. G. Lie and J. P. Carbotte, Solid State Commun. , 511(1978); , 599 (1980); , 127 (1980). B. Mitrovi`c and J. P. Carbotte, Can. J. Phys. , 784(1983). E. Schachinger, M. G. Greeson, and J. P. Carbotte, Phys.Rev. B , 406 (1990). Y. Yokoya and Y. O. Nakamura, Solid State Commun. ,133 (1996); Physica C , 187 (1996). Y. Yokoya, Phys. Rev. B , 6107 (1997). A. D. Caviglia, S. Gariglio, N. Reyren, D. Jaccard, T.Schneider, M. Gabay, S. Thiel, G. Hammerl, J. Mannhart,and J. M. Triscone, Nature , 624 (2008). K. Ueno, S. Nakamura, H. Shimotani, A. Ohtomo, N.Kimura, T. Nojima, H. Aoki, Y. Iwasa, and M. Kawasaki,Nat. Mater. , 855 (2008) J. T. Ye, S. Inoue, K. Kobayashi, Y. Kasahara, H. T. Yuan,H. Shimotani, and Y. Iwasa, Nat. Mater. , 125 (2010). K. Ueno, S. Nakamura, H. Shimotani, H. T. Yuan, N.Kimura, T. Nojima, H. Aoki, Y. Iwasa, and M. Kawasaki,Nat. Nanotechnol. , 408 (2011). K. Taniguchi, A. Matsumoto, H. Shimotani, and H. Takagi,Appl. Phys. Lett. , 042603 (2012). J. T. Ye, Y. J. Zhang, R. Akashi, M. S. Bahramy, R. Arita,and Y. Iwasa, Science , 1193 (2012). P. B. Allen and B. Mitrovi´c, in
Solid State Physics , editedby H. Ehrenreich, F. Seitz, and D. Turnbull (Academic,New York, 1982), Vol. 37, p. 1. We note that the condition ǫ i ≫ T must be retained whenwe take limit N ǫ → ∞ . We found that it is necessary to multiply p ( ξ ′ ± ω ) alsofor the terms originating from J ; otherwise, numerically I − J does not become zero at ξ =0 and consequently thekernel diverges for ξ → R. Akashi, K. Nakamura, R. Arita, and M. Imada, Phys.Rev. B , 054513 (2012). C. Felser and R. Seshadri, J. Mater. Chem. , 459 (1999). The non-monotonic dependence of T asymc and T symc on E F is due to our approximation that λ and µ are constantregardless of E F , and so we do not focus on that. S. Massidda, F. Bernardini, C. Bersier, A. Continenza, P.Cudazzo, A. Floris, H. Glawe, M. Monni, S. Pittalis, G.Profeta, A. Sanna, S. Sharma, and E. K. U. Gross, Super-cond. Sci. Technol. , 034006 (2009). To calculate exchange-correlation kernels and electronicenergy eigenvalues, we performed electronic structure andphonon calculations within the local-density approxima-tion using ab initio plane-wave pseudopotential cal-culation codes
Quantum Espresso . The pseudopo-tential for Al was generated in the configuration of(3 s ) . (3 p ) . (3 d ) . , whereas that for Nb was generated inthe configuration of (4 s ) . (4 p ) . (4 d ) . . The plane-waveenergy cutoff was set to 20 and 80 Ry for Al and Nb,respectively. The charge density was calculated with the16 × × k points in the Monkhorst-Pack grid for Aland Nb. The optimized lattice parameters used for thecalculations were 7.483 and 6.099 bohr for Al (fcc) andNb (bcc). The dynamical matrices for Al and Nb werecalculated on the 8 × × q points, where the gaussian ofwidth 0.200 Ry was used for the Fermi-surface integra-tion. The phonon linewidths for Al (Nb) were integrated using supplementary 32 × ×
32 (48 × × k points andthe gaussian of width 0.020 Ry (the first-order Hermite-gaussian function of width 0.020 Ry). The Eliashbergfunctions were calculated using the tetrahedron interpola-tion method . Using the output of pw code, the electrondielectric functions ε were calculated within the random-phase approximation on the 7 × × × × q pointsusing the tetrahedron interpolation method with the Rath-Freeman treatment considering 30 (40) bands for Al(Nb). The SCDFT gap equation was solved for Al (Nb)with 15000 (18000) k points for bands crossing the Fermilevel and 500 (600) points for the other bands, where weconsidered 12 (26) bands. D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. , 566(1980). J. P. Perdew and A. Zunger, Phys. Rev. B , 5048 (1981). P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car,C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococ- cioni, I. Dabo, A. Dal Corso, S. Fabris, G. Fratesi, S.de Gironcoli, R. Gebauer, U. Gerstmann, C. Gougous-sis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari,F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L.Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P.Seitsonen, A. Smogunov, P. Umari, and R. M. Wentz-covitch, J. Phys.: Condens. Matter N. Troullier and J. L. Martins, Phys. Rev. B , 1993(1991). M. Methfessel and A. T. Paxton, Phys. Rev. B , 3616(1989). G. Lehmann and M. Taut, Phys. Stat. Sol. , 469 (1972). M. S. Hybertsen and S. G. Louie, Phys. Rev. B , 5585(1987); , 5602 (1987). J. Rath and A. J. Freeman, Phys. Rev. B11