Spin-orbit coupling, minimal model and potential Cooper-pairing from repulsion in BiS$_2$-superconductors

We develop the realistic minimal electronic model for recently discovered BiS$_2$ superconductors including the spin-orbit coupling based on a first-principles band structure calculations. Due to strong spin-orbit coupling, characteristic for the Bi-based systems, the tight-binding low-energy model necessarily includes $p_x$, $p_y$, and $p_z$ orbitals. We analyze a potential Cooper-pairing instability from purely repulsive interaction for the moderate electronic correlations using the so-called leading angular harmonics approximation (LAHA). For small and intermediate doping concentrations we find the dominant instabilities to be $d_{x^2-y^2}$-wave, and $s_{\pm}$-wave symmetries, respectively. At the same time, in the absence of the sizable spin fluctuations the intra and interband Coulomb repulsion are of the same strength, which yields the strongly anisotropic behaviour of the superconducting gaps on the Fermi surface in agreement with recent ARPES findings. In addition, we find that the Fermi surface topology for BiS$_2$ layered systems at large electron doping can resembles the doped iron-based pnictide superconductors with electron and hole Fermi surfaces with sufficient nesting between them. This could provide further boost to increase $T_c$ in these systems.


I. INTRODUCTION
Recent discovery of BiS 2 -layered superconductors with T c up to 10K has attracted significant attention due to striking similarities of their crystal structure with cuprate and iron-based high-T c superconductors [1][2][3][4][5][6][7][8][9][10][11] . According to the density functional theory (DFT) calculations, the parent compounds of the BiS 2 -based superconductors are semiconductors and the metallic behavior is achieved by electron doping of the conduction band, which mainly consists of Bi 6p orbitals 1,12 .
Although the origin of superconductivity is not yet clear in these compounds, it was attributed to the electron-phonon interaction [13][14][15] . This seems reasonable in view of the low superconducting transition temperature and weaker electronic correlations of p-orbitals than in their 3d-counterparts. However, in the recent neutron scattering experiment, the observed almost unchanged low-energy modes indicated that the electron phonon coupling could be weaker than expected 16 . In addition, the large ratio 2∆/T c may suggest that the pairing mechanism is unconventional 17,18 . In addition, no isotope effect on superconducting transition temperature was found 19 .
Very recently, using 'ab-initio' approach the strength of the electron-phonon coupling in LaO 0.5 F 0.5 BiS 2 was estimated to be λ < 0.5 20 , which is too low to explain the superconducting transition temperature in this system. Electron-electron correlations can, in principle, be also responsible for the Copper pairing and there were several theoretical studies about possible pairing symmetries in these new superconductors arising due to repulsive in-teractions [21][22][23][24] . However, as the nominal compositions of the superconducting materials referred to the significant electron doping 3,11 , all these previous studies concentrated on the high electron doping region where the electronic structure was featured with large Fermi surfaces in close vicinity to Van Hove singularity. However, very recently, it was reported that the actual electronic filling in these systems is much lower than the nominal one and there are only two small electron pockets around the X points of the Brillouin Zone 25,26 corresponding to x = 0.22 27 . Most importantly, the recent observation of the strongly anisotropic superconducting gaps on the electron pockets in NdO 0.71 F 0.29 BiS 2 using laser angle resolved photoemission spectroscopy (laser ARPES) suggests that the pairing mechanism for this system is likely an unconventional one and was attributed to be a result of the multiple paring interactions 28 .
Here, based on the fully relativistic ab-initio band structure calculations, we develop the realistic minimal electronic model for recently discovered BiS 2 superconductors including the spin-orbit coupling, relevant for Bi-based systems. We show that the tight-binding lowenergy model necessarily includes all three p-orbitals, p x , p y , and p z . Using angular harmonics approximation (LAHA) 29 we analyze a potential Cooper-pairing instability from purely repulsive interaction concentrating on the case of weak to intermediate electronic correlations strength for small and intermediate doping levels. Similar to the previous spin fluctuations based studies 12,30 we find global d x 2 −y 2 (B 1g )-wave, and anisotropic s-wave (A 1g ) symmetries, respectively, to be the dominant symmetries. However, our results show that for weak corre-lations, i.e. when the intraband and the interband repulsion are of similar strengths, the gap cannot be described by the single harmonics of the corresponding functions of the A 1g or B 1g irreducible representations typical for the case of the strong spin fluctuations. Instead, each gap acquires a strong angular dependence in the form of the multiple harmonics on each of the Fermi surface pockets allowed by the global symmetry. This is a consequence of the fact that the electronic system reduces the net repulsion. We anaylze the form of these harmonics and show that they are in agreement with recent experimental ARPES observation 28 . We further investigate the Fermi surface topology for BiS 2 layered systems at large electron doping and show its similarity with the doped iron-based pnictide superconductors with electron and hole Fermi surfaces with sufficient nesting between them. This could provide further boost to increase T c in these systems.

II. TIGHT-BINDING MODEL
The fully relativistic band structure calculation of the LaOBiS 2 system was performed within DFT using the Perdew-Burke-Ernzerhof exchange-correlation functional as implemented in WIEN2K program 31 . The muffintin radius of each atom R MT was chosen such that its product with the maximum modulus of reciprocal vectors K max become R MT K max = 7.0. The lattice parameters were taken from Ref. 32 and the corresponding Brillouin zone was sampled using a 12×4 k−mesh. From this band calculation, as shown in Fig.1, we constructed a 12-band tight binding (TB) model using the maximally localized Wannier functions [33][34][35] with Bi 6p (p x , p y , and p z ) orbitals as projection centers. The resulting TB model was further reduced to a 3-band model plus spin-orbit coupling by ignoring the inter-layer hopping between the BiS 2 layers. Here, t ij;νν = t(x i − x j , y i −y j ; νν ) are the nearest and next-nearest neighbor hopping parameters obtained from projecting the results of ab initio calculations. Their values are shown in Table  I where ν, and ν refer to the orbital indices. Thus, in the Momentum space the Hamiltonian can be written aŝ with the hopping integrals T νν =t ν 0,0 + 2t ν 1,0 cos k x + 2t ν 0,1 cos k y + 2t ν 1,1 cos(k x + k y ) + 2t ν 1,−1 cos(k x − k y ), T νν =T ν ν = t νν 0,0 + 2t νν 1,0 cos k x + 2t νν 0,1 cos k y + 2t νν 1,1 cos(k x + k y ) + 2t νν 1,−1 cos(k x − k y ).
Here t ν x,y and t νν x,y denote intra-and interorbital electron hopping, as determined in Table I. On top of that, the on-site spin-orbit (SO) coupling Hamiltonian, H SO = λ L · S, written for the bismuth p-orbitals have the following form: hereafter, we set λ = 0.874eV as a fitting parameter to match the results of the effective model to the DFT band structure. The full electronic structure of the resulting tight-binding Hamiltonian is shown in Fig.2(a). For the low filling the band structure can be indeed reproduced roughly by taking account only the p x and p y orbitals and the parameters we find are similar to the ones found previously 12 . We should note that the proper inclusion of the spin-orbit coupling does require to involve the p zorbital into consideration. This is valid even for the small fillings that always there is an admixture of the p z orbitals to the lowest band due to spin-orbit coupling. For the small fillings the Fermi surface of the BiS 2layer consists of the two electron pockets centered near (±π, 0) and (0, ±π) points of the BZ, which resembles some of the iron-based electron-doped chalcogenide superconductors. Upon doping the electron pockets become bigger and at x = 0.5 the system undergoes the transition from the small electronic Fermi surface pockets to the large Fermi surface, similar to Ref. 12 . Due to spin-orbit coupling the contribution of p z -orbital to the conducting band increases and becomes quite significant for x > 0.5. Studying the evolution of the electron structure for larger doping we observe that for the doping range n > 1.5, the Fermi surface topology resembles strongly the one found in the iron pnictide superconductors (see Fig.2(c)). That corresponds to the weakly nested electron and hole Fermi surfaces centered near the corresponding points of the BZ. Such a nesting tends to boost the inter-band electron-electron scattering and favor non-phononic mechanism of superconductivity.

III. LINDHARD RESPONSE FUNCTION
As mentioned in the introduction, the origin of superconductivity in BiS 2 layers is still debated. Nevertheless, given the similarity of the BiS 2 -systems and the ironbased superconductors regarding crystal structure and Fermi surface topology, it is tempting to apply the concept of Cooper-pairing from repulsion to these systems  TABLE I. Hopping parameters t(xi − xj, yi − yj; νν ) = t(∆x, ∆y; νν ) (in meV) for the simplified two-dimensional model. Note that t(∆x, ∆y; νν ) =t(-∆x, −∆y; νν ) and t(∆x, ∆y; νν ) =t(∆x, ∆y; ν ν). as well, given much weaker correlations strength in Bibased materials. Previously, assuming spin fluctuations scenario and the presence of the disconnected Fermi surfaces for the small electron doping the d x 2 −y 2 -wave symmetry for the BiS 2 systems was obtained 12,30 . At the same time, recent ARPES experiments point towards strongly anisotropic superconducting gap 28 , which would not agree with the simple spin-fluctuation mediated scenario with strong interband repulsion. Here, we concentrate on the more realistic situation of the weak spin fluctuations and analyze the potential Cooper-pairing instabilities for the case of disconnected Fermi surface pockets using the leading angular harmonics approximation (LAHA), developed previously and specially designed for the systems when the spin fluctuations are either intermediate or weak. 29,36,37 . To search for the potential instabilities of the electron system, we perform the calculations of the Lindhard response function, which gives information on the potential instabilities of electronic system. For the multiorbital systems Based on the Lindhard spin response function theory, the bare spin susceptibility is defined by where ν i runs over orbitals, andŜ s are the spin operators given byŜ  whereσ s are Pauli matrices. Therefore, by applying Wick's theorem and neglecting contractions leading to q = 0, one can find that the different components of the bare spin susceptibility, χ ss 0 , in the frequency domain is with the Green's functions Introduced by here we set iω = iω n + iΩ. Since the Green's functions are not diagonal in the orbital basis it is convenient to move to the band representation, and find where µ being the band index, and the matrix elements ζ µ kνσ , are connecting the band to the orbital basis: c kνσ (τ ) = µ ζ µ kνσ b kµ (τ ). By performing the Matsubara frequency sum over ω n , yields defined the transverse susceptibility as follows and by consideringσ = −σ, we find for the longitudinal susceptibility. The results of the calculations for the two representative concentrations of the model are present in Fig3 (a),(b). For n = 0.2 and n = 1.4 the effect of the spin orbit-coupling is an introduction of the weak easy-plane anisotropy, i.e. χ +− > χ zz , which appears to be an overall prefactor. Therefore we discuss only the behaviour of χ +− and the situation for χ zz is similar. For n = 0.2 we clearly observe the intraband (small q) and interband (large Q) scattering of the electron pockets of similar strength. This occurs due to similar orbital content of the pockets. Nevertheless, the clear separation of the peaks allows for the purely electronic mechanism of the Cooper-pairing even in the presence of weak or moderate correlations as we show later. At the same time, much stronger electronic response is found for n = 1.4 where one clearly finds the strong scattering between electron and hole pockets near Q = (π, 0)[(0, π)]. It dominates the instability of the electronic gas at this filling arising from the logarithmic divergence of the Lindhard response function near this wave vector due to near nesting of the electron an hole bands separated by Q. This is somewhat reminiscent of the electronic structure of the iron-based superconductors and worth studying further.

IV. RANDOM PHASE APPROXIMATION
In order to study the superconducting state, we introduce an intra-site Hubbard model, including intra and inter-orbital interactions (U and U , respectively), as well as the Hund's rule J and pair hopping energy J , to compute the total electronic response within the random Phase approximation (RPA) approach. In this respect, the real space picture of the interacting Hamiltonian is given by 23 Therefore, in the band representation the total Hamiltonian readŝ where the quartic terms describe the scattering of pairs (k ↑, −k ↓) on the pocket µ to (k ↑, −k ↓) on the pocket µ , with subsequent scattering amplitude: The peaks in the bare susceptibility are enhanced within the RPA only for the larger value of n. If the incommensurate peak at (π/3, π/3) and the one close to (π, π) for n = 0.2 are only slightly enhanced, the peaks close to the antiferromagnetic wave-vector Q = (π, 0) for n = 1.4 are enhanced significantly due to the fact the the pockets have a different character. While the pockets near the Γ and M pockets are hole like, the pockets near X and Y -pints are electron like. Momentum ranges are given by (−π ≤ q x/y ≤ π).
in which, the vertex Γ sq tp (k 1 , k 2 ) describes the amplitude of the scattering in orbital basis, and it is a linear combination of the interactions: U , U , J, and J . Here, ζ p kµ = ζ µ * kp are the matrix elements connecting the orbital and band basis, and they represent the orbital contributions to the Fermi surface, see Fig. 2(c&d).
The final RPA result of the susceptibility can be addressed in the form of Dyson equations as follows: where (U µ ) u,v w,z are the coefficients of the interacting Hamiltonian arranged into the three orbital basis matrix, and the repeated indices are summed over. For the transversal susceptibility the above equation results in where the non-zero components of U µ are given by however, for the longitudinal susceptibility we find where and finally Despite the fact that the correlations should be moderate in BiS 2 systems, we still include the RPA corrections assuming the strength of the interactions to be U < t. Nevertheless, as the Lindhard response function shows, there are several intraband and interband scattering peaks, well separated in momentum space. Therefore, we expect them to be enhanced differently once correlations are included in the particle-hole channel. In particular, we see that the RPA spin susceptibilities, despite of moderate interactions, further strengthen the features, present in the Lindhard response function. For n = 0.2 the interband and intraband scattering are only slightly enhanced, which indicates very weak instability of the electron gas in the magnetic channel. The interband and intraband magnetic fluctuations are by far not enough to generate the true magnetic instability. They are, however, sufficiently strong to generate the so-called nodeless d x 2 −y 2 -wave superconductivity with sign change of the gap between electron pockets. This type of superconductivity was originally proposed for relatively small electron filling within spin fluctuation scenario 12 . The crucial difference is, however, that this gap possesses strong angular dependence on the Fermi surface sheets making it difficult to identify the global d x 2 −y 2 −wave symmetry. Furthermore, for large n = 1.4 the instability with respect to the (π, 0) or (0, π) density wave type ordering is stronger, which is somewhat similar to the iron-based superconductors. Such an instability is known to favor the so-called s ± -wave symmetry of the superconducting order parameter, driven by short-range magnetic fluctuations. We analyze these instabilities quantitatively in the section using the so-called leading angular harmonic approximation, introduced previously 29,36,37 . As we show later, even in this case there are higher angular harmonics in the superconducting gap on the Fermi surface.

V. LEADING ANGULAR HARMONICS APPROXIMATION
The LAHA approximation was developed for the ironbased superconductors 29,36,37 as a method to analytically solve and characterize the superconducting gap equation. This method is particularly well suited once the interactions are moderate, which do not give rise to a pronounced enhancement of the spin fluctuations as it is the case in LiFeAs 29 . In this section we extend it to the BiS 2 superconductors. The basic idea of LAHA is that the Γ µµ , defined in eq. (14), are dependent on the angles along different Fermi surface sheets, which are well separated in the BZ. Therefore, they can be defined as simple functions of the momenta k and k 29,36,37 . In particular, one can decompose the Γ µµ (k, k ) into representations of the tetragonal space group A 1g (s-wave) and B 1g (d-wave), i.e., with Ψ s and Ψ d being the basis functions of the A 1g and B 1g representations, respectively. Finally, the Γ (A1g) and Γ (B1g) can be expanded to model the angular dependence of the pair-scattering. In this work, the LAHA approximation will be adapted and implemented to analyze the superconducting instabilities of BiS 2 superconductors in two limiting cases: the small (n < 0.5), and the large (n > 1.5) electron fillings. The Fermi Surface topology is characterized in both cases by the presence of disconnected electron and/or hole pockets centered near the high-symmetry points of the first BZ. Furthermore, we expect the spin singlet state to dominate the actual pseudo-spin-singlet Cooper-pairing and that the Cooperpairs can be further characterized as even-parity spin singlet solutions. This is due to the fact that the spin-orbit coupling only introduces the overall difference between the longitudinal and the transverse components, which is weakly momentum dependent. Thus, one could consider the Cooper-pairing can be still determined as an even parity wave functions in the psudospin basis. Furthermore, as the admixture of the p z orbitals is also often weak, the spin remains in most of the cases a good quantum number.
Here the upper (lower) sign corresponds to pocket 1 (2), φ kei (φ k hi ) represents the angles along the i th electron (hole) pocket, and Γ s and Γ d refer to the extended s-wave and d x 2 −y 2 -wave Cooper-pairing vertex, respectively. As the correlations are weak, it is instructive to consider the LAHA projections not for the RPA but only restricting to the second order perturbation for the interactions. They are shown in Table II.   Observe that the constant part of the interaction in the A 1g -pairing channel is larger than in the B 1g -one. Nevertheless, the constant part of the interaction is always repulsive for the Fermi surface topology for small n, while the interband part of the repulsion is pair-building for the B 1g -symmetry, due to the change of the sign of the order parameter between the (π, 0) and the (0, π) pockets. Nevertheless, the inter and intraband repulsions are the same due to the fact that both bands are electronlike. An inclusion of the RPA corrections in the particlehole channel due to interband nesting helps in achieving U e1e2 >Ũ e1e1 ,Ũ e2e2 . Nevertheless, all interactions are still very similar to each other as the spin fluctuations are not strongly enhanced, see Fig.3(c). As a result the solution for the superconducting gap acquires a strong angular dependence in the B 1g -channel in the form ∆ e1 (φ) = ∆ e +∆ e cos 2φ +∆ e cos 4φ + ...
where ... stands for the higher harmonics of cos 6θ and cos 8θ, which we have not explicitly included. Furthermore, it is important to realize that for the intraband and interband repulsive interaction of similar strength, the superconductivity occurs only if the angular dependent harmonics are of the same magnitude as the constant parts of the gap 29 , which yields nearly nodal behaviour of the superconducting gap on each of the Fermi surface sheets, i.e. we find ∆ e ∼∆ e ∼∆ e . This is in qualitative agreement with the experimental data 28 . Nevertheless, the most important prediction of the superconducting gap to belong to the B 1g -irreducible representation is the antiphase character of the constant magnitudes and cos 4θ harmonics on the two electron pockets, which still needs to be tested experimentally. Observe, that near-nodal d x 2 −y 2 -wave solution dominates for the Fermi topology consisting of the electron pockets and the angular harmonics are increasing in magnitudes with increasing n. Importantly, this is not altered by the inclusion of the spin-orbit coupling and this symmetry appears to be still the most dominant one for n < 0.4. For larger doping the Fermi surface undergoes a so-called Lifshitz transition and changes its topology 12,23 . Although for n > 0.4 the unconventional superconducting mechanism cannot be fully excluded, the superconducting mechanism cannot be purely repulsive as the momentum structure of the response function is less pronounced. Therefore it will likely lose against simple s-wave symmetry, driven by the electron-phonon interaction. For a larger doping level n, we find a richer and more complicated Fermi Surface: the electron pockets described in the previous section have now merged giving rise to two hole pockets located at Γ and M . On the other hand, elliptical electron pockets from the upper energy band appear at X and Y . All pockets have a similar radius, and are separated by a large wave vector Q = (π, 0) or (0, π).
Due to the presence of hole pockets, we need to take into account hole-hole Γ hihj , and electron-hole Γ eihj scattering vertices. On the other hand, all four pockets are highly symmetric now. Therefore, we neglect the harmonics higher than cos 2φ on the electron pockets for simplicity:   have, Finally, the electron-hole vertices read here, upper and lower sign corresponds to Γ h1e1 and Γ h2e2 respectively. We present the results of the projections for the bare interactions in Table III. As one clearly sees, the largest repulsive interaction in the Cooper-pairing channel occurs in the A 1g −channel. However, in contrast to the case of smaller n, here the Fermi surface topology allows to have the A 1g -solution due to the appearance of the hole pockets near the Γ and M points of the Brillouin Zone. As a result, the repulsion between the electron and the hole pockets becomes effectively pair-building for the so-called s ± -wave symmetry due to the change of sign of the superconducting gap between electron and hole pockets. Furthermore, due to the different character of the bands, the interaction between electron and hole pockets will be enhanced even for the moderate renormalization within RPA, as seen from Fig.3(d). Therefore, this state wins over d x 2 −y 2 -wave for these doping concentrations. In particular, we find for the superconducting on the electron and hole pockets in the A 1g -channel where the signs of the constant gap on the electron and hole pockets are opposite. It is interesting to note that even for the moderate correlations within RPA the interband repulsion between the electron and the hole pockets becomes larger than the intraband or the interband within pockets of the same character, i.e. h 1 h 2 or e 1 e 2 .
As a result the gaps on the hole pockets can be considered as constants to a good approximation. On the electron pockets the gaps still acquire some significant angle dependence with∆ e < ∆ e . Analysing further the behaviour of the interactions and solving the linear version of the BCS gap equation we found that this part of the doping phase diagram is dominated by the s +− −wave symmetry with sign changing gap on the electron and hole pockets. Although such a filling factor was not achieved for BiS 2 -systems, yet it would be interesting to investigate its potential realization in these and similar systems.

VI. CONCLUSIONS
To conclude, we develop the effective low-energy Hamiltonian of BiS 2 lyered superconductors based on the fully relativistic ab-initio band structure calculations, including the spin-orbit coupling, relevant for Bi-based systems. The model consists of the tight-binding Hamiltonian, based on all three p-orbitals, p x , p y , and p z . The spin orbit coupling introduces weak spin anisotropy, which does not modify the spin structure of the Cooperpairing. In particular, due to relatively weak contribution of the p z orbitals to the Fermi surface, the Cooper-pairing can be still regarded as even parity spin singlet state. Despite of the weakness of correlations for the p−electrons we find that the purely repulsive interaction still yields unconventional superconductivity as soon as the interand and intraband interactions can be separated in the momentum space. In particular, for n ∼ 0.3 the superconducting gap possesses the global d x 2 −y 2 −wave symmetry yet it acquires strong angular dependence forming accidental nodes or deep minima on the electron pockets, located at X and Y points of the BZ. This generally agrees with recent ARPES experiments 28 and indicates that the mechanism of the Cooper-pairing is still repulsive in nature. Despite the near-nodal or nodal behavior and strong angular dependence of the gaps on the Fermi surface, the global symmetry remains d x 2 −y 2 -wave. This is reflected by the opposite signs of the cos 4φ harmonics on the two electron pockets, respectively. This interesting predication needs to be further verified experimentally. Furthermore, with increasing n we find another doping region where the interaction has a well-defined momentum structure with a clear separation of the inetrband and intraband scattering. The topology of the Fermi surface is then similar to the case of iron-based superconductors with electron and hole pockets. In this case the nesting between electron and hole bands promotes the nodeless A 1g -symmetry representation to be the dominant solution in this case even for the weakly correlated case.