Pressure-induced superconductivity in kagome metal CsCr₃Sb₅: Role of spin–orbit coupling and inter-orbital spin fluctuations

Kagome lattice, a network of corner-sharing triangles, has been extensively studied in condensed matter physics for its ability to realize unconventional magnetic and electronic states.^[1–10] Recently, significant attention has been drawn to the newly discovered quasi-two-dimensional kagome metals AV₃Sb₅ (A = K, Rb, Cs) owing to their rich phenomena.^[11–16] Reported properties include a giant anomalous Hall effect,^[17] an unconventional charge density wave (CDW) phase accompanied by time-reversal symmetry breaking,^[18] pair density wave,^[19] electronic nematicity,^[20] and superconductivity (SC).^[20–24]

The newly synthesized kagome metal CsCr₃Sb₅ [Fig. 1(b)], a structural analogue of the AV₃Sb₅ family, exhibits a distinct phase diagram.^[25] At ambient pressure, it undergoes a transition near 55 K that simultaneously involves CDW, spin density wave (SDW), and orbital nematicity.^[25,26] Under applied pressure, these density-wave (DW) phases are progressively suppressed and coexist with SC around 3.65–3.80 GPa before being eliminated, giving rise to a superconducting dome that peaks at 4.2 GPa with a transition temperature of 6.4 K and vanishes above 10 GPa.^[25] Moreover, the absence of a soft acoustic phonon mode suggests that the SC phase is unconventional.^[27]

Unlike AV₃Sb₅, where the Fermi level resides near van Hove singularities,^[16] CsCr₃Sb₅ hosts a flat band derived from the Cr d orbitals near the Fermi level and van Hove singularities with binding energy of ∼ 0.5 eV, observed in both angle-resolved photoemission spectroscopy experiments^[28–30] and theoretical calculations.^[31–34] A random phase approximation (RPA) study^[35] employing a 30-orbital tight-binding model for CsCr₃Sb₅ suggests an s^±-wave pairing instability with competing (d_xy, d_x²-y²)-wave symmetry. However, due to the significantly smaller differences in the on-site energies between d orbitals in CsCr₃Sb₅ compared to AV₃Sb₅^[36] (see Appendix A), spin–orbit coupling (SOC) in CsCr₃Sb₅ may play a more prominent role than in AV₃Sb₅. Indeed, introducing a small SOC can lead to a qualitative change in the Fermi surface (FS) structure [see Figs. 1(c) and 1(d)]. Therefore, similar to some iron-based superconductors,^[37–40] although the low-energy electronic structure is dominated by 3d electrons, SOC could have a significant impact on the superconducting pairing symmetry in CsCr₃Sb₅.

In this paper, using an effective four-orbital tight-binding model with local multi-orbital Hubbard interactions and an RPA analysis, we investigate the effect of SOC on the superconducting pairing symmetry in pressurized CsCr₃Sb₅. Our results reveal that the E_2g symmetry is the dominant pairing symmetry both with and without SOC, except in the regime of very small Hund’s coupling without SOC, where the A_2g symmetry becomes stabilized. The enhancement of the inter-orbital spin fluctuations stabilizes the E_2g symmetry. When the SOC is included, spin fluctuations remain the primary driving force for pairing, with orbital fluctuations playing a subsidiary role. In both cases, the dominant pairing interaction originates from the β band, primarily derived from the d_x²-y² orbital.

To study the superconductivity of CsCr₃Sb₅, we first perform the first-principles calculation (see details in Appendix A) on CsCr₃Sb₅ under external pressure of 5 GPa. Structural relaxation yielded a hexagonal lattice [Fig. 1(b)] with parameters a = 5.3235 Å and c = 8.6439 Å. Figure 1(a) presents the resulting electronic band structure and projected density of states, revealing the dominant contributions near the Fermi energy from the Cr-d_xz, Cr-d_yz, Cr-d_x²-y², and Sb₁-p_z orbitals. Furthermore, the bands originating from the Cr d orbitals exhibit minimal dispersion along the Γ–A path, resulting in the flat bands. To capture these features, we employed the effective four-orbital tight-binding model H_t established in our previous work,^[34] which incorporates the Cr-d_xz, Cr-d_yz, Cr-d_x²-y², and Sb₁-p_z orbitals, to fit the DFT-calculated bands in the k_z = 0 plane at 5 GPa. The hopping parameters for CsCr₃Sb₅ under pressure are given in Appendix A. As shown in Fig. 1(a), the effective tight-binding model qualitatively reproduces the bands in the DFT calculations near the Fermi level.

To investigate the influence of SOC on superconducting pairing function, we restrict the SOC to on-site terms, defined as HSOC=λSOC∑ili⋅si, where s_i and l_i denote the spin and orbital angular momentum operators, respectively, for the i-th Cr atom (see Appendix A). The total Hamiltonian is then expressed as H = H_t + H_SOC + H_I, where H_I corresponds to the local multi-orbital Hubbard-type interaction for Cr d-electrons (see Appendix B), parametrized by the Hubbard U, inter-orbital Coulomb interaction U′ and Hund’s coupling J_H. Based on the single-particle Hamiltonian H₀ = H_t + H_SOC, we present the orbital-resolved FSs in Figs. 1(c) and 1(d), calculated without and with SOC, respectively. In both cases, the bands denoted by α originate predominantly from the Sb₁-p_z and Cr-d_yz orbitals near the Fermi level, while the bands β and γ are primarily contributed by the Cr-d_x²-y² and Cr-d_xz orbitals, respectively.

To explore the superconducting pairing symmetry in CsCr₃Sb₅ under pressure, we solve the linearized gap equation projected onto the FS, formulated as an eigenvalue problem:^[40–43]

The pairing kernel is given by

where |vF(kF′)| corresponds to the Fermi velocity at the Fermi momentum kF′, V_BZ denotes the volume (or area in two dimensions) of the Brillouin zone, and the indices ξ label the (pseudo)spin-singlet (ξ=0) and triplet (ξ = x,y,z) pairing channels. Γξξ′(kF,kF′) represents the effective pairing interaction projected onto the FS (see Appendix B). The linearized gap equation (Eq. (1)) yields eigenvalues λ and eigenvectors Δ_ξ(k_F), which correspond to the superconducting gap functions. A leading eigenvalue reaching λ = 1 indicates the onset of a superconducting instability, unless it is suppressed by a competing instability in the particle–hole channel. To examine the possibility of such particle–hole instabilities within the RPA framework, we diagonalize the matrix −V^χ^0(q,0), where χ^0(q,0) is the bare particle–hole susceptibility and V^ is the bare interaction vertex (see Appendix B). An eigenvalue of −V^χ^0(q,0) attaining unity, corresponding to the Stoner criterion, signals the emergence of a particle–hole instability.

Since the spin-fluctuation-mediated pairing interaction Γξξ′(kF,kF′) is directly related to the particle–hole susceptibility χ^(q=kF−kF′,0) (see Appendix B), it is crucial to analyze which particle–hole fluctuations provide the dominant contribution. To this end, we diagonalize the static particle–hole susceptibility χ^(q,0) and consider the leading eigenvalues λχmax(q). When the eigenvalue λχmax(q) is strongly enhanced at a characteristic wave vector Q, the scattering amplitude Γξξ′(kF,kF′) of a Cooper pair from (kF′,−kF′) to (kF′+Q,−kF′−Q) is dominated by the particle–hole fluctuation associated with the corresponding eigenvector at Q. To clarify the degrees of freedom underlying this fluctuation, we further decompose the static particle–hole susceptibility χ^(q,0) into orbital-resolved spin and charge components. The charge contribution is negligible compared to the spin channel and therefore not shown. The total spin susceptibility is defined as χS=∑μνχμνS(q), with orbital indices μ and ν. The individual component is given by

where l and l′ denote sublattice indices. The spin susceptibilities of the d_x²-y²-orbital and {d_xz,d_yz} orbitals are written as χdx20y2S=χdx20y2,dx20y2S(q) and χdxz+dyzS=∑μν∈{dxz,dyz}χμνS(q), respectively. Similarly, the pairing kernel M^n′n′ξξ′(kF,kF′) can be written as^[40]

where g_R,ξ(k_F) and gL,ξ′(kF′) are the right and left eigenvectors of Mξξ′(kF,kF′), respectively. Symbol λ₁ denotes the largest eigenvalue of the pairing kernel. The ellipsis represents contributions from subleading eigenvalues. The singlet component of the pairing kernel projected onto the leading superconducting channel is given by λ1gR,0(kF)gL,0(kF′). In the following calculations, since the triplet component of the pairing kernel is found to be negligible compared to the singlet component, we focus solely on the singlet channel. In this work, we employ a 100 × 100 q-mesh and set the temperature to k_BT = 1/β = 0.001 eV.

Figure 2 exhibits the evolution of the leading instability channels as a function of the Hubbard interaction U and Hund’s coupling J_H. In Figs. 2(a) and 2(b), the left solid line marks the boundary separating the normal metal phase from the superconducting phase, determined by the condition that the leading eigenvalue λ of the gap equation (1) attains unity. The right solid line indicates the phase boundary between the superconducting and DW phases, determined by the Stoner criterion that the largest eigenvalue of −V^χ^0(q,0) reaches unity. Within the normal metal phase where λ < 1, that is, the region to the left of the superconducting boundary, the dominant particle–particle instability is additionally presented. In the absence of SOC, the leading superconducting pairing initially favors A_2g symmetry for small U and J_H/U. However, at larger values of J_H/U, the leading pairing symmetry changes from A_2g to E_2g with the increase of U. With further enhancement of U, the system ultimately enters a SDW. When SOC is introduced, the leading superconducting pairing favors instead the E_2g symmetry, and at sufficiently large U, a particle–hole instability emerges as the dominant ordering tendency.

To elucidate the superconducting gap functions, we first analyze the case without SOC at two representative parameter sets: J_H/U = 0.2, U = 0.9964 eV and J_H/U = 0.02, U = 1.0853 eV [see Figs. 2(a) and 2(b)]. As illustrated in Fig. 2(c), these values of interaction parameters correspond to points where the leading eigenvalue λ of the linearized gap equation reaches unity. Figure 3(a) displays the gap function with A_2g symmetry, while Figs. 3(b) and 3(c) exhibit the two degenerate solutions with E_2g symmetry. The characteristic wave vectors Q₁ = (0.2252,−0.3900) and Q₂ = (0.4070,−0.2350), which connect regions with opposite signs of the gap function, coincide with the peak positions of the leading eigenvalues λχmax(q) of particle–hole susceptibility χ^(q,0) [Figs. 3(d) and 3(g)]. In the spin-fluctuation-mediated superconducting pairing scenario, the peak of eigenvalue λχmax(q) at wave vector Q (=Q₁, Q₂) indicates that the pairing scattering is dominated by particle–hole fluctuation associated with the corresponding eigenvector at Q, and that the pairing interaction satisfies Γ00(kF,kF′)>0 where kF−kF′=Q. To maximize the eigenvalue λ of the gap equation (1), the gap function must change sign between two FS points connected by these wave vectors (note the minus sign in Eq. (2)),^[40,42,43] i.e., Δ0(kF′+Q)Δ0(kF′)<0. For both sets of interaction parameters, the primary peaks for λχmax(=λχmax(q)) are located at Q₁ and Q₂, as shown in Figs. 3(e) and 3(h). These peaks are predominantly associated with spin fluctuations, as evidenced by the significant enhancement of the total spin susceptibility χ^S [Figs. 3(e) and 3(h)]. The spin fluctuations at Q₁ and Q₂ originate primarily from the d_x²-y² orbital, as indicated by the pronounced peaks in χdx20y2S [Figs. 3(e) and 3(h)]. These peaks are closely tied to the nesting of the β band, which is mainly derived from the d_x²-y² orbital [Fig. 1(c)]. The Fermi-surface-nesting origin of these peaks is further confirmed by examining the properties of the bare spin susceptibility χ0,dx20y2S, whose characteristics are entirely determined by the FS, with no contribution from interaction effects. We find that χ0,dx20y2S exhibits clear maxima at Q₁ and Q₂ [inset of Fig. 3(e)], whereas χ0,dxz+dyzS lacks corresponding peaks. At small Hund’s couplings, such as J_H/U = 0.02, the peak at Q₁ exceeds that at Q₂, stabilizing the A_2g pairing symmetry. In contrast, at larger Hund’s couplings, such as J_H/U = 0.2, the Q₂ peak becomes dominant, favoring the E_2g pairing symmetry. This evolution reflects the interplay between intra- and inter-orbital spin fluctuations. For weak Hund’s coupling, the Hubbard interaction U enhances intra-orbital spin fluctuations, and the bare spin susceptibility in the dominant d_x²-y²-orbital channel exhibits a stronger peak at Q₁ than that at Q₂ [inset of Fig. 3(e)], leading to the dominance of Q₁ in λχmax. By contrast, at stronger Hund’s coupling, the inter-orbital spin fluctuations are increased, as evidenced by the enhanced peaks for χdxz+dyzS [Fig. 3(h)]. This shifts the leading peak to Q₂. To shed light on the microscopic process behind the enhancement of inter-orbital spin fluctuations, we analyze the largest eigenvalue λχmax(q) of χ^(q,0) together with its eigenvector v_lν(q). The largest eigenvalue can be written as

where D(q)=∑ll′l1ν1μνvlν*(q)[V^]l1ν1lν[χ^0]l′μl1ν1(q,0)vl′μ(q). The corresponding particle–hole fluctuation mode is represented by the operator O(q)=∑lνvlν(q)Slνz(q), with the restriction to the S^z component enforced by spin-rotation symmetry without SOC. As shown in Fig. 5(a), the dominant weight of the eigenvector at Q₁ originates from the d_x²-y² orbital on both A and B sublattices (with |vA,x20y2(q)|2=|vB,x20y2(q)|2), followed by subdominant contributions from the corresponding d_xz orbitals. At Q₂, the weight is primarily from the B-sublattice d_x²-y² orbital, with the next-largest and third-largest contributions from B-sublattice d_xz and d_yz orbitals, respectively. Accordingly, we restrict the summation in D(q) to these dominant degrees of freedom. With increasing Hund’s coupling, the weight shifts from the d_x²-y² orbital to the d_xz and d_yz orbitals. Simultaneously, −D(q) increases and approaches unity, with a notably faster increase at Q₂ than at Q₁ at large Hund’s coupling, as shown in Fig. 5(a). This enhancement primarily stems from the contribution of the inter-orbital spin channel to −D(q). These results demonstrate that Hund’s coupling selectively enhances inter-orbital spin fluctuations, with a dominant effect at Q₂, which consequently causes the leading peak in λχmax to shift from Q₁ to Q₂. In addition to the above analysis indicating that the pairing interaction originates primarily from the scattering channels within the β band, this conclusion is also supported by the singlet component of the pairing kernel. As demonstrated in Figs. 3(f) and 3(i), the pairing interaction within the β band is the strongest.

Next, we discuss the cases with SOC. Since the results for λ_SOC = 25–40 meV^[44] are qualitatively similar, owing to the unchanged Fermi-surface topology, we present the results only for the case of λ_SOC = 25 meV in Fig. 4. When SOC is included, the leading superconducting instability favors E_2g symmetry. The two degenerate gap functions in this channel are illustrated in Figs. 4(a) and 4(b) for J_H/U = 0.2 and U = 0.98675 eV. The characteristic wave vectors Q₃=(0.3984,−0.2300), Q₄=(0.3984,0.000), and Q4′=(0.1867,−0.5748), linking Fermi-surface regions with opposite signs of the gap function, align with the maxima of λχmax [Fig. 4(c)]. Among these vectors, the leading and subleading peaks in λχmax occur at Q₃ and Q₄, respectively [Fig. 4(e)]. At Q₃, the dominant peak of λχmax is predominantly driven by spin fluctuations, consistent with the significant enhancement of the total spin susceptibility χ^S at this wave vector. The subleading peak at Q₄ reflects an increased contribution from orbital fluctuations. This is evident from the fact that although the peak of χ^S at Q₅=(0.4503,0.0000) is larger than that at Q₄, the situation is reversed in λχmax, suggesting an additional orbital contribution at Q₄. To further confirm the presence of additional orbital fluctuations, we compute the orbital susceptibility χ^O, defined as χO=∑ll′s1s2∑μνμ′ν′[χ^]l′ν′s2,l′μ′s2lμs1,lνs1(q,0)OνμOν′μ′ (see Appendix B). As shown in Fig. 5(b), the dressed orbital susceptibility χ^O exhibits a pronounced peak at Q₄, in contrast to its intensity at Q₅, verifying the existence of additional orbital fluctuations. Consistent with the case without SOC, the dominant pairing interaction remains within the β band, as confirmed by the singlet component of the pairing kernel [Fig. 4(d)]. However, compared to the scenario without SOC, the inter-band pairing interaction between the β and α bands is enhanced, aligning with the above analysis that additional orbital fluctuations at Q₄ play a role with increasing SOC. The interband pairing interaction between the β and γ bands is also enhanced, as shown in Fig. 4(d). One representative channel is the Q4′-mediated scattering, where Q4′ connects the β and γ bands [Fig. 4(a)], as illustrated in Fig. 4(c).

In summary, we have theoretically investigated the superconducting pairing symmetry in pressurized CsCr₃Sb₅ using an effective four-orbital tight-binding model combined with RPA analysis. Our results reveal that the strong Hund’s coupling enhances the inter-orbital spin fluctuations, and results in a shift in the maximum of spin fluctuations from the Γ–K direction to the Γ–M direction. The significantly enhanced inter-orbital spin fluctuations promote a superconducting gap function with E_2g symmetry. In the case of a small Hund’s coupling, where the inter-orbital spin fluctuations are weak, the leading pairing symmetry is A_2g symmetry. The inclusion of SOC enhances the orbital fluctuations and stabilizes the E_2g symmetry. In both scenarios with and without SOC, the β band, which is mainly derived from the d_x²-y² orbital, provides the prominent contribution to the pairing interaction.

The first-principles calculations were performed using Vienna ab initio simulation package (VASP)^[45] with the projector augmented-wave method.^[46] The exchange–correlation interactions were modeled using the Perdew–Burke–Ernzerhof (PBE) form of the generalized gradient approximation (GGA).^[47] A Γ-centered k-point grid of dimensions 12 × 12 × 5 was implemented in the self-consistent calculations. To avoid double counting of the SOC, the SOC was not included in these calculations.

Guided by symmetry analysis under the D_6h point Group assumption, we derive an effective tight-binding model, which can be represented as a 10 × 10 matrix^[34]

where A, B, and C denote the sublattice indices. The fermionic operator Ψ^σ(k) is defined as

where the subscript σ indicates the spin indices. The orbital indices are xz, yz, x²–y² for the Cr-d orbitals and p_z for the Sb₁-p_z orbital. The hopping matrices T^AA, T^AB, T^Ap, and T^pp are given by^[34]

where k1=k⋅a1,k2=k⋅a2, and k3=k⋅(a2−a1). a₁ and a₂ are the lattice vectors. The hopping matrices T^BB(T^Bp) and T^CC(T^Cp) can be obtained by changing the parameters (k₁,k₂,k₃) in T^AA(T^Ap) to (−k₂,−k₃,k₁) and (k₃,−k₁,−k₂), respectively. The explicit hopping parameters were determined through a least-squares optimization procedure that minimizes the deviation between the DFT-calculated bands and those derived from the effective tight-binding model. The hopping parameters for CsCr₃Sb₅ at 5 GPa are listed in Table 1.

In the chosen orbital basis d_xz, d_yz, and d_x²-y², the orbital angular momentum operator l_i of the Cr ions can be expressed in matrix form as

The bare particle–hole susceptibility χ^0 is defined as

where ω_n represents a fermionic Matsubara frequency and T_τ denotes the time-ordering operator with respect to the imaginary time. The average 〈⋯〉₀ depends on e^−βH₀ with inverse temperature β = 1/k_BT. H₀ is the quadratic part of the Hamiltonian H. The subscript α_n = (l_n,μ_n,σ_n) denotes the sublattice l_n, orbital μ_n, and spin σ_n indices. The Fourier-transformed fermionic operators are written as

where r_i denotes the lattice vector corresponding to the ith unit cell and δ_l represents an intra-unitcell vector with sublattice l. In the normal state, the bare particle–hole susceptibility χ^0 becomes

where unα(k)=〈α|kn〉 is the projection of eigenvector onto α state, corresponding to the eigenvalue ε_k,n of H₀. The function f(ε_k,n1) represents the Fermi distribution function. Within the RPA framework, the particle–hole susceptibility χ^ is given by

where V^(q) is the bare particle–hole vertex. For the multi-orbital Hubbard-type interaction, defined as

the bare particle–hole vertex V^(q) satisfies

where, in this instance, the q-independent interaction vertex V^(q) can be expressed as a combination of spin (U^s) and charge (U^c) vertices,

The non-zero elements of the charge and spin vertices are respectively given by^[40,42,43]

In this work, we adopt the relationships U′ = U – 2J_H and J_P = J_H.

The spin or charge susceptibilities could be obtained from the generalized particle–hole susceptibility χ^(q,iωn) by projecting onto the spin or charge channels. Specifically, the total spin susceptibility is given by

where σ^a=x,y,z represents the Pauli matrix. To obtain the other physical susceptibilities, such as orbital susceptibility, similar projections can be performed. For instance, the orbital susceptibility χ^O discussed in the main text is defined as

where the orbital matrix O is expressed in the basis of d_xz, d_yz, and d_x²-y² orbitals as

To examine the possible particle–hole instabilities, we consider the static susceptibility χ^(q,0). The instability occurs when det(1+V^χ^0(q,0))=0. Equivalently the emergence of a particle–hole instability is indicated when the largest eigenvalue of the matrix −V^χ^0(q,0) reaches unity. The computation is restricted to the submatrix of −V^χ^0(q,0) defined by the set of indices (α₁α₂,α₃α₄), where [V^]α3α4α1α2≠0, since these terms provide the dominant contributes near the critical point. A similar restriction is applied in evaluating the largest eigenvalue λχmax of the particle–hole susceptibility χ^(q,0).

For the analysis of superconductivity, we focus on the pairing function with zero center momentum. The effective pairing interaction by spin-fluctuations within RPA can be written as^[40,41]

Due to the prevalence of scattering events occurring in proximity to the FS, attention is directed towards the pairing interactions among the states located on the FS. The effective pairing interaction on the FS is given by

where subscript n₁ and ς₁ denote the band and pseudo-spin space, respectively. If the band indices n₁ or n₂ do not correspond to a band at the Fermi surface, the Fermi-surface projected pairing interaction Γξξ′(kF,kF′) necessarily vanishes. Therefore, the band indices need not be explicitly written in the pairing interaction Γξξ′(kF,kF′). The γ^ξ matrices are defined as

Here, σ^ξ=0 represents the 2 × 2 identity matrix. The symbol γ^0 signifies the pseudo-spin singlet, whereas γ^ξ denotes the pseudo-spin triplet, with ξ belonging to the set {x,y,z}. Owing to the gauge freedom inherent in the eigenstate resolution process, we implement the symmetry-adapted methodology outlined in Ref. [40] to circumvent non-unique gauge selections.

In order to categorize the pairing functions, the projection operator is utilized to map these functions onto the irreducible representation space associated with the corresponding point group. The projection operator is given by

where Dⁱ(R) is the i-th irreducible representation of point group G with an order of g. The notation P_R denotes the symmetry operator corresponding to the element R in the symmetry group G. l_i denotes the dimension of the i-th irreducible representation space. The crystal structure of CsCr₃Sb₅ possesses D_6h point group symmetry, which comprises twelve irreducible representations: A_1g, A_2g, B_1g, B_2g, E_1g, E_2g, A_1u, A_2u, B_1u, B_2u, E_1u, and E_2u.