Pairing states of spin-3/2 fermions: Symmetry-enforced topological gap functions

We study the topological properties of superconductors with paired $j=\frac{3}{2}$ quasiparticles. Higher spin Fermi surfaces can arise, for instance, in strongly spin-orbit coupled band-inverted semimetals. Examples include the Bi-based half-Heusler materials, which have recently been established as low-temperature and low-carrier density superconductors. Motivated by this experimental observation, we obtain a comprehensive symmetry-based classification of topological pairing states in systems with higher angular momentum Cooper pairing. Our study consists of two main parts. First, we develop the phenomenological theory of multicomponent (i.e., higher angular momentum) pairing by classifying the stationary points of the free energy within a Ginzburg-Landau framework. Based on the symmetry classification of stationary pairing states, we then derive the symmetry-imposed constraints on their gap structures. We find that, depending on the symmetry quantum numbers of the Cooper pairs, different types of topological pairing states can occur: fully gapped topological superconductors in class DIII, Dirac superconductors and superconductors hosting Majorana fermions. Notably, we find a series of nematic fully gapped topological superconductors, as well as double-Dirac superconductors with quadratic dispersion. Our approach, applied here to the case of $j=\frac{3}{2}$ Cooper pairing, is rooted in the symmetry properties of pairing states, and can therefore also be applied to other systems with higher angular momentum and high-spin pairing. We conclude by relating our results to experimentally accessible signatures in thermodynamic and dynamic probes.


I. INTRODUCTION
In condensed matter physics, the study of superconductors has traditionally been guided by two defining characteristics of a bulk superconductor: the nature of the pairing order parameter and the mechanism of Cooper pairing [1][2][3][4]. Recent years, however, have witnessed great progress in understanding phases of quantum matter from the perspective of topology. In particular, in the case of superconductors, it has become clear that a global property of the Cooper pair wavefunction, encoded in its topology, constitutes a third defining characteristic. Nontrivial topology leads to the presence of quasiparticle excitations on surfaces and edges [5][6][7][8][9]. Specifically, the class of topological superconductors-in much the same way as topological insulators and topological semimetals-can be distinguished from ordinary superconductors by gapless quasiparticle excitations on the boundary, protected by the bulk superconducting gap structure. The latter is a manifestation of the bulkboundary correspondence, which establishes an inherent link between surface properties and bulk topology.
Topological superconductors with a bulk pairing gap are defined by a gap structure which cannot be adiabatically deformed into an s-wave superconductor without closing the pairing gap. Evidently, this implies that the question of pairing symmetry and bulk topology are inti- * jwfv@sas.upenn.edu mately related. Indeed, time-reversal invariant topological superconductors in class DIII are known to require odd-parity pairing [10,11]. The close connection between unconventional pairing symmetry and bulk topology is also manifest in topological nodal superconductors, i.e., superconductors with topologically protected nodal degeneracies in the bulk quasiparticle spectrum and distinctive gapless excitations at the surface [12]. A famous example of the latter are the topological bulk point nodes of the superfluid 3 He A-phase, which originate from the time-reversal breaking chiral pairing [13]. Therefore, superconductors with unconventional pairing symmetry generally inspire the question whether they realize topological pairing states, and thus have protected Andreev surface states.
In this paper we address this question for superconductors with a semimetallic normal state characterized by quadratically dispersing spin-orbit coupled j = 3 2 bands. An important motivation for this undertaking is the observation of superconductivity in the class of Bi-based half-Heusler materials APtBi and APdBi, where A can be a rare-earth element or Y/Lu [14][15][16][17][18][19][20][21][22][23][24]. Experimental evidence, in particular peneration depth measurements reported in Ref. 24, has given indications that the pairing in YPtBi is unconventional.
An additional incentive to consider the interplay of unconventional pairing and topology in spin-orbit coupled j = 3 2 bands is the possibility of high-spin Cooper pairing. This was recognized in important papers focusing on a specific set of fully gapped pairing states [25,26] and on-site pairings [27]. Reference 27 in particular has set the stage for studying superconductivity in the half-Heusler compounds [28][29][30][31][32][33][34]; the present authors have investigated the pairing instabilities in the p-wave pairing channels [31].
Whereas previous work has focussed primarily on the question of pairing symmetry, specifically in the context of materials such as YPtBi, the aim of this paper is to provide a comprehensive topological gap structure classification of spin j = 3 2 pairing states. Such classification, which encompasses all pairing channels, is desirable for the practical purpose of interpreting ongoing and future experiments, and stands to enable important progress in identifying the nature of the pairing order parameter in j = 3 2 systems. We proceed in two main steps. First, for multicomponent pairing channels, i.e., channels of Cooper pairing with nonzero total angular momentum, we obtain the stationary points of the free energy within a Ginzburg-Landau expansion using a symmetry-based strategy developed for 3 He [35]. Since only pairing states corresponding to stationary points can be minima of the free energy, any analysis of gap structure can be limited to this set of possible superconducting ground states determined by energetics. This, in practice, is a significant simplification. The second step is then to systematically analyze the topology of gap structures of each stationary pairing state by deriving the constraints imposed on the superconducting gap function by discrete (e.g., timereversal, inversion, mirror) and rotational symmetries.
Notably, we find a series of fully gapped topological superconductors which spontaneously break rotation symmetry and have a nematic axis [36,37]. In addition, we obtain different classes of point nodal superconductors, hosting low-energy Dirac or Majorana bulk quasiparticles with dispersion relations which depend on topological properties of the point node. Importantly, despite starting from a normal state with full rotational symmetry (emergent at low-energies), our formalism includes pairing states with discrete spatial symmetry and can therefore be viewed as including crystal anisotropy effects. Furthermore, since our approach relies on symmetry arguments, the results of our work are relevant to a broad range of spin-orbit coupled systems with higher angular momentum pairing. We begin by introducing the Hamiltonian of the normal state electronic structure close to the semi-metallic touching point at the zone center. We assume that other electron or hole pockets are absent. The normal state Hamiltonian is expressed as where c k = (c k 3 2 , c k 1 2 , c k,− 1 2 , c k,− 3 2 ) T are the j = 3 2 quasiparticle annihilation operators and h k takes the isotropic Luttinger form [38,39] h k = (κ 1 + 5 4 Here, m is an effective mass, µ is the chemical potential, and S = (S x , S y , S z ) T are the three spin matrices.
(Explicit expressions of the spin matrices are provided in Appendix A.) The Luttinger Hamiltonian describes a touching of quadratically dispersing bands at Γ which are spin-orbit split by the term (k · S) 2 . As a consequence of both time-reversal (Θ) and inversion symmetry (P ) the bands remain twofold degenerate at each momentum k. The Luttinger Hamiltonian can be diagonalized and brought into the form where the energies measured with respect to the chemical potential are given by and the operators f † k = (f † k↑ , f † k↓ ) and d † k = (d † k↑ , d † k↓ ) create quasiparticles in the energy eigenstates. The twofold degeneracy of each band defines an effective pseudospin degree of freedom, which we denote by ↑, ↓. Equation (4) shows that the coefficients κ 1,2 directly relate to the band curvatures. In this work we will particularize to the regime where κ 2 > κ 1 > 0. These conditions ensure that one pair of degenerate bands is electron-like and curving upward, the ε c k solution, and the other pair is hole-like and curving downward, the ε v k solution. We refer to these bands as the conduction band (c) and valence band (v), respectively. Furthermore, these conditions imply that the valence band states are | 3 2 , m j = ± 3 2 angular momentum states and the conduction band states are | 3 2 , m j = ± 1 2 states. (This may be seen by considering h k along k z .) In this way, by tuning the chemical potential, we have access to a valence band Fermi surface consisting of pseudospin ± 3 2 states and a conduction band Fermi surface consisting of pseudospin ± 1 2 states. Due to the different axial angular momentum of states on these Fermi surfaces pairing is expected to affect them differently. Therefore, in our study of pairing gap structures we clearly distinguish between valence band and conduction band Fermi surfaces. In Sec. IV, where we present the detailed analysis of gap structures, we will focus on both these cases, with a special emphasis on the more intriguing case of a pseudospin ± 3 2 Fermi surface. The operators f † k create quasiparticles in valence band eigenstates, i.e., f † kµ |0 = |k, µ; v , and similarly for d † k . It would be desirable for these pseudospin operators to transform as ordinary spin under spatial and timereversal symmetries. It is not guaranteed that such a basis for the band eigenstates exists. We can, however, choose a basis such that the pseudospin states, |k, ↑; v and |k, ↓; v (in case of the valence band), transform as canonical Kramers partners under time-reversal and inversion symmetry. This justifies the (pseudo)spin labeling ↑, ↓, and implies that, when considering the pairing on the Fermi surface, we can speak of pseudospin-singlet and pseudospin-triplet pairing.
In the Luttinger Hamiltonian of Eq. (2) we have neglected the terms originating from the crystal field. As a result, h k is invariant under continuous joint spatial and spin rotations. This approximation may be justified at low energies close to the touching point when crystal anisotropy effects can be considered small. Furthermore, and perhaps more importantly, for our purpose of a gap structure classification of stationary pairing states it is natural to choose a starting point of higher symmetry. In fact, as will be shown in Secs. III and IV, a gap structure classification developed on the basis of a rotationally symmetric model naturally includes the analysis of pairing states with discrete crystal symmetry, since the symmetry group which leaves stationary points of the free energy invariant may in principle be any subgroup of the full rotation group. Therefore, our gap structure classification also applies to pairing states arising in cubic models (see Sec. IV G). This implies, for instance, that our study bears a direct connection to the Bi-based half-Heusler superconductors, in particular YPtBi. We do note, however, that in this work we consider inversion symmetric systems with a twofold degenerate Fermi surface.
We note in passing that in this paper we exclusively focus on superconductivity, and assume the presence of a Fermi surface due to hole or electron doping; other ordering instabilities, relevant at the touching point, have been addressed in Refs. 40-43.

B. Pairing channels and their symmetry
The first step towards an analysis of pairing states is the identification of distinct pairing channels. Given the symmetry group G of the normal state material, the irreducible pairing channels are classified by the representations of G. In the present case, the Luttinger Hamiltonian of Eq. (2) has both time-reversal and inversion symmetry, and is invariant under joint rotations of spatial and spin degrees of freedom. As a result, including U (1) charge conservation, the symmetry group can be written as G = U (1) × SO(3) × P × Θ. The irreducible pairing channels can be distinguished by the angular momentum quantum numbers of the Cooper pairs. The total angular momentum J is the sum of the Cooper pair orbital angular momentum L and spin angular momentum S; as a consequence of spin-orbit coupling the symmetry quantum numbers of the Cooper pairs are (L, S; J, M J ), where M J is the magnetic quantum number describing the axial angular momentum.
Pairing channels with nonzero J have 2J + 1 independent components, transforming as partners under rotations, and are called multicomponent channels. The components are degenerate right at the superconducting transition temperature T c : the transition temperature is a property of the channel and symmetry requires the symmetry-related components to have the same T c . Our study of quasiparticle spectra and gap structures will require explicit expressions for these gap function components; they can be obtained using the standard L, S-coupling scheme for addition of angular momenta [25,31], as we will now briefly describe.
In the case of j = 3 2 quasiparticles the total spin of the Cooper pair can take the values S = 0, 1, 2, 3. Cooper pairs in a total spin S = 0 and S = 1 state are conventionally called singlet and triplet pairing states; by analogy S = 2, 3 states can be called quintet and septet pairings. We define Π † SM S (k) as the creation operator of a pair of quasiparticles with momenta k and −k, in a state with total spin S and magnetic quantum number Here, the matrices S SM are the multipole matrices of spin j = 3 2 fermions [26] and the anti-symmetric matrix T = e iπSy plays the role of ≡ is y familiar from spin-1 2 pairing.
The S = 0 matrix is proportional to the identity, i.e., S 00 = 1/2, and corresponds to a rotationally invariant spin-singlet pairing. The S = 1 matrices transform as a magnetic dipole (i.e., pseudovector) and are given by linear combinations of the spin matrices S. The S = 2, 3 matrices are higher order multipole matrices, describing spin quadrupolar and octupolar pairing, respectively, and transform as rank-2 and rank-3 tensors. Together these matrices span the space of Hermitian 4 × 4 matrices. (A more detailed discussion of the multipole matrices, including an explicit construction, can be found in Appendix A.) The internal spatial structure of the Cooper pair is captured by the orbital part of the Cooper pair wave function and is given by the spherical harmonics Y LM L (k), wherek = k/|k|. In the familiar nomenclature, superconductors with orbital angular momentum L = 0, 1, 2, 3 are referred to as s-, p-, d-, and f -wave, respectively. The spin and orbital angular momenta of Cooper pairs are constrained by Fermi statistics: since the matrices S SM S are symmetric (anti-symmetric) for even (odd) S, S and L must either both be even or both be odd. To form irreducible pairings, let Π † JM J (k) be the operator which creates a pair of electrons in a state of total angular momentum J and M J . Such irreducible pair creation operators take the general form where now the momentum-dependent matrices J JM J (k) are linear combination of Y LM L (k) and the spin matrices S SM S . The appropriate linear combinations are uniquely determined by the Clebsch-Gordan coefficients; specifically one has J Combinations of (L, S) such that L + S = J The complex expansion coefficients ∆ M J define the multicomponent pairing order parameter; superconductors with total angular momentum J are described by a 2J+1component order parameter. Different order parameter configurations with equal norm generally define different pairing states, with different symmetry properties (unless, clearly, the two configurations are related by a global rotation). A phenomenological theory of the superconducting order parameter will be developed in the next section, where ground state solutions of the free energy and their symmetry breaking patterns are discussed. In Table I we list a number of (L, S; J) pairing multiplets for orbital angular momenta up to L = 2. The pairing functions J JM J (k) corresponding to the latter pairing channels up to L = 1 are tabulated in Ref. 31. For wellscreened short-ranged interactions, the s-wave (L = 0) and p-wave (L = 1) pairing channels are expected to have highest T c , with higher angular momentum L channels being suppressed. The relative strength of s-wave and p-wave pairing instabilities is a question of high current interest, which we will not address here. For the purpose of this work we will take the position that both s-wave and p-wave pairing are likely to be relevant for the experimental systems under study.
It is useful to consider the symmetries of the pairing matrices J JM J (k) in more detail. In particular, the transformation properties under Θ and P will be of interest. First note that Θ and P act on the quasiparticle-operators as where T = e iπSy is the time-reversal matrix of Eq.
The spin multipole matrices transform under timereversal as ΘS SM S Θ −1 = (−1) S+M S S S,−M S , which can be derived from T S * T † = −S; for the spherical harmonics one has ΘY LM L (k) These symmetry properties, in combination with the transformation properties under SO(3) rotations, are at the heart of both the phenomenological theory of multicomponent pairing and the subsequent gap structure analysis. We note that the pairings J JM J (k) are eigenstates of rotations about the z-axis. If U θz is the j = 3 2 spinor representation of a rotation C θz about the z-axis, . We conclude this section with two remarks. The first concerns the symmetry of the normal state. When cubic crystal anisotropy effects are important, the irreducible pairing channels are labeled by representations of the cubic point group. The number of distinct (i.e., orthogonal) pairing channels in systems with discrete crystal symmetry is finite. Since the dimension of cubic representations is at most three, the pairings of Eq. (6) are generally split into pairings with distinct cubic symmetry. For instance, assuming an inversion symmetric normal state, the J = 2 channel is split into J = 2 → E g,u + T 2g,u , where g/u denotes the parity of the pairing channel, i.e., even/odd under inversion; the pairing channels J = 3, 4 of Table I are split as J = 3 → A 2g,u + T 1g,u + T 2g,u and J = 4 → A 1g,u + E g,u + T 1g,u + T 2g,u . The more precise splitting of the pairing components of Eq. (7) into pairing functions transforming as partners of cubic representations is tabulated in Table VII of Appendix F. For pairing channels up to L = 1 such splitting has been worked out in 31.
Furthermore, since the number of cubic representations is finite, pairings from distinct channels (L, S; J) will collapse onto the same cubic channel. For instance, as may be seen from the splitting of the J = 2, 3, 4 channels, the gap functions of cubic T 2g,u pairing can have contributions from all three isotropic channels. As a result, gap functions of cubic pairings will be linear combinations of symmetry-allowed terms with coefficients not determined by symmetry, and in general can be quite complicated.
Crucially, however, the symmetry group of a given pairing state, i.e., the subgroup of the normal state symmetry group which leaves a pairing state invariant, is manifest and independent of material-specific details. The symmetry group can be used to establish universal properties of gap functions which are independent of their specific form.
Second, we note that in an analysis of pairing one may choose to focus exclusively on the Fermi surface band (either valence or conduction band), and project out the "high-energy" band. Since the projected Fermi surface pseudspin operators only admit pseudospin-singlet and pseudospin-triplet pairing, S = 0 and S = 2 states cannot be distinguished on the Fermi surface, affecting the classification of irreducible pairings [31]. For instance, (L, S) = (0, 2) pairing effectively collapses onto (L, S) = (2, 0) pairing on the Fermi surface [27,31]. In this work, in order to correctly capture the essential features of pairing gap structures, it is important to include the effect of pairing-induced coupling of conduction and valence band.

III. PHENOMENOLOGY OF MULTICOMPONENT PAIRING
The Hamiltonian of Eq. (8) describes Cooper pairing within an irreducible pairing channel (L, S; J). A particular pairing state is specified by the superconducting order parameter ∆ M , which carries all information on its symmetry properties. (We drop the subscript J, as M = M J will always correspond to J in what follows.) For multicomponent superconductors, the order parameter not only has an overall amplitude and phase, but also internal structure: different order parameter configurations generally correspond to distinct pairing states. These pairing states can be sharply distinguished by symmetry; multicomponent superconductors, in addition to U (1) charge conservation, break symmetries such as time-reversal symmetry or rotation symmetry.
The spontaneous breaking of symmetry due to the selection of a specific pairing state occurs below T c . (Right at T c all states in a channel are degenerate). Consider, for instance, the simple example of a two-component pwave order parameter (p x , p y ) in two dimensions. Below T c the system will select either one of two states as its ground state: a time-reversal even but anisotropic state, given by cos θp x + sin θp y (where θ parametrizes a family of states), or a chiral pairing state of the form p x ± ip y .
To find the order parameter configuration of the superconducting ground state, one minimizes the free energy of the superconductor, denoted F J . At temperatures below but close to T c , a phenomenological Ginzburg-Landau theory (GL) is applicable, and one may expand F J in powers of the order parameter and its gradients. The phenomenological GL expansion parameters determine the superconducting state below T c . (At zero temperature an expansion is no longer valid, and the minimiza-tion must rely on the full free energy.) In general, for multicomponent orders, analytical minimization of the energy functional can become challenging or even untractable, as the number of interaction parameters increases with the number of components. It is not clear that a full analytical solution can be found when the number of components becomes large.
A powerful and elegant alternative strategy to obtain the free energy minima is based on the observation that solutions corresponding to free energy minima typically have residual symmetry, i.e., they do not fully break the symmetry of the normal state. This has motivated the expectation that states with residual symmetry are primary ground state candidates. In fact, by deriving all states invariant under a subgroup of the full symmetry group (which may be continuous or discrete subgroups), it was demonstrated that it is possible to systematically identify stationary states of the energy functional [35,44], which may then simply be compared by directly computing the energy. Even though there is no general proof that this delivers all stationary states, in all known cases where an analytical solution is available, the result matches the energy comparison of stationary states [35,45,46,53].
Two kinds of stationary states with residual symmetry exist. Inert states are stationary points of the free energy independent of the precise form of F J [47][48][49][50][51], whereas noninert states depend on the interaction parameters of the energy functional, requiring knowledge of its precise form. Following the method presented in Ref. 51 in the context of spinor Bose-Einstein condensates, both types of stationary states can be obtained from a symmetry classification of order parameter configurations.
For the purpose of studying gap structures of multicomponent pairing states, which is the focus of this work, these considerations lead to the important conclusion that we can restrict to studying the class of stationary states. This is a significant simplification, since in practice (i.e., for paring channels with nonzero but small J), the set of stationary states is rather tractable. Furthermore, insofar as inert states are concerned, details of the energy functional are unimportant. In this section our aim is to describe how the stationary solutions of the free energy can be derived using symmetry principles, and address their symmetry properties. This will provide the foundation for our gap structure analysis in Sec. IV. We pay particular attention to the cases J = 1, 2, 3, 4. In this section, we also develop a Ginzburg-Landau theory for multicomponent pairing. This allows us to energetically compare the stationary solutions, and study how degeneracy lifting and spontaneous symmetry breaking can occur immediately below T c .

A. Symmetry properties and stationary pairing states
The symmetry of the superconducting order parameter ∆ M is determined by the pairings defined in Eq. (7). For instance, it follows from Eqs. (10) and (11) that the order parameter transforms under Θ and P as The transformation properties under SO(3) rotations are uniquely fixed by the total angular momentum J.
The symmetry of the order parameter can be made more transparent by adopting a representation for the pairing states which exploits the analogy with angular momentum states. Specifically, given the order parameter ∆ M , we may write the pairing state |∆ as In this state vector representation, the states |J, M are identified with the pairings J JM (k) of Eq. All pairing states given by |J, M have the special property that they are inert states of the free energy [49][50][51]: they are stationary points of the energy independent of its precise form. The states |J, M have a continuous isotropy group, where the isotropy group is defined as the subgroup of total symmetry group G which leaves the state invariant. In the case of the states |J, M the isotropy group is isomorphic to SO(2), i.e., the group of rotations about the quantization axis. Members of SO(3) not part of the isotropy group generate different but symmetry-equivalent and energetically degenerate states. We may therefore take the z-axis as the axis of continuous rotations, without loss of generality.
Additional inert states can be obtained by considering discrete isotropy subgroups of SO(3); the discrete subgroups of SO(3) are C n (cyclic group rotations by an angle 2π/n about a special axis), D n (dihedral group of C n and an additional orthogonal two-fold axis), O (point group of the octahedron), T (point group of the tetrahedron), Y (point group of the icosahedron). An example of such a state is given by |∆ D4 = |2, 2 + |2, −2 , which is a pairing state of a J = 2 superconductor with D 4 symmetry.
In the process of constructing pairing states with discrete symmetry of a specific pairing channel one may find that some states are not uniquely determined. To obtain a stationary state of the free energy one minimizes the free energy over the manifold parametrizing these states. Such a stationary state is a noninert state, as it depends on the form of the energy.
Before we proceed to the GL theory, it is worth pointing out that instead of the magnetic basis used in Eq. (13), one may write the pairing state |∆ in a "real" basis as Here, |J, a are chosen such that time-reversal simply acts as |J, a → |J, a , which implies for the order parameter Θ : ∆ a → ∆ * a . For instance, for J = 1 one has a ∈ {x, y, z} and hence ∆ = (∆ x , ∆ y , ∆ z ) T ; for J = 2 one has a ∈ {x 2 − y 2 , 3z 2 − r 2 , xz, yz, xy}. In this basis, timereversal symmetry breaking pairing states are defined by order parameters which are not equal to their complex conjugates (up to a phase), which connects to standard treatments of multicomponent superconductivity. Note also that the rank-J tensor structure of ∆ is particularly transparent in this basis.
B. Ginzburg-Landau theory for general J For general angular momentum J, the free energy F J of the superconductor is an integral over the free energy density f J , In the GL regime (i.e., in the vicinity of T c ), where the strength of pairing is small, the free energy density f J can be expanded in powers of the order parameter ∆ and its gradients. At given order, the expansion consists of all independent terms invariant under the symmetries of the normal state. For our purposes it is sufficient to consider the homogeneous part of the free energy density and ignore contributions from spatial inhomogeneities. Up to fourth order in ∆, the free energy density can be expressed in the following general form where where n(ε F ) is the density of states at the Fermi energy.) Clearly, the first two terms only depend on the overall magnitude of ∆, and do not depend on the internal structure of the order parameter. (|∆| 2 is the only symmetry-allowed term at second order.) The third term is a sum over the magnitudes of the subsidiary order parameters I KN , where K is an angular momentum and N the axial angular momentum. Subsidiary order parameters are bilinears (i.e., composites) of the superconducting order parameter ∆ and capture the broken symmetries of the superconducting state. In the present case, the subsidiary orders I KN describe the magnetic multipole moments of the superconductor. For instance, superconductors with nonzero I 1,N =1,0,−1 have a magnetic dipole moment and thus have a chirality; superconductors with nonzero I 2,N =2,...,−2 have a magnetic quadrupole moment. It follows that the subsidiary order parameters encode the symmetry properties of the superconductor, and are thus sensitive to the internal structure of ∆. As a result, the GL coefficients v K are responsible for energetically discriminating different pairing states below T c , and they favor (or disfavor) pairing states with a certain structure of multipole moments. As a result, the GL analysis of multicomponent superconductors is an analysis of the terms with interaction parameters v K . We can write the subsidiary orders I KN as where I KN are the corresponding multipole matrices of an angular momentum J, whose dimension (2J + 1) × (2J + 1) therefore depends on J. The structure of the matrices I KN is discussed in more detail in Appendix B.
(Note that the I KN are gauge invariant, as is required for magnetic multipole order parameters.) As is clear from Eq. (17), in the case of superconductors with Cooper pair angular momentum J, one can form 2J distinct subsidiary order parameters. Importantly, however, the sum over K in Eq. (16) can be restricted to K = 1, . . . , J.
A proof of this, which builds on Ref. 47, is provided in Appendix C. (Note that K = 0 can be excluded as it would simply give another term ∝ |∆| 4 .) Since the subsidiary order parameters fully encode the symmetry properties of multicomponent pairing states, they can be used to uniquely distinguish classes of pairing states. More precisely, if two order parameter config-urations represent the same pairing state, they are characterized by the same pattern of subsidiary orders, up to a global rotation. These states have the same structure of multipole moments. If, on the other hand, two stationary states of the free energy have distinct isotropy groups, and are therefore different pairing states, this will be reflected in their multipole moment signature; they will be associated with a different structure of subsidiary order.
To establish a more concrete connection between the stationary states and subsidiary orders, consider, for instance, the pairing states |J, M . The chiral pairing states |J, M = 0 have a nonzero magnetic dipole moment proportional to M along the z-axis. To see this, let us define as the magnetic dipole moments along the x, y, z axes, which is often referred to as chirality. These are related to {I 11 , I 10 , I 1−1 } as I z = I 10 and I 1±1 = ∓(I x ±iI y )/ √ 2. Then, the chiral pairing states |J, M = 0 are characterized by I z ∝ M .
Clearly, the pairing states |J, M = 0 have vanishing chirality. In fact, all time-reversal invariant pairing states must have vanishing magnetic multipole moments with odd K. Pairing states which do break time-reversal symmetry but have vanishing chirality belong to an exotic class of states with nonzero higher order odd-K multipole moment. The rotational symmetry breaking of the pairing states |J, M = 0 is reflected in a nonzero quadrupole moment (and, in general, higher order even-K multipole moments). The quadrupole moment must be invariant under the rotations about the polar or nematic axis.
We conclude the general discussion of GL theory with a remark regarding remaining degeneracies of the functional (16). In Eq. (16) free energy density is expanded up to fourth order. While this is often sufficient to determine the ground state below T c , it is possible that degeneracies (i.e., degeneracies of states with distinct symmetry) remain at this order, which only get lifted at sixth or higher order. If two symmetry distinct states are found to have the same energy at fourth order, one must take the GL expansion to the next order to find the ground state. Within the phenomenological GL theory, higher order terms (which are not simply products of lower order terms) are systematically and straightforwardly constructed using the multipole moment subsidiary order parameters. For instance, sixth order invariants are simply obtained by considering products I K1N1 I K2N2 I K3N3 with N 1 + N 2 + N 3 = 0 and summing with the appropriate Clebsch-Gordan coefficients to generate total singlets. C. Examples: Application to J = 1, 2, 3, 4 Let us consider some examples, starting with the simplest case J = 1. Writing the order parameter in the real basis of Eq. (14) as ∆ = (∆ x , ∆ y , ∆ z ) T , the free energy density f J=1 is given by The superconducting state below T c is controlled by a single GL interaction coefficient v 1 , giving rise to two possible ground states [52]. The first is a non-magnetic polar state favored when v 1 > 0; the second, favored when v 1 < 0, is chiral and maximizes the magnetic dipole moment.
To see this more clearly, consider the chirality I of Eq.
The matrix elements of I x,y,z have a very simple form, given by As a result, one finds I = (I x , I y , In units where the pairing amplitude is set to 1, i.e., |∆| 2 = 1, the solutions of Eq (19) are given by I = 0 (v 1 > 0) and |I| = 1 (v 1 < 0). In terms of Eq. (13) these solutions are simply expressed as the polar pairing state |1, 0 and the chiral pairing state |1, 1 . (Note that |1, −1 is related to |1, 1 by a twofold rotation about the x-axis.) Next, we turn to the five-component J = 2 superconductor. The corresponding five-component order parameter is defined as and the free energy density f J=2 takes the form where now v 1,2 are two GL interaction coefficients. This GL energy functional was analytically solved by Mermin [53], who studied a purely orbital L = 2 d-wave order parameter, and was reconsidered by Sauls and Serene in the context of an (L, S; J) = (1, 1; 2) (or 3 P 2 ) order parameter for massive neutron stars [54].
To consider the free energy minima, it is convenient for the present purpose to follow the symmetry classification approach of Ref. 51. Of the stationary states of f J=2 only a subset of four corresponds to minima; the remaining states are saddle points. Two of these minima are inert states with a continuous symmetry group SO(2). They are given by |∆ 2 = |2, 2 and |∆ 0 = |2, 0 , or alternatively, in the real basis of Eq. (21) by and Whereas the former is chiral and has I = (0, 0, I z ) T = (0, 0, 2) T , the latter state is nematic with vanishing dipole moment.
The two remaining states corresponding to free energy minima have discrete symmetry: the state |∆ T = (|2, 2 + |2, −2 )/2 + i|2, 0 / √ 2 has tetrahedral symmetry and the state |∆ D4 = (|2, 2 + |2, −2 )/ √ 2 has dihedral D 4 symmetry. In the real basis these take the form and respectively. Neither of these two states has a magnetic dipole moment: I = 0 for both states. The tetrahedral state nevertheless breaks time-reversal symmetry, as is signaled by the relative phase in Eq. (25), and this is reflected in a nonzero octupole moment I 3N . The dihedral state is time-reversal invariant and its lowest nonzero multipole moment is quadrupolar. Even though they have different symmetry, the nematic state (24) and the dihedral state are similar in the sense that they both are time-reversal invariant and have nonvanishing quadrupole moment I 2N . In fact, these states are known as the uniaxial and biaxial nematic states, and they remain degenerate in energy to fourth order in the GL expansion [53]. The lifting of this degeneracy occurs at higher order in the expansion. That such a degeneracy lifting should occur is expected from the gap structures of these states, as we will demonstrate in the next section. Within weak-coupling BCS theory, it was found that the uniaxial nematic pairing is favored below T c [32].
Finally, we briefly discuss the cases J = 3 and J = 4. The free energy densities follow directly from Eq. (16) and are straightforward generalizations of Eqs. (19) and (22). Here, we will not quote their expressions explicitly, but instead focus our attention on the set of stationary states in these higher angular momentum channels. Clearly, the pairing states |3, M and |4, M are stationary states with a continuous symmetry group. In addition to these, a number of stationary states with discrete symmetry and total angular momentum J = 3 and J = 4 exist. To illustrate this, let us consider the inert stationary states with discrete symmetry. For J = 3 there are two such states, with octahedral O and dihedral D 6 symmetry, given by In the case of J = 4 pairing, there are one octahedral, one tetrahedral, and three dihedral states, given by (Note that we have not been concerned with normalization.) The tetrahedral state has the same symmetry as the state of Eq. (25). Except for this tetrahedral state, all these inert states correspond to time-reversal invariant pairing states. Apart from inert states, the J = 3, 4 pairing channels admit a number noninert states, which we will not list exhaustively here. An example will be considered in Sec. IV, where we present our gap structure analysis and classification; in particular, the gap structures of all inert states listed here will be considered. Insofar as the energetics of these J = 3, 4 pairing states with discrete symmetry is concerned, we make a general observation. Since octahedral (i.e., cubic) symmetry forbids a quadrupole moment, and time-reversal symmetry forbids both a dipole and an octupole moment, the octahedral states (16)]. For the case J = 3, this implies that the octahedral pairing state minimizes the GL free energy when v 1,2,3 > 0. Similarly, for the case J = 4 it is possible to show that the octahedral pairing state minimizes the GL free energy when v 1,2,3 > 0 and v 4 < 0, since the octahedral state maximizes the total hexadecapole moment N |I 4N | 2 . This observation can be viewed as an example of the general utility of expressing the energy functional in terms of quantities directly reflecting the pairing symmetry.

IV. QUASIPARTICLE GAP STRUCTURES AND TOPOLOGY
Based on the analysis of possible pairing ground states, in this section we turn to a detailed analysis of their quasiparticle gap structures, where we focus in particular on the associated topological quantum numbers. In three dimensions, four generic types of pairing gap structures can be distinguished: (i) a full pairing gap; (ii) bulk point nodes (codimension-3 nodes); (iii) bulk line nodes (codimension-2 nodes); and (iv) Bogoliubov Fermi surfaces (codimension-1 nodes). Gap structures of the latter kind have recently been introduced in Ref. 30, where it was shown that these nodal degeneracies of codimension-1 are topologically stable in parity-even multiband superconductors with spontaneously broken time-reversal symmetry.
Bulk point nodes correspond to Berry curvature monopoles in momentum space and must therefore come in pairs of opposite monopole charge [55,56]. In super-conductors particle-hole symmetry (Ξ) imposes the constraint that a point nodal degeneracy at momentum K on the Fermi surface must have a partner at −K (e.g., the antipodal point on a spherical Fermi surface) with opposite monopole charge. If the quasiparticle spectrum consists of a single pair of point nodes, or more generally multiple non-degenerate pairs, the low-energy gapless quasiparticles obey the Majorana equation of motion and realize itinerant Majorana fermions in three dimensions [37,57,58].
Different nodal gap structures arise when symmetries force pairs of point nodes to be degenerate. For instance, when time-reversal symmetry Θ is present each point node must be degenerate with a node of opposite Berry monopole charge. Such point nodes can be called Dirac points, by analogy with Dirac semimetals, realizing Dirac superconductors [59]. A second kind of degenerate point nodes occurs when the degenerate nodes have the same monopole charge, as is the case in the canonical example of the superfluid 3 He-A phase [13].
These general considerations demonstrate that the topological properties of the quasiparticle spectrum are inextricably linked to the symmetry of superconducting state, as symmetries can put constraints on the gap structure. In particular, the parity of the superconducting state plays a crucial role: time-reversal invariant topological superconductors (with a full pairing gap) must have odd-parity pairing [10,11]. Similarly, the parity of the pairing state is known to determine whether line nodes are stable degeneracies [60,61].
An analysis of gap structure topology must therefore clearly discriminate superconducting states with different symmetry. Accordingly, our derivation and classification of topological pairing states is built on the symmetry classification of stationary pairing states presented in the previous section.
The organization of this section reflects this. We begin by both reviewing and establishing a number of general implications of symmetry-mandated constraints on gap structures. Armed with these, we then carefully examine the gap structures of: single-component J = 0 superconductors (Sec. IV C), multicomponent M = 0 pairing states (Sec. IV D), multicomponent chiral pairing states (Sec. IV E), and, finally, pairing states with discrete symmetry (Sec. IV F).
To describe and study the quasiparticle gap structures of pairing states we adopt the mean-field formalism and define the Nambu spinor The superconducting mean-field Hamiltonian then takes the form with H k given by Here, h k is the Luttinger Hamiltonian of Eq. (2), ∆ k is the pairing potential, and we have introduced a set of Pauli matrices τ z and τ ± = (τ x ±iτ y )/2 acting on Nambu space. The pairing potential ∆ k follows from the pairing Hamiltonian in Eq. (8) and is given by As stated earlier, we will focus on the gap structures of order parameter configurations ∆ = (∆ J , . . . , ∆ −J ) T corresponding to the possible mean-field ground states which were obtained in the previous section; pairing states that do not correspond to free energy extrema are not considered. At this stage, it is useful to consider the discrete symmetry properties of the Hamiltonian H k . H k possesses a manifest particle-hole symmetry defined as Ξ = CK, where K is complex conjugation and C is a unitary matrix given by Specifically, the Hamiltonian satisfies Furthermore, depending on the parity of the orbital angular momentum L, the pairing potential is either even or odd under inversion, i.e., where P acts as the identity. For odd-parity pairing states, the inversion can be redefined as P acting as τ z , such that the Hamiltonian H k is inversion-symmetric, PH k P † = H −k . This implies, however, that P and Ξ do not commute, but instead satisfy the anticommutation relation {Ξ, P} = 0. A time-reversal symmetric pairing potential satisfies Since the pairing potential also obeys Fermi statistics, expressed as ∆ k = ∆ T −k , time-reversal invariance implies that the pairing potential is Hermitian: As a result, time-reversal invariant superconductors in three dimensions (Altland-Zirnbauer class DIII) admit a Z topological classification in terms of a winding number [5,62]. Any improper spatial symmetry, i.e., an inversion or mirror symmetry, forces the winding number to be zero, and this has lead to the important insight that time-reversal invariant topological superconductors in three dimensions must have odd-parity pairing symmetry [10,11]. In particular, when the Fermi surface (or, more generally, the set of disconnected Fermi surfaces) enclose an odd number of time-reversal invariant momenta, a fully gapped odd-parity superconductor is a topological superconductor. This is a powerful corollary which we can directly apply to the present case where we consider a single (either valence or conduction band) Fermi surface around Γ.
The mean-field Hamiltonian of Eq. (36) is expressed in the orbital basis; since we are interested in pairing on the Fermi surface it is advantageous rewrite it in the band basis, defined by the f k and d k operators. The quasiparticle operators c k and c † k can then be expressed in terms of f k and d k as where the matrices V k and W k contain the eigenvectors of the valence band and conduction band states, respectively. We choose a basis such that the Fermi surface pseudospin degrees of freedom transform under Θ and P as Θf kµ Θ −1 = µν f −kν and P f kµ P −1 = f −kµ . (See appendix D for explicit expressions.) Using Eq. (41) we rewrite H in the band basis as where ψ k and φ k are Nambu spinors for the band operators, i.e., The projection of H k onto the valence band, denoted H vv k , is given by and H cc k is simply obtained from the substitutions v ↔ c and V k ↔ W k . The off-diagonal blocks, which represent a pairing-induced coupling of the valence and conduction bands, take the form To describe pairing within the valence one can simply project out the conduction band and take Eq. (44). This ignores the effects of coupling to the conduction band captured by Eq. (45) and potentially misses qualitative features of the gap structure with topological origin [30]. Let us therefore take more formal approach which can account for all constraints imposed by the symmetry of the system. The resolvent corresponding to H k takes the form Using the properties of the resolvent, its valence band block is obtained as where the effective HamiltonianH vv k (ω) is given bỹ The poles of (47) still give the exact eigenenergies as long as the corresponding eigenstates have nonzero support on the valence band states. Since pairing is typically assumed to involve states on the Fermi surface, and one is interested in small energies ω compared to the Fermi energy, the effective Hamiltonian can be expanded to lowest order in ω/ε c , where ε c is the energy of the conduction band at momenta on the Fermi surface. Importantly, a number of properties of the effective valence band pairing Hamiltonian (48) can be established by simply invoking symmetry arguments. First, note that symmetries of the normal state, by definition, do not mix conduction and valence band states, implying that their action is block-diagonal in the (ψ k , ϕ k ) T basis. Second, note that Eq. (48) together with the poles of G vv (k, ω), implicitly defines the full effective pairing potential∆ k of the valence band, wherẽ ∆ k is a 2 × 2 matrix in pseudospin space. The symmetry properties of ∆ k [Eq. (37)] carry over to∆ k . In particular, given our choice of pseudospin basis, for even-/oddparity pairing one has∆ k = ±∆ −k , and time-reversal symmetric pairing implies ∆ * k T =∆ −k . In combination with Fermi statistics time-reversal symmetry leads to a Hermitian pairing potential∆ † k =∆ k . Then, in the case of odd-parity pairing,∆ k may be expanded in pseudospin Pauli matrices s x,y,z as∆ k = g(k) · s, where g(k) = g * (k) = −g(−k) is real. We would like to consider constraints imposed on g(k) by mirror symmetries. To this end we must establish how the valence band pseudospin states transform under mirror symmetry. Using the pseudospin basis presented in Appendix D, we find that the pseudospin matrix representations O Mx and O My of the mirror symmetries M x : x → −x and M y : y → −y take the simple form This proves that a pseudospin basis exists such that the pseudspin transforms as an ordinary spin under mirror symmetry, and immediately implies that mirror symmetry imposes constraints on the gap function. Specifically, on the mirror plane g must be normal to the mirror plane. Moreover, since the equation g x (k) = 0 has a one-parameter family of solutions on a yz mirror plane, the valence band gap structure has mirror symmetry protected point nodes where the solutions of g x (k) = 0 intersect the Fermi surface [60]. These points nodes must be degenerate due to time-reversal symmetry.
We have thus obtained the result that time-reversal invariant odd-parity pairing states with a mirror symmetry generically have degenerate point nodes (with opposite Berry monopole charge). The point nodes are protected by mirror symmetry [59,60,63].

A. Pairing states with a rotation axis
We have learned from Sec. III that pairing states of multicomponent superconductors generically have special axes of rotation symmetry. These may be principal axes of continuous rotations or (a set of equivalent) discrete rotation axes. At Fermi surface momenta which lie on the rotation axis, rotation symmetry can give rise to constraints on the gap structure, leading to point nodal degeneracies and gapless quasiparticle excitations [58]. Therefore, to study the gap structure of superconductors with a rotation axis, we develop a symmetry-based theory for the quasiparticle spectrum in the vicinity of the rotation-invariant Fermi surface momenta, which we denote ±K, see Fig. 1. In most cases we will be able to choose K along the z-axis, without loss of generality, in which case K = k F v,cẑ (Fig. 1 A). Here, K = k F v,c is the Fermi momentum of a valence or conduction band Fermi surface, and is given by In the case of a valence band Fermi surface (a case we will often consider as an example), we expand the Nambu spinor ψ k in small momenta q around ±K and define the spinor Ψ qµ as where we have introduced the label 1, 2 for the nodal degree of freedom ±K. Recall that for the valence band µ =↑, ↓ refers to the ± 3 2 pseudospin states. Here, the quantization axis is chosen along the rotation axis defined by K, such that under rotations by an angle θ one has C θ f K↑,↓ C −1 θ = e ±i3θ/2 f K↑,↓ . Similarly, in the case of a conduction band Fermi surface, we collect the conduction band degrees of freedom close to ±K in the spinor Φ qµ , given by Note that now, however, the pseudospin label µ =↑, ↓ refers to the ± 1 2 states, such that under rotations one has C θ d K↑,↓ C −1 θ = e ±iθ/2 d K↑,↓ . Expanded in these degrees of freedom, the Luttinger Hamiltonian of Eq. (2) takes the form depending on whether one is considering a valence band or conduction band Fermi surface. Here, h v,c q are given by with ε v,c ±q ≡ ε v,c K±q . Note that we have used the inversion symmetry of the normal state: ε v,c k = ε v,c −k . Now, let us specifically consider a valence band Fermi surface. The dispersion ε v ±q can be expanded in small q as whereK = K/k F v and v F v is the valence band Fermi velocity. If ±K is along the z-axis, this reduces to Note that the conduction band constitutes a high-energy subspace, located at energy ε c K = µ(κ 1 + κ 2 )/(κ 1 − κ 2 ). Next, consider the pairing at and close to ±K. We distinguish two cases: even-parity pairing and odd-parity pairing. Even-parity pairing states have total spin S = 0, 2 (see Sec. II B), which, when projected onto the valence band, implies pseudospin-singlet pairing. To describe even-parity pairing in the vicinity of ±K we define the pseudospin-singlet operator F † s (q) as Odd-parity pairing states have total spin S = 1, 3, and consequently have pseudospin-triplet structure when projected onto the valence band. In accordance, we define the three pseudospin-triplet operators as In addition to these pairing operators, in order to capture the full low-energy structure of pairing in the valence band, we define the following effective Zeeman-type spinsplit operator whereas for odd-spin pairing it is given by It follows from (59) that δ ∼ |∆| 2 /ε c K . The general program for the remaining parts of this section which pertain to pairings with a rotation axis is to derive constraints imposed by symmetry on effective low-energy gap functions. In appendix E we show how to obtain these gap functions from any particular pairing potential ∆ k [Eq. In the familiar case of spin j = 1 2 pairing, Cooper pairs can be in a spin-singlet (S = 0) or spin-triplet (S = 1) state (assuming parity is a good quantum number). The quasiparticle spectrum of a spin-singlet superconductor is manifestly two-fold spin-degenerate. In contrast, spin-triplet superconductors are either unitary or non-unitary. Non-unitary superconductors necessarily break time-reversal symmetry and have the property that the two spin species have different quasiparticle spectra, i.e., the spectrum is not manifestly degenerate. The latter can have important implications for the gap structure, as it is a necessary condition for non-degenerate point nodes to exist. We note in passing that the converse is not true: time-reversal symmetry breaking does not necessarily imply non-unitary pairing. In spin-orbit coupled systems, however, time-reversal symmetry breaking generically leads to non-unitary pairing, since symmetryallowed terms in the gap function (which reflect spinorbit coupling) generically give different pairing for spinup and spin-down.
The notion of degenerate and non-degenerate quasiparticle spectra can be generalized to pairing states of j = 3 2 fermions. We will call pairing states with a nondegenerate quasiparticle spectrum spin-selective pairing states, and states with a manifestly two-fold degenerate spectrum spin-degenerate pairing states. (Note that the distinction 'unitary' versus 'non-unitary' is specific to spin-1 2 pairing.) Clearly, spin-selective versus spindegenerate pairing pertains to the spin structure of the Cooper pairs. Spin-selective pairing states should be understood as states described by a pairing potential which contains S SM S and S S,−M S in an asymmetric way.
Since time-reversal symmetry guarantees a two-fold degenerate spectrum, spin-selective pairing states must break time-reversal symmetry. In general, the converse is certainly not true. As in the aforementioned case of j = 1 2 systems, however, spin-orbit coupling generically leads to spin-selective pairing when time-reversal symmetry is broken. To see this in the present context, one may consider the irreducible spin-orbit coupled pairings given by Eq. (7). A given pairing J JM with M = 0 is an asymmetric sum over orbital and spin angular momentum. For instance, the pairing J 31 = c 1 Y 11 S 30 +c 2 Y 10 S 31 +c 3 Y 1,−1 S 32 (c 1,2,3 are unimportant Clebsch-Gordan coefficients) contains S 31 and S 32 but neither S 3,−1 nor S 3,−2 . This implies different pairing for spin species related by timereversal symmetry and thus constitutes spin-selective pairing. DIII and may be viewed as higher spin analogs 3 He Bphase. They are, however, different from the 3 He B-phase in a crucial way, which depends on the structure of the Fermi surface. Both Refs. 25 and 26 have reported that in the case of a pseudospin ± 3 2 Fermi surface (a valence band Fermi surface, in the present context) the Z winding number characterizing the odd-parity J = 0 topological superconductor is ±3, making it topologically distinct form the 3 He B-phase. This distinction is expressed through the surface state spectrum.
The even-parity J = 0 pairing states are topologically trivial; Table II  As discussed in Sec. III A, the pairing states |J, 0 are time-reversal invariant (up to a phase) and therefore spindegenerate. As far as spatial symmetries are concerned, the states |J, 0 can be distinguished by two symmetry quantum numbers: the parity of L (i.e., even/odd under inversion) and the parity of J. The parity of J determines whether |J, 0 is even or odd with respect to twofold rotation about an axis perpendicular to the quantization axis. Without loss of generality we take the latter to be the z-axis. One then has that |J, 0 is even (odd) under a twofold rotation about an axis in the xy-plane when J is even (odd). This is why we may call even-J states nematic and odd-J states polar.
Importantly, since mirror reflections in a plane that contains the z-axis are the product of inversion and perpendicular twofold rotations, the parity of L and J also determine the mirror symmetry properties of |J, 0 . The mirror symmetry properties have important consequences for the gap structures of both even-parity and odd-parity pairing states.
We first consider odd-parity pairing (i.e., L odd). Then, the pairing states |J, 0 are even (odd) under mirror reflections in planes perpendicular to the xy-plane when J is odd (even). For instance, the state |2, 0 is odd under mirror symmetry. This directly leads us to an important observation: odd-parity pairing states |J, 0 with even J have neither an inversion symmetry nor a mirror symmetry, and, as a result, there are no constraints on the gap function which might enforce nodal degeneracies. We conclude that these pairing states have a full pairing gap on the Fermi surface and are time-reversal invariant topological pairing states in class DIII. Notably, these topological superconductors are different from the J = 0 superconductors mentioned in Sec. IV C, since the former break rotation symmetry and have a nematic axis.
Due to the presence of a mirror symmetry, the gap structure of odd-parity pairing states |J, 0 with odd J is different. As discussed earlier in this section, see also Eq. (49), mirror symmetries force the pairing gap to vanish at points on the Fermi surface. In the present case, since the |J, 0 have a continuous rotation symmetry about the z-axis, the point nodes are located on the z-axis, i.e., at ±K = k F vẑ , see Sec. IV A. To demonstrate this in more detail, we rearrange the spinor components of Ψ q↑,↓ in Eq. (51) and define Ψ q1,2 as The low-energy pairing hamiltonian near ±K, defined in Eq. (61), can be expressed as where we recall that τ ± = (τ x ± iτ y )/2, and ∆ q is given by Here, as before, s z and s ± = (s x ± is y )/2 are Pauli matrices in pseudospin space. Rotation symmetry and mirror symmetry impose constraints on the three functions {∆ q+ , ∆ q0 , ∆ q− }. [Note that the action of mirror symmetries on the pseudospin degrees of freedom was determined in Eq. (49).] In particular, under rotations C θz by an angle θ the spin matrices transform as C θz : s ± → e ±3iθ s ± ; s z is left invariant. As a result, to lowest linear order in q one finds ∆ q0 = 0 and ∆ q± ∝ ∆(q x ∓ iq y ) 3 , where ∆ is the pairing amplitude. This not only shows that the pairing gap vanishes on the rotation z-axis, but also proves that the quasiparticle dispersion is cubic in momenta q x,y in directions orthogonal to the rotation axis. Similarly, one may consider Eq. (64) for the pseudospin ± 1 2 Fermi surface, in which case the Hamiltonian (63) should be expressed in terms of spinors Φ q1,2 defined in analogy with Eq. (62). One then finds that ∆ q0 = 0 and ∆ q± ∝ ∆(q x ∓ iq y ), implying that the quasiparticle dispersion is linear in all directions away from ±K.
We can express this in terms of a Hamiltonian H q for the low-energy quasiparticles at ±K. Introducing a set of Pauli matrices σ x,y,z for the nodal degree of freedom, where q ± = (q x ± iq y ). A Hamiltonian of this form describes Dirac quasiparticles with cubic dispersion. We thus draw the important conclusion that the pseudospin ± 3 2 gap structure of odd-parity pairing states |J, 0 with odd J is given by triple Dirac points on the rotation axis.
For the pseudospin ± 1 2 states, on the other hand, the Hamiltonian for the low-energy gapless quasiparticles is obtained as where v ∆ is an effective velocity in the x, y directions, set by the pairing strength. Hamiltonian (66) shows that the low-energy quasiparticles are governed by a Dirac equation with linear dispersion.
As an example of an odd-parity pairing state |J, 0 which gives rise to Dirac points on the rotation axis one may consider the state |3, 0 given by Here, Y 1,±1 = Y 1,±1 (k) are spherical harmonics. Since We now come to even-parity |J, 0 pairing states. Once more, mirror reflections prove to be important. The evenparity pairing states are odd (even) under mirror reflections in planes perpendicular to the xy-plane when J is odd (even). First, consider odd J. States with odd J are odd under mirror reflections and this forces the gap function to vanish on any of these mirror planes. (Note that even-parity pairing states are pseudospin singlets.) Since the states |J, 0 have a continuous rotation symmetry about the z-axis, this implies that the pairing gap vanishes on the entire Fermi surface. As a result, the even-parity pairing states |J, 0 with odd J remain fully gapless.
Finally, the even-parity states |J, 0 with even J are mirror symmetric and generically have line nodes. The latter simply follows from the fact that the pseudospinsinglet gap function must have sign changes on the Fermi surface.
These results are summarized in Table III. It is important to point out that these results are independent of whether a valence band or conduction band Fermi surface is considered. In particular, Eqs. (63) and (64), and subsequent analysis, remain valid when applied to a conduction band Fermi surface. The key symmetry property of the chiral states is their eigenvalue of rotation about the quantization axis, which depends on the axial angular momentum M . As discussed in Sec. III A, we can take the z-axis as the rotation axis. Then, focusing on ±K = ±k F vẑ , in this section the aim is to derive the constraints on the gap functions ∆ qs and {∆ q+ , ∆ q0 , ∆ q− } of Eqs. (60) and (61) imposed by rotation symmetry.
Rotation symmetry imposes the constraint that the orbital angular momentum of the low-energy gap functions ∆ qs and {∆ q+ , ∆ q0 , ∆ q− } must match the axial angular momentum M . Here, we anticipate differences for a valence band and conduction band Fermi surface. The pseudospin-triplet operators F † t± of Eq. (58) carry angular momentum ±3, whereas in the case of a conduction band Fermi surface they carry angular momentum ±1, affecting the matching conditions.
As a result of time-reversal symmetry breaking, the chiral pairing states are generically spin-selective have nonzero δ in Eqs. (60) and (61). Within the class of spinselective chiral pairing states we can formulate a more precise constraint on δ by considering the full pairing potential ∆ k of Eq. (37) at ±K. If ∆ k vanishes at ±K, i.e., ∆ ±K = J JM J (±K) = 0, then δ must be zero. Since only Y L0 (±K) = 0, the pairing potential ∆ ±K must vanish whenever the corresponding Clebsch-Gordan coefficient vanishes, i.e., For the L = 1 multiplets this only occurs for (L, S; J) = (1, 1; 2) when M = 2.

Pseudospin-singlet Hamiltonian from symmetry
The pseudospin-singlet pairing Hamiltonian defined in Eq. (60) can be combined with normal state Hamiltonian Eq. (54) to obtain the full Hamiltonian of the low-energy quasiparticle degrees of freedom. Rather than Ψ q↑,↓ , it is convenient to rearrange the operators and form the spinors Ψ q± given by The full Hamiltonian is block diagonal in this basis, i.e., H ± = 1 2 q Ψ † q± H q± Ψ q± , with the Hamiltonian matrices H q± are given by The constraint on the gap function ∆ qs is that the angular momentum of the pseudospin-singlet pairing matches M . More specifically, the gap function must be of the form The Hamiltonian of Eq. (72) can be recast using Pauli matrices σ z = ±1 for the nodal degree of freedom ±K.
Equation (73) shows that ∆ q is either even or odd under q → −q. When the gap function is even, Eq. (72) can be expressed as  and when it is odd we find The spectrum takes the same general form on both cases. We find the four energy eigenvalues The most significant feature of these solutions is that they generically describe two nodes of codimension 1, one enclosing K and one enclosing −K. These nodes are defines by the solutions of the equation Bogoliubov Fermi surfaces of this type were first described in Ref. [30], which noted that they may be viewed as inflated point nodes. This interpretation naturally follows from the picture presented here. That Fermi surfaces may be topologically stable features of a gap structure follows from a topological Z 2 invariant associated with nodes of codimension 1 in even-parity time-reversal symmetry broken superconductors [61,65]. Additional Bogoliubov Fermi surfaces generically occur on the Fermi surface equator, i.e., in the vicinity k z = 0, of even-parity chiral pairings with odd M . For odd M , |J, M is odd under twofold rotation about the z-axis. Since we have assumed even parity, |J, M is also odd under a mirror reflection in the xy-plane. This would imply a line node on the equator, however, since the pairing is spin-selective these line nodes generically will be inflated to nodes of codimension 1.

Pseudospin-triplet Hamiltonian from symmetry
Pseudospin-triplet pairing must be of odd-parity type, and therefore nodes of codimension 1 (i.e., surfaces) are not topologically stable [61,65]. Since chiral pairing states are generically spin-selective, the effective pseudopsin-splitting δ is nonzero, implying that point nodes on the rotation axis (if they exist) are nondegenerate. This may be compared to non-degenerate nodal degeneracies in ferromagnetic superconductors, where the Zeeman splitting originates from ferromagnetic order rather than spin-selective pairing [66].
The splitting of pseudospin-↑ and ↓ implies that, in order to determine the symmetry-mandated low-energy dispersion of quasiparticles on the rotation z-axis, we need to examine the gap functions ∆ q+ and ∆ q− , since these correspond to ↑↑and ↓↓-pairing. In the spirit of Refs. 58, 67, and 68 one finds that constraints derived from rotational symmetry (i.e., the angular momentum quantum numbers) fully determine the form of ∆ q+ and ∆ q− . We demonstrate this by considering ∆ q+ . The gap function ∆ q+ can be expanded in momenta q x , q y perpendicular to the rotation axis as where and are nonnegative integers defining orbital axial angular momentum quantum numbers; A are coefficients. In terms of the quantum numbers , the orbital angular momentum of ∆ q+ is given by − . Furthermore, in the case of the (valence band) pseudospin ± 3 2 states the pseudospin angular momentum of ↑↑-pairing is 3 2 + 3 2 = 3. (The latter equals 1 2 + 1 2 = 1 for the conduction band ± 1 2 pseudospin states.) Since the total axial angular momentum of the pairing state |J, M is M , the sum of orbital and pseudospin angular momentum must be equal to M , and we thus arrive at the matching condition Since in the expansion of Eq. (77) we are interested in the lowest order terms, we only consider solutions of (78) for which either or is zero. For instance, in the case when M = 2, the matching condition gives ( , ) = (0, 1). Clearly, a matching condition similar to Eq. (78) exists for ∆ q− , in which case the pseudospin angular momentum is −3. Applying these matching conditions to the cases M = 1, 2, 3, we arrive at Table IV. An analogous analysis can be straightforwardly performed for a pseudospin ± 1 2 conduction band Fermi surface, in which case the pseudospin-triplet operators carry angular momentum ±1 and one should replace 3 with 1 in (78). The low-energy behavior of pseudospin ± 1 2 gap functions ∆ q± is summarized in Table IV as well. Table IV shows that special cases arise when M = 3 (for pseudospin 3 2 ) and M = 1 (for pseudospin 1 2 ) , i.e., when the angular momentum of the pairing state matches the angular momentum of the operator for pseudospin-↑↑ pairing. In this case, the pseudospin-↑ quasiparticles can pair at ±K and develop a pairing gap. The pseudospin-↓ quasiparticles must remain gapless, however. Therefore, when the angular momentum of the pairing state matches the angular momentum of the Cooper pair f † 1µ f † 2µ (µ =↑ or ↓), one of the two nodes along the rotation axis will be gapped out, leaving a single 3D Majorana fermion behind.
The low-energy spectral properties of chiral pairing states |J, M = 0 summarized in Table IV have been established based on symmetry arguments which take into account the pseudospin splitting δ implicitly. The low-energy quasiparticle gap structure of odd-parity chiral pairing states may also be obtained from an explicit calculation based on the Hamiltonian of Eq. (61). Specifically, using the spinors defined in Eq. (62), the pairing Hamiltonian H ∆ of Eq. (61) (which explicitly includes the pseudospin splitting proportional to δ) is expressed as where have suppressed the contribution from Ψ q2 since all spectral information is contained in (79). As in Eq. (63), the pairing potential ∆ q contains the three gap functions {∆ q+ , ∆ q0 , ∆ q− }. These gap functions as well as δ can be determined using the perturbative approach detailed in Appendix E. Upon including the normal state part of the Hamiltonian, to the lowest order given by ε v q τ z , we obtain the four branches E ±± q of the quasiparticle spectrum as where Λ q is defined as It is straightforward to establish that Λ q is only nonzero for chiral states and must be zero when time-reversal symmetry is present. To see this, note that time-reversal symmetry requires ∆ * q0 = ∆ q0 , ∆ * q+ = −∆ q− , ∆ * q− = −∆ q+ , and, as discussed above, δ = 0. It then simply follows that Λ q = 0 in this case.
Even though Eqs. (80) and (81) appear rather complicated, they describe low-energy gap structures whose essential properties have been rigorously determined from symmetry arguments, and are given by Table IV. What Eqs. (80) and (81) nevertheless serve to illustrate is the importance of the energy scale set by δ. In particular, as mentioned in Sec. IV A, δ describes a pairing-induced splitting of pseudospin states proportional to |∆| 2 /ε c K . Correspondingly, as may be checked directly from Eqs. (80) and (81), the two point nodes are separated by a momentum ∼ |∆| 2 /(v F v ε c K ) along the k z axis. The emergence of this energy scale is important for the potential observation of the nodal structure through thermodynamic probes, since it determines the temperature range over which the characteristic temperature dependence of thermodynamic quantities is accessible.
We conclude the discussion of odd-parity chiral pairing states by illustrating the general considerations with simple examples. We focus our attention on the chiral pairing states with M = 1, 2, 3 in a J = 3 channel with angular momentum quantum numbers (L, S) = (1, 3). The corresponding pairings J 3M are given by where we have suppressed the momentum dependence of the spherical harmonics. These pairings generate the lowenergy nodal structures listed in Table IV. To gain a better understanding of these example pairings, consider the terms proportional to Y 10 ∼ k z , which are nonzero on the rotation z-axis. In case of the J 31 pairing with M = 1, the spin matrix S 31 connects states which differ by one unit of angular momentum. Therefore, S 31 does not directly connect the pseudospin-3 2 states, which have relative angular momentum 3, but does connect the pseudospin-1 2 states. As a result, both pseudospin-3 2 species remain gapless, whereas the pseudospin + 1 2 states can pair, leaving only the pseudospin − 1 2 gapless. This qualitatively explains the first row of Table IV. A similar argument can be made in case of J 33 pairing: the spin matrix S 33 connects the pseudospin- 3 2 states, such that a pseudospin-↑↑ pairing can form on the rotation axis. This corresponds to the aforementioned special case of M = 3; see also Table IV. In contrast, the matrix S 32 does not directly connect any of the states within a pseudospin sector. In all these three cases, the terms proportional to Y 10 are responsible for finite δ, giving rise to the pseudospin splitting.

F. Gap structures of pairing states with discrete symmetry
In the final part of this section, we consider gap structures of pairing states with discrete symmetry. As mentioned in Sec. III A, even when the normal state has full SO(3) rotational symmetry, there can be-and typically will be-stationary states of the free energy with discrete spatial symmetry. In Sec. III C we have discussed a number of examples of such pairing states, focusing in particular on the inert states. Here, we examine their gap structures.
In order to do so, it is necessary to comment on the precise structure of the isotropy groups of these pairing states. We have briefly mentioned the definition of isotropy groups in Sec. III C; they are subgroups of the full symmetry group G which leave the state invariant. Importantly, this implies that elements of the isotropy group may be composites of spatial transformations and U (1) gauge factors. (We have implicitly made use of this in the case of the continuous isotropy groups of the |J, M states.) In the case of discrete isotropy groups such as O, T , and D n , it is particularly important to properly account for phase factors associated with spatial symmetries. Consider, for instance, the pairing state |2, 2 + |2, −2 of a J = 2 superconductor, which has D 4 symmetry. The two generators of the isotropy group D 4 are given by {e iπ C 4z , C 2x }, showing that the fourfold rotation leaves the state invariant only in combination with the phase factor e iπ . The significance of this for our purposes is that the precise structure of the isotropy group can depend on the total angular momentum of the superconductor. In particular, two pairing states with the same discrete symmetry may still have different gap structure due to a different realization of the isotropy group.
We further note that the isotropy groups are taken to be subgroups of U (1) × SO(3). For superconductors, the symmetry of the pairing state under inversion, and more generally under improper rotations, is fixed by the parity of the pairing state. Therefore, we treat even-and oddparity pairing separately.
Our goal in this section is to illustrate our classification by focusing primarily on inert pairing states with O, T , and D n symmetry; one example of noninert pairing states will be explicitly discussed. Furthermore, we will consider specific pairing states with these symmetry groups up to total angular momentum J = 4, the highest total angular momentum up to p-wave order, see Sec. II B. Generalization to higher angular momentum channels is straightforward and will be mentioned where appropriate.
where we have indicated the angular momentum J as |∆ O J . (Here, we will not be concerned with the normalization of these states.) The two generators of the isotropy groups are given by {e iπ C 4z , e iπ C 2,z+x } and {C 4z , C 2,z+x }, respectively, where C 2,z+x is a twofold rotation about the (101) axis. We further observe that both states are time-reversal invariant pairing states, and are thus spin-degenerate. The state |∆ O 4 is interesting, since it is invariant under all rotation symmetries of the cube. This implies in particular that odd-parity pairing states |∆ O 4 do not possess any mirror symmetry. Based on the discussion of mirror symmetries presented in the introductory part of this section, and in light of similar considerations in Sec. IV D, we conclude that odd-parity |∆ O 4 states realize fully gapped topological superconductors in class  Table I. DIII. Intriguingly, whereas the nematic topological superconductors described in Sec. IV D have a nonzero quadrupole moment [in the sense of Eq. (17)], the symmetry of the |∆ O 4 states does not allow a quadrupole moment. In fact, the highest nonzero multipole moment (i.e., subsidiary order) is a hexadecapole moment. The even-parity |∆ O 4 states clearly do have mirror symmetry; the planes perpendicular to any of the twofold axes are mirror planes. In the even-parity case this does not mandate degeneracies and this generically leads to a full pairing gap for even-parity |∆ O 4 states. These gapped even-parity superconductors are topologically trivial.
We turn to the J = 3 states |∆ O 3 . Due to the phase factor associated with the twofold rotations the odd-parity |∆ O 3 states are invariant under mirror reflection in planes perpendicular to the (110) axis (and equivalent axes). In contrast, the even-parity |∆ O 3 states are odd under mirror reflection in planes perpendicular to the (110) axis (and equivalent axes). For the even-parity states this implies line nodes on the Fermi surface. Instead, for the odd-parity pairing states the mirror symmetries, in combination with time-reversal invariance, leads to point nodes on the Fermi surface, located along the (001) as well as the (111) directions (and all equivalent directions).
To establish the dispersion of the low-energy gapless quasiparticles at the nodes, we proceed in the same way as in Sec. IV D. We first treat the fourfold axis along (001), i.e., the z-direction. The pseudospin-triplet pairing ∆ q matrix was defined in Eqs. (63) and (64), and takes the form Under fourfold rotation one has C 4z : s ± → e ±3iπ/2 s ± ; s z is invariant. Since the pairing must be odd under fourfold rotation one finds where A ± and B ± are expansion coefficients. One of the mirror planes is perpendicular to (110) and, invoking the same arguments which led to Eq.
which, in combination with the normal state contribution given in Eq. (56), gives rise to a Dirac Hamiltonian for the gapless low-energy quasiparticles, with linear dispersion to lowest order [c.f. Eq. (66)]. We note that the analysis is similar for the case of a pseudospin-1 2 conduction band Fermi surface, giving rise to linear dispersion. Now, consider the threefold axis along the (111) direction. It is convenient to apply a global rotation to the pairing state |∆ O 3 such that the threefold axis is oriented along the z-direction. Alternatively, one may view this as choosing local coordinates (q x , q y ) perpendicular to the rotation axis, see Fig. 1. We furthermore imagine that the state has been rotated such that the mirror symmetry is given by M y : y → −y. The group of symmetries which leave the intersection of the threefold axis and the Fermi surface K invariant is given by C 3v , i.e., the threefold rotations and three equivalent mirror reflections. Importantly, it follows from group theory that the symmetry group C 3v does not protect degeneracies for pseudospin- 3 2 states. In the present case this implies that no point nodes exist along the (111) direction for a pseudospin- 3 2 Fermi surface. In the odd-parity pairing state |∆ O 3 the (valence band) ± 3 2 Fermi surface only exhibits point nodes along the (001) axis (and equivalent axes). This is indeed different for pseudospin-1 2 states: a pseudospin- 1 2 Fermi surface exhibits points along (111) direction. Using the threefold rotations and mirror symmetries (in the rotated basis) we obtain the low-energy pairing matrix ∆ q given by where we have defined q ± = q x ±iq y . This defines another set of Dirac points, in addition to the Dirac points along the (001) direction.

Pairing states with tetrahedral T symmetry
Next, we consider pairing states with tetrahedral symmetry. As compared to the octahedral states, these pairing states lack a fourfold rotation axis; two examples are given by The two generators of the respective isotropy groups are given by {e ∓i2π/3 C 3n , C 2z }, where −, + applies to J = 2, 4, andn = (1, 1, 1) T / √ 3 is a unit vector along the threefold axis. Importantly, the tetrahedral states break time-reversal symmetry but do not have a chirality [Eq. (18)], which follows directly from Eqs. (89) and (90). This is true for general tetrahedral pairing states. As a result, tetrahedral pairing is spin-selective.
Given that the tetrahedral pairing states have a threefold axis along the (111) and equivalent directions, we can invoke the arguments of Secs. IV A and IV E to study the low-energy gap structure at Fermi momenta ±K defined along the (111) rotation axis. Consider oddparity pairing first. As is the case for chiral pairing states (see Sec. IV E) we can focus on the gap functions ∆ q± for pseudospin-↑ and -↓ pairing. The pseudospin-triplet operators have angular momentum ±3 and are therefore invariant under threefold rotations. Consequently, the orbital angular momentum of ∆ q± must match the rotation eigenvalue of the pairing state. This implies that both ∆ q± ∝ (q x − iq y ) for |∆ T 2 , and, similarly, ∆ q± ∝ (q x +iq y ) for |∆ T 4 , giving rise to linearly dispersing Majorana fermions along the rotation axes in each pseudospin sector.
In the case of a conduction band Fermi surface, i.e., when the pseudospin-triplet operators carry angular momentum ±1 and transform as e ±i2π/3 under threefold rotation, either ∆ q+ or ∆ q− can acquire a constant nonq-dependent part. Therefore, only a single pseudospin species of Majorana fermions exists on the rotation axis: the Majorana fermions are fully spin-polarized. A realization of such Majorana fermions were theoretically found in the tetrahedral pairing state of 3 P 2 superfluids [69].
Majorana fermions on the threefold rotation axis are a generic property of tetrahedral pairing states, and the angular momentum of the paired electrons determines whether Majorana fermions of a single or both pseudospin species is present.
Finally, we note that even-parity pairing states with tetrahedral symmetry will generically have Z 2 Fermi surfaces enclosing the Fermi momenta along the threefold rotation axes. This follows from the arguments presented in Sec. IV E. In addition, since the even-parity pairing have a mirror plane orthogonal to the twofold rotation, they generically feature line nodal degeneracies as well.

Pairing states with dihedral Dn symmetry
At last, we turn to the class of dihedral pairing states with isotropy groups D n . As a first example, consider the case n = 8 and the pairing state |∆ D8 4 = |4, 4 − |4, −4 of Eq. (31). This is a time-reversal invariant pairing state with isotropy group {e iπ C 8z , e iπ C 2x }. For oddparity pairing states, this implies four symmetry-related mirror planes (e.g., the xz and yz planes are both mirror planes). As a result, point nodes must appear along the eightfold axis, i.e., the z-axis. To obtain the lowenergy gap structure of the point nodes at momenta ±K along rotation axis, we must require that ∆ q is odd under eightfold rotation and respects all mirror symmetries.
Using that s ± transform as C nz : s ± → e ±6iπ/n s ± for pseudospin ± 3 2 states, we find Instead, for pseudospin ± 1 2 fermions transforming as C nz : s ± → e ±2iπ/n s ± we find that the low-energy quasiparticle dispersion takes the form As before, in both cases the significance lies in the second term. The low-energy gap structure describes quasiparticles with linear and cubic Dirac dispersion, respectively. For the ± 1 2 pseudospin states we thus find a another new type of low-energy quasiparticle: Dirac fermions with cubic dispersion.
Next, consider n = 6 with the pairing states are given by [see Eqs. (28) and (32)] The structure of the isotropy group is the same in both cases and given by {e iπ C 6z , e iπ C 2x }. In addition, the two dihedral states are time-reversal invariant. The latter is true for all pairing states with D 6 symmetry. Clearly, referring earlier arguments, even-parity |∆ D6 3 and |∆ D6 4 pairing states are odd under certain mirror reflections and must therefore have line nodes. For the odd-parity states we make the following observations. The structure of the isotropy group generators implies that the odd-parity |∆ D6 states have a set of three vertical mirror planes given by the mirror operation M x : x → −x and its two equivalents related by threefold rotation. In addition, the odd-parity states are invariant under mirror reflection in the xy plane, since the isotropy group contains the element e iπ C 2z . These constraints have different implications for the pseudospin ± 3 2 and ± 1 2 Fermi surfaces. In particular, in the case of a pseudospin- 3 2 Fermi surface, no points nodes are present along the sixfold rotation axis, i.e., the (001) direction. This follows from the requirement that ∆ q must be invariant under the subgroup C 3v and odd under the sixfold rotations, which is not sufficient to force ∆ q=0 to vanish for ± 3 2 doublets. In contrast, a pseudospin-1 2 Fermi surface has symmetry-protected point nodes along sixfold rotation axis. More specifically, we find that the lowenergy gap structure∆ q is given by This shows that odd-parity |∆ D6 states can realize double Dirac points: low-energy gapless quasiparticles with quadratic dispersion in the x and y directions. The presence of mirror symmetry implies that the pseudospin- 3 2 Fermi surface must have point nodes somewhere, even if they are not located along the sixfold rotation axis (where they might be expected). The location of these point nodes can be determined with the help of the mirror plane perpendicular to the sixfold axis. Indeed, the intersection of the three vertical mirror planes and xy mirror plane defines six points on the Fermi surface, located along the (010) and equivalent directions, which remain gapless. The dispersion of the gapless quasiparticles is linear, As a third example of dihedral states, consider states with D 4 symmetry. Two examples are given by The generators of the isotropy group are given by {e iπ C 4z , C 2x }, and all pairing states with D 4 symmetry are time-reversal invariant. The odd-parity pairing states have a mirror symmetry M (110) : (x, y) → (y, x) and its equivalent related by twofold rotation. This implies point nodes along the fourfold axis. The dispersion of the nodal quasiparticles is derived in the same way as before; we obtain the gap structure ∆ q for momenta ±K along the fourfold z-axis as This shows that the gapless low-energy quasiparticles of odd-parity |∆ D4 pairing states have linear dispersion to lowest order and are yet another realization of Dirac superconductors. (Note that pseudospin-± 1 2 pairing also gives rise to linear dispersion.) Finally, as an example of dihedral pairing states which break time-reversal symmetry, consider the states with D 3 symmetry. Here, x ± ≡ √ 1 ± x and x is a parameter which will depend on the details of the free energy. (These are therefore noninert states.) The generators of the isotropy group of both states are given by {C 3z , e iπ C 2x }. Even though these states do not have a chirality, they break time-reversal symmetry and hence define spin-selective pairing states. Let us focus on the odd-parity realizations of these D 3 pairing states. We then notice that along the threefold axis, the gap functions ∆ q± of pseudospin-3 2 triplet pairing can have a constant part, i.e., ∆ q± ∝ 1, since the corresponding pseudospin-triplet pairing operators are invariant under threefold rotation. This implies a full pairing gap along the rotation axis. In contrast, for ∆ q± ∝ q ∓ for pseudospin-1 2 triplet pairing, giving rise to Majorana fermions with linear dispersion on the rotation axis.

G. Application: cubic crystal anisotropy
Following the detailed exposition of pairing states with discrete symmetry, we conclude this section by demonstrating how the gap structure classification may be directly applied to systems with a normal state exhibiting crystal anisotropy. We focus the discussion on the cubic group, since one of the main motivations of this work are the half-Heusler materials. (We recall that the splitting of the isotropic channels in terms of cubic channels is listed in Table VII.) As discussed in Sec. II B, when crystal anisotropy effects reduce the spatial symmetry group of the (spin-orbit coupled) normal state to the crystal point group, pairing channels are labeled by representations of the crystal point group, and are necessarily finite dimensional. (Recall that the cubic representations have dimension one, two, or three.) In manner fully analogous to Secs. II B and III, one may determine the set of stationary pairing states within each channel using symmetry arguments. The isotropy groups, which can be taken as a definition of distinct stationary pairing states, are necessarily discrete, since they must be subgroups of the normal state symmetry group. A complete list of cubic stationary states and their isotropy groups has been given by Volovik and Gorkov [45].
For the purpose of deriving symmetry-enforced constraints on the gap structure of stationary pairing states, only the isotropy group is needed. The explicit form of the gap function is not required. This is important, since gap functions can be rather complicated in crystal systems due to the infinitely many symmetry-allowed terms within a representation (i.e., the available symmetry quantum numbers are greatly reduced in crystal point groups). In this regard, as far as the question of manifest (symmetry-enforced) gap structure properties is concerned, the question whether the normal state has full rotational symmetry or discrete crystal symmetry is secondary. What matters is the structure of the symmetry group of the pairing state; if it is discrete, as must be case with a cubic normal state, the theory of Sec. IV F applies.
To make this more specific, consider the cubic normal state (spatial) symmetry group O h = O × P . Pairing states can be distinguished based on the parity eigenvalue, and for definiteness here we restrict to odd-parity pairing states. Then, there are five distinct pairing channels, labeled by the cubic representations A 1 , A 2 , E, T 1 , and T 2 (see Sec. II B; here we suppress the oddparity designation). The A 1 and A 2 pairing channels are single-component channels, and therefore give rise to free energies with one unique stationary point. The symmetry of the pairing states then follows directly from the representations; the isotropy groups are generated by {C 4z , C 2,z+x } and {e iπ C 4z , e iπ C 2,z+x }, respectively. These symmetry groups may be recognized as the groups of the two octahedral states of Sec. IV F 1 (see also Table  VII in appendix F), which implies that the gap structure is identical.
The E pairing channel is two-dimensional and the corresponding order parameter can be written as ∆ = (∆ 3z 2 −r 2 , ∆ x 2 −y 2 ) T in the basis of Eq. (14). The GL functional for the two-component order parameter has three minima; two of them are given by ∆ = (1, 0) T and ∆ = (0, 1) T . Both pairing states have dihedral symmetry group D 4 , but in the former case it is generated by {C 4z , C 2x }, whereas in the latter case D 4 is generated by {e iπ C 4z , C 2x }. As a result, the state ∆ = (0, 1) T has the same isotropy group and gap structure as the |∆ D4 states of IV F 3 (see Table VII). In contrast, as discussed in Sec. IV F, odd-parity pairing states with isotropy groups generated by pure rotations (i.e., no phase factors) are fully gapped due to the absence of constraints deriving from mirror symmetry. This directly applies to ∆ = (1, 0) T .
As a final example, consider the three-component pairing channel T 2 . A superconducting order parameter can be defined as ∆ = (∆ yz , ∆ zx , ∆ xy ) T [again in the basis of Eq. (14)]. In total, four distinct pairing states can arise in the T 2 channel. For the purpose of illustration, here we just consider two: the time-reversal invariant state ∆ = (1, 1, 1) T and chiral state ∆ = (1, i, 0) T . The former has dihedral isotropy group D 3 generated by {C 3,x+y+z , C 2,x−y }, which per Sec. IV F implies a full (topological) pairing gap. The time-reversal odd pairing state is left invariant under the group generated by e iπ/2 C 4z , implying that ∆ = (1, i, 0) T is chiral and has (axial) angular momentum +1 along the threefold axis. The gap structure may then be obtained by applying the arguments of Sec. IV E 2 to discrete n-fold rotations.
To summarize these considerations, even in cases where the normal state has discrete crystal symmetry, the gap structure classification developed in this section can be directly applied, since pairing states with discrete symmetry are naturally included. We have demonstrated this explicitly using the example of the cubic group, but the conclusion holds for any other crystal point group.

V. DISCUSSION AND CONCLUSION
In this work we have presented a comprehensive topological gap structure classification of j = 3 2 pairing states, obtained through a systematic analysis of the constraints enforced by symmetry. Our analysis of multicomponent pairing states demonstrates that in strongly spin-orbit coupled systems with higher total angular momentum pairing, and in particular in systems with high-spin pairing, topological pairing states form a significant subset of the class of possible superconducting ground states. Four broad classes of topological pairing states should be distinguished: fully gapped time-reversal invariant topological superconductors, nodal Dirac superconductors, nodal superconductors hosting Majorana fermions, and superconductors with Z 2 protected Bogoliubov Fermi surfaces.
Within each class, further distinctions can be made. For instance, fully gapped topological superconductors can be either isotropic or nematic. Nematic superconductors spontaneously break rotation symmetry and have an anisotropic pairing gap [36]. The latter provides a useful experimental diagnostic, as it does not require phase sensitive probes. Within the class of superconductors with bulk nodal gapless excitations, pairing states can be distinguished based on the Berry monopole charge of the point nodes, where the monopole charge is directly related to the dispersion of the low-energy quasiparticles. For instance, Dirac or Majorana quasiparticles with linear dispersion are different from quasiparticles with quadratic dispersion, and define distinct superconducting states.
Our work shows that Majorana fermions generically occur in spin-orbit coupled j = 3 2 superconductors with multicomponent odd-parity pairing which spontaneously breaks time-reversal symmetry. As discussed in Sec. IV B, spin-orbit coupled superconductors with broken time-reversal symmetry are generically spin-selective, allowing for an effective pairing-induced pseudospin splitting. This splitting is responsible for the lifting of pseudospin degeneracies of nodal points on rotation axes, thereby giving rise to non-degenerate point nodes. As far as multicomponent even-parity superconductors are concerned, the spin-selectiveness of the pairing implies that gap structures generically feature the Z 2 surface degeneracies [30].
The gap structure classification we establish in this work provides a useful framework for interpreting ongoing and future experiments which target bulk properties of superconductors. In particular, thermodynamic probes such as specific heat, penetration depth, or NMR spin relaxation time measurements are sensitive to the nature of low-energy excitations [70]. The low-temperature behavior of these quantities directly reflects the density of low-energy quasiparticle states. More specifically, whereas fully gapped superconductors exhibit exponentially activated temperature dependence, nodal superconductors exhibit a power-law dependence at temperatures T T c . The power-law exponent is directly related to the low-energy quasiparticle density of states, and therefore allows to distinguish nodes with different codimension. Notably, however, the density of states of point nodes depends on the Berry monopole charge, which can give rise to low-temperature behavior expected for nodes of different codimension. As a notable example, point nodes with quadratic dispersion can masquerade as line nodes. Our classification is therefore directly useful for the purpose of assigning candidate pairing states to experimentally observed behavior. It is also worth pointing out that for odd-parity time-reversal symmetry breaking pairing states which host non-degenerate Majorana fermions, a further experimental signature is NMR spin relaxation time anisotropy [58].
The symmetry properties of pairing states are fundamental to our classification of gap structures. These symmetries properties are uniquely encoded in the subsidiary order parameters associated with the superconducting state, which take the form magnetic multipole moments (e.g., dipole moment or chirality; quadrupole moment). Therefore, information on the nature of the superconducting state becomes accessible by probing the structure of the multipole moments. Time-reversal symmetry and rotation symmetry breaking, for instance, can be determined by polar Kerr effect measurements and thermal conductivity or specific heat measurements as function of magnetic field direction, respectively.
The defining physical manifestation of bulk topology are the gapless excitations on the boundary of the material. Topological superconductors with a full pairing gap host two-dimensional gapless Majorana fermions on their surfaces. The existence of these surface Majorana fermions does not depend on surface termination. However, in the case of nematic superconductors, the precise form of the surface quasiparticle dispersion is expected to be anisotropic and sensitive to surface termination with respect to the nematic axis. Bulk nodal superconductors are characterized by gapless Majorana arc surface states, which connect the projections onto the surface Brillouin zone of bulk nodes with opposite monopole charge [57]. As a result, their structure is inherently surface termination dependent. In Dirac superconductors, which are time-reversal invariant and possess a mirror symmetry, these Majorana arcs must come in pairs: Majorana-Kramers pairs [59,63].
The gapless surface excitations can be probed using tunneling microscopy experiments, which couple to the surface density of states. An interesting direction for future work is to study the surface tunneling spectra for different pairing states and different surface terminations.
We conclude this paper by pointing out two important implications of our work. First, since the formalism of our gap structure classification includes pairing states with discrete symmetry, it encompasses the pairing ground states which can arise when crystal anisotropy effects are accounted for. As a result, insofar as the question of quasiparticle gap structure is concerned-the primary interest of this work-the application of our classification is not limited to superconductors with full rotational symmetry. Furthermore, despite our focus on j = 3 2 pairing in the Luttinger model, our topological gap structure classification is directly relevant to higher angular momentum pairing in a more general setting, in particular other systems with strong spin orbit coupling. Symmetry arguments are the work horse of our approach and we therefore expect our analysis of multicomponent topological pairing states to find broad application.
-Note added. After finalization of this manuscript we became aware of a preprint which also considers topological superconducting states in the Luttinger models [71]. momentum, can be obtained as follows. First, notice that the matrices S 1M are proportional to the spin matrices S ± and S z of Eq. (A1). Specifically, one has The higher order multipole matrices S SM , where S = 2, 3, can be obtained applying the recursive formula to the highest weight matrix with M = S. For each S, the highest weight matrix is obtained by setting S SS ∝ S S 11 and requiring that the normalization of the matrices S SM is such that the sum rules The matrices S SM encode the spin multipole structure of the Cooper pair. Since the total spin of the Cooper pair can be S = 1, 2, 3 (apart from S = 0), Cooper pairs can have spin dipole, quadrupole and octupole moments. To highlight the interpretation of spin multipole moments, we take S = 2 as an and construct the multipole components contained in the set S 2M explicitly. The spin quadruple matrices are defined by a rank-2 symmetric traceless tensor Q ab , where a, b ∈ {x, y, z}, given by Symmetric and traceless tensors such as Q ab have five independent components, which precisely matches the number of S = 2 matrices S 2M . The explicit linear correspondence between the five components of Q ab and S 2M is presented in Table VI. Note that Q † ab = Q ba and therefore the quadruple components are real.  is always integer (and not half-odd integer). The matrices I KN have dimensions (2J + 1) × (2J + 1), and can be constructed in the same way as S SM S . The dipole matrices Recall that the highest multipole possible is K = 2J, implying that for a total angular momentum J superconductor 2J distinct subsidiary order parameters can be defined.
(conduction band). Here, µ labels the pseudospin degree of freedom of the two bands, ± 3 2 and ± 1 2 , denoted as µ =↑, ↓. We require that the basis for this pseudopsin is chosen such that |k, µ transform under Θ and P as an a usual spin. This implies P |k, µ = | − k, µ (D1) Θ|k, µ = µν | − k, ν . (D2) Taking the valence band as an example, the matrix which relates the operators f k and c k is defined as V k and can be explicitly represented as where v 1,2 are the vectors of the |k, µ states in the basis of c † kα |0 . Note that this makes V k a 4 × 2 matrix. The relation between f k and c k then reads as Here,P v is the projection operator onto the valence band states, i.e., it projects out operators of the conduction band. Now, the symmetry requirements of Eqs. (D1) and (D2) can be formulated in terms of the eigenvector matrix U k . Using that c k transforms under time-reversal as Θc kα Θ −1 = T αβ c kβ , we find Eq. (D2) implies where ≡ iσ y . Note that * = and T = − . The requirement of inversion symmetry is simply V k = V −k , which is trivially satisfied since the Hamiltonian is even under inversion. Similarly, the matrix of conduction band eigenvectors is defined as W k , and we have The quasiparticle operators c k and c † k can then be expressed in terms of f k and d k as A set of basis vectors v 1,2 and w 1,2 may be found by diagonalizing (k · S) 2 and choosing the eigenvectors such that the requirements of Eqs. Observe that even though the L = 1 spherical harmonics are odd under inversion, these states satisfy v 1,2 (k) = v 1,2 (−k), and therefore the matrix V k of Eq. (D3) constructed from these states trivially satisfies V k = V −k . In the same way, for w 1,2 one has where N = 2Y 2 10 + |Y 11 | 2 . To consider the action of spatial symmetries on the pseudospin operators f k and d k we denote an element of O(3) as R. The spin j = 3 2 quasiparticle operators c † k transform under R aŝ where U R is the j = 3 2 matrix representation of R. We define the matrix representation of R on the pseudospin degree of freedom f † kµ as O R (k), which in general will depend on k. Then, we find that O R (k) is related to U R as from which we obtain O R (k) as Clearly, a similar relation holds for W k .
Appendix E: Low-energy quasiparticle gap structure: Explicit projection In this appendix we derive general expressions for the projected pairing. Specifically, given a pairing potential ∆ k in Eq. (36) we project onto the low-energy Fermi surface degrees of freedom at special momenta ±K. We start by decomposing the quasiparticle operators c k at ±K in terms of valence band and conduction band operators f k and d k ; we find where V ≡ V K = V −K and W ≡ W K = W −K are the matrices of eigenvectors, see Eq. (41). With the help of these relations we expand the pairing Hamiltonian in the vicinity of ±K as where Ψ q and Φ q were defined in Eqs. (51) and (52). The Hamiltonian components H vv q and H vc q (v and c label the valence and conductions bands, respectively) are given by the matrix expressions where we have defined ∆ ±q ≡ ∆ K±q as in the main text, and The Hamiltonian blocks H cc q and H cv q are simply obtained from Eqs. (E3) and (E4), respectively, by substituting V ↔ W . In these expressions ± applies to even-parity (+) and odd-parity (−) pairing states.
Since we are assuming a valence band Fermi surface, the conduction band defines a high-energy manifold. To project the pairing onto the valence band subspace close to ±K we apply perturbation theory. An effective Hamiltonian H vv eff (q) for the valence band subspace is given by an expression similar to Eq. (48) as and we can expand (H cc q ) −1 in powers of (ε c q=0 ) −1 . Here, ε c q=0 = ε c K is the conduction band energy at the Fermi momentum K. The structure of Eq. (E3) shows that H cc q is the sum of the normal state part and the pairing part, and can be expressed as where we have neglected the q-dependence of the normal state contribution. The matrix X q describes the pairing part and ∆ is the overall amplitude of the superconducting order parameter, which we may take to be real. With this we may expand (H cc q ) −1 as (E7) Note that the expansion parameter ∆/ε c K is typically small. Substituting this expansion into Eq. (E5) we obtain Here, P c K ≡ W W † is the matrix projector onto the conduction band states at K (and hence −K). In this expression, the first term can be recognized as a particlehole term and is responsible for the effective Zeeman-type pseudospin splitting, see Eq. (59). To obtain δ one sets q = 0 in the this first term. The first term may also contain a renormalization of the single-particle energies.
The second term describes a contribution to the pairing of valence band due to coupling to the conduction band. As a result, it is smaller by one order of ∆/ε c K and not expected to be of significance.
Appendix F: Connection to cubic symmetry TABLE VII. Splitting of pairing channels in cubic systems. In the presence of cubic crystal anisotropy, the isotropic pairing channels labeled by angular momentum J are split into cubic pairing channels labeled by cubic representations, as explained in Sec. II B. This Table lists the splitting of the isotropic pairing components given in Eq. (7) into pairing functions transforming as partners of the cubic representations; this is shown in the third column, using the notation of Eq. (13). We list the splitting of angular momentum channels J = 1, 2, 3, 4. Note that some cubic representations appear multiple times; the corresponding pairing functions are "degenerate" in cubic symmetry. Note further that the parity g, u depends on the quantum numbers (L, S). A number of pairing states which can arise in cubic systems (and thus necessarily have discrete symmetry, see IV G) are discussed in Sec. IV F; these are listed in the fourth column. Column five refers to the specific equation.