Multiorbital singlet pairing and d + d superconductivity

Abstract

Recent experiments in multiband Fe-based and heavy-fermion superconductors have challenged the long-held dichotomy between simple s- and d-wave spin-singlet pairing states. Here, we advance several time-reversal-invariant irreducible pairings that go beyond the standard singlet functions through a matrix structure in the band/orbital space, and elucidate their naturalness in multiband systems. We consider the sτ₃ multiorbital superconducting state for Fe-chalcogenide superconductors. This state, corresponding to a d + d intra- and inter-band pairing, is shown to contrast with the more familiar d + id state in a way analogous to how the B- triplet pairing phase of ³He superfluid differs from its A- phase counterpart. In addition, we construct an analog of the sτ₃ pairing for the heavy-fermion superconductor CeCu₂Si₂, using degrees-of-freedom that incorporate spin-orbit coupling. Our results lead to the proposition that d-wave superconductors in correlated multiband systems will generically have a fully-gapped Fermi surface when they are examined at sufficiently low energies.

Introduction

Unconventional superconductivity of strongly correlated systems is centrally important in condensed matter physics, with the symmetry of the superconducting order parameter being a key issue of the field. This question appears to have reached a consensus in some notable instances. An example is the d-wave symmetry for the Cooper pairs in the well-studied Cu-based superconductors (SCs)^1,2. However, the pairing symmetry remains enigmatic in other classes of strongly correlated materials. For singlet superconductivity, the long-held dichotomy is between fully-gapped s- and nodal d-wave pairing states. However, it has been increasingly recognized that multi-band/orbital systems are inherently richer for pairings^3,4. A canonical setting for multiorbital spin-singlet pairings is the Fe-based superconductors^{5,6,7,8,9,10,11}, especially for the Fe-chalcogenide cases. Here, the discovery of an orbital-selective Mott crossover in the normal state^12,13 motivated the notion of orbital-selective superconductivity¹⁴. The latter opens up the possibilities for a variety of orbital-dependent pairing states, which have been studied in recent years both theoretically^{15,16,17,18,19,20} and experimentally^21,22,23. In addition, heavy fermion SCs, a class that includes about fifty members, have emerged as another prominent setting for singlet pairing states beyond the traditional possiblities²⁴.

Recent experiments have directly challenged the conventional s- and d-wave dichotomy. In alkaline Fe-selenides, inelastic neutron scattering^25,26 revealed signatures of in-gap spin resonances, whose characteristic wavevectors qualify them as a typical indicator of sign-changing d-wave order parameters^7,9,27,28,29. By contrast, ARPES studies have indicated fully-gapped superconductivity^30,31,32,33, even for a Fermi pocket near the center of the two-dimensional Brillouin zone (BZ), which appears consistent with s- wave symmetry. Understanding the Fe-chalcogenide SC is crucially important, since the Fe-based superconductivity with the highest superconducting transition temperature (T_c) occurs in this category.

A similar situation has emerged in the heavy-fermion system CeCu₂Si₂ (Ref. ³⁴). A host of properties, including the inelastic neutron-scattering spectrum²⁸, have traditionally been interpreted in terms of a sign-changing d-wave pairing state, yet recent specific heat³⁵ and London penetration depth^34,36 results at very low temperatures pointed toward a fully gapped Fermi-surface (FS).

It is surprising that the SC phases exhibit s- and d-wave characters simultaneously. One possible origin is s + id pairing, which breaks the point-group (PG) and time-reversal (TR) symmetries. While TR symmetry breaking may develop in special instances in the bulk³⁷ or on the surface³⁸, FeSCs typically preserve TRS. Especially for the alkaline Fe-selenides, there is no evidence for either TR symmetry-breaking or two-stage phase transitions as the temperature is lowered. Thus, it is important to identify an alternate candidate pairing. For the Fe-chalcogenide SC, one candidate pairing state was named sτ₃¹⁷. It has the s−waveform factor, but τ₃, a Pauli matrix in the xz, yz 3d-electron orbital basis, turns the pairing state into d − wave-like; indeed, in the band basis, the sτ₃ pairing has the intra- and inter-band d + d form. That both intra- and inter-band pairings can play a role is to be expected in this type of model³⁹ and other cases^4,40. However, this d + d form is highly unusual, thereby raising the question of both its naturalness and generality.

With the stage set by the above, the present work makes two advances. First, we demonstrate that the d + d pairing state belongs to matrix singlet pairing order parameters with non-trivial orbital structures that are natural and likely common-place in multi- orbital/band systems. As the orbital degrees-of-freedom (DOF) transform non-trivially under PG operations, these matrices can be chosen as one of the irreducible representations of the same group.

We make a case for the matrix singlet pairing’s naturalness by presenting an in-depth analysis of the sτ₃ pairing state. Written in the band basis, the sτ₃ pairing has the intra- and inter-band d + d form, but it remains an irreducible B_1g representation of the (tetragonal D_4h) PG. The unusual d + d pairing state is to be contrasted with its more commonly discussed d + id counterpart. Nonetheless, it is well defined. We demonstrate this point by showing that the d + d singlet pairing state can be compared and contrasted with the more familiar d + id state in analogy with how, in the case of superfluid ³He, the well-defined B-phase is measured against the equally well-known A-phase. The latter are spin-triplet pairing states that have an inherent matrix structure—in spin space—even for single-band cases.

Second, we illustrate the matrix singlet pairing’s generality by constructing this type of state in other multiband systems, for the case of heavy fermion superconductor CeCu₂Si₂. This is an important undertaking, given that CeCu₂Si₂ is the first-ever discovered unconventional SC⁴¹, and also recognizing that heavy fermion systems represent a prototype setting for strong correlations and unconventional superconductivity in general^42,43,44. Using DOF that incorporate spin-orbit coupling (SOC), we introduce an sΓ₃ state. This provides the theoretical basis for the excellent description of the experimental results in CeCu₂Si₂ in terms of the d + d pairing order parameter^24,34.

We will for the most part direct our analysis towards the effect of multiple orbitals/bands on the nature of the pairing states. Therefore, the issue of what drives such pairing states in the multi-orbital/band settings will only be briefly considered. Where this is done, our emphasis is on the short-range spin exchange interactions that are themselves induced by the underlying Coulomb (Hubbard and Hund’s) interactions. What we have achieved, from these microscopic calculations, is to demonstrate the relevance of the general considerations given above. We expect that our calculation will motivate further microscopic studies that include additional microscopic physics, such as orbital fluctuations, or are based on other approaches to the electron-correlation effects.

The emphasis of the present work is on singlet pairing states. Triplet pairing already has a matrix form that transforms nontrivially in spin space, even for single-band systems such as ³He. However, candidate solid state systems for triplet pairing often involve multiple orbitals/bands and strong correlations^{45,46,47,48,49,50}. Thus, the type of matrix pairing structure in the orbital/band space we consider here may produce triplet superconducting states^51,52 and the associated excitations that are of potential interest to quantum computing.

The remainder of the paper is organized as follows. We begin the first subsection of our results by discussing some of the most relevant general properties of non-trivial matrix pairing in the context of the Fe-based SCs. We subsequently define the sτ₃ pairing, and discuss the unusual properties of this state and show how it can be stabilized in an s- to d-wave degeneracy regime. We also support our discussion with numerical results for the pertinent five-orbital models of the Fe-based SCs. Furthermore, we consider sτ₃ in the band basis and illustrate how it is analogous to ³He -B. In the second subsection, we contrast the multiband d + d intra- and inter-band pairing state with the single-band d + id pairing state, and argue that these two cases are the analogs of ³He B and A. We show how they can be stabilized in a t − J₁ − J₂ model. In the third subsection, we extend the notion of non-trivial orbital structure beyond the Fe-based compounds by discussing a candidate analogous to sτ₃ for the heavy-fermion SC CeCu₂Si₂. In order to clearly see these results, we also present the irreducible representations of the D_4h point group in the context of CeCu₂Si₂. The Methods section contains additional accounts of the numerical results which support the stability of sτ₃ pairing for the alkaline Fe-selenides. We also discuss the t − J model and its solutions which illustrate the case of d + id pairing. Additional important aspects of matrix pairing are discussed in the Supplementary Notes. There, for completeness, we outline the role of the matrix-pairing functions in the various phases of superfluid ³He, where spin provides the analog of the orbital DOF. We highlight the lessons we believe can be applied to the case of multiorbital pairing in unconventional SCs. In addition, we illustrate how s- and d-wave states can coexist without breaking either PG or TR-symmetries via a general Landau-Ginzburg analysis. The band-basis representation of sτ₃ pairing and an illustration of the effects of damping on Bogoliubov-de Gennes (BdG) quasiparticles are also presented in the Supplementary Notes.

Results

d + d matrix singlet pairing as an analog of ³He -B

In solid-state systems, electrons inherit the orbital structure of the underlying ions which form the crystalline lattice. The set of local DOF must include the additional orbital structure. In turn, Cooper pairs formed out of the same electrons are naturally characterized by these additional local, orbital DOF.

Consider the concrete case of an electronic system on a lattice with D_4h tetragonal point-group symmetry. Further, assume that the dominant contribution to the lowest-lying bands is due to xz and yz orbital local DOF. For simplicity, we ignore SOC. In general, the pairing interactions $V{({\bf{k}},{\bf{k}}^{\prime} )}_{\alpha \beta ;\Gamma \delta }$ depend on the momenta as well as the orbital and spin indices of the two electrons. This two-dimensional space turns out to be relevant for Fe-based SCs ^53,54,55, and we will first define the sτ₃ pairing state in this space. The pairing is orbitally selective in that it is intraorbital and its amplitude is orbital dependent. We will then consider the stability of the matrix singlet pairing state in more realistic five-orbital models. Through the d + d representation in the band basis, we present an intriguing analogy of the singlet pairing state as an analog of ³He -B.

Matrix pairing in multiorbital Fe-based SCs

A spin-singlet pairing restricted in the orbital space to the xz, yz cubic harmonics will have the general form

$$\hat{{\Delta }}({\bf{k}})=\Delta \hat{g}{({\bf{k}})}_{\alpha \beta }i{\sigma }_{2}.$$

(1)

The even-parity matrix $\hat{g}({\bf{k}})$ denotes the components of the pairing in the four-dimensional space spanned by the tensor products of the two orbital DOF. These tensor-product states are analogs to the spin-1/2 product states in triplet ³He (see Supplementary Note 1). Likewise, they depend on the relative momentum of the pair. Finally, iσ₂ denotes spin-singlet pairing. We do not consider this additional structure since it plays no essential role in the subsequent discussion.

The pairing matrix can be decomposed into components which transform according to one of the five even-parity irreducible representations of the D_4h point group. This allows for additional separation of the DOF as

$${\hat{g}}^{(i)}{({\bf{k}})}_{\alpha \beta }={g}^{(i)}({\bf{k}}){\hat{\tau }}_{\alpha \beta }^{(i)},$$

(2)

where i labels one of the five, even-parity A_1g, A_2g, B_1g, B_2g, and E_g irreducible representations of D_4h. The functions g⁽ⁱ⁾(k) can likewise be chosen to belong to one of these representations. To illustrate, s-wave states such as ${s}_{{x}^{2}+{y}^{2}}({\bf{k}})$ and ${s}_{{x}^{2}{y}^{2}}({\bf{k}})$ both belong to A_1g. Standard d-wave states such as ${d}_{{x}^{2}-{y}^{2}}({\bf{k}})$ and d_xy(k) are B_1g and B_2g representations, respectively. The xz, yz orbital doublet transforms as the two-component E_g representation. The ${\tau }_{\alpha \beta }^{(i)}$ identity and Pauli matrices describe linear combinations of the tensor-product states which transform according to one of four irreducible representations contained in the E_g × E_g = A_1g + A_2g + B_1g + B_2g decomposition of the tensor-product space of the two E_g orbital DOF. By analogy to the total S = 1 spin states of ³He, these matrix-elements play the role of effective Clebsch-Gordan coefficients. The τ₀, τ₁, and τ₃ matrices transform according to A_1g, B_2g, and B_1g, respectively. In this work, we consider parity-even spin-singlet pairings belonging to one-dimensional irreducible representations of D_4h. This naturally excludes pairing states involving τ₂ matrices, which would be parity-odd.

These arguments point to an important aspect. In ³He, the relative angular momentum and local (spin) DOF transform independently under separate groups. In the present case, g⁽ⁱ⁾(k) and orbital matrix parts (τ⁽ⁱ⁾) are necessarily coupled since they both transform under the same PG. In effect, this constitutes an inherent SOC-like locking of the different spatial DOF of the Cooper pair.

We note that the single-component representation pairings in Eqs. (1) and (2) are unitary such that

$${\hat{\Delta }}^{\dagger }({\bf{k}})\hat{\Delta }({\bf{k}})={\Delta }^{2}{g}^{2}({\bf{k}}){\tau }_{0}.$$

(3)

Of particular relevance to our discussion is the fact that pairing with non-trivial matrix structure in general allows for several inequivalent representations of the PG. The problem of determining the stability of the different pairings, including those with non-trivial structure, is a challenging task, which is typically treated numerically on a case-by-case basis. We illustrate this point further below in this section, within a five-orbital t-J₁-J₂ model.

s τ ₃ pairing state

Of interest here is the sτ₃ pairing. In terms of Eqs. (1) and (2), it corresponds to

$$g({\bf{k}})={s}_{{x}^{2}{y}^{2}}({\bf{k}})$$

(4)

$${\hat{\tau }}_{\alpha \beta }={\tau }_{3,\alpha \beta }.$$

(5)

It transforms as the B_1g representation due exclusively to the τ₃ matrix. Because of the orbital struture, it represents neither simple s- nor d-wave states. However, sτ₃ pairing preserves both PG- and TR-symmetries of the normal state.

To illustrate the properties of the sτ₃ pairing, we first consider a simplified two-orbital model⁵³ and neglect any possible subleading channels. The TB and pairing parts of the BdG Hamiltonian in the orbital basis read³⁹

$${\hat{H}}_{\text{BdG}}({\bf{k}})={\hat{H}}_{\text{TB}}({\bf{k}})+{\hat{H}}_{\text{Pair}}({\bf{k}})$$

(6)

$${\hat{H}}_{\text{TB}}({\bf{k}})=\left[\left({\xi }_{0}({\bf{k}})-\mu \right){\tau }_{0}+{\xi }_{1}({\bf{k}}){\tau }_{1}+{\xi }_{3}({\bf{k}}){\tau }_{3}\right]\otimes {\gamma }_{3}$$

(7)

$${\hat{H}}_{\text{Pair}}({\bf{k}})=\Delta (k){s}_{{x}^{2}{y}^{2}}(\hat{{\boldsymbol{k}}}){\tau }_{3}\otimes {\gamma }_{1}.$$

(8)

The γ Pauli matrices act in Nambu space. To simplify the expressions, we discuss one of the two spin-sectors. With singlet pairing, the Hamiltonian for the other sector can be obtained in straightforward fashion. Note that, from the perspective of point-group symmetry classification, sτ₃ transforms in the same B_1g representation as the diagonal-in-orbital-space ${d}_{{x}^{2}-{y}^{2}}{\tau }_{0}$ pairing, as has been discussed in this type of model³⁹ and related settings^56,57,58,59. What distinguishes the sτ₃ pairing is the nontrivial commutation relation between the corresponding pairing and kinetic parts of the Hamiltonian¹⁷.

It is instructive to recognize that the ξ₁τ₁ and ξ₃τ₃ terms of ${\hat{H}}_{\text{TB}}$ play a role similar to a Rashba SOC. The bands corresponding to the normal-state dispersion are

$${\epsilon }_{\pm }({\bf{k}})={\xi }_{0}({\bf{k}})\pm \sqrt{{\xi }_{1}^{2}({\bf{k}})+{\xi }_{3}^{2}({\bf{k}})},$$

(9)

reflecting the space-group allowed varying orbital-content and splitting of the Fermi surfaces (FSs). We refer the reader to “Pairing channels of the five-orbital t − J₁ − J₂ model” subsection in Methods for detailed expressions of the ξ’s. The FS corresponding to this effective model has electron pockets centered at the ( ± π, 0) and (0, ± π) points of an one-Fe unit cell.

In Ref. ¹⁷, we showed that the general BdG dispersion is always gapped along the FS. Nodes away from the FS can appear for larger band splitting^17,60,61, reflecting the corrections to the gap term due to the non-commuting TB and pairing parts. However, in alkaline-Fe selenides, where sτ₃ pairing is competitive, the small band splitting at the center of the Brillouin zone precludes the appearance of nodes. Even in the cases when the nodes were to appear in the BdG spectrum away from the Fermi surface, it will not affect our conclusion. The point is that, in strongly correlated systems, only nodal excitations on the Fermi energy are long lived and, thus, sharply defined. For states away from the Fermi energy, any putative nodal excitations will necessarily involve a large damping caused by the underlying electron correlations, which obviates the distinction between nodal and gapped excitations. We illustrate how this can occur in Supplementary Note 5.

Another important characteristic of such a gapped sτ₃ state is its sign change under a π/2 rotation. Such a sign-change leads to the formation of an in-gap spin resonance. sτ₃ is then a pairing state which reconciles a fully-gapped FS with the presence of a spin-resonance, typically associated with a d − wave gapless order parameter.

Although we focus on a simplified two-orbital model in order to illustrate the salient properties of sτ₃ pairing, the latter can also be stabilized in more general five-orbital models of the alkaline Fe-selenide class of SCs. The pairing matrix in the t − J₁ − J₂ model can be decomposed into all the symmetry-allowed channels. The complex coefficients of these components have both amplitude and phase. The illustrative results are shown in Fig. 1a, b. The zero-temperature pairing amplitudes of all symmetry-allowed pairing channels have been determined in a five-orbital t − J₁ − J₂ model with nearest and next-nearest exchange couplings. This model and its solution method are discussed in subsection “Five-orbital t − J₁ − J₂ model and solution method” in Methods section.

Fig. 1: Zero-temperature results for a five-orbital t − J₁ − J₂ model.

a Dimensionless pairing amplitudes of the leading B_1g channels as compared to that of the B_1g sτ₃ channel as functions of J₁/J₂ for a five-orbital t − J₁ − J₂ model of the alkaline Fe-selenides. The numerically-determined pairing amplitudes are the weights of each of the PG symmetry-allowed channels. 1_xy denotes the trivial 1 × 1 matrix in the d_xy orbital sector. See subsections “Pairing channels of the five-orbital t − J₁ − J₂ model” and “Five-orbital t − J₁ − J₂ model and solution method” in Methods for the details of the calculation. The sτ₃ pairing with non-trivial orbital structure is dominant in the 0.8 ≤ J₁/J₂ ≤ 1.0 window. b Phases of the leading B_1g channels relative to the sτ₃ channel as functions of the tuning parameter. These are obtained from the difference in the phases of each symmetry-allowed channel which are determined from the self-consistent solution. In the [0.8, 1] interval where sτ₃ is dominant, these relative phases are either zero or ± π. Here, the amplitudes of the coexisting B_1g channels are comparable to that of sτ₃. This illustrates that the sτ₃ pairing which is equivalent to d + d, effectively preserves TR and PG symmetries.

The TB part and the associated FS are chosen to be consistent with LDA studies¹⁷. The dominant pairing amplitudes are intra-orbital. The pairing state is orbital selective in the sense that the pairing amplitude and its phase are orbital sensitive. We focus on the case where the pairing amplitude is largest for the xz, yz orbitals while also allowing inter-orbital pairing. This reflects orbital-selective correlation effects in the normal state. The J₁/J₂ ratio controls the symmetry of the dominant pairing channel with ${s}_{{x}^{2}{y}^{2}}({\bf{k}}){\tau }_{0}$ and ${d}_{{x}^{2}-{y}^{2}}({\bf{k}}){\tau }_{0}$ states for small and large values of the ratio, respectively. The sτ₃ pairing is dominant near the transition separating order-parameters belonging to A_1g to B_1g representations for a finite range of the control parameter about the point where J₁/J₂ ≈ 1. A ${d}_{{x}^{2}-{y}^{2}}{\tau }_{0}$ with trivial orbital structure provides the subleading pairing with comparable amplitude. See subsection “Pairing channels of the five-orbital t − J₁ − J₂ model” of the Methods for more details.

It is important to put the results of microscopic studies in a more general perspective. Our calculations indicate that a subleading ${d}_{{x}^{2}-{y}^{2}}{\tau }_{0}$ pairing of comparable amplitude is present in the regime where sτ₃ is dominant. While we have focused on the properties of the dominant sτ₃ pairing alone, a more realistic picture would involve coexisting sτ₃ and ${d}_{{x}^{2}-{y}^{2}}{\tau }_{0}$ in the vicinity of s- to d-wave phase transition. This superposition of pairing states with different orbital structure which belong to the same B_1g irreducible representation preserves both PG- and TR-symmetries. In Supplementary Note 3, we present a Landau-Ginzburg analysis to show that, generically, the pairing state involves a linear superposition of these two components and there is only one superconducting transition at a single T_c.

Intra- and inter-band d + d pairing and its analogy with the ³He B-phase

It is instructive to consider the sτ₃ pairing in a band basis:

$${\hat{H}}_{\text{Pair}}({\bf{k}})={\Delta }_{3}({\bf{k}}){\alpha }_{3}+{\Delta }_{1}({\bf{k}}){\alpha }_{1},$$

(10)

where α_1,3 are Pauli matrices in the two-band space and where the form factors Δ_3,1 transform as ${d}_{{x}^{2}-{y}^{2}}$ and d_xy, respectively. The details of the transformation from orbital to band basis are discussed in Supplementary Note 4. There, we also show that the α_3,1 matrices are equivalent to A_1g and A_2g representations, respectively, by applying the inverse transformation from band to orbital basis. The same conclusion can be reached by requiring that each of the two terms in Eq. (10) transforms as B_1g. Because the overall pairing is in the irreducible B_1g channel, it is natural that the intraband α₃ part has the ${d}_{{x}^{2}-{y}^{2}}$ form factor. Likewise, the interband α₁ component has the d_xy form factor. Thus, the sτ₃ pairing is equivalent to a d + d pairing.

When the pairing matrix is squared, the intra- and inter-band d-waves add in quadrature as ${\Delta }_{1}^{2}({\bf{k}})+{\Delta }_{2}^{2}({\bf{k}})$ to produce a gap which does not close along the FSs centered on the BZ center corresponding to the two bands, as shown in Fig. 2. This is due to the anti-commuting nature of the two Pauli matrices α₃ and α₁ which denote intra- and inter-band pairing, respectively. As in the orbital basis, corrections to this gap are present due to the splitting of the FSs. As discussed in the previous subsection, in cases relevant to our discussion, these additional effects are typically small and consequently do not close the gap; and generically, the FS is always fully gapped.

Fig. 2: The gapping of an illustrative FS by d + d pairing.

The blue, dashed line indicates the nodes of the intra-band component Δ₃ (Eq. (10)), which transforms as a B_1g representation of D_4h. The red, dotted line shows the nodes of the inter-band component Δ₁ (also in Eq. (10)), which transforms as a B_2g representation. The two components add in quadrature to produce a nonzero gap everywhere on the FS. Note that possible nodes of a common s-wave form factor (Eqs. 28-29 of Supplementary Note 4 for sτ₃) are not shown here as they are irrelevant to our argument.

The band basis reveals a pairing structure which is very similar to that in ³He -B. Referring to Supplementary Notes 1 and 2, the matrix order-parameter in that case is typically expressed as ${\hat{\Delta }}_{{\,}^{3}\text{He - B}}({\bf{k}}) \sim ({\bf{k}}\cdot {\boldsymbol{\sigma }})i{\sigma }_{2}$. This amounts to a linear superposition of p-wave states, p_x, p_y and p_z, together with a matrix structure made possible due to spin-triplet pairing as represented by the σ Pauli matrices. The anti-commuting nature of these matrices ensures that three p-waves add in quadrature to produce a full gap. The situation clearly mirrors the case of sτ₃ in the band basis, where two d-wave states likewise produce a finite gap. The sτ₃ pairing thus provides a remarkable example where a phase which is similar to ³He -B via a structure in the band-basis is stabilized in a solid-state SC model.

Along with this similarity between the sτ₃ pairing state and the B phase of ³He, it is important to also note on the ways in which they differ. The distinctions are due primarily to the continuous rotation symmetries of ³He as contrasted with the discrete nature of the PG in the inter- and intra-band d-wave case. The latter belong to a single irreducible representation of a PG involving discrete operations. As such, they break no symmetries of the normal state with the trivial exception of a global phase rotation due to pairing. By contrast, ³He -B breaks the SO(3)_L × SO(3)_S symmetry of the normal state down to SO(3)_L+S, via a relative spin-orbit symmetry breaking^62,63. Specifically, the invariance of the normal state under continuous and independent rotations of angular-momentum and spin, respectively, is broken down to simultaneous rotations in both sectors. In spite of this additional symmetry-breaking, we note that the B phase has the largest residual symmetry of all superfluid ³He phases. In this respect, it still resembles intra- and inter-band d-wave pairing which preserves both PG and TR symmetries.

d + d and d + id pairing: Analogy with ³He -B vs.³He -A

We have seen that an orbital basis is convenient for classifying the pairing states according to symmetry and for solving microscopic models. Physically, the equivalent band basis is more natural in connecting with experiments. We have also seen that the non-trivial sτ₃ pairing is equivalent to simultaneous intra- and inter-band ${d}_{{x}^{2}-{y}^{2}}$ and d_xy pairings (Eq. (10)). These add in quadrature to produce a full gap on the FS on either of the two bands and their sign-changing factors allow for the formation of in-gap spin resonances. For simplicity, we consider only unitary pairings. The intra- and inter-band terms are consequently associated with α₃ and α₁ Pauli matrices, respectively. Importantly, a d + d pairing does not break either PG or TR symmetries of the Hamiltonian. This amounts to associating both d-wave components with a single irreducible representation in an orbital basis.

We have shown that the d + d pairing is a well-defined pairing state, through an analogy with the B phase of ³He. To further elucidate the naturalness of this unusual pairing state, we compare and contrast it with the more familiar d + id pairing. We show that d + d vs. d + id pairing is analogous to the B- vs. A-phases of ³He.

d + d in a multiorbital model vs. d + id in a single-orbital model

An intra-band d + id pairing, where the two components are ${d}_{{x}^{2}-{y}^{2}}$- and d_xy-waves, respectively, is a natural competitor to the intra- and inter-band d + d. Here, we show how the intra-band d + id can be stabilized in a one-band t − J₁ − J₂ model in the vicinity of the J₁ ≈ J₂ point.

In the first subsection, we discussed how the sτ₃ orbital non-trivial pairing channel becomes dominant for a finite range of the J₁/J₂ tuning parameter in a realistic five-orbital t − J₁ − J₂ model for the alkaline Fe-selenides. The details of the calculations are given in the Methods section. We showed how sτ₃ pairing is equivalent to a d + d intra- and inter-band pairing. To further understand the nature of the sτ₃ -dominated state, we plot the phases of the leading B_1g channels relative to sτ₃ as functions of J₁/J₂ in Fig. 1 (b). The leading B_1g channels have relative phases wrt sτ₃ which are closely centered around 0 or π for the entire range of the tuning parameter. In the [0.8, 1] interval where sτ₃ is dominant the relative phases are either zero or ± π. Here the amplitudes of the subleading ${d}_{{x}^{2}-{y}^{2}}{\tau }_{0}$ and ${s}_{{x}^{2}+{y}^{2}}{\tau }_{3}$B_1g channels are comparable to that of the leading sτ₃. Therefore, this regime corresponds to a pairing state with non-trivial orbital structure which preserves TR and PG symmetries. We note that all A_1g and B_2g channels are strongly suppressed in the regime where sτ₃ is dominant which we consider here (Fig. 3).

Fig. 3: Zero-temperature leading pairing amplitudes (dimensionless) as functions of ${J}_{1}^{xz/yz}/{J}_{2}^{xz/yz}$ for a five-orbital t − J₁ − J₂ model of the alkaline Fe-selenides.

${J}_{2}^{xz}={J}_{2}^{yz}=1/2$ in units of the bandwidth. The exchange interactions for the d_xy orbital are five times smaller while the exchange couplings for all remaining orbitals are zero. Please see Ref. ¹⁷ for a detailed account of the model and solution. For J₁/J₂ ≤ 0.7 A_1g pairing channels are dominant with leading ${s}_{{x}^{2}{y}^{2}}{\tau }_{0}$ in the xz/yz sector. This pairing has trivial orbital structure. There is a narrow region of coexistence between finite A_1g and B_1g channels for 0.7 ≤ J₁/J₂ ≤ 0.8. Beyond this range, B_1g channels dominate with leading ${s}_{{x}^{2}{y}^{2}}{\tau }_{3}$ in the xz/yz sector which has non-trivial orbital structure. At even larger values, an orbital-trivial ${d}_{{x}^{2}-{y}^{2}}{\tau }_{0}$ phase dominates.

We next discuss the orbital-trivial d + id pairing. For simplicity, we consider a single d_xy orbital t − J₁ − J₂ model on a square lattice. We choose the tight-binding (TB) parameters and chemical potential to be consistent with a circular hole pocket at the center of the BZ. The details of the model are discussed in subsection “Single-orbital t − J₁ − J₂ model” of the Methods. The model is solved using a self-consistent decomposition of the exchange interactions as in the five-orbital cases discussed previously^17,64. The resulting zero-temperature pairing amplitudes for J₂ = 1/2 in units of the bandwidth, and for a finite range of the ratio J₁/J₂ are shown in Fig. 4a. For J₁/J₂ < 0.8, the only significant pairing occurs in the d_xy, B_2g channel. For higher values of J₁/J₂, the amplitude of a ${d}_{{x}^{2}-{y}^{2}},{B}_{1g}$ pairing becomes finite. These two remain finite up to J₂/J₁ ≈ 2.1, where the d_xy component vanishes continuously. Beyond this point, two additional ${s}_{{x}^{2}+{y}^{2}}$ and ${s}_{{x}^{2}{y}^{2}}$ order-parameters emerge. A similar conclusion has been reached in a related model⁶⁵.

Fig. 4: Zero-temperature results for a single-orbital t − J₁ − J₂ model close to half-filling.

See subsection “Single-orbital t − J₁ − J₂ model” of the Methods for details of the model. (a) Dimensionless pairing amplitudes as functions of the ratio J₁/J₂. When the tuning parameter is less than 0.8, only the d_xy, B_2g channel has finite amplitude. In the 0.8 ≤ J₁/J₂ ≤ 2.1 interval, d_xy coexists with a ${d}_{{x}^{2}-{y}^{2}},{B}_{1g}$ channel. For larger values of the tuning parameter, the d_xy channel is suppressed and ${s}_{{x}^{2}+{y}^{2}}$ and ${s}_{{x}^{2}{y}^{2}}$A_1g channels emerge. (b) Relative phase of the two d-wave channels modulo π as a function of J₁/J₂. When both d-waves have finite amplitudes, a π/2 relative phase is clearly visible. When one of the two is suppressed, the relative phase is essentially arbitrary.

To illustrate that the two coexisting d-wave components are locked into a d + id state, we plot their relative phases mod π in Fig. 4b. The relative phases are obtained from the difference in the phases of each symmetry-allowed channel which are determined from the self-consistent solution. While these relative phases are essentially arbitrary whenever one of the d-waves vanishes, a π/2 relative phase is clear in the interval J₁/J₂ ∈ [0.8, 2.1] where both coexist. Although these results were obtained for a single-orbital model, they do demonstrate how a d + id pairing with trivial orbital structure can become stable in similar two-orbital models.

d + d pairing vs. d + id: Analogy with B- vs. A-phases of ³He superfluid

In the first subsection, we showed that the d + d pairing is closely analogous to the B phase of ³He, where the pairing is a superposition of p_x,y,z-waves corresponding to equal- and opposite-spin pairing as illustrated by Eqs. 2 and 3 of Supplementary Notes 1 and 2. Just like the B phase, the d + d pairing is an irreducible representation, here of the PG, and preserves the TR symmetry of the normal state by construction.

By contrast, the intra-band (d + id)α₀ pairing, where α₀ is the identity matrix in the band basis, is a linear superposition of two irreducible representations.

In general, the α₀ matrix in band space would correspond to an identity τ₀ matrix in orbital space. The d + id pairing spontaneously breaks both PG and TR symmetries. Therefore, it is a natural analog of ³He -A, which is typically described in terms of equal-spin p_x + ip_y pairing. This phase spontaneously breaks both rotational and TR-symmetries of the normal state as illustrated by Eq. 5 of Supplementary Note 2. The band and spin matrices in the ³He -A and d + id cases are analogous.

In the first subsection, we discussed how d + d differs from ³He -B due mainly to the discrete versus continuous symmetries of the two, respectively. This kind of difference also exists between d + id and ³He -A. The latter spontaneously breaks both angular momentum and spin continuous rotational symmetries down to a ${U}_{{L}_{z}-\phi }\times {U}_{{S}_{z}}$ group of independent rotations in each sector about preferred axes (Supplementary Note 2). When dipole-dipole interactions are negligible, the directions of either axes are arbitrarily chosen. By contrast, d + id involves a superposition of pairings belonging to two irreducible representations of a discrete PG corresponding to fixed symmetry axes. Additionally, in ³He -A, the two components, p_x and p_y are exactly degenerate, and there is only a single transition temperature T_c. By contrast, in d + id, the two d-components are generically non-degenerate, and two stages of phase transitions are to be expected when the temperature is varied.

In spite of clear differences, the formal similarities between ³He -B and d + d and likewise between ³He -A and d + id, which are due to the presence of non-trivial matrix structure, are intriguing. In this sense, ³He provides both a well-established parallel and a prototype for the emergence and description of the effects of non-trivial matrix structure in unconventional singlet SCs.

Given the venerable status of superfluid ³He -B, we believe that revealing the above connections elevates the status of the d + d spin-singlet pairing state. In turn, this connection motivates the consideration of such d + d spin-singlet pairing beyond the context of Fe-based superconductors. Indeed, this leads us to the second part of our work, which is to propose a microscopic pairing state that is capable of understanding the heavy fermion superconductor CeCu₂Si₂.

Matrix singlet pairing with spin-orbit coupling: CeCu₂Si₂

Another class of multiband superconductors arises in heavy fermion systems, in which quasi-localized f electrons hybridize with dispersive spd-conduction (c) electrons. These include CeCu₂Si₂, which is the first-ever discovered unconventional SC⁴¹ and one of the best-studied heavy-fermion SCs. For most of its history, this compound was believed to have a conventional d-wave order parameter. Such a conclusion has been supported by inelastic neutron scattering experiments which revealed a spin-resonance peak in the SC state²⁸ together with angle-resolved resistivity measurements of the upper critical field H_c2⁶⁶, among others. Remarkably, recent measurements of the specific heat³⁵ and London penetration depth^34,36 down to lower temperatures indicated a fully-gapped SC state. The apparent contradiction between these different experimental probes is reminiscent of the situation in the Fe-chalcogenide SCs. In those cases, we argued that the fully-gapped but sign-changing sτ₃ provide a natural resolution. An analogous proposal for CeCu₂Si₂ is clearly of great interest. Note that a d + d inter- and intra-band pairing directly inspired from the Fe-based cases provides a good fit to the the superfluid-density and specific-heat results in CeCu₂Si₂^24,34.

Objective and outline of the subsection

Here, we construct the analog of the sτ₃ pairing for CeCu₂Si₂. For reasons that will become clear, we shall refer to this state as sΓ₃ to indicate the associated non-trivial pairing matrix, as in the case of the Fe-based SCs.

We consider the pairing between the composite heavy quasiparticles in terms of simultaneous f-f, f-c, and c-c pairing in the original electron basis prior to hybridization. Of the three, f-f pairing is expected to be the strongest, reflecting the more localized nature of the heavy bands. The albeit weaker f-c and c-c pairing terms will be important, especially when they are involved in creating a pairing component that opens a gap.

In contrast to the case of the Fe-based SC, an important ingredient for constructing pairing states in a heavy fermion metal such as CeCu₂Si₂ is that SOC plays a 0th-order role. The local orbital and spin DOFs transform simultaneously under PG operations. This imposes additional constraints on any matrix associated with the local DOF. Due to the large SOC, the local f-electron manifold splits as a consequence of the crystal field. The resulting multiplets, which are labeled according to the irreducible representations of D_4h, play a role analogous to that of the d_xz/yz orbitals in the Fe-based SC case.

A number of experiments^67,68 as well as LDA+DMFT studies⁶⁹ have indicated that one of the Γ₇ doublets of the crystal-field split ²F_5/2 local electron is the dominant contribution to the heavy FS sheets. The lowest-lying excited states of the f electron are composed of a Γ₆ and another Γ₇ doublet. Our analysis will also involve Γ₆ Wannier orbitals of the conduction electron states near the Fermi energy, and these Wannier orbitals will hybridize with the excited Γ₆ f-level and thereby makes it a small but nonzero component in the ground-state manifold.

In this subsection, we will use these DOFs to advance a matrix pairing state, sΓ₃, which transforms in B_1g under D_4h. To clarify the involved DOFs, we also discuss the character table of the D_4h point group and its irreducible representations and construct the conduction-electron Γ₆ Wannier orbitals from the Cu 3d orbitals.

Spin-orbital coupled local states

Our aim is to propose a minimal symmetry-allowed candidate for CeCu₂Si₂ which has properties similar to that of sτ₃ in the Fe-based cases. By construction, such a state must belong to one of the single-component double-valued irreducible representations of D_4h, as required by strong SOC. To illustrate, the even-parity double-valued irreducible representation ${\Gamma }_{1}^{+}$ is the analog of A_1g, while ${\Gamma }_{3}^{+}$ is the analog of B_1g. The latter is our prime candidate.

Either f or c electrons originate from an odd-spin state and therefore transform as either Γ₆ or Γ₇ representations of the PG. A minimal structure for combined local orbital-spin DOF is determined by a 2 × 2 matrix Σ. This matrix must belong to a non-trivial irreducible representation of the PG; e.g., it must change sign under a C_4z rotation. To ensure that the rotation properties are determined exclusively by the local orbital-spin DOF, the pairing must be a product between the non-trivial orbital-spin matrix and a form-factor belonging to the identity representation. In addition to the matrix structure of the local DOF, the pairing matrix must also incorporate c, f indices.

Thus, a minimal order-parameter is a 4 × 4 matrix. We consider singlet, parity-even pairing exclusively. Hence, candidate pairing matrices must be odd under exchange and TR-invariant. For simplicity, we restrict our discussion to pairings which are even under f-c exchange. This necessarily implies that Σ is anti-symmetric. Furthermore, as the Σ matrix can transform under inversion, we only consider pairings between electrons belonging to irreducible representations of identical parity. Following the notation used previously, possible candidates are chosen to be of the form

$$\hat{\Delta }({\bf{k}})=\Delta g({\bf{k}})\hat{\Sigma }\otimes \hat{\Xi }.$$

(11)

The components of the local-DOF multiplet are determined by the 2 × 2 Σ matrix while the f, c nature of the paired electrons is given by the 2 × 2 Ξ matrix. As in the more familiar case of full spin rotational symmetry, the matrix elements of the Σ matrices are effective Clebsch–Gordan coefficients adapted to the cases of discrete PG symmetry⁷⁰.

Conventional B _1g pairing

We first consider candidates on the ${\Gamma }_{7}^{-}$ ground-state doublet. The superscript denotes odd parity under inversion. Although these naturally correspond to f-f pairing involving the Γ₇ ground-state local multiplet, they also cover possible f-c pairings with c electrons which belong to the same representation. In the latter case, the c electrons would correspond to p-type linear-superposition of Wannier orbitals. The tensor product of two such doublets decomposes into the irreducible representations of D_4h as⁷⁰

$${\Gamma }_{7}^{-}\times {\Gamma }_{7}^{-}={\Gamma }_{1}^{+}+{\Gamma }_{2}^{+}+{\Gamma }_{5}^{+}.$$

(12)

Here, ${\Gamma }_{1,2}^{+}$ are one-dimensional representations which are analogous to the A_1g and B_2g in the absence of SOC⁷⁰. The two-dimensional ${\Gamma }_{5}^{+}$ is analogous to the xz, yz (E_g) doublet. Following Ref. ⁷⁰, the matrices corresponding to each of the three ${\Gamma }_{1,2,5}^{+}$ representations are:

$${{{\Sigma }}}_{{\Gamma }_{1}}=\frac{\,\text{i}\,}{\sqrt{2}}{\sigma }_{2}$$

(13)

$${{{\Sigma }}}_{{\Gamma }_{2}}=-\frac{\,\text{i}\,}{\sqrt{2}}{\sigma }_{1}$$

(14)

$${{{\Sigma }}}_{{\Gamma }_{5},x}^{(5)}=\frac{\,\text{i}\,}{\sqrt{2}}{\sigma }_{3}$$

(15)

$${{{\Sigma }}}_{{\Gamma }_{5},y}^{(5)}=\frac{1}{\sqrt{2}}{\sigma }_{0}.$$

(16)

The σ’s are standard Pauli matrices. Recall that for f-f or symmetric f-c pairings, we require Σ to be anti-symmetric. The only choice is ${\Sigma }_{{\Gamma }_{1}}$ which transforms as the trivial representation. This matrix is the analog of simple singlet-pairing in the standard BCS case and is invariant under all PG operations. It is clear that f-f or symmetric f-c singlet pairing on the ${\Gamma }_{7}^{-}$ manifold does not support any non-trivial structure in the local DOF. This contrasts with the Fe-based case, where the absence of SOC allowed for all τ matrices in the xz, yz manifold.

We can construct a standard d-wave pairing belonging to the Γ₃ representation which is analogous to a B_1g representation without SOC. We do so by choosing $g({\bf{k}})={d}_{{x}^{2}-{y}^{2}}({\bf{k}})$ and $\hat{\Sigma }={\Sigma }_{{\Gamma }_{1}}$ in Eq. (11). Likewise, $\hat{{{\Xi }}}$ can be chosen to be proportional to either Ξ₁ or (1/2)(Ξ₀ − Ξ₃), where Ξ₀ and Ξ_1,3 denote identity and Pauli matrices, respectively. The two cases correspond to f-c and f-f pairing, respectively.

Matrix B _1g pairing

We next consider pairing between electrons belonging to distinct ${\Gamma }_{7}^{-}$ and ${\Gamma }_{6}^{-}$ manifolds. This can correspond to f-f pairing between electrons belonging to the ${\Gamma }_{7}^{-}$ ground-state and f electrons belonging to the excited ${\Gamma }_{6}^{-}$ manifolds, respectively. Alternately, it can denote f-c pairing between ${\Gamma }_{7}^{-}$f-electrons and ${\Gamma }_{6}^{-}$ conduction c electrons. Further below, we illustrate how intra-unit cell linear combinations of Cu 3d states in the presence of SOC can form bases for ${\Gamma }_{6}^{-}$ conduction electrons. The product states decompose as⁷⁰

$${\Gamma }_{6}^{-}\times {\Gamma }_{7}^{-}={\Gamma }_{3}^{+}+{\Gamma }_{4}^{+}+{\Gamma }_{5}^{+}.$$

(17)

The corresponding matrices are⁷⁰

$${\Sigma }_{{\Gamma }_{3}}=\frac{\,\text{i}\,}{\sqrt{2}}{\chi }_{2}$$

(18)

$${\Sigma }_{{\Gamma }_{2}}=-\frac{\,\text{i}\,}{\sqrt{2}}{\chi }_{1}$$

(19)

$${\Sigma }_{{\Gamma }_{5},x}=\frac{1}{\sqrt{2}}{\chi }_{0}$$

(20)

$${\Sigma }_{{\Gamma }_{5},y}=\frac{\,\text{i}\,}{\sqrt{2}}{\chi }_{3}.$$

(21)

The χ’s are Pauli matrices. Note however that they represent different DOF and thus transform differently under the PG. Therefore, one should not confuse the meaning of the χ Pauli matrices defined in this case with those of the Γ₇ − Γ₇ case discussed previously. The only anti-symmetric matrix is ${{{\Sigma }}}_{{\Gamma }_{3}}$. It transforms as the ${\Gamma }_{3}^{+}$ irrep of D_4h and is an analog of the τ₃ matrix in the Fe-based cases. Moreover, this matrix is invariant under TR. We conclude that a counterpart of the sτ₃ pairing for CeCu₂Si₂ is an s-wave pairing belonging to the sign-changing Γ₃ representation, or sΓ₃ pairing:

$$\hat{\Delta }({\bf{k}})=\Delta s({\bf{k}})i{\chi }_{2}\otimes \hat{\Xi },$$

(22)

where s(k) corresponds to a s-wave form factor which transforms according to the ${\Gamma }_{1}^{+}$ trivial irrep. sΓ₃ pairing, which involves electrons belonging to different irreducible representations due to the Γ₃ matrix, is necessarily non-local, and thus vanishes when r_Relative = 0, where r_Relative is the distance between two paired electrons. Therefore, we do not restrict the s-wave form factor to be of sign-changing form. The form of the Ξ matrix differs depending on either f-f or f-c pairings. In the f-c case it can be chosen to be proportional to a Ξ₁ Pauli matrix. In the f − f case, it can be made proportional to a Ξ₀ − Ξ₃ matrix. In either case, the gap is determined by the amplitude and form factor only similarly to what happens for sτ₃. In a multiband model of CeCu₂Si₂³⁴, this pairing produces a full gap.

On general grounds, the non-trivial Γ₇-Γ₆ pairing in either f-f or f-c cases is likely weaker than the Γ₇-Γ₇ f-f pairing. However, there are cases where such Γ₇-Γ₆ contributions can be important. Consider a dominant Γ₇-Γ₇ f-f pairing corresponding to a d-wave state with nodes along the FS. An admixture of non-trivial pairing either from f-c or from f-f involving the excited local manifold can open a gap. While we can also consider non-trivial pairing terms in the c-c sector, these are expected to be weaker than their f-f and f-c counterparts. Likewise, other candidates with non-trivial orbital-spin structure can be obtained if we relax some of our assumptions such as the symmetry of the f-c terms under exchange. We reserve a detailed analysis of these cases for future work.

Our candidate Γ₇-Γ₆, sΓ₃ pairing represents an sτ₃ analog for CeCu₂Si₂. As in the Fe-based cases, the structure of the local DOF allows a natural interpolation between simple s − and d − wave states. Such a pairing can in principle reconcile the difficulties in interpreting the more recent experimental results.

We note that the ${\Gamma }_{6}^{-}$ conduction electrons which enter the matrix B_1g pairing likely originate from Cu 3d orbitals (see below). Indeed, several experiments^24,71,72 have indicated that the strongest suppression of T_c occurs upon substituting Cu by non-magnetic impurities. Our matrix B_1g pairing candidate, which involves ${\Gamma }_{6}^{-}$ conduction electrons from 3d Cu states, is naturally consistent with these findings.

Similar to the Fe-chalcogenide case, for unitary pairing we expect the sΓ₃ pairing in the band basis to contain the intraband α₃ and interband α₁ components. Each must be in d-wave state, with the form factor of the intraband α₃ being ${d}_{{x}^{2}-{y}^{2}}$. Thus, the sΓ₃ pairing realizes a d + d multiband pairing. Importantly, the d + d pairing does not break either PG or TR symmetries of the Hamiltonian.

As discussed in the first subsection, we expect that the sΓ₃ matrix-pairing will coexist with a conventional d-wave pairing below T_c since they both belong to the same Γ₃ irreducible representation of D_4h. The admixture between these will ensure that the SC state preserves both PG and TR-symmetries but also exhibits a gap which is finite everywhere along the FS.

We stress that our analysis is distinguished from the well-known symmetry-based procedure typically considered in the context of heavy-fermion SCs⁷³. The latter does not explicitly treat possible non-trivial matrix structures associated with the local DOF. Instead, the order-parameters are generically classified according to the irreducible representations of the various PGs in the context of a LG analysis. In our case, we have focused on the non-trivial role of the local DOF.

Irreducible representations of D _4h

To expound on the local DOFs, we now turn to the character table for the double-valued representations of the tetragonal D_4h point group, showing it in Table 1.

Table 1 Character table for the double-valued representations of the tetragonal D_4h point group.

We follow the conventions of Ref. ⁷⁰. E represents the identity, C_n are rotations about z by 2π/n, while ${C}_{2}^{\prime}$ and ${C}_{2}^{^{\prime\prime} }$ are π rotations about x, y and (−x, y), (x, y) axes, respectively. I denotes inversion. S₄ indicates a C₄ rotation about z followed by a reflection in the xy plane perpendicular to this axis. σ_h, σ_v, and σ_d are reflections through the xy, xz and yz, and diagonal-z planes, respectively.

In the presence of SOC, the double-valued, even-parity irreducible representations ${\Gamma }_{1-5}^{+}$ correspond to the single-valued A_1g, A_2g, B_1g, B_2g, and E_g representations, respectively in the absence of SOC. Likewise, the odd-parity ${\Gamma }_{1-5}^{-}$ correspond to the A_1u, A_2u, B_1u, B_2u, and E_u representations. ${\Gamma }_{6}^{+}$ and ${\Gamma }_{7}^{+}$ transform as spinors and product of spinors and linear combinations of ${\Gamma }_{1-4}^{+}$ and similarly for the odd-parity ${\Gamma }_{6}^{-}$ and ${\Gamma }_{7}^{-}$.

Typical bases for representations relevant to this work are ${\Gamma }_{3}^{+}:{x}^{2}-{y}^{2}$, ${\Gamma }_{4}^{+}:xy$, ${\Gamma }_{5}^{+}:xz,yz$. For odd-parity representations we mention ${\Gamma }_{1}^{-}:({x}^{2}-{y}^{2})xyz$, ${\Gamma }_{3}^{-}:xyz$, ${\Gamma }_{6}^{-}:{\Gamma }_{6}^{+}\times {\Gamma }_{1}^{-}$ and ${\Gamma }_{7}^{-}:{\Gamma }_{7}^{+}\times {\Gamma }_{3}^{-}$.

${\Gamma }_{6}^{-}$ conduction electrons

Previously, we discussed matrix B_1g pairing between ${\Gamma }_{6}^{-}$c -electrons and ${\Gamma }_{7}^{-}$ f-electrons. In this subsection, we illustrate how linear combinations of intra-unit-cell Cu 3d orbitals in the presence of SOC provide ${\Gamma }_{6}^{-}$ conduction electron states. We consider crystal field-split ${d}_{{x}^{2}-{y}^{2}}$ Cu orbitals for simplicity although similar constructions are possible for other orbitals.

Consider one of the two Cu planes in the unit cell of CeCu₂Si₂⁶⁹. The Cu sites are located halfway along the edges of a plaquette as illustrated in Fig. 5.

Fig. 5: Single Cu plane in the unit cell of CeCu₂Si₂.

The four sites labeled (1) − (4) correspond to Cu ${d}_{{x}^{2}-{y}^{2}}$ orbitals in the plane. The dashed-line circles represent the Ce sites projected onto the Cu-plane.

The linear combinations of the four orbitals at the Cu sites

$${p}_{x}={d}_{{x}^{2}-{y}^{2}}^{(4)}-{d}_{{x}^{2}-{y}^{2}}^{(2)}$$

(23)

$${p}_{y}={d}_{{x}^{2}-{y}^{2}}^{(1)}-{d}_{{x}^{2}-{y}^{2}}^{(3)}$$

(24)

transform as an (x, y) doublet under the D_4h point group. Consequently, (p_x, p_y) belong to a two-component ${\Gamma }_{5}^{-}$ irreducible representation. We further include the local spin-1/2 DOF which belongs to a ${\Gamma }_{6}^{+}$ irreducible representation⁷⁰. From the direct-product states of p-orbital linear combinations and spinor states ϕ_±1/2 we can construct states which belong to ${\Gamma }_{6}^{-}$ representation⁷⁰:

$${{{\Psi }}}_{{\Gamma }_{6}^{-};1/2}=\frac{\,\text{i}\,}{\sqrt{2}}[{p}_{x}+\,\text{i}\,{p}_{y}]{\phi }_{-1/2}$$

(25)

$${{{\Psi }}}_{{\Gamma }_{6}^{-};-1/2}=\frac{\,\text{i}\,}{\sqrt{2}}[{p}_{x}-\,\text{i}\,{p}_{y}]{\phi }_{1/2},$$

(26)

where SOC was taken into account. Similar states can be constructed in the remaining Cu plane. The symmetric linear combination between ${\Gamma }_{6}^{-}$ states in both Cu planes likewise belongs to ${\Gamma }_{6}^{-}$ doublet.

Discussion

Recent experiments in multiband Fe-based and heavy-fermion SCs are inconsistent with either simple s- or d-wave pictures, with no conclusive evidence for time reversal symmetry breaking. We argued for alternatives which can interpolate between the two simple cases without breaking the PG and TR symmetries via pairings with non-trivial matrix-structure in the orbital DOF. We discussed how matrix singlet pairings can emerge in unconventional SCs.

To support our general arguments, we considered the specific context of the Fe-based SCs. We present microscopic results showing that the phase difference of the intra-band ${d}_{{x}^{2}-{y}^{2}}$ and inter-band d_xy pairing components to be either 0 or π. This d + d pairing is the band basis equivalent of the sτ₃ form in the orbital basis, and is an irreducible B_1g representation of the (tetragonal D_4h) PG. We demonstrate that this d + d singlet pairing state is well defined, by showing that it can be compared and contrasted with the more familiar d + id state in a way analogous to how the well-defined B-phase in the case of superfluid ³He is measured against the equally well-known A-phase. The d + d pairing state allows for the reconcillation between seemingly contradictory experimental observations.

Non-trivial orbital structure can be relevant to unconventional SCs beyond the Fe-based family. To illustrate this, we constructed a pairing analogous to sτ₃ for the heavy-fermion CeCu₂Si₂ using general group-theoretical arguments. This sΓ₃ pairing state is also expected to have a d + d pairing structure in the band basis. It provides a natural theoretical basis to understand the striking low-temperature properties recently measured in the superconducting state of CeCu₂Si₂.

In these d + d pairing states, the anti-commuting nature of the two pairing components leads to their contributing to the single-particle excitation spectrum through an addition in quadrature, making it a fully-gapped superconducting state. The formation of the gap is connected to the energetic stabilization of such a state over a range of microscopic parameters. These results lead us to suggest that d − wave superconductors of strongly correlated multiorbital systems will inherently have a fully-gapped Fermi surface, even though the gap can be very small.

Note added

During the reviewing process of this manuscript, Ref. ⁷⁴ appeared with the results of recent x-ray spectroscopy experiments in CeCu₂Si₂ that support the sΓ₃ matrix pairing proposed here for CeCu₂Si₂. The sΓ₃ pairing includes paired Γ₇ f-electrons and Γ₆ conduction electrons. As the latter must hybridize with the excited Γ₆ f-electron states, a small but nonzero mixture of Γ₆ f-electrons is expected in the ground-state manifold. This mixture was shown for CeCu₂Si₂ in Ref. ⁷⁴.

Methods

Pairing channels of the five-orbital t − J ₁ − J ₂ model

We present our results for the five-orbital t − J₁ − J₂ model of the alkaline Fe-selenides¹⁷. The leading pairing amplitudes at zero-temperature are shown in Fig. 3 as functions of ${J}_{1}^{xz/yz}/{J}_{2}^{xz/yz}$. Exchange interactions in the xz/yz sector are identical for the two orbitals. ${J}_{2}^{xz/yz}=1/2$ in units of the band-width while the exchange interactions for d_xy orbital are 5 times smaller. Interactions in the remaining orbitals are ignored. J₁ and J₂ refer to their values for the xz/yz sector. For small and large values of the tuning parameter, ${s}_{{x}^{2}{y}^{2}}{\tau }_{0},{A}_{1g}$ and ${d}_{{x}^{2}-{y}^{2}}{\tau }_{0},{B}_{1g}$ orbital-trivial pairings are dominant. In the interval 0.8 ≤ J₁/J₂ ≤ 1, the ${s}_{{x}^{2}{y}^{2}}{\tau }_{3}$ pairing with non-trivial orbital structure is dominant with sub-leading ${d}_{{x}^{2}-{y}^{2}}{\tau }_{0}$ channel. The abrupt change around J₁/J₂ ≈ 0.75 is due to a transition from dominant A_1g to B_1g channels which become quasi-degenerate in this region. For the FS considered here, large J₂ favors dominant s_x2y2τ₀, A_1g while large J₁ favors d_x2−2y2τ₀, B_1g channels, respectively.

The form factors for the pairing terms have the standard expressions

$${s}_{{x}^{2}+{y}^{2}}({\bf{k}})=\frac{1}{2}\left[\cos ({k}_{x}a)+\cos ({k}_{y}a)\right]$$

(27)

$${s}_{{x}^{2}{y}^{2}}({\bf{k}})=\cos ({k}_{x}a)\cos ({k}_{y}a)$$

(28)

$${d}_{{x}^{2}-{y}^{2}}({\bf{k}})=\frac{1}{2}\left[\cos ({k}_{x}a)-\cos ({k}_{y}a)\right]$$

(29)

$${d}_{xy}({\bf{k}})=\sin ({k}_{x}a)\sin ({k}_{y}a).$$

(30)

Five-orbital t − J ₁ − J ₂ model and solution method

The pairing instabilities in the different symmetry channels of the alkaline Fe-selenides were obtained via a five-orbital t − J₁ − J₂ model:

$$\begin{array}{l}H=-{\sum \limits_{i\,{<}\,j}}({t}_{ij}^{\alpha \beta }{c}_{\alpha }^{\dagger }{c}_{\beta }+\,\text{h.c.}\,)+{\sum \limits _{i,\alpha }}\left({\epsilon }_{i\alpha }-\mu \right){n}_{i}+{\sum \limits _{\left\langle ij\right\rangle ,\alpha ,\beta }}{J}_{1}^{\alpha \beta }\left({{\bf{S}}}_{i\alpha }\cdot {{\bf{S}}}_{j\beta }-\frac{1}{4}{n}_{i\alpha }{n}_{j\beta }\right)\\ +{\sum \limits _{\left\langle \left\langle ij\right\rangle \right\rangle ,\alpha ,\beta }}{J}_{2}^{\alpha \beta }\left({{\bf{S}}}_{i\alpha }\cdot {{\bf{S}}}_{j\beta }-\frac{1}{4}{n}_{i\alpha }{n}_{j\beta }\right),\end{array}$$

(31)

where i, j indices cover all of the sites of a two-dimensional square lattice and α, β ∈ {1, …5} represent the ${d}_{xz},{d}_{yz},{d}_{{x}^{2}-{y}^{2}},{d}_{xy}$, and ${d}_{{z}^{2}}$ orbitals, respectively. The parameters of the model are specified in Ref. ⁶⁴. Different orbitals exhibit varying degrees of correlations, such that the exchange couplings are orbital dependent. More specifically, intra-orbital exchange couplings for the ${d}_{{x}^{2}-{y}^{2}}$ and ${d}_{{z}^{2}}$ orbitals are set to zero. Both NN and next-NN (NNN) exchange couplings are equal in the d_xz/yz sector and are larger by a factor of 5 than the exchange couplings in the d_xy sector. Inter-orbital exchanges have a small effect¹⁷ and are neglected here.

The interactions are decoupled in the particle-particle channel and the model is solved at T = 0 within a self-consistent approach. The double-occupancy constraint is introduced via an effective band renormalization^17,64.

We calculate the intra-orbital, NN and NNN pairing bonds, driven by J₁ and J₂ exchange couplings respectively, along $\hat{{\bf{x}}}$, $\hat{{\bf{y}}}$ and $\hat{{\bf{x}}}+\hat{{\bf{y}}}$ and $\hat{{\bf{x}}}-\hat{{\bf{y}}}$ respectively as

$${\Delta }_{{\bf{e}},\alpha \alpha }=\frac{1}{2}\left\langle {c}_{{{\bf{R}}}_{i}\alpha \uparrow }^{\dagger }{c}_{{{\bf{R}}}_{i}+{\bf{e}}\alpha \downarrow }^{\dagger }-{c}_{{{\bf{R}}}_{i}\alpha \downarrow }^{\dagger }{c}_{{{\bf{R}}}_{i}+{\bf{e}}\alpha \uparrow }^{\dagger }\right\rangle$$

(32)

where ${\bf{e}}\in \{\hat{{\bf{x}}},\hat{{\bf{y}}},\hat{{\bf{x}}}+\hat{{\bf{y}}},\hat{{\bf{x}}}-\hat{{\bf{y}}}\}$, R_i is a site vector, and α is an orbital index. The NN and NNN pairing bonds enter the pairing part of a Nambu Hamiltonian via

$${\Delta }_{{\bf{k}},\alpha \alpha }={\sum \limits_{{\bf{e}}}}{J}_{{\bf{e}}}\cos \left({\bf{k}}\cdot {\bf{e}}\right)$$

(33)

The dimensionless pairing amplitudes reported in the Results section are obtained by taking appropriate linear combinations of the NN and NNN pairing bonds. As such, the procedure does not bias towards any particular pairing channel. In addition, there are no approximations for the shape of the FS, and the pairing bonds are determined via averages where the momentum summation is over the entire Brillouin zone. The calculation is initiated with a random set of NN and NNN pairing bonds for all of the orbitals and subsequently allowed to converge. The procedure is repeated until a set of 300 converged solutions are obtained. From this set of converged solutions, we select the one which corresponds to the absolute minimum in the associated free-energy. This solution, again obtained without any superfluous conditions, corresponds to the physical solution reported in the manuscript.

Single-orbital t − J ₁ − J ₂ model

The Hamiltonian of the single d_xy orbital t − J₁ − J₂ model defined on a 2D square lattice is

$$H={H}_{\text{TB}}+{J}_{1}{\sum \limits _{\left\langle ij\right\rangle }}{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}+{J}_{2}{\sum \limits _{\left\langle \left\langle ij\right\rangle \right\rangle }}{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j},$$

(34)

where i, j label the lattice sites. The spin-density operators are defined as ${{\bf{S}}}_{i}=(1/2){\sum }_{ab}{c}_{a}^{\dagger }({{\bf{R}}}_{i}){{\boldsymbol{\sigma }}}_{ab}{c}_{b}({{\bf{R}}}_{j})$, where a, b are spin indices.

The TB part is determined by

$$\begin{array}{l}{H}_{TB}=-{t}_{1}{\sum \limits_{\left\langle ij\right\rangle }}{\sum \limits _{a}}{c}_{a}^{\dagger }({{\bf{R}}}_{i}){c}_{a}({{\bf{R}}}_{j})\\ -{t}_{3}{\sum \limits_{\left\langle \left\langle ij\right\rangle \right\rangle }}{\sum \limits_{a}}{c}_{a}^{\dagger }({{\bf{R}}}_{i}){c}_{a}({{\bf{R}}}_{j})-\mu {\sum \limits_{i}}{\sum \limits _{a}}{c}_{a}^{\dagger }({{\bf{R}}}_{i}){c}_{a}({{\bf{R}}}_{i}).\end{array}$$

(35)

The band is determined by

$$\begin{array}{l}\epsilon ({\bf{k}})=-2{t}_{1}\left[\cos ({k}_{x}a)+\cos ({k}_{y}a)\right]\\ -4{t}_{3}\cos ({k}_{x}a)\cos ({k}_{y}a)-\mu ,\end{array}$$

(36)

where a is the NN distance.

The TB coefficients are chosen as t₁ = 2t₃ = −0.5. The resulting band is shown in Fig. 6. Near half-filling we take μ ≈ −0.3 to obtain the FS shown in Fig. 7.

Fig. 6: Dispersion corresponding to Eq. (36) in units of 2∣t₁∣ with t₁ = 2t₃ = −0.5.

μ ≈ −0.3 ensures the FS shown in Fig. 7 near half-filling.

Fig. 7: Hole-like FS for a single-orbital d_xy model.

Here, the system is close to half-filling.

We implicitly take into account the renormalization of the bandwidth near half-filling by considering a large, fixed effective J₂ = −2t₁ = 1 while J₁ is allowed to vary. We decouple the exchange interactions in the pairing channels. The model is solved using the methods of Refs. ^17,64 near half-filling.