Direct evaluation of the isotope effect within the framework of density functional theory for superconductors

Within recent developments of density functional theory, its numerical implementation and of the superconducting density functional theory is nowadays possible to predict the superconducting critical temperature, , with sufficient accuracy to anticipate the experimental verification. In this paper we present an analytical derivation of the isotope coefficient within the superconducting density functional theory. We calculate the partial derivative of with respect to atomic masses. We verified the final expression by means of numerical calculations of isotope coefficient in monatomic superconductors (Pb) as well as polyatomic superconductors (CaC6). The results confirm the validity of the analytical derivation with respect to the finite difference methods, with considerable improvement in terms of computational time and calculation accuracy. Once the critical temperature is calculated (at the reference mass(es)), various isotope exponents can be simply obtained in the same run. In addition, we provide the expression of interesting quantities like partial derivatives of the deformation potential, phonon frequencies and eigenvectors with respect to atomic masses, which can be useful for other derivations and applications.


Introduction
The discovery of isotope effect (IE) in superconductors, that is the variation of the superconducting critical temperature (T c ) upon isotope substitution of the constituent atoms of the material, has had a fundamental historical role. In fact, after coefficient, equal to 1 2 . Although, this is the results of the many approximations within BCS theory, the prediction was quali tatively confirmed for many superconductors (Zn, Cd, Sn, Hg, Pb, Tl) [6]. However, the observation of reduced isotope effect (γ < 1 2 ) in many materials (Ru, Os, Mo) [6] posed more strin gent limits in a microscopic theory of superconductivity.
Beyond its historical importance, the measure of isotope effect also represents a key experimental tool in the study of superconductivity. Generally, IE is considered a measure of the phonon contribution on the electron pairing and its lack, as a possible indication of other pairing mechanisms. For example, in recent years, IE played a fundamental role con firming phonon mediated superconductivity in SH 3 [7], MgB 2 [8] and in CaC 6 [9], and several studies of isotope shifts on T c have been carried in almost all known cuprates [10] and iron superconductors [11].
Strictly speaking the value of the isotope coefficient has only a relative importance, without a reliable theory that interprets it. In fact, isotope coefficient is a 'number' that in principle contains all the relevant physics of the normal (non super conducting) phase and its final value can hinder peculiar and interesting effects like anharmonic contrib ution to explain the inverse isotope effect in superconducting Palladium-Hydride compounds [12]. On the contrary, even without anharmonic contributions, its final value accounts for all the interactions responsible for the superconducting phase: electron-phonon, Coulomb interactions, spin fluctuations, etc. A theory able to predict the isotope coefficient must be able to deal with all these relevant interactions.
Due to its fundamental importance, the straightforward comparison with the experiments and to complete the form ulation of the SCDFT, in this paper we present an analytical derivation of the isotope coefficient within the SCDFT. This new approach will allow to obtain isotope coefficient as a post processing with a considerable gain in terms of computational time and precision. In addition, the analytical calculation of partial derivatives of the interaction functionals gives access to many interesting quantities like derivative of deformation potential, phonon frequencies and eigenvectors with respect to atomic masses.
In the following we present the analytical derivation of IE (section 2) then, we performed the partial derivatives (section 3) and then present a numerical implementation and discuss results on Pb and CaC 6 (section 4)

Basic definitions
Considering a superconductor with critical temperature T c , the isotope effect is described by the dependence of T c on the mass of the constituent atoms. In the case of oneatom per unit cell: where M is the atomic mass. In case of several atoms per unit cell, this can be generalized to or making use of a reference mass, M 0 , with reference critical temperature, T 0 c : For small changes in the mass we obtain The partial isotopecoefficient, γ α , can be obtained from or, if we introduce the critical inverse temperature β c = 1/(k B T c ),

Dependence of the critical temperature on the masses
Within the SCDFT [14,15] the superconducting critical temperature is determined by the linearized gap equa tion (lin. gap. equation) (refer to [14,15] for the meaning of the symbols) 8 : The interaction kernel, K kk , contains the sum of two contrib utions: the electron-phonon (el-ph) interaction (K ph kk ) and the Coulomb electron-electron (el-el) interaction (K el kk ). First, we rewrite the linearized gap equation as a matrix equation with a symmetric kernel: The linearized gap equation is only valid at the critical temper ature, or equivalently β c , which enters the equation via the kernel M, which also depends-through the el-ph matrix ele ments [14,15]-on the nuclear masses.
In order to calculate the dependence of β c on the nuclear masses, we first rewrite the lin. gap equation as an eigenvalue equation, so that it is valid for all temperatures.
The critical temperature is then defined by the lowest eigen value being equal to unity Since the eigenvalue equation (12) depends separately on β and on the spectral function, also the eigenvalue λ 9 depends (independently) both on β and the masses M α , while β c depends, on the masses M α through equation (13).
We now take the derivative of equation (13) with respect to the masses: From this we obtain: Next, we need to calculate the partial derivatives of the eigen value λ. To this end we take the partial derivatives of equa tion (12) with respect to β and M α . Using the normalization of ∆ (for each β), we get Again, due to the normalization, we find This leads us to the desired result: that, from the definition in equation (7) is proportional to the partial isotope coefficient.

Partial derivatives of the kernels
In order to calculate equation (19), we need to evaluate the partial derivatives of the interaction kernels with respect to the atomic masses ( M α ) and inverse temperature (β), which represent the main part of the present paper. From equation (11) we find (21) where we used that the electronic contribution (K el ) does not depend on the nuclear masses.
So, for the numerator and denominator of equation (19) we obtain:

Analytical derivatives
In order to calculate the above terms, we need to perform par tial derivatives of the renormalization term ( Z k ) and of inter action kernel (K ph kk ) with respect to both nuclear masses ( M α ) and inverse temperature (β). In this work we adopt the kernels (density functionals) of [14] 10 : 9 Not to be confused with the electron-phonon parameter. 10 In [14] the Coulomb interaction K el kk was considered in the Thomas-Fermi approximation. However, the formalism is equally valid for any other type of screening as the RPA [20,27,82], therefore with respect to [14] we drop here the TF label of K el .
The partial derivative with respect to β does not pose rele vant problems, because the temperature is present only within the Fermi and Bose factors, which have a clear analytical expression. The complete expression of the partial derivatives of Z k and K ph kk with respect to β is reported in appendix A.
3.1.1. The kernel derivatives. The partial derivatives with respect to the masses of the kernels are more involved. They can be cast as: The above expressions involve the calculation of the deriv ative of the both phonon frequencies and electron-phonon coupling constant with respect to the atomic masses, which we afford in the next section.

Derivatives of the phonon frequencies.
For a mode λ at wavevector q the polarization vectors ζ λq αµ (with α the atom index and µ the cartesian coordinate) and phononfrequencies Ω λq are determined as eigenvectors and eigenvalues of the dynamical matrix [44]: where each matrix element does not depend on the nuclear masses (the sum is over the atoms of the crystal at cell position R i , with respect the refer ence cell R 0 ). From equation (30), and using, again, the (massinde pendent) normalization of the eigenvectors we obtain (see appendix B for the details of the calculations): 3.1.3. Derivatives of the coupling constants. The electronphonon coupling constants are given by the matrix element of the derivative of the selfconsistent potential (V λq (r)) with respect to the phonon mode λ at wavevector q between elec tron wavefunctions at wavevectors k and k = k + q: The selfconsistent potential can be written as [45]: where χ(r, r ) is the full response function, and the derivative of the bare potential. Since the full response function is not explicitly known, the response equation (34) is normally rewritten as which involves only the Kohn-Sham response function (χ 0 (r , r )), which is explicitly known, but it has to be solved selfconsistently with respect to V λq (r).
In order to calculate the partial derivative of the coupling constant with respect to the masses we need to calculate the partial derivative of the potential, We write V 0 λq (r) as: The same change of basis can be applied to V λq (r): Thus equation (34) can be written as: This form of the displacement potential now allows us to take the derivatives with respect to the nuclear masses without additional solutions of the response equation, since the response screening is performed on potentials which do not depend on the masses, and the eigenvectors, frequencies, etc are pulled outside the response function. Therefore, At this point we only need the derivative of the phonon fre quencies and eigenvectors with respect to the masses, which are explicitly calculated in appendices B and C. At this point, we transform also the matrix element through withg αµ;q kk = d 3 rϕ * k (r)Ṽ αµ;q (r)ϕ k (r).
Using equations (42) in (44) and including the derivatives obtained in appendices B and C, we obtain: 11 As general rule we decided to choose any temperature in which ∆ 2 k becomes linear with T. The BCStype second order phase transition requires that near T c the superconducting gap may be expressed as ∆ ∼ (1 − T Tc ) 1 2 . Anyway, we always verified that this last choice ensures that γ does not depend on the temperature.
Since the eigenvectors of the dynamical matrix form a unitary matrix, the g αµ;q kk can also be generated from the firstprinci ples calculated g λq kk through the transformatioñ Finally, the required derivative of the modulus squared is (48) and this completes our derivation of equation (7). The complete analytic expression of the isotope coefficients is clearly very complicated and we are not able to extract exact features from it, as an upper limit for the isotope coeffi cients. However it can be easily evaluated numerically for any given material.

Implementation and test cases: Pb and CaC 6
We test the validity of the analytical approach by direct tests with usual finite difference method (FDM). In this last case, in order to calculate γ α we need to calculate the critical temper ature for two different isotopic masses of the same αatom. The two approaches must give the same result. Evaluation of γ α (equation (7)), requires the calculation of the oneelectron energies near the Fermi level (ξ k ), phonon frequencies (Ω λq ), the phonon eigenvectors (ζ λq αµ ) and elec tron-phonon matrix elements (g λq kk ). All these quantities have been calculated within density functional theory perturbation theory [44,[46][47][48] using a pseudopotential implementation [49]. In addition, we solved the gap equation to find T c and obtain the values of ∆ k sufficiently close to T c (where the gap vanishes by definition) 11 . Thus, the main steps of the calcul ation can be summarized as: (i) calculate the el-ph matrix ele ments and phonons frequencies and eigenvectors on a regular grid for the electron wavevector (k) and phonon wavevectors (q) in the Irriducible Brillouin zone. (ii) Obtain ∆ k on a non uniform [15] mesh of k and k' near the Fermi level and solve the superconducting gap equation to find the value of T c . (iii) Calculate γ α using the same nonuniform mesh.
Pb crystallizes in a fcc lattice and its Fermi surface is com posed of three sheets arising from the Pb p orbitals; at 7.2 K it becomes superconductor showing a multigap superconduc tivity [43].
The calculation was performed using normconserving pseudopotential [57] with 32 Ry of energy cutoff. Converged Brillouin zone integration of the charge density was achieved with 26 3 Monkhorst-Pack (MP) [58] grid. Phonon frequen cies and electron-phonon matrix elements are calculated on a 8 3 MP q grid and 26 3 MP k grid. The selfconsistent gap equation was solved on a nonuniform kmesh in the IBZ of 6000 (500) points for bands crossing (not crossing) the Fermi level. The results are identical to those reported in [43] and not discussed here. We focus, instead, on the isotope effect. First, we determine the isotope coefficient calculating the crit ical temperature T 1 c and T 2 c for two different isotopic masses, M 1 = 207.21 and M 2 = 208 respectively. Thus γ is obtained as: The error bar accounts for the indeterminacy in the deter mination of T c .
Application of the analytical procedure, using the same sets of k and q points, gives γ = 0.48. The experimental value is 0.48 ± 0.01 [59,60]. The agreement is excellent.
In principle, the only value that can be compared with the experimental one, is the one obtained in the same conditions, that is, with a finite variation of the isotopic mass. The analyt ical derivation, assumes an infinitesimal variation around the reference mass. However, the dependence of γ on the refer ence mass itself is weak and irrelevant for practical purposes.
CaC 6 represents a stringent test for our analytic derivation. In fact, the multiatom case is far from trivial, as it accounts for many different interaction between different masses. In order to be consistent with calculations of [65], we used the same computational details. Even in this case, our analyt ical procedure reproduces (within 0.01) the finitedifference values (obtained as described in the case of Pb). In particular, we find γ Ca = 0.24 and γ C = 0.26. These values in agreement with previous estimations [63,65] outline the important role played by both C and Ca atoms in the superconducting phase transition and the high coupling of low energy Ca phonons, in agreement with recent ARPES measurements [70]. However the Ca isotope coefficient disagrees with the values measured in [9], that estimates a value between 0.4 and 0.56. This point, remains up to now unresolved [81].

Conclusions
This paper presents a development of the superconducting density functional theory presented in [14,15]. We derived an analytical expression of the isotope coefficient in the framework of SCDFT which has many advantages with respect to the finite difference method. In particular: (i) in the finite difference method, a numerically stable value of γ, requires an evaluation of T c with a high enough acc uracy (less than 1%). This is a formidable task by itself. In fact, close to T c the selfconsistent gap is approaching zero, making the selfconsistency very difficult and time con suming to achieve. In addition, in order to obtain a reliable estimate of T c the calculation of the superconducting gap has to be performed for many different temperatures close to T c . In case of multiband and polyatomic systems this pro cedure can represent a serious limitation. The advantage of the analytical approach is that it does not require extra self consistent solutions of the gap equation. (ii) The determina tion of γ needs (at least) two independent calculations of T c for two different atomic masses. The new approach needs only one calculation of T c at the reference mass. (iii) In case of polyatomic systems, evaluation of γ α for each species is an independent calculation. In the analytical approach γ α is calculated within the same run for each atomic species (α). Therefore the overall calculation of isotope coefficients is reduced by a factor 2 × N scf × N s , where N scf is the number of selfconsistent cycles in the solution of the gap equa tion and N s the number of atomic species in the system.
We want to underline that the present analytical approach (even if at the moment rather complicated) is suitable to analyze the dependence of the isotope coefficients on many different normal states properties (electron eigenenergies, phonon frequencies and eigenvectors, superconducting gap, Coulomb potentials, etc.).
In conclusion, the present formulation of the isotope effect, completes the formal derivation of SCDFT. As for the BCS formulation, we were able to derive an analytical expression of isotope effect. This allows a relevant time saving, higher accuracy with respect to the usual numerical calculation and an analytical framework to analyze the com putational results.

Appendix A. Partial derivatives of interaction kernels with respect β
The partial derivative of Z k and K ph kk with respect to β are defined as: According to [14,15], the function I is defined as: and J as: In the following sections we will derive analytical expression of the derivative of I and J with respect to β.

Appendix B. Partial derivative of phonon frequencies with respect to atomic masses
From equation (30) we can write: which ultimately have to be inserted in equation (42).

Appendix C. Partial derivative of the phonon eigenvectors
We will use perturbation theory in order to calculate the partial derivative of the phonon eigenvectors. First, we linearize the dynamical matrix in terms of the masses: and then we will calculate the derivative as:  (42), we obtain (D.1)