Topological $p_{z}$-wave nodal-line superconductivity with flat surface bands in the AH_xCr₃As₃ (A=Na, K, Rb, Cs) superconductors

In recent years, the quasi-one-dimensional (Q1D) Cr-based superconductors have attracted great research interests. There are two Q1D Cr-based superconductors families, i.e. the A₂Cr₃As₃ (abbreviated as the 233) family (A=Na, K, Rb, Cs) [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29] and the ACr₃As₃ (abbreviated as the 133) family [33, 34, 35, 36, 37, 39, 30, 32, 31, 38], with both sharing the same Q1D alkaline-cations-separated [(Cr₃As₃)^2-]_∞ double-walled subnanotubes (DWSN) structure. The 233 family were synthesized earlier, and have shown unconventional pairing nature [2, 20, 22, 3, 4, 5, 21, 7, 40, 23, 25] with possible triplet pairing symmetry[2, 20, 22, 21, 40, 25, 41, 18]. The 133 family were synthesized slightly later, obtained by deintercalating half of the alkali metal atoms with ethanol[33, 34, 35, 36]. The primitive ACr₃As₃ have a spin glass ground state [33, 34, 39], which is consistent with the relaxation behavior of density functional theory (DFT) calculations[42, 39, 37, 38] and the muon spin-rotation ($\mu$SR) experiments [5]. However, when enough amount of hydrogen atoms are intercalated in the material, superconductivity (SC) emerges [43, 44, 45, 46], with highest $T_{c}$ around 8 K[30, 31]. Therefore, the correct chemical formula for“ACr₃As₃” should be AH_xCr₃As₃. It is the difference in H content that leads to the difference in ground state properties between superconducting and non-superconducting samples [43, 44, 45, 46]. The proximity of the magnetic ordered state to the SC in the phase diagram implies spin fluctuation-driven pairing mechanism, similar to the cuprates and the iron-based superconductors.

The density-functional theory (DFT) based calculations suggest that the low-energy degrees of freedom for the AH_xCr₃As₃ is the Cr-3d orbitals[43, 44, 45, 46, 47], similar with those of the 233 family[2, 20, 22, 3, 6, 24, 13, 25, 15, 26, 27, 14, 28, 16, 17, 29], implying the strong electron-electron interactions can account for the origin of the SC in the system. Furthermore, the presence of H induces effective electron doping, which confirms that H doping has the effect of raising Fermi level [43, 47, 44]. A remarkable property of the band structure of the AH_xCr₃As₃ is the presence of the Lifshitz transition under moderate hole doping, during which its 3D FSs are changed to two disconnected Q1D FS sheets [44, 47]. It’s clarified in our previous work[47] based on the random-phase approximation (RPA) [48, 49, 50, 51, 52, 54, 53, 55, 56] that the presence of the type-II VHS[57, 58, 59, 60, 61] during the Lifshitz transition doping regimes has driven the triplet $p_{z}$-wave pairing symmetry. It’s interesting to investigate the topological properties of this triplet pairing state. However, as the spin-orbit coupling (SOC) has not been considered in that work, the three spin components of the triplet pairing are exactly degenerate, which leads to uncertainty of the superconducting ground state.

In this paper, we report our RPA study on the pairing symmetry and the topological properties of the KH_xCr₃As₃, based on our first-principle DFT band structure considering SOC. Our main results are as follow. Due to the inversion symmetry of the material, the double-folded spin degeneracy survives the SOC. The main effect of the SOC on the band structure lies in the elimination of the double degeneracy along the $\Gamma$-$A$ and the $K$-$H$ lines. Starting from our previous six-orbital tight-binding (TB) model[47] without SOC, we add the dominating on-site SOC conserving the spin-U(1) symmetry (SUS), leading to a band structure well consistent with the DFT one. Then we add the multi-orbital Hubbard interactions into this model and perform the RPA calculations. The leading pairing symmetry obtained from our RPA calculations is triplet $p_{z}$-wave pairing in the $S_{z}=0$ channel, i.e. the $\uparrow\downarrow+\downarrow\uparrow$. This pairing state belongs to the SUS-protected time-reversal-invariant (TRI) topological nodal-line SC defined in the Ref[62], which includes the 233 family as a good example. Determined by the momentum-dependent topological winding number, the whole surface Brillouin zone (BZ) is covered with topological flat bands, with different regimes covered with different numbers of flat bands, as shown by our surface spectrum in the $(001)$ surface BZ. Our results appeal more experimental access into this intriguing family.

The crystal structure of KHCr₃As₃ is shown in Fig. 1(a). This compound consists of an infinite linear chain of DWSN [(Cr₃As₃)^2-]_∞ extended along the c-axis, forming a Q1D structure, which is interconnected by K⁺ cations. The hydrogen atoms locate at the center of the Cr₆ octahedrons in such a DWSN structure[43, 44, 45, 46, 47]. Compared with the crystal structure of K₂Cr₃As₃ [2, 20, 22, 3, 6, 24, 13, 25, 15, 26, 27, 14, 28, 16, 17, 29] and KCr₃As₃ [56, 33, 34, 35, 39, 43, 36, 37, 30, 32, 31, 38], KHCr₃As₃ is more similar to KCr₃As₃, which has a centrosymmetric structure with space group $P6_{3}/m$ (No.176)) and point group $C_{6h}$.

The band structure of KHCr₃As₃ is calculated using DFT theory, adopting the ultra-soft pseudopotential (USPP) of Perdew-Burke-Ernzerhof (PBE) [63] in QUANTUM ESPRESSO (QE) [64]. We use $5\times 5\times 11$ sparse k-points in the entire Brillouin Zone (BZ) for self-consistent calculation. We set the truncated kinetic energy of the wave function to 60Ry and the charge density to 400Ry. The experimental lattice parameters of the crystal structure of KHCr₃As₃ were adopted: $a=9.09481$ Å  and $c=4.17703$ Å [43], with only the atomic positions within the unit cell to be determined by optimization.

The obtained band structure of KHCr₃As₃ without considering SOC is shown in Fig. 1(b) by the red curves. It can be seen that there are five bands near the Fermi level, labeled as $\alpha,\beta,\gamma,\delta$ and $\epsilon$, respectively. The three bands across the Fermi level are the $\alpha$-, $\beta$- and $\gamma$- bands. These bands form two Q1D Fermi surfaces ($\alpha$ and $\beta$) and a 3D Fermi surface ($\gamma$), as shown in Fig. 1(c). Note that the $\beta$ and $\gamma$- Fermi surfaces (FSs) intersect at the $k_{x}=k_{y}=0,k_{z}=k_{F}$ point on the $\Gamma$-$A$ line. Compared with the band structure of KCr₃As₃ [56, 33, 34, 35, 39, 43, 36, 37, 30, 32, 31, 38], it can be found that the Fermi level rises when H is inserted, which reflects that H is used as an electron donor in KHCr₃As₃, consistent with Ref[43].

After the band structure calculation, a six-band low-energy tight-binding model has been constructed in our previous work[47]. Briefly, the relevant orbital degrees of freedom include the Cr-3d${}_{z^{2}}$, 3d${}_{x^{2}-y^{2}}$ and 3d_xy orbitals from the two inequivalent Cr atoms within each unit cell, and the effective hopping parameters are obtained through the calculations of maximal localization for these orbitals using Wannier90 [65]. The obtained TB Hamiltonian with $C_{6h}$ point-group symmetry can be expressed in the momentum space as,

Here $\mu,\nu=1,\cdots,6$ indicate the orbital-sublattice indices, containing the $d_{z^{2}}$, $d_{x^{2}-y^{2}}$ and $d_{xy}$ orbitals of A or B sublattices. The elements of the $h(\bm{k})$ matrix is obtained through Fourier transformation using the real-spacing hopping-integral data provided in our previous work[47]. The band structure thus obtained is well consistent with that of the DFT at low energy near the Fermi level, see Ref[47].

Using DFT calculations, we obtain the band structure with SOC. The band structures with and without SOC are compared in Fig.  1(b). From the comparison, we find the two band structures are overall similar except that the degeneracies along the high-symmetric lines are lift up by the SOC. Without SOC, the $\beta$- and $\gamma$- bands are degenerate at the $\Gamma$- and $K$- points, or in general along the $\Gamma$-$A$ and $K$-$H$ lines. Similarly, the $\delta$- and $\varepsilon$- bands are also degenerate along these high-symmetric lines. However, with SOC, the energies of the $\beta$- and $\gamma$- bands at the $\Gamma$- and $K$- points are split by about 0.02eV and 0.05eV, and the energies of the $\delta$- and $\varepsilon$- bands at the two momenta are split by about 0.08eV and 0.006eV, respectively. Note that the spin degeneracy is maintained with SOC, due to the combination of the inversion and time-reversal (TR) symmetries. The FS with SOC is shown in Fig. 1(d). Now the $\beta$- and $\gamma$- FSs no longer intersect at the $(0,0,k_{F})$ point on the $\Gamma$-$A$ line.

In the following studies involving electron-electron interactions, we need a simple TB formula for the SOC part of the Hamiltonian. For simplicity, we only consider the dominating part of the SOC, i.e. the on-site SOC term. From combination of the rotation symmetry, the mirror-reflection symmetry about the $xy$-plane, and the TR symmetry of the system, the only possible on-site SOC takes the following formula[62, 27],

This SOC Hamiltonian also satisfies the inversion symmetry $\hat{R}$ which dictates $c^{\dagger}_{\mathbf{i}\mu\sigma}\to c^{\dagger}_{\left(\hat{R}\mathbf{i}\right)\mu\sigma}$.

Our full model band Hamiltonian with SOC is given by

Setting $\lambda_{\text{SOC}}=10$ meV, the obtained band structure and FS for this model are shown in Fig. 2 (a) and (b) respectively. Comparing Fig.  2 (a) with Fig.  1 (b), it’s found that this band structure is qualitatively consistent with the DFT one with SOC. Firstly, the spin degeneracy is kept in this band structure, i.e. $\varepsilon_{\mathbf{k}\uparrow}^{\alpha}=\varepsilon_{\mathbf{k}\downarrow}^{\alpha}\equiv\varepsilon_{\mathbf{k}}^{\alpha}$, due to the combined TRS and inversion symmetry, where $\varepsilon_{\mathbf{k}\sigma}^{\alpha}$ denotes the band energy for the spin $\sigma$ and momentum $\mathbf{k}$ in the $\alpha$-th band. Secondly, the degeneracies between the $\beta$- and $\gamma$- bands and those between the $\delta$- and $\varepsilon$- bands at the $\Gamma$- and $K$- points and along the $\Gamma$-$A$ lines (as well as along the unshown $K$-$H$ lines) are lift up. Furthermore, the FS for the up or down spin shown in Fig. 2 (b) for this model is also qualitatively consistent with the DFT FS with SOC shown in Fig. 1(d).

According to our previous RPA result for this material without considering SOC[47], the triplet $p_{z}$-wave pairing is the leading pairing symmetry whose three spin components are exactly degenerate. Here we consider the SOC to lift up such degeneracy. We adopt the following extended Hubbard model Hamiltonian in our study,

Here, the interaction parameters $U$, $V$, and $J_{H}$ denote the intra-orbital, inter-orbital Hubbard repulsion, and the Hund’s rule coupling (as well as the pair hopping) respectively, which satisfy the relation $U=V+2J_{H}$. Note that this Hamiltonian conserves the SUS, i.e. the total $S_{z}$ is a good quantum number.

The pairing symmetry of this model is obtained through the standard multi-orbital RPA approach considering SOC[27]. The procedure of this approach is briefly introduced as follow. Firstly, we need to calculate the susceptibility tensor as function of the momentum. As the Hamiltonian (4) conserves the SUS, the spin indices of the nonzero susceptibility tensor elements can only be $\uparrow\uparrow\uparrow\uparrow$, $\uparrow\uparrow\downarrow\downarrow$, $\downarrow\downarrow\uparrow\uparrow$, $\downarrow\downarrow\downarrow\downarrow$, $\uparrow\downarrow\downarrow\uparrow$, $\downarrow\uparrow\uparrow\downarrow$. The renormalized susceptibility tensor can be obtained from the bare one via the Dyson equation in the RPA level. Secondly, we need to calculate the effective pairing interaction mediated by exchanging the particle-hole excitations whose propagators are represented by these susceptibility tensor elements. Note that with SOC here, the effective pairing interactions in the three spin channels of the triplet pairing are no longer degenerate. Finally, we need to solve the linearized mean-field (MF) gap equation near the SC $T_{c}$ involving the obtained effective pairing interaction function to obtain all possible pairing eigenvalues $\lambda$ and corresponding pairing eigenvectors. The pairing eigenvectors corresponding to the largest pairing eigenvalue $\lambda$ determines the pairing symmetry, and the $T_{c}$ is related to $\lambda$ via $T_{c}\sim W_{\text{band}}e^{-1/\lambda}$, with $W_{\text{band}}$ denoting the low-energy band width.

The doping $x$ dependences of the pairing eigenvalues $\lambda$ for the five leading pairing symmetries are shown in Fig.3 in the doping regime $x\in(0.4,1)$ wherein definite experimental evidences for SC have been detected. Note that for the triplet pairings, due to the SUS, the $S_{z}$ of the Cooper pair is a good quantum number, which suggests that we can get the $S_{z}=1$ (i.e. the $\uparrow\uparrow$), $S_{z}=-1$ (i.e. the $\downarrow\downarrow$) or $S_{z}=0$ (i.e. the $\uparrow\downarrow+\downarrow\uparrow$) pairing channels in the presence of the SOC. While the $S_{z}=1$ and $S_{z}=-1$ pairing channels are still degenerate due to the spin degeneracy, they are not degenerate with the $S_{z}=0$ pairing channel. As a result, Fig.3 shows that the triplet $p_{z}$-wave pairing in the $S_{z}=0$ (i.e. the $\uparrow\downarrow+\downarrow\uparrow$) channel is the leading pairing symmetry. This pairing phase diagram can be compared with that without SOC provided in Ref[47].

The BCS-mean-field (MF) Hamiltonian describing the triplet $p_{z}$-wave pairing in the $S_{z}=0$ channel reads,

Here the real gap function $\Delta_{\mathbf{k}}^{\alpha}$ obtained from our RPA calculations satisfies the $p_{z}$-wave pairing symmetry. This symmetry requires the gap function to be zero on the $k_{z}=0$ plane. Therefore, the FS cut on this plane gives the gap nodal lines. It’s interesting to investigate the topological properties of such nodal-line SC.

It’s easily checked that the MF Hamiltonian (5) satisfies the combined TRS and SUS. The topological classification of nodal-line superconductors hosting the two symmetries can be attributed to the following momentum-dependent topological invariant $Z(k_{x},k_{y})$ defined as[62]

Here $\theta_{\mathbf{k}\alpha}$ denotes the complex phase angle of $\varepsilon_{\mathbf{k}}^{\alpha}+i\Delta_{\mathbf{k}}^{\alpha}$, and the integer $I_{\alpha}(k_{x},k_{y})\equiv\frac{1}{2\pi}\int^{\pi}_{-\pi}d\theta_{\mathbf{k}\alpha}$ represents the winding number of $\theta_{\mathbf{k}\alpha}$ along a closed path in the BZ which passes through $(k_{x},k_{y},0)$ and is perpendicular to the $k_{z}=0$ plane. As topological invariant, $2Z(k_{x},k_{y})$ counts the number of gapless surface modes at the momentum $(k_{x},k_{y})$ in the $(001)$-surface BZ.

As an example, we illustrate how to calculate $I_{\gamma}(k_{x},k_{y})$ for two typical $(k_{x},k_{y})$ points at the doping level $x=0.85$. The $\gamma$- FS for $x=0.85$ is shown in Fig. 4(a). This FS consists of two Q1D FS sheets nearly parallel to the $k_{z}=0$ plane, connected by six hollow tubes (only three are visible in the figure) nearly perpendicular to the $k_{z}=0$ plane. The six hollow tubes cut the $k_{z}=0$ plane to form six pockets, as shown in Fig. 4(b). The $(k_{x1},k_{y1})$ locates near the $\Gamma$- point, passed through by the closed vertical path $L_{1}$, and the $(k_{x2},k_{y2})$ locates within one of the six pockets, passed through by the path $L_{2}$. The $k_{z}$-dependence of $\Delta^{\alpha}_{\mathbf{k}}$ averaged on the FSs is shown in Fig. 4(c), where the sign of $\Delta^{\alpha}_{\mathbf{k}}$ follows that of $k_{z}$. This character is also satisfied for the $\Delta^{\alpha}_{\mathbf{k}}$ with fixed $(k_{x},k_{y})$, as the Q1D $p_{z}$-wave pairing gap function value doesn’t obviously depend on $(k_{x},k_{y})$.

Let’s first evaluate $I_{\gamma}$ for $(k_{x1},k_{y1})$. Fig. 4(a) shows the locations of five typical momenta $A_{1}\sim E_{1}$ along $L_{1}$, which intersects with the FS twice. The energies $\varepsilon_{\mathbf{k}}^{\gamma}$ for $C_{1}$ and $D_{1}$ are positive, and those for the remaining three are negative. The gap functions $\Delta_{\mathbf{k}}^{\gamma}$ for $D_{1}$ and $E_{1}$ are positive, and those for $B_{1}$ and $C_{1}$ are negative. Due to the $p_{z}$-wave symmetry, the $\Delta_{\mathbf{k}}^{\gamma}$ for $A_{1}$ is zero. Fig. 4(d) shows the locations of the $\varepsilon_{\mathbf{k}}^{\gamma}+i\Delta_{\mathbf{k}}^{\gamma}$ for the five momenta $A_{1}\sim E_{1}$ in the complex plane expanded by the $\varepsilon_{\mathbf{k}}^{\gamma}$- axis and the $\Delta_{\mathbf{k}}^{\gamma}$- axis, which form a closed trajectory circling around the original point, suggesting that the phase winding number of $\varepsilon_{\mathbf{k}}^{\gamma}+i\Delta_{\mathbf{k}}^{\gamma}$ along $L_{1}$ is 1. Therefore, we have $I_{\gamma}(k_{x1},k_{y1})=1$. Then let’s investigate $I_{\gamma}(k_{x2},k_{y2})$. As $L_{2}$ locates within one of the six hollow tubes, it does not intersect with the $\gamma$-FS. Consequently, the $\varepsilon_{\mathbf{k}}^{\gamma}$ for the five momenta $A_{2}\sim E_{2}$ are all negative. On the other hand, the signs for the $\Delta_{\mathbf{k}}^{\gamma}$ for the five momenta are the same as those for $A_{1}\sim E_{1}$. As a result, the trajectory formed by the locations of the $\varepsilon_{\mathbf{k}}^{\gamma}+i\Delta_{\mathbf{k}}^{\gamma}$ for $A_{2}\sim E_{2}$ in the complex plane doesn’t circle around the original point, as shown in Fig. 4(e), which leads to zero winding number of the phase angle of $\varepsilon_{\mathbf{k}}^{\gamma}+i\Delta_{\mathbf{k}}^{\gamma}$ along $L_{2}$. Thus, we have $I_{\gamma}(k_{x2},k_{y2})=0$.

Besides the $\gamma$- FS, the compound with $x=0.85$ hosts two other Q1D FSs, i.e. the $\alpha$- and $\beta$- FSs. These two FSs take simple shapes, with each hosting two Q1D FS sheets almost parallel to the $k_{z}=0$ plane, similar as those for $x=1$ shown in Fig. 1 (d) and Fig. 2 (b). For arbitrary $(k_{x},k_{y})$, the vertical path $L$ passing through it always intersects with the $\alpha$- and the $\beta$- FS twice. On the other hand, the gap functions for these two bands are similar with that for the $\gamma$ band in that the sign of $\Delta_{\mathbf{k}}^{\alpha/\beta}$ follows that of $k_{z}$. Following the above illustrated procedure, we can easily obtain $I_{\alpha/\beta}(k_{x},k_{y})=1$ for an arbitrary $(k_{x},k_{y})$. The total topological invariant $Z(k_{x},k_{y})$ is the sum of $I_{\alpha}(k_{x},k_{y})$, $I_{\beta}(k_{x},k_{y})$ and $I_{\gamma}(k_{x},k_{y})$, which takes different nonzero values in different regimes of the $(k_{x},k_{y})$ plane. For example, for the two typical momenta shown in Fig. 4 (b), we have $Z(k_{x1},k_{y1})=3$ and $Z(k_{x2},k_{y2})=2$.

The nontrivial topological invariant $Z(k_{x},k_{y})$ leads to exactly-flat surface bands on the $(001)$ surface[62]. Concretely, for nonzero $Z(k_{x},k_{y})$, there would be $2Z(k_{x},k_{y})$ surface bands which are exactly flat when they pass through the momentum $(k_{x},k_{y})$ in the surface BZ[62]. Such knowledge guides us to study the surface spectrum on the $(001)$ surface. For this purpose, we adopt open boundary condition along the $z$- axis and periodic ones along the $x$- and $y$- axes, under which we Fourier transform the BCS-MF Hamiltonian (5) to the $\left(k_{x},k_{y},i_{z}\right)$ space. Then we diagonalize the Bogoliubov-de-Genes equation for the fixed $(k_{x},k_{y})$ to get the corresponding energy spectrum. The obtained energy spectrum as function of $k_{x}$ and $k_{y}$ is shown in Fig. 5 along the high symmetric lines in the $(k_{x},k_{y})$ plane. To enhance visibility, we have enlarged the gap amplitude by an order of magnitude. As shown in Fig. 5, some regime in the $(k_{x},k_{y})$-plane is covered by six flat bands while the remaining regime is covered by four flat bands. We have checked that these flat bands are surface bands as the corresponding wave functions are localized on the $(001)$ surface. Comparing Fig. 5 with Fig. 4 (b), it’s found that the regime within the six Fermi pockets containing the typical $(k_{x2},k_{y2})$ point wherein $Z(k_{x},k_{y})=2$ is covered by four flat bands, and the remaining regime containing the typical $(k_{x1},k_{y1})$ point wherein $Z(k_{x},k_{y})=3$ is covered by six flat bands.

The FS topology of KH_xCr₃As₃ can be drastically changed upon tuning $x$ which is experimentally accessible[44]. On the other hand, Fig. 3 shows that the triplet $p_{z}$-wave pairing in the $S_{z}=0$ channel is always the leading pairing symmetry, belonging to the SUS-protected TRI topological nodal-line SC with flat surface bands. Consequently, the distribution of the number of flat bands in the surface BZ will be easily engineered through tuning the hydrogen doping. In Fig. 6, the FSs and surface spectra for $x=1$, $x=0.61$ and $x=0.45$ are shown. Clearly, the different FS topologies for the three different hydrogen-doping levels lead to different distributions of the number of flat bands in the surface BZ. For all the hydrogen-doping levels we studied, the whole surface BZ is covered with flat bands. Particularly in the case of $x=0.61$, due to the simple FS shape with three Q1D FSs, the whole surface BZ is covered with six flat bands.

The SUS-protected TRI topological nodal-line SC with flat surface bands predicted here for the AH_xCr₃As₃ family is similar with that predicted for the A₂Cr₃As₃ family, which can be detected in the surface conductance spectrum in the scanning tunneling microscope (STM) experiment by the sharp zero-bias peak[62]. There can be inter-site spin-flipping SOC in the two families, which however should be much weaker than the on-site SOC adopted here. In the presence of such weak spin-flipping SOC, the SUS would be slightly broken, leading to slightly broadening of the zero-bias peak in the STM, which is still well detectable[62].

In conclusion, we have performed a RPA-based study on the pairing symmetry and topological properties of the KH_xCr₃As₃, based on our DFT band structure considering SOC. Our results yield triplet $p_{z}$-wave pairing in the $(\uparrow\downarrow+\downarrow\uparrow)$ channel to be the leading pairing symmetry. This pairing state belongs to the SUS-protected TRI topological nodal-line SC hosting flat surface bands. Determined by the $\mathbf{k}$-dependent topological invariant $Z(k_{x},k_{y})$, the whole $(001)$-surface BZ is covered with topological flat bands, with different regimes covered with different numbers of flat surface bands, as shown by our surface spectra, which can be verified by experiments.