Pseudospin-Triplet Pairing in Iron-Chalcogenide Superconductors

We study superconductivity of electron systems with both spin and pseudospin-1/2 degrees of freedom. By solving linearized gap equations, we derive a weak coupling criterion for the even-parity spin-singlet pseudospin-triplet pairing. It can generally mix with the on-site s-wave pairing since both of them belong to the same symmetry representation ($A_{1g}$) and their mixture could naturally give rise to anisotropic intra-band pairing gap functions with or without nodes. This may directly explain why some of the iron-chalcogenide superconductors are fully gapped (e.g. FeSe thin film) and some have nodes (e.g. LaFePO and LiFeP). We also find that the anisotropy of gap functions can be enhanced when the principal rotation symmetry is spontaneously broken in the normal state such as nematicity, and the energetic stabilization of pseudospin-triplet pairings indicates the coexistence of nematicity and superconductivity. This could be potentially applied to bulk FeSe, where gap anisotropy has been experimentally observed.


Introduction
The symmetry principle is one of the most powerful tools to diagnose low-energy electronic band structures, lattice vibrations, and linear responses [1], and is also valuable to explore various symmetry-breaking ordered phases such as magnetism, charge/spin density-wave, nematicity and superconductivity [2].The crystal symmetry of a solid-state system dictates the normal band structures it hosts near the Fermi level, which could in turn determine the most favorable superconducting pairing symmetry [3,4].This symmetry principle for superconductors (SC) is recently extended to investigate multi-band unconventional superconductivity [5][6][7].Interestingly, the orbital-independent and orbital-dependent pairings that belong to the same symmetry representation may coexist with each other [8].Such orbital-dependent pairings have been studied in a wide variety of systems with multiband character, including Sr 2 RuO 4 [9], iron-chalcogenide SCs [10][11][12][13], Cu-doped Bi 2 Se 3 [14] and half-Heusler compounds [15][16][17][18], from which the guiding principle by symmetry has been shown to be crucial to understanding the nature of unconventional superconductivity.
A few specific systems can be effectively characterized by a general normal-state model Hamiltonian that contains both spin ({↑, ↓}) and pseudospin ({1, 2}) degrees of freedom, where pseudospin could originate from two atomic orbitals, two sublattices, two layers, or two valleys [6].We start from a spin-singlet centrosymmetric SC to explore the existence of even-parity pseudospin-triplet * hu.lunhui.zju@gmail.compairings, for example, , and further investigate their valuable roles in tailoring anisotropic pairing gap functions with or without nodes [19].Different from spintriplet pairings, spin-singlet pseudospin-triplet pairings have not been much explored in real materials since such pairings are usually considered to be energetically unfavorable.This is partly due to the common belief that the double degeneracy of the two orbitals is lifted by orbital hybridization so that the orbital-dependent pairing would be severely suppressed under crystal field splitting or electron-electron repulsive interaction.One aim of this work is concerned with the possible condition for the existence of even-parity spin-singlet orbital-dependent pairings, and possible applications to real materials.
On the other hand, the effects of symmetry breaking in unconventional SCs is an important topic that has attracted tremendous interest.The symmetry could be broken explicitly by external fields or strain, or be broken spontaneously from many-body interactions.Two typical examples are rotational symmetry breaking [20,21] and time-reversal-symmetry (TRS) breaking [22][23][24][25][26]. Besides, the interplay between nematicity and superconductivity is yet to be fully understood in some real materials, such as FeSe [27,28], where gap functions can be highly anisotropic.These systems are all multi-band SCs, while symmetry-reducing signatures are experimentally observed above the superconducting transition temperature, which is mainly caused by both crystal field splittings and interaction-induced order parameters (e.g.nematicity).Thus, discovering the coexistence of nematicity and superconductivity in these multi-band systems arXiv:2109.06039v4[cond-mat.supr-con]5 Dec 2023 can shed new light on understanding the underlying favorable pairing symmetries.
The main finding of this work is that the anisotropic gap functions with or without nodes could be attributed to the mixing of isotropic s-wave pairing and even-parity spin-singlet pseudospin-triplet pairing, even though both of them belong to the A 1g symmetry representation.For technical conveniences, we adopt an orbital d o (k)-vector notations [11] to describe the pairing matrix and similarly a g o (k)-vector for orbital hybridization in the twoorbital subspace ({1, 2}).Solving linearized gap equations, we show that the presence of g o -vector generally suppresses the superconductivity with orbital d o -vector except for d o (k) ∥ g o (k), which is consistent with the concept of superconducting fitness [6].This sets up weakcoupling criteria for A 1g -type orbital-dependent pairings that could naturally give rise to anisotropic gap functions in real superconducting materials.Moreover, we reveal a deep connection between two-orbital nematic SC and pseudospin-triplet pairings.Within the mean-field theory for electron-electron repulsive interactions, the nematic order develops in the orbital subspace at T < T nem , which also contributes to the total orbital hybridization, g tot = g o + g nem .This leads to the stabilization of a nematic orbital d o -vector for d o (k) ∥ g tot (k), indicating the coexistence of nematicity and superconductivity.The direct applications to FeSe [27,28] are also discussed.We also generalize it to a two-valley system with C 6 breaking terms (e.g., Kekulé distortion).In the end, we also predict an orbital-polarized superconducting state.

Results
Classification of Spin-singlet Orbital-triplet pairings.To explore the weak-coupling criterion for the energetically favorable even-parity spin-singlet pseudospintriplet pairing, we consider the mean-field pairing Hamiltonian, where F † s1a,s2b (k) = c † s1a (k)c † s2b (−k) is the creation operator of Cooper pairs, s 1 , s 2 are indices for spins and a, b are for pseudospins (e.g., two orbitals {1, 2}).A general pairing potential of a two-band model is a four-by-four matrix [6].In particular, the spin-singlet pairing function ∆ a,b s1,s2 (k) = f (k)M a,b (k)(iσ 2 ) s1,s2 consists of the angular form factor f (k) and M a,b (k) in the orbital channel.The spin-singlet pairings are not mixed with spin-triplet pairings in the absence of spin-orbit coupling (SOC).In analogy to spin-triplet SCs, for the technical convenience, we then use an orbital d o (k)-vector for the spin-singlet orbital-dependent pairing potential [11], (2) where ∆ s and ∆ o are pairing strengths in orbitalindependent and orbital-dependent channels, respectively.Here τ and σ are Pauli matrices acting on the orbital and spin subspace, respectively, and τ 0 is a 2-by-2 identity matrix.When both ∆ s and ∆ o are real, a real orbital d o (k)-vector preserves TRS while a complex one spontaneously breaks TRS (T = iτ 0 σ 2 K with K being complex conjugate).The Fermi statistics requires Ψ s In other words, d 2 o (k) describes oddparity spin-singlet orbital-singlet pairings and the other two are for even-parity spin-singlet orbital-triplet pairings.Moreover, we provide an alternative definition of orbital d o -vectors in Appendix A. Even though the orbitalindependent part Ψ s (k) is also "orbital-triplet" by statistics, it is completely trivial.Hereafter, we only refer to d 1 o (k) and d 3 o (k) as orbital-triplet pairings [29].In addition, the basis functions for both Ψ s (k) and orbital d o (k)-vectors in Eq. ( 2) could be classified by crystalline symmetry.
Under the action of an n-fold rotation operator C n about the z-axis, the pairing potential ∆(k) transforms as where D[C n ] is the corresponding matrix representation, J is the orbital angular momentum quantum number, and also labels the irreducible representations of the C n point group.For example, J = 0 is for A representation and J = 2 is for B representation.Firstly, the TRS requires the coexistence of ∆J and ∆−J with equal weight.If the rotation symmetry C n is further imposed, then J and −J have to be equivalent modulo n, i.e.J ≡ −J mod n.
At the mean-field level, the Bogoliubov de-Gennes (BdG) Hamiltonian is given by where H 0 (k) represents a two-band normal-state Hamiltonian with both spin and pseudospin degrees of freedom.
In general, the BdG Hamiltonian is also invariant under the C n rotation symmetry, i.e., Here we assume both inversion and time-reversal symmetries are preserved.To be specific, we consider a SOCfree Hamiltonian, where the basis is Classification of spin-singlet pairing potentials for Eq. ( 2).Here we consider a spin-singlet two-orbital superconductors with {dxz, dyz}-orbitals.Based on the n-fold rotation symmetry Cn about z-axisd and TRS, we have J = −J mod n, which leads to all the pairing channels with orbital-independent Ψs(k) and orbital-dependent do(k)-vector in Eq. ( 2).Here, for J = 0 pairing subspace of C3, the (d relative to the chemical potential µ, m is the effective mass, λ o represents the orbital hybridization and . And the g 3 -component leads to the different effective masses of different orbitals.
As mentioned earlier, this vector notation is just for the technical convenience.Besides, the g 1 and g 2 components are determined by symmetries.For example, TRS requires , which is the same as the constraint for the orbital d o -vector.The more explicit form of g o (k) is determined by other crystal symmetries.
In general, the pseudospin-triplet (i.e.orbital-triplet) pairing state shares some similarities with the spin-triplet pairing state [31].To show that, we first discuss the superconducting quasi-particle spectrum of orbital-triplet SCs in the absence of band-splitting caused by orbital hybridizations, i.e., g o (k) = 0 for Eq.(5).In this case, the superconducting gaps on the Fermi surface are for the ∆ s = 0 limit.This indicates that there are two distinct gaps if TRS is spontaneously broken.In the following, we mainly focus on the time-reversal-invariant superconducting states, i.e., real d o -vectors, for which the classification of pairing potentials is shown in Table (I) based on Eq. ( 3).We will show the interplay between ∆ s and ∆ o can lead to anisotropic superconducting gaps on different Fermi surfaces.Moreover, its stability against orbital-hybridization, electron-electron interactions, and applications to real materials will be discussed in detail as follows.We will also briefly comment on the effects of TRS-breaking in the end.
We apply the weak-coupling scheme [6] for spin-singlet orbital-triplet pairings against crystal field splittings, which cause orbital hybridizations [i.e. the g o (k) term in Eq. ( 5)].We analytically calculate the superconduc-tivity instability for the orbital d o -vector by BCS decoupling scheme.The superconducting transition temperature T c of orbital-dependent pairing channels is calculated by solving the linearized gap equation, where the Matsubara Green's function for electrons with ω n = (2n + 1)π/β and G h (k, iω n ) = −G * e (k, iω n ).We expand the attractive interactions as V s1a,s2b Here Γ labels the irreducible representation with l = 1, 2, ..., Dim Γ.In this work we focus on 1d representations, i.e.Dim Γ = 1, which already include many interesting cases and are sufficient for the applications discussed in later sections.Due to the possible existence of multiple pairing channels belonging to different representations, each channel has its own critical temperature T Γ c , the largest of which becomes the actual critical temperature of the system.In the weak-coupling theory, T Γ c follows the standard BCS form and is solely determined by the corresponding pairing interaction v Γ in that particular channel.To the leading order of ), the equation for T c for the channel Γ reads (see details in the Methods section), ln where T c0 is the critical temperature for λ o = 0, Ω is the solid angle of k, ĝo = g o (k)/|g o (k)| are normalized vectors.Here we take , where ψ (0) (z) is the digamma function.
We now discuss its implications.In general, the λ oterm describes a pair-breaking term, since C 0 (T c ) ≤ 0 and it monotonically decreases as λ o increases, hence the right-hand side of Eq. ( 8) suppresses T c in general.However, if we focus on one-dimensional representations, i.e.Dim Γ = 1, it is straightforward to see that d Γ o ∥ g o can lead to T c = T c0 for any value of λ o , which indicates that the orbital d o -vector that is parallel with g o is unaffected by the orbital hybridizations.It is worth mentioning that due to the possible suppression of T c , depending on the relation between d Γ o and g o the leading instability channel at λ o = 0 could be suppressed more than some of the other coexisting channels and may eventually become sub-leading.This interesting behavior is discussed further in Appendix C. For notional simplicity, we will drop the representation index Γ when there is no danger of confusion.Choosing g o (k) = (2k x k y , 0, k 2 x − k 2 y ), the numerical results are shown in Fig. 1.The black line confirms that T c is unaffected as , which is the unconventional A 1g pairing.However, T c for other d ovectors are severely suppressed.The light-blue line is for d o (k) = 1 √ 2 (1, 0, 1), and the light-orange line for . Therefore, we conclude that the orbital d o -vector could exist in SCs with two active orbitals that are not fully degenerate.This is similar to spin-triplet SCs, where the A 1g -type spin d svector could exist in noncentrosymmetric SCs because d s ∥ g s is optimally satisfied [4,6].
It is worth mentioning that the results presented above is using a continuum form of the Hamiltonian based on k • p theory.For real materials, given the interaction on the lattice, the components of the interaction in terms of the basis functions of the representations might not be exactly the same with the form of the vector g o .As a result, the parallel condition presented above may not be exactly satisfied.However, the theory developed in this work is generally applicable and the extend to which the parallel condition holds can still be a useful criterion for the most favorable pairing.
Next, we include ∆ s , and investigate the coupling between Ψ s and d o .Solving the coupled linearized gap equations up to (λ o k 2 F /µ) 2 order (see details in Appendix C, we find that the results from Eq. ( 8) are still correct.Besides, the magnitude of orbital d o -vectors might be determined as d o (k) = Ψ s (k)ĝ o (k).It implies that Ψ s and d o belong to the same representation of crystalline groups.Therefore, the stability of orbital d o -vector by Eq. ( 8) indicates the symmetry principle for spin-singlet orbital-triplet pairings.
We now explain Eq. ( 8) from the band picture.
Within the band basis, the pairing potential in the orbital subspace becomes ∆band with the normal band dispersion The intra-orbital pairing naturally gives rise to the (1,0,1) Stability of orbital do-vectors vs orbital hybridization λ0 in Eq. (5).Shown are the transition temperature Tc/Tc0 as a function of λok 2 F /kBTc0 for go(k) = (2kxky, 0, k 2 x − k 2 y ).Tc0 is Tc at λ0 = 0.The curves from top to bottom correspond to do(k intra-band pairing.However, it is different for orbitaldependent pairings.To show that, we decompose the orbital We find that the d ∥ -part gives rise to the intra-band pairing, while the d ⊥ -part leads to the inter-band pairing (see Appendix D).If the band splitting is much larger than the pairing gap (λ o k 2 F ≫ ∆ o ), the inter-band pairing is not energetically favorable in the weak-coupling pairing limit.It means that the inter-band pairing will be severely suppressed if we increase the orbital hybridization λ o , consistent with Eq. ( 8) and results in Fig. 1.Now if we again include the orbital-independent pairing part ∆ s Ψ s (kτ 0 iσ 2 ), the relation between d o and Ψ s (k) obtained previously from solving the coupled linearized gap equation (see Appendix C) can also be reproduced in the band picture by considering the maximization of the condensation energy.The total condensation energy per volume and per spin of the two intra-band pairings is given by where N ± are the density of states on the two Fermi surfaces (E ± ).And ∆ s Ψ s (k) ± ∆ o d ∥ (k) are the pairing gaps on these two Fermi surfaces.In order to maximize δE, we have Applications to superconductors with/without nodes.As a consequence of the mixing of the orbitalindependent pairing (∆ s ) and orbital-dependent pairing (∆ o ) discussed in the previous section, there could be a nodal SC.In this section, we apply the results of the previous section to study superconductors with two orbitals, where ∆ s and ∆ o coexist.It is shown that the anisotropic gap functions with/without nodes depend on the ratio of ∆ s and ∆ o superconducting order parameters.Our weak-coupling theory might have potential applications to some of the nodal/nodeless SCs in the ironchalcogenides family.For example, the angle-resolved photoemission spectroscopy (ARPES) measurements indicate a nontrivial superconducting gap anisotropy for the monolayer FeSe thin film [32].The penetration depth measurements on both LaFePO [33] and LiFeP [34] show a linear dependence on T , suggesting the presence of superconducting gap nodes.
As an example, we consider the pairing potential in Eq. ( 2) for monolayer FeSe, where there is no hole pocket around the Γ-point, and a two-spin two-orbital model has been shown to be a good approximation around the electron pockets near the M point of the Brillouin zone (two Fe unit cell).The density functional theory calculations show that there are four bands around the M point, giving rise to only two electron pockets.In the one Fe unit cell, there is one pocket near the X and Y points, respectively.After folding with respect to the unit cell with two Fe, we obtain two pockets around the M point.Considering spin degrees of freedom, it naturally resembles a C 4z -invariant two-orbital model [35], where ϵ(k) = (k 2 x + k 2 y )/(2m) − µ with m > 0 the effective mass, A leads to the anisotropic effective mass (i.e., orbital hybridization), and v so represents SOC that still preserves inversion symmetry.These four states are degenerate at the M point since they form the fourdimensional representation of the space group No. 129 (P 4/nmm) [36].We take the parameters for the FeSe thin film as µ = 55 meV, 1/(2m) = 1375 meV• Å2 , A = 600 meV• Å2 and v so ≤ 15 meV• Å [35].The SOC is very weak to open a tiny gap along the k x = 0 and k y = 0 lines, shown in Fig. 2 (a).As what we expect, it shows two C 4z rotational-invariant Fermi surfaces, and the maximal gap, which is induced by the z-component of the g o vector, is around 12 meV along the (11) and (1 1) directions.This is larger than the typical superconducting gaps in iron-chalcogenide SCs (∼ 4 meV), implying that the effect of the orbital hybridization on the pairing symmetries should not be neglected.
We now use the criteria derived above (Eq.( 8)) to examine the superconducting states.Specifically, the weakcoupling criterion indicates that the most favorable pair-ing to characterize the anisotropic superconducting gap is the A 1g -type s-wave pairing symmetry, The ratio between ∆ s and ∆ o determines the superconducting nodal structure.To simplify the analysis, we turn off the weak SOC.In the band basis, the dispersion of Here ± label the band index.Projecting ∆(k) onto the bands leads to ∆ ± = ∆ s ± ∆ o |k x k y |.Given that ∆ s , ∆ o > 0, nodal points can only appear for ∆ − on the "−" band.The nodal condition would be |k x k y | = ∆ s /∆ o has solution on the FS given by ϵ − (k) = 0.By using the mathematical inequality k 2 x + k 2 y ≤ 2|k x k y |, it can be shown that the nodal condition is given by, which is shown in Fig. 2 (b).In general, the ratio ∆ s /∆ o should depend on both interaction strength in each pairing channel and the orbital hybridization strength.This gives rise to the condition of nodal A 1g -type s-wave superconducting states.Therefore, it could not only explain the anisotropic gap functions observed in the FeSe thin film (fully gapped) but also the nodal superconductivity in LaFePO and LiFeP.Around one linear Dirac node, the effective Hamiltonian up to linear-k can be mapped out as where k 1 , k 2 are liner combinations of k x and k y .All the other Dirac nodes are related to this one by reflection symmetries.Then, we only need to focus on H D , which is a Dirac Hamiltonian with topological charge (winding number) ±2, whose node is protected by the chiral symmetry (i.e., the product of time-reversal symmetry and particle-hole symmetry).The 2Z winding number is due to the presence of inversion symmetry and time-reversal symmetry.To analytically show the topology of Dirac nodes, we apply perturbation analysis with respect to PT symmetry (i.e., the product of time-reversal symmetry and inversion symmetry) and Chiral symmetry.Note that the PT symmetry can be also C 2z T symmetry for a 2D or quasi-2D SC.The projected symmetry representations are given by PT = σy τ0 and C = σy τx .As expected, the PT symmetry commutes with H D , while the Chiral symmetry anti-commutes with H D .Then, local perturbations preserving PT and Chiral are where m 1 and m 2 represent perturbation strengths or mass terms.The spectrum of H D + H ′ D are given by Nodal SC Fully gapped which indicates that the Dirac nodes are movable but not removable.For example, k 1 = 510.7kx + 76.5k y and k 2 = −14.7kx − 40.9k y around one Dirac node.Then, turning on the SOC v so = 15 meV• Å, we numerically confirm the nodal SC phase with ∆ s = 3 meV and ∆ o = 200 meV• Å2 , shown in Fig. 2 (c), where the logarithm of superconducting gaps are plotted.The eight dark red points are the linear Dirac nodes.Based on the topology-protection argument, the interplay between intra-and inter-orbital pairings for nodal superconductivity is robust against local perturbations.Note that our results are different from a previous work [35], in which the d-wave pairing symmetry induced nodal SC.In experiments, the nodal gap structure could be detected by measuring the temperature dependence of physical quantities like specific heat and penetration depth at low temperatures.A power law dependence usually indicates the existence of nodal structures (point nodes or line nodes), whereas exponential dependence implies the SC is fully gapped [3].
Applications to superconductors with nematic order.In addition to the crystal field splitting, the many-body electron-electron interactions may also lead to orbital hybridization, such as the nematic ordering in the normal states (See Appendix E for details).The rotational symmetry reduction could either be from interaction-induced spontaneous symmetry breaking or from explicit symmetry breaking from, say, adding external strain.Then the natural question to ask is whether it is still possible to have an orbital-dependent pairing order characterized by some d o -vector.Interestingly, we find that the orbital-dependent pairing can coexist with the electronic nematic ordering as long as d o is parallel to the g tot , which is an effective orbital-hybridization vector that also contains the nematic order.This establishes a deep connection between SCs with nematic order and spin-singlet orbital-triplet pairings.In the following, we study two typical examples.
• For case A [two-orbital system], we apply the theory to fit the anisotropic superconducting gap of the hole pocket in the bulk FeSe measured by the quasiparticle interference imaging [28].
• For case B [two-valley system], we use a toy model to demonstrate the possible existence of s + d-like nematic nodal superconductor in two-valley systems on a honeycomb lattice.We also show the transition between U-shaped and V-shaped quasiparticle density-of-state by tuning the chemical potential.
Case A: Two-orbital model for the bulk FeSe SC.
We discussed the possible anisotropic A 1g -type s-wave pairing states for the C 4 -symmetric iron-chalcogenide SCs including fully gapped FeSe thin film and nodal SC in LiFeP and LaFePO.Here we investigate the C 4 -breaking nematic SC in bulk FeSe.Let us revisit the iron-based SC with a well-established nematic ordering.We consider , where ni is electron density operator for the i-atomic orbital.If ) is spontaneously broken down to C 2 and we have the nematic order.The intra-orbital interaction does not alter the mean-field results for nematic orders (see details in Appendix E).The total inter-orbital hybridization contains two parts, where g o (k) is caused by the crystal field splitting and g nem = (0, 0, Φ) is induced due to the nematicity Φ = v 1 ⟨n 1 − n2 ⟩, which is momentum-independent if translation symmetry is to be preserved.Hereafter, we focus on the hole pockets around the Γ point to fit the experimental data of superconducting gap functions [28].We  18).All the other parameters used are the same [37], including the chemical potential, effective mass, orbital hybridization, and nematic order.
will see that even this simplified weak-coupling model, where the coupling between the hole pockets at the Γ point and the electron pockets at the M point is ignored, can produce a descent fit the experimental data.A similar result is expected for the electron pockets near the M point.Replacing g o with g tot in Eq. ( 5), we can still use Eq. ( 8) to investigate the interplay between superconductivity and nematic order, thus the orbital d o -vector satisfying d o ∥ g tot leads to the nematic superconductivity.Thus, it generally shows the A 1g -type s-wave spin-singlet orbital-triplet pairings in nematic SCs.This scenario can be adopted to study the quasi-two dimensional bulk FeSe, where superconductivity (T c ∼ 8 K) emerges inside a well-developed nematic phase (transition temperature T nem ∼ 90 K [38]), shown in Fig. 3 (a).For a minimal two-band model [39] for the bulk FeSe with {d xz , d yz }-orbitals, g o = (2k x k y , 0, k 2 x − k 2 y ) and g nem = (0, 0, Φ) [37,40].Therefore, the nematic orbital d o -vector with d o ∥ g tot breaks C 4 (see Appendix E for more details).The projected pairing gap function on the large Fermi surface is given by that is in the isotropic limit.The presence of Φ is the driving force for the anisotropy of ∆ FS (k).When the nematicity Φ is strong enough, the orbital d o -vector will be pinned along the z-axis, resulting in the so-called orbitalselective pairing states.We adopt the realistic parameters for the bulk FeSe SC from Ref. [37] to calculate the superconducting gap measured by the quasiparticle interference imaging [28].In Fig. 3 (b), we show the angular dependence of the pairing gap around the hole pocket at the Γ-point of FeSe in the presence of nematic order.
Our theory provides an equally decent fit to recent experimental data [28] as the intra-orbital s + d-pairing theory proposed by Kang et al. [37], even though our work uses a simplified model without considering the coupling to the other two electron pockets.Our theory shows more clearly the role of nematic order on the pairing symmetries.Therefore, the theory developed in this work may alternatively explain the experimental evidence of orbital-selective pairings of the FeSe SC in Refs.[27,28], and reveal a deep connection between nematic SC and spin-singlet orbital-triplet pairings.It has to be mentioned that here we only focused on the hole pockets around the Γ point and discussed the nematicity-induced gap anisotropy around the hole FS.There are other possible mechanisms for gap anisotropy in Fe-based SCs.For example, a previous work [41] discussed, among other things, the anisotropy/isotropy of the SC gap around the electron pockets at the M point, where the degree of anisotropy depends on the J 1 -J 2 magnetic frustration in the proposed five-orbital t-J 1 -J 2 microscopic model.
Case B: Two-valley system superconductivity.Similar to the two-orbital systems considered above, we discuss in this section superconductivity in two-valley systems, like single layer graphene SC [42] or transition metal dichalcogenide (TMD) [43], where the pairing can be between opposite valleys K ± .The spin-singlet pairing is merely characterized by the orbital d o -vector with ∆ s = 0 in Eq. ( 2).For the single-particle Hamiltonian, the inter-valley hopping is naturally forbidden by translational symmetry, namely, λ o = 0 in Eq. ( 5).Then, we consider the inter-valley scattering Hamiltonian, k).It generates the inter-valley coupling g int by defining the order parameter Φ ⟩ that spontaneously breaks the translational symmetry, where Both d 1 (k) and d 2 (k) are real to preserve TRS.As for the interactioninduced g int , T and I require g int,1 (k) = g int,1 (−k) and g int,2 (k) = 0.By symmetry, there are two general possibilities.One is g int,1 (k) = 1, so C 3 × I = C 6 is preserved, and it describes the charge-density-wave order [44,45].The other one is that spontaneously breaks C 6 down to C 2 , forming a nematic order.This is experimentally possible for the strain-induced Kekul'e distortion (i.e., √ 3 × √ 3 type).We next discuss superconductivity in the presence of inter-valley couplings, by replacing the g o -vector in Eq. ( 5) with the interaction-induced g int .As a result, Eq. ( 8) is still applicable.It is similar to a recent work [46] where the charge order coexists with a sublattice-selective non-unitary pairing state.
The nematic inter-valley coupling is represented as Here the normalization factor has been dropped without changing the essential physics.The system is fully gapped if the s-wave gap is dominant ).As a concrete toy model, we look at superconductivity on a generic honeycomb lattice with two valleys K ± , with the Hamiltonian around the two valleys given by, where the two-valley basis used here is given by ψ ) and ϵ(k) takes the same form as in Eq. ( 5).The parameter α determines the C 3 anisotropy of the continuum model around each valley.This Hamiltonian was used as an effective model [47] to study twisted bilayer graphene.
Including the inter-valley scattering effects, the oneband model is given by where the λ int determines the strength of the intervalley scattering.In the absence of inter-valley scattering (λ int = 0), the Fermi surfaces (FSs) around the two K ± valleys are plotted in Fig. 4 (a).As expected, with a fully symmetric s-wave pairing characterized by d o = (1, 0, 0), a fully gapped or U-shaped quasi-particle density-ofstates (DOS) is obtained and shown in Fig. 4 (d).Then we include the aforementioned inter-valley scattering g int that breaks C 6 down to C 2 .As a result, our theory implies that the effective nematicity generated will favor a nematic pairing characterized by d o ∥ g int .Consider the generic form x − k 2 y ), 0, 0), the resulting C 6 -breaking FS are shown in Fig. 4 (b) and (c), where the nodal lines of the pairing are also shown.By tuning the chemical potential µ, the FSs and nodal lines can go from non-intersecting in Fig. 4 (b) to inter-secting in Fig. 4 (c), leading to the corresponding evolution from the gapped U-shaped DOS in Fig. 4 (e) to the gapless V-shaped DOS in Fig. 4 (f).

Discussions
We briefly discuss the difference between our theory and the previous studies [21] for nematic SCs.One example is a pairing state belonging to a 2D irreducible representation (Irrep), e.g., the E-pairing in Cu or Nb-doped Bi 2 Se 3 [48,49] and UPt 3 [50,51].A real order parameter vector (∆ E,1 , ∆ E,2 ) spontaneously breaks C 3 , leading to nematic superconductivity.Alternatively, a nematic SC can be formed by mixing two different 1D-Irrep-pairing channels.In FeSe [52,53], the nematic order breaks the C 4 down to C 2 , which mixes the s-wave and d-wave pairing channels.However, T c of the (s + d) orbitalindependent pairing state could be generally affected by increasing nematicity, because of the significant change in the density of states at the Fermi energy.In our theory, the (s + d)-like nematic d o -vector coexists with the nematic order, so T c is almost unaffected by increasing nematicity.Therefore, it may help to distinguish our results from previous proposals in experiments, where one may use the chemical or physical pressures to tune the nematicity and measure T c as a function of pressure [54].Nevertheless, more efforts are necessary to test the results established in this work for nematic SCs.
In addition to the above discussions for the timereversal-invariant superconducting states, we also comment on the effects of the spontaneous TRS-breaking, where a complex orbital d o -vector generates the orbital orderings as M o = −iγ 1 /α M (d × d * ), of which only the y-component breaks TRS, as illustrated in Fig. 5 (a).Alternatively, the corresponding quasi-particle spectrum in Fig. 5 (b) shows the two distinct gaps, similar to the range given by Eq. ( 6).More explicitly, we schematically plot the atomic orbital-polarized density of states (DOS) by defining |±⟩ = |1⟩ + i|2⟩ for complex orbitals, where D + ̸ = D − at finite energy clearly indicates that the DOS is orbital-polarized, which is consistent with the Ginzburg-Landau theory, shown in Appendix B.Moreover, we also find that the orbital-spin conversion would lead to the spin-polarized DOS [55].
The above result for orbital-triplet pairings is similar to the superconducting gaps for non-unitary spin-triplet SCs [3].By symmetry, the Gizburg-Landau free energy is the same.To show the similarity, for the singleband spin-triplet SCs [56], the spin-triplet pairing potential is generally given by ∆ , where ∆ 0 is the pairing strength and σ are Pauli matrices in the spin subspace.Due to the Fermi statistics, the spin d s (k)-vector has to satisfy d s (k) = −d s (−k).The d s -vector formalism is firstly developed for He 3 superfluid [57].And it also occurs in noncentrosymmetric SCs, the spin d s (k)-vector is usually pinned along a certain crystal axis since superconductivity is nonsuppressed only for d s (k) ∥ g s (k), where g s (k) represents the Rashba spin-orbit coupling (SOC) [4,6].Besides, It is an equal-spin pairing so that σ 3 is conserved, and nonzero M s leads to two distinct superconducting gaps of the quasi-particle spectrum [58], shown in Fig. 5 (d).In addition, the density of states (DOS) is spin-polarized, namely, D ↑ ̸ = D ↓ at finite energy ω, as illustrated in Fig. 5 (d).
To summarize, we have derived a general weakcoupling criterion to investigate the spin-singlet orbitaltriplet pairings in nematic SCs.For technical convenience, we adopt the orbital d o -vector to describe the spin-singlet orbital-dependent pairing states and the g ovector for the orbital hybridizations.The main results of this work include, first, we demonstrate that an orbital d o -vector that is parallel with g o -vector for orbital hybridizations is possible to be realized in real superconducting materials.Second, the interplay between intraorbital and orbital-dependent pairings that belong to the same symmetry representation can explain the observation of robust Dirac nodes in the quasi-2D iron-based SCs.Remarkably, we find that d o -vectors could even coexist with many-body interaction-induced nematic orders or charge-density-wave orders when d o ∝ g tot = g o + g nem (or g int ).Moreover, our theory discovers the important role of nematic orders in SC pairing symmetry, which builds a possible bridge between repulsive interaction-induced nematic orders and nematic superconductivity and also reveals a deep connection between spin-singlet orbital-triplet pairings in nematic SCs.Our results may be helpful in understanding the nematic superconductivity in bulk FeSe.Our work will motivate more theoretical and experimental efforts to search for spin-singlet orbital-triplet SCs, even for topological superconductivity, which might contribute to further understanding the effects of spontaneous symmetry breaking on superconductivity.

Methods
In this section, we present the derivation for the main result Eq. ( 8), which is first order in λ o , by solving the linearized gap equation.The second order results are presented in Appendix C. The general k • p normal Hamiltonian considered in the main text reads, where the electronic basis is made of {1, 2}-orbitals x +k 2 y )/2m−µ is the band energy measured relative to the chemical potential µ, λ o represents the orbital hybridization and Besides, we set λ o > 0 without loss of generality.The Matsubara Green's function for electrons is Here β = 1/k B T and ω n = (2n + 1)π/β with n integer.Therefore, where We expand the attractive interactions as where v Γ 0 > 0 is the interaction strength of the irreducible representation channel Γ of the crystalline group, and l = 1, 2, ..., Dim Γ.Each pairing channel Γ gives rise to a SC critical temperature T Γ c , and the actual transition temperature of the system is given by the largest of these critical temperatures.In our work, we main focus on the case where Dim Γ = 1, which is sufficient for the applications discussed in the main text.The coupling between orbital-dependent pairings and orbital-independent pairings will be discussed in details later.The transition temperature T Γ c of orbital-dependent pairing channels is calculated by solving the linearized gap equation, which is reduced to v Γ 0 χ Γ (T ) − 1 = 0 with the superconductivity susceptibility χ Γ (T ) in the channel Γ defined as, where α, β ∈ {+, −}.For notional simplicity, the superscript Γ will be dropped when there is no danger of confusion.Firstly, let us calculate the trace part.In the following calculate, we will use And, Therefore, we arrive at Then we have Next, we calculate the integration for k,ωn by using, where N 0 is the density of states at Fermi surface and Ω is the solid angle of k on Fermi surfaces.Then, On one hand, where the approximation is done at low temperature when β → ∞.
On the other hand, we could find a series representation for χ 0 , which also applies to the case where λ o ̸ = 0, so that χ 0 ≡ χ(λ o = 0) and χ(λ o ̸ = 0) can be related by a simple relation.The way to do it is to perform the integration in ϵ first.More precisely, where the low temperature limit is again assumed and the integration is done using residue theorem.In the same spirit, we have, Now by introducing the digamma function defined on the complex plane, we have the following relation, where Therefore, Now we can proceed to calculate χ(T ) given in Eq. (33), In the calculation, we use normalized gap functions with where T c0 is T c for λ o = 0 case by solving v 0 χ 0 (T c0 ) = 1.This is the Eq. ( 8) in the main text.In general, the righthand side of Eq. (C13) suppress T c .It clearly indicates that T c would not be suppressed by orbital hybridization once d o ∥ g o for all k.So we conclude that the orbital d o -vector is possible to be stabilized in materials.
A. Two definitions for the orbital do-vector In the main text, we take the general pairing potential of a two-orbital SC, where ∆ s and ∆ o are pairing strengths in orbital-independent and orbital-dependent channels, respectively.Here τ are Pauli matrices acting on the orbital subspace and τ 0 is a 2-by-2 identity matrix.In the absence of band-splitting caused by spin-orbital couplings, the gap function on the Fermi surface is where This expression is mathematically similar to the superconducting gap of non-unitary spin-triplet SCs.
At this point, it is a good place to comment on the other possible way to defining the orbital d o -vector.Different from the one used in the main text, this definition groups the pairing term into orbital-singlet and orbital-triplet parts.In the form of Eq. (A1), Ψ s (k) and d 1,3  o (k) are even in k, but d 2 o (k) is odd in k due to Fermi statistics.By regrouping the terms based on the parity in k, we have which contains do • τ with the new do -vector redefined in terms of the original amplitudes and form factors.In this form, the first part is odd in k, which is the orbital-singlet part, and the second part is even in k and gives orbital-triplet state.Table II gives a detailed comparison between the two definitions of the orbital d o -vector.The spin d s -vector is also presented for completeness.It shows that the definition of orbital d o -vector used in the main text is more convenient to discuss the spontaneous TRS-braking pairing states.
Comparison between the two possible definitions of the orbital do-vector in spin-singlet SCs, together with the spin ds-vector of spin-triplet SCs.The parity properties are obtained from Fermi statistics.The TRS row gives the transformation properties in order to preserve TRS.Both the atomic orbital polarization (AOP) and the spin polarization (SP) take the same form in terms of their respective d-vectors.

B. Classification of spin singlet pairing states with Cn and TRS
In this section, we classify the possible spin-singlet pairing states constrained by C n about z-axis and TRS.And we also discuss the spontaneous time-reversal symmetry breaking pairings and the induced orbital polarized densityof-states.

Classification of pairings
The pairing potential ∆(k) transforms under the rotation C n as where J labels the irreducible representations of the C n point group.For example, J = 0 is for A representation and J = 2 is for B representation.Firstly, the TRS requires the coexistence of ∆J and ∆−J with equal weight.If the rotation symmetry C n is further imposed, then J and −J have to be equivalent modulo n, i.e.J ≡ −J mod n.
The results for the basis functions of Ψ s (k) and d o (k) are summarized in Table (1) in the main text.Here, the k z -dependent pairing symmetries are also presented for completeness.However, such pairings are neglected in the main text where we mainly focus on 2D systems.When inversion symmetry is also present, it leads to the following constraints for different orbital basis, 1.)If the inversion symmetry is I = τ 0 σ 0 , two atomic orbitals have the same parity, it require that d 2 o = 0. 2.) If the inversion symmetry is I = τ 3 σ 0 , two atomic orbitals have opposite parities, it require that d 1 o = 0. 3.) If the inversion symmetry is I = τ 1 σ 0 , two orbitals are the valley indexes, it require that d 3 o = 0.

Spontaneous TRS-breaking orbital-polarization
Next, we study spontaneous TRS-breaking and its consequences for a two-band SC with {d xz , d yz }-orbitals.I = τ 0 σ 0 constrains the orbital d o -vector to be ( , the set of superconducting order parameters are given by {∆ s , ∆ o , d ≜ (d 1 , 0, d 3 )}.Furthermore, the orbital orderings can be characterized by where , determining the gap strengths by using with β, β ′ > 0, which determines d.In addition, there are two possibilities to achieve the spontaneously TRS-breaking states, which are described respectively by F 3 and F 4 , The F 3 term in Eq. (B3) helps to develop a relative phase difference between ∆ s and ∆ o of being ±π/2 when b 1 = 0 and b 2 > 0 [59].As for the b 2 < 0 case, the TRS-breaking is caused solely by the F 4 term in Eq. (B4).For example, γ 0 (T ) = γ 0 (T /T ′ c − 1) and T ′ c < T c , where T ′ c is the critical temperature for the spontaneous TRS-breaking inside the superconducting states.When T < T ′ c , the orbital d o -vector becomes complex, then it generates the orbital orderings as M o = −iγ 1 /α M (d × d * ), of which only the y-component breaks TRS, as illustrated in the main text (see Fig. 1).More precisely, M y o ∝ k,σ ⟨n σ,+ (k) − nσ,− (k)⟩.Here we define |±⟩ = |1⟩ + i|2⟩ for complex orbitals, thus M y o ̸ = 0 indicates the TRS-breaking orbital-polarization (OP), similar to the time-reversal-odd polarization of the Cooper pairs discussed in Ref. [60,61].
We next solve the Bogoliubov-de-Gennes Hamiltonian, The quasi-particle spectrum is plotted in the main text (see Fig. 1), where two distinct gaps appear.Then, we calculate the atomic orbital-polarized density of states (DOS), where u n κ,σ = u n dxz,σ − iκu n dyz,σ with κ = ± for d xz ± id yz orbitals, and δ(x) is the delta function.In Fig. 6, the numerical results helps to confirm a two-gap feature due to the spontaneous breaking of TRS, compared with the quasi-particle spectrum.Moreover, D + ̸ = D − at finite energy clearly indicates that the DOS is orbital-polarized, which is consistent with the GL analysis.The orbital-spin conversion would lead to the spin-polarized DOS [55].

C. The stability of orbital do-vector under crystal fields
In the Method section of the main text, we derived our main result up to first-order in the coupling λ o .Here, we first address the situation where multiple pairing channels with possibly different pairing strengths coexist.Then we show the second-order result in λ o .1. First-order result applied to multiple coexisting pairing channels We have the following first-order result, ln which is derived with the assumption that there is only one pairing channel.Here, we elaborate on a subtlety mentioned in the main text that might arise due to coexisting multiple pairing channels belonging to different 1d irreducible representations.In the weak-coupling theory and without orbital hybridization, the critical temperature for a particular channel Γ is simply obtained by solving the linearized gap equation v Γ χ 0 (T c0 ) = 1 and the solution is given by T This means that the critical temperature in each channel is solely determined by the strength of the pairing interaction in that particular channel.The leading instability channel has the largest pairing interaction, which determines the T c .However, when orbital hybridization is considered, the story could change.Depending on the relation between d Γ o and g o , some pairing channels will be suppressed more than the others.Therefore, the previous leading instability channel could become sub-leading in the presence of orbital hybridization.In Figure 7 we take the J = 0 and J = 2 representations under C 4 as an example, with the assumption that v J=0 > v J=2 .We see that T

Second-order approximated results
In this subsection, we consider the coupling between orbital-independent pairing (Ψ s (k)-part in Eq. (A1)) and orbital-dependent pairing (d o (k)-part in Eq. (A1)), and study the second-order approximated results for the above conclusion.
The attractive interaction is now decomposed as where v 0 is the interaction strength in the orbital-dependent channel and v 1 is the interaction strength in the orbitalindependent channel.And they belong to the same representation of symmetry groups, leading to the coupled linearized gap equation, where It leads to Considering the v 0 > v 1 case firstly, then, the bare T c of orbital-dependent pairings are larger than that of orbitalindependent pairings, we have from which, we define the total superconductivity susceptibility as, where χ(T ) has been calculated in the above subsection (see Eq. ( 45)), and the second part is the second-order correction.After tracing out the spin degrees of freedom, we have χ os (T ) = χ so (T ).Following the same procedure as in the first-order case, we have which would vanish if λ o = 0, i.e. no orbital hybridization, based on the definitions of G +/− e/h , which in turn, reproduces the first-order calculation above.At non-zero, but small λ o (λ o k 2 F < µ), χ os (T ) will also be small but non-zero.For convenience of discussion, we define where δ(T, λ o ) ∼ λ o k 2 F /µ would vanish at leading order (see Eq. ( 35)).Then we have With this, the total superconductivity susceptibility in Eq. (C7) becomes where δχ(T ) is the second-order correction due to the coupling between orbital-independent pairings (Ψ s (k)-part in Eq. (A1)) and orbital-dependent pairing (d o (k)-part in Eq. (A1)), Since we assumed v 0 > v 1 , i.e.T c0 > T s , then the actual transition temperature would be T c ∼ T c0 > T s , giving log(T c /T s ) > 0. As a result, the correction to the susceptibility is positive: δχ(T ) > 0. Then following the same procedure as in the previous section, we have ln where the first part is the first-order result (see Eq. (C13) or Eq. ( 5) in the main text; order as O(λ o k 2 F /µ)), and the second part is the second-order result (order as O((λ o k 2 F /µ) 2 )).Therefore, we conclude that, • The first part: it determines the direction of orbital d o -vector to be d o ∥ ĝo .Because of C 0 (T c0 ) ≤ 0 , once d o ∥ ĝo at any momentum k, the first part vanishes.
• The second part: it relates the form Ψ s (k) to the orbital d o -vector: Thus, the second part becomes maximum, leading to the maximal increase of T c0 .
Here we have used the fact, 2 reaches 1 only when i = j.
Next, we briefly discuss the case where v 0 < v 1 .The same result can be similarly argued.In this case, the dominant pairing channel is the orbital-independent pairing (T s > T c0 ), which can induce the orbital d o -vector via their couplings.Similar to Eq. (C6), we define the total superconductivity susceptibility for orbital-independent pairings, which leads to • The numerator: it determines the local magnitude of orbital d o -vector to be d o ∝ Ψ s (k)ĝ o .Thus, the numerator becomes maximum, leading to the increase of T c0 maximally.
Therefore, according to both Eq.(C13) and Eq.(C17), we conclude that the orbital d o -vector that is parallel with orbital hybridization g o -vector could be generally stabilized in real materials.And we find that which is shown in Eq. ( 6) in the main text.However, it has also a Z 2 phase ±, which can be pinned by taking higher order corrections into account.
D. Formation of the pairing near Fermi surface in band picture Here, we provide another perspective on the pairing in orbital channel near the Fermi surface (FS) by looking at the total free energy of the system in band picture, where the pairing amplitude is treated perturbatively.
In the presence of orbital hybridization or nematic order, the Hamiltonian without pairing is given by where λ is taken to be positive.The degeneracy in the orbital channel will be lifted, whereas the spin channel still has the double-degeneracy.Effectively, the vector g acts as a "Zeeman field" in the orbital space, and the pseudo-spin will be parallel or anti-parallel to the field for the two splitting levels.And we notice that More precisely, the two eigenstates of H 0 can be denoted by |E ± ⟩ ≡ |ĝ; ±⟩, where +/− means parallel/anti-parallel (eigenvalues of the symmetry operator ĝ • τ ).Please note that ĝ = g/|g|.And with λ > 0. Setting E ± = 0, it gives rise to two FSs (labeled as FS ± ) with energy splitting as 2λ|g(k)|, which is approximately as ∼ λk 2 F with respect to ϵ k F = µ.Now we consider the spin-singlet pairing part in the original basis, The BdG Hamiltonian is then given by where γ i are the Pauli matrices in particle-hole channel.Next, we consider weak-coupling limit (infinitesimal pairing strength, namely, ∆ o → 0) and we use the band picture to study the pairing Hamiltonian.For this purpose, we rewrite H 0 as where τ = ± is the band index and s is for spin.The unitary transformation matrix U (k) in the orbital subspace leads to the diagonalization of H 0 , In the limit λk 2 F ≫ ∆ o (i.e., band splitting is much larger than the pairing gap), the inter-band pairing is not energetically favorable in the weak-coupling pairing theory (i.e., attractive interaction is infinitesimal small).It means the inter-band pairing will be severely suppressed if we increase the orbital hybridization λ, consistent with the calculation in the main text (see Fig. From the above analysis, we conclude that the orbital d o -vector should be parallel with the orbital hybridization g. To determine the magnitude of the orbital d o -vector, we turn on the orbital-independent pairing ∆ s ̸ = 0. We also assume both pairing channels are small.Assuming the two FSs have DOS N ± near the FS (ignoring the momentumdependence if the FS is almost isotropic), then the total condensation energy per volume and per spin of the two intra-band pairings is given by

G. Spin and orbital magnetizations: Ms and Mo
Here we present the definitions of spin and orbital magnetizations in our mean-field analysis.The spin magnetization is defined as, See Appendix D for details).Even though the intra-orbital pairing and the orbital-triplet pairing belong to the same symmetry representation, the different k-dependencies of Ψ s (k) and d ∥ (k) can naturally lead to the anisotropic superconducting gap on the Fermi surface observed in experiments.

FIG. 2 .
FIG. 2. The application to iron-chalcogenide superconductors with/without linear Dirac nodes.In (a), the two-electron pockets around the M point.For zero spin-orbit coupling, vso = 0, (b) shows the phase diagram as a function of the intra-orbital pairing ∆s and the inter-orbital pairing ∆o.For the gap parameters represented by the green dot in (b), the nodal superconductor is exhibited in (c), where the eight dark red points represent the chiral symmetry-protected Dirac nodes.

15 FIG. 4 .
FIG.4.Fermi surfaces (FSs) at the K± valleys and the quasi-particle density of states (DOS).The C6 symmetric FS without inter-valley scattering is shown in (a) and its DOS with an isotropic s-wave pairing is given in (d).(b) shows C6-breaking FSs due to the inter-valley scattering, together with the nodal lines of nematic pairing.There are no nodes on the FSs and the corresponding DOS is shown in (e).(c) is similar to (b) but with chemical potential µ adjusted so that the nodal lines intersect the FSs, hence a V-shaped DOS is obtained as in (f).Parameters used are the following, α = 0.2, t1 = 0.15, t2 = 0.25, ∆o = 0.5.For (a) and (d) µ = 1.5, gint = 0; for (b) and (e) µ = 1.5, gint = 0.7; for (c) and (f) µ = 2.7, gint = 0.7.

FIG. 5 .
FIG. 5. Schematic diagrams for the TRS-breaking effects for non-unitary orbital-triplet superconductors (SCs) and spintriplet SCs.As for orbital-triplet SCs, (a) shows a complex orbital do-vector that spontaneously breaks TRS and results in the TRS-breaking orbital-polarization with Mo ∝ ido × d * o ; and (b) shows the quasi-particle spectrum along kx and the orbital-polarized density of states (DOS) D± with |±⟩ representing |1⟩ ± i|2⟩.For spin-triplet SCs, (c) shows the superconductivity-induced spontaneous spin-polarization with Ms ∝ ids × d * s ; and (d) shows the two distinct gaps of the quasi-particle spectrum along kx and the spin-polarized density of state Dσ with σ = {↑, ↓}.The gapped spectrum is plotted for kz ̸ = 0 and the node in DOS profile is due to the nodal line at kz = 0.
J=0 c starts higher than T J=2 c , but since d J=2 o is parallel g o and d J=0 o is not, T J=2 c is not suppressed by λ o whereas T J=0 c is suppressed and eventually becomes lower than T J=2 c .

3 𝒅 3 0FIG. 9 .
FIG.8.The coexistence of orbital do-vector and nematic order.In (a), each tnem corresponds to a particular form of the d-vector, which determines Tc based on Eq. (C13).The three curves correspond to three different values for the nematic order Φ.The Tc is not suppressed by the nematic order as long as do ∥ gtot, i.e. tnem = 1.The tnem = 1 case is further illustrated in (b), where it is shown that the magnitude of the nematic order does not change Tc (up to the order of approximation made in Eq. (C13)).(c) shows non-zero nematic order breaks the original C4 (red line) down to C2 (blue line).Here go = (k 2x − k 2 y , 0, 3kxky).

2 ,
⟨c † s1 (k)σ s1s2 c s2 (k)⟩, (G1)and the orbital magnetization is defined as,M o = k,s,a,b ⟨c † s,a (k)τ ab c s,b (k)⟩.(G2)More specifically, the orbital magnetization has the following components, c s,dyz + c † s,dyz c s,dxz ⟩, c s,dyz − c † s,dyz c s,dxz ⟩ = 1 2 k,s ⟨n s,dxz+idyz − ns,dxz−idyz ⟩, c s,dxz − c † s,dyz c s,dyz ⟩, (G5) where M x,z o breaks the C 4 rotation symmetry and M y o breaks the TRS.In our work, we focus on spontaneous TRS breaking, with C 4 preserved.M y o (k) gives the local density difference in mementum space for d xz ± id yz orbitals.A non-zero M y o (k) implies TRS breaking at k, because the two orbitals are TR partners.However, it has to be noted that, similar to the spin-triplet case, local TRS breaking does not necessarily imply the global TRS is broken.To obtain the total overall magnetization, we still need to sum over momentum around the FS.In the Ginzburg-Landau formalism, the Im[d o × d * o ], whose only non-zero component is the y-component based on symmetry constraints in our formalism, is coupled to the induced magnetization M y o so that TRS is still retained at the Lagrangian level.Therefore, we have the following relevant terms αM y o Im[(d o × d * o ) y ] + β (M y o ) which upon functional derivative with respect to the induced magnetization would give M y o ∝ Im[(d o × d * o ) y ] = −i (d o × d * o ) y .