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τ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 + id state in a way analogous to how the B- triplet pairing phase of 3He superfluid differs from its A- phase counterpart. In addition, we construct an analog of the sτ3 pairing for the heavy-fermion superconductor CeCu2Si2, 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 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][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 24 .
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 2 Si 2 (Ref. 34 ). A host of properties, including the inelastic neutron-scattering spectrum 28 , have traditionally been interpreted in terms of a sign-changing d-wave pairing state, yet recent specific heat 35 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 37 or on the surface 38 , 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τ 3 17 . It has the s−waveform factor, but τ 3 , 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τ 3 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 39 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τ 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 structurein 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 2 Si 2 . This is an important undertaking, given that CeCu 2 Si 2 is the first-ever discovered unconventional SC 41 , 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Γ 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 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 multiorbital/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 electroncorrelation 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 3 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τ 3 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τ 3 in the band basis and illustrate how it is analogous to 3 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 3 He B and A. We show how they can be stabilized in a t − J 1 − J 2 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τ 3 for the heavy-fermion SC CeCu 2 Si 2 . In order to clearly see these results, we also present the irreducible representations of the D 4h point group in the context of CeCu 2 Si 2 . The Methods section contains additional accounts of the numerical results which support the stability of sτ 3 pairing for the alkaline Feselenides. 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 3 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τ 3 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 3 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ðk; k 0 Þ αβ;Γδ depend on the momenta as well as the orbital and spin indices of the two electrons. This twodimensional space turns out to be relevant for Fe-based SCs 53-55 , and we will first define the 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 analog of 3 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 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 analogs to the spin-1/2 product states in triplet 3 He (see Supplementary Note 1). 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ŝ g ðiÞ ðkÞ αβ ¼ g ðiÞ ðkÞτ ðiÞ αβ ; ( 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-wave states such as s x 2 þy 2 ðkÞ and s x 2 y 2 ðkÞ both belong to A 1g . Standard d-wave states 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 two-component E g representation. The τ ðiÞ αβ identity and Pauli matrices describe linear combinations of the tensorproduct 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. 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 τ 2 matrices, which would be parity-odd.
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. E.M. Nica and Q. Si 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 nontrivial 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 1 -J 2 model. 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 s-nor d-wave states. However, sτ 3 pairing preserves both PG-and TR-symmetries of the normal state.
To illustrate the properties of the sτ 3 pairing, we first consider a simplified two-orbital model 53 H Pair ðkÞ ¼ ΔðkÞs x 2 y 2 ðkÞτ 3 γ 1 : 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 pointgroup symmetry classification, sτ 3 transforms in the same B 1g representation as the diagonal-in-orbital-space d x 2 Ày 2 τ 0 pairing, as has been discussed in this type of model 39 and related settings [56][57][58][59] . What distinguishes the sτ 3 pairing is the nontrivial commutation relation between the corresponding pairing and kinetic parts of the Hamiltonian 17 .
It is instructive to recognize that the ξ 1 τ 1 and ξ 3 τ 3 terms ofĤ TB play a role similar to a Rashba SOC. The bands corresponding to the normal-state dispersion are 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 1 − J 2 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. 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 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τ 3 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τ 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 also be stabilized in more general five-orbital models of the alkaline Fe-selenide class of SCs. The pairing matrix in the t − J 1 − J 2 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 1 − J 2 model with nearest and next-nearest exchange couplings. This model and its solution method are discussed in subsection "Five-orbital t − J 1 − J 2 model and solution method" in Methods section.
The TB part and the associated FS are chosen to be consistent with LDA studies 17 . The dominant pairing amplitudes are intra-orbital. The sτ 3 pairing with non-trivial orbital structure is dominant in the 0.8 ≤ J 1 /J 2 ≤ 1.0 window. b Phases of the leading B 1g channels relative to the sτ 3 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 selfconsistent solution. In the [0.8, 1] interval where sτ 3 is dominant, these relative phases are either zero or ± π. Here, the amplitudes of the coexisting B 1g 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 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 orbitalselective correlation effects in the normal state. 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 subsection "Pairing channels of the five-orbital t − J 1 − J 2 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 τ 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 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 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 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 α 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 inter-band dwaves 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, as shown in Fig. 2. 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. 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.
The band basis reveals a pairing structure which is very similar to that in 3 He -B. Referring to Supplementary Notes 1 and 2, 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-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τ 3 in the band basis, where two d-wave states 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 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, 3 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 3 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 3 He -B vs. 3 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 nontrivial sτ 3 pairing is equivalent to simultaneous intra-and interband 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 α 3 and α 1 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  In the first subsection, 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 the Methods section. 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 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 either zero or ± π. 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).
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 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 1 − J 2 model" of the Methods. The model is solved using a selfconsistent decomposition of the exchange interactions as in the five-orbital cases discussed previously 17,64 . The resulting zerotemperature 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. 4a. For J 1 /J 2 < 0.8, the only significant pairing occurs in the d xy , B 2g channel. For higher values of J 1 /J 2 , the amplitude of a d x 2 Ày 2 ; B 1g pairing becomes finite. These two remain 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 65 .
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 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.   17 for a detailed account of the model and solution. For J 1 /J 2 ≤ 0.7 A 1g pairing channels are dominant with leading s x 2 y 2 τ 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 1 /J 2 ≤ 0.8. Beyond this range, B 1g channels dominate with leading s x 2 y 2 τ 3 in the xz/yz sector which has non-trivial orbital structure. At even larger values, an orbital-trivial d x 2 Ày 2 τ 0 phase dominates. When the tuning parameter is less than 0.8, only the d xy , B 2g channel has finite amplitude. In the 0.8 ≤ J 1 /J 2 ≤ 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 1 /J 2 . 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.

E.M. Nica and Q. Si
By contrast, 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 analog 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. 5 of Supplementary Note 2. The band and spin matrices in the 3 He -A and d + id cases are analogous.
In the first subsection, 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 (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 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 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 3 He -B, we believe that revealing the above connections elevates the status of the d + d spinsinglet 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 2 Si 2 .
Matrix singlet pairing with spin-orbit coupling: CeCu 2 Si 2 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 2 Si 2 , which is the first-ever discovered unconventional SC 41 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 spinresonance peak in the SC state 28 together with angle-resolved resistivity measurements of the upper critical field H c2 66 , among others. Remarkably, recent measurements of the specific heat 35 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τ 3 provide a natural resolution. An analogous proposal for CeCu 2 Si 2 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 2 Si 2 24,34 .
Objective and outline of the subsection. Here, we construct the analog of the sτ 3 pairing for CeCu 2 Si 2 . For reasons that will become clear, we shall refer to this state as sΓ 3 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 2 Si 2 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 69 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 lowest-lying excited states of the f electron are composed of a Γ 6 and another Γ 7 doublet. Our analysis will also involve Γ 6 Wannier orbitals of the conduction electron states near the Fermi energy, and these Wannier orbitals will hybridize with the excited Γ 6 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Γ 3 , 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 Γ 6 Wannier orbitals from the Cu 3d orbitals.
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 analog of A 1g , while Γ þ 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 Γ 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 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 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 70 .

E.M. Nica and Q. Si
Conventional B 1g 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 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 70 Here, Γ þ 1;2 are one-dimensional representations which are analogous to the A 1g and B 2g in the absence of SOC 70 . The twodimensional Γ þ 5 is analogous to the xz, yz (E g ) doublet. Following Ref. 70 , 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 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 Γ À 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.
Matrix B 1g 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. Further below, we illustrate how intra-unit cell linear combinations of Cu 3d states in the presence of SOC can form bases for Γ À 6 conduction electrons. The product states decompose as 70 The corresponding matrices are 70 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 analog 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 s-wave form factor which transforms according to the Γ þ 1 trivial irrep. sΓ 3 pairing, which involves electrons belonging to different irreducible representations due to the Γ 3 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 Ξ 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 34 , 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 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 Γ 7 -Γ 6 , sΓ 3 pairing represents an sτ 3 analog 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 − wave states. Such a pairing can in principle reconcile the difficulties in interpreting the more recent experimental results.
We note that the Γ À 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 Γ À 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Γ 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 the first subsection, we expect that the sΓ 3 matrix-pairing will coexist with a conventional d-wave pairing 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 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 73 Table 1.
We follow the conventions of Ref. 70 . E represents the identity, C n are rotations about z by 2π/n, while C 0 2 and C 00 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.
Γ À 6 conduction electrons. Previously, we discussed matrix B 1g pairing between Γ À 6 c -electrons and Γ À 7 f-electrons. In this subsection, we illustrate how linear combinations of intra-unit-cell Cu 3d orbitals in the presence of SOC provide Γ À 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 2 Si 2

69
. The Cu sites are located halfway along the edges of a plaquette as illustrated in Fig. 5.
The linear combinations of the four orbitals at the Cu sites x 2 Ày 2 À d ð2Þ x 2 Ày 2 (23) x 2 Ày 2 À d ð3Þ transform as an (x, y) doublet under the D 4h point group.
Consequently, (p x , p y ) belong to a two-component Γ À 5 irreducible representation. We further include the local spin-1/2 DOF which belongs to a Γ þ 6 irreducible representation 70 . From the direct-product states of p-orbital linear combinations and spinor states ϕ ±1/2 we can construct states which belong to Γ À 6 representation 70 : where SOC was taken into account. Similar states can be constructed in the remaining Cu plane. The symmetric linear combination between Γ À 6 states in both Cu planes likewise belongs to Γ À 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τ 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 lowtemperature 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 singleparticle 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. 74 appeared with the results of recent x-ray spectroscopy experiments in CeCu 2 Si 2 that support the sΓ 3 matrix pairing proposed here for CeCu 2 Si 2 . The sΓ 3 pairing includes paired Γ 7 f-electrons and Γ 6 conduction electrons. As the latter must hybridize with the excited Γ 6 f-electron states, a small but nonzero mixture of Γ 6 f-electrons is expected in the ground-state manifold. This mixture was shown for CeCu 2 Si 2 in Ref. 74 .

METHODS
Pairing channels of the five-orbital t − J 1 − J 2 model 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  Interactions in the remaining orbitals are ignored. J 1 and J 2 refer to their values for the xz/yz sector. For small and large values of the tuning parameter, s x 2 y 2 τ 0 ; A 1g and d x 2 Ày 2 τ 0 ; B 1g orbital-trivial pairings are dominant. In the interval 0.8 ≤ J 1 /J 2 ≤ 1, the s x 2 y 2 τ 3 pairing with non-trivial orbital structure is dominant with sub-leading d x 2 Ày 2 τ 0 channel. The abrupt change around J 1 /J 2 ≈ 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 2 favors dominant s x2y2 τ 0 , A 1g while large J 1 favors d x2−2y2 τ 0 , B 1g channels, respectively.
Five-orbital t − J 1 − J 2 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 1 − J 2 model: 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. 64 . 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 17 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 doubleoccupancy constraint is introduced via an effective band renormalization 17,64 .
We calculate the intra-orbital, NN and NNN pairing bonds, driven by J 1 and J 2 exchange couplings respectively, alongx,ŷ andx þŷ andx Àŷ respectively as Δ e;αα ¼ 1 2 c y Ri α" c y Ri þeα# À c y Ri α# c y Ri þeα" where e 2 fx;ŷ;x þŷ;x Àŷg, 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 Δ k;αα ¼ X e J e cos k Á e ð Þ 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 1 − J 2 model The Hamiltonian of the single d xy orbital t − J 1 − J 2 model defined on a 2D square lattice is where i, j label the lattice sites. The spin-density operators are defined as E.M. Nica and Q. Si The band is determined by ϵðkÞ ¼ À2t 1 cosðk x aÞ þ cosðk y aÞ Â Ã À4t 3 cosðk x aÞ cosðk y aÞ À μ; 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. 6. Near half-filling we take μ ≈ −0.3 to obtain the FS shown in Fig. 7.
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 model is solved using the methods of Refs. 17,64 near halffilling.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.