Weak phonon-mediated pairing in BiS$_2$ superconductor from first principles

TL;DR

This study uses first-principles calculations to investigate the electron-phonon pairing mechanism in LaO$_{0.5}$F$_{0.5}$BiS$_2$, concluding that phonons alone cannot account for its high superconducting transition temperature.

Contribution

The paper provides a more accurate calculation of electron-phonon coupling in BiS$_2$ superconductors, challenging previous assumptions about phonon-mediated pairing.

Findings

Charge density wave instability confirmed

Electron-phonon coupling strength is lower than previous estimates

Phonon-mediated pairing cannot explain high T_c in this material

Abstract

Superconductivity in novel bismuth-sulphur superconductors has attracted large research efforts, both experimental and theoretical, but a consensus on the nature of superconductivity in these materials has yet to be reached. Using density functional theory for superconductors, we study the electron-phonon pairing mechanism in LaO$_{0.5}$F$_{0.5}$BiS$_2$. We first confirm the presence of a commensurate charge density wave instability, in accordance with previous studies. Using a recently developed integration scheme for the electron-phonon coupling, we found that its strength is much lower than previously calculated, due to improved density of state calculations. We finally conclude that conventional phonon-mediated pairing cannot explain the high superconducting transition temperatures observed in this material.

Tables

Table 1. Table 1: Superconducting transition temperatures and effective Coulomb repulsion factors in LaO 0.5 F 0.5 BiS 2 calculated using density functional theory for superconductors.

	static	dynamical
$T_c$	0.2 K	0.4 K
${\mu }^{\ast }$	0.18	0.16

Equations

$${\Delta }_{n\text{k}}=-{\mathcal{Z}}_{n\text{k}}{\Delta }_{n\text{k}}-\frac{1}{2}\sum _{n^{\prime }{\text{k}}^{\prime }}{\mathcal{K}}_{n\text{k}{n}^{\prime }{\text{k}}^{\prime }}\frac{\tanh [(\beta /2)E_{n^{\prime }{\text{k}}^{\prime }}]}{E_{n^{\prime }{\text{k}}^{\prime }}}{\Delta }_{n^{\prime }{\text{k}}^{\prime }}$$

$${\alpha }^{2}F(\omega )=$$

$$\times \delta ({ϵ}_{n^{\prime }\text{k}}-E_F)\delta (\omega -{\omega }_{\nu \text{q}})$$

$$\lambda =2\int d\omega \frac{{\alpha }^{2}F}{\omega }$$

$${\omega }_{ln}=\exp \left(\frac{\int d\omega \frac{{\alpha }^{2}F(\omega )}{\omega }\ln (\omega )}{\int d\omega \frac{{\alpha }^{2}F(\omega )}{\omega }}\right)$$

$$T_c=\frac{{\omega }_{ln}}{1.2}\exp \left(-\frac{1.04(1+\lambda )}{\lambda -{\mu }^{\ast }(1+0.62\lambda )}\right)$$

$${\alpha }^2{F}_{\text{resc}}(\omega )=$$

$$\times \frac{\tilde{N}({\tilde{E}}_F)}{N(E_F)}\widetilde{{\alpha }^{2}F}(\omega )$$

$${\alpha }^2{F}_{\text{resc}}(\omega )\equiv$$

$$\times \widetilde{{\alpha }^{2}F}(\omega )$$

Full text

I Introduction

Superconductivity was discovered in two layered materials sharing the same layer of bismuth and sulphur atoms Mizuguchi et al. (2012a, b); Awana et al. (2013), which has been shown to be the superconducting layer Usui et al. (2012); Usui and Kuroki (2015); Morice et al. (2015). The material having reached the highest superconducting transition temperature ($T_{c}$) of 10.6 K so far is La(O,F)BiS2 Mizuguchi et al. (2012b). Its parent compound is insulating and becomes superconducting upon doping, which raises the Fermi level in a rigid band fashion Usui et al. (2012). Various variations of this composition have been attempted, of which the most studied have been the replacement of lanthanum with a different lanthanide Xing et al. (2012); Jha and Awana (2013); Demura et al. (2013); Jha et al. (2013); Zhai et al. (2014); Yazici et al. (2013).

Interest in this family of superconductors was sparked by the theoretical prediction of the presence of a charge density wave instability in La(O,F)BiS2 Wan et al. (2013); Yildirim (2013), and its observation in EuFBiS2 Zhai et al. (2014). Moreover, ferromagnetism was measured in the superconducting phase of Ce(O,F)BiS2, in agreement with ab-initio calculations Xing et al. (2012); Morice et al. (2016). Finally, quantum critical fluctuations of the magnetic moments were measured in CeOBiS2 Higashinaka et al. (2015).

La(O,F)BiS2 diplays a large range of superconducting transition temperatures between 2.5 and 11.5 K Mizuguchi et al. (2014); Zhang et al. (2016). Its change with doping was measured, and gives a dome structure centred on $x=0.5$ Deguchi et al. (2013). The synthesis method has been shown to be crucial in changing $T_{c}$: the crystals synthesised using a simple solid state reaction at ambient pressure (AP) have $T_{c}\sim 3$ K, while the ones having gone through an extra high pressure annealing stage (HP) have $T_{c}\sim 10$ K Mizuguchi et al. (2012b). Interestingly, this change is revertible: HP samples annealed again at ambient pressure recover the structure and properties of AP samples Kajitani et al. (2014). The difference in crystal structure between AP and HP samples is hard to characterise, as X-ray peaks of HP samples are much broader than the ones of AP samples Mizuguchi et al. (2012b); Deguchi et al. (2013). This difference has however been linked to strain along the c-axis Kajitani et al. (2014), and to a change from a larger positive Hall resistance in AP samples to a smaller negative one in HP samples Pallecchi et al. (2014).

The evolution of La(O,F)BiS2 under pressure has been measured for both AP and HP samples. $T_{c}$ in AP samples exhibit a sharp increase under pressure from around 3 K to around 9.5 K at 0.7 GPa Wolowiec et al. (2013); Tomita et al. (2014). In comparison, in HP samples, $T_{c}$ rises very slightly with pressure Kotegawa et al. (2012); Wolowiec et al. (2013). This change in $T_{c}$ has been shown to coincide with a structural phase transition from tetragonal to monoclinic Tomita et al. (2014). However relating this to the crystal structure of the HP samples is difficult since X-ray and neutron measurements both gave a tetragonal structure for HP samples Mizuguchi et al. (2012b); Lee et al. (2013)

The prediction of a charge density wave instability in La(O,F)BiS2 Wan et al. (2013); Yildirim (2013) led to in-depth studies of its vibrational modes using neutron diffraction Lee et al. (2013); Athauda et al. (2015). No phonon anomaly in either doping or temperature change was detected, which suggested that the electron-phonon coupling could be much weaker than the calculated values Lee et al. (2013). Furthermore, neutron diffraction data yielded the presence of local charge fluctuations Athauda et al. (2015).

Angle-resolved photoemission spectroscopy (ARPES) was performed on La(O,F)BiS2, Ce(O,F)BiS2 and Nd(O,F)BiS2 Zeng et al. (2014); Ye et al. (2014); Terashima et al. (2014); Sugimoto et al. (2015). Large discrepancies with electronic structure calculations were obtained, particularly the absence of the central nested Fermi pocket Zeng et al. (2014); Ye et al. (2014); Sugimoto et al. (2015). This was attributed to a level of doping lower than expected Ye et al. (2014); Sugimoto et al. (2015). This was confirmed by measurements on precisely calibrated samples Terashima et al. (2014). Finally, electronic correlations were found to be weak in these systems Ye et al. (2014); Terashima et al. (2014)

The diversity in the phenomenology of these materials drove theoretical studies to probe the possibility of unconventional pairing. Spin-fluctuation mediated pairing has been calculated using the random-phase approximation, Martins et al. (2013); Zhou and Wang (2013); Wu et al. (2014), numerics Liang et al. (2014); Liu and Feng (2014) and functional renormalisation group Yang et al. (2013).

Ab-initio calculations outputted large electron-phonon couplings, between 0.8 and 0.85 Wan et al. (2013); Yildirim (2013); Li et al. (2013). Two calculations found unstable phonon modes around ($\pi$,$\pi$,0) Wan et al. (2013); Yildirim (2013), leading to a commensurate charge density wave Wan et al. (2013); Yildirim (2013). The phonon dispersion along the $c$ axis was found to be small, highlighting the small magnitude of the interlayer coupling across the buffer layers Wan et al. (2013). All these studies concluded that superconductivity in La(O,F)BiS2 is strongly coupled and conventional Wan et al. (2013); Yildirim (2013); Li et al. (2013). The same conclusion has been reached in La(O,F)BiSe2 Feng et al. (2014)

In this paper, we study conventional superconductivity in LaO0.5F0.5BiS2, using density functional theory for superconductors Lüders et al. (2005); Marques et al. (2005). We first confirm the arising of a commensurate charge density wave instability, in accordance with previous studies Wan et al. (2013); Yildirim (2013). Using a novel integration scheme for the electron-phonon coupling, we find that its strength is much lower than previously calculated, due to our improved convergence of density of state calculations. We finally conclude that conventional phonon-mediated pairing cannot explain the high $T_{c}$s observed in this material.

II Methods

Electronic structures, dynamical matrices and electron-phonon coupling coefficients were calculated using norm-conserving pseudopotentials and a plane wave basis set as implemented in quantum espresso Giannozzi et al. (2009). The plane-wave energy cutoff for the wavefunctions was set to 90 Ry. We also used the all-electron full-potential linearised augmented plane-wave code Elk Elk (2004). Electron dielectric functions were calculated within the frequency-dependent random-phase approximation.

We employed the exchange-correlation potential implemented by Perdew and Zunger based on the local density approximation Perdew and Zunger (1981). The relativistic effects are mostly treated at the scalar-relativistic level Koelling and Harmon (1977), whereas the spin-orbit coupling was explicitly included for re-examining the charge instability (see below).

We used density functional theory for superconductors to calculate superconducting transition temperatures ($T_{c}$) from first principles Lüders et al. (2005); Marques et al. (2005). It is an ab-initio theory which allows one to calculate the conventional strong-coupling $T_{c}$ within the level of the Migdal-Eliashberg theory Migdal (1958); Eliashberg (1960, 1961) and without ajustable parameters.

It consists in solving self-consistently the gap equation:

[TABLE]

where $n$ and $n^{\prime}$ run over bands, k and $\textbf{k}^{\prime}$ over crystal momentum, $\beta$ is the inverse temperature, $\Delta_{n\textbf{k}}$ is the superconducting gap function, and $E_{n\textbf{k}}=\sqrt{(\epsilon_{n\textbf{k}}-\mu)^{2}+\Delta_{n\textbf{k}}^{2}}$ where $\epsilon_{n\textbf{k}}$ is the Kohn-Sham one-particle energy and $\mu$ is the chemical potential for the normal state. In subsequent calculations, $\mu$ is approximated to the Fermi level for the normal state $E_{F}$, where the temperature dependence is ignored. $\mathcal{Z}$ and $\mathcal{K}$ are exchange-correlation kernels. $\mathcal{Z}$ represents the mass renormalisation of the normal-state band structure due to phonon exchange, while $\mathcal{K}=\mathcal{K}^{\text{ph}}+\mathcal{K}^{\text{el}}$ represents electron-phonon and electron-electron interactions.

This gap equation is derived by using an extended version of the standard finite-temperature density-functional procedure Mermin (1965) consisting in minimising a grand canonical potential function of three densities: the density of electrons, of nuclei, and the superconducting order parameter Lüders et al. (2005). Choosing a system of non-interacting electrons but interacting nuclei as the Kohn-Sham system, one obtains the Kohn-Sham equations in the Bogoliubov-de Gennes form. Neglecting the feedback effect of the superconducing gap on the normal electronic states finally gives the gap equation Lüders et al. (2005).

The kernels of the gap equation, defined by the previous procedure as functional derivatives of terms of the grand canonical potential, are derived perturbatively using many-body theory Görling and Levy (1994); Lüders et al. (2005). We employed the phononic kernels $\mathcal{K}^{\text{ph}}$ and $\mathcal{Z}^{\text{ph}}$ averaged over $k$-point defined in equation (23) in Marques et al. (2005) and equation (40) in Akashi and Arita (2013a), respectively. The phonon and electron-phonon coupling properties were calculated using density functional perturbation theory Baroni et al. (2001); Giannozzi et al. (2009). The pairing and mass-renormalisation effects are thus treated at the level of the Migdal-Eliashberg theory with varying density of states Pickett (1982).

The electronic kernel $\mathcal{K}^{\text{el}}$ is calculated using the frequency-dependent dielectric function Akashi and Arita (2013b, 2014); Akashi et al. (2015), thus including dynamical Coulomb effects Takada (1978).

III Results

III.1 Structural instability

We calculated the phonon dispersion in the tetragonal unit cell using density functional perturbation theory on a sampling of the Brillouin zone of $4\times 4\times 2$ q-points with quantum espresso. Hereafter we refer to the phonon wave number by q. The obtained dispersion features a structural instability which extends all along the $\Gamma$-M line, matching previous studies.

In order to test the robustness of this instability, both to other approximations and to the inclusion spin-orbit coupling, we recalculated this dispersion using an all-electron calculation with and without spin-orbit coupling, and a q-grid of $2\times 2\times 1$ points. This yielded two dispersions, which are similar, with a slightly reduced depth of the imaginary frequencies from $-20$ meV to $-14$ meV. This confirms that this instability is robust.

We obtained the fully relaxed structure by using a unit cell of dimensions $\sqrt{2}a\times\sqrt{2}a$ where $a$ is the unit cell parameter of the tetragonal structure. For convenience we call this new crystal structure the $\sqrt{2}\times\sqrt{2}$ structure.

We calculated the phonon dispersion in this $\sqrt{2}\times\sqrt{2}$ structure (Figure 1). The obtained dispersion does not feature the instability 111We find imaginary modes in the small region around the Gamma point, which could be due to insufficient energy cutoff for charge density.. We therefore use this structure to calculate electron-phonon coupling coefficients.

III.2 Electron-phonon coupling

The effect of electron-phonon coupling on conventional superconductivity is well represented by the Eliashberg function Scalapino (1969):

[TABLE]

where $\omega$ is the frequency, $n$ and $n^{\prime}$ run over electronic bands, $\nu$ over phononic bands, k and q over the Brillouin zone, and $g$ are the electron-phonon coupling coefficients. Integrating the Eliashberg function over frequency gives the electron-phonon coupling factor $\lambda$:

[TABLE]

The characteristic phonon frequency $\omega_{ln}$ is defined by:

[TABLE]

The effects of the electron-phonon coupling on $T_{c}$ are well quantified by the McMillan-Allen-Dynes formula McMillan (1968); Allen and Dynes (1975):

[TABLE]

where $\mu^{*}$ is an adjustable pseudopotential parameter representing the effective electron-electron Coulomb repulsion suppressing the pairing instability Morel and Anderson (1962).

In practice, the integration of equation 1 is carried out by approximating the singular delta function by a smeared function $\tilde{\delta}$ such as the Hermite-Gaussian Methfessel and Paxton (1989). However, this treatment can introduce two types of quantitative errors. First, the integrated value converged with respect to the $k$-point density suffers from a systematic error due to the width of the smearing function. Second, the convergence with respect to the $k$-point density tends to be very difficult to achieve, especially when the smearing width is small.

The magnitude of the error is well quantified by examining the electronic density of states $N(E)=\sum_{n\textbf{k}}\delta(E-\epsilon_{n\textbf{k}})$. We define $\tilde{N}$ as the density of states calculated using this smearing scheme. Here and hereafter, the functions and values with tilde denote the ones calculated with the smearing scheme. Equation 1 indicates that $\tilde{\alpha}^{2}F$ scales with $\tilde{N}(\tilde{E}_{F})^{2}$ in the limit where the electron-phonon coupling coefficients $g$ are independent of $n$, $n^{\prime}$, k and q 222In the Phonon package of quantum espresso, the prefactor $1/N(E_{F})$ in equation 1 is replaced by $1/\tilde{N}(E_{F})$, hence in this case we only obtain a $\tilde{N}(E_{F})$ dependence when summing the delta functions.. To estimate the impact of the above quantitative errors, we examined the dependence of the electronic density of states at the Fermi energy $\tilde{N}(\tilde{E}_{F})$ with respect to the smearing width and $k$-point grid (Figure 3). We obtained a well-converged value for $N(E_{F})$ using a very dense $k$-point grid and the tetrahedron integration method, which we refer to as the “exact” $N(E_{F})$ hereafter. We found that $\tilde{N}(\tilde{E}_{F})$ differs from the exact value by approximately 20% at the Gaussian width 0.5 eV, which is within a standard range of the width. Moreover, we found slowly convergent behaviour for widths lower than 0.1 eV. Notably, the values of $\tilde{N}(\tilde{E}_{F})$ are systematically larger than the exact one, indicating that $\tilde{\alpha}^{2}F$ can suffer from systematic overestimation in this material.

In order to mitigate the above-mentioned overestimation we modified the $\alpha^{2}F$ formula in practice, following recent technical developments Koretsune and Arita (2016), and obtained the rescaled Eliashberg function 333In the case of quantum espresso, because the program uses the prefactor $1/\tilde{N}(E_{F})$ in equation 1, we need to use the following formula:

$\displaystyle\alpha^{2}F_{\text{resc}}(\omega)=$ $\displaystyle\frac{\sum_{\textbf{q}}\sum_{nn^{\prime}\textbf{k}}\delta(\epsilon_{n\textbf{k}+\textbf{q}}-E_{F})\delta(\epsilon_{n^{\prime}\textbf{k}}-E_{F})}{\sum_{\textbf{q}}\sum_{nn^{\prime}\textbf{k}}\tilde{\delta}(\epsilon_{n\textbf{k}+\textbf{q}}-\tilde{E}_{F})\tilde{\delta}(\epsilon_{n^{\prime}\textbf{k}}-\tilde{E}_{F})}$

$\displaystyle\times\frac{\tilde{N}(\tilde{E}_{F})}{N(E_{F})}\tilde{\alpha^{2}F}(\omega)$

(5) :

[TABLE]

The numerator is replaced to the exact $N(E_{F})^{2}$ since it corresponds to the dense k and q-points limit. Note that regardless of the k-point density, this modified formula gives us the exact value in the limit where the electron-phonon coupling coefficients $g$ are independent of $n$, $n^{\prime}$, k and q. With this formula, as shown in figure 4, we obtained values of $\lambda$ almost independent of the smearing width.

We confirmed these results by calculating $\alpha^{2}F$ with the optimized tetrahedron method Kawamura et al. (2014) (Figure 2). We used four different q-point meshes, two with $3\times 3\times 3$ points (and a $18\times 18\times 6$ k-point mesh) and two with $4\times 4\times 2$ points (and a $24\times 24\times 4$ k-point mesh), one shifted along the $z$ axis and the other along all the axes 444The unshifted q-mesh cannot be employed for the tetrahedron integration since it suffers from numerical divergence at $\textbf{q}=0$. Note that this is an artefact due to the discretisation and can be avoided if the q-point integral around $\textbf{q}=0$ is calculated analytically 555We have also employed the rescaling formula (Equation (6)) for $\alpha^{2}F_{tetra}$ to accelerate the k-point convergence. In this case, the denominator of the prefactor is calculated with the k and q-point mesh employed for the integration in equation (1).. We name $\alpha^{2}F_{\text{tetra}}$ the Eliashberg function obtained using the all-shifted $4\times 4\times 2$ q-point mesh. We obtained $\lambda=0.39$ and $\lambda=0.58$ for the first sampling, with $z$-shifted and all-shifted meshes, respectively, and $\lambda=0.48$ and $\lambda=0.47$ for the second sampling. The values obtained with $3\times 3\times 3$ points differ clearly, indicating a lack of convergence, while the ones obtained with $4\times 4\times 2$ points are both very close and between the two previous ones. This indicates a strong convergence. Moreover, the values obtained with the $4\times 4\times 2$ q-point mesh match closely the values obtained by integrating $\alpha^{2}F_{\text{resc}}$ in the limit of small smearing (Figure 4), which is the exact limit. The small deviation between the values of $\lambda$ indicates that the converged $\lambda$ is about or below 0.5.

The results obtained for $\omega_{ln}$ are less conclusive: indeed the values obtained for the two $3\times 3\times 3$ q-point grids are 116 K and 153 K, while the ones obtained using the two $4\times 4\times 2$ q-point grids are 98 K and 143 K. Let us note that $T_{c}$ is affected exponentially by $\lambda$ whereas it is affected only linearly by $\omega_{ln}$ (Equation (3)), which means that the convergence of $\lambda$ gives us the order of magnitude of $T_{c}$. In the following subsection, we show that this order of magnitude is clearly different from the one observed experimentally.

III.3 Superconducting transition temperature

We calculated the superconducting transition temperature ($T_{c}$) of LaO0.5F0.5BiS2 using density functional theory for superconductors 666We sampled the Brillouin zone and the different bands using a weighted random set of points Marques et al. (2005). The standard error on the value of the gap was found to be 18.3%.. We ran a “static” calculation which only uses the zero frequency dielectric function, and a “dynamical” calculation which uses the full frequency-dependent dielectric function, both calculated using a $4\times 4\times 2$ k-point grid. We used $\alpha^{2}F_{\text{tetra}}$ in both calculations. The results are summarised in table 1.

It is interesting to compare this result to the output of the McMillan-Allen-Dynes formula McMillan (1968). We calculated $T_{c}$ with this formula using the parameters corresponding to $\alpha^{2}F_{\text{tetra}}$ for a large range of values of $\mu^{*}$. We then can obtain the $\mu^{*}$ corresponding to the value of $T_{c}$ calculated using density functional theory for superconductors. We obtain values of $\mu^{*}$ larger than in most conventional superconductors Carbotte (1990); Akashi et al. (2012) (Table 1). This is reminiscent of other ab-initio studies of superconductors Akashi et al. (2012).

The value of $T_{c}$ calculated here are very distant from all the values of $T_{c}$ measured in this compound, which are between 2.5 and 11.5 K Mizuguchi et al. (2014). We therefore conclude that conventional phonon-coupling superconductivity cannot explain the experimental data on LaO0.5F0.5BiS2.

IV Conclusion

We calculated phonon frequencies in the original unit cell of LaO0.5F0.5BiS2, using two different approximations and spin-orbit coupling, and found a structural instability close to the M point of the Brillouin zone. This confirms previous findings Wan et al. (2013); Yildirim (2013).

The electron-phonon coupling was calculated in a relaxed $\sqrt{2}\times\sqrt{2}$ structure using the standard smearing integration scheme and obtained values for $\lambda$ with a strong smearing width dependence (Figure 4). These results match previous calculations well Yildirim (2013).

We then used a cutting-edge rescaling technique Koretsune and Arita (2016) to mitigate the influence of errors in the calculation of the density of states, which resulted in a much smaller dependence of $\lambda$ on smearing width (Figure 4). The obtained $\lambda$ matches closely results obtained with a novel tetrahedron method Kawamura et al. (2014). We conclude that $\lambda<0.5$.

Finally, using the Eliashberg function calculated using the tetrahedron method, we calculated $T_{c}$ using density functional theory for superconductors Lüders et al. (2005); Marques et al. (2005). We obtained $T_{c}=0.4$ K. This corresponds to a pair-breaking electron-electron repulsion $\mu^{*}=0.16$, which is larger than conventional values. We thus conclude that conventional phonon-mediated pairing cannot explain the high $T_{c}$s observed in this material. This is an important step toward the clarification of the origin of superconductivity in BiS2 superconductors.

Acknowledgments

We would like to thank K. Kuroki and H. Ikeda for interesting discussions. We acknowledge support from EPSRC, Corpus Christi College, and JST, PRESTO and JSPS KAKENHI grants number JP15H03696, JP16H00924, JP16H06345 and JP15K20940. This work was supported by MEXT Element Strategy Initiative to Form Core Research Center in Japan. Part of the calculations were performed at the Supercomputer Center at the Institute for Solid State Physics in the University of Tokyo.