Topological nodal i-wave superconductivity in PtBi2

DFT calculations and BdG model

We performed DFT calculations using the full-potential local orbital code FPLO³⁸ within the generalized gradient approximation³⁹, using the tetrahedron method with 12³ points for the Brillouin integration. Subsequently, a maximally projected symmetry-conserving Wannier model⁴⁰ containing Wannier orbitals for each Bi 6p and Pt 6s, 5d basis orbital was constructed. The model is mapped onto a semi-infinite slab with a surface block consisting of three Bi₆Pt₃ layers on which only a nonzero gap function is added. In detail, the gap function reads

with \(D(k)=1,020{k}_{x}{k}_{y}(3{k}_{x}^{4}-10{k}_{x}^{2}{k}_{y}^{2}+3{k}_{y}^{4})\) being a scaled Taylor expansion of sin(6ϕ). The scaling factor was chosen to fulfil D(k) = 1 at k = (0.4, 0.0325) Å⁻¹. The Bloch spectral density for a penetration depth of three surface blocks is obtained by solving the BdG equations using Green function recursion^2,41 in this semi-infinite geometry.

In Extended Data Fig. 1, we show results for the SC gap on the (001) surface as a function of the distance from the node for different values of V₀ (compare with Fig. 3b). Also, for this termination, we observe the same features: a V-shaped gap that increases with V₀.

Symmetry-allowed SC states

As noted in the main text, the symmetry of possible SC order parameters can be classified in terms of the irreps A₁, A₂ and E of the point group C_3v of PtBi₂. For trigonal PtBi₂, the mirror planes contain the Γ–K lines and not the Γ–M lines, and thus they do not include the centre points of the Fermi arcs, at which the apparent nodes are located. However, time-reversal symmetry acts like twofold rotation symmetry about the z-axis for momenta k = (k_x, k_y). Consequently, any time-reversal-symmetric function of k that is even under all mirror reflections of the lattice is also even, and that which is odd is also odd, under mirror reflections with respect to vertical planes through the arc centres.

Owing to the lack of inversion symmetry, spin-singlet and spin-triplet pairing generically mix, and it is necessary to consider 2 × 2 pairing matrices Δ(k) appearing in the BdG Hamiltonian

To construct possible pairing matrices Δ(k), we have to consider the symmetry properties of k-dependent form factors and of matrices acting on spin space.

The k-dependent factors are relevant only on the Fermi lines. To parameterize them, we start from the polar angle ϕ of k. Electronic bands are continuous, so that the Fermi arcs are connected by the Fermi surfaces of bulk states. As the Fermi arcs are horseshoe-shaped, the Fermi surface is not convex, and the points on the arcs are not uniquely labelled by ϕ. This can be resolved by deforming the parameterization without changing its symmetry and is irrelevant for our analysis. The lowest-order basis functions of ϕ together with their irrep and their sign under time reversal are listed in Extended Data Table 1. Basis functions of higher order modulate only solutions that can already be constructed from them.

A basis of the space of 2 × 2 matrices is given by the identity matrix σ₀ and the Pauli matrices σ_x, σ_y and σ_z. These transform as irreducible tensor operators of the irreps, as shown in Extended Data Table 2.

By multiplying the form factors and the basis matrices, we can obtain all possible SC states. Using standard rules for products of irreps, we can choose them to belong to specific irreps. However, only products that are even under time reversal satisfy fermionic antisymmetry Δ^T(−k) = −Δ(k) (ref. ⁴²).

The SC state of full symmetry, that is, belonging to A₁, was considered in ref. ²⁷. The irrep A₁ is even under all mirror reflections and rotations, and thus A₁ symmetry does not impose any gap nodes.

The irrep E is 2D, leading to a two-component SC order parameter. The first component, by itself as well as any symmetry-related SC states, is odd under some mirror reflections. This imposes nodes at some arc centres but not at all of them, thereby breaking the threefold rotation symmetry. The second component, by itself and symmetry-related states, does not impose any nodes. This follows from the observation that any order parameter belonging to the second component of E can be constructed from an A₁ (full symmetry) order parameter by multiplication by cos 2ϕ, which does not have zeros on the arcs. The gap magnitude is not the same at all arcs, so that such a state also breaks rotation symmetry. We do not find experimental indications for this symmetry breaking. There are also time-reversal-symmetry-breaking E states, constructed by superposition of the two components with a phase shift of ±π/2. These states do not have nodes. Moreover, there is no experimental evidence for broken time-reversal symmetry.

The remaining irrep A₂ is odd under all mirror reflections. This fact imposes nodes at all arc centres. In particular, it is odd under ϕ ↦ −ϕ. Owing to the preserved rotation and time-reversal symmetries, the gap profile is the same for all arcs. Hence, only A₂ pairing symmetry is consistent with nodes at the arc centres.

Constructing the possible pairing matrices as described above, we obtain

where the coefficients f₁, f₃, f₅ and f₆ can be chosen real because of time-reversal symmetry. U_T = iσ_y is the unitary part of the time-reversal operator. The four terms describe p-wave (l = 1) in-plane spin-triplet pairing, f-wave (l = 3) out-of-plane spin-triplet pairing, h-wave (l = 5) in-plane spin-triplet pairing and i-wave (l = 6) spin-singlet pairing, respectively. The SC energy gap on the Fermi arcs is a real function of ϕ. As the BdG Hamiltonian with the pairing matrix in equation (2) preserves threefold rotation and time-reversal symmetries, so does the energy gap. Moreover, we have shown above that it must be odd under reflection at the arc centres. Hence, the energy gap is a basis function belonging to the irrep A₂, proportional to sin 6ϕ plus higher harmonics of the same symmetry, that is, it has i-wave form.

Anomalous topological superconductivity

Two-dimensional superconducting systems that obey time-reversal symmetry (class DIII in the Altland–Zirnbauer table³¹) may host gapless Majorana cones in their 2D Brillouin zone. These gapless points are in many respects analogous to the Weyl points of 3D crystals, with the important exception that they require chiral symmetry to remain protected. They are characterized by an integer topological invariant (here, a winding number instead of the Chern number associated with Weyl points) and always occur in pairs. Furthermore, time-reversal symmetry relates cones with opposite momenta and the same integer invariant, just as with Weyl cones. Thus, as the total topological invariant associated with all band crossings must vanish in a periodic Brillouin zone, Majorana and Weyl cones must come in multiples of four as long as time-reversal symmetry is preserved.

In PtBi₂, instead, our results indicate the presence of six Majorana cones on a single superconducting surface, violating the above requirement. Moreover, all six cones have the winding numbers of the same sign, such that their sum does not add up to zero. The resolution of this apparent paradox is that other Majorana cones occur on the opposite surface of the crystal. Taking both surfaces into account, the total number of Majorana cones is 12 (a multiple of four), and the sum of all winding numbers vanishes.

This is the sense in which we state that the 2D superconductor forming on the surface of PtBi₂ is anomalous. Given their multiplicity and winding number, their Majorana cones cannot occur in a purely 2D, standalone superconducting state, but only as surface modes of a higher-dimensional, 3D system. This is analogous to the unidirectional edge modes of quantum Hall systems, which cannot be realized as standalone 1D systems but only as edge states.

PtBi₂ tight-binding model calculations

We explore the consequences of surface Majorana cones using the toy model recently introduced in ref. ²⁷, the properties of which we briefly summarize below. It is defined on a trigonal lattice, with Bravais vectors a₁ = (0, 1, 0), \({{\bf{a}}}_{2}=(\sqrt{3}/2,-1/2,0)\), a₃ = (0, 0, 1), and consists of two spinful orbitals per unit cell. Setting k_i = a_i ⋅ k, the momentum-space Hamiltonian reads

with Γ₁ = τ_zσ₀, Γ₂ = τ_xσ_x, Γ₃ = τ_yσ₀, where Pauli matrices τ_x, τ_y and τ_z encode the orbital and Pauli matrices σ_x, σ_y and σ_z denote spin, and

where

represents threefold rotations around a₃, that is, around the z-axis. Besides the threefold rotation,

the model also reproduces the other lattice symmetries of PtBi₂. There is a mirror symmetry along the k₁ = −2k₂ plane of the Brillouin zone, with M₁ = iτ₀σ_x,

and the two other mirror planes are obtained by applying the threefold rotation symmetry. Furthermore, the model obeys time-reversal symmetry \(T={\rm{i}}{\tau }_{0}{\sigma }_{y}{\mathcal{K}}\) with the complex conjugation \({\mathcal{K}}\) such that TH(k₁, k₂, k₃)T⁻¹ = H^*(−k₁, −k₂, −k₃).

We choose the hopping amplitude t = 1 as the energy unit and express all other energy scales relative to it. In all numerical simulations, we have used the parameters μ = 2, β = −0.75, λ = 3, α = 0.75 and γ = 0.5, for which the band structure hosts 12 well-separated Weyl cones, thus reproducing the behaviour of PtBi₂.

The numerical results shown in Extended Data Fig. 2 are obtained for systems that have a finite number of unit cells in the z-direction (either in a slab geometry or in a prism geometry), by adding SC pairing to the two top-most and the two bottom-most unit cells. The BdG Hamiltonian takes the same block form as equation (1), with the upper-diagonal block given by equation (3). As all pairing terms with A₂ symmetry produce nodes along the Γ–M directions, we choose the simplest, p-wave term in the toy model, resulting in

where

Linearizing equation (8) produces the same type of pairing matrix as the p-wave term of equation (2). As mentioned above, the amplitude f₁(z) = 2 for the two top-most and bottom-most unit cells, whereas f₁ = 0 otherwise. We chose such a large value for the pairing term to enhance the gap along the Fermi arc, thus reducing finite-size effects and enabling us to visualize Majorana states for numerically accessible system sizes.

Finally, to better differentiate between the Fermi arcs on the top and bottom surfaces, we shift the value of μ on the bottom-most unit cell (z = 0) from μ = 2 to μ = 1.7. This causes the top and bottom Fermi arcs to occur at different momenta in the slab Brillouin zone, such that they can be more easily visualized.

Extended Data Fig. 2a,b are obtained in a slab geometry consisting of 80 unit cells in the z-direction, in which the Cartesian momentum directions are defined as k_y = k₁ and \({k}_{x}=({k}_{1}+2{k}_{2})/\sqrt{3}\). The gap is always computed as the absolute value of the eigenenergy closest to the Fermi level, E_F = 0.

For each surface Majorana cone, we determine the winding number by using the standard approach of rotating the Hamiltonian to a block off-diagonal form. We use the full slab Hamiltonian, enabling us to determine the winding numbers of Majorana cones on both the top and the bottom surfaces.

To minimize finite-size effects, Extended Data Fig. 2c is obtained for a slab consisting of 1,000 unit cells in the z-direction and shows the gap opening in the Majorana cone on the top surface, indicated by an orange arrow in Extended Data Eig. 2b. The Zeeman field is included by adding an onsite term V_zτ₀σ_z to the Hamiltonian equation (3).

Extended Data Fig. 4d is reproduced from Fig. 3c and shows the coexistence of bulk Weyl cones, surface Majorana cones and Majorana hinge modes. To reduce the finite-size gap that would otherwise be present in these features, the panel includes results from two separate simulations. The green, red and blue points are obtained in a slab geometry (infinite along both k_x and k_y), and consist of 320 unit cells along the z-direction. Eigenvalues are computed for different values of k_x, in steps of 10⁻³, and plotted as a function of k_y ∈ [1.05, 1.45]. The colour represents the state probability density summed over the bottom half of the slab, meaning for values z < 160. The apparently sharp transition between the different colours is a consequence of our plotting choice. Multiple points of different colours are plotted on top of each other, and the only visible colours correspond to those points that are plotted last. The hinge states, shown in black, are then superimposed on the slab band structure plot. The hinge modes are obtained in a prism geometry, infinite along y, 40 unit cells along z and 300 unit cells along x. This is also the geometry of the system that is used in Fig. 3e, where the colour scale denotes the probability density summed over the four states closest to E = 0 at k_y = 1.3, as marked by the green arrow in Extended Data Fig. 4d.

Leading edge, peak and trailing edge positions of the energy distribution curves along the arc of PtBi₂

As mentioned in the main text, compared with the leading edge, the energy distribution curve peak, and more significantly, the trailing edge positions at higher binding energies are more accurate for determining the gap function near the node. In Extended Data Fig. 3a, the trailing edge positions are overlaid on the Fermi surface map for visualization. For comparison, the leading edge, peak and trailing edge positions of the energy distribution curves are plotted as a function of distance from the node in Extended Data Fig. 3d. It is to be noted that the SC gap is plotted with respect to the distance from the node in the calculations (Fig. 3b), whereas the gap function is presented with respect to θ in Fig. 2. To draw one-to-one correspondence, we have plotted the trailing edge as a function of both θ and distance from the node in Extended Data Fig. 3c. The dip in the trailing edge position near the node is much sharper compared with the leading edge position and agrees more with the gap function obtained from calculations near the node (Fig. 3a).

Leading edge, peak and trailing edge positions for BSCCO

To demonstrate how the LEG curves are affected by the Fermi function, we have plotted the leading edge, peak and trailing edge positions for BSCCO (bismuth strontium calcium copper oxide), which is a well-established nodal superconductor. Extended Data Fig. 4a–c shows the overall good quality of the collected ARPES data. Similar to PtBi₂, the trailing edge position of BSCCO harbours a much sharper feature near the node compared with the leading edge and peak positions (Extended Data Fig. 4d).