Multiorbital singlet pairing and $d+d$ superconductivity

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\tau_{3}$ 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 + i d$ state in a way analogous to how the B- triplet pairing phase of $^{3}$He superfluid differs from its A- phase counterpart. In addition, we construct an analogue of the $s\tau_{3}$ pairing for the heavy-fermion superconductor CeCu$_{2}$Si$_{2}$, 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 be fully gapped when they are examined at sufficiently low energies.


I. 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 wellstudied Cu-based superconductors [1,2]. However, the pairing symmetry remains enigmatic in other classes of strongly correlated materials. For singlet superconductivity, the longheld dichotomy is between fully-gapped sand 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 (SC) [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 [14]. 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] and experimentally [20,21]. 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 [22].
Recent experiments have directly challenged the conventional sand d-wave dichotomy. In alkaline Fe-selenides, inelastic neutron scattering [23,24] revealed signatures of ingap spin resonances, whose characteristic wavevectors qualify them as a typical indicator of sign-changing d-wave order parameters [25]. By contrast, ARPES studies have indicated fully-gapped superconductivity [26][27][28][29], even for a Fermi pocket near the center of the two-dimensional Brillouin zone (BZ), which appears consistent with swave symmetry.
A similar situation has emerged in the heavy-fermion sys- * Corresponding author: enica@asu.edu tem CeCu 2 Si 2 (Ref. 30). A host of properties, including the inelastic neutron-scattering spectrum [31], have traditionally been interpreted in terms of a sign-changing d-wave pairing state, yet recent specific heat [32] and London penetration depth [30,33] results at very low temperatures pointed toward a fully gapped Fermi-surface (FS). It is surprising that the SC phases exhibit sand d-wave characters simultaneously. One possible origin is s + id pairing, which breaks the point-group (PG) and time-reversal (TR) symmetries. However, in the absence of conclusive evidence pointing to either TR symmetry-breaking or two-stage phase transitions as the temperature is lowered, an alternate candidate pairing must be identified. For the Fe-chalcogenide SC, one candidate pairing state was named sτ 3 [17]. It has the s−wave form factor, but τ 3 , a Pauli matrix in the xz, yz 3delectron orbital basis, turns the pairing state into d−wave-like. 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. Yet, the sτ 3 pairing has the intra-and inter-band d + d form, which is highly unusual, thereby raising the question of both its naturalness and generality.
Here, we argue that the d+d pairing state belongs to matrix singlet pairing order parameters with non-trivial orbital structure 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τ 3 pairing state. Written in the band basis, the sτ 3 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 3 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.
We also 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 2 Si 2 . This is an important understaking, given that CeCu 2 Si 2 is the firstever discovered unconventional SC [34], and also recognizing that heavy fermion systems represent a prototype setting for strong correlations and unconventional superconductivity in general [35][36][37]. Using DOF that incorporate spin-orbit coupling (SOC), we introduce an sΓ 3 state. This provides the theoretical basis for the excellent description of the experimental results in CeCu 2 Si 2 in terms of the d + d pairing order parameter [22,30].
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.
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 3 He . However, candidate solid state systems for triplet pairing often involve multiple orbitals/bands and strong correlations [38][39][40][41][42][43]. Thus, the type of matrix pairing structure in the orbital/band space we consider here may produce new types of triplet superconducting states [44,45] and the associated excitations that are of potential interest to quantum computing.
The remainder of the paper is organized as follows. We begin Sec. II 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τ 3 pairing, and discuss the unusual properties of this state and show how it can be stabilized in an sto d-wave degeneracy regime in Sec. II B. We also support our discussion with numerical results for realistic models of the Fe-based SCs. In Sec. II C, we consider sτ 3 in the band basis and illustrate how it is analogous to 3 He -B. In Sec. III, 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 analogues of 3 He B and A. We show how they can be stabilized in a t − J 1 − J 2 model. In Sec. IV, we extend the notion of nontrivial orbital structure beyond the Fe-based compounds by discussing a candidate analogous to sτ 3 for the heavy-fermion SC CeCu 2 Si 2 . We summarize our results and discuss their implications in Sec. V. In Appendix A, for completeness, we outline the role of the matrix-pairing functions in the various phases of superfluid 3 He , where spin here provides the analogue of the orbital DOF. We highlight the lessons we believe can be applied to the case of multiorbital pairing in unconventional SCs. In Appendix B, we illustrate how sand d-waves can coexist without breaking either PG or TR-symmetries via a general Landau-Ginzburg analysis. Appendix C contains an additional account of the numerical results which support the stability of sτ 3 pairing for the alkaline Fe-selenides. In Appendix D, we discuss the t − J model and its solution which illustrates the case of d + id pairing. Appendix E contains a brief review of the D 4h PG and its irreducible representations.
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 lowest-lying bands is due to xz and yz orbital local DOF. For simplicity, we ignore SOC. In general, the pairing interactions V (k, k ) αβ;γδ 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 [46][47][48][49], and we will first define sτ 3 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 analogue of 3 He -B.

A. 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 The even-parity matrixĝ(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 analogues to the spin-1/2 product states in triplet 3 He (see Appendix A). Likewise, they depend on the relative momentum of the pair. Finally, iσ 2 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 aŝ 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 (i) (k) can likewise be chosen to belong to one of these representations. To illustrate, s-waves such as s x 2 +y 2 (k) and s x 2 y 2 (k) both belong to A 1g . Standard d-waves such as d x 2 −y 2 (k) and d xy (k) are B 1g and B 2g representations, respectively. The xz, yz orbital doublet transforms as the twocomponent E g representation. The τ (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 3 He , these matrix-elements play the role of effective Clebsch-Gordan coefficients. The τ 0 , τ 1 , and τ 3 matrices transform according to A 1g , B 2g , and B 1g , respectively.
These arguments point to an important aspect. In 3 He , the relative angular momentum and local (spin) DOF transform independently under separate groups. In the present case, g (i) (k) and orbital matrix parts (τ (i) ) 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 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 fiveorbital t-J 1 -J 2 model.

B. sτ3 pairing state
Of interest here is the sτ 3 pairing. In terms of Eqs. 1 and 2, it corresponds to It transforms as the B 1g representation due exclusively to the τ 3 matrix. Because of the orbital struture, it represents neither simple snor d-wave. However, sτ 3 pairing preserves both PG-and TR-symmetries of the normal state.
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 the ξ 1 τ 1 and ξ 3 τ 3 terms ofĤ T B play a role similar to a Rashba SOC. The bands corresponding to the normalstate dispersion are reflecting the space-group allowed varying orbital-content and splitting of the FSs. We refer the reader to Refs. 17 and 46 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. 17, we showed that the general BdG dispersion is always gapped along the FS. Nodes away from the FS can appear for larger band splitting, reflecting the corrections to the gap term due to the non-commuting TB and pairing parts. However, in alkaline-Fe selenides, where sτ 3 pairing is competitive, the small band splitting at the center of the Brillouin zone precludes the appearance of nodes.
Another important characteristic of such a gapped sτ 3 state is its sign change under a π/2 rotation. Such a sign-change leads to the formation of an in-gap spin resonance. sτ 3 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τ 3 pairing, the latter can be stabilized in realistic five-orbital models of the alkaline Fe-selenide class of SCs. The zero-temperature pairing amplitudes of all symmetry-allowed pairing channels have been determined in a five-orbital t − J 1 − J 2 model with nearest and next-nearest exchange couplings. The tight-binding (TB) part and the associated FS are chosen to be consistent with LDA studies [17] and orbital-selective pairing was taken into account by tuning the strength of the pairing interactions in the xz, yz versus xy orbitals. We focus on the case where pairing mainly involved the xz, yz orbitals while also allowing inter-orbital pairing. In the following, J 1 and J 2 refer to their values for the xz/yz sector. The J 1 /J 2 ratio controls the symmetry of the dominant pairing channel with s x 2 y 2 (k)τ 0 and d x 2 −y 2 (k)τ 0 states for small and large values of the ratio, respectively. The sτ 3 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 1 /J 2 ≈ 1. A d x 2 −y 2 τ 0 with trivial orbital structure provides the subleading pairing with comparable amplitude. See Appendix C for more details.
A dominant sτ 3 pairing emerges in the vicinity of the J 2 /J 1 ≈ 1 point. In our model, the J 2 exchange interaction equally projects [50] pairings with s x 2 y 2 and d xy form factors. The J 1 exchange equally prefers [50] s x 2 +y 2 -and d x 2 −y 2 -wave order parameters. The inter-orbital pairings [51] involving τ 1 are suppressed, even when inter-orbital exchange interactions are included [17,50].
In addition, the FS with electron pockets at the BZ edge suppresses a conventional s x 2 +y 2 pairing. Finally, d xy pairings which have vanishing amplitude on the the same pockets are likewise not favored. Thus, in the J 1 ≈ J 2 regime, interactions strongly favor s x 2 y 2 and d x 2 −y 2 with either τ 0 or τ 3 .
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 τ 0 pairing of comparable amplitude is present in the regime where sτ 3 is dominant. While we have focused on the properties of the dominant sτ 3 pairing alone, a more realistic picture would involve coexisting sτ 3 and d x 2 −y 2 τ 0 in the vicinity of sto 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 Appendix B, 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 .

C. Intra-and inter-band d + d pairing and its analogy with the 3 He B-phase
It is instructive to consider the sτ 3 pairing in a band basis: where α 1,3 are Pauli matrices in the two-band space. Because the overall pairing is in the irreducible B 1g channel, it is natural that the intraband α 3 part has the d x 2 −y 2 form factor. Likewise, the interband α 1 component has the d xy form factor. Thus, the sτ 3 pairing is equivalent to a d + d pairing.
When the pairing matrix is squared, the intra-and interband d-waves add in quadrature as ∆ 2 1 (k) + ∆ 2 2 (k) to produce a gap which does not close along the FSs centered on the BZ center corresponding to the two bands. This is due to the anti-commuting nature of the two Pauli matrices α 3 and α 1 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. In cases relevant to our discussion, these additional effects are typically small and consequently do not close the gap.
The band basis reveals a pairing structure which is very similar to that in 3 He -B. Referring to Appendix A, the matrix order-parameter in that case is typically expressed aŝ ∆3 He−B (k) ∼ (k·σ)iσ 2 . This amounts to a linear superposition of p-waves 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τ 3 in the band basis, where two d-waves likewise produce a finite gap. The sτ 3 pairing thus provides a remarkable example where a phase which is similar to 3 He -B via a structure in the band-basis is stabilized in a solid-state SC model.
Along with this similarity between the sτ 3 pairing state and the B phase of 3 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 3 He as contrasted with Pairing Amplitude Phase Wrt s In the [0.8, 1] interval where sτ3 is dominant, these relative phases are zero. Here, the amplitudes of the coexisting B1g channels are comparable to that of sτ3 . This illustrates that the sτ3 pairing which is equivalent to d + d, effectively preserves TR and PG symmetries.
the discrete nature of the PG in the inter-and intra-band dwave 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, 3 He -B breaks the SO(3) L × SO(3) S symmetry of the normal state down to SO(3) L+S . 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 3 He phases. In this respect, it still resembles intra-and inter-band d-wave pairing which preserves both PG and TR symmetries.
We have seen that 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 experiment. We have also seen that the non-trivial sτ 3 pairing is equivalent to simulta-neous 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 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 α 3 and α 1 Pauli matrices, respectively. Importantly, we assume that 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 3 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 3 He .
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 1 − J 2 model in the vicinity of the J 1 ≈ J 2 point.
In Sec. II B, we discussed how the sτ 3 orbital non-trivial pairing channel becomes dominant for a finite range of the J 1 /J 2 tuning parameter in a realistic five-orbital t − J 1 − J 2 model for the alkaline Fe-selenides. The details of the calculations are given in Appendix C. We showed how sτ 3 pairing is equivalent to a d + d intra-and inter-band pairing. To further understand the nature of the sτ 3 -dominated state, we plot the phases of the leading B 1g channels relative to sτ 3 as functions of J 1 /J 2 modulo π in Fig. 1 (b). The leading B 1g channels have relative phases wrt sτ 3 which are closely centered around 0 or π for the entire range of the tuning parameter. In the [0.8, 1] interval where sτ 3 is dominant the relative phases are zero. Here the amplitudes of the subleading d x 2 −y 2 τ 0 and s x 2 +y 2 τ 3 B 1g channels are comparable to that of the leading sτ 3 . 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τ 3 is dominant which we consider here (Fig. 3 in Appendix C).
We next discuss the orbital-trivial d + id pairing. For simplicity, we consider a single d xy orbital t − J 1 − J 2 model on a square lattice. We choose the 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 Appendix D. The model is solved using a mean-field decomposition of the exchange interactions as in the five-orbital cases discussed previously [17,52]. The resulting zero-temperature pairing amplitudes for J 2 = 1/2 in units of the bandwidth, and for a finite range of the ratio J 1 /J 2 are shown in Fig. 2 (a). Pairing Amplitude Results for a single-orbital t − J1 − J2 model close to halffilling. Please see Appendix D for details of the model. (a) Pairing amplitudes as functions of the ratio J1/J2. When the tuning parameter is less than 0.8, only the dxy, B2g channel has finite amplitude.
For larger values of the tuning parameter, the dxy channel is suppressed and s x 2 +y 2 and s x 2 y 2 A1g channels emerge. (b) Relative phase of the two d-wave channels modulo π as a function of J1/J2. 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.
finite up to J 2 /J 1 ≈ 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 [53].
To illustrate that the two coexisting d-wave components are locked into a d+id state, we plot their relative phases mod π in Fig. 2 (b). While these relative phases are essentially arbitrary whenever one of the d-waves vanishes, a π/2 relative phase is clear in the interval J 1 /J 2 ∈ [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. In Sec. II C, we showed that the d+d pairing is closely analogous to the B phase of 3 He , where the pairing is a superposition of p x,y,z -waves corresponding to equal-and oppositespin pairing as illustrated by Eqs. A2 and A3 of Appendix A. 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 constrast, the intra-band (d + id)α 0 pairing, where α 0 is the identity matrix in the band basis, is a linear superposition of two irreducible representations.
In general, the α 0 matrix in band space would correspond to an identity τ 0 matrix in orbital space. The d + id pairing spontaneously breaks both PG and TR symmetries. Therefore, it is a natural analogue of 3 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. A5 of Appendix A. The band and spin matrices in the 3 He -A and d + id cases are analogous.
In Sec. II C, we discussed how d + d differs from 3 He -B due mainly to the discrete versus continuous symmetries of the two, respectively. This kind of difference also exists between d + id and 3 He -A. The latter spontaneously breaks both angular momentum and spin continuous rotational symmetries down to a U Lz−φ × U Sz group of independent rotations in each sector about preferred axes (Appendix A). 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 3 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 3 He -B and d + d and likewise between 3 He -A and d + id, which are due to the presence of non-trivial matrix structure, are intriguing. In this sense, 3 He provides both a wellestablished parallel and a prototype for the emergence and description of the effects of non-trivial matrix structure in unconventional singlet SCs.

IV. MATRIX SINGLET PAIRING WITH SPIN-ORBIT COUPLING: CECU2SI2
Another class of multiband superconductors arises in heavy fermion systems, in which quasi-localized f electrons hybridize with dispersive spdconduction (c) electrons. These include CeCu 2 Si 2 , which is the first-ever discovered unconventional SC [54] 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 [31] together with angle-resolved resistivity measurements of the upper critical field H c2 [55], among others. Remarkably, recent measurements of the specific heat [32] and London penetration depth [30,33] 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τ 3 provide a natural resolution. An analogous proposal for CeCu 2 Si 2 is clearly of great interest. Note that a d + d interand intra-band pairing directly inspired from the Fe-based cases provided a good fit to the the London penetration depth measurements in CeCu 2 Si 2 [22,30].
In contrast to the case of the Fe-based SC, SOC is important in CeCu 2 Si 2 . The local orbital and spin DOF transform simultaneously under PG operations. This imposes additional constraints on any matrix associated with the local DOF. A number of experiments [56,57] as well as LDA+DMFT studies [58] have indicated that one of the Γ 7 doublets of the crystal-field split 2 F 5/2 local electron is the dominant contribution to the heavy FS sheets. The low-lying excited states of the f electron are composed of Γ 6 and Γ 7 doublets. The

A. Spin-orbital coupled local states
Our aim is to propose a minimal symmetry-allowed candidate for CeCu 2 Si 2 which has properties similar to that of sτ 3 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 Γ + 1 is the analogue of A 1g , while Γ + 3 is the analogue 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 Γ 6 or Γ 7 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 TRinvariant. 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 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 [59].

B. Conventional B1g pairing
We first consider candidates on the Γ − 7 ground-state doublet. The superscript denotes odd parity under inversion. Although these naturally correspond to f -f pairing involving the Γ 7 ground-state local multiplet, they also cover possible f -c pairing 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 [59] Here, Γ + 1,2 are one-dimensional representations which are analogous to the A 1g and B 2g in the absence of SOC [59]. The two-dimensional Γ + 5 is analogous to the xz, yz(E g ) doublet. Following Ref. 59, the matrices corresponding to each of the three Γ + 1,2,5 representations are: 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 Σ Γ1 which transforms as the trivial representation. This matrix is the analogue of simple singletpairing 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 Γ − 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 Γ 3 representation which is analogous to a B 1g representation without SOC. We do so by choosing g(k) = d x 2 −y 2 (k) andΣ = Σ Γ1 in Eq. 11. Likewise,Ξ can be chosen to be proportional to either Ξ 1 or (1/2)(Ξ 0 − Ξ 3 ), where Ξ 0 and Ξ 1,3 denote identity and Pauli matrices, respectively. The two cases correspond to f -c and f -f pairing, respectively.

C. Matrix B1g pairing
We next consider pairing between electrons belonging to distinct Γ − 7 and Γ − 6 manifolds. This can correspond to f -f pairing between electrons belonging to the Γ − 7 ground-state and f electrons belonging to the excited Γ − 6 manifolds, respectively. Alternately, it can denote f -c pairing between Γ − 7 f -electrons and Γ − 6 conduction c electrons. As before, c electrons would most likely have p-orbital character. The product states decompose as [59] The corresponding matrices are [59] 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 Γ 7 − Γ 7 case discussed previously. The only anti-symmetric matrix is Σ Γ3 . It transforms as the Γ + 3 irrep of D 4h and is an analogue of the τ 3 matrix in the Fe-based cases. Moreover, this matrix is invariant under TR. We conclude that a counterpart of the sτ 3 pairing for CeCu 2 Si 2 is an s-wave pairing belonging to the sign-changing Γ 3 representation, or sΓ 3 pairing: where s(k) corresponds to a sign-changing s-wave form factor which transforms according to the Γ + 1 trivial irrep. As in the sτ 3 case, the sign-changing nature of the form factor is important since the wavefunction changes signs and thus must vanish when r Relative = 0. 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 Ξ 1 Pauli matrix. In the f − f case, it can be made proportional to a Ξ 0 −Ξ 3 matrix. In either case, the gap is determined by the amplitude and form factor only similarly to what happens for sτ 3 . In a multiband model of CeCu 2 Si 2 [30], this pairing produces a full gap.
On general grounds, the non-trivial Γ 7 -Γ 6 pairing in either f -f or f -c cases is likely weaker than the Γ 7 -Γ 7 f -f pairing . However, there are cases where such Γ 7 -Γ 6 contributions can be important. Consider a dominant Γ 7 -Γ 7 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 ff 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 orbitalspin 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 Γ 7 -Γ 6 , sΓ 3 pairing represents an sτ 3 analogue for CeCu 2 Si 2 . As in the Fe-based cases, the structure of the local DOF allows a natural interpolation between simple s− and d−waves. Such a pairing can in principle reconcile the difficulties in interpreting the more recent experimental results.
Similar to the Fe-chalcogenide case, for unitary pairing we expect the sΓ 3 pairing in the band basis to contain the intraband α 3 and interband α 1 components. Each must be in d-wave state, with the form factor of the intraband α 3 being d x 2 −y 2 . Thus, the sΓ 3 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 Sec. II A, we expect that the sΓ 3 matrixpairing will coexist with a conventional d-wave pairing (Sec. IV B) below T c since they both belong to the same Γ 3 irreducible representation of D 4h . The admixture between these will ensure that the SC state preserves both PG and TRsymmetries but also exhibits a gap which is finite everywhere along the FS.
We stress that our analysis is distinguished from the wellknown symmetry-based procedure typically considered in the context of heavy-fermion SCs [60]. The latter do 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.

V. CONCLUSIONS
Recent experiments in multiband Fe-based and heavyfermion SCs are inconsistent with either simple sor 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 matrixstructure 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τ 3 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 3 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τ 3 for the heavy-fermion CeCu 2 Si 2 using general group-theoretical arguments. This sΓ 3 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 2 Si 2 .
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 gapped spectrum, even though the gap can be very small. The p-wave spin-triplet pairing in 3 He leads to a great variety of superfluid phases. These states are captured by a nontrivial matrix structure of the pairing in spin space. Here, we review some of the most important superfluid phases of 3 He paying particular attention to the form of the pairing matrix. In the main text (Sections II C, III A, and III B), we draw analogies between superfluid 3 He and spin-singlet multiband SCs where the matrix structure is associated with local orbital DOF. These analogies can be made in spite of such obvious differences as the real-space symmetry which is continuous in 3 He and discrete in the SC cases and the normal-state bands which are typically more complex in SCs than in 3 He case.
Cooper pairing is the basic feature of all SCs. Electrons are bound in pairs to form a state which saves ground-state energy. Even in the simplest BCS description of Cooper pairing, the spin-1/2 DOF of the two electrons play a crucial role. The bound state has a net zero total spin. Paired electrons accordingly behave like bosons which effectively undergo condensation into the superconducting ground state. In turn, this makes it possible to understand the emergence of superconductivity as macroscopic breaking of a U (1) symmetry associated with the net condensate. This macroscopic breaking of symmetry underlies most of the remarkable properties observed in experiments [61,62]. The success of the basic picture of super-conductivity is due in large part to the fundamental concepts of spontaneous symmetry-breaking with little reliance on the detailed microscopic knowledge of the interactions responsible for the Cooper instability.
In addition to the net spin, Cooper pairs can also have net relative angular momentum (AM).
We restrict the discussion to homogeneous SC states. The symmetry of the Cooper pairs in that case is captured within a BdG formalism by the pairing function is a pairing amplitude. The matrix elements of g Relative AM denote the weights of the tensor-product states of the local DOF. These weights are in general functions of the relative AM of the Cooper pair. In the simplest spin-1/2 singlet-pairing case, these matrix elements are the Clebsch-Gordan coefficients of the total spin S = 0 of the Cooper pair.
In the absence of detailed microscopic models, symmetry provides a powerful tool in understanding the nature of the transition to the paired state as well as a host of the latter's experimental signatures. The case of superfluid 3 He is illustrative. The basic microscopic constituents are essentially featureless spin-1/2 atoms. Normal 3 He is a Fermi liquid which preserves rotational invariance in real space, spin space, and a global phase [62]. Mathematically, this defines a group G = SO(3) L ⊗ SO(3) S ⊗ U (1). Interactions which lead to Cooper pairing preserve this symmetry. In the weak-coupling regime, the pairing amplitude is dwarfed by an energy scale which characterizes the normal liquid [63]. To a good approximation, the interactions and the pairing function can be projected onto sectors belonging to one of the irreducible representations of G [62,63]. Each representation is assigned a critical temperature T c determined by the eigenvalues of the linearized gap equation. The irreducible representation with the highest T c determines the nature of the paired state. Instead of going into the microscopic analyses [64], we will focus on symmetry-based analyses for the superfluid 3 He . The p-wave spin-triplet channel with relative AM L = 1 and total spin S = 1 was singled-out as the most likely candidate [63].
The only essential local DOF in 3 He are the spin-1/2 moments of the pairing atoms. In the physically-relevant triplet sector the order parameter can be written as [60] where the σ Pauli matrices act on the spin-1/2 tensor-product states. Odd-parity, spin-exchange symmetric pairing is restricted to the three complex components of d, which transforms as a vector under spin rotations [62]. The physical significance of this vector is important. For unitary states, where id × d * = 0, it can be shown [62,65] that the vector d is real within ak-dependent overall phase. In these cases, d points along the direction of total S = 0.
Due to the separate SO(3) symmetry of both relative AM and spin sectors, the order parameter has a large degeneracy. The superfluid phases of 3 He can be understood in great part via the breaking of the symmetry of G in order to relieve the inherent degeneracy. The local spin DOF play a crucial role in the nature of the resulting phases, as we now proceed to review.

Contrasting B and A phases
In zero magnetic field and for most pressures, 3 He condenses directly into the superfluid B phase. At low temperatures it is the most stable state. The B phase came to be identified with a unitary, "Balian-Werthamer" (BW) [66] order-parameter. In it's simplest form, the BW vector orderparameter is The orientation of the spins and of the relative AM are locked together, reflecting the breaking of symmetry down to a group G = SO(3) L+S of simultaneous rotations. The simplest BW state then corresponds to a unique J = 0 irreducible representation of G . The gap, given bŷ is finite everywhere along the Fermi surface (FS) for both spin-1/2 quasiparticles. The emergence of the uniform gap can formally be traced to the anti-commuting Pauli matrices. Therefore, it is clear that the local spin structure of this pairing plays a very important role. Due to the uniform gap, the BW order-parameter is the most favored configuration at weakcoupling [66].
In a restricted region of the pressure-temperature phase diagram, 3 He condenses directly into the A phase described by an unitary "Anderson-Brinkman-Morel" (ABM) orderparameter vector [62,63] d ABM =d k ·m + ik ·n . (A5) The system spontaneously selects preferred quantization axes for the spin and relative AM of the pair, respectively. Equivalently, the symmetry is spontaneously broken to a group G = U Lz−φ × U Sz , where φ denotes a global phase rotation. This is a unitary order parameter and the vectord thus determines a fixed direction in spin-space along which the pairs have a vanishing dipole-moment S Pair = 0. For example, when the vectord lies in the xy plane, the system exhibits equal-spin pairing in each of the independent sectors with total S z = ±1 while pairing does not occur for opposite-spin partners. The relative AM likewise acquires a preferred direction which is determined by the triad of mutually-perpendicular unit vectorsm,n andl. The latter determines the direction in k-space where the gap vanishes and thus where nodes occur. This phase also breaks time-reversal (TR) symmetries [67]. Due to the appearance of nodes, the ABM order-parameter is not stable at weak-coupling where the BW is always preferred [62,63]. The ABM pairing is stabilized when feedback effects are taken into account [62]. The P − T phase diagram reflects this, as the A phase only occurs in a restricted region and always gives way to the B phase at lower temperatures [62].
The phases of 3 He in magnetic fields H are also instructive. The order parameters in these cases turn out to be non-unitary. These phases also manifest orientation effects due to preferred axes for the spins of the paired electrons. To illustrate, we focus on the A phase. Near T c and under an applied field, the A phase gives way to a non-unitary A 1 phase. Pairing occurs in only one of the sectors with total S z = ±1, as a consequence of the lifting of the degeneracy of the spin-up and -down FSs. For lower temperatures, the non-unitary A 2 phase stabilizes. It is described by the complex vector Pairing now occurs in both of the sectors of total S z = ±1. However, due to the split of the FS under applied field, the two coefficients A(B) ∼ ∆ ↑↑ ± ∆ ↓↓ reflect different pairing amplitudes for the two sectors. Both A 1,2 phases are non-unitary pairing states where ∆ † (k)∆(k) is not simply proportional to the identity matrix. In such cases, it can be shown that a pair at k acquires a finite spin dipole-moment [62,65]: In the simplest cases, this spin-dipole moment points along the applied field H. The vectorsd andê are perpendicular to the applied field, while the local moment is parallel to a mutually-perpendicular directionf =d ×ê. A comparison of A and B phases of superfluid 3 He is important for our analysis of SC phases with non-trivial orbital local DOF. The B phase is the more symmetric of the two in the sense that it's residual-symmetry is still enhanced with respect to that of the A phase. The spin and relative AM DOF are locked to produce a maximal uniform gap. This phase is consequently always preferred at lower temperatures. The more constrained A phase, which further breaks the rotation and TR symmetries has nodes and is only stabilized at higher pressure and temperature. We believe that these salient properties are not restricted to superfluid 3 He but can also manifest in solid-state SC with non-trivial local DOF. Thus, we expect that pairing can naturally take advantage of additional structure due to local orbital DOF to induce the strongest and most uniform gap along the FS, as in 3 He -B.
Appendix B: Landau-Ginzburg theory for coexisting inequivalent representations Here we construct a generalized LG theory following the arguments of Sec. II A. We consider pairing withing a twoorbital xz, yz model. We restrict our discussion to intraorbital pairing interactions where the Greek indices denote the two orbitals. Such pairing interactions can originate from a t−J 1 −J 2 model [17,50,52]. For these types of interactions, inter-orbital pairing characterized by the off-diagonal τ 1 matrix are suppressed. Furthermore, we consider that s x 2 y 2 and d x 2 −y 2 channels are equallyfavored. This would correspond to J 1 ≈ J 2 in a t − J 1 − J 2 model. As discussed in Sec. II A, we expect that orbital structure plays an important role near this point. Hence, we consider two, nearly-degenerate s x 2 y 2 τ 3 and d x 2 −y 2 τ 0 order parameters. Each of these transforms as the single-component B 1g irreducible representation of D 4h . We denote the order-parameters associated with s x 2 y 2 τ 3 and d x 2 −y 2 τ 0 by ∆ 1 and ∆ 2 , respectively. The general LG free-energy for these two candidates is Terms involving α 1,2 and β 1−4 are also allowed in the case of almost-degenerate but orbital-trivial s and d order parameters [68]. Couplings involving α 3 , β 5 , and β 6 are made possible by the orbital structure. We assume that α 3 is independent of temperature.
In any SC phase, we require thatα < 0,β > 0. We assume that near the superconducting transition, quartic terms can essentially be neglected. The minimum of F LG is then determined by the minimum ofα. For given α 1−3 , we find the extrema ofα as a function of θ, φ for coexisting channels with θ = 0, π/2. We obtain: The second condition implies φ = 0 or π. Using double-angle formulas, the first condition is For simplicity we consider α 3 > 0. To see what the extrema determined above imply, we consider two limiting cases. First, we consider almost degenerate channels with |(α 1 −α 2 )/2α 3 | 1. The two extrema correspond to leading order to for φ = π and 0, respectively. It is straightforward to show that the solution corresponding to φ = π is always smaller than either α 1 or α 2 which correspond to a pairing in either s x 2 y 2 τ 3 or d x 2 −y 2 τ 0 channels below T c . This implies that the bilinear terms in the LG theory always prefer a coexisting solution immediately below T c . The other limiting cases occur for |(α 1 − α 2 )/2α 3 | 1. In this case, θ is close to either 0 or π/2. Below T c one channel will will dominate with a very small admixture of the other.
In general, we expect that no additional transitions occur below T c although this can occur in some instances as illustrated by the analogous case of 3 He -A 1 and 3 He -A 2 .
These arguments illustrate that coexisting sτ 3 and dτ 0 channels are in general favored when pairing interactions are similar in the two cases. We present our results for the five-orbital t − J 1 − J 2 model of the alkaline Fe-selenides [17]. The leading pairing amplitudes at zero-temperature are shown in Fig. 3 where i, j label the lattice sites. The spin-density operators are defined as S i = (1/2) ab c † a (R i )σ ab c b (R j ), where a, b are spin indices.
The TB part is determined by The band is determined by where a is the NN distance. The TB coefficients are chosen as t 1 = 2t 3 = −0.5. The resulting band is shown in Fig. 4. Near half-filling we take µ ≈ −1.3 to obtain the FS shown in Fig. 5. We implicitly take into account the renormalization of the bandwidth near half-filling by considering a large, fixed effective J 2 = −2t 1 = 1 while J 1 is allowed to vary. We decouple the exchange interactions in the pairing channels. The Pairing Amplitude

Appendix E: Irreducible representations of D 4h
We include the character table for the double-valued representations of the tetragonal D 4h point group.
We follow the conventions of Ref. 59. E represents the identity, C n are rotations about z by 2π/n, while C 2 and C 2 are π rotations about x, y and (−x, y), (x, y) axes, respectively. I denotes inversion. S 4 indicates a C 4 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 Γ + 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 Γ − 1−5 correspond to the A 1u , A 2u , B 1u , B 2u , and E u representations. Γ combinations of Γ + 1−4 and similarly for the odd-parity Γ − 6 and Γ − 7 .