Three-dimensional Majorana fermions in chiral superconductors.

Using a systematic symmetry and topology analysis, we establish that three-dimensional chiral superconductors with strong spin-orbit coupling and odd-parity pairing generically host low-energy nodal quasiparticles that are spin-nondegenerate and realize Majorana fermions in three dimensions. By examining all types of chiral Cooper pairs with total angular momentum J formed by Bloch electrons with angular momentum j in crystals, we obtain a comprehensive classification of gapless Majorana quasiparticles in terms of energy-momentum relation and location on the Fermi surface. We show that the existence of bulk Majorana fermions in the vicinity of spin-selective point nodes is rooted in the nonunitary nature of chiral pairing in spin-orbit-coupled superconductors. We address experimental signatures of Majorana fermions and find that the nuclear magnetic resonance spin relaxation rate is significantly suppressed for nuclear spins polarized along the nodal direction as a consequence of the spin-selective Majorana nature of nodal quasiparticles. Furthermore, Majorana nodes in the bulk have nontrivial topology and imply the presence of Majorana bound states on the surface, which form arcs in momentum space. We conclude by proposing the heavy fermion superconductor PrOs4Sb12 and related materials as promising candidates for nonunitary chiral superconductors hosting three-dimensional Majorana fermions.

INTRODUCTION

Chiral superconductors exhibit Cooper pairing with finite angular momentum, thus spontaneously breaking time-reversal symmetry (). Two-dimensional (2D) chiral superconductors have been extensively studied in the context of Sr2RuO4 (). They are generally expected to have a full superconducting gap and support topologically protected quasiparticles at the edge and in the vortex core. In contrast, 3D chiral superconductors generally have nodes. A well-known example is superfluid 3He in the p + ip paired A phase, which has two point nodes on the Fermi surface along the p axis (). Quasiparticles near these nodes are spin-degenerate and correspond to Weyl fermions (, ). When spin-orbit coupling is present, chiral superconductors with odd-parity (for example, p-wave) pairing may have nonunitary gap structures and spin-selective point nodes (). In this case, despite the fact that the Fermi surface is spin-degenerate, only states of one spin polarization at the nodal points are gapless in the superconducting state, whereas states of the opposite spin polarization are gapped. Consequently, low-energy nodal quasiparticles arise from pairing within states of the same spin. These quasiparticles are identical to their antiparticles and thus are the solid-state realization of 3D Majorana fermions.

Majorana fermions in condensed matter have recently attracted much attention (, ). So far, most of the studies have focused on localized Majorana fermion zero modes in quantum devices. In contrast, 3D Majorana fermions that are naturally occurring as itinerant quasiparticles in bulk chiral superconductors are not well studied. In particular, it has been unclear what distinctive properties these Majorana quasiparticles have and in what materials they are likely to be found.

Here, we develop a systematic approach to classifying different types of Majorana quasiparticles around spin-selective point nodes in chiral superconductors. We present the criterion for these Majorana nodes on a high-symmetry axis based on the symmetry of the superconducting order parameter and the band symmetry in the normal state. We further infer the presence of Majorana nodes away from the high-symmetry axis from the topology of gap structures in the momentum space. We show that the Majorana nature of nodal quasiparticles gives rise to a strongly anisotropic spin relaxation rate depending on the spin direction, which can be directly measured in the nuclear magnetic resonance (NMR) experiment. Similar to Weyl fermions in topological semimetals, the presence of Majorana quasiparticles in chiral superconductors leads to a nodal topological superconductor phase, which exhibits Majorana fermion surface states. As we demonstrate explicitly, zero-energy Majorana surface states form arcs in the surface Brillouin zone, which end at the bulk Majorana nodes. Finally, we propose the heavy fermion superconductor PrOs4Sb12 as a promising candidate for chiral superconductor hosting Majorana quasiparticles.

RESULTS

Symmetry analysis of quasiparticle gap structures

We start with a general symmetry-based analysis of superconducting gap nodes in chiral superconductors with strong spin-orbit coupling and inversion symmetry. We assume that time-reversal symmetry is present in the normal state and is spontaneously broken in the superconducting state due to chiral pairing. We assume that chiral Cooper pairs carry a nonzero total angular momentum (including both orbital and spin) J along a crystal axis of n-fold rotation C, which acts jointly on the electron’s coordinate and spin. Here, n can only be 2, 3, 4, or 6 for discrete rotation symmetry of crystals, and J is only defined mod n. Moreover, because ±J corresponds to time-reversed chiral states, it suffices to consider positive integers J = 1,..., n/2 for n = 2,4,6 and J = 1 for n = 3.

Here, we address the gap structure associated with the points on the Fermi surface along the n-fold axis (hereafter denoted as z), whose momenta are given by ±K = ±kFẑ, where kF is the Fermi momentum. Our approach to deriving the gap structures relies on both symmetry and topological arguments. On the basis of a systematic symmetry analysis, we show that the form of the gap structure at ±K (that is, the C-invariant Fermi surface momenta) is entirely determined by the total angular momentum J of the Cooper pair and the angular momentum of energy bands at ±K in the normal state. Following the symmetry analysis, we invoke a topological constraint on the nodal structure of the quasiparticle spectrum to deduce the full low-energy gap structure, both at and away from ±K. Using these two complementary methods, we will demonstrate the existence of two types of point nodes, located on and off the C axis, respectively.

In the presence of both time-reversal (Θ) and inversion (P) symmetries, spin-orbit–coupled energy bands remain twofold-degenerate at each momentum, and we label the degenerate bands by a pseudospin index, α = ↑, ↓. For simplicity, we will refer to α as spin. The presence of C, Θ, and P symmetries guarantees that one can choose a basis for Bloch states at ±K such that (i) the state with α = ↑ (↓) has an angular momentum j (−j), that iswhere j is a positive half-integer; (ii) PcP−1 = c−; and (iii) ΘcΘ−1 = εαβc−, where εαβ is the Levi-Civita symbol.

Having specified the angular momentum J of the chiral Cooper pair and the angular momentum ±j of Bloch electrons, we are ready to deduce the gap structure near ±K by symmetry analysis. Only pseudospin-triplet pairings, which have odd-parity symmetry, may generate spin-dependent superconducting gaps necessary for 3D Majorana fermions. There are three triplet pairing operators between states near ±K, denoted bywhere Γ1,2,3 at q = 0 carry angular momenta 2j, −2j, and 0, respectively. In general, the pairing potential near ±K is a mixture of these three pairing operators, with corresponding form factors

Because we are interested in the gap structure near q = 0, it suffices to expand Δ(q) to the leading order in qwhere we defined q± = q ± iq, and (q, q) is the momentum tangential to the Fermi surface at ±K. The exponents a, b are integers greater than or equal to zero. When a ≠ b, the smaller of the two determines the leading-order behavior of the gap function, whereas the other can be neglected. When a = b, both terms are equally important and should be kept together.

The form of Δ(q) is constrained by the requirement that the pairing term H carries the angular momentum J. This completely determines the exponents a, b [that is, the analytic form of Δ(q) at small q], allowing us to deduce the gap structures in the vicinity of ±K.

First, consider the case J = 0 mod n, that is, when the superconducting order parameter has effectively zero angular momentum with respect to the C rotation axis. In this case, the triplet pairing component with zero angular momentum Γ3 is allowed at ±K, that is, Δ3(q) is finite at q = 0, creating a full pairing gap without any low-energy quasiparticles.

Second, consider nonzero (mod n) J. If J ≠ 2j mod n, none of the three triplet pairing terms can be finite at ±K, that is, Δ → 0 as q → 0 for all i = 1,2,3. This implies that both spin ↑ and ↓ electrons are gapless at ±K, resulting in spin-degenerate nodes at ±K and non-Majorana nodal quasiparticles. The low-energy Hamiltonian for these gapless quasiparticles can be determined from Eqs. 3 and 4.

Majorana nodes on rotation axis

The type of chiral pairing that gives rise to the Majorana nodal quasiparticles—the focus of this work—corresponds to J = 2j mod n. This implies odd J, and we list all these cases in Table 1. Except for two cases and (to be addressed separately later), we have 2j ≠ −2j mod n. Under this condition, the spin ↑ states that carry angular momentum j are allowed to (and generally will) pair up and form Cooper pairs that carry total angular momentum 2j, whereas the spin ↓ states remain gapless at ±K because of the angular momentum mismatch. The resulting nodal quasiparticles are therefore spin-nondegenerate Majorana fermions.

Table 1

Classification of pairing potentials.

Summary of the classification of pairing potentials Δ≡Δ2, of the spin ↓ states to lowest order in (q+, q−), with q± = q ± iq. The potentials are classified for a given combination of (n, j), where n describes an n-fold rotation axis and j is the spin angular momentum. The chiral superconductor has total angular momentum 2j, and the effective orbital angular momentum of Δ is given by l.

Cn	j	J = 2j	l (mod n)	Pairing Δq
n = 2	j=12	J = 1	l = 0	∝ 1
n = 3	j=12	J = 1	l = −1	∝ q−
n = 4	j=12	J = 1	l = −2,2	∝q−2,q+2
	j=32	J = 3	l = −2,2	∝q−2,q+2
n = 6	j=12	J = 1	l = 2	∝q+2
	j=32	J = 3	l = 0	∝ 1
	j=52	J = 5	l = −2	∝q−2

The low-energy Hamiltonian for these quasiparticles is given bywhere we have defined c ≡ c± and Δ ≡ Δ2,. In addition, ξ ≡ ε − μ, where ε is the single-particle energy-momentum relation and μ is the chemical potential. For small q, we have ξ = vFq, where vF = kF/m is the Fermi velocity in the ẑ direction.

It is instructive to write H in Nambu space by introducing the four-component fermion operator so that H can be expressed aswith the 4 × 4 matrix H(q) taking the general form

The four-component quantum field Ψ satisfies the same reality condition as Majorana fermions in high-energy physics, which reads as in momentum space or, equivalently, in real space, where the Pauli matrix τ acts on Nambu space and Ψ is the transpose of Ψ. This reality condition demonstrates that the low-energy quasiparticles can be regarded as Majorana fermions in three dimensions.

At small q, the pairing term Δ in Eq. 8 can be expanded in powers of q+ or q−. The exponent is determined by the mismatch between the angular momentum of the Cooper pair J = 2j and that of the spin ↓ pairing operator Γ2 at q = 0, which is equal to − 2j. Hence, one finds that

The smallest allowed integer |l| gives the form of Δ to the leading order. For any given (n,j) and with J = 2j fixed, the corresponding l is listed in Table 1 (additional details can be found in the Supplementary Materials).

From Table 1, we find three types of pairing terms Δ with different l’s, which give rise to two types of Majorana fermions with different energy-momentum relations. First, for , one has l = 1 mod n; hence, |Δ| ∝ q⊥, where we defined . This implies that the quasiparticles near the nodes ±K disperse linearly with q in all directions, as governed by the following effective Hamiltonian to first order in qwhere σ = ±1 denotes the two nodes ±K, and vΔ is defined via |Δ| = vΔq⊥ + O(q2). Except for the velocity anisotropy, H(q) is identical to the relativistic Hamiltonian for Majorana fermions in particle physics.

Second, we find several cases for which l = ±2 mod n. According to Eq. 9, this implies that the gapless quasiparticles disperse quadratically in q, q and linearly in q (see Table 1), as governed by the following effective Hamiltonian H(q) to second order in qwhere mΔ is an effective mass defined by .

In the case of fourfold rotational symmetry, that is, n = 4, both and terms, with angular momenta of l = 2 and −2, respectively, are allowed in Δ. As a result, the Hamiltonian H(q) takes a more involved form, which is discussed in the Supplementary Materials.

The above cases of chiral pairing with |l| = 1 and 2 both give rise to gapless Majorana quasiparticles at ±K. According to Table 1, there are two remaining cases that both have l = 0 mod n: and . The property l = 0 mod n implies that spin ↓ states at ±K are allowed to pair and form a Cooper pair carrying the same angular momentum 2j = −2j mod n as the spin ↑ Cooper pair . As a result, both Cooper pairs coexist in the superconducting state and generate a full gap at ±K.

Spin-orbit coupling and nonunitary pairing.

It is clear from our derivation of the Majorana nodal quasiparticles that these can only be present in chiral superconductors with nonunitary gap structures, that is, with a spin-nondegenerate quasiparticle spectrum, such that spin ↑ and ↓ states have different gaps (). So far, nonunitary superconductors have received much less attention than their unitary counterparts. Although nonunitary pairing states have been discussed in relation to UPt3 (–), to Sr2RuO4 (, , ), and recently to LaNiGa2 (, ), the only established example of nonunitary pairing is superfluid 3He in high magnetic fields (), known as the A1 phase. However, from a symmetry point of view, nonunitary pairing is generic and more natural (in a theoretical sense) in chiral superconductors with strong spin-orbit coupling. This is a consequence of the lack of spin-rotation symmetry, replaced by the symmetry of combined spin and momentum rotation under the crystal point group. In these cases, there are typically more than one pairing components with different spin S or orbital angular momentum L but the same total angular momentum J = L + S. As a result, the full bulk gap function Δ of the chiral superconductor, defined with , is generally a mixture of these pairing components, which all belong to the same irreducible representation of the point group. Specifically, Δ can be written aswhere are pairing components (that is, crystal harmonics) with total angular momentum J but different L and S, and λ are dimensionless coefficients describing the admixture of these different components. For each pairing channel J, the set of allowed pairing components depends on both the point group symmetry of the crystal and the spin angular momentum j. In Table 2, we present a full list of gap function components for trigonal (C3), tetragonal (C4), and hexagonal (C6) superconductors and for general spin angular momentum j. Table 2 thus generalizes standard gap function classifications for Bloch electrons () and applies to energy bands of, for instance, electrons, such as reported in the half-Heusler superconductors YPtBi and LuPtBi (, ). In addition, the heavy fermion superconductor UPt3 has recently been proposed to have bands ().

Table 2

Complete set of gap functions for chiral spin-orbit–coupled superconductors.

List of allowed gap function components of Eq. 12 for the chiral pairing channels J = 1,2,3, (pseudo)spin angular momentum , and crystal rotation symmetries C with n = 3,4,6. For each combination (J, j), a complete set of components is given; any other allowed gap function component is generated by multiplying with fully point group symmetry invariant functions (). Because angular momenta are only defined mod n, some entries in the table are equivalent, for example, under C3 symmetry, where is the time-reversed partner of . Recall that s± = s ±is and s are Pauli matrices acting on the Bloch electron (pseudo)spin.

(J, j)	Trigonal (C3)	Tetragonal (C4)	Hexagonal (C6)
(1,12)	k_+s_z,k_zs_+,k_−s_−,k_zk+2s_−k_zk−2s_z,(k+3−k−3)s_+	k_+s_z,k_zs_+,k_z(k+4−k−4)s_+k_zk+2s_−,k_zk−2s_−,k−3s_z	k_+s_z,k_zs_+,k_zk+2s_−k−5s_z,k_zk−4s_−,k_zk−6s_+
(1,32)	k+sz, k+s+, k+s−k_zk−2s_z,k_zk−2s_+,k_zk−2s_−	k_+s_z,k_zs_−,k−3s_zk_zk−2s_+,k_zk+2s_+,k_zk+4s_−	k_+s_z,k_zk−2s_+,k_zk−2s_−k−5s_z,k_zk+4s_+,k_zk+4s_−
(1,52)	k+sz, k−s+, kzs−k+3s_−,k_zk−2s_z,k_zk+2s_+	k_+s_z,k−3s_z,k_zk+4s_+k_zs_+,k_zk+2s_−,k_zk−2s_−	k_+s_z,k−5s_z,k_zs_−k_zk+2s_+,k_zk−4s_+,k_zk+6s_−
(2,12)	≅(−1,12)	k_+s_+,k_−s_−,k_zk+2s_zk+3s_−,k−3s_+,k_zk−2s_z	k_+s_+,k_zk+2s_z,k_zk−4s_zk+3s_−,k−3s_−,k−5s_+
(2,32)	≅(−1,32)	k_+s_−,k_zk+2s_z,k+3s_+k_−s_+,k_zk−2s_z,k−3s_−	k_zk+2s_z,k_−s_+,k_−s_−k_zk−4s_z,k+5s_+,k+5s_−
(2,52)	≅(−1,52)	k_zk+2s_z,k_zk−2s_z,k_+s_+k−3s_+,k_−s_−,k+3s_−	k_zk+2s_z,k_zk−4s_z,k_+s_−k−5s_−,k−3s_+,k+3s_+
(3,12)	≅(0,12)	≅(−1,12)	k+3s_z,k_zk+2s_+,k_zk+4s_−k−3s_z,k_zk−2s_−,k_zk−4s_+
(3,32)	≅(0,32)	≅(−1,32)	k_zs_+,k_zs_−,k_zk+6s_+k+3s_z,k−3s_z,k_zk−6s_−
(3,52)	≅(0,52)	≅(−1,52)	k+3s_z,k−3s_z,k_zk−2s_+k_zk+4s_+,k_zk+2s_−,k_zk−4s_−

As an example of Δ in Eq. 12, consider the following gap function of a J = 1 superconductor of electrons, consisting of two pairing components with (L, S) = (1, 0) and (L, S) = (0, 1), respectivelywhere we defined k± = k ± ik and s± = s ± is. It is straightforward to verify that the pairing is nonunitary as a result of the second (L, S) = (0, 1) term. The first term of Eq. 13 corresponds to the gap function of the A phase of superfluid 3He (). The spin-degenerate quasiparticle spectrum of the 3He-A phase hosts Weyl fermions at low energies (, ), which can be viewed as a complex quantum field made up of two degenerate Majorana fields. The admixture of the (L, S) = (0, 1) component, which is enabled by spin-orbit coupling, gaps out the spin ↑ states at ±K and gives rise to gapless spin ↓ excitations governed by Hamiltonian (Eq. 11) (Supplementary Materials).

This example illustrates a general feature of spin-orbit–coupled chiral superconductors: The lack of spin-rotation symmetry naturally leads to nonunitary pairing, which serves as a spin-selective gapping mechanism and creates spin-nondegenerate nodal excitations, obeying the Majorana reality condition.

As we pointed out earlier, the Majorana condition makes the quantum field Ψ a four-component real field. One may be tempted to rewrite Eq. 5 in terms of a two-component complex quantum field, : , which is invariant under the U(1) transformation, f† → f†e. However, this U(1) symmetry is only present in the presence of translational symmetry and broken by impurity-induced potential scattering between the nodes. To see this, consider the spin-conserving internode scattering termwhere M is the scattering amplitude at the momentum 2kF. In terms of the complex field , H involves terms and thus removes the emergent U(1) symmetry in the clean limit. In terms of the four-component Majorana field Ψ, H is given by

Including Eq. 14 in the Majorana Hamiltonian (Eq. 8), the energy E of the Majorana nodes with linear dispersion, Eq. 10, is given by . Thus, the effect of internode scattering is to generate a mass term without U(1) symmetry, that is, a Majorana mass term. As a result, the fundamental quantum field describing the gapless quasiparticles is a four-component real field, that is, a field obeying the Majorana condition.

We note that spin-nondegenerate point nodes also occur when the Fermi surface in the normal state is already spin-split due to magnetism—as theoretically shown in magnetic topological insulator-superconductor heterostructures () and ferromagnetic p-wave superconductors () and in the mixing of chiral d- and p-waves ()—or spin-orbit coupling in noncentrosymmetric superconductors (–). This is different from our case, where the spin-selective point nodes occur via nonunitary pairing in a spin-degenerate normal state. Also, the present case of nonunitary pairing does not assume any special feature in the band structure and should be distinguished from theoretical surveys of possible pairing states in Weyl and Dirac semimetals (–).

Off-axis point nodes

A key result of our symmetry analysis presented in the first part of this section is the presence of nodal Majorana excitations at the rotationally invariant Fermi surface momentum K on the principal rotation axis for nonzero l mod n (see Eq. 9). Rotational symmetry further dictates the form of the energy-momentum dispersion of these on-axis Majorana quasiparticles. We find that spin-orbit–coupled odd-parity chiral superconductors can have additional point nodes located at generic Fermi surface momenta away from the north and south poles, that is, off-axis Majorana nodes. We will first illustrate the presence of these point nodes using examples and then explain their topological origin.

As a first example, let us consider a chiral superconductor with C3 symmetry and angular momentum J = 1. We now show that its gap structure can exhibit nodes at Fermi surface momenta other than ±K. The full pairing potential Δ of Eq. 12 is given to leading p-wave order by (see the Supplementary Materials for details)where λ are three real admixture coefficients. In addition to the on-axis nodes at k = ±K, the quasiparticle spectrum corresponding to Δ of Eq. 16 exhibits six nodes located at off-axis Fermi surface momenta. Writing the Fermi momenta as , the location of these nodes can be expressed as the relations and (here, we have assumed λλ > 0. If λλ < 0, one should substitute ). There are three nodes on the northern Fermi surface hemisphere, which are related by threefold rotation C3, and each of these nodes has a partner on the southern hemisphere related by inversion P, as shown in Fig. 1C.

Fig. 1

Schematic structure of Majorana point nodes of spin-orbit–coupled chiral superconductors with total angular momentum J = 1 with an n-fold (n = 2, 3, 4, 6) rotation axis along z.

Two types of Majorana nodes are shown: on- and off-axis nodes. Whereas the former are pinned to the rotation axis (that is, ±K), the latter appear at generic Fermi surface momenta. (A) C6-symmetric case with double Majorana nodes at ±K. (B) C4-symmetric case. (C and D) C3- and C2-symmetric cases, respectively, including a view from the top (projection on the xy plane). The gap structure of the C3-symmetric superconductor has both on- and off-axis nodes, whereas that of the C2-symmetric superconductor only has off-axis nodes. Nodes with a positive (negative) monopole charge C (see Eq. 18) are indicated by solid black (white) dots, with the monopole charge (that is, C = ±1, ±2) explicitly given. In case of C4 symmetry, the sign of the Majorana node monopole charge at ±K depends on microscopic details (see the Supplementary Materials).

As a second example, consider a J = 1 superconductor with a twofold rotational symmetry C2. The pairing potential Δ is composed of all terms that are odd under C2, and to p-wave order in spherical harmonics is given by

At the Fermi surface momentum K, one has Δ = Δ0(λs+ + λs−), which implies a pairing gap for both c± and c±, in agreement with the result shown in Table 1. Although no nodes exist on the twofold z axis, it is straightforward to verify that for general nonzero admixture coefficients, two pairs of nodes are located on the Fermi surface, each pair related by twofold rotation with the partners of a pair related by the inversion P. For instance, taking λ = 0 (λ = 0), the nodes are located on the intersection of the Fermi surface with the yz (xz) plane, as shown in Fig. 1D.

In both these examples, all point nodes located at generic rotation noninvariant Fermi surface momenta are spin-nondegenerate, and the gapless quasiparticles are therefore Majorana fermions. We call these nodes off-axis Majorana nodes. The presence of the off-axis Majorana nodes in these examples motivates the question of whether these are accidental or whether the off-axis nodes are related to the Majorana nodes on the rotation axis at a deep level.

We now address this question by focusing on the topological nature of the Majorana point nodes. We show that the classification of different types of low-energy Majorana quasiparticles in terms of location on the Fermi surface (that is, on- or off-axis) and energy-momentum dispersion (that is, linear or quadratic) is linked to a topological property of point nodes in momentum space.

In band theory, crossing points of (quasiparticle) energy bands are endowed with an integer topological quantum number given by the Chern number C defined as

Here, F = ∇ × A(k) is the Berry curvature, that is, the field strength of the momentum space gauge Berry connection A(k), which is integrated over a surface Ω, enclosing the point node. Hence, C quantifies the Berry curvature monopole strength of the point node. Because monopoles cannot be removed unless they annihilate with a monopole of equal but opposite strength, point nodes, which carry a nonzero monopole charge, are topologically protected.

It is straightforward to show that the Majorana nodes with linear dispersion described in Eq. 10 have a monopole charge C = ∓1 at ±K, similar to Weyl fermions in topological semimetals (, ), and that the Majorana nodes in Eq. 11, which disperse quadratically tangential to the Fermi surface, have a monopole charge C = ±2, similar to double-Weyl fermions (). Therefore, the former may be called single Majorana nodes and the latter double Majorana nodes. These single and double on-axis Majorana nodes, corresponding to C3 and C4,6 symmetry, respectively, are schematically shown in Fig. 1, with the monopole charge C explicitly indicated.

From the perspective of topology, the presence of the off-axis Majorana nodes in superconductors with C3 or C2 can be understood from monopole charge conservation. For instance, the C3-symmetric superconductor with the gap function 16 can be viewed as a descendent of the C6-symmetric superconductor with the gap function 13, obtained from lowering the symmetry. As a result of lower symmetry, additional gap functions can mix in with the J = 1 channel, and according to our symmetry analysis, this transforms the C = 2 double Majorana node at K of a C6 superconductor into a C = −1 Majorana node at K of a C3 superconductor. Because the total monopole charge must be conserved and the gap structure must preserve the C3 symmetry, three additional point nodes with a monopole charge of C = +1 must exist. This effective “splitting” of a C = 2 Majorana node into four single nodes agrees with the explicit analysis in Eq. 16 and is schematically shown in Fig. 1. An analogous argument explains the existence of two C = +1 off-axis Majorana nodes in the case of the C2-symmetric superconductor, which has a full pairing gap at K. The nodal structure of the C2 superconductor is depicted in Fig. 1. It is therefore natural to think of the off-axis Majorana nodes as originating from the on-axis double Majorana nodes in a J = 1 superconductor, obtained by lowering C6 or full rotational symmetry to trigonal (C3) and orthorhombic (C2) symmetry.

These topological arguments demonstrate that the chiral superconductors discussed here are topological nodal superconductors. Topological nodal superconductors are superconductors with topologically protected gapless quasiparticle excitations in the bulk, in analogy with the protection of Weyl fermions in Weyl semimetals (see also the “Surface Andreev bound states: Majorana arcs” section) (, ). In particular, the gapless Majorana quasiparticles are topologically protected by monopole charge conservation.

Detecting 3D Majorana fermions

Here, we explore ways to detect the 3D Majorana nodal quasiparticles in chiral spin-orbit–coupled superconductors. Because the Majorana fermions arise as a result of a spin-selective gapping mechanism associated with nonunitary pairing, we first address the experimentally observable consequences of nonunitary pairing from a general perspective and then turn to a specific probe sensitive to the spin-polarized low-energy Majorana quasiparticles—the NMR spin relaxation rate.

Signatures of nonunitary chiral superconductors

The nonunitary gap structure of chiral spin-orbit–coupled superconductors gives rise to a number of distinctive experimental signatures. The most prominent characteristic of nonunitary pairing, the nondegenerate quasiparticle excitation spectrum, leads to a spin-dependent density of states: The densities of states N↑,↓(E) for (pseudo)spin ↑,↓ excitations are unequal [that is, N↑(E) ≠ N↓(E)], and in particular, the low-energy branch of the spectrum consists of fully spin-polarized states near the point nodes. As a consequence of the nondegenerate spectrum, the total density of states ∑αNα(E) can exhibit two distinct peaks, rather than a single peak at Δ0, which is characteristic of conventional s-wave superconductors. This can lead to a two-gap–like feature in the specific heat, as is demonstrated more explicitly in simple examples in the Supplementary Materials. Therefore, experimental signatures that are commonly attributed to multiband superconductivity may actually originate from nonunitary pairing in a single band.

A well-known and discriminating property of nodal superconductors is the characteristic temperature dependence of a diverse set of dynamic and thermodynamic quantities, such as the electronic part of the specific heat, the London penetration depth, and the NMR spin relaxation rate () (). The low-energy branch of the quasiparticle spectrum of chiral nonunitary superconductors consists of C = ±1 or C = ±2 point nodes, which implies that the density of states N(E) in the regime E ≪ Δ0 takes the form N(E)/N0 ~ (E/Δ0), where N0 is the normal-state density of states. The exponent n depends on the nature of the point node: It equals n = 2 for C = ±1 and n = 1 for C = ±2. The form of the density of states at the nodes is responsible for the typical power law temperature dependence of the specific heat, the penetration depth, and the spin relaxation, which probe the density of quasiparticle states, at temperatures T ≪ Tc ~ Δ0.

In the next section, we consider the NMR spin relaxation rate in more detail. In the case of chiral superconductors with low-energy Majorana quasiparticles, not only should the spin relaxation rate exhibit power law temperature dependence at low temperatures, but the temperature dependence of is also expected to crucially depend on the direction of the nuclear spin polarization. This follows from the fact that the only quasiparticle excitations available at low energy are spin ↓ states, which is intimately related to the nonunitary nature of the superconducting state. Therefore, we derive below the theory of NMR spin relaxation in nonunitary superconductors, which serves as a powerful tool to identify Majorana nodal quasiparticles.

NMR: Spin relaxation rate.

The measurement of the NMR spin-lattice relaxation rate 1/T1 at low temperature is a well-established experimental technique to probe the gap structure of superconductors. The temperature dependence of 1/T1 at temperatures T ≪ Tc ~ Δ0 can be used as a measure of the density of low-energy quasiparticle states and allows to distinguish fully gapped, point nodal gap, and line nodal gap structures ().

The coupling of quasiparticle states to the nuclear spin originates from the hyperfine interaction between the nuclear spin and itinerant electrons. The hyperfine coupling Hamiltonian Hhf is given bywhere Ŝ are the components of the nuclear spin operator and s are Pauli matrices representing the electron pseudospin (summation of repeated spin indices α, β is implied). Furthermore, γ is the gyromagnetic ratio of the nuclear spin, Ahf is the hyperfine coupling constant, and g(k, k′) is a momentum-dependent tensor describing the coupling of the nuclear spin to the pseudospin of electrons in spin-orbit–coupled materials. The form of this tensor can be complicated and material-specific. However, because only quasiparticles around the nodes contribute to the spin relaxation at low temperature, it suffices to consider g(k, k′) at the nodes, a key simplification that enables us to find universal features of the spin relaxation rate below.

Our aim is to derive 1/T1 for chiral nonunitary superconductors hosting 3D Majorana fermions. Because the Majorana nodes are nondegenerate with definite pseudospin, the low-energy quasiparticles couple anisotropically to the nuclear spin. To demonstrate this, consider the case of gapless Majorana particles pinned to ±K. Projecting the Hamiltonian (Eq. 19) into the space of Bloch states near ±K (P projects onto the low-energy Hilbert space), one findswhere g1 = g(K, K), g2 = g(−K, −K), and g3 = g(K, −K) are the g tensors evaluated at the nodes. Only the z component of the nuclear spin enters the low-energy Hamiltonian due to the rotational symmetry of the crystal around the z axis. As a result, we expect that the nuclear spin relaxation time 1/T1 is highly direction-dependent and, to leading-order approximation, diverges when the nuclear spin is initially polarized along the z direction.

As we show below, the strongly anisotropic spin relaxation rate is a generic consequence of spin-selective nodes, whereas the divergence for nuclear spin polarization along the nodal direction is an artifact of the leading-order Hamiltonian (Eq. 20). We now go beyond this leading-order approximation and expand the full form factor g(k, k′) into crystal spherical harmonics. Keeping the lowest-order s-wave component, we have g(k, k′) ~ δ, reducing the hyperfine coupling to

Using Eq. 21, we proceed to explicitly calculate the NMR relaxation rate 1/T1 for a nonunitary superconductor. For simplicity, we consider a nuclear spin of S = 1/2. The spin relaxation rate 1/T1 is expressed through the transverse spin susceptibility χ−+(k, ω) (that is, transverse to nuclear spin direction) and reads as (, )

When evaluating 1/T1, it is important to distinguish unitary and nonunitary pairing states. In the case of former, it has been shown that the spin relaxation rate does not depend on the direction of polarization of the nuclear spin (reviewed in the Supplementary Materials). The case of nonunitary pairing requires separate and more careful treatment. For concreteness, we first consider the C6-symmetric J = 1 superconductor, which only has double nodes at ±K. Other nodal nonunitary superconductors are discussed toward the end of this section.

The gap structure of the J = 1 superconductor with C6 symmetry is given by Eq. 13. The full BdG Hamiltonian is diagonalized in terms of Bogoliubov quasiparticle operators, which we define as a. Writing the hyperfine interaction Hamiltonian (Eq. 21) in terms of the Bogoliubov quasiparticle operators, the spin susceptibility can be readily evaluated. Because only low-energy excitations contribute to the relaxation rate at low temperatures (that is, ~T ≪ Δ0 min{λ,λ}), we can restrict to Bogoliubov quasiparticle states with small momenta q relative to ±K and only keep the low-energy gapless states a corresponding to energies (see Eq. 11), where . Then, to the leading order in small momentum q, the Hamiltonian (Eq. 21) reads aswith the form factors and G(q, q′) defined asand the momentum-dependent factors . Here, p, p′ = +1 (−1) for the north (south) node.

In Eq. 23, the anomalous terms aa and have been omitted because they do not contribute to 1/T1 due to energy conservation. The effective mass mΔ of the Bloch states c± and energy associated with the Bloch states c± are defined in terms of gap function parameters in the Supplementary Materials. Equation 23, which is the projection of Hamiltonian (Eq. 19) into the space of low-energy Bogoliubov quasiparticles, should be compared to the projection into the space of low-energy Bloch states given in Eq. 20. The former contains a coupling to Ŝ± = Ŝ ± iŜ, originating from the nonzero support of the Bogoliubov quasiparticle states on the Bloch electron states c±. The support is vanishingly small near the nodes as a result of , indicating a suppression of the spin relaxation rate for a nuclear spin initially polarized along the z direction. In the limit that this coupling is strictly absent, as in Eq. 20, the spin initially polarized along the z direction does not relax.

Using the projection of the hyperfine coupling into the space of low-energy Bogoliubov quasiparticles given by Eq. 23, it is straightforward to calculate the NMR relaxation rate given by Eq. 22 in the low-temperature limit T ≪ Δ0 min{λ,λ} (Supplementary Materials). We find 1/T1 to be given bywhere S⊥ and S are the projections of the nuclear spin polarization on the xy plane and z axis, respectively, and the coefficients D⊥, are given by

Equation 25 proves that there is a strong anisotropy of spin relaxation depending on the polarization of the nuclear spin. For a nuclear spin polarized perpendicular to the z axis, the spin relaxation behaves as 1/T1 ~ T3, as expected for C = ±2 point nodes with quadratic dispersion and linear dependence of density of states on energy. Instead, for a nuclear spin polarized along the z axis, the spin relaxation rate is suppressed by a factor of , that is, the zeroth-order term in an expansion in is absent.

The strong relaxation rate anisotropy can be intuitively understood from a simple physical picture. The hyperfine interaction leading to nuclear spin relaxation is given by Eq. 21, and consequently, nuclear spin relaxation occurs simultaneously with an electron spin flip. However, at low energies, only the c± Bloch states are available, implying that a nuclear spin initially polarized along z is vanishingly improbable to relax. This physical picture is captured by Eq. 20. The qualitative difference of this result compared to the case of unitary pairing () can be understood in a similar way. In the case of the latter, the quasiparticle energy spectrum is doubly degenerate, implying that both spin species are present and, consequently, leading to the same temperature dependence of 1/T1 for all nuclear spin polarizations.

A qualitatively similar anisotropic spin relaxation rate has been predicted for 2D Majorana fermions, which live on the surface of a 3D topological superfluid, that is, the 3He-B phase (, ). In this case, where the 2D surface Majorana modes are described by a two-component real quantum field, the anisotropy in the spin relaxation arises because one can only construct the Ising spin operator, which points perpendicular to the surface and takes the form of a Majorana mass term. This is different from our case, that is, with the 3D Majorana fermions, where the spin operator constructed from the low-energy quasiparticles does not correspond to a mass term.

On the basis of the result for the J = 1 superconductor with C6 symmetry and quadratic point nodes at ±K (Eq. 25), we now comment on other chiral nonunitary superconductors. First, consider the J = 1 superconductor, with C4 symmetry shown in Fig. 1 and discussed in more detail in the Supplementary Materials. Because the C4-symmetric superconductor only has on-axis Majorana nodes with C = ±2, Eq. 25 remains valid, and a significant suppression of 1/T1 for the nuclear spin polarized along the z axis is expected.

Next, consider the J = 1 superconductor with C3 symmetry with gap function 16. In this case, the low-energy gap structure consists of eight single Majorana nodes (that is, linear dispersion), two of which are located at ±K, and six at off-axis Fermi surface momenta (see Fig. 1). This nodal structure complicates the explicit derivation of an analytical expression for 1/T1, but the final result, however, can be inferred from the C6 symmetric case. For a weak trigonal anisotropy and intermediate temperatures given by the condition Δ0λ < T < Δ0 min{λ, λ}, the trigonal λ term can be neglected, and one can expect the same behavior as in the hexagonal case with the relaxation rate given by Eq. 25.

However, at the lowest temperatures given by T ≪ Δ0 min{λ, λ, λ}, the linear dispersion of the nodes comes into play: The density of states becomes a quadratic function of energy, leading to 1/T1 ~ T5. Whereas at the on-axis nodes only quasiparticles with definite spin ↓ are available at low energies, the low-energy quasiparticles at the off-axis nodes are mixtures of spin ↑ and ↓. As a result, the temperature dependence of 1/T1 will exhibit the same power law behavior, that is, ~T5 for all nuclear spin polarizations. However, the numerical prefactors will reflect a directional anisotropy, which may be large.

Superconductors with C2 symmetry, due to the four off-axis linear nodes, are expected to show behavior similar to C3 crystals, exhibiting strong anisotropy in the functional temperature dependence at intermediate temperatures and point node power law behavior at the lowest temperatures for all nuclear spin polarizations. Again, numerical prefactors will generically be different.

Surface Andreev bound states: Majorana arcs

Chiral superconductors with point nodal quasiparticle excitations are topological nodal superconductors due to the nonzero Berry monopole charge of the point nodes (), as discussed following Eq. 18. Through the bulk-boundary correspondence, the topological nature of the bulk superconducting state is reflected on a surface boundary separating the superconductor from the vacuum or, equivalently, a gapped s-wave superconductor. The surface Andreev bound states of topological nodal superconductors take the form of arcs in surface momentum space, connecting the projections of the topological bulk nodes onto surface momentum space, as schematically shown in Fig. 2. A canonical example of surface states in nodal topological systems is the surface Fermi arcs of the superfluid 3He-A phase, which originate from and terminate at projections of the bulk Weyl nodes (, ). A similar surface state structure has been explored in the context of UPt3 ().

Fig. 2

Schematic representation of Majorana arc surface Andreev bound states of nodal superconductors.

For a given surface termination, the projections of the bulk Majorana nodes onto the surface momentum space (transparent gray planes) are connected by the surface Majorana arcs (thick blue lines). The surface Majorana arcs must start and terminate at nodes with opposite monopole charge. (A) Arc structure on a side surface of the A phase of 3He. (B) Schematic arc structure of C3-symmetric J = 1 chiral superconductor for a side surface in the y direction (see also Fig. 3). The projection of bulk Majorana nodes (coming from northern Fermi surface hemisphere) on the top surface is also shown (compare Fig. 3B).

Fig. 3

Majorana arc surface states.

Plots of the zero-energy (E = 0) surface Majorana arc states in surface momentum space for chiral J = 1 superconductors with C6 symmetry (A) and C3 symmetry (B to D) and gap functions 13 and 16, respectively. (A), (C), and (D) show the Majorana arc states of a surface boundary in the xz plane (that is, semi-infinite superconductor at y > 0), whereas (B) shows the Majorana arc states of surface boundary in the xy plane (superconductor z > 0). As all panels show, the surface Majorana arcs connect the projections of the bulk Majorana nodes. The dashed circle shows the radius of the Fermi surface projection. In (A), the straight dashed blue line denotes the Fermi arcs of superfluid 3He-A for comparison. The parameters used are given by λΔ0/μ = 0.013, λΔ0/μ = 0.01, and λΔ0/μ = 0.004, 0.009 in (B), (C), and (D), respectively.

Weyl nodal fermions have recently attracted a great deal of attention in semimetallic materials, referred to as Weyl semimetals (, , –), which are semimetals with nondegenerate point nodal touchings of bulk energy bands, associated with nonzero Berry monopole charge, and Fermi arcs on the surface. The Weyl semimetal state has recently been predicted and observed in the TaAs materials class (–) and photonic crystals (–). The superconducting analog of Weyl semimetals was first considered by Meng and Balents () based on a topological insulator–s-wave–superconductor heterostructure model. The surface arcs connect the projections of the nondegenerate bulk nodes in the Bogoliubov quasiparticle spectrum, and because of the redundancy built into the BdG mean-field description, these surface arcs are Majorana arcs.

In general, nodal touchings of energy bands can only occur when at least one of two symmetries, time-reversal Θ or inversion P symmetry, is broken (). The experimentally found Weyl semimetals in TaAs and related materials all break inversion symmetry, and although much effort has been devoted to looking for time-reversal breaking Weyl semimetals (, , , , –), their conclusive observation in materials remains an open challenge.

The chiral superconductors in this work break time-reversal symmetry and have odd-parity pairing, implying that they preserve a Z2 symmetry, given by τP (see the Supplementary Materials). Hence, the only symmetry manifest at any surface is particle-hole symmetry. Because the bulk point nodes are Majorana nodes, the surface states of chiral nonunitary superconductors are Majorana arcs.

In this section, we calculate and study the structure of the Majorana arcs of chiral nonunitary superconductors, with a focus on J = 1 superconductors with C6 and C3 symmetry with gap functions 13 and 16, respectively. The profile of surface Majorana arcs depends on the projections of the bulk Majorana nodes onto the surface and therefore depends on boundary geometry. Here, we will consider two different semi-infinite geometries: (i) a boundary in the xz plane, separating the vacuum (y < 0) and the superconductor (y > 0), and (ii) a boundary in the xy plane (that is, vacuum z < 0 and superconductor z > 0).

Starting with a boundary in the xz plane, the first quantized Hamiltonian is given by HBdG(k, r) = HBdG(k)θ(y), where HBdG(k) is the (first quantized) BCS-BdG Hamiltonian of the superconductor with pairing potential 13. We solve the equation HBdG(−i∇, r)Ψ(r) = EΨ(r) for zero-energy solutions (E = 0), localized at the boundary, that is, wave functions decaying exponentially at y → ∞ and satisfying the Dirichlet boundary condition Ψ| = 0. In this geometry, k and k remain good quantum numbers, and we substitute k → −i∂ in HBdG(k).

The gap function of a C6-symmetric superconductor, given by Eq. 13, gives rise to nondegenerate C = ±2 bulk nodes at ±K. To solve for the zero-energy states of the Majorana arcs, we look for a general solution of the form Ψ(y) ∝ exp(αy), with the condition Re α < 0. The equation HBdG(−i∇)Ψ = 0 then translates into a polynomial in an α of degree eight, whose roots determine the wave function solutions. The roots can be found explicitly and are given bywhere we defined and , and the eigenvalues with Re α1,2 < 0 are implied because only they can be used to satisfy the localization condition Ψ| → 0. These four eigenvalues are then used to construct the solution that satisfies the Dirichlet boundary condition at y = 0, Ψ| = 0. This condition results in the implicit equation for the zero-energy Majorana arc profile in the k − k plane. The implicit equation reads as

The resulting profile of (zero-energy) Majorana arc states is bow tie–shaped in surface momentum space and is shown in Fig. 3A. The zero-energy solutions are nondegenerate, apart from the surface momentum (k, k) = (0, 0), but are related by particle-hole symmetry. This profile should be compared to the fully degenerate surface Majorana arcs of superfluid 3He-A, corresponding to λ = 0 in Eq. 13 and shown with a dashed line in Fig. 3A (see also Fig. 2). As a result of the degeneracy of the two sheets of Majorana arcs, they effectively form a single complex Fermi arc.

In the vicinity of (k, k) = (0,0), the structure of the Majorana arcs of the C6-symmetric chiral superconductor can be obtained within the semiclassical or Andreev approximation [we follow Park et al. ()]. Defining as the semiclassical momentum perpendicular to the boundary and assuming that k⊥ ≫ kF(Δ0/εF), one obtains the Andreev Hamiltonian from the BdG Hamiltonian as HBdG(k, r) → H⊥(k⊥, −i∇) + H∥(k∥), with k∥ ≡ (k, k) and H⊥,∥ given in the Supplementary Materials. In obtaining the surface Majorana arc Hamiltonian, one first solves H⊥(k⊥, − i∇) for zero-energy solutions at k∥ = 0 and projects H∥(k∥) into the subspace of these solutions. At k∥ = 0, we find two solutions from which we construct the second quantized operators and , satisfying the Majorana condition . Projecting H∥(k∥) into the subspace of , we obtainwhere labels the surface Majorana degree of freedom and I2 is a 2 × 2 identity matrix. The profile of zero-energy states is simply obtained as |λk| = |λk|, in agreement with Eq. 27 and Fig. 3.

Next, we consider the case of the C3-symmetric superconductor with gap function (16). As discussed in the “Off-axis point nodes” section, the gap structure consists of eight single Majorana bulk nodes shown in Fig. 1 (also in Fig. 2, including projections onto surface momentum space). An analytical expression for the zero-energy mode profile analogous to Eq. 27 cannot be obtained in this case, and we numerically solve the characteristic polynomial. The result is shown in Fig. 3 (C and D).

Similar to the case of C6 superconductors, the Hamiltonian of the surface Majorana arcs can be constructed in the vicinity of k∥ = 0 within the Andreev approximation. Again, we obtain two surface state Majorana operators . The Hamiltonian of the Majorana arcs reads aswhere the coefficients A and A are given by and , respectively. The zero-energy states of the Majorana arcs in the vicinity of k∥ = 0 are then given by the equation and exist only provided . Only in this range of parameters is the Hamiltonian (Eq. 30) meaningful (otherwise, there is no low-energy excitations near k∥ = 0). This result is in total correspondence with the exact numerical solution shown in Fig. 3 (C and D).

To conclude, we consider a surface boundary in the xy plane. In this geometry, the bulk C = ±2 nodes of the C6 superconductor project to a single surface momentum k = k = 0 such that no well-defined Majorana arcs exist in this case. In contrast, the projections of the six off-axis C = ±1 nodes of the C3 superconductor do not coincide and are connected by Majorana arc states, as shown in Fig. 3B. As is clear from Fig. 3B, threefold rotational symmetry is preserved.

The above analysis shows that spin-orbit–coupled chiral superconductors are topological nodal systems with characteristic surface arc states. Chiral superconductors hosting Majorana fermions in the bulk have spin-nondegenerate Majorana arcs on surface boundaries, unlike complex fermions in Weyl semimetals or spin-degenerate boundary states in superfluid 3He-A.

Candidate materials

The purpose of this section is to connect the general theory of spin-orbit–coupled odd-parity chiral superconductors presented in this study to reported experimental evidence for chiral and nodal pairing in certain known superconductors. Because chiral nonunitary pairing critically relies on spin-orbit coupling, we focus the search for candidate materials hosting Majorana fermions on materials with spin-orbit coupling.

Heavy fermion materials are typically strongly spin-orbit–coupled, and the vast majority of known heavy fermion superconductors are believed to have unconventional pairing symmetry. Of particular interest to this study are two heavy fermion compounds with filled skutterudite structure: PrOs4Sb12 (, ) and PrPt4Ge12 (). For both materials, signatures consistent with point nodes have been observed, although the determination of the pairing state is not yet definitive and the Majorana nature of nodal quasiparticles remains to be tested.

The skutterudite superconductor PrOs4Sb12 has tetrahedral crystal structure with the point group Th. Thermal transport measurements are indicative of a superconducting phase with point nodes (). The presence of point nodes has been further corroborated by the temperature dependence of the specific heat (, ), the penetration depth and NMR spin relaxation rate (), and especially Sb-NQR (), finding spin relaxation rate proportional to T5 at temperatures considerably below Tc. In addition, μSR measurements have been interpreted as supporting time-reversal symmetry breaking in the superconducting state (), and Knight shift measurements are suggestive of triplet pairing (). Very recent Kerr angle measurements provide even more support for time-reversal symmetry breaking in the superconducting state ().

A number of theoretical studies have proposed unconventional pairing symmetries as possible descriptions of the superconducting phases in PrOs4Sb12 (, ). It has been argued that, assuming broken time-reversal symmetry and the existence of point nodes, to best fit experiments, the phenomenological order parameter should be of the three-component Tu symmetry (). The time-reversal symmetry broken phase then corresponds to the chiral combination (that is, phase difference e±) of two components with different amplitude. Within this framework, the resulting chiral superconductor is a nonunitary pairing state with twofold rotational symmetry and nondegenerate point nodes not pinned at a twofold axis. This quasiparticle spectrum can be captured by the pairing gap function given by Eq. 17. As a result, the heavy fermion superconductor PrOs4Sb12 is a promising candidate for realizing the off-axis gapless Majorana fermions.

Another member of the family of filled skutterudites with the same crystal structure as PrOs4Sb12 is the material PrPt4Ge12. Superconductivity has been observed in PrPt4Ge12 (), and penetration depth in combination with specific heat measurements has provided evidence of point-like nodes (, ). Furthermore, a subsequent μSR study has reported a spontaneous magnetization below Tc, which suggests time-reversal symmetry breaking (). These results point toward PrPt4Ge12 as a second candidate to host gapless Majorana fermions. However, an NQR spin relaxation study has provided support for a weakly coupled BCS superconductor with an anisotropic s-wave pairing gap with point nodes (). In addition, it should be noted that some experimental studies have found evidence for two-band superconductivity in PrPt4Ge12 (–), similar to the case of PrOs4Sb12 (–). As discussed when we considered the experimental detection of Majorana fermions, this may still be consistent with nonunitary superconductivity.

DISCUSSION

A central pillar of this paper is the fact that time-reversal breaking pairing in crystals with strong spin-orbit coupling is generically nonunitary, leading to a (spin-)nondegenerate quasiparticle gap structure. In spin-orbit coupling systems, chiral pairing channels are labeled by the total angular momentum J of the Cooper pairs, and the gap function is a general linear superposition of spherical harmonics degenerate in the pairing channel. As shown in the “Symmetry analysis of quasiparticle gap structures” section, the resulting gap function is typically nonunitary.

A consequence of the nonunitary nature of the pairing is the presence of nondegenerate point nodes, which satisfy the same Majorana reality condition as in high-energy particle physics and therefore constitute Majorana fermion quasiparticles in three dimensions. Depending on the (discrete) n-fold rotational symmetry of the crystal and the angular momentum j of the Bloch electrons, these Majorana nodes can be single nodes with linear dispersion or double nodes with quadratic dispersion tangential to the Fermi surface. In addition, symmetry may pin the nodes to the rotation axis, in which case we call them on-axis nodes. When the rotational symmetry of the crystal is low, the on-axis nodes can be split and the off-axis nodes (that is, nodes at generic nonrotation invariant Fermi surface momenta) can appear. The splitting of Majorana nodes and the appearance of off-axis nodes are determined by the conservation of the topological monopole charge associated with point nodes. Hence, it is an analog of the trigonal warping of 2D Dirac fermions in bilayer graphene (). Here, we find both trigonal and orthorhombic warping of Majorana nodes.

We note that our symmetry-based approach is general and complete in the sense of treating the general case of Bloch electrons with angular momentum j, forming Cooper pairs with total angular momentum J in C-symmetric crystals. In particular, our analysis applies to superconductors where the pairing is not between spin- electrons but more generally between spin-j electrons, such as as in half-Heusler compounds (, ) or ().

Experimental manifestations of nonunitary pairing, and consequently of Majorana nodal fermions, can be searched by means of probes sensitive to the difference in spin ↑ and ↓ gap structure. For instance, the density of states clearly reflects the nondegeneracy of the quasiparticle spectrum, giving rise to two-gap features reminiscent of multiband superconductors. We demonstrated that the NMR spin relaxation rate in nonunitary superconductors with spin-selective low-energy quasiparticle excitations shows a marked anisotropic dependence on the polarization of the nuclear spin. We find that for a nuclear spin initially polarized along z, the relaxation rate is significantly suppressed, which serves as a discriminating feature of nonunitary superconductivity.

Chiral superconductors hosting Majorana fermions belong to the class of topological nodal superconductors. Topological nodal superconductors are analogs of topological semimetals called Weyl semimetals. Similar to the latter, topological nodal superconductors have special surface states: Majorana arcs in momentum space connecting projections of bulk nodes. Whereas the Majorana arcs come in (spin-)degenerate pairs (effectively forming Fermi arcs) in unitary superconductors, nonunitary superconductors have nondegenerate Majorana arcs at surface boundaries, as demonstrated in the previous section.

Most superconductors extensively studied in the literature are examples of unitary pairing states. Comparatively, nonunitary superconductors have received less attention. As discussed in the first part of Results, from a conceptual standpoint, relying on symmetry arguments, nonunitary pairing is natural and potentially widespread in spin-orbit–coupled systems. We propose the heavy fermion superconductor PrOs4Sb12 as a promising candidate for this nonunitary pairing and consequently as a realization of Majorana fermions in three dimensions.

MATERIALS AND METHODS

Details on symmetries and BCS-BdG mean-field theory

This section introduces the BCS-BdG mean-field theory and the implementation of symmetries from a more formal perspective. In crystals with strong spin-orbit coupling and in the presence of both time-reversal (Θ) and parity (P) symmetries, all electronic bands remain twofold-degenerate, and the bands are labeled by an effective pseudospin α = ↑, ↓. In the presence of Θ, P, and crystal rotation symmetry C, one can choose a basis such that the electrons c transform under Θ and P asrespectively, and under n-fold rotation aswhere k* = Ck. Furthermore, θ = 2π/n is the angle of rotation, and ±j is the total angular momentum of the Bloch electrons c. Note that as a result, α = ↑,↓ labels the general angular momentum ±j states.

The normal-state Hamiltonian is given by , with H0(k) = ξδαβ. The energy relative to the chemical potential, given by ξ = ε − μ, is a scalar function of momentum composed of terms invariant under the crystal symmetry group. The invariance of the normal-state Hamiltonian under an n-fold rotation C is explicitly expressed as .

In a BCS-BdG mean-field theory formulation, the pairing Hamiltonian HΔ is expressed aswhere the pairing matrix Δ is the momentum-dependent gap function and s are Pauli spin matrices acting on α = ↑,↓. We define the Nambu spinor Φ in the canonical basis with the following

In terms of the Nambu spinor, the mean-field theory BdG Hamiltonian takes the form with

The pairing Hamiltonian transforms as ΘHΔΘ−1 under time reversal, which implies that the gap function transforms as . Chiral superconductors break time-reversal symmetry, and one has . Similarly, odd-parity pairing defined as PHΔP−1 = −HΔ implies Δ− = −Δ for the pairing gap function. As a result of the latter, the BdG Hamiltonian has an effective Z2 symmetry given by τP such that τPHBdG(k)(τP)−1 = HBdG(−k). In addition to these symmetries, the BdG Hamiltonian manifestly obeys a particle-hole symmetry ΞHBdGΞ−1 = −HBdG, which implies for HBdG(k): .

Details on chiral nonunitary superconductors with angular momentum J

This paper studies odd-parity chiral superconductors in which the Cooper pairs have a total angular momentum J = L + S. Because of spin-orbit coupling, only total angular momentum is a good quantum number. The pairing gap function of angular momentum J chiral superconductors is defined through the n-fold rotation C, which we assume to be a rotation about the z axis. Concretely, the pairing gap function satisfies

In a crystal with discrete C rotation symmetry, angular momentum is only a defined mod n, which was manifested in Eq. 36. As a consequence of Eq. 36, total angular momentum J = L + S labels the different pairing channels of the chiral superconductors. Specifically, the gap function takes the general form of Eq. 12, and all functions with combinations (L, S), such that L + S = J, are allowed to mix in with coefficients λ.

The superconducting gap function of Eq. 12 can be explicitly expanded in the spin matrices s as Δ = d(k)⋅s, where d(k) is the momentum-dependent vector. One then finds for

When d* × d = 0, the pairing is said to be unitary, and the quasiparticle spectrum of unitary pairing states is manifestly twofold-(spin-)degenerate. In contrast, when d* × d ≠ 0, the pairing is said to be nonunitary, and the quasiparticle energies are given by . Consequently, nonunitary pairing states are characterized by a nondegenerate quasiparticle energy dispersion, with different gap structures for the spin ↑ and spin ↓ electrons.

As a consequence of spin-orbit coupling and the lack of spin-rotation symmetry, the generic gap function given by Eq. 12 corresponds to nonunitary pairing states. This is easily seen using Table 2.

Details on NMR spin relaxation calculation

This section briefly recapitulates how the NMR relaxation rate is obtained. We calculate the NMR relaxation rate using Fermi’s golden rule, which could be shown to be equivalent to Eq. 22. Defining the nuclear spin coherent state |S〉 for a spin initially polarized along S = (sin θ cos φ, sin θ sin φ, cos θ) as |S〉 = [cos(θ/2), esin(θ/2)] (we considered nuclear spin for simplicity), Fermi’s golden rule for the NMR relaxation takes the formwhere a are Bogoliubov quasiparticles, f is the Fermi-Dirac distribution function, and E1,2(k) are eigenenergies of BdG Hamiltonian.

Details on Majorana arcs in the Andreev approximation

This section explains how the Majorana arcs are calculated within the Andreev approximation. To derive the Majorana arc surface states within the semiclassical or Andreev approximation [following Park et al. ()], we solve the BdG Hamiltonian (Eq. 35) in the presence of a spatially dependent pairing potential Δ(r). Specifically, we assume that Δ(r) is given by Δ(r) = ΔΘ(y), that is, a surface boundary in the xz plane. Substituting k → −i∇, the first-quantized BdG equation reads as

Here, Ψ(r) is a (first-quantized) spinor wave function, which we further decomposed into , where Ψ0 is a spinor, k∥ = (k, k) is the momentum parallel to the boundary surface, which is a good quantum number, and the functions take the general form

Here, N is a normalization constant, the parallel coordinates are given by r∥ = (x, z), and is defined as the semiclassical momentum perpendicular to the boundary surface. χ(y) is a scalar function. We demand that the wave functions satisfy the Dirichlet boundary condition at y = 0 and y → ∞, that is, Ψ(r)| = 0 and Ψ(r)| = 0. The function χ(y) will always be such that the latter is satisfied, and to satisfy the former, we take superpositions of incident and reflected waves.

In the semiclassical approximation, defined by the condition k⊥ ≫ kF(Δ0/εF), the BdG equation (Eq. 39) takes the form of the Andreev equation for χ(y), given by

Here, the Hamiltonian H0(k∥) depends only on the momenta k∥, and H⊥(k⊥,−i∇) is a function of and the semiclassical momentum k⊥. Our strategy will be to solve Eq. 41 for the case k∥ = 0 and E = 0 and then obtain the effective Hamiltonian of the Majorana arcs in the vicinity of k∥ = 0 by projecting H0(k∥) into the space of zero-energy solutions of H⊥(k⊥,−i∇).