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

Recent discovery of BiS₂-layered superconductors with T_c up to 10 K has attracted significant attention due to striking similarities of their crystal structures with cuprate and iron-based high-T_c superconductors [1–11]. According to the density functional theory (DFT) calculations, the parent compounds of the BiS₂-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–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.5F_0.5BiS₂ 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 on possible pairing symmetries in these new superconductors arising due to repulsive interactions [21–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.71F_0.29BiS₂ 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 a realistic minimal electronic model for recently discovered BiS₂ superconductors including the spin–orbit (SO) coupling, relevant for Bi-based systems. We show that the tight-binding low-energy model necessarily includes all three p-orbitals, p_x, p_y, and p_z. Using leading 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 correlations, 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 representation 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 analyze the form of these harmonics and show that they are consistent with those found in recent experimental ARPES data [28]. We further investigate the Fermi surface topology for BiS₂ 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.

The fully relativistic band structure calculation of the LaOBiS₂ system was performed within DFT using the Perdew–Burke–Ernzerhof exchange-correlation functional as implemented in WIEN2K program [31]. The muffin-tin radius of each atom R_MT was chosen such that its product with the maximum modulus of reciprocal vectors K_max become R_MTK_max = 7.0. The lattice parameters were taken from [32] and the corresponding Brillouin zone was sampled using a 12 × 12 × 4 k-mesh. From this band calculation, as shown in figure 1, we constructed a 12-band tight-binding (TB) model using the maximally localized Wannier functions [33–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, expressed as

$$\begin{eqnarray}&&{\hat{{ \mathcal H }}}_{0}=\displaystyle \sum _{{ij}}\displaystyle \sum _{\nu \nu ^{\prime} \sigma }{t}_{{ij};\nu \nu ^{\prime} }{c}_{i\nu \sigma }^{\dagger }{c}_{j\nu ^{\prime} \sigma }^{}+\displaystyle \sum _{i\nu \sigma }{\varepsilon }_{\nu }{c}_{i\nu \sigma }^{\dagger }{c}_{i\nu \sigma }^{},\end{eqnarray} \tag{ 1 }$$

plus SO coupling by ignoring the inter-layer hopping between the BiS₂ layers. Here, ${t}_{{ij};\nu \nu ^{\prime} }=t({x}_{i}-{x}_{j},{y}_{i}-{y}_{j};\nu \nu ^{\prime} )$ are the nearest and next-nearest neighbor hopping parameters obtained from projecting the results of ab initio calculations. Their values are shown in table 1 where ν, and ν' refer to the orbital indices. Thus, in the momentum space the Hamiltonian can be written as

$$\begin{eqnarray}&&{\hat{{ \mathcal H }}}_{0}=\displaystyle \sum _{{\bf{k}}\nu \nu ^{\prime} \sigma }{T}_{\nu \nu ^{\prime} }({\bf{k}}){c}_{{\bf{k}}\nu \sigma }^{\dagger }{c}_{{\bf{k}}\nu ^{\prime} \sigma }^{},\end{eqnarray} \tag{ 2 }$$

with the hopping integrals

$$\begin{eqnarray*}{T}_{\nu \nu } & = & {t}_{0,0}^{\nu }+2{t}_{1,0}^{\nu }\cos {k}_{x}+2{t}_{0,1}^{\nu }\cos {k}_{y}+2{t}_{1,1}^{\nu }\cos ({k}_{x}+{k}_{y})+2{t}_{1,-1}^{\nu }\cos ({k}_{x}-{k}_{y}),\\ {T}_{\nu \nu ^{\prime} } & = & {T}_{\nu ^{\prime} \nu }={t}_{0,0}^{\nu \nu ^{\prime} }+2{t}_{1,0}^{\nu \nu ^{\prime} }\cos {k}_{x}+2{t}_{0,1}^{\nu \nu ^{\prime} }\cos {k}_{y}+2{t}_{1,1}^{\nu \nu ^{\prime} }\cos ({k}_{x}+{k}_{y})+2{t}_{1,-1}^{\nu \nu ^{\prime} }\cos ({k}_{x}-{k}_{y}).\end{eqnarray*}$$

Here ${t}_{x,y}^{\nu }$ and ${t}_{x,y}^{\nu \nu ^{\prime} }$ denote intra- and interorbital electron hopping, as determined in table 1.

Zoom In Zoom Out Reset image size

On top of that, the on-site SO coupling Hamiltonian, ${H}_{{\rm{SO}}}=\lambda \,{\bf{L}}\cdot {\bf{S}}$, written for the bismuth p-orbitals have the following form:

$$\begin{eqnarray}{{ \mathcal H }}_{{SO}}=\displaystyle \frac{\lambda }{2}\left(\begin{array}{cccccc}0 & -{\rm{i}} & 0 & 0 & 0 & 1\\ {\rm{i}} & 0 & 0 & 0 & 0 & -{\rm{i}}\\ 0 & 0 & 0 & -1 & {\rm{i}} & 0\\ 0 & 0 & -1 & 0 & {\rm{i}} & 0\\ 0 & 0 & -{\rm{i}} & -{\rm{i}} & 0 & 0\\ 1 & {\rm{i}} & 0 & 0 & 0 & 0\end{array}\right),\end{eqnarray} \tag{ 3 }$$

hereafter, we set λ = 0.874 eV 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 figure 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 SO coupling does require to involve the p_z-orbital into consideration. This is valid even for the small fillings where there always exists an admixture of the p_z orbitals to the lowest band due to SO coupling.

Zoom In Zoom Out Reset image size

For the small fillings the Fermi surface of the BiS₂-layer 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 [12]. Due to SO 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 regime 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 figure 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 interband electron–electron scattering and favor non-phononic mechanism of superconductivity.

As mentioned in the introduction, the origin of superconductivity in BiS₂ layers is still debated. Nevertheless, given the similarity of the BiS₂-systems and the iron-based superconductors regarding crystal structure and Fermi surface topology, it is tempting to apply the concept of Cooper-pairing from repulsion to these systems as well, given much weaker correlations strength in Bi-based 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₂ 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 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 the bare spin susceptibility is defined by

$$\begin{eqnarray}&&{({\chi }^{{ss}^{\prime} })}_{{\nu }_{1}{\nu }_{2}}^{{\nu }_{3}{\nu }_{4}}({\boldsymbol{q}})=\langle T{\hat{S}}_{{\nu }_{1}{\nu }_{2}}^{s}({\boldsymbol{q}},\tau ){\hat{S}}_{{\nu }_{3}{\nu }_{4}}^{s^{\prime} }(-{\boldsymbol{q}},0)\rangle ,\end{eqnarray} \tag{ 4 }$$

where ν_i runs over orbitals, and ${\hat{S}}^{s}$ are the spin operators given by

$$\begin{eqnarray}&&{\hat{S}}_{\mu \nu }^{s}({\boldsymbol{q}},\tau )=\displaystyle \frac{1}{2}\displaystyle \sum _{{\boldsymbol{k}}\sigma \sigma ^{\prime} }{c}_{{\boldsymbol{k}}\mu \sigma }^{\dagger }(\tau ){\hat{\sigma }}_{\sigma \sigma ^{\prime} }^{s}{c}_{{\boldsymbol{k}}+{\boldsymbol{q}}\nu \sigma ^{\prime} }^{}(\tau ),\end{eqnarray} \tag{ 5 }$$

where ${\hat{\sigma }}^{s}$ are Pauli matrices. Therefore, by applying Wick's theorem and neglecting contractions leading to ${\boldsymbol{q}}={\bf{0}}$, one can find that the different components of the bare spin susceptibility, ${\chi }_{0}^{{ss}^{\prime} }$, in the frequency domain is

$$\begin{eqnarray}&&{({\chi }_{0}^{{ss}^{\prime} })}_{{\nu }_{3}{\nu }_{4}}^{{\nu }_{2}{\nu }_{1}}({\boldsymbol{q}},{\rm{i}}{\rm{\Omega }})=\displaystyle \frac{-T}{4{N}}\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{{\boldsymbol{k}},{\rm{i}}{\omega }_{{n}}}{\gamma \delta \gamma ^{\prime} \delta ^{\prime} }}{\hat{\sigma }}_{\gamma \delta }^{s}{\hat{\sigma }}_{\gamma ^{\prime} \delta ^{\prime} }^{s^{\prime} }{G}_{{\nu }_{4}{\nu }_{1},\delta ^{\prime} \gamma }^{{\boldsymbol{k}},{\rm{i}}{\omega }_{{n}}}{G}_{{\nu }_{2}{\nu }_{3},\delta \gamma ^{\prime} }^{{\boldsymbol{k}}+{\boldsymbol{q}},{\rm{i}}{\omega }_{{n}}^{{\prime} }},\end{eqnarray} \tag{ 6 }$$

with the Green's functions introduced by

$$\begin{eqnarray}&&{G}_{\nu \nu ^{\prime} ,\sigma \sigma ^{\prime} }^{{\boldsymbol{k}},{\rm{i}}{\omega }_{{n}}}=-{\int }_{0}^{\beta }{\rm{d}}\tau {{\rm{e}}}^{{\rm{i}}{\omega }_{{n}}\tau }\langle {T}_{\tau }\,{c}_{{\boldsymbol{k}}\nu \sigma }(\tau ){c}_{{\boldsymbol{k}}\nu ^{\prime} \sigma ^{\prime} }^{\dagger }(0)\rangle ,\end{eqnarray} \tag{ 7 }$$

here we set ${\rm{i}}\omega ^{\prime} ={\rm{i}}{\omega }_{{n}}+{\rm{i}}{\rm{\Omega }}$. Since the Green's functions are not diagonal in the orbital basis it is convenient to move to the band representation, and find

$$\begin{eqnarray}&&{G}_{\nu \nu ^{\prime} ,\sigma \sigma ^{\prime} }^{{\boldsymbol{k}},\tau }=\displaystyle \sum _{\mu }{\zeta }_{{\boldsymbol{k}}\nu \sigma }^{\mu }{\zeta }_{{\boldsymbol{k}}\nu ^{\prime} \sigma ^{\prime} }^{\mu * }{G}_{\mu }^{{\boldsymbol{k}},\tau },\end{eqnarray} \tag{ 8 }$$

where μ being the band index, and the matrix elements ${\zeta }_{{\boldsymbol{k}}\nu \sigma }^{\mu }$, are connecting the band to the orbital basis: ${c}_{{\boldsymbol{k}}\nu \sigma }(\tau )={\sum }_{\mu }{\zeta }_{{\boldsymbol{k}}\nu \sigma }^{\mu }{b}_{k\mu }(\tau )$. By performing the Matsubara frequency sum over ω_n, yields

$$\begin{eqnarray}&&{({\chi }_{0}^{{ss}^{\prime} })}_{{\nu }_{3}{\nu }_{4}}^{{\nu }_{2}{\nu }_{1}}({\boldsymbol{q}},\omega )=-\displaystyle \frac{T}{4N}\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{{\boldsymbol{k}},\mu \mu ^{\prime} }{\gamma ,\delta ,\gamma ^{\prime} ,\delta ^{\prime} }}{\hat{\sigma }}_{\gamma \delta }^{s}{\hat{\sigma }}_{\gamma ^{\prime} \delta ^{\prime} }^{s^{\prime} }({\zeta }_{{\boldsymbol{k}}{\nu }_{4}\delta ^{\prime} }^{\mu }{\zeta }_{{\boldsymbol{k}}{\nu }_{1}\gamma }^{\mu * }{\zeta }_{{\boldsymbol{k}}+{\boldsymbol{q}}{\nu }_{2}\delta }^{\mu ^{\prime} }{\zeta }_{{\boldsymbol{k}}+{\boldsymbol{q}}{\nu }_{3}\gamma ^{\prime} }^{\mu ^{\prime} * })\displaystyle \frac{f({E}_{{\boldsymbol{k}}+{\boldsymbol{q}},\mu ^{\prime} })-f({E}_{{\boldsymbol{k}},\mu })}{\omega +{\rm{i}}{0}^{+}+{{E}}_{{\boldsymbol{k}}+{\boldsymbol{q}},\mu ^{\prime} }-{{E}}_{{\boldsymbol{k}},\mu }},\end{eqnarray} \tag{ 9 }$$

where f(...) are Fermi functions. The straightforward calculations based on the Pauli matrices relation:

$$\begin{eqnarray*}{\hat{\sigma }}_{\gamma \delta }^{+}{\hat{\sigma }}_{\gamma ^{\prime} \delta ^{\prime} }^{-}=\left\{\begin{array}{lcc}4 & \mathrm{if} & \gamma =\delta ^{\prime} \ne \gamma ^{\prime} =\delta \\ 0 & & \mathrm{otherwise}\end{array}\right.\end{eqnarray*}$$

and

$$\begin{eqnarray*}{\hat{\sigma }}_{\gamma \delta }^{z}{\hat{\sigma }}_{\gamma ^{\prime} \delta ^{\prime} }^{z}=\left\{\begin{array}{lcc}1 & \mathrm{if} & \gamma =\delta =\gamma ^{\prime} =\delta ^{\prime} \\ -1 & \mathrm{if} & \gamma =\delta =-\gamma ^{\prime} =-\delta ^{\prime} \\ 0 & & \mathrm{otherwise}\end{array}\right.\end{eqnarray*}$$

defined the transverse susceptibility as follows

$$\begin{eqnarray}&&{({\chi }_{0}^{\pm })}_{{\nu }_{3}{\nu }_{4}}^{{\nu }_{2}{\nu }_{1}}({\boldsymbol{q}},\omega )=\displaystyle \frac{-T}{4N}\displaystyle \sum _{{\boldsymbol{k}}\mu \mu ^{\prime} }({\zeta }_{{\boldsymbol{k}}{\nu }_{4}\uparrow }^{\mu }{\zeta }_{{\boldsymbol{k}}{\nu }_{1}\uparrow }^{\mu * }{\zeta }_{{\boldsymbol{k}}+{\boldsymbol{q}}{\nu }_{2}\downarrow }^{\mu ^{\prime} }{\zeta }_{{\boldsymbol{k}}+{\boldsymbol{q}}{\nu }_{3}\downarrow }^{\mu ^{\prime} * })\displaystyle \frac{f({E}_{{\boldsymbol{k}}+{\boldsymbol{q}},\mu ^{\prime} })-f({E}_{{\boldsymbol{k}},\mu })}{\omega +{\rm{i}}{0}^{+}+{{E}}_{{\boldsymbol{k}}+{\boldsymbol{q}},\mu ^{\prime} }-{{E}}_{{\boldsymbol{k}},\mu }},\end{eqnarray} \tag{ 10 }$$

and by considering $\tilde{\sigma }=-\sigma$, we find

$$\begin{eqnarray}&&{({\chi }^{{zz}})}_{{\nu }_{3}{\nu }_{4}}^{{\nu }_{2}{\nu }_{1}}({\boldsymbol{q}},\omega )=\displaystyle \frac{-T}{4N}\displaystyle \sum _{{\boldsymbol{k}}\mu \mu ^{\prime} \sigma }({\zeta }_{{\boldsymbol{k}}{\nu }_{4}\sigma }^{\mu }{\zeta }_{{\boldsymbol{k}}{\nu }_{1}\sigma }^{\mu * }{\zeta }_{{\boldsymbol{k}}+{\boldsymbol{q}}{\nu }_{2}\sigma }^{\mu ^{\prime} }{\zeta }_{{\boldsymbol{k}}+{\boldsymbol{q}}{\nu }_{3}\sigma }^{\mu ^{\prime} * }-{\zeta }_{{\boldsymbol{k}}{\nu }_{4}\sigma }^{\mu }{\zeta }_{{\boldsymbol{k}}{\nu }_{1}\tilde{\sigma }}^{\mu * }{\zeta }_{{\boldsymbol{k}}+{\boldsymbol{q}}{\nu }_{2}\tilde{\sigma }}^{\mu ^{\prime} }{\zeta }_{{\boldsymbol{k}}+{\boldsymbol{q}}{\nu }_{3}\sigma }^{\mu ^{\prime} * })\displaystyle \frac{f({E}_{{\boldsymbol{k}}+{\boldsymbol{q}},\mu ^{\prime} })-f({E}_{{\boldsymbol{k}},\mu })}{\omega +{\rm{i}}{0}^{+}+{{E}}_{{\boldsymbol{k}}+{\boldsymbol{q}},\mu ^{\prime} }-{{E}}_{{\boldsymbol{k}},\mu }},\end{eqnarray} \tag{ 11 }$$

for the longitudinal susceptibility.

The results of the calculations for the two representative concentrations of the model are present in figures 3(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. ${\chi }^{+-}\gt {\chi }^{{zz}}$, which appears to be an overall prefactor. Therefore we discuss only the behavior of ${\chi }^{+-}$ 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 RPA approach. In this respect, the real space picture of the interacting Hamiltonian is given by [23]

$$\begin{eqnarray}&&{{ \mathcal H }}_{{\rm{int}}}=\displaystyle \sum _{{\boldsymbol{i}},\nu }{{Un}}_{{\boldsymbol{i}}\nu \uparrow }{n}_{{\boldsymbol{i}}\nu \downarrow }+\displaystyle \sum _{{\boldsymbol{i}},\nu \ne \nu ^{\prime} }\left[\left(U^{\prime} -\displaystyle \frac{J}{2}\right){n}_{{\boldsymbol{i}}\nu }{n}_{{\boldsymbol{i}}\nu ^{\prime} }-2J{{\boldsymbol{S}}}_{{\boldsymbol{i}}\nu }\cdot {{\boldsymbol{S}}}_{{\boldsymbol{i}}\nu ^{\prime} }+J^{\prime} ({c}_{{\boldsymbol{i}}\nu \uparrow }^{\dagger }{c}_{{\boldsymbol{i}}\nu \downarrow }^{\dagger }{c}_{{\boldsymbol{i}}\nu ^{\prime} \downarrow }{c}_{{\boldsymbol{i}}\nu ^{\prime} \uparrow }+{\rm{h}}{\rm{.}}{\rm{c}}.)\right],\end{eqnarray} \tag{ 12 }$$

where ${n}_{{\boldsymbol{i}}\nu }={\sum }_{\sigma }{n}_{{\boldsymbol{i}}\nu \sigma }={\sum }_{\sigma }{c}_{{\boldsymbol{i}}\nu \sigma }^{\dagger }{c}_{{\boldsymbol{i}}\nu \sigma }^{}$, and ${{\boldsymbol{S}}}_{{\boldsymbol{i}}\nu }={\sum }_{\sigma ,\sigma ^{\prime} }{c}_{{\boldsymbol{i}}\nu \sigma }^{\dagger }{\hat{\sigma }}_{\sigma \sigma ^{\prime} }^{}{c}_{{\boldsymbol{i}}\nu \sigma ^{\prime} }^{}$. Therefore, in the band representation the total Hamiltonian reads

$$\begin{eqnarray}&&\hat{{ \mathcal H }}=\displaystyle \sum _{\mu ,{\boldsymbol{k}}}{\epsilon }_{\mu }({\boldsymbol{k}}){b}_{\mu {\boldsymbol{k}}}^{\dagger }{b}_{\mu {\boldsymbol{k}}}^{}+\displaystyle \sum _{\mu \mu ^{\prime} ,{\boldsymbol{k}}{\boldsymbol{k}}^{\prime} }{{\rm{\Gamma }}}_{\mu \mu ^{\prime} }({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} ){b}_{\mu {\boldsymbol{k}}}^{\dagger }{b}_{\mu ,-{\boldsymbol{k}}}^{\dagger }{b}_{\mu ^{\prime} {\boldsymbol{k}}^{\prime} }{b}_{\mu ^{\prime} ,-{\boldsymbol{k}}^{\prime} },\end{eqnarray} \tag{ 13 }$$

where the quartic terms describe the scattering of pairs $({\boldsymbol{k}}\uparrow ,\,-{\boldsymbol{k}}\downarrow )$ on the pocket μ to $({\boldsymbol{k}}^{\prime} \uparrow ,\,-{\boldsymbol{k}}^{\prime} \downarrow )$ on the pocket μ', with subsequent scattering amplitude:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mu \mu ^{\prime} }({{\boldsymbol{k}}}_{1},{{\boldsymbol{k}}}_{2})=\displaystyle \sum _{{ptsq}}{\zeta }_{-{{\boldsymbol{k}}}_{1}\mu }^{p* }{\zeta }_{{{\boldsymbol{k}}}_{1}\mu ^{\prime} }^{t* }[{{\rm{\Gamma }}}_{{tp}}^{{sq}}({{\boldsymbol{k}}}_{1},{{\boldsymbol{k}}}_{2})]{\zeta }_{{{\boldsymbol{k}}}_{2}\mu }^{s}{\zeta }_{-{{\boldsymbol{k}}}_{2}\mu ^{\prime} }^{q},\end{eqnarray} \tag{ 14 }$$

in which, the vertex $[{{\rm{\Gamma }}}_{{tp}}^{{sq}}({{\boldsymbol{k}}}_{1},{{\boldsymbol{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, ${\zeta }_{{\boldsymbol{k}}\mu }^{p}={\zeta }_{{\boldsymbol{k}}p}^{\mu * }$ are the matrix elements connecting the orbital and band basis, and they represent the orbital contributions to the Fermi surface, see figures 2(c) and (d).

The final RPA result of the susceptibility can be described in the form of Dyson equations as follows:

$$\begin{eqnarray}&&{({\chi }_{{\rm{RPA}}})}_{s,t}^{p,q}={({\chi }_{0})}_{s,t}^{p,q}+{({\chi }_{{\rm{RPA}}})}_{u,v}^{p,q}{({U}_{\mu })}_{w,z}^{u,v}{({\chi }_{0})}_{s,t}^{w,z},\end{eqnarray} \tag{ 15 }$$

where $({U}_{\mu }{)}_{w,z}^{u,v}$ 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

$$\begin{eqnarray}&&{({\chi }_{{\rm{RPA}}}^{\pm })}_{s,t}^{p,q}={({\chi }_{0}^{\pm })}_{s,t}^{p,q}+{({\chi }_{{\rm{RPA}}}^{\pm })}_{u,v}^{p,q}{({U}_{\mu }^{\pm })}_{w,z}^{u,v}{({\chi }_{0}^{\pm })}_{s,t}^{w,z},\end{eqnarray} \tag{ 16 }$$

where the non-zero components of U_μ are given by

$$\begin{eqnarray*}{({U}_{\mu }^{\pm })}_{a,a}^{a,a} & = & \tilde{U};\ \ \ \ \ \ {({U}_{\mu }^{\pm })}_{b,b}^{a,a}=\tilde{J};\\ {({U}_{\mu }^{\pm })}_{a,b}^{a,b} & = & \tilde{V};\ \ \ \ \ \ {({U}_{\mu }^{\pm })}_{b,a}^{a,b}=\tilde{J}^{\prime} .\end{eqnarray*}$$

however, for the longitudinal susceptibility we find

$$\begin{eqnarray}&&{({\chi }_{{\rm{RPA}}}^{{zz}})}_{s,t}^{p,q}=\displaystyle \frac{1}{4}\left[\displaystyle \sum _{\sigma }{({\chi }_{{\rm{RPA}}}^{\sigma \sigma })}_{s,t}^{p,q}-\displaystyle \sum _{\sigma \ne \sigma ^{\prime} }{({\chi }_{{\rm{RPA}}}^{\sigma \sigma ^{\prime} })}_{s,t}^{p,q}\right],\end{eqnarray} \tag{ 17 }$$

where

$$\begin{eqnarray}&&{({\chi }_{{\rm{RPA}}}^{\sigma \sigma ^{\prime} })}_{s,t}^{p,q}={({\chi }_{0}^{\sigma \sigma ^{\prime} })}_{s,t}^{p,q}+{({\chi }_{{\rm{RPA}}}^{\sigma \sigma ^{\prime} })}_{u,v}^{p,q}{({U}_{s}^{\sigma ^{\prime} \sigma })}_{w,z}^{u,v}{({\chi }_{0}^{\sigma \sigma ^{\prime} })}_{s,t}^{w,z},\end{eqnarray} \tag{ 18 }$$

with

$$\begin{eqnarray*}{({U}_{\mu }^{\sigma \sigma ^{\prime} })}_{a,a}^{a,a} & = & \tilde{U};\ \ \ \ \ \ {({U}_{\mu }^{\sigma \sigma ^{\prime} })}_{b,b}^{a,a}=\tilde{V};\\ {({U}_{\mu }^{\sigma \sigma ^{\prime} })}_{a,b}^{a,b} & = & \tilde{J};\ \ \ \ \ \ {({U}_{\mu }^{\sigma \sigma ^{\prime} })}_{b,a}^{a,b}=\tilde{J}^{\prime} ,\end{eqnarray*}$$

and finally

$$\begin{eqnarray}&&{({\chi }_{{\rm{RPA}}}^{\sigma \sigma })}_{s,t}^{p,q}={({\chi }_{0}^{\sigma \sigma })}_{s,t}^{p,q}+{({\chi }_{{\rm{RPA}}}^{\sigma \sigma })}_{u,v}^{p,q}{({U}_{\mu }^{\sigma \sigma })}_{w,z}^{u,v}{({\chi }_{0}^{\sigma \sigma })}_{s,t}^{w,z},\end{eqnarray} \tag{ 19 }$$

with

$$\begin{eqnarray*}{({U}_{\mu }^{\sigma \sigma })}_{a,a}^{a,a} & = & 0;\ \ \ \ \ \ {({U}_{\mu }^{\sigma \sigma })}_{b,b}^{a,a}=\tilde{V}+\tilde{J};\\ {({U}_{\mu }^{\sigma \sigma })}_{a,b}^{a,b} & = & \tilde{J}+\tilde{V};\ \ \ \ \ \ {({U}_{\mu }^{\sigma \sigma })}_{b,a}^{a,b}=0.\end{eqnarray*}$$

Despite the fact that the correlations should be moderate in BiS₂ 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}} \mbox{-}$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 LAHA, 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.

The LAHA approximation was developed for the iron-based 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₂ superconductors. The basic idea of LAHA is that the ${{\rm{\Gamma }}}_{\mu \mu ^{\prime} }$, defined in equation (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 ${\boldsymbol{k}}$ and ${\boldsymbol{k}}^{\prime}$ [29, 36, 37]. In particular, one can decompose the ${{\rm{\Gamma }}}_{\mu \mu ^{\prime} }({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )$ into representations of the tetragonal space group A_1g (s-wave) and B_1g (d-wave), i.e.,

$$\begin{eqnarray*}&&{{\rm{\Gamma }}}^{({A}_{1g})}({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )={{\rm{\Gamma }}}^{s}({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )=\displaystyle \sum _{m,n}{A}_{{mn}}^{s}{{\rm{\Psi }}}_{m}^{s}({\boldsymbol{k}}){{\rm{\Psi }}}_{n}^{s}({\boldsymbol{k}}^{\prime} ),\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{{\rm{\Gamma }}}^{({B}_{1g})}({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )={{\rm{\Gamma }}}^{d}({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )=\displaystyle \sum _{m,n}{A}_{{mn}}^{d}{{\rm{\Psi }}}_{m}^{d}({\boldsymbol{k}}){{\rm{\Psi }}}_{n}^{d}({\boldsymbol{k}}^{\prime} ),\end{eqnarray*}$$

with Ψ^s and Ψ^d being the basis functions of the A_1g and B_1g representations, respectively. Finally, the ${{\rm{\Gamma }}}^{({A}_{1g})}$ and ${{\rm{\Gamma }}}^{({B}_{1g})}$ 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₂ 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 pseudospin-singlet Cooper-pairing and that the Cooper-pairs can be further characterized as even-parity spin singlet solutions. This is due to the fact that the SO 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 pseudospin 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.

5.1. Low doping limit: n = 0.2

The Fermi surface of the system for small doping limit consists of two electron pockets at (π, 0), e₁, and (0, π), e₂. In contrast to the Fe-based superconductors, the topology of these pockets is rather non-elliptical (figure 2), implying that higher harmonics (up to cos 4ϕ) need to be considered to correctly describe the scattering processes:

$$\begin{eqnarray*}{{\rm{\Gamma }}}_{{e}_{1}{e}_{2}}^{s} & = & {U}_{{ee}}[1+2{\alpha }_{{ee}}(\cos 2{\phi }_{{k}_{e1}}-\cos 2{\phi }_{{k}_{e2}^{{\prime} }})-4{\beta }_{{ee}}\cos 2{\phi }_{{k}_{e1}}\cos 2{\phi }_{{k}_{e2}^{{\prime} }}+2{\gamma }_{{ee}}(\cos 4{\phi }_{{k}_{e1}}+\cos 4{\phi }_{{k}_{e2}^{{\prime} }})\\ & & +2{\delta }_{{ee}}(\cos 2{\phi }_{{k}_{e1}}\cos 4{\phi }_{{k}_{e2}^{{\prime} }}-\cos 2{\phi }_{{k}_{e2}^{{\prime} }}\cos 4{\phi }_{{k}_{e1}})+4{\eta }_{{ee}}\cos 4{\phi }_{{k}_{e1}}\cos 4{\phi }_{{k}_{e2}^{{\prime} }}],\\ {{\rm{\Gamma }}}_{{e}_{1}{e}_{2}}^{d} & = & {\tilde{U}}_{{ee}}[-1-2{\tilde{\alpha }}_{{ee}}(\cos 2{\phi }_{{k}_{e1}}-\cos 2{\phi }_{{k}_{e2}^{{\prime} }})+4{\tilde{\beta }}_{{ee}}\cos 2{\phi }_{{k}_{e1}}\cos 2{\phi }_{{k}_{e2}^{{\prime} }}-2{\tilde{\gamma }}_{{ee}}(\cos 4{\phi }_{{k}_{e1}}+\cos 4{\phi }_{{k}_{e2}^{{\prime} }})\\ & & \mp 2{\tilde{\delta }}_{{ee}}(\cos 2{\phi }_{{k}_{e1}}\cos 4{\phi }_{{k}_{e2}^{{\prime} }}-\cos 2{\phi }_{{k}_{e1}^{{\prime} }}\cos 4{\phi }_{{k}_{e2}})-4{\tilde{\eta }}_{{ee}}\cos 4{\phi }_{{k}_{e1}}\cos 4{\phi }_{{k}_{e2}^{{\prime} }}],\\ {{\rm{\Gamma }}}_{{e}_{i}{e}_{i}}^{s} & = & {U}_{{ee}}[1\pm 2{\alpha }_{{ee}}(\cos 2{\phi }_{{k}_{{ei}}}+\cos 2{\phi }_{{k}_{{ei}}^{{\prime} }})+4{\beta }_{{ee}}\cos 2{\phi }_{{k}_{{ei}}}\cos 2{\phi }_{{k}_{e1}^{{\prime} }}+2{\gamma }_{{ee}}(\cos 4{\phi }_{{k}_{{ei}}}+\cos 4{\phi }_{{k}_{{ei}}^{{\prime} }})\\ & & \pm 2{\delta }_{{ee}}(\cos 2{\phi }_{{k}_{{ei}}}\cos 4{\phi }_{{k}_{{ei}}^{{\prime} }}+\cos 2{\phi }_{{k}_{{ei}}^{{\prime} }}\cos 4{\phi }_{{k}_{{ei}}})+4{\eta }_{{ee}}\cos 4{\phi }_{{k}_{{ei}}}\cos 4{\phi }_{{k}_{{ei}}^{{\prime} }}],\\ {{\rm{\Gamma }}}_{{e}_{i}{e}_{i}}^{d} & = & {\tilde{U}}_{{ee}}[1\pm 2{\tilde{\alpha }}_{{ee}}(\cos 2{\phi }_{{k}_{{ei}}}+\cos 2{\phi }_{{k}_{{ei}}^{{\prime} }})+4{\tilde{\beta }}_{{ee}}\cos 2{\phi }_{{k}_{{ei}}}\cos 2{\phi }_{{k}_{{ei}}^{{\prime} }}+2{\tilde{\gamma }}_{{ee}}(\cos 4{\phi }_{{k}_{{ei}}}+\cos 4{\phi }_{{k}_{{ei}}^{{\prime} }})\\ & & \pm 2{\tilde{\delta }}_{{ee}}(\cos 2{\phi }_{{k}_{{e}_{i}}}\cos 4{\phi }_{{k}_{{ei}}^{{\prime} }}+\cos 2{\phi }_{{k}_{{ei}}^{{\prime} }}\cos 4{\phi }_{{k}_{{ei}}})+4{\tilde{\eta }}_{{ee}}\cos 4{\phi }_{{k}_{{ei}}}\cos 4{\phi }_{{k}_{{ei}}^{{\prime} }}].\end{eqnarray*}$$

Here the upper (lower) sign corresponds to pocket 1 (2), ${\phi }_{{k}_{{ei}}}$ (${\phi }_{{k}_{{hi}}}$ ) represents the angles along the ith 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 2.

Table 2.  Parameters of the LAHA projection of the electronic interactions projected on the extended s-wave and ${d}_{{x}^{2}-{y}^{2}}$-wave symmetry representations for n = 0.2 and with U = 0.6, V = 0.36 and J = J' = 0.12 (all in eV).

Uee	αee	βee	γee	δee	ηee
0.393	−1.463 × 10−4	1.438 × 10−5	8.866 × 10−3	1.009 × 10−4	6.172 × 10−4
${\tilde{U}}_{{ee}}$	${\tilde{\alpha }}_{{ee}}$	${\tilde{\beta }}_{{ee}}$	${\tilde{\gamma }}_{{ee}}$	${\tilde{\delta }}_{{ee}}$	${\tilde{\eta }}_{{ee}}$
0.041	0.0764	5.912 × 10−3	0.0516	7.885 × 10−3	2.660 × 10−3

Observe that the constant part of the interaction in the A_1g-pairing channel is larger than that 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 electron like. An inclusion of the RPA corrections in the particle–hole channel due to interband nesting helps in achieving ${\tilde{U}}_{e1e2}\gt {\tilde{U}}_{e1e1},{\tilde{U}}_{e2e2}$. Nevertheless, all interactions are still very similar to each other as the spin fluctuations are not strongly enhanced, see figure 3(c). As a result the solution for the superconducting gap acquires a strong angular dependence in the B_1g-channel in the form

$$\begin{eqnarray}{{\rm{\Delta }}}_{{e}_{1}}(\phi ) & = & {{\rm{\Delta }}}_{e}+{\bar{{\rm{\Delta }}}}_{e}\cos 2\phi +{\tilde{{\rm{\Delta }}}}_{e}\cos 4\phi +...\\ {{\rm{\Delta }}}_{{e}_{2}}(\phi ) & = & -{{\rm{\Delta }}}_{e}+{\bar{{\rm{\Delta }}}}_{e}\cos 2\phi -{\tilde{{\rm{\Delta }}}}_{e}\cos 4\phi \,+\,...,\end{eqnarray} \tag{ 20 }$$

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 behavior of the superconducting gap on each of the Fermi surface sheets, i.e. we find ${{\rm{\Delta }}}_{e}\sim \bar{{{\rm{\Delta }}}_{e}}\sim \tilde{{{\rm{\Delta }}}_{e}}$. This is similar to the experimental findings [28] (see equation (1) in [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 SO 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.

5.2. Intermediate doping limit: n = 1.4

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 ${{\rm{\Gamma }}}_{{h}_{i}{h}_{j}}$, and electron–hole ${{\rm{\Gamma }}}_{{e}_{i}{h}_{j}}$ 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:

$$\begin{eqnarray*}&&{{\rm{\Gamma }}}_{{e}_{i}{e}_{i}}^{s}={U}_{{ee}}[1\pm 2{\alpha }_{{ee}}(\cos 2{\phi }_{{k}_{{ei}}}+\cos 2{\phi }_{{k}_{{ei}}^{{\prime} }})+4{\beta }_{{ee}}\cos 2{\phi }_{{k}_{{ei}}}\cos 2{\phi }_{{k}_{{ei}}^{{\prime} }}],\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{{\rm{\Gamma }}}_{{e}_{i}{e}_{i}}^{d}={\tilde{U}}_{{ee}}[1\pm 2{\tilde{\alpha }}_{{ee}}(\cos 2{\phi }_{{k}_{{ei}}}+\cos 2{\phi }_{{k}_{{ei}}^{{\prime} }})+4{\tilde{\beta }}_{{ee}}\cos 2{\phi }_{{k}_{{ei}}}\cos 2{\phi }_{{k}_{{ei}}^{{\prime} }}]\end{eqnarray*}$$

and the hole pockets can be reasonably well approximated by circle i.e. by a constant term. Therefore we have,

$$\begin{eqnarray*}{{\rm{\Gamma }}}_{{h}_{i}{h}_{i}} & = & {{\rm{\Gamma }}}_{{h}_{i}{h}_{i}}^{s}+{{\rm{\Gamma }}}_{{h}_{i}{h}_{i}}^{d}={U}_{{h}_{i}{h}_{i}}+{\tilde{U}}_{{h}_{i}{h}_{i}}\cos 2{\phi }_{{k}_{{hi}}}\cos 2{\phi }_{{k}_{{hj}}^{{\prime} }},{{\rm{\Gamma }}}_{{h}_{1}{h}_{2}}={{\rm{\Gamma }}}_{{h}_{1}{h}_{2}}^{s}+{{\rm{\Gamma }}}_{{h}_{1}{h}_{2}}^{d}\\ & = & -{U}_{{h}_{1}{h}_{2}}-{\tilde{U}}_{{h}_{1}{h}_{2}}\cos 2{\phi }_{{k}_{{hi}}}\cos 2{\phi }_{{k}_{{hj}}^{{\prime} }}.\end{eqnarray*}$$

Finally, the electron–hole vertices read

$$\begin{eqnarray*}&&{{\rm{\Gamma }}}_{{h}_{i}{e}_{j}}={{\rm{\Gamma }}}_{{h}_{i}{e}_{j}}^{s}+{{\rm{\Gamma }}}_{{h}_{i}{e}_{j}}^{d}=\pm {U}_{{h}_{i}{e}_{j}}[1\pm {\alpha }_{{h}_{i}{e}_{j}}\cos (2{\phi }_{{k}_{{ej}}^{{\prime} }})]+{\tilde{U}}_{{h}_{i}{e}_{j}}\cos (2{\phi }_{{k}_{{hi}}})[1\pm {\tilde{\alpha }}_{{h}_{i}{e}_{j}}\cos (2{\phi }_{{k}_{{ej}}^{{\prime} }})].\end{eqnarray*}$$

Here, upper and lower sign corresponds to ${{\rm{\Gamma }}}_{{h}_{1}{e}_{1}}$ and ${{\rm{\Gamma }}}_{{h}_{2}{e}_{2}}$ respectively. We present the results of the projections for the bare interactions in table 3. 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 figure 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

$$\begin{eqnarray}{{\rm{\Delta }}}_{{h}_{1}}(\phi ) & = & {{\rm{\Delta }}}_{{h}_{1}},\\ {{\rm{\Delta }}}_{{h}_{2}}(\phi ) & = & {{\rm{\Delta }}}_{{h}_{2},}\\ {{\rm{\Delta }}}_{{h}_{3}}(\phi ) & = & {{\rm{\Delta }}}_{{h}_{3}},\\ {{\rm{\Delta }}}_{{e}_{1}}(\theta ) & = & {{\rm{\Delta }}}_{e}+{\bar{{\rm{\Delta }}}}_{e}\cos 2\theta ,\\ {{\rm{\Delta }}}_{{e}_{2}}(\theta ) & = & {{\rm{\Delta }}}_{e}-{\bar{{\rm{\Delta }}}}_{e}\cos 2\theta ,\end{eqnarray} \tag{ 21 }$$

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₁h₂ or e₁e₂. 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 ${\bar{{\rm{\Delta }}}}_{e}\lt {{\rm{\Delta }}}_{e}$. Analyzing further the behavior 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₂-systems, yet it would be interesting to investigate its potential realization in these and similar systems. The overall phase diagram as a function of doping is the shown in figure 4.

Table 3.  Parameters of the LAHA projection of the electronic interactions projected onto the extended s-wave and ${d}_{{x}^{2}-{y}^{2}}$-wave symmetry representations for n = 1.4 with U = 0.6, V = 0.36 and J = J' = 0.12 (all in eV).

Uee	αee	βee	Uh1h1	Uh2h2	Uh1h2	Uh1e	Uh2e	αh1e	αh2e
0.405	−7.909 × 10−4	3.156 × 10−5	0.303	0.321	0.309	0.237	0.238	9.928 × 10−4	−2.806 × 10−3
${\tilde{U}}_{{ee}}$	${\tilde{\alpha }}_{{ee}}$	${\tilde{\beta }}_{{ee}}$	${\tilde{U}}_{h1h1}$	${\tilde{U}}_{h2h2}$	${\tilde{U}}_{h1h2}$	${\tilde{U}}_{h1e}$	${\tilde{U}}_{h2e}$	${\tilde{\alpha }}_{h1e}$	${\tilde{\alpha }}_{h2e}$
0.0262	0.0894	8.144 × 10−3	0.0187	2.644 × 10−3	6.816 × 10−3	0.0221	8.143 × 10−3	0.179	0.182

To conclude, we develop the effective low-energy Hamiltonian of BiS₂ layered superconductors based on the fully relativistic ab initio band structure calculations, including the SO coupling, relevant for Bi-based systems. The model consists of a tight-binding Hamiltonian, based on all three p-orbitals, p_x, p_y, and p_z. The sizeable SO coupling introduces relatively weak spin anisotropy, which does not modify the spin structure of the Cooper-pairing. 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. This validates previous approach, based on the two-orbital model only [12]. Furthermore, this hints towards a different microscopic mechanism of superconducting in the other Bi-based compounds. Despite the weakness of correlations for the p-electrons we find that the purely repulsive interaction still yields unconventional superconductivity as soon as the interband and intraband interactions can be separated in the momentum space. In particular, for n ∼ 0.3, which corresponds to the experimentally achieved values, 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 and an application of the uniaxial strain could be also interesting to study. Furthermore, with increasing n we find another doping region where the interaction has a well-defined momentum structure with a clear separation of the interband 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.

We acknowledge helpful discussions with F Ahn, D Altenfeld, A Donkov and J Knolle. SC and AA wish to acknowledge the Korea Ministry of Education, Science and Technology, Gyeongsangbuk-Do and Pohang City for Independent Junior Research Groups at the Asia Pacific Center for Theoretical Physics. The work by AA was supported through NRF funded by MSIP of Korea (2015R1C1A1A01052411), and by Max Planck POSTECH/KOREA Research Initiative (No. 2011-0031558) programs through NRF funded by MSIP of Korea. IE acknowledges support by the Ministry of Education and Science of the Russian Federation in the framework of Increase Competitiveness Program of NUST MISiS (K2-2017-085). MSB gratefully acknowledges the support from the CREST, JST (No JPMJCR16F1).