# New type of helical topological superconducting pairing
at finite excitation energies

Masoud Bahari<EMAIL_ADDRESS>Institute for Theoretical
Physics and Astrophysics, University of Würzburg, D-97074 Würzburg, Germany
Würzburg-Dresden Cluster of Excellence ct.qmat, Germany  Song-Bo Zhang
Department of Physics, University of Zürich, Winterthurerstrasse 190, 8057,
Zürich, Switzerland  Chang-An Li  Institute for Theoretical Physics and
Astrophysics, University of Würzburg, D-97074 Würzburg, Germany  Würzburg-
Dresden Cluster of Excellence ct.qmat, Germany  Sang-Jun Choi  Institute for
Theoretical Physics and Astrophysics, University of Würzburg, D-97074
Würzburg, Germany  Würzburg-Dresden Cluster of Excellence ct.qmat, Germany
Carsten Timm  Würzburg-Dresden Cluster of Excellence ct.qmat, Germany
Institute of Theoretical Physics, Technische Universität Dresden, 01062
Dresden, Germany  Björn Trauzettel  Institute for Theoretical Physics and
Astrophysics, University of Würzburg, D-97074 Würzburg, Germany  Würzburg-
Dresden Cluster of Excellence ct.qmat, Germany

(September 3, 2024)

###### Abstract

We propose a new type of helical topological superconductivity away from the
Fermi surface in three-dimensional time-reversal-symmetric odd-parity
multiband superconductors. In these systems, pairing between electrons
originating from different bands is responsible for the corresponding
topological phase transition. Consequently, a pair of helical topological
Dirac surface states emerges at finite excitation energies. These helical
Dirac surface states are tunable in energy by chemical potential and strength
of band-splitting. They are protected by pseudospin rotation symmetry. We
suggest concrete materials in which this phenomenon could be observed.

Introduction.—The search for intrinsic topological superconductors has been an
active research area in the condensed matter community for many years
(Bernevig, 2013; Hasan  _et al._ , 2015; Ando and Fu, 2015; Sato and Ando,
2017; Shen, 2017). In weakly interacting systems, the essential ingredient for
intrinsic topological phase transitions (TPTs) is a weak attractive electron-
electron interaction appearing close to the Fermi surface. This results in
particular bulk properties including topological superconducting gaps (Fu and
Kane, 2008; Qi  _et al._ , 2009; Sato, 2009; Fu and Berg, 2010; Hor  _et al._
, 2010; Wray  _et al._ , 2010; Qi and Zhang, 2011; Fang  _et al._ , 2015; Yang
_et al._ , 2016; Wang  _et al._ , 2018; Khalaf, 2018; Călugăru  _et al._ ,
2019; Zhang  _et al._ , 2020), possibly with nodes (Joynt and Taillefer,
2002b; Sato, 2006; Béri, 2010; Kobayashi  _et al._ , 2014; Schnyder and
Brydon, 2015; Agterberg  _et al._ , 2017; Kawakami  _et al._ , 2018; Tamura
_et al._ , 2020; Nayak  _et al._ , 2021). Due to particle-hole symmetry, the
bulk-boundary correspondence allows for emergent Majorana boundary states
protected by band topology.

Finite energy (FE) Cooper pairing is a particular type of pairing that can
happen away from the Fermi surface (Moreo  _et al._ , 2009; Komendová  _et
al._ , 2015; Triola and Balatsky, 2017; Linder and Balatsky, 2019; Tang  _et
al._ , 2021; Bahari  _et al._ , 2022; Kanasugi and Yanase, 2022; Chakraborty
and Black-Schaffer, 2022; Bahari  _et al._ , 2022; Rüßmann  _et al._ , 2023).
It arises from multiband effects in superconductors (Binnig  _et al._ , 1980;
Shan  _et al._ , 2011; Lin  _et al._ , 2014; Kittaka  _et al._ , 2014; Nomoto
_et al._ , 2016; Nomoto and Ikeda, 2016b; Seo  _et al._ , 2017; Zhao  _et al._
, 2018; Singh  _et al._ , 2018; Menke  _et al._ , 2019; Kanasugi and Yanase,
2019; Xu  _et al._ , 2019; Lado and Sigrist, 2019; Hu  _et al._ , 2020; Timm
and Bhattacharya, 2021; Shang  _et al._ , 2021; Wang  _et al._ , 2021a; Sevik
_et al._ , 2022). Dips in the density of states at FEs is a particular
spectral feature of such pairing found in hybrid structures based on the
superconducting proximity effect of $\text{NbSe}_{2}$ in the topological
insulators $\text{\text{BiSbTe}}_{1.25}\text{Se}_{1.75}$ (Banerjee  _et al._ ,
2018) and $\text{Bi}_{2}\text{Se}_{3}$ (Li  _et al._ , 2017). In this Letter,
we show that TPTs can be exclusively induced by FE Cooper pairing emerging
from time-reversal-symmetric odd-parity multiband superconductors. The
interplay of spin-orbit coupling and electron-phonon interaction allows for
the formation of unconventional FE Cooper pairing with topologically
nontrivial order. Bulk band crossings are partially gapped out in momentum
space at FEs due to odd-parity superconductivity. This leads to TPTs as
evidenced by helical topological Dirac surface states at finite excitation
energies. Such surface states are composed of electrons and holes with
different magnetic quantum numbers. They are distinct from conventional
Majorana boundary states due to broken particle-hole symmetry away from the
Fermi energy. The helical Dirac points are (i) protected topologically by
pseudospin rotation symmetry originating from parity-time-reversal symmetry as
we specify below, and (ii) tunable in energy by chemical potential and
strength of band splitting.

To observe FE Cooper pairing, the normal state should have at least two energy
bands with the same sign of curvature close to the Fermi energy. The bands
need to have different effective masses. Therefore, our results are relevant
for a broad range of multiband superconductors as we elaborate on below.

To capture the underlying physics, we develop a theory for generic
superconducting multiband systems with four-valued local degrees of freedom.
For concreteness, we work with a model Hamiltonian describing the band
structure of electrons with effective angular momentum $j=3/2$ near the Fermi
energy. However, the formalism is not restricted to this particular choice of
angular momentum quantum numbers.

Model.—We start with the Luttinger-Kohn Hamiltonian $\hat{\mathcal{H}}_{{\bf
k}}=\alpha|{\bf k}|^{2}+\beta({\bf k}\cdot\mathbf{J})^{2}-\mu$, which
describes $j=3/2$ electrons within the $\Gamma_{8}$ bands (Agterberg  _et al._
, 2017; Brydon  _et al._ , 2016; Luttinger and Kohn, 1955; Dresselhaus, 1955).
The kinetic, symmetric spin-orbit-coupling, and chemical potential terms are
parameterized by $\alpha$, $\beta$, and $\mu$, respectively. ${\bf
k}=(k_{x},k_{y},k_{z})$ is the three-dimensional (3D) momentum and
$\mathbf{J}=(\hat{J}_{x},\hat{J}_{y},\hat{J}_{z})$ are the $4\times 4$ spin
matrices in $j=3/2$ formalism. The normal state has a pair of doubly
degenerate quadratic energy bands given by $E_{{\bf
k}}^{+}=(\alpha+9\beta/4)|{\bf k}|^{2}-\mu$ and $E_{{\bf
k}}^{-}=(\alpha+\beta/4)|{\bf k}|^{2}-\mu$. The chemical potential $\mu$
should fulfill $|\mu|<N\hbar\omega_{D}$ where $\omega_{D}$ is the Debye
frequency and $N=1.25+\alpha/\beta$ is a material dependent number. While
$\hat{\mathcal{H}}_{{\bf k}}$ preserves O(3) symmetry, our predictions are
generic and can be applied to systems with relative point group symmetry in
the normal state (Com, j).

Figure 1:  (a) BdG energy spectra along the [001] direction for three
different values of the chemical potential. Solid (dashed) lines depict
twofold degenerate particle (hole) bands. Energy bands plotted with distinct
colors (blue and black) are decoupled due to pseudospin rotation symmetry. In
the third panel, the fade gray lines denote $\Delta=0$. (b) Spectrum versus
$\mu$ for a (001) slab with $80$ layers. The color bar indicates the inverse
participation ratio (Com, b). The gray and blue colors represent bulk states,
and red color denotes edge states. (c) Topological phase diagram induced by FE
pairing. Region III (II) indicates that the helical Dirac points are located
(coexist with) within the low energy gap (bulk states). (d) Topological
invariant corresponding to panel (b). $A_{2u}$ pairing is responsible for
superconductivity in all panels. Other parameters are $\beta=0.3\alpha$,
$\alpha=-\Delta$ and $k_{0}=a^{-1}$ with $a$ being the lattice constant in
corresponding tight-binding calculations.

The superconducting properties of the system can be captured by the
Bogoliubov-de Gennes (BdG) Hamiltonian

$\hat{H}({\bf k})=\left(\begin{array}[]{cc}\hat{\mathcal{H}}_{{\bf
k}}&\hat{\Delta}_{{\bf k}}\\\ \hat{\Delta}_{{\bf
k}}^{\dagger}&-\hat{\mathcal{H}}_{-{\bf k}}^{T}\end{array}\right),$ (1)

where $\hat{\mathcal{H}}_{{\bf k}}$ is the inversion-symmetric normal part and
$\hat{\Delta}_{{\bf k}}$ is a $4\times 4$ pairing matrix describing pairing of
electrons with four-valued degrees of freedom. Generally, the explicit form of
$\hat{\Delta}_{{\bf k}}$ depends on the basis of the system. For concreteness,
we focus on systems with cubic point group $O_{h}$. In this case,
$\hat{\Delta}_{{\bf k}}$ can be expanded into irreducible tensor operators of
irreducible representations (irreps) of $O_{h}$ (Savary  _et al._ , 2017;
Venderbos  _et al._ , 2018; Yu and Liu, 2018; Bahari  _et al._ , 2022).

Topology at finite energies.—The prerequisite for TPTs away from the Fermi
energy is that the $E_{{\bf k}}^{+}$ electron and $-E_{{\bf k}}^{-}$ hole
bands cross each other. These band crossings become anti-crossings in the
presence of odd-parity FE Cooper pairing as we explain in detail in the
Supplemental Material (SM) (Sup, ). For concreteness, we consider the $A_{2u}$
($p$-wave septet) pairing given by $\hat{\Delta}_{{\bf
k}}=\Delta(\mathbf{k}\cdot\hat{\mathbf{T}})\hat{\mathcal{R}}$, where $\Delta$
is the pairing strength, $\hat{\mathcal{R}}=e^{i\pi\hat{J}_{y}}$ is the
fermionic antisymmetry factor, and
$\hat{T}_{i}=\\{\hat{J}_{i},\hat{J}_{i+1}^{2}-\hat{J}_{i+2}^{2}\\}$ with
$i+1=y$ if $i=x$, etc., cyclically (Brydon  _et al._ , 2016; Savary  _et al._
, 2017). Note that $p$-wave septet pairing has been argued to be the favorable
pairing symmetry in systems with cubic symmetry (Brydon  _et al._ , 2016; Kim
_et al._ , 2018; Ishihara  _et al._ , 2021). Although we focus on this
specific pairing channel, our qualitative results are general. They apply to
several odd-parity pairing channels summarized in Table. 1.

$O_{h}$ |  | ${\bf k}$ | $\hat{\Delta}_{{\bf k}}^{+-}=\Delta\times$
---|---|---|---
$E_{u}$ | $E_{u}^{(2)}$ | $(0,0,k_{z})$ | $k_{z}\hat{\tau}_{0}$
$T_{2u}$ | $T_{2u}^{(1)}$ | $(0,0,k_{z})$ | $k_{z}\hat{\tau}_{z}$
| $T_{2u}^{(2)}$ | $(0,k_{y},0)$ | $k_{y}\hat{\tau}_{x}$
| $T_{2u}^{(3)}$ | $(k_{x},0,0)$ | $ik_{x}\hat{\tau}_{x}$
$A_{2u}$ |  | $(k_{x},0,0)$ | $-k_{x}\hat{\tau}_{0}$
|  | $(0,k_{y},0)$ | $-ik_{y}\hat{\tau}_{z}$
|  | $(0,0,k_{z})$ | $k_{z}\hat{\tau}_{0}$
$T_{2u}$ | $T_{2u}^{(1)}$ | $(0,0,k_{z})$ | $k_{z}\hat{\tau}_{z}$
|  | $(\pm k,\pm k,0)$ | $k(\pm\hat{\tau}_{0}\pm i\hat{\tau}_{z})$
| $T_{2u}^{(2)}$ | $(k_{x},0,0)$ | $k_{x}\hat{\tau}_{x}$
|  | $(0,\pm k,\pm k)$ | $k(\hat{\tau}_{y}\pm\hat{\tau}_{z}/3)$
| $T_{2u}^{(3)}$ | $(0,k_{y},0)$ | $k_{y}\hat{\tau}_{x}$
|  | $(\pm k,0,\pm k)$ | $k(\pm\hat{\tau}_{0}/3-i\hat{\tau}_{y})$
Table 1: TPTs at FEs in $p$-wave pairing channels distinguished by cubic
irreps (first column) and their components (second column). The direction of
TPTs is denoted in the third column. The last column indicates the FE pairing
potential responsible for TPTs. The pairing matrices for the given irreps are
given in the SM (Sup, ; Com, j).

Without loss of generality, we illustrate our predictions for the [001]
direction (i.e., $k_{x}=k_{y}=0$). We observe in Fig. 1(a) that the bulk
bands, i.e., electron $E_{k_{z}}^{+}$ and hole $-E_{k_{z}}^{-}$ branches, do
not touch each other for $\mu=-5.7\Delta$ close to $k_{z}=\pi$ away from the
Fermi surface. In this case, the system exhibits a topologically trivial
phase. Increasing slightly the Fermi energy to $\mu=-5.5\Delta$, this makes
the two bands touch each other. For $\mu=-5.2\Delta$, the presence of FE
pairing, i.e., pairing between electrons with different magnetic quantum
numbers ($|m_{j}|=3/2$ with $|m_{j}|=1/2$), opens a band gap indicating a TPT.
This TPT is accompanied by the emergence of helical topological surface
states. To illustrate the appearance of helical surface states (HSSs), the
spectrum of a finite system is depicted in Fig. 1(b) (Com, c). The localized
fourfold degenerate HSSs are marked by red color indicating large inverse
participation ratio (IPR) (Com, b). They emerge within the range of Fermi
energies $\mu/\Delta\in[-5.5,0]$, in which a nonzero topological invariant,
see Fig. 1(d), can be found. We describe below how the topological invariant
is defined. Generally, when the energy bands have downward (upward) curvature,
i.e., $\text{sgn}(\alpha)=\text{sgn}(\beta)=-(+)1$, the topological phase
occurs in the range $\mu/\Delta\in[\mu_{c},0]\ ([0,\mu_{c}])$ with the
critical value $\mu_{c}=4\alpha+5\beta$. The analysis of the spectrum along
other directions reveals that these HSSs are helical topological Dirac points
occurring at finite excitation energies as shown in Fig. 2. Notably, the
helical surface points do not hybridize with the bulk states due to the
presence of pseudospin rotation symmetry (Com, d).

According to the topological phase diagram, depicted in the plane of
$(\mu/\alpha,\beta/\alpha)$ in Fig. 1(c), the FE Cooper pairing in odd-parity
multiband superconductors is intrinsically topological (regions II and III)
within a certain range of Fermi energies.

Theory for topological phase transition.—The nontrivial topology at FEs
originates from unconventional pairing of FE electrons. To understand the
underlying mechanism, we project the BdG Hamiltonian onto the doubly
degenerate eigenspinors of the normal state $\hat{V}_{{\bf k}}^{\pm}$ as the
pseudospin basis. In this basis, the inter-band representation of the BdG
Hamiltonian becomes

$\hat{\mathscr{H}}({\bf k})=\Bigg{(}\begin{array}[]{cc}\hat{\mathsf{H}}_{{\bf
k}}^{+-}&\hat{\Delta}_{{\bf k}}^{\text{intra}}\\\ \big{(}\hat{\Delta}_{{\bf
k}}^{\text{intra}}\big{)}^{\dagger}&\hat{\mathsf{H}}_{{\bf
k}}^{-+}\end{array}\Bigg{)},$ (2)

where $\hat{\mathsf{H}}_{{\bf k}}^{+-}$ and $\hat{\mathsf{H}}_{{\bf k}}^{-+}$
denote the inter-band parts of the BdG Hamiltonian and $\hat{\Delta}_{{\bf
k}}^{\text{intra}}$ describes intra-band pairing at the Fermi energy. The
subblock matrices in Eq. $\eqref{InterBandH}$ can be written as

$\\!\\!\\!\hat{\mathsf{H}}_{{\bf
k}}^{+-}=\left(\\!\\!\\!\begin{array}[]{cc}E_{{\bf k}}^{+}&\hat{\Delta}_{{\bf
k}}^{+-}\\\ (\hat{\Delta}_{{\bf k}}^{+-})^{\dagger}&-E_{{\bf
k}}^{-}\end{array}\\!\\!\\!\right),\ \ \hat{\Delta}_{{\bf
k}}^{\text{intra}}\\!=\\!\left(\\!\\!\\!\begin{array}[]{cc}0&\hat{\Delta}_{{\bf
k}}^{++}\\\ (\hat{\Delta}_{{\bf
k}}^{--})^{\dagger}&0\end{array}\\!\\!\\!\right),\\!\\!$ (3)

where $\hat{\Delta}_{{\bf k}}^{\nu\nu^{\prime}}$ with
$\nu,\nu^{\prime}\in\\{+,-\\}$ denotes the pairing potential projected onto
the normal state bands defined by $\hat{\Delta}_{{\bf
k}}^{\nu\nu^{\prime}}=\hat{V}_{{\bf k}}^{\nu\dagger}\hat{\Delta}_{{\bf
k}}(\hat{V}_{-{\bf k}}^{\nu^{\prime}\dagger})^{T}$ where $\nu\neq\nu^{\prime}$
($\nu=\nu^{\prime}$) indicates inter- (intra-) band pairing between electrons
having different (identical) magnetic quantum numbers in magnitude.

Figure 2:  Spectra of a slab with (001) surfaces as a functions of $k_{y}$
indicating (a) topologically trivial phase, (b) topological phase transition,
and (c) topological phase at FEs. The thickness is 200 layers. (d) 3D energy
dispersion versus $k_{x}$ and $k_{y}$. (e) Helical surface cones at positive
and negative energies related by particle-hole symmetry. (d,e) Yellow arrows
indicate left- (right-) handed spin texture in the upper (lower) cone.
$E_{0}=0.1\Delta$ and other parameters are the same as in Fig. 1.

In our model, the most convenient way to explain TPTs at FEs is to make
$\hat{\mathscr{H}}({\bf k})$ block diagonal. Thus, TPTs occur in directions in
momentum space ,where the odd-parity inter-band pairing is finite, whereas the
intra-band pairing vanishes (nodes at the Fermi energy), i.e.,
$\hat{\Delta}_{{\bf k}}^{\nu\nu^{\prime}}\neq 0$ and $\hat{\Delta}_{{\bf
k}}^{\text{intra}}=0$. In this case, Eq. $\eqref{InterBandH}$ becomes block
diagonal and $\hat{\mathsf{H}}_{{\bf k}}^{\nu\nu^{\prime}}$ is further
reducible. To block diagonalize $\hat{\mathsf{H}}_{{\bf
k}}^{\nu\nu^{\prime}}$, we define a rotation symmetry operator in inter-band
basis

$\mathcal{\hat{D}}_{\mathbf{n}}(\theta,\epsilon)\equiv\text{diag}\big{(}e^{i\theta/2(\mathbf{n}\cdot\hat{\boldsymbol{\tau}})},\epsilon
e^{-i\theta/2(\mathbf{n}\cdot\hat{\boldsymbol{\tau}})^{*}}),$ (4)

where $\theta$ is the rotation angle about an arbitrary unit vector
$\mathbf{n}$. The vector of Pauli matrices in pseudospin basis is
$\hat{\boldsymbol{\tau}}=(\hat{\tau}_{x},\hat{\tau}_{y},\hat{\tau}_{z})$ and
$\epsilon=\pm 1$ stands for two unitarily equivalent representations.
Interestingly, the superconducting Hamiltonian of each block preserves a
symmetry under pseudospin rotation by $\theta=\pi$, i.e.,
$[\mathcal{\hat{D}}_{\mathbf{n}}(\pi,\epsilon),\hat{\mathsf{H}}_{{\bf
k}}^{\nu\nu^{\prime}}]=0$. This allows us to label the corresponding energy
bands by the sign $\pm$ corresponding to the eigenvalues $\lambda=\pm i$ of
$\mathcal{\hat{D}}_{\mathbf{n}}(\pi,\epsilon)$. This pseudospin rotation
symmetry fulfills the property
$\mathcal{\hat{D}}_{\mathbf{n}}^{2}(\pi,\epsilon)=-\hat{\mathbbm{1}}$. It
appears when the system has at least four (or more) degrees of freedom
combined with inversion and time-reversal symmetries Suh  _et al._ (2023). The
basis elements of $\hat{\mathsf{H}}_{{\bf k}}^{+-}$ under $\pi$-rotation about
the $x$ ($y$)-axis flips the orientation of the pseudospin indices
$\\{\uparrow,\downarrow\\}\leftrightarrow\\{\downarrow,\uparrow\\}$ while the
momentum remains intact. In our model, the pseudospin up (down) is composed of
a linear combination of magnetic quantum numbers $m_{j}\in\\{\pm 1/2,\pm
3/2\\}$. This symmetry results in the conservation of pseudospin quantum
numbers.

$\hat{\mathsf{H}}_{{\bf k}}^{\nu\nu^{\prime}}$ can be reduced to two
independent subblock matrices in the eigenspaces of
$\mathcal{\hat{D}}_{\mathbf{n}}(\pi,\epsilon)$. The corresponding similarity
transformation is $\hat{\mathcal{U}}^{\dagger}\hat{\mathsf{H}}_{{\bf
k}}^{\nu\nu^{\prime}}\hat{\mathcal{U}}=\text{diag}\big{(}\hat{h}_{{\bf
k},-}^{\nu\nu^{\prime}},\hat{h}_{{\bf k},+}^{\nu\nu^{\prime}}\big{)}$, where
$\hat{\mathcal{U}}$ is the unitary matrix defined by the eigenvectors of
$\mathcal{\hat{D}}_{\mathbf{n}}(\pi,\epsilon)$. Then, $\hat{h}_{{\bf
k},\lambda}^{\nu\nu^{\prime}}$ becomes

$\displaystyle\hat{h}_{{\bf k},\lambda}^{\nu\nu^{\prime}}$
$\displaystyle=\left(\begin{array}[]{cc}-E_{{\bf
k}}^{\nu^{\prime}}&\big{(}\delta_{{\bf
k},\lambda}^{\nu\nu^{\prime}}\big{)}^{*}\\\ \delta_{{\bf
k},\lambda}^{\nu\nu^{\prime}}&E_{{\bf k}}^{\nu}\end{array}\right),$ (7)

where $\delta_{{\bf k},\lambda}^{\nu\nu^{\prime}}$ is the FE pairing
represented in the eigenbasis of
$\mathcal{\hat{D}}_{\mathbf{n}}(\pi,\epsilon)$ (Com, f). Eventually, the
transformed BdG Hamiltonian, Eq. (2), decouples into four $2\times 2$
subblocks,

$\hat{\mathscr{H}}^{\prime}({\bf k})=\text{diag}\big{(}\hat{h}_{{\bf
k},+}^{+-},\hat{h}_{{\bf k},-}^{+-},\hat{h}_{{\bf k},+}^{-+},\hat{h}_{{\bf
k},-}^{-+}\big{)}.$ (8)

The topological properties at FEs are encoded in the subblocks of
$\hat{\mathscr{H}}^{\prime}({\bf k})$. Generically, when the pairing term in
$\hat{h}_{{\bf k},\lambda}^{\nu\nu^{\prime}}$ is odd in momentum, the FE bulk
band gaps between electron and hole branches cross or anticross at the parity-
invariant momenta. This happens through the tuning of the normal state
parameters, i.e., $(\alpha,\beta,\mu)$, and leads to a topologically
nontrivial phase as described below.

The topological classification of our 3D model in class DIII is captured by a
$\mathbb{Z}$ topological index at low energies (Com, g; Fu and Kane, 2007;
Schnyder  _et al._ , 2008; Hughes  _et al._ , 2011; Morimoto and Furusaki,
2013; Chiu  _et al._ , 2016). We are, however, interested in a different type
of topological classification at finite excitation energies. This
classification is based on the sum of four $\mathbb{Z}_{2}$ indices stemming
form a quantization of the Berry phase defined for each block presented in Eq.
(8). Note that $\hat{h}_{{\bf k},-}^{+-}$ is essentially a replica of
$\hat{h}_{{\bf k},+}^{+-}$, owing to the presence of
$\mathcal{\hat{D}}_{\mathbf{n}}(\pi,\epsilon)$ symmetry. Additionally,
$\hat{h}_{{\bf k},\pm}^{-+}$ represent the particle-hole counterparts of
$\hat{h}_{{\bf k},\pm}^{+-}$. Therefore, the decoupled sectors in Eq. (8) are
topologically identical. The topological invariant is given by
$Z=\sum_{\lambda\nu\nu^{\prime}}N_{\lambda}^{\nu\nu^{\prime}}$ with
$\mathbb{Z}_{2}$ index $N_{\lambda}^{\nu\nu^{\prime}}=(1/\pi)\oint_{c}d{\bf
k}\cdot i\langle u_{{\bf
k},\lambda}^{\nu\nu^{\prime}}|\boldsymbol{\nabla}_{{\bf k}}|u_{{\bf
k},\lambda}^{\nu\nu^{\prime}}\rangle$, where $|u_{{\bf
k},\lambda}^{\nu\nu^{\prime}}\rangle$ is the eigenspinor of $\hat{h}_{{\bf
k},\lambda}^{\nu\nu^{\prime}}$ associated with the lower energy band (Berry,
1984; Zak, 1989; Resta, 1994; Xiao  _et al._ , 2010). In the topological
phase, $N_{\lambda}^{\nu\nu^{\prime}}$ is quantized due to parity symmetry
implied by $\hat{\tau}_{z}\hat{h}_{{\bf
k},\lambda}^{\nu\nu^{\prime}}\hat{\tau}^{-1}_{z}=\hat{h}_{{-\bf
k},\lambda}^{\nu\nu^{\prime}}$. We obtain analytically that $\hat{h}_{{\bf
k},\lambda}^{\nu\nu^{\prime}}$ is topologically nontrivial with
$N_{\lambda}^{\nu\nu^{\prime}}=1$, when $\mu/\Delta\in[\mu_{c},0]\
([0,\mu_{c}])$ for $\text{sgn}(\alpha)=\text{sgn}(\beta)=-(+)1$. The
topological phase diagram in Fig. 1(c) illustrates this analysis Com (h). The
topological index becomes a net value of $Z=4$ indicating the emergence of two
pairs of helical Dirac states on the surface of superconductor, cf. Fig. 2(e).
One helical pair appears at positive excitation energy while the other pair
emerges at negative excitation energy (Com, i).

2D surface states.— The HSSs disperse linearly in the vicinity of the $\Gamma$
point of the 2D surface Brillouin zone as illustrated in Fig. 2(c) and 2(d).
At larger momenta, the surface dispersion connects to the bulk states. This is
due to the intra-band pairing at the Fermi energy, i.e., $\hat{\Delta}_{{\bf
k}}^{\text{intra}}$. It is finite for nonvanishing $k_{x}$ or $k_{y}$. Hence,
$\hat{\mathscr{H}}^{\prime}({\bf k})$ is no longer block-diagonal and the
pseudospin rotation symmetry is broken. This lifts the twofold degeneracy
except at the Dirac points, which is the reason for the emergence of helical
Dirac surface cones at FEs.

To gain insight into the dispersion of the HSSs, we consider an open surface
perpendicular to the $z$ direction and derive an effective Hamiltonian for the
HSSs near the Dirac points. We consider $k_{x}$ and $k_{y}$ small and
decompose the BdG Hamiltonian as $\hat{H}({\bf
k})=\hat{H}(0,0,k_{z})+\delta\hat{H}(k_{x},k_{y},0)$. First, we derive the
Dirac point energy at $k_{x}=k_{y}=0$. At this point, $\hat{H}(0,0,k_{z})$
becomes block diagonal through the unitary transformation $\hat{\mathscr{U}}$
constructed from the pseudospin symmetry
$\mathcal{\hat{D}}_{\mathbf{n}}(\pi,\epsilon)$ in the inter-band basis. This
leads to
$\hat{\mathscr{H}}^{\prime}(0,0,k_{z})=\text{diag}(\hat{h}_{k_{z},+}^{+-},\hat{h}_{k_{z},-}^{+-},\hat{h}_{k_{z},+}^{-+},\hat{h}_{k_{z},-}^{-+})$,
which is identical to Eq. (8). We consider a semi-infinite system in the half
space $z\geq 0$ with open boundary conditions. Therefore, $k_{z}$ is no longer
conserved and we treat the $z$ direction in real space, where the momentum
operator is $-i\partial_{z}$. The eigenvalue equation for each subblock reads
$\hat{h}_{-i\partial_{z},\lambda}^{\nu\nu^{\prime}}\hat{\Phi}(\xi,z)=\mathscr{E}_{\text{DP}}\hat{\Phi}(\xi,z)$,
where $\nu,\nu^{\prime}\in\\{+,-\\}$ and
$\hat{\Phi}(\xi,z)=(u,v)^{T}\text{exp}(\xi z)$ is the trial eigenspinor of the
Dirac point with $\xi$ being its penetration length. $|u|^{2}$ ($|v|^{2}$) are
probability weights for electron (hole) bands having different band indices.
We derive the energy of the helical topological Dirac surface points induced
by unconventional FE pairing as

$\mathscr{E}_{\text{DP}}=\pm\mu\frac{m-m^{\prime}}{m+m^{\prime}},$ (9)

where $m=\alpha+\beta/4$ and $m^{\prime}=\alpha+9\beta/4$ are the masses for
the normal state bands. Notably, the Dirac points are tunable in energy by
chemical potential $\mu$ and spin-orbit-coupling strength
$m-m^{\prime}=-2\beta$ (Com, i).

Projecting the BdG Hamiltonian onto the basis of surface states results in an
effective Hamiltonian for the FE helical topological surface states of the
form (Sup, )

$\displaystyle\hat{\mathtt{H}}(k_{x},k_{y})$
$\displaystyle=(\mathscr{E}_{\text{DP}}-\varsigma_{1}{\bf
k}_{||}^{2})\hat{\sigma}_{0}+\varsigma_{2}(k_{y}\hat{\sigma}_{x}-k_{x}\hat{\sigma}_{y}),$
(10)

where ${\bf k}_{||}^{2}=k_{x}^{2}+k_{y}^{2}$, and $\varsigma_{1(2)}$ is the
group velocity of the HSSs. To leading order, the HSSs exhibit the Dirac
dispersion
$\mathcal{E}_{\pm}\approx\mathscr{E}_{\text{DP}}\pm\varsigma_{2}|\mathbf{k}_{||}|$.
In the limit, where either $\mu$ or $m-m^{\prime}$ vanishes, the Dirac points
shift to zero energy. Then, the FE helical topological surface states become
dispersive helical Majorana modes (Bernevig, 2013; Hasan  _et al._ , 2015;
Ando and Fu, 2015; Sato and Ando, 2017; Shen, 2017). In this case,
$\hat{h}_{k_{z},\pm}^{+-}$ ($\hat{h}_{k_{z},\pm}^{-+}$) resembles an ordinary
$p$-wave superconducting Hamiltonian (Kitaev, 2001). However, for finite $\mu$
and $m-m^{\prime}$, topological HSSs emerge at FE.

To select a proper material to observe such phenomenon, certain conditions
need to be satisfied. In the normal state: (i) relevant energy bands should
have the same signs of curvature close to the Fermi energy, i.e., both curving
downward or upward, (ii) parity-time-reversal symmetry should be present.
Examples are antiperovskite oxides $\text{Sr}_{3-x}\text{SnO}$ (Oudah  _et
al._ , 2016) and half-Heusler materials RPdBi with R$\,\in${Y,Dy,Tb,Sm}
(Nakajima  _et al._ , 2015; Com, a). In the superconducting state: (i) pairing
should be odd-parity, (ii) nodes should be present in the BdG spectrum at the
Fermi energy such as in $\text{UTe}_{2}$ (Metz  _et al._ , 2019), (iii)
nonvanishing pairing at FEs should be allowed. In addition to the materials
mentioned above, hybrid structures such as
$\text{\text{BiSbTe}}_{1.25}\text{Se}_{1.75}/\text{NbSe}_{2}$ (Banerjee  _et
al._ , 2018; Xi  _et al._ , 2016), and
$\text{Bi}_{2}\text{Te}_{3}/\text{NbSe}_{2}$ (Zhang  _et al._ , 2009; Xi  _et
al._ , 2016), and X/SBi with S$\,\in${YPt,LuPd} (Brydon  _et al._ , 2016; Kim
_et al._ , 2018; Ishihara  _et al._ , 2021), and X$\,\in${Si,Ge} are suitable
platforms. The excitations in silicon, germanium (Grundmann, 2006), and SBi
(Brydon  _et al._ , 2016; Ishihara  _et al._ , 2021) have $j=3/2$ character
near the Fermi level. Notably, the $\Gamma_{8}$ bands of Si and Ge have
suitable curvature. The pairing order parameter in SBi is believed to have
$p$-wave septet symmetry (Kim  _et al._ , 2018; Ishihara  _et al._ , 2021).
Therefore, the superconducting proximity effect of SBi to Si or Ge should
induce FE pairing signaled by the emergence of HSSs.

Conclusions.— We have shown that a new type of helical topological
superconducting pairing in multiband time-reversal-symmetric odd-parity
superconductors can emerge due to unconventional pairing away from the Fermi
energy. This leads to the appearance of tunable helical topological Dirac
surface states at finite energies protected by pseudospin rotation symmetry.
Our work opens new avenues to search for helical topological superconductivity
at finite excitation energies in multiband systems. Promising experimental
probes are (spin-polarized) angle-resolved photoemission (ARPES) (Hsieh  _et
al._ , 2009; Zhang  _et al._ , 2011; Evtushinsky  _et al._ , 2014; Sie  _et
al._ , 2019; Thompson  _et al._ , 2020; Liu  _et al._ , 2020, 2021; Lin  _et
al._ , 2021; Chen  _et al._ , 2022) and scanning tunneling (STS)
spectroscopies (Bode, 2003; Wiebe  _et al._ , 2004; Alpichshev  _et al._ ,
2010; Fischer  _et al._ , 2007; Oka  _et al._ , 2014; Neupane  _et al._ ,
2016; Cornils  _et al._ , 2017; Choi  _et al._ , 2017; Fan  _et al._ , 2021;
Wang  _et al._ , 2021b; Schneider  _et al._ , 2021). Both methods have pros
and cons to observe the predicted helical Dirac surface states at finite
excitation energies related to working temperature and momentum resolution.
Therefore, the method should be chosen appropriately for a given material. Our
predictions trigger several interesting future research directions, for
instance, how topological finite energy pairing manifests in mesoscopic
quantum transport.

###### Acknowledgements.

The work was supported by the DFG (SPP 1666, SFB 1170 ToCoTronics, and SFB
1143, project A04, Project-Id 247310070), and the Würzburg-Dresden Cluster of
Excellence ct.qmat, EXC 2147, Project-Id 390858490. We thank the Bavarian
Ministry of Economic Affairs, Regional Development and Energy for financial
support within the High-Tech Agenda Project “Bausteine für das Quanten
Computing auf Basis topologischer Materialen”.

## References

  * Bernevig (2013) B. A. Bernevig,  _Topological insulators and topological superconductors_ (Princeton University Press, 2013).
  * Hasan  _et al._ (2015) M. Z. Hasan, S.-Y. Xu, and G. Bian, Phys. Scr. T164, 014001 (2015).
  * Ando and Fu (2015) Y. Ando and L. Fu, Annu. Rev. Condens. Matter Phys. 6, 361 (2015).
  * Sato and Ando (2017) M. Sato and Y. Ando, Rep. Prog. Phys. 80, 076501 (2017).
  * Shen (2017) S. Shen,  _Topological Insulators: Dirac Equation in Condensed Matter_, Springer Series in Solid-State Sciences (Springer Singapore, 2017).
  * Fu and Kane (2008) L. Fu and C. L. Kane, Phys. Rev. Lett. 100, 096407 (2008).
  * Qi  _et al._ (2009) X.-L. Qi, T. L. Hughes, S. Raghu, and S.-C. Zhang, Phys. Rev. Lett. 102, 187001 (2009).
  * Sato (2009) M. Sato, Phys. Rev. B. 79, 214526 (2009).
  * Fu and Berg (2010) L. Fu and E. Berg, Phys. Rev. Lett. 105, 097001 (2010).
  * Hor  _et al._ (2010) Y. S. Hor, A. J. Williams, J. G. Checkelsky, P. Roushan, J. Seo, Q. Xu, H. W. Zandbergen, A. Yazdani, N. P. Ong, and R. J. Cava, Phys. Rev. Lett. 104, 057001 (2010).
  * Wray  _et al._ (2010) L. A. Wray, S.-Y. Xu, Y. Xia, Y. S. Hor, D. Qian, A. V. Fedorov, H. Lin, A. Bansil, R. J. Cava, and M. Z. Hasan, Nat. Phys. 6, 855 (2010).
  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Fang  _et al._ (2015) C. Fang, B. A. Bernevig, and M. J. Gilbert, Phys. Rev. B 91, 165421 (2015).
  * Yang  _et al._ (2016) W. Yang, Y. Li, and C. Wu, Phys. Rev. Lett. 117, 075301 (2016).
  * Wang  _et al._ (2018) Y. Wang, M. Lin, and T. L. Hughes, Phys. Rev. B 98, 165144 (2018).
  * Khalaf (2018) E. Khalaf, Phys. Rev. B 97, 205136 (2018).
  * Călugăru  _et al._ (2019) D. Călugăru, V. Juričić, and B. Roy, Phys. Rev. B 99, 041301 (2019).
  * Zhang  _et al._ (2020) R.-X. Zhang, Y.-T. Hsu, and S. Das Sarma, Phys. Rev. B 102, 094503 (2020).
  * Joynt and Taillefer (2002b) R. Joynt and L. Taillefer, Rev. Mod. Phys. 74, 235 (2002b).
  * Sato (2006) M. Sato, Phys. Rev. B. 73, 214502 (2006).
  * Béri (2010) B. Béri, Phys. Rev. B. 81, 134515 (2010).
  * Kobayashi  _et al._ (2014) S. Kobayashi, K. Shiozaki, Y. Tanaka, and M. Sato, Phys. Rev. B. 90, 024516 (2014).
  * Schnyder and Brydon (2015) A. P. Schnyder and P. M. R. Brydon, J. Condens. Matter Phys. 27, 243201 (2015).
  * Kawakami  _et al._ (2018) T. Kawakami, T. Okamura, S. Kobayashi, and M. Sato, Phys. Rev. X 8, 041026 (2018).
  * Tamura  _et al._ (2020) S.-T. Tamura, S. Iimura, and S. Hoshino, Phys. Rev. B 102, 024505 (2020).
  * Nayak  _et al._ (2021) A. K. Nayak, A. Steinbok, Y. Roet, J. Koo, G. Margalit, I. Feldman, A. Almoalem, A. Kanigel, G. A. Fiete, B. Yan, Y. Oreg, N. Avraham, and H. Beidenkopf, Nat. Phys. 17, 1413 (2021).
  * Agterberg  _et al._ (2017) D. F. Agterberg, P. M. R. Brydon, and C. Timm, Phys. Rev. Lett. 118, 127001 (2017).
  * Moreo  _et al._ (2009) A. Moreo, M. Daghofer, A. Nicholson, and E. Dagotto, Phys. Rev. B. 80, 104507 (2009).
  * Komendová  _et al._ (2015) L. Komendová, A. V. Balatsky, and A. M. Black-Schaffer, Phys. Rev. B 92, 094517 (2015).
  * Triola and Balatsky (2017) C. Triola and A. V. Balatsky, Phys. Rev. B 95, 224518 (2017).
  * Linder and Balatsky (2019) J. Linder and A. V. Balatsky, Rev. Mod. Phys. 91, 045005 (2019).
  * Tang  _et al._ (2021) G. Tang, C. Bruder, and W. Belzig, Phys. Rev. Lett. 126, 237001 (2021).
  * Bahari  _et al._ (2022) M. Bahari, S.-B. Zhang, and B. Trauzettel, Phys. Rev. Research 4, L012017 (2022).
  * Kanasugi and Yanase (2022) S. Kanasugi and Y. Yanase, Commun. Phys. 5, 39 (2022).
  * Chakraborty and Black-Schaffer (2022) D. Chakraborty and A. M. Black-Schaffer, Phys. Rev. B 106, 024511 (2022).
  * Rüßmann  _et al._ (2023) P. Rüßmann, M. Bahari, S. Blügel, and B. Trauzettel, (2023), arXiv:2307.13990 .
  * Binnig  _et al._ (1980) G. Binnig, A. Baratoff, H. E. Hoenig, and J. G. Bednorz, Phys. Rev. Lett. 45, 1352 (1980).
  * Shan  _et al._ (2011) L. Shan, Y.-L. Wang, B. Shen, B. Zeng, Y. Huang, A. Li, D. Wang, H. Yang, C. Ren, Q.-H. Wang, S. H. Pan, and H.-H. Wen, Nature Physics 7, 325 (2011).
  * Lin  _et al._ (2014) X. Lin, G. Bridoux, A. Gourgout, G. Seyfarth, S. Krämer, M. Nardone, B. Fauqué, and K. Behnia, Phys. Rev. Lett. 112, 207002 (2014).
  * Kittaka  _et al._ (2014) S. Kittaka, Y. Aoki, Y. Shimura, T. Sakakibara, S. Seiro, C. Geibel, F. Steglich, H. Ikeda, and K. Machida, Phys. Rev. Lett. 112, 067002 (2014).
  * Nomoto  _et al._ (2016) T. Nomoto, K. Hattori, and H. Ikeda, Phys. Rev. B. 94, 174513 (2016).
  * Nomoto and Ikeda (2016b) T. Nomoto and H. Ikeda, Phys. Rev. Lett. 117, 217002 (2016b).
  * Seo  _et al._ (2017) K. Seo, J. D. Sau, and S. Tewari, Phys. Rev. B. 95, 205107 (2017).
  * Zhao  _et al._ (2018) Y. Zhao, S. Zeng, C. Lian, Z. Dai, S. Meng, and J. Ni, Phys. Rev. B. 98, 134514 (2018).
  * Singh  _et al._ (2018) J. Singh, A. Jayaraj, D. Srivastava, S. Gayen, A. Thamizhavel, and Y. Singh, Phys. Rev. B. 97, 054506 (2018).
  * Menke  _et al._ (2019) H. Menke, C. Timm, and P. M. R. Brydon, Phys. Rev. B. 100, 224505 (2019).
  * Kanasugi and Yanase (2019) S. Kanasugi and Y. Yanase, Phys. Rev. B. 100, 094504 (2019).
  * Xu  _et al._ (2019) C. Q. Xu, B. Li, J. J. Feng, W. H. Jiao, Y. K. Li, S. W. Liu, Y. X. Zhou, R. Sankar, N. D. Zhigadlo, H. B. Wang, Z. D. Han, B. Qian, W. Ye, W. Zhou, T. Shiroka, P. K. Biswas, X. Xu, and Z. X. Shi, Phys. Rev. B. 100, 134503 (2019).
  * Lado and Sigrist (2019) J. L. Lado and M. Sigrist, Phys. Rev. Research 1, 033107 (2019).
  * Hu  _et al._ (2020) L.-H. Hu, P. D. Johnson, and C. Wu, Phys. Rev. Research 2, 022021 (2020).
  * Timm and Bhattacharya (2021) C. Timm and A. Bhattacharya, Phys. Rev. B. 104, 094529 (2021).
  * Shang  _et al._ (2021) T. Shang, W. Xie, J. Z. Zhao, Y. Chen, D. J. Gawryluk, M. Medarde, M. Shi, H. Q. Yuan, E. Pomjakushina, and T. Shiroka, Phys. Rev. B. 103, 184517 (2021).
  * Wang  _et al._ (2021a) Z. Wang, S. Zeng, Y. Zhao, X. Wang, and J. Ni, Phys. Rev. B. 104, 174519 (2021a).
  * Sevik  _et al._ (2022) C. Sevik, J. Bekaert, M. Petrov, and M. V. Milošević, Phys. Rev. Materials 6, 024803 (2022).
  * Banerjee  _et al._ (2018) A. Banerjee, A. Sundaresh, R. Ganesan, and P. S. A. Kumar, ACS Nano 12, 12665 (2018).
  * Li  _et al._ (2017) H. Li, T. Zhou, J. He, H.-W. Wang, H. Zhang, H.-C. Liu, Y. Yi, C. Wu, K. T. Law, H. He, and J. Wang, Phys. Rev. B. 96, 075107 (2017).
  * Luttinger and Kohn (1955) J. M. Luttinger and W. Kohn, Phys. Rev. 97, 869 (1955).
  * Dresselhaus (1955) G. Dresselhaus, Phys. Rev. 100, 580 (1955).
  * Brydon  _et al._ (2016) P. M. R. Brydon, L. Wang, M. Weinert, and D. F. Agterberg, Phys. Rev. Lett. 116, 177001 (2016).
  * Com (j) The O(3) symmetry-broken case in the normal state is analyzed in the Supplemental Material (Sup, )
  * Savary  _et al._ (2017) L. Savary, J. Ruhman, J. W. F. Venderbos, L. Fu, and P. A. Lee, Phys. Rev. B. 96, 214514 (2017).
  * Venderbos  _et al._ (2018) J. W. F. Venderbos, L. Savary, J. Ruhman, P. A. Lee, and L. Fu, Phys. Rev. X. 8, 011029 (2018).
  * Yu and Liu (2018) J. Yu and C.-X. Liu, Phys. Rev. B. 98, 104514 (2018).
  * (64) See Supplemental Material for additional information about the reasons for TPTs at FE, properties of pseudospin rotation symmetry, global symmetries and an alternative topological index in terms of parity symmetry, O(3) symmetry-broken case, derivation of helical Dirac surface states at FEs, the effective surface Hamiltonian, and the pairing matrices in cubic symmetry. The Supplemental Material includes Refs. (Agterberg  _et al._ , 2017; Schnyder  _et al._ , 2008; Chiu  _et al._ , 2016; Fang  _et al._ , 2015; Yang  _et al._ , 2016; Tamura  _et al._ , 2020; Kawakami  _et al._ , 2018; Savary  _et al._ , 2017; Venderbos  _et al._ , 2018; Yu and Liu, 2018; Bahari  _et al._ , 2022; Fu and Kane, 2007; Hughes  _et al._ , 2011; Timm and Bhattacharya, 2021; Brydon  _et al._ , 2016; Morimoto and Furusaki, 2013).
  * Kim  _et al._ (2018) H. Kim, K. Wang, Y. Nakajima, R. Hu, S. Ziemak, P. Syers, L. Wang, H. Hodovanets, J. D. Denlinger, P. M. R. Brydon, D. F. Agterberg, M. A. Tanatar, R. Prozorov, and J. Paglione, Sci. Adv. 4, eaao4513 (2018).
  * Ishihara  _et al._ (2021) K. Ishihara, T. Takenaka, Y. Miao, Y. Mizukami, K. Hashimoto, M. Yamashita, M. Konczykowski, R. Masuki, M. Hirayama, T. Nomoto, R. Arita, O. Pavlosiuk, P. Wiśniewski, D. Kaczorowski, and T. Shibauchi, Phys. Rev. X. 11, 041048 (2021).
  * Com (c) To evaluate the tight-binding calculations, we used standard lattice regularizations $k_{z}^{2}\rightarrow 2(1-\text{cos}(k_{z}a))/a^{2}$ and $k_{z}\rightarrow\text{sin}(k_{z}a)/a$ where $a$ is the lattice constant.
  * Com (b) The IPR is defined by $I_{E}=\text{ln}(\sum_{i=1}^{8N}|\langle i|\psi\rangle|^{4}/\text{ln}(8N))$ where $N$ is the total number of sites and $|\psi\rangle$ indicates the eigenvector associated to $E$. The localized (extended) boundary (bulk) states have $I_{E}\rightarrow 0\ (-1)$.
  * Com (d) Such coexistence establishes topological bound states in the continuum (BIC). This has been addressed in different quantum mechanical systems (Plotnik  _et al._ , 2011; Mur-Petit and Molina, 2014; Baum  _et al._ , 2015; Hsu  _et al._ , 2016; Ricco  _et al._ , 2016; Xiao  _et al._ , 2017; Mukherjee  _et al._ , 2018; Bahari and Hosseini, 2019; Takeichi and Murakami, 2019; Chen  _et al._ , 2019; Jangjan and Hosseini, 2020; Benalcazar and Cerjan, 2020; Cerjan  _et al._ , 2020; Gong  _et al._ , 2021). However, the emergent topological BIC induced intrinsically by FE superconducting correlations has not been reported so far.
  * Suh  _et al._ (2023) H. G. Suh, Y. Yu, T. Shishidou, M. Weinert, P. M. R. Brydon, and D. F. Agterberg, arxiv.2301.11458 (2023).
  * Com (j) Our table is slightly different from the results given in Ref. (Savary  _et al._ , 2017). The details of our calculations are given explicitly in the SM (Sup, ).
  * Com (f) Note that $\hat{h}_{{\bf k},\lambda}^{-+}$ can be easily derived by substitution of $+\leftrightarrow-$ in Eq. (7). Note that $\delta_{{\bf k},\lambda}^{+-}$ fulfills the symmetry relation $\delta_{{\bf k},\lambda}^{-+}=-\delta_{{\bf-k},\lambda}^{+-}$.
  * Com (g) The BdG Hamiltonian preserves time-reversal $\hat{T}$, particle-hole $\hat{\mathscr{P}}$, chiral $\hat{\mathcal{C}}$, and inversion $\hat{\mathcal{P}}_{u}$ symmetries. The symmetry operators fulfill the properties $\hat{T}^{2}=-1$ and $\hat{\mathscr{P}}^{2}=\hat{\mathcal{P}}^{2}_{u}=1$.
  * Fu and Kane (2007) L. Fu and C. L. Kane, Phys. Rev. B 76, 045302 (2007).
  * Schnyder  _et al._ (2008) A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Phys. Rev. B. 78, 195125 (2008).
  * Hughes  _et al._ (2011) T. L. Hughes, E. Prodan, and B. A. Bernevig, Phys. Rev. B 83, 245132 (2011).
  * Morimoto and Furusaki (2013) T. Morimoto and A. Furusaki, Phys. Rev. B 88, 125129 (2013).
  * Chiu  _et al._ (2016) C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Rev. Mod. Phys. 88, 035005 (2016).
  * Berry (1984) M. V. Berry, Proc. R. Soc. Lond. A. Math. Phys. Sci. 392, 45 (1984).
  * Zak (1989) J. Zak, Phys. Rev. Lett. 62, 2747 (1989).
  * Resta (1994) R. Resta, Rev. Mod. Phys. 66, 899 (1994).
  * Xiao  _et al._ (2010) D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).
  * Com (h) The FE topological Dirac surface states are not restricted to the $A_{2u}$ channel. Rather, they also exist for the other pairing channels summarized in Table 1. The TPTs at FEs happen along the $C_{2}$ axes for the $E_{u}$ and $T_{2u}$ channels, and along the $C_{2}^{\prime}$ axes for the $T_{2u}$ channel.
  * Com (i) The Dirac points are two-fold degenerate at one surface at positive and at negative excitation energies due to the conservation of pseudospin rotation symmetry. Considering two surfaces, we obtain the same set of Dirac points on each surface. Thus, the degeneracy is doubled as compared to the case of single surface.
  * Com (i) Although, spin-orbit coupling is an essential ingredient for the formation of such surface states in our model, the effective masses in systems with nearly vanishing spin-orbit coupling can stem from different orbitals giving rise to non-vanishing values for $m-m^{\prime}$.
  * Kitaev (2001) A. Y. Kitaev, Phys-usp. 44, 131 (2001).
  * Oudah  _et al._ (2016) M. Oudah, A. Ikeda, J. N. Hausmann, S. Yonezawa, T. Fukumoto, S. Kobayashi, M. Sato, and Y. Maeno, Nat. Commun. 7, 13617 (2016).
  * Nakajima  _et al._ (2015) Y. Nakajima, R. Hu, K. Kirshenbaum, A. Hughes, P. Syers, X. Wang, K. Wang, R. Wang, S. R. Saha, D. Pratt, J. W. Lynn, and J. Paglione, Sci. Adv. 1, e1500242 (2015).
  * Com (a) The inversion breaking anti-symmetric spin-orbit-coupling (ASOC) in half-Heusler material is weak compared to symmetric spin-orbit coupling.
  * Metz  _et al._ (2019) T. Metz, S. Bae, S. Ran, I.-L. Liu, Y. S. Eo, W. T. Fuhrman, D. F. Agterberg, S. M. Anlage, N. P. Butch, and J. Paglione, Phys. Rev. B 100, 220504 (2019).
  * Xi  _et al._ (2016) X. Xi, Z. Wang, W. Zhao, J.-H. Park, K. T. Law, H. Berger, L. Forró, J. Shan, and K. F. Mak, Nature Physics 12, 139 (2016).
  * Zhang  _et al._ (2009) H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang, Nature Physics 5, 438 (2009).
  * Grundmann (2006) M. Grundmann, _The Physics of Semiconductors: An Introduction Including Devices and Nanophysics_ (Springer, Berlin Heidelberg, 2006).
  * Hsieh  _et al._ (2009) D. Hsieh, Y. Xia, L. Wray, D. Qian, A. Pal, J. H. Dil, J. Osterwalder, F. Meier, G. Bihlmayer, C. L. Kane, Y. S. Hor, R. J. Cava, and M. Z. Hasan, Science 323, 919 (2009).
  * Zhang  _et al._ (2011) Y. Zhang, L. X. Yang, M. Xu, Z. R. Ye, F. Chen, C. He, H. C. Xu, J. Jiang, B. P. Xie, J. J. Ying, X. F. Wang, X. H. Chen, J. P. Hu, M. Matsunami, S. Kimura, and D. L. Feng, Nat. Mater. 10, 273 (2011).
  * Evtushinsky  _et al._ (2014) D. V. Evtushinsky, V. B. Zabolotnyy, T. K. Kim, A. A. Kordyuk, A. N. Yaresko, J. Maletz, S. Aswartham, S. Wurmehl, A. V. Boris, D. L. Sun, C. T. Lin, B. Shen, H. H. Wen, A. Varykhalov, R. Follath, B. Büchner, and S. V. Borisenko, Phys. Rev. B. 89, 064514 (2014).
  * Sie  _et al._ (2019) E. J. Sie, T. Rohwer, C. Lee, and N. Gedik, Nat. Commun. 10, 3535 (2019).
  * Thompson  _et al._ (2020) J. J. P. Thompson, D. Pei, H. Peng, H. Wang, N. Channa, H. L. Peng, A. Barinov, N. B. M. Schröter, Y. Chen, and M. Mucha-Kruczyński, Nat. Commun. 11, 3582 (2020).
  * Liu  _et al._ (2020) W. Liu, L. Cao, S. Zhu, L. Kong, G. Wang, M. Papaj, P. Zhang, Y.-B. Liu, H. Chen, G. Li, F. Yang, T. Kondo, S. Du, G.-H. Cao, S. Shin, L. Fu, Z. Yin, H.-J. Gao, and H. Ding, Nat. Commun. 11, 5688 (2020).
  * Liu  _et al._ (2021) C. Liu, R. P. Day, F. Li, R. L. Roemer, S. Zhdanovich, S. Gorovikov, T. M. Pedersen, J. Jiang, S. Lee, M. Schneider, D. Wong, P. Dosanjh, F. J. Walker, C. H. Ahn, G. Levy, A. Damascelli, G. A. Sawatzky, and K. Zou, Nat. Commun. 12, 4573 (2021).
  * Lin  _et al._ (2021) C. Lin, M. Ochi, R. Noguchi, K. Kuroda, M. Sakoda, A. Nomura, M. Tsubota, P. Zhang, C. Bareille, K. Kurokawa, Y. Arai, K. Kawaguchi, H. Tanaka, K. Yaji, A. Harasawa, M. Hashimoto, D. Lu, S. Shin, R. Arita, S. Tanda, and T. Kondo, Nat. Mater. 20, 1093 (2021).
  * Chen  _et al._ (2022) S.-D. Chen, M. Hashimoto, Y. He, D. Song, J.-F. He, Y.-F. Li, S. Ishida, H. Eisaki, J. Zaanen, T. P. Devereaux, D.-H. Lee, D.-H. Lu, and Z.-X. Shen, Nature. 601, 562 (2022).
  * Bode (2003) M. Bode, Rep. Prog. Phys 66, 523 (2003).
  * Wiebe  _et al._ (2004) J. Wiebe, A. Wachowiak, F. Meier, D. Haude, T. Foster, M. Morgenstern, and R. Wiesendanger, Rev. Sci. Instrum. 75, 4871 (2004).
  * Alpichshev  _et al._ (2010) Z. Alpichshev, J. G. Analytis, J.-H. Chu, I. R. Fisher, Y. L. Chen, Z. X. Shen, A. Fang, and A. Kapitulnik, Phys. Rev. Lett. 104, 016401 (2010).
  * Fischer  _et al._ (2007) O. Fischer, M. Kugler, I. Maggio-Aprile, C. Berthod, and C. Renner, Rev. Mod. Phys. 79, 353 (2007).
  * Oka  _et al._ (2014) H. Oka, O. O. Brovko, M. Corbetta, V. S. Stepanyuk, D. Sander, and J. Kirschner, Rev. Mod. Phys. 86, 1127 (2014).
  * Neupane  _et al._ (2016) M. Neupane, N. Alidoust, M. M. Hosen, J.-X. Zhu, K. Dimitri, S.-Y. Xu, N. Dhakal, R. Sankar, I. Belopolski, D. S. Sanchez, T.-R. Chang, H.-T. Jeng, K. Miyamoto, T. Okuda, H. Lin, A. Bansil, D. Kaczorowski, F. Chou, M. Z. Hasan, and T. Durakiewicz, Nat. Commun. 7, 13315 (2016).
  * Cornils  _et al._ (2017) L. Cornils, A. Kamlapure, L. Zhou, S. Pradhan, A. A. Khajetoorians, J. Fransson, J. Wiebe, and R. Wiesendanger, Phys. Rev. Lett. 119, 197002 (2017).
  * Choi  _et al._ (2017) D.-J. Choi, C. Rubio-Verdú, J. de Bruijckere, M. M. Ugeda, N. Lorente, and J. I. Pascual, Nat. Commun. 8, 15175 (2017).
  * Fan  _et al._ (2021) P. Fan, F. Yang, G. Qian, H. Chen, Y.-Y. Zhang, G. Li, Z. Huang, Y. Xing, L. Kong, W. Liu, K. Jiang, C. Shen, S. Du, J. Schneeloch, R. Zhong, G. Gu, Z. Wang, H. Ding, and H.-J. Gao, Nat. Commun. 12, 1348 (2021).
  * Wang  _et al._ (2021b) D. Wang, J. Wiebe, R. Zhong, G. Gu, and R. Wiesendanger, Phys. Rev. Lett. 126, 076802 (2021b).
  * Schneider  _et al._ (2021) L. Schneider, P. Beck, J. Wiebe, and R. Wiesendanger, Sci. Adv. 7, eabd7302 (2021).
  * Plotnik  _et al._ (2011) Y. Plotnik, O. Peleg, F. Dreisow, M. Heinrich, S. Nolte, A. Szameit, and M. Segev, Phys. Rev. Lett. 107, 183901 (2011).
  * Mur-Petit and Molina (2014) J. Mur-Petit and R. A. Molina, Phys. Rev. B. 90, 035434 (2014).
  * Baum  _et al._ (2015) Y. Baum, T. Posske, I. C. Fulga, B. Trauzettel, and A. Stern, Phys. Rev. Lett. 114, 136801 (2015).
  * Hsu  _et al._ (2016) C. W. Hsu, B. Zhen, A. D. Stone, J. D. Joannopoulos, and M. Soljavcic, Nat. Rev. Mater. 1, 16048 (2016).
  * Ricco  _et al._ (2016) L. S. Ricco, Y. Marques, F. A. Dessotti, R. S. Machado, M. de Souza, and A. C. Seridonio, Phys. Rev. B. 93, 165116 (2016).
  * Xiao  _et al._ (2017) Y.-X. Xiao, G. Ma, Z.-Q. Zhang, and C. T. Chan, Phys. Rev. Lett. 118, 166803 (2017).
  * Mukherjee  _et al._ (2018) S. Mukherjee, J. Gomis-Bresco, P. Pujol-Closa, D. Artigas, and L. Torner, Phys. Rev. A. 98, 063826 (2018).
  * Bahari and Hosseini (2019) M. Bahari and M. V. Hosseini, Phys. Rev. B. 99, 155128 (2019).
  * Takeichi and Murakami (2019) M. Takeichi and S. Murakami, Phys. Rev. B. 99, 035128 (2019).
  * Chen  _et al._ (2019) Z.-G. Chen, C. Xu, R. Al Jahdali, J. Mei, and Y. Wu, Phys. Rev. B. 100, 075120 (2019).
  * Jangjan and Hosseini (2020) M. Jangjan and M. V. Hosseini, Sci. Rep. 10, 14256 (2020).
  * Benalcazar and Cerjan (2020) W. A. Benalcazar and A. Cerjan, Phys. Rev. B. 101, 161116 (2020).
  * Cerjan  _et al._ (2020) A. Cerjan, M. Jürgensen, W. A. Benalcazar, S. Mukherjee, and M. C. Rechtsman, Phys. Rev. Lett. 125, 213901 (2020).
  * Gong  _et al._ (2021) Z. Gong, J. Serafini, F. Yang, S. Preble, and J. Yao, Phys. Rev. Applied 16, 024059 (2021).

Supplemental Material

In this section, we address why TPTs at FEs occur only in odd-parity multiband
superconductors with time-reversal symmetry. We show that the FE pairing
potential is odd under the parity operation. To do so, we start by examining
the parity and time-reversal symmetry properties of the normal state and the
pairing potential. A fermionic state $|{\bf k},m_{j}\rangle$, with ${\bf k}$
as the 3D momentum, and $m_{j}$ being the magnetic quantum number of total
angular momentum $j$, transforms under inversion $P$ and time-reversal $T$
operations, respectively, as

$\displaystyle P|{\bf k},m_{j}\rangle$ $\displaystyle=$ $\displaystyle|-{\bf
k},m_{j}\rangle,$ (11) $\displaystyle T|{\bf k},m_{j}\rangle$ $\displaystyle=$
$\displaystyle(-1)^{j+m_{j}}|-{\bf k},-m_{j}\rangle.$ (12)

Inversion symmetry acts only on the momenta and the magnetic quantum number
remains intact. In this case, the matrix form of the inversion operator is a
$4\times 4$ identity matrix. The matrix representation of the anti-unitary
time-reversal operator takes the form $\hat{T}=\hat{\mathcal{R}}\mathcal{K}$
with $\mathcal{K}$ being the complex conjugation and the unitary part
$\hat{\mathcal{R}}=e^{i\pi\hat{J}_{y}}=i\hat{\sigma}_{x}\otimes\hat{\sigma}_{y}$.

The normal state preserves inversion and time-reversal symmetries described by

$\displaystyle\hat{P}\hat{\mathcal{H}}_{{\bf k}}\hat{P}^{-1}$
$\displaystyle=\hat{\mathcal{H}}_{-{\bf k}},\ \ \
\hat{\mathcal{R}}\hat{\mathcal{H}}_{{\bf
k}}^{*}\hat{\mathcal{R}}^{-1}=\hat{\mathcal{H}}_{-{\bf k}}.$ (13)

Moreover, the normal state preserves the combination of these symmetries given
by $(\hat{P}\hat{T})\hat{\mathcal{H}}_{{\bf
k}}(\hat{P}\hat{T})^{-1}=\hat{\mathcal{H}}_{{\bf k}}$ leading to doubly
degenerate eigenstates.

The parity of a pairing potential is distinguished by the orbital angular
momentum of Cooper pairing, i.e., $L=2n+1$ $(L=2n)$ denotes odd- (even-)
parity angular momenta with $n$ being non-negative integers. In this case, the
pairing potential is odd $u$ (even $g$) under the parity operation given by

$\hat{P}\hat{\Delta}_{{\bf k}}\hat{P}^{-1}=-(+)\hat{\Delta}_{-{\bf k}}.$ (14)

We can show that the BdG Hamiltonian, given in Eq. (1) of the main paper,
satisfies inversion symmetry for odd- (even-) parity pairing potential given
by

$\hat{\mathcal{P}}_{u,g}\hat{H}({\bf
k})\hat{\mathcal{P}}_{u,g}^{-1}=\hat{H}(-{\bf k}),$ (15)

where $\hat{\mathcal{P}}_{u}=\hat{\sigma}_{z}\otimes\hat{P}$ and
$\hat{\mathcal{P}}_{g}=\hat{\sigma}_{0}\otimes\hat{P}$. Furthermore, a pairing
potential is time-reversal symmetric described by

$\displaystyle\hat{\mathcal{R}}\hat{\Delta}_{{\bf
k}}^{*}\hat{\mathcal{R}}^{-1}$ $\displaystyle=\hat{\Delta}_{-{\bf
k}}=\pm\hat{\Delta}_{{\bf k}},$ (16)

where $-(+)$ corresponds to odd- (even-) parity pairing potentials. In our
investigations, we focus on odd-parity pairing potentials which preserve time-
reversal symmetry.

Due to inversion symmetry in the normal state, the pseudospin operator
$\hat{V}_{{\bf k}}=\\{\hat{V}_{{\bf k}}^{+},\hat{V}_{{\bf k}}^{-}\\}$ is an
even-parity matrix function, i.e., $\hat{V}_{{\bf k}}^{\pm}=\hat{V}_{-{\bf
k}}^{\pm}$. It is given explicitly by

$\displaystyle\hat{V}_{{\bf k}}^{+}$ $\displaystyle=\frac{1}{2}\frac{|{\bf
k}_{||}|}{|\mathbf{k}|}\left(\begin{array}[]{cc}2k_{z}k_{-}/k_{+}^{2}&k_{-}/k_{+}\\\
\sqrt{3}k_{-}/k_{+}&0\\\ 0&\sqrt{3}\\\ 1&-2k_{z}/k_{-}\end{array}\right),$
(21)

and

$\displaystyle\hat{V}_{{\bf k}}^{-}$ $\displaystyle=\frac{1}{\Gamma_{{\bf
k}}^{-}}\left(\begin{array}[]{cc}2\sqrt{3}k_{z}k_{-}{\bf
k}_{||}^{2}/k_{+}^{2}&-\sqrt{3}k_{-}^{2}\\\ -(|\mathbf{k}|^{2}+3k_{z}^{2}){\bf
k}_{||}^{2}/k_{+}^{2}&0\\\ 0&|\mathbf{k}|^{2}+3k_{z}^{2}\\\ \sqrt{3}{\bf
k}_{||}^{2}&2\sqrt{3}k_{z}k_{+}\end{array}\right),$ (26)

where ${\bf k}_{||}^{2}=k_{x}^{2}+k_{y}^{2}$, $\Gamma_{{\bf
k}}^{-}=2|\mathbf{k}|\sqrt{|\mathbf{k}|^{2}+3k_{z}^{2}}$,
$|\mathbf{k}|=\sqrt{{\bf k}_{||}^{2}+k_{z}^{2}}$, and $k_{\pm}=k_{x}\pm
ik_{y}$.

The parity of the pairing potential projected onto the pseudospin basis, i.e.,
$\hat{\Delta}_{{\bf k}}^{\nu\nu^{\prime}}$ with $\nu\nu^{\prime}\in\\{+,-\\}$,
is the same as the parity of the unprojected pairing potential. We prove this
by the symmetry relation

$\displaystyle\hat{\mathfrak{P}}\hat{\Delta}_{{\bf
k}}^{\nu\nu^{\prime}}\hat{\mathfrak{P}}^{-1}$
$\displaystyle=\hat{\mathfrak{P}}\hat{V}_{{\bf
k}}^{\nu\dagger}\hat{\Delta}_{{\bf k}}(\hat{V}_{-{\bf
k}}^{\nu^{\prime}\dagger})^{T}\hat{\mathfrak{P}}^{-1}$
$\displaystyle=\hat{\mathfrak{P}}\hat{V}_{{\bf
k}}^{\nu\dagger}(\pm\hat{P}^{-1}\hat{\Delta}_{-{\bf k}}\hat{P})(\hat{V}_{-{\bf
k}}^{\nu^{\prime}\dagger})^{T}\hat{\mathfrak{P}}^{-1}$
$\displaystyle=\pm\hat{V}_{-{\bf k}}^{\nu\dagger}\hat{\Delta}_{-{\bf
k}}(\hat{V}_{+{\bf k}}^{\nu^{\prime}\dagger})^{T}$
$\displaystyle=\pm\hat{V}_{{\bf k}}^{\nu\dagger}\hat{\Delta}_{-{\bf
k}}(\hat{V}_{-{\bf k}}^{\nu^{\prime}\dagger})^{T}$
$\displaystyle=\pm\hat{\Delta}_{-{\bf k}}^{\nu\nu^{\prime}},$ (27)

where $\hat{\mathfrak{P}}$ is the inversion operator defined by a $2\times 2$
identity matrix in pseudospin basis. Note that to derive Eq. (27), we consider
the even-parity property of the eigenspinors. This shows that the inversion
symmetry of the normal state is crucial to the determine the parity of the
projected pairing. Therefore, the FE pairing potential, which is responsible
for TPTs at FEs, has odd- (even-) parity if the orbital angular momentum of
Cooper pairing is odd (even), i.e.,

$\hat{\mathfrak{P}}\hat{\Delta}_{{\bf
k}}^{+-}\hat{\mathfrak{P}}^{-1}=-(+)\hat{\Delta}_{-{\bf k}}^{+-}.$ (28)

In the lattice representation, we replace even (odd) momenta in the finite
energy pairing potential $\hat{\Delta}_{k_{\nu}}^{+-}$ by

$\displaystyle\\!k_{\nu}^{2n+1}\\!\rightarrow\\![2(1-\text{cos}(k_{\nu}))]^{n}\text{sin}(k_{\nu}),\\!$
$\displaystyle\\!\\!\text{for }L$ $\displaystyle=\text{odd},\\!\\!$ (29)
$\displaystyle\\!k_{\nu}^{2n}\\!\rightarrow\\![2(1-\text{cos}(k_{\nu}))]^{n},\\!$
$\displaystyle\\!\\!\text{for }L$ $\displaystyle=\text{even},\\!\\!$ (30)

where $\nu\in\\{x,y,z\\}$ and we assume the lattice constant to be unity. The
odd-parity pairing channels contain odd powers of momenta in
$\hat{\Delta}_{{\bf k}}^{+-}$ along the TPTs direction. Thus, the
$\text{sin}(k_{\nu})$ term in Eq. (29) forces
$\hat{\Delta}_{k_{\nu}=\pm\pi}^{+-}=0$ at the parity-time-reversal-invariant
momenta (PTRIM), i. e., $k_{\nu}=\pm\pi$. This allows to close and reopen the
gap, through the manipulation of the normal state parameters at FEs.
Therefore, helical Dirac surface states emerge only for the odd-parity pairing
channels.

In contrast, bulk bands for even-parity pairing channels at FEs never close at
$k_{\nu}=\pm\pi$ , i.e., $\hat{\Delta}_{k_{\nu}=\pm\pi}^{+-}=4^{n}\neq 0$.
This prohibits the quantization of the topological index. Therefore, the
system remains in the topologically trivial phase.

We illustrate the aforementioned general argument by two examples. Since we
focus on the $p$-wave pairing channel, the momenta appear in linear order,
i.e., $L=1$ ($n=0$). Consider the $A_{2u}$ pairing along the direction ${\bf
k}=(0,0,k_{z})$. The BdG Hamiltonian is block diagonal through the unitary
transformation constructed from $\mathcal{\hat{D}}_{\mathbf{n}}(\pi,\epsilon)$
symmetry in the interband basis, i.e.,
$\hat{\mathscr{H}}^{\prime\prime}(0,0,k_{z})=\text{diag}(\hat{h}_{k_{z},+}^{-+},\hat{h}_{k_{z},+}^{+-},\hat{h}_{k_{z},-}^{+-},\hat{h}_{k_{z},-}^{-+})$.
The explicit matrix form of $\hat{h}_{k_{z},\pm}^{+-}$ in the lattice
representation is

$\hat{h}_{k_{z},\pm}^{+-}=\left(\begin{array}[]{cc}-E_{k_{z}}^{-}&\Delta\text{sin}(k_{z})\\\
\Delta\text{sin}(k_{z})&E_{k_{z}}^{+}\end{array}\right),$ (31)

where the off-diagonal terms correspond to the FE pairing. In Eq. (31),
$-E_{k_{z}}^{-}$ ($E_{k_{z}}^{+}$) is the hole (electron) band given by

$\displaystyle E_{k_{z}}^{-}$
$\displaystyle=2(\alpha+\beta/4)(1-\text{cos}(k_{z}))-\mu,$ (32)
$\displaystyle E_{k_{z}}^{+}$
$\displaystyle=2(\alpha+9\beta/4)(1-\text{cos}(k_{z}))-\mu.$ (33)

The TPT occurs at the PTRIM ${\bf k}=(0,0,\pm\pi)$ where the spectrum of Eq.
(31) becomes

$\displaystyle\varepsilon_{1}$ $\displaystyle=\mu-(4\alpha+\beta),\ \
\varepsilon_{2}=4\alpha+9\beta-\mu.$ (34)

We assume the normal state energy bands to have identical sign of curvature,
i.e., $\text{sgn}(\alpha)=\text{sgn}(\beta)=\pm 1$. In this case, the system
exhibits a gap at ${\bf k}=(0,0,\pm\pi)$, i.e.,
$|\varepsilon_{1}-\varepsilon_{2}|=2|4a+5b-\mu|>0$. The TPT occurs when
$|\varepsilon_{1}-\varepsilon_{2}|=0$. This results in the phase transition
relation $\mu=4\alpha+5\beta$. Hence, the TPT at nonzero excitation energies
can occur through the interplay between odd-parity superconductivity and
normal state crossings at finite excitation energies.

In realistic materials, $\alpha$ and $\beta$ are fixed and the system can
exhibit a topological phase for a certain range of chemical potentials. This
is illustrated by the topological phase diagram in Fig. 1(c) in the main text.

In contrast, the FE pairing never vanishes at PTRIM for even-parity pairing
channels. To illustrate this, consider the even-parity s-wave $E_{g}$ pairing
given by the pairing matrix $\hat{\Delta}_{{\bf
k}}=(\Delta/3)(2\hat{J}_{z}^{2}-\hat{J}_{x}^{2}-\hat{J}_{y}^{2})\hat{\mathcal{R}}$
(Agterberg  _et al._ , 2017). This pairing matrix exhibits vanishing
(nonvanishing) intra- (inter-) band pairing at momenta
$2k_{z}^{2}=k_{x}^{2}+k_{y}^{2}$, i.e., $\hat{\Delta}_{{\bf
k}}^{\text{intra}}=0$ and $\hat{\Delta}_{{\bf k}}^{+-}\neq 0$. Therefore, this
instability channel is a candidate for the TPTs at FEs. The FE pairing is even
under the parity exchange, i.e., $\hat{\mathfrak{P}}\hat{\Delta}_{{\bf
k}}^{+-}\hat{\mathfrak{P}}^{-1}=\hat{\Delta}_{-{\bf k}}^{+-}$ with
$\hat{\Delta}_{{\bf k}}^{+-}=\Delta(k_{+}/k_{-})\hat{\tau}_{x}$ where
$k_{\pm}=k_{x}\pm ik_{y}$. In this case, the interband superconducting
Hamiltonian becomes

Figure 3: Helical topological Dirac cones at finite excitation energies
induced by FE Cooper pairing. A pair of helical Dirac cones (pink) emerge on a
surface and the degenerate partner appears on the other surface (green). Black
arrows illustrate the spin texture of the surface states exhibiting spin-
momentum locking due to the helical property. $\mathscr{E}_{\text{DP}}$
denotes energy of the Dirac points and $\mu$ is the Fermi energy of the normal
state.

$\hat{\mathsf{H}}_{{\bf k}}^{+-}=\left(\begin{array}[]{cc}E_{{\bf
k}}^{+}\hat{\tau}_{0}&\Delta(k_{+}/k_{-})\hat{\tau}_{x}\\\
\Delta(k_{-}/k_{+})\hat{\tau}_{x}&-E_{{\bf
k}}^{-}\hat{\tau}_{0}\end{array}\right),$ (35)

where

$\displaystyle E_{{\bf k}}^{+}$
$\displaystyle=(3/8)(4\alpha+9\beta)(k_{x}^{2}+k_{y}^{2})-\mu,$ (36)
$\displaystyle E_{{\bf k}}^{-}$
$\displaystyle=(3/8)(4\alpha+\beta)(k_{x}^{2}+k_{y}^{2})-\mu.$ (37)

Equation (35) preserves a rotational symmetry along the z-axis in pseudospin
basis with rotation angle $\theta=\pi$ given by
$\mathcal{\hat{D}}_{\mathbf{n}_{z}}(\pi,1)=\text{diag}(i\hat{\tau}_{z},-i\hat{\tau}_{z})$.
Therefore, $\hat{\mathsf{H}}_{{\bf k}}^{+-}$ becomes diagonal in the
eigenbasis of $\mathcal{\hat{D}}_{\mathbf{n}_{z}}(\pi,1)$ labeled by
eigenvalues $\lambda=\pm i$ obtained by
$\hat{\mathcal{U}}^{\dagger}\hat{\mathsf{H}}_{{\bf
k}}^{+-}\hat{\mathcal{U}}=\text{diag}(\hat{h}_{{\bf k},+i}^{+-},\hat{h}_{{\bf
k},-i}^{+-})$ with

$\hat{h}_{{\bf k},\pm i}^{+-}=\left(\begin{array}[]{cc}-E_{{\bf
k}}^{-}&\Delta(k_{-}/k_{+})\\\ \Delta(k_{+}/k_{-})&E_{{\bf
k}}^{+}\end{array}\right),$ (38)

where the matrix of eigenvectors for
$\mathcal{\hat{D}}_{\mathbf{n}_{z}}(\pi,1)$ is given by

$\hat{\mathcal{U}}=\left(\begin{array}[]{cccc}0&1&0&0\\\ 0&0&0&1\\\ 0&0&1&0\\\
1&0&0&0\end{array}\right).$ (39)

The spectrum of Eq. (38) is always gapped at FEs. Hence, TPTs at FEs are not
possible.

##  Appendix A Pseudospin rotation symmetry properties

In this section, we explain the constraint imposed by pseudospin rotation
symmetry on the FE pairing potential. The motivation behind such analysis is
to show the explicit form of a general FE pairing matrix fulfilling the
pseudospin rotation symmetry. To do so, we present the FE pairing matrix in a
general form

$\hat{\Delta}_{{\bf k}}^{+-}\equiv\boldsymbol{\mathfrak{g}}_{{\bf
k}}\cdot\hat{\boldsymbol{\tau}},$ (40)

where the four-component $\boldsymbol{\mathfrak{g}}_{{\bf k}}$ vector is given
by $\boldsymbol{\mathfrak{g}}_{{\bf
k}}=(\mathfrak{g}_{0},\mathfrak{g}_{x},\mathfrak{g}_{y},\mathfrak{g}_{z})$
with the components being momentum dependent functions (this dependency is
dropped). The four-component vector of Pauli matrices is defined in the inter-
band basis by
$\hat{\boldsymbol{\tau}}=(\hat{\tau}_{0},\hat{\tau}_{x},\hat{\tau}_{y},\hat{\tau}_{z})$
with $\hat{\tau}_{0}$ being a $2\times 2$ identity matrix.

The matrix representation of a FE pairing potential converts to a pairing
function $\delta_{{\bf k},\lambda}^{+-}$ (see Eq. (5) in the main text) in the
presence of $\mathcal{\hat{D}}_{\mathbf{n}}(\pi,\epsilon)$ symmetry where the
rotation angle is about the arbitrary unit vector
$\mathbf{n}=(n_{x},n_{y},n_{z})$. In this case, the relation between
$\hat{\Delta}_{{\bf k}}^{+-}$and $\delta_{{\bf k},\lambda}^{+-}$ can be
distinguished by the components of the vector $\boldsymbol{\mathfrak{g}}_{{\bf
k}}$. Interestingly, the pseudospin rotation symmetry about the x(y)[z] axis,
defined by

$\mathbf{n}_{x}\equiv(1,0,0),\ \mathbf{n}_{y}\equiv(0,1,0),\
\mathbf{n}_{z}\equiv(0,0,1),$ (41)

enforces $\delta_{{\bf k},\lambda}^{+-}$ to select up to two components of the
vector $\boldsymbol{\mathfrak{g}}_{{\bf k}}$. After some algebra, such
constraints for the $\epsilon=-1$ representation can be obtained by

$\displaystyle\\!\\![\mathcal{\hat{D}}_{\mathbf{n}_{z}}(\pi,-1),\hat{\mathsf{H}}_{{\bf
k}}^{+-}]$ $\displaystyle\\!=\\!0\rightarrow\delta_{{\bf
k},\lambda}^{+-}=\mathfrak{g}_{0}+\lambda\mathfrak{g}_{z},$ (42)
$\displaystyle\\!\\![\mathcal{\hat{D}}_{\mathbf{n}_{y}}(\pi,-1),\hat{\mathsf{H}}_{{\bf
k}}^{+-}]$ $\displaystyle\\!=\\!0\rightarrow\delta_{{\bf
k},\lambda}^{+-}=\lambda i\mathfrak{g}_{x}-\mathfrak{g}_{z},$ (43)
$\displaystyle\\!\\![\mathcal{\hat{D}}_{\mathbf{n}_{x}}(\pi,-1),\hat{\mathsf{H}}_{{\bf
k}}^{+-}]$ $\displaystyle\\!=\\!0\rightarrow\delta_{{\bf
k},\lambda}^{+-}=\mathfrak{g}_{0}+\lambda\mathfrak{g}_{x}.$ (44)

To interpret the above relations, we consider Eq. (42) as an example. In this
case, $\delta_{{\bf k},\lambda}^{+-}$ preserves the pseudospin rotation
symmetry under the rotation angle $\theta=\pi$ about the $z$-axis if either
$\mathfrak{g}_{0}$ or $\mathfrak{g}_{z}$, or both of them are finite, and the
other components vanish, i.e., $\boldsymbol{\mathfrak{g}}_{{\bf
k}}=(\mathfrak{g}_{0},0,0,\mathfrak{g}_{z})$. This is a useful result since we
can directly ascertain the pseudospin rotation symmetry of a FE pairing
potential in multiband superconductors by checking only the components of the
vector $\boldsymbol{\mathfrak{g}}_{{\bf k}}$ .

The dot product of the vector $\boldsymbol{\mathfrak{g}}_{{\bf k}}$ and the
Pauli matrices in Eq. (40) allows us to find the relation between
$\hat{\Delta}_{{\bf k}}^{+-}$ and a proper
$\mathcal{\hat{D}}_{\mathbf{n}}(\pi,\epsilon)$ symmetry given by

$\displaystyle\\!\\!\mathcal{\hat{D}}_{\mathbf{n}_{z}}(\pi,-1):$
$\displaystyle\ \hat{\Delta}_{{\bf
k}}^{+-}=\mathfrak{g}_{0}\hat{\tau}_{0}+\mathfrak{g}_{z}\hat{\tau}_{z},$ (45)
$\displaystyle\\!\\!\mathcal{\hat{D}}_{\mathbf{n}_{y}}(\pi,-1):$
$\displaystyle\ \hat{\Delta}_{{\bf
k}}^{+-}=\mathfrak{g}_{x}\hat{\tau}_{x}+\mathfrak{g}_{z}\hat{\tau}_{z},$ (46)
$\displaystyle\\!\\!\mathcal{\hat{D}}_{\mathbf{n}_{x}}(\pi,-1):$
$\displaystyle\ \hat{\Delta}_{{\bf
k}}^{+-}=\mathfrak{g}_{0}\hat{\tau}_{0}+\mathfrak{g}_{x}\hat{\tau}_{x}.$ (47)

Importantly, we observe that the component $\mathfrak{g}_{y}$ in Eqs. (42-47)
is absent. However, the pseudospin rotation symmetry with $\epsilon=+1$
representation reveals the allowed symmetry form of $\delta_{{\bf
k},\lambda}^{+-}$ including the $y$-component of
$\boldsymbol{\mathfrak{g}}_{{\bf k}}$,

$\displaystyle\\!\\![\mathcal{\hat{D}}_{\mathbf{n}_{z}}(\pi,1),\hat{\mathsf{H}}_{{\bf
k}}^{+-}]$ $\displaystyle\\!=\\!0\rightarrow\delta_{{\bf
k},\lambda}^{+-}=\mathfrak{g}_{x}-i\lambda\mathfrak{g}_{y},$ (48)
$\displaystyle\\!\\![\mathcal{\hat{D}}_{\mathbf{n}_{y}}(\pi,1),\hat{\mathsf{H}}_{{\bf
k}}^{+-}]$ $\displaystyle\\!=\\!0\rightarrow\delta_{{\bf
k},\lambda}^{+-}=\mathfrak{g}_{0}+\lambda\mathfrak{g}_{y},$ (49)
$\displaystyle\\!\\![\mathcal{\hat{D}}_{\mathbf{n}_{x}}(\pi,1),\hat{\mathsf{H}}_{{\bf
k}}^{+-}]$ $\displaystyle\\!=\\!0\rightarrow\delta_{{\bf
k},\lambda}^{+-}=-i\lambda\mathfrak{g}_{y}-\mathfrak{g}_{z}.$ (50)

In this case, the $\mathcal{\hat{D}}_{\mathbf{n}}(\pi,\epsilon)$ symmetry
forces the FE pairing matrix to have the explicit form

$\displaystyle\\!\\!\mathcal{\hat{D}}_{\mathbf{n}_{z}}(\pi,1)$
$\displaystyle:\ \hat{\Delta}_{{\bf
k}}^{+-}=\mathfrak{g}_{x}\hat{\tau}_{x}+\mathfrak{g}_{y}\hat{\tau}_{y},$ (51)
$\displaystyle\\!\\!\mathcal{\hat{D}}_{\mathbf{n}_{y}}(\pi,1)$
$\displaystyle:\ \hat{\Delta}_{{\bf
k}}^{+-}=\mathfrak{g}_{0}\hat{\tau}_{0}+\mathfrak{g}_{y}\hat{\tau}_{y},$ (52)
$\displaystyle\\!\\!\mathcal{\hat{D}}_{\mathbf{n}_{x}}(\pi,1)$
$\displaystyle:\ \hat{\Delta}_{{\bf
k}}^{+-}=\mathfrak{g}_{y}\hat{\tau}_{y}+\mathfrak{g}_{z}\hat{\tau}_{z}.$ (53)

Taking into account Eqs. (45-47) and (51-53), we can easily understand the
explicit form of the $\mathcal{\hat{D}}_{\mathbf{n}}(\pi,\epsilon)$ symmetry
along the TPT direction by looking at the components of $\hat{\Delta}_{{\bf
k}}^{+-}$, e.g., the last column of Table. I in the main text. For instance,
one of the TPT directions for the $T_{2u}$ irrep is $\mathbf{k}\in(\pm k,0,\pm
k)$ and the FE pairing potential is given by $\hat{\Delta}_{k}^{+-}=\Delta
k(\pm\hat{\tau}_{0}/3-i\hat{\tau}_{y})$. According to Eq. (52),
$\hat{\Delta}_{k}^{+-}$ preserves $\mathcal{\hat{D}}_{\mathbf{n}_{y}}(\pi,1)$
symmetry with $\mathfrak{g}_{0}=\pm\Delta k/3$ and $\mathfrak{g}_{y}=\pm
i\Delta k$.

Note that when $\hat{\Delta}_{{\bf k}}^{+-}$ contains only one component of
$\boldsymbol{\mathfrak{g}}_{{\bf k}}$ along the TPT direction, the rotation
axis of the $\mathcal{\hat{D}}_{\mathbf{n}}(\pi,\epsilon)$ is not unique.
Importantly, the momentum dependency of components of the
$\boldsymbol{\mathfrak{g}}_{{\bf k}}$ along the TPT direction is proportional
to the orbital angular momentum of Cooper pairing. For instance,
$\boldsymbol{\mathfrak{g}}_{{\bf k}}$ is proportional to linear order of
momenta at TPT directions since we focus on $p$-wave pairing, i, e., $L=1$.

Figure 4:  Enlarged view of the spectra given in Fig. 2(c) of the main text.
The thickness is 180 layers and $E_{0}=\Delta$/2. The high-symmetry points are
$X=(0,\pi/a,0)$ and $\Gamma=(0,0,0)$ where $a$ is the lattice constant in
tight binding calculations. The surface (bulk) states are colorful (dark).
Blue arrows indicate the helical Dirac points induced by topological FE
pairing. The gray dashed line indicates the Fermi energy. (a) [(b)] Surface
states inside (outside) of the bulk states before (after) a Lifshitz
transition for $\mu=-5.2E_{0}$ ($\mu=-E_{0}$). (b) The helical surface states
resemble a butterfly shape, and establish dispersive Majorana modes away from
the high symmetry points at the Fermi energy. The other parameters are
$\alpha=-\Delta/2$ and $\beta=0.3\alpha$.

The reason for having two representations for the
$\mathcal{\hat{D}}_{\mathbf{n}}(\pi,\epsilon)$ symmetry can be seen by the
unitary transformation $\hat{\mathcal{V}}$ between the eigenbases for
$\epsilon=+1$ and $\epsilon=-1$,

$\hat{\phi}_{+1}=\hat{\mathcal{V}}\hat{\phi}_{-1},\ \ \
\hat{\mathcal{V}}=\hat{\mathcal{V}}^{-1}=\hat{\mathcal{V}}^{\dagger}=\left(\begin{array}[]{cccc}0&0&1&0\\\
0&1&0&0\\\ 1&0&0&0\\\ 0&0&0&1\end{array}\right),$ (54)

where $\hat{\phi}_{\epsilon=\pm 1}=\hat{\mathcal{U}}^{-1}\hat{\varphi}_{{\bf
k}}^{+-}$ with $\hat{\varphi}_{{\bf k}}^{+-}$being the $4\times 1$ column as
the basis of FE superconducting Hamiltonian $\hat{\mathsf{H}}_{{\bf k}}^{+-}$,
and $\hat{\mathcal{U}}$ is the eigenvector matrix for
$\mathcal{\hat{D}}_{\mathbf{n}}(\pi,\epsilon)$ operator.

##  Appendix B Global symmetries

The explicit matrix representations for time-reversal and particle-hole
symmetry operators in BdG formalism are given by
$\hat{\mathcal{T}}=i\hat{\varsigma}_{0}\otimes\hat{\sigma}_{x}\otimes\hat{\sigma}_{y}\mathcal{K}$
and
$\hat{\mathscr{P}}=\hat{\varsigma}_{x}\otimes\hat{\sigma}_{0}\otimes\hat{\sigma}_{0}\mathcal{K}$,
respectively. Here, $\mathcal{K}$ is the complex-conjugate operator,
$\hat{\varsigma}_{i}$ and $\hat{\sigma}_{i}$ with $i\in\\{0,x,y,z\\}$ are
Pauli matrices acting on particle-hole and spin subspaces, respectively. The
symmetry relation reads $\hat{\mathcal{O}}\hat{H}^{*}({\bf
k})\hat{\mathcal{O}}^{-1}=\epsilon\hat{H}(-{\bf k})$ with $\epsilon=+1(-1)$
for $\hat{\mathcal{O}}=\hat{\mathcal{T}}(\hat{\mathscr{P}})$. Additionally, we
can construct a chiral symmetry operator
$\hat{\mathcal{C}}=\hat{\mathcal{T}}\hat{\mathscr{P}}$ fulfilling
$\hat{\mathcal{C}}\hat{H}({\bf k})\hat{\mathcal{C}}^{-1}=-\hat{H}({\bf k})$
(Schnyder  _et al._ , 2008; Morimoto and Furusaki, 2013; Chiu  _et al._ ,
2016). The symmetry operators have the properties $\hat{\mathcal{T}}^{2}=-1$
and $\hat{\mathscr{P}}^{2}=1$. In addition, the inversion symmetry operator is
defined by $\hat{\mathcal{P}}_{u}=\hat{\sigma}_{z}\otimes\hat{P}$ fulfilling
$\hat{\mathcal{P}}_{u}^{2}=1$ .

We complete our discussion on topological phase transitions induced by FE
pairing by proposing an alternative approach to calculate the topological
index due to the presence of inversion symmetry (Fu and Kane, 2007; Hughes
_et al._ , 2011). The BdG Hamiltonian commutes with $\hat{\mathcal{P}}_{u}$
and $\hat{T}$ operators at eight PTRIM given by

$\displaystyle\\!\\!\mathbf{k}_{\text{p}}\\!$
$\displaystyle=\\!(0,\\!0,\\!0),\mathbf{k}_{\text{p}}\\!=\\!(0,\\!0,\\!\pi),\mathbf{k}_{\text{p}}\\!=\\!(0,\\!\pi,\\!0),\mathbf{k}_{\text{p}}\\!=\\!(\pi,\\!0,\\!0),\\!\\!$
$\displaystyle\\!\\!\mathbf{k}_{\text{p}}\\!$
$\displaystyle=\\!(0,\\!\pi,\\!\pi),\mathbf{k}_{\text{p}}\\!=\\!(\pi,\\!\pi,\\!0),\mathbf{k}_{\text{p}}\\!=\\!(\pi,\\!0,\\!\pi),\mathbf{k}_{\text{p}}\\!=\\!(\pi,\\!\pi,\\!\pi).\\!\\!$
(55)

Generally, the FE pairing arising form odd-parity pairing channels can induce
topological phase transition along the directions $\mathbf{K}$ connecting the
PTRIM. The criteria rely on the presence and absence of finite and low energy
pairings along such directions, respectively. In this case, the BdG
Hamiltonian becomes block diagonal

| $\hat{\mathcal{T}}$ | $\hat{\mathscr{P}}$ | $\hat{\mathcal{C}}$ | $\hat{\mathcal{P}}_{u}$
---|---|---|---|---
$\hat{\mathcal{T}}$ | Commute | Commute | Commute | Commute
$\hat{\mathscr{P}}$ |  | Commute | Commute | Anticommute
$\hat{\mathcal{C}}$ |  |  | Commute | Anticommute
$\hat{\mathcal{P}}_{u}$ |  |  |  | Commute
Table 2: Commutation/anticommutation relations between the discrete
symmetries.

$\hat{\mathscr{H}}^{\prime}(\mathbf{K})=\text{diag}\big{(}\hat{h}_{\mathbf{K},+}^{+-},\hat{h}_{\mathbf{K},+}^{+-},\hat{h}_{\mathbf{K},-}^{-+},\hat{h}_{\mathbf{K},-}^{-+}\big{)},$
(56)

where

$\hat{h}_{\mathbf{K},\lambda}^{+-}=\left(\begin{array}[]{cc}-E_{\mathbf{K}}^{-}&\big{(}\delta_{\mathbf{K},\lambda}^{+-}\big{)}^{*}\\\
\delta_{\mathbf{K},\lambda}^{+-}&E_{\mathbf{K}}^{+}\end{array}\right).$ (57)

In Eq. (56), $\hat{h}_{\mathbf{K},\lambda}^{-+}$ is the particle-hole partner
of $\hat{h}_{\mathbf{K},\lambda}^{+-}$ implied by
$\hat{\tau}_{y}(\hat{h}_{\mathbf{K},\lambda}^{+-})^{*}\hat{\tau}_{y}^{-1}=\hat{h}_{\mathbf{K},\lambda}^{-+}$.
Note that each block $\hat{h}_{\mathbf{k},\lambda}^{\pm\mp}$ preserves an
effective time-reversal symmetry $\hat{T}=\mathcal{K}$ with $\hat{T}^{2}=+1$.
Particle-hole symmetry and the conventional time-reversal symmetry with
$\hat{T}^{2}=-1$ are broken due to the different diagonal entries in the
symmetry blocks. Also, each block in Eq. (56) satisfies inversion symmetry
implied by
$\hat{\tau}_{z}\hat{h}_{\mathbf{K},\lambda}^{\nu\nu^{\prime}}\hat{\tau}_{z}^{-1}=\hat{h}_{-\mathbf{K},\lambda}^{\nu\nu^{\prime}}$
with $\nu,\nu^{\prime}\in\\{+,-\\}$. In this case,
$\hat{h}_{\mathbf{K},\lambda}^{\nu\nu^{\prime}}$ commutes with
$\hat{\tau}_{z}$, and the negative parity of the eigenstates associated to the
lower energy band in sector $\hat{h}_{\mathbf{K},\lambda}^{\nu\nu^{\prime}}$
determines the topological nature of the phase transition. Note that only one
of the directions connecting the PTRIM is sufficient to capture the nontrivial
topological phase induced by FE Cooper pairing. Therefore, to characterize the
topological phases, we define a topological index

$\mathcal{N}=|n_{\Gamma}-n_{\mathbf{p}}|,$ (58)

where $n_{\Gamma}$ and $n_{\mathbf{p}}$ are the number of negative eigenvalues
of the parity operator at the $\Gamma$ point
$\mathbf{k}_{\text{p}}=(0,\\!0,\\!0)$, and other PTRIM. Note that
$\mathcal{N}$ is a $\mathbb{Z}_{2}$ topological index taking two values, i.e.,
$\mathcal{N}=1(0)$ in the topologically nontrivial (trivial) phase. The full
topological index can be derived by the summation of the topological indices
associated to the decoupled blocks in Eq. (56) as

$Z=\sum_{\lambda=\pm}(\mathcal{N}_{\lambda}^{+-}+\mathcal{N}_{\lambda}^{-+}),$
(59)

where $\mathcal{N}_{\lambda}^{\nu\nu^{\prime}}$ corresponds to the block
$\hat{h}_{\mathbf{K},\lambda}^{\nu\nu^{\prime}}$. In the topologically
nontrivial phase, conservation of parity leads to the quantization of the
topological index and the emergence of surface states. Pseudospin rotation
symmetry ensures that the surface states come in pairs establishing FE helical
Dirac points. Hence, the appearance of a pair of helical surface states at
positive (negative) excitation energies is signaled by
$\sum_{\lambda=\pm}\mathcal{N}_{\lambda}^{+-}=2$
($\sum_{\lambda=\pm}\mathcal{N}_{\lambda}^{-+}=2$). Considering two surfaces,
we observe four Dirac surface states on each surface (two for positive
excitation energies and two for negative energies) due to parity-time-reversal
symmetry, cf. Fig. 3.

A Lifshitz transition can move the FE helical Dirac surface states, shown in
Fig. 4(a), to the low energies. In this case, the surface states establish
Majorana modes at the Fermi energy, see Fig. 4(b).

Note that the point group symmetry of the normal-state Hamiltonian for any
real system is lower than the O(3) symmetry of the Luttinger-Kohn Hamiltonian
discussed in the main text. The Luttinger-Kohn Hamiltonian serves as a
valuable approximation for energy bands close to high-symmetry points.
However, the symmetry of the superconducting (BdG) Hamiltonian is relevant for
our analysis. It belongs to the $O_{h}$ symmetry group, as the O(3) symmetry
of the normal state transforms into cubic point group symmetry due to the
pairing matrices derived through irreducible representations of the cubic
point group symmetry. Importantly, the pseudospin rotation symmetry mentioned
above is present for multiband system with parity-time-reversal symmetry. Its
presence is independent of the O(3) symmetry of the normal state Hamiltonian.
Nonetheless, we examine our predictions when the normal state is influenced by
O(3) symmetry-breaking terms. To do so, we incorporate the cubic spin-orbit
coupling term to the normal state Hamiltonian as

Figure 5: Spectra of a slab with (001) surfaces in the topologically
nontrivial regime. The thickness is 100 layers. The cubic spin-orbit coupling
term is chosen as $\gamma=\alpha$. Hence, O(3) symmetry is broken in the
normal state. In panels (a) and (c), $k^{\prime}$ defines the direction in
momentum space with $k_{x}=k_{y}$. The chemical potential is chosen as (a)
$\mu=-5.2E_{0}$, (b) $\mu=-5.1E_{0}$, and (c,d) $\mu=-E_{0}$. The other
parameters are $\beta=0.3\alpha$, $\alpha=-(3/2)\Delta$, and
$E_{0}=(2/3)\Delta$. The high-symmetry points are $X=(0,\pi/a,0)$ and
$M=(\pi/a,\pi/a,0)$.

$\\!\\!\mathcal{\hat{H}}_{{\bf k}}\\!=\alpha|{\bf
k}|^{2}\hat{I}_{4}\\!+\\!\beta\sum_{i}k_{i}^{2}\hat{J}_{i}^{2}\\!+\\!\gamma\sum_{i\neq
j}k_{i}k_{j}\hat{J}_{i}\hat{J}_{j}-\mu,\\!$ (60)

where $\gamma$ parametrizes the strength of the cubic spin-orbit coupling
term. This Hamiltonian refers to the O(3) symmetry-broken case in the normal
state if $\gamma\neq\beta$. We adopt the tight-binding regularization of Eq.
(60) in the BdG form. The spectral results within the topologically non-
trivial phase are illustrated in Fig. 5. Evidently, we can also observe the
emergence of helical Dirac surface states with particle-hole character away
form the Fermi energy due to the unconventional finite energy Cooper pairing
in the case with broken O(3) symmetry, see Fig. 5(a,b). In addition, we
investigate a Lifshitz transition around the $\Gamma$ point by a variation of
the chemical potential as illustrated in Fig. 5(c,d). In this case, the
surface states are shifted towards low energies. If the superconducting gap
possesses nodes at zero excitation energy, these surface states connect the
nodal points, see Fig. 5(c). If the superconducting gap is fully developed at
zero excitation energy, these surface states connect to the bulk states at
finite excitation energy, see Fig. 5(d).

##  Appendix C Helical Dirac surface states at finite energies

In this section, we provide details for the derivation of the helical Dirac
surface states described in the main text. The BdG Hamiltonian along the [001]
direction becomes

$\hat{H}(k_{z})=\left(\begin{array}[]{cccc}\hat{H}_{1}&0&\hat{H}_{3}&0\\\
0&\hat{H}_{2}&0&\hat{H}_{3}\\\ \hat{H}_{3}&0&-\hat{H}_{1}&0\\\
0&\hat{H}_{3}&0&-\hat{H}_{2}\end{array}\right),$ (61)

with

$\hat{H}_{1}=\text{diag}(E_{k_{z}}^{+},E_{k_{z}}^{-}),\
\hat{H}_{2}=\text{diag}(E_{k_{z}}^{-},E_{k_{z}}^{+}),$ (62)

where $E_{k_{z}}^{+}=(\alpha+9\beta/4)k_{z}^{2}-\mu$,
$E_{k_{z}}^{-}=(\alpha+\beta/4)k_{z}^{2}-\mu$,
$\hat{H}_{3}\\!\\!=\\!\\!\boldsymbol{\Delta}k_{z}\hat{\sigma}_{x}$ , and
$\boldsymbol{\Delta}=\sqrt{3}\Delta/2$. To be able to apply our theory,
presented in Eqs. (2-6) of the main text, we should represent $\hat{H}(k_{z})$
in the pseudospin basis where each diagonal block contains a pair of doubly
degenerate bands. This can be done through the unitary transformation
$\hat{\mathscr{U}}_{1}$

$\displaystyle\\!\hat{H}_{1}(k_{z})\\!$
$\displaystyle=\\!\hat{\mathscr{U}}_{1}^{\dagger}H(k_{z})\hat{\mathscr{U}}_{1}\\!$
(63)
$\displaystyle=\\!\\!\left(\\!\\!\begin{array}[]{cccc}E_{k_{z}}^{+}&\\!0&0&\\!\boldsymbol{\Delta}k_{z}\\\
0&\\!E_{k_{z}}^{-}&\boldsymbol{\Delta}k_{z}&\\!0\\\
0&\\!\boldsymbol{\Delta}k_{z}&-E_{k_{z}}^{+}&\\!0\\\
\boldsymbol{\Delta}k_{z}&\\!0&0&-E_{k_{z}}^{-}\end{array}\\!\\!\right)\\!\otimes\\!\hat{\sigma}_{0},\\!\\!\\!\\!\\!\\!\\!$
(68)

where the unitary matrix
$\hat{\mathscr{U}}_{1}^{\dagger}=\hat{\mathscr{U}}_{1}^{-1}$ is given by

$\hat{\mathscr{U}}_{1}=\left(\begin{array}[]{cccccccc}1&0&0&0&0&0&0&0\\\
0&0&1&0&0&0&0&0\\\ 0&0&0&1&0&0&0&0\\\ 0&1&0&0&0&0&0&0\\\ 0&0&0&0&1&0&0&0\\\
0&0&0&0&0&0&1&0\\\ 0&0&0&0&0&0&0&1\\\ 0&0&0&0&0&1&0&0\end{array}\right).$ (69)

In the next step, Eq. (68) should be represented in the inter-band basis. This
can be done through another unitary transformation $\hat{\mathscr{U}}_{2}$ as

$\displaystyle\hat{\mathscr{H}}(k_{z})$
$\displaystyle=\hat{\mathscr{U}}_{2}^{\dagger}\hat{H}_{1}(k_{z})\hat{\mathscr{U}}_{2}=(\hat{\mathscr{U}}_{1}\hat{\mathscr{U}}_{2})^{\dagger}H(k_{z})\hat{\mathscr{U}}_{1}\hat{\mathscr{U}}_{2}$
(70)
$\displaystyle=\left(\begin{array}[]{cc}\hat{\mathsf{H}}_{k_{z}}^{+-}&0\\\
0&\hat{\mathsf{H}}_{k_{z}}^{-+}\end{array}\right),$ (73)

where

$\displaystyle\\!\hat{\mathsf{H}}_{k_{z}}^{+-}\\!$
$\displaystyle=\\!\left(\\!\begin{array}[]{cc}E_{k_{z}}^{+}&\boldsymbol{\Delta}k_{z}\\\
\boldsymbol{\Delta}k_{z}&-E_{k_{z}}^{-}\end{array}\\!\right)\\!\otimes\\!\hat{\sigma}_{0},\\!\\!$
(76) $\displaystyle\\!\hat{\mathsf{H}}_{k_{z}}^{-+}\\!$
$\displaystyle=\\!\left(\\!\begin{array}[]{cc}E_{k_{z}}^{-}&\boldsymbol{\Delta}k_{z}\\\
\boldsymbol{\Delta}k_{z}&-E_{k_{z}}^{+}\end{array}\\!\right)\\!\otimes\\!\hat{\sigma}_{0}.\\!\\!$
(79)

The transformation $\hat{\mathscr{U}}_{2}$ fulfills the unitary property
$\hat{\mathscr{U}}_{2}^{\dagger}=\hat{\mathscr{U}}_{2}^{-1}$. Its explicit
matrix form is given by

$\hat{\mathscr{U}}_{2}=\left(\begin{array}[]{cccccccc}0&1&0&0&0&0&0&0\\\
1&0&0&0&0&0&0&0\\\ 0&0&0&0&0&1&0&0\\\ 0&0&0&0&1&0&0&0\\\ 0&0&0&0&0&0&0&1\\\
0&0&0&0&0&0&1&0\\\ 0&0&0&1&0&0&0&0\\\ 0&0&1&0&0&0&0&0\end{array}\right).$ (80)

Note that $\hat{\mathscr{H}}(k_{z})$, given in Eq. (73), is identical to Eq.
(2) of the main text. In this case, $\boldsymbol{\Delta}k_{z}\hat{\sigma}_{0}$
is the FE pairing potential in Eqs. (76) and Eqs. (79). The zeros on the off-
diagonal blocks of Eq. $\text{(\ref{FEblocks}})$ indicate vanishing intra-band
pairing giving rise to the presence of nodes at the Fermi energy. The diagonal
blocks in $\hat{\mathscr{H}}(k_{z})$ exhibit pseudospin-$\pi$ rotation
symmetry along the z(x)-axis, i.e., the symmetry relations given in Eqs. (45)
and (47). The explicit matrix form for such an operator can be obtained as

$\hat{\mathcal{D}}_{\mathbf{n}_{z}}(\pi,-1)=\left(\begin{array}[]{cccc}0&i&0&0\\\
i&0&0&0\\\ 0&0&0&i\\\ 0&0&i&0\end{array}\right).$

Thus, Eq. $(\text{\ref{FEblocks}})$ is further reducible in the eigenspace of
$\hat{\mathcal{D}}_{\mathbf{n}_{z}}(\pi,-1)$. This can be done through the
unitary matrix $\hat{\mathscr{U}}_{3}$

$\displaystyle\hat{\mathscr{H}}^{\prime}(k_{z})$
$\displaystyle=\hat{\mathscr{U}}_{3}^{\dagger}\hat{\mathscr{H}}(k_{z})\hat{\mathscr{U}}_{3}$
(81)
$\displaystyle=(\hat{\mathscr{U}}_{1}\hat{\mathscr{U}}_{2}\hat{\mathscr{U}}_{3})^{\dagger}H(k_{z})\hat{\mathscr{U}}_{1}\hat{\mathscr{U}}_{2}\hat{\mathscr{U}}_{3}$
(82)
$\displaystyle=\text{diag}(\hat{h}_{k_{z},+}^{+-},\hat{h}_{k_{z},+}^{+-},\hat{h}_{k_{z},+}^{-+},\hat{h}_{k_{z},-}^{-+}),$
(83)

with

$\\!\hat{h}_{k_{z},\pm}^{+-}\\!=\\!\left(\\!\begin{array}[]{cc}-E_{k_{z}}^{-}&\boldsymbol{\Delta}k_{z}\\\
\boldsymbol{\Delta}k_{z}&E_{k_{z}}^{+}\end{array}\\!\right)\\!,\
\hat{h}_{k_{z},\pm}^{-+}\\!=\\!\left(\\!\begin{array}[]{cc}-E_{k_{z}}^{+}&\boldsymbol{\Delta}k_{z}\\\
\boldsymbol{\Delta}k_{z}&E_{k_{z}}^{-}\end{array}\\!\right)\\!,\\!$ (84)

and $\hat{\mathscr{U}}_{3}=\text{diag}(\hat{\mathscr{V}},\hat{\mathscr{V}})$
where $\hat{\mathscr{V}}$ is the matrix of eigenvectors of
$\hat{\mathcal{D}}_{\mathbf{n}_{z}}(\pi,-1)$ given by

$\hat{\mathscr{V}}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cccc}0&1&0&-1\\\
0&1&0&1\\\ 1&0&-1&0\\\ 1&0&1&0\end{array}\right).$

We solve the eigenvalue problem for one of the decoupled subblocks in Eq.
(83). This problem is given by
$\hat{h}_{k_{z},\pm}^{-+}\hat{\Phi}(\xi,z)=\mathscr{E}_{\text{DP}}\hat{\Phi}(\xi,z)$
where $\hat{\Phi}(\xi,z)=(u,v)^{T}\text{exp}(\xi z)$ is the ansatz for the
decaying eigenspinor. $\xi$ denotes the localization factor, $|u|^{2}$
($|v|^{2}$) is the probability weight for electron (hole) states with
different magnetic quantum number $m_{j}=\pm 3/2$ ($m_{j}=\pm 1/2$). We
consider a semi-infinite system in $z\geq 0$ space. Therefore, $k_{z}$ is no
longer conserved and we use its real space representation
$k_{z}=k_{z}^{\dagger}=-i\partial_{z}$. In this case, the secular equation
$(\hat{h}_{-i\partial_{z},\pm}^{-+}-\mathscr{E}_{\text{DP}})\hat{\Phi}(\xi,z)=0$
can be evaluated by setting its determinant to zero,

$\Bigg{|}\begin{array}[]{cc}m^{\prime}\xi^{2}+\mu-\mathscr{E}_{\text{DP}}&-i\boldsymbol{\Delta}\xi\\\
-i\boldsymbol{\Delta}\xi&-(m\xi^{2}+\mu)-\mathscr{E}_{\text{DP}}\end{array}\Bigg{|}=0,$
(85)

where $m=\alpha+\beta/4$ and $m^{\prime}=\alpha+9\beta/4$.red The solution of
the secular equation yields

$\displaystyle\xi_{\pm}=$
$\displaystyle\sqrt{\frac{1}{2mm^{\prime}}\Big{(}\Lambda\pm\sqrt{\Lambda^{2}+4mm^{\prime}(\mathscr{E}_{\text{DP}}^{2}-\mu^{2})}\Big{)}},$
(86)

where
$\Lambda=\boldsymbol{\Delta}^{2}-\mathscr{E}_{\text{DP}}\vartheta^{-}-\mu\vartheta^{+}$
and $\vartheta^{\pm}=(m^{\prime}\pm m)/2$. In the absence of pairing, i.e.,
for $\boldsymbol{\Delta}=0$, the localization length becomes purely imaginary
leading to extended states. However, $\xi_{\pm}$ obtains a real component in
the topological phase induced by FE pairing. This leads to proper surface
state solutions.

To specify the components of $\hat{\Phi}(\xi,z)$, we use the secular equation,
and obtain

$\displaystyle u\equiv u_{\iota}=-i\boldsymbol{\Delta}\xi_{\iota}\rightarrow
v\equiv v_{\iota}$
$\displaystyle=m^{\prime}\xi_{\iota}^{2}+\mu-\mathscr{E}_{\text{DP}},$ (87)
$\displaystyle v\equiv V_{\iota}=-i\boldsymbol{\Delta}\xi_{\iota}\rightarrow
u\equiv U_{\iota}$
$\displaystyle=m\xi_{\iota}^{2}+\mu+\mathscr{E}_{\text{DP}}.$ (88)

Consequently, we have a pair of eigenspinors given by

$\displaystyle\\!\hat{\Phi}(\xi_{\iota},z)$
$\displaystyle=\Bigg{(}\begin{array}[]{c}u_{\iota}\\\
v_{\iota}\end{array}\Bigg{)}e^{\xi_{\iota}z},\ \
\hat{\Phi}^{\prime}(\xi_{\iota},z)=\Bigg{(}\begin{array}[]{c}U_{\iota}\\\
V_{\iota}\end{array}\Bigg{)}e^{\xi_{\iota}z}\\!.$ (93)

We can construct two sets of wave functions by the superposition of
eigenspinors $\hat{\Phi}(\xi_{\iota},z)$ $(\hat{\Phi}^{\prime}(\xi,z))$,

$\hat{\mathsf{\Psi}}_{1}(z)\\!\\!=\\!\\!\sum_{\iota=\pm}C_{\iota}\hat{\Phi}(\xi_{\iota},z),\
\ \
\hat{\mathsf{\Psi}}_{2}(z)\\!\\!=\\!\\!\sum_{\iota=\pm}Q_{\iota}\hat{\Phi}^{\prime}(\xi_{\iota},z),$
(94)

where the summations run over the decay factors and the coefficients of the
expansion are denoted by $C_{\iota}$ and $Q_{\iota}$. To have surface state
solutions, the wave functions and their first derivatives must vanish at the
interface of the system and far away from the interface, i.e.,

$\hat{\mathsf{\Psi}}(\infty)=\hat{\mathsf{\Psi}}^{\prime}(\infty)=0,\ \
\hat{\mathsf{\Psi}}(0)=\hat{\mathsf{\Psi}}^{\prime}(0)=0.$ (95)

Note that we assume just one interface in our analytical calculations. The
boundary conditions at $z=0$ give rise to two pairs of equations

$\displaystyle C_{-}$ $\displaystyle\Bigg{(}\begin{array}[]{c}u_{-}\\\
v_{-}\end{array}\Bigg{)}+C_{+}\Bigg{(}\begin{array}[]{c}u_{+}\\\
v_{+}\end{array}\Bigg{)}=0,$ (100) $\displaystyle Q_{-}$
$\displaystyle\Bigg{(}\begin{array}[]{c}U_{-}\\\
V_{-}\end{array}\Bigg{)}+Q_{+}\Bigg{(}\begin{array}[]{c}U_{+}\\\
V_{+}\end{array}\Bigg{)}=0.$ (105)

Re-arranging Eqs. (100) and (105), we arrive at

$\displaystyle\frac{C_{-}}{C_{+}}$
$\displaystyle=-\frac{\xi_{+}}{\xi_{-}}=-\frac{m^{\prime}\xi_{+}^{2}\\!+\\!\mu\\!-\mathscr{E}_{\text{DP}}}{m^{\prime}\xi_{-}^{2}\\!+\\!\mu\\!-\mathscr{E}_{\text{DP}}},$
(106) $\displaystyle\frac{Q_{-}}{Q_{+}}$
$\displaystyle=-\frac{\xi_{+}}{\xi_{-}}=-\frac{m\xi_{+}^{2}\\!+\\!\mu\\!+\mathscr{E}_{\text{DP}}}{m\xi_{-}^{2}\\!+\\!\mu\\!+\mathscr{E}_{\text{DP}}}.$
(107)

Combining Eqs. (106) and (107) results in the explicit formula for the energy
of the helical Dirac surface points $\mathscr{E}_{\text{DP}}$ given in Eq. (7)
of the main text.

Choosing $C_{-}=Q_{-}=\xi_{+}$ and $C_{+}=Q_{+}=-\xi_{-}$, this allows us to
derive the general eigenfunction corresponding to $\mathscr{E}_{\text{DP}}$,

$\hat{\mathsf{\Psi}}(z)=\mathcal{C}\xi_{+}\Bigg{(}\begin{array}[]{c}i\boldsymbol{\Delta}\xi_{-}\\\
m^{\prime}\xi_{-}^{2}+\mu-\mathscr{E}_{\text{DP}}\end{array}\Bigg{)}(e^{-\xi_{-}z}-e^{-\xi_{+}z}),$
(108)

where $\mathcal{C}$ is the normalization factor

$\mathcal{C}=\frac{1}{\sqrt{\left(|\kappa_{1}|^{2}+|\kappa_{2}|^{2}\right)}}\frac{1}{\sqrt{\int_{0}^{\infty}dz\
|f(z)|^{2}}}.$ (109)

##  Appendix D Effective 2D helical surface Hamiltonian

In this section, we derive the effective Hamiltonian for the 2D helical
surface states given in Eq. (8) of the main text. To do so, we need to project
the bare BdG Hamiltonian onto the helical Dirac surface states basis. Note
that $\hat{\mathsf{\Psi}}(z)$ in Eq. (108) is the eigenfunction corresponding
to the subblock matrix $\hat{h}_{{\bf k},+}^{-+}$, and we have defined
$f(z)\equiv(e^{-\xi_{-}z}-e^{-\xi_{+}z})$, $\kappa_{1}\equiv
i\boldsymbol{\Delta}\xi_{+}\xi_{-}$, and
$\kappa_{2}\equiv\xi_{+}(m^{\prime}\xi_{-}^{2}+\mu-\mathscr{E}_{\text{DP}})$.
To have a proper projection basis, we also need the eigenfunction for the
subblock $\hat{h}_{{\bf k},+}^{+-}$. It is given by
$\hat{\varphi}(z)=\mathcal{C}(\kappa_{3},\kappa_{4})^{T}\gamma(z)$ where
$\kappa_{3}\equiv\xi_{+}(m\xi_{-}^{2}+\mu+\mathscr{E}_{\text{DP}})$,
$\kappa_{4}\equiv-i\varDelta\xi_{+}\xi_{-}$, and $\gamma(z)=-f(z)$.

$\hat{\mathsf{\Psi}}(z)$ and $\hat{\varphi}(z)$ are $2\times 1$ column
vectors. In order to use them for the projection method, we convert them to
$8\times 1$ representation since the BdG Hamiltonian is a $8\times 8$ matrix.
This can be done through the transformation made by the $8\times 2$ columns of
the matrix
$\hat{\mathscr{U}}^{-1}=\\{\hat{\gamma}_{1},\hat{\gamma}_{2},\hat{\gamma}_{3},\hat{\gamma}_{4}\\}$.
The first and forth subblock matrices in Eq. (83) are identical corresponding
to $\hat{h}_{{\bf k},+}^{-+}$ and $\hat{h}_{{\bf k},-}^{-+}$, respectively.
Thus, they correspond to doubly degenerate helical surface states with
eigenvalue $\mathscr{E}_{\text{DP}}$. Therefore, their $8\times 1$
representations take the form

$\displaystyle\hat{\Gamma}_{1}^{-}$
$\displaystyle\equiv\hat{\gamma}_{1}\hat{\mathsf{\Psi}}(z)=\mathcal{C}f(z)(0,\kappa_{2},0,0,\kappa_{1},0,0,0)^{T},$
(110) $\displaystyle\hat{\Gamma}_{2}^{-}$
$\displaystyle\equiv\hat{\gamma}_{4}\hat{\mathsf{\Psi}}(z)=\mathcal{C}f(z)(0,0,\kappa_{2},0,0,0,0,\kappa_{1})^{T}.$
(111)

We repeat the above steps for $\hat{\varphi}(z)$ to obtain the proper basis
for the sectors $\hat{h}_{{\bf k},+}^{+-}$ and $\hat{h}_{{\bf k},-}^{+-}$,

$\displaystyle\hat{\Gamma}_{1}^{+}$
$\displaystyle\equiv\hat{\gamma}_{2}\hat{\varphi}(z)=\mathcal{C}\gamma(z)(\kappa_{4},0,0,0,0,\kappa_{3},0,0)^{T},$
(112) $\displaystyle\hat{\Gamma}_{2}^{+}$
$\displaystyle\equiv\hat{\gamma}_{3}\hat{\varphi}(z)=\mathcal{C}\gamma(z)(0,0,0,\kappa_{4},0,0,\kappa_{3},0)^{T}.$
(113)

We use Eqs. (110-113) as the proper orthonormal set of eigenfunctions to
project the bulk superconducting Hamiltonian to the surface. Note that the
orthonormality condition reads
$\int_{0}^{\infty}dz[\hat{\mathscr{Y}}^{\nu}]^{\dagger}\hat{\mathscr{Y}}^{\nu}=\hat{\sigma}_{0}$
where
$\hat{\mathscr{Y}}^{\pm}=\\{\hat{\Gamma}_{1}^{\pm},\hat{\Gamma}_{2}^{\pm}\\}$.

To derive the effective Hamiltonian for the 2D helical surface states, we
consider the conserved wave vectors $k_{x}$ and $k_{y}$ in the BdG Hamiltonian
to be small close to the $\Gamma$ point. Then, we project
$\hat{H}(k_{x},k_{y},-i\partial_{z})$ onto the basis of the helical Dirac
surface states. Thus, the effective Hamiltonian at finite excitation energies
becomes

$J$ | $O_{h}$ | $\eta$ | $\hat{\mathcal{O}}_{\eta}(\hat{J})$ | $\hat{\eta}$ | $(L,S)$ | $(L,S)$
---|---|---|---|---|---|---
$0$ | $A_{1g,u}$ | $A_{1g,u}$ | $\frac{1}{2}\hat{I}_{4}$ | $\hat{A}_{1g,u}\\!=\\!\hat{\mathcal{N}}_{0,0}$ | $(0,0)$ | $(1,1)$
$1$ | $T_{1u}$ | $T_{1u}^{(1)}$ | $\frac{1}{\sqrt{5}}\hat{J}_{z}$ | $\hat{T}_{1u}^{(1)}\\!=\\!\hat{\mathcal{N}}_{1,0}$ | $\times$ | $(1,1)$
|  | $T_{1u}^{(2)}$ | $\frac{1}{\sqrt{5}}\hat{J}_{y}$ | $\hat{T}_{1u}^{(2)}\\!=\\!\frac{i}{\sqrt{2}}(\hat{\mathcal{N}}_{1,-1}+\hat{\mathcal{N}}_{1,1})$ | $\times$ | $(1,1)$
|  | $T_{1u}^{(3)}$ | $\frac{1}{\sqrt{5}}\hat{J}_{x}$ | $\hat{T}_{1u}^{(3)}\\!=\\!\frac{1}{\sqrt{2}}(\hat{\mathcal{N}}_{1,-1}-\hat{\mathcal{N}}_{1,1})$ | $\times$ | $(1,1)$
$2$ | $E_{g,u}$ | $E_{g,u}^{(1)}$ | $\frac{1}{6}(3\hat{J}_{z}^{2}-\hat{\mathbf{J}})$ | $\hat{E}_{g,u}^{(1)}\\!=\\!\hat{\mathcal{N}}_{2,0}$ | $(0,2)$ | $(1,1),(1,3)$
|  | $E_{g,u}^{(2)}$ | $\frac{1}{2\sqrt{3}}(\hat{J}_{x}^{2}-\hat{J}_{y}^{2})$ | $\hat{E}_{g,u}^{(2)}\\!=\\!\frac{1}{\sqrt{2}}(\hat{\mathcal{N}}_{2,-2}+\hat{\mathcal{N}}_{2,2})$ | $(0,2)$ | $(1,1),(1,3)^{*}$
| $T_{2g,u}$ | $T_{2g,u}^{(1)}$ | $\frac{1}{\sqrt{3}}(\hat{J}_{x}\hat{J}_{y}+\hat{J}_{y}\hat{J}_{x})$ | $\hat{T}_{2g,u}^{(1)}\\!=\\!\frac{i}{\sqrt{2}}(\hat{\mathcal{N}}_{2,-2}-\hat{\mathcal{N}}_{2,2})$ | $(0,2)$ | $(1,1),(1,3)^{*}$
|  | $T_{2g,u}^{(2)}$ | $\frac{1}{\sqrt{3}}(\hat{J}_{z}\hat{J}_{x}+\hat{J}_{x}\hat{J}_{z})$ | $\hat{T}_{2g,u}^{(2)}\\!=\\!\frac{i}{\sqrt{2}}(\hat{\mathcal{N}}_{2,-1}-\hat{\mathcal{N}}_{2,1})$ | $(0,2)$ | $(1,1),(1,3)^{*}$
|  | $T_{2g,u}^{(3)}$ | $\frac{1}{\sqrt{3}}(\hat{J}_{y}\hat{J}_{z}+\hat{J}_{z}\hat{J}_{y})$ | $\hat{T}_{2g,u}^{(3)}\\!=\\!\frac{i}{\sqrt{2}}(\hat{\mathcal{N}}_{2,-1}+\hat{\mathcal{N}}_{2,1})$ | $(0,2)$ | $(1,1),(1,3)^{*}$
$3$ | $A_{2u}$ | $A_{2u}$ | $\frac{1}{\sqrt{3}}(\hat{J}_{x}\hat{J}_{y}\hat{J}_{z}+\hat{J}_{z}\hat{J}_{y}\hat{J}_{x})$ | $\hat{A}_{2u}\\!=\\!\frac{i}{\sqrt{2}}(\hat{\mathcal{N}}_{2,-2}-\hat{\mathcal{N}}_{2,2})$ | $\times$ | $(1,3)^{*}$
| $T_{1u}$ | $T_{1u}^{(1)}$ | $\frac{4}{\sqrt{365}}\hat{J}_{z}^{3}$ | $\hat{T}_{1u}^{(1)}\\!=\\!\frac{12}{5\sqrt{73}}\hat{\mathcal{N}}_{3,0}$ | $\times$ | $(1,3)$
|  | $T_{1u}^{(2)}$ | $\frac{4}{\sqrt{365}}\hat{J}_{y}^{3}$ | $\hat{T}_{1u}^{(2)}\\!=\\!\frac{3}{5\sqrt{73}i}[\sqrt{5}(\hat{\mathcal{N}}_{3,-3}+\hat{\mathcal{N}}_{3,3})+\sqrt{3}(\hat{\mathcal{N}}_{3,1}+\hat{\mathcal{N}}_{3,-1})]$ | $\times$ | $(1,3)$
|  | $T_{1u}^{(3)}$ | $\frac{4}{\sqrt{365}}\hat{J}_{x}^{3}$ | $\hat{T}_{1u}^{(3)}\\!=\\!\frac{3}{5\sqrt{73}}[\sqrt{5}(\hat{\mathcal{N}}_{3,-3}-\hat{\mathcal{N}}_{3,3})+\sqrt{3}(\hat{\mathcal{N}}_{3,1}-\hat{\mathcal{N}}_{3,-1})]$ | $\times$ | $(1,3)$
| $T_{2u}$ | $T_{2u}^{(1)}$ | $\frac{1}{\sqrt{3}}\lceil\hat{J}_{z}(\hat{J}_{x}^{2}-\hat{J}_{y}^{2})\rfloor$ | $\hat{T}_{2u}^{(1)}\\!=\\!\frac{1}{\sqrt{2}}(\hat{\mathcal{N}}_{3,-2}+\hat{\mathcal{N}}_{3,2})$ | $\times$ | $(1,3)^{*}$
|  | $T_{2u}^{(2)}$ | $\frac{1}{\sqrt{3}}\lceil\hat{J}_{x}(\hat{J}_{y}^{2}-\hat{J}_{z}^{2})\rfloor$ | $\hat{T}_{2u}^{(2)}\\!=\\!\frac{1}{4}[\sqrt{3}(\hat{\mathcal{N}}_{3,3}-\hat{\mathcal{N}}_{3,-3})+\sqrt{5}(\hat{\mathcal{N}}_{3,1}-\hat{\mathcal{N}}_{3,-1})]$ | $\times$ | $(1,3)^{*}$
|  | $T_{2u}^{(3)}$ | $\frac{1}{\sqrt{3}}\lceil\hat{J}_{y}(\hat{J}_{z}^{2}-\hat{J}_{x}^{2})\rfloor$ | $\hat{T}_{2u}^{(3)}\\!=\\!\frac{1}{4i}[\sqrt{3}(\hat{\mathcal{N}}_{3,3}+\hat{\mathcal{N}}_{3,-3})-\sqrt{5}(\hat{\mathcal{N}}_{3,1}+\hat{\mathcal{N}}_{3,-1})]$ | $\times$ | $(1,3)^{*}$
Table 3: Decomposition of total angular momentum $J$ (first column) in the
irreducible representation (irrep) of $O_{h}$ symmetry (second column). The
dimension of an irrep is distinguished by the number of components $\eta$
given in the third column. The fourth column denotes the normalized
irreducible basis matrices for a component of $O_{h}$ symmetry in $j=3/2$
representation. The fifth column demonstrates the correspondence between the
matrix form for components of a cubic irrep $\hat{\eta}$ and the components of
$SO(3)$ symmetry, namely $\hat{\mathcal{N}}_{J,m_{j}}$ being the total angular
momentum tensor matrices. The last two columns indicate the spin $S$ and
orbital $L$ angular momenta of Cooper pairs associated to $J$ and components
of a given irrep. The Fermi statistics forces both $L$ and $S$ to be even
($g$) or odd ($u$). The symmetrization of the basis matrices is denoted by
$\lceil\hat{A}\hat{B}\hat{C}\rfloor=(\hat{A}\hat{B}\hat{C}+\hat{A}\hat{C}\hat{B}+\hat{B}\hat{C}\hat{A}+\hat{B}\hat{A}\hat{C}+\hat{C}\hat{A}\hat{B}+\hat{C}\hat{B}\hat{A})/3!$.
The pairing channels satisfying TPT at FEs are marked by $(*)$.

$\displaystyle\hat{\mathtt{H}}(k_{x},k_{y})$
$\displaystyle=\int_{0}^{\infty}dz\hat{\mathscr{Y}}^{\nu\dagger}(z)\hat{H}(k_{x},k_{y},-i\partial_{z})\hat{\mathscr{Y}}^{\nu}(z)$
$\displaystyle=\left(\begin{array}[]{cc}A_{1,1}&A_{1,2}\\\
A_{2,1}&A_{2,2}\end{array}\right),$ (116)

where the matrix elements of $\hat{\mathtt{H}}(k_{x},k_{y})$ are given by

$\displaystyle A_{i,j}$
$\displaystyle=\int_{0}^{\infty}dz\hat{\Gamma}_{i}^{\dagger}(z)\hat{H}(k_{x},k_{y},-i\partial_{z})\hat{\Gamma}_{j}(z),$
(117)

with $i,j\in\\{1,2\\}$. $\hat{\mathtt{H}}(k_{x},k_{y})$ should be Hermitian,
thus, the matrix elements in Eq. (116) must fulfill the relations

$\displaystyle A_{1,1}$ $\displaystyle=A_{1,1}^{*}=A_{2,2}=A_{2,2}^{*},$ (118)
$\displaystyle A_{1,2}$ $\displaystyle=A_{2,1}^{*},\ \ A_{2,1}=A_{1,2}^{*}.$
(119)

After straightforward algebra, we arrive at the 2D effective Hamiltonian for
the helical surface states given in Eq. (8) of the main text. Note that the
group velocity of the helical topological Dirac surface states at finite
excitation energies takes the form

$\varsigma_{2}=\frac{\sqrt{3}}{4}\frac{\Delta}{(|\kappa_{1}|^{2}+|\kappa_{2}|^{2})}\text{Im}\left[\left(\kappa_{2}\kappa_{1}^{*}-\kappa_{1}\kappa_{2}^{*}\right)\right].$
(120)

It is clear that $\varsigma_{2}$ depends on the components of the wave
function with direct proportionality to the pairing strength $\Delta$.

##  Appendix E Pairing matrices

In this section, we demonstrate the method for obtaining the pairing matrices
in cubic symmetry, given in Table I of the main paper. The cubic pairing
matrices with dependency of momenta up to the linear order can be derived
through the relation between SO(3) (Yang  _et al._ , 2016; Fang  _et al._ ,
2015; Savary  _et al._ , 2017) and cubic (PairingMatrices-2011; Savary  _et
al._ , 2017; Venderbos  _et al._ , 2018; Yu and Liu, 2018; Kawakami  _et al._
, 2018) symmetries as we elaborate in the following. In addition, in Ref.
(Timm and Bhattacharya, 2021), a general group-theoretical method to derive
the symmetry-allowed pairing matrices for any local degrees of freedom (spin,
orbital, and basis site) is introduced. The method relies on the reduction of
product representation of the crystallographic point group.

The irreducible representations of SO(3) symmetry are labeled by the total
angular momentum $J$, which combines orbital $L$ and spin $S$ angular momenta.
The Cooper pairs formed by two electrons with $j=3/2$ total angular momentum
can have singlet ($J=0$), triplet ($J=1$), quintet ($J=2$), or septet $(J=3)$
angular momenta. Considering orbital angular momentum up to the p-wave
channel, i.e., $L\in\\{0,1\\}$, the intrinsic spin angular momentum can take
values $S\in\\{0,1,2,3\\}$ according to relation $|L-S|\leq J\leq|L+S|$. The
correspondence between the components of a cubic irrep and the SO(3) symmetry
is given by the relation (Tamura  _et al._ , 2020),

$\\!\\!\hat{\eta}=\\!\\!\\!\\!\sum_{m_{j_{1}},m_{j_{2}}}\\!\sum_{m_{J}}\\!\big{[}\big{(}\hat{\mathcal{O}}_{\eta}(\hat{J})\hat{\mathcal{R}}\big{)}_{m_{j_{1}},m_{j_{2}}}\\!\langle
J,m_{J}|m_{j_{1}},m_{j_{2}}\rangle\big{]}\hat{\mathcal{N}}_{J,m_{J}},\\!\\!$
(121)

where $\langle J,m_{J}|m_{j_{1}},m_{j_{2}}\rangle$ are the Clebsch-Gordan
coefficients, $m_{j_{1}},m_{j_{2}}\in\\{\pm 3/2,\pm 1/2\\}$,
$\hat{\mathcal{R}}=e^{i\pi\hat{J}_{y}}$, $\hat{\mathcal{N}}_{J,m_{J}}$ denotes
the multipole matrix for the total angular momentum labeled with magnetic
quantum number $m_{J}\in\\{-J,\ldots,J\\}$, and $\hat{\eta}$ is the matrix
representation for the component of an irrep with the normalized basis matrix
$\hat{\mathcal{O}}_{\eta}(\hat{J})$ given in the third and fourth columns of
Table 3, respectively. Note that the fifth column of Table 3 is obtained
according to Eq. (121). The explicit matrix form of $\hat{\mathcal{N}}_{J,m}$
can be derived by the expansion

$\hat{\mathcal{N}}_{J,m}=\sum_{m_{L},m_{S}}\langle m_{L},m_{S}|J,m_{J}\rangle
Y_{m_{L}}^{L}(\mathbf{k})\hat{\mathcal{S}}_{m_{S}}^{S},$ (122)

where $\langle m_{L},m_{S}|J,m_{J}\rangle$ denotes the Clebsch-Gordan
coefficient, $Y_{m_{L}}^{L}(\mathbf{k})$ and $\hat{\mathcal{S}}_{m_{S}}^{S}$
are spherical harmonics and irreducible spin tensor matrices labeled with
$m_{L}\in[-L,\ldots,L]$ and $m_{S}\in[-S,\ldots,S]$ as the orbital axial
angular momentum and spin magnetic quantum number, respectively. Since we are
interested in the odd-parity pairings up to linear momenta, the orbital
angular momentum is $L=1$ ($p$-wave) with the relative spherical harmonics
$Y_{1,\pm 1}(\mathbf{k})=\mp\sqrt{3/8\pi}k_{\pm}/|{\bf k}|$ and
$Y_{1,0}(\mathbf{k})=\sqrt{3/4\pi}k_{z}/|{\bf k}|$. Moreover, the derivation
method for obtaining $\hat{\mathcal{S}}_{m_{S}}^{S}$ follows Refs. (Venderbos
_et al._ , 2018; Yu and Liu, 2018; Bahari  _et al._ , 2022). Using the fifth
columns of Table 3 and Eq. (122), the cubic pairing matrices can be derived
straightforwardly, giving

$\hat{\Delta}_{\mathbf{k}}=|\mathbf{k}|\,\hat{\eta}\,\hat{\mathcal{R}}.$ (123)

We start with the $A_{2u}$ irrep to find its pairing matrix. $A_{2u}$ is a
one-dimensional irrep appearing only for the $J=3$ irrep corresponding to
$L=1$ and $S=3$. The matrix representation for $A_{2u}$ is given by

$\displaystyle\hat{A}_{2u}\\!$
$\displaystyle=\\!\frac{i}{\sqrt{2}}(\hat{\mathcal{N}}_{2,-2}-\hat{\mathcal{N}}_{2,2})\\!=\\!\frac{\sqrt{2}}{12i}[Y_{1}^{1}(3\hat{\mathcal{S}}_{-3}^{3}-\sqrt{15}\hat{\mathcal{S}}_{1}^{3})$
$\displaystyle+$
$\displaystyle\sqrt{12}Y_{0}^{1}(\hat{\mathcal{S}}_{-2}^{3}+\hat{\mathcal{S}}_{2}^{3})\\!+\\!Y_{-1}^{1}(3\hat{\mathcal{S}}_{3}^{3}-\sqrt{15}\hat{\mathcal{S}}_{-1}^{3})],\\!\\!\\!$
(124)

where we dropped the momentum dependency of spherical harmonics. Spin
multipole matrices, which fulfill the relation
$\hat{\mathcal{S}}_{-m_{S}}^{S}=(-1)^{m_{S}}(\hat{\mathcal{S}}_{m_{S}}^{S})^{T}$,
are given by

$\\!\\!\hat{\mathcal{S}}_{0}^{3}\\!=\frac{\sqrt{5}}{10}\left(\\!\\!\begin{array}[]{cccc}1&\\!\\!0&0&0\\\
0&\\!\\!-3&0&0\\\ 0&\\!\\!0&3&0\\\
0&\\!\\!0&0&-1\end{array}\\!\\!\right)\\!,\hat{\mathcal{S}}_{1}^{3}\\!=\\!\frac{\sqrt{5}}{5}\\!\left(\\!\\!\begin{array}[]{cccc}0&-1&0&\\!\\!0\\\
0&0&\sqrt{3}&\\!\\!0\\\ 0&0&0&\\!\\!-1\\\
0&0&0&\\!\\!0\end{array}\\!\\!\\!\right)\\!\\!,\\!\\!\\!$ (125)

and

$\hat{\mathcal{S}}_{2}^{3}=\frac{\sqrt{2}}{2}\left(\begin{array}[]{cccc}0&0&1&0\\\
0&0&0&-1\\\ 0&0&0&0\\\ 0&0&0&0\end{array}\right)\\!,\
\hat{\mathcal{S}}_{3}^{3}=\left(\begin{array}[]{cccc}0&0&0&-1\\\ 0&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\end{array}\right).$ (126)

Eventually, the pairing matrix for the $A_{2u}$ irrep takes the form (Brydon
_et al._ , 2016)

$\displaystyle\hat{\Delta}_{{\bf k}}$
$\displaystyle=|\mathbf{k}|\,\hat{A}_{2u}\,\hat{\mathcal{R}}$
$\displaystyle=\Delta\left(\begin{array}[]{cccc}\frac{3}{4}k_{-}&\frac{\sqrt{3}}{2}k_{z}&\frac{\sqrt{3}}{4}k_{+}&0\\\
\frac{\sqrt{3}}{2}k_{z}&\frac{3}{4}k_{+}&0&-\frac{\sqrt{3}}{4}k_{-}\\\
\frac{\sqrt{3}}{4}k_{+}&0&-\frac{3}{4}k_{-}&\frac{\sqrt{3}}{2}k_{z}\\\
0&-\frac{\sqrt{3}}{4}k_{-}&\frac{\sqrt{3}}{2}k_{z}&-\frac{3}{4}k_{+}\end{array}\right).$
(131)

In the next step, we consider the three-dimensional $T_{2u}$ irrep. It appears
for the $J=3$ and $J=5$ irreps. Due the the nature of $j=3/2$ electron, only
the septet component is allowed. In this case, the first component of the
$T_{2u}$ irrep up to linear order of momenta takes the form

$\displaystyle\hat{T}_{2u}^{(1)}\\!$
$\displaystyle=\\!\frac{1}{\sqrt{2}}(\hat{\mathcal{N}}_{3,-2}+\hat{\mathcal{N}}_{3,2})\\!=\\!\frac{1}{6\sqrt{2}}[Y_{-1}^{1}(\sqrt{15}\hat{\mathcal{S}}_{-1}^{3}+3\hat{\mathcal{S}}_{3}^{3})$
$\displaystyle+$
$\displaystyle\sqrt{12}Y_{0}^{1}(\hat{\mathcal{S}}_{2}^{3}-\hat{\mathcal{S}}_{-2}^{3})-Y_{1}^{1}(3\hat{\mathcal{S}}_{-3}^{3}+\sqrt{15}\hat{\mathcal{S}}_{1}^{3})].$
(132)

The second component is

$\displaystyle\hat{T}_{2u}^{(2)}=$
$\displaystyle\frac{1}{4}[\sqrt{3}(\hat{\mathcal{N}}_{3,3}-\hat{\mathcal{N}}_{3,-3})+\sqrt{5}(\hat{\mathcal{N}}_{3,1}-\hat{\mathcal{N}}_{3,-1})]$
$\displaystyle=$
$\displaystyle\frac{1}{24}[\sqrt{3}Y_{-1}^{1}((5\hat{\mathcal{S}}_{2}^{3}-3\hat{\mathcal{S}}_{-2}^{3})-\sqrt{30}\hat{\mathcal{S}}_{0}^{3})$
$\displaystyle+\sqrt{3}Y_{1}^{1}((5\hat{\mathcal{S}}_{-2}^{3}-3\hat{\mathcal{S}}_{2}^{3})-\sqrt{30}\hat{\mathcal{S}}_{0}^{3})$
$\displaystyle+Y_{0}^{1}(9(\hat{\mathcal{S}}_{-3}^{3}+\hat{\mathcal{S}}_{3}^{3})+\sqrt{15}(\hat{\mathcal{S}}_{-1}^{3}+\hat{\mathcal{S}}_{1}^{3}))],$
(133)

and the third component

$\displaystyle\hat{T}_{2u}^{(3)}=$
$\displaystyle\frac{1}{4i}[\sqrt{3}(\hat{\mathcal{N}}_{3,3}+\hat{\mathcal{N}}_{3,-3})-\sqrt{5}(\hat{\mathcal{N}}_{3,1}+\hat{\mathcal{N}}_{3,-1})]$
$\displaystyle=$
$\displaystyle\frac{i}{24}(\sqrt{3}Y_{1}^{1}(3\hat{\mathcal{S}}_{2}^{3}-5\hat{\mathcal{S}}_{-2}^{3}-\sqrt{30}\hat{\mathcal{S}}_{0}^{3})$
$\displaystyle+\sqrt{3}Y_{-1}^{1}(5\hat{\mathcal{S}}_{2}^{3}-3\hat{\mathcal{S}}_{-2}^{3}+\sqrt{30}\hat{\mathcal{S}}_{0}^{3})$
$\displaystyle+Y_{0}^{1}(9(\hat{\mathcal{S}}_{-3}^{3}-\hat{\mathcal{S}}_{3}^{3})+\sqrt{15}(\hat{\mathcal{S}}_{1}^{3}-\hat{\mathcal{S}}_{-1}^{3}))).$
(134)

Finally, the explicit pairing matrices for the $T_{2u}$ irrep become

$\displaystyle\hat{\Delta}_{{\bf k}}$
$\displaystyle=|\mathbf{k}|\,\hat{T}_{2u}^{(1)}\,\hat{\mathcal{R}}$
$\displaystyle=\\!\Delta\\!\left(\\!\\!\begin{array}[]{cccc}\frac{3}{4}k_{-}&\frac{\sqrt{3}}{2}k_{z}&\frac{\sqrt{3}}{4}k_{+}&0\\\
\frac{\sqrt{3}}{2}k_{z}&\frac{3}{4}k_{+}&0&\frac{1}{4}\sqrt{3}k_{-}\\\
\frac{\sqrt{3}}{4}k_{+}&0&\frac{3}{4}k_{-}&-\frac{\sqrt{3}}{2}k_{z}\\\
0&\frac{\sqrt{3}}{4}k_{-}&-\frac{\sqrt{3}}{2}k_{z}&\frac{3}{4}k_{+}\end{array}\\!\\!\right),\\!\\!$
(139) $\displaystyle\hat{\Delta}_{{\bf k}}$
$\displaystyle=|\mathbf{k}|\,\hat{T}_{2u}^{(2)}\,\hat{\mathcal{R}}$
$\displaystyle=\Delta\left(\begin{array}[]{cccc}3k_{z}&\frac{1}{\sqrt{3}}k_{-}^{\prime}&\frac{1}{\sqrt{3}}k_{z}&ik_{y}\\\
\frac{1}{\sqrt{3}}k_{+}^{\prime}&k_{z}&3ik_{y}&\frac{1}{\sqrt{3}}k_{z}\\\
\frac{1}{\sqrt{3}}k_{z}&3ik_{y}&k_{z}&\frac{-1}{\sqrt{3}}k_{+}^{\prime}\\\
ik_{y}&\frac{1}{\sqrt{3}}k_{z}&\frac{-1}{\sqrt{3}}k_{+}^{\prime}&3k_{z}\end{array}\right),$
(144) $\displaystyle\hat{\Delta}_{{\bf k}}$
$\displaystyle=|\mathbf{k}|\,\hat{T}_{2u}^{(3)}\,\hat{\mathcal{R}}$
$\displaystyle=\Delta\left(\begin{array}[]{cccc}3k_{z}&\frac{-1}{\sqrt{3}}k_{-}^{\prime\prime}&\frac{-1}{\sqrt{3}}k_{z}&-k_{x}\\\
\frac{-1}{\sqrt{3}}k_{-}^{\prime\prime}&-k_{z}&-3k_{x}&\frac{1}{\sqrt{3}}k_{z}\\\
\frac{-1}{\sqrt{3}}k_{z}&-3k_{x}&k_{z}&\frac{-1}{\sqrt{3}}k_{+}^{\prime\prime}\\\
-k_{x}&\frac{1}{\sqrt{3}}k_{z}&\frac{-1}{\sqrt{3}}k_{+}^{\prime\prime}&-3k_{z}\end{array}\right),$
(149)

where $k_{\pm}^{\prime}=4k_{x}\pm ik_{y}$ and $k_{\pm}^{\prime\prime}=k_{x}\pm
4ik_{y}$.

$E_{u}$ $(T_{2u})$ is a two- (three-) dimensional irrep. It appears for the
$J=2$ and $J=4$ irreps. In the $j=3/2$ representation, only the $J=2$ channel
is allowed. Note that the second component of $E_{u}$ and all the components
of $T_{2u}$ can undergo TPTs at FEs. Thus, the expansion for
$\hat{E}_{u}^{(2)}$ is

$\displaystyle\\!\\!\hat{E}_{u}^{(2)}$
$\displaystyle\\!=\frac{1}{\sqrt{2}}(\hat{\mathcal{N}}_{2,-2}+\hat{\mathcal{N}}_{2,2})\\!=\frac{\sqrt{14}}{42}[\sqrt{3}Y_{1}^{1}(\sqrt{15}\hat{\mathcal{S}}_{-3}^{3}+\hat{\mathcal{S}}_{1}^{3})$
$\displaystyle+\sqrt{3}$ $\displaystyle
Y_{-1}^{1}(\hat{\mathcal{S}}_{-1}^{3}+\sqrt{15}\hat{\mathcal{S}}_{3}^{3})-\sqrt{15}Y_{0}^{1}(\hat{\mathcal{S}}_{-2}^{3}+\hat{\mathcal{S}}_{2}^{3})].$
(150)

The pairing matrix for $E_{u}^{(2)}$ becomes

$\displaystyle\hat{\Delta}_{{\bf k}}$
$\displaystyle=|\mathbf{k}|\,\hat{E}_{u}^{(2)}\,\hat{\mathcal{R}}$
$\displaystyle=\Delta\left(\begin{array}[]{c<EMAIL_ADDRESS>-5k_{z}\hfil\hskip-1.0&-\sqrt{3}k_{+}\hfil\hskip-1.0&0&k_{-}\\\
-k_{+}\hfil\hskip-1.0&0\hfil\hskip-1.0&\sqrt{3}k_{-}&-5k_{z}\\\
0\hfil\hskip-1.0&k_{-}\hfil\hskip-1.0&-5k_{z}&-5\sqrt{3}k_{+}\end{array}\right).$
(155)

Likewise, the pairing matrix for the first component of the $T_{2u}$ irrep can
be straightforwardly derived

$\displaystyle\\!\\!\hat{\Delta}_{{\bf k}}$
$\displaystyle=|\mathbf{k}|\,\hat{T}_{2u}^{(1)}\,\hat{\mathcal{R}}$
$\displaystyle=\Delta\\!\left(\\!\\!\begin{array}[]{cccc}-5\sqrt{3}k_{-}&\\!\\!5k_{z}&\\!\\!k_{+}&\\!\\!0\\\
5k_{z}&\\!\\!\sqrt{3}k_{+}&\\!\\!0&\\!\\!k_{-}\\\
k_{+}&\\!\\!0&\\!\\!\sqrt{3}k_{-}&\\!\\!-5k_{z}\\\
0&\\!\\!k_{-}&\\!\\!-5k_{z}&\\!\\!-5\sqrt{3}k_{+}\end{array}\\!\\!\right)\\!,\\!\\!\\!\\!$
(160)

and the second component takes the form

$\displaystyle\\!\\!\hat{\Delta}_{{\bf k}}$
$\displaystyle=|\mathbf{k}|\,\hat{T}_{2u}^{(2)}\,\hat{\mathcal{R}}$
$\displaystyle=\Delta\left(\\!\\!\begin{array}[]{cccc}0&-5k_{-}&4k_{z}&\sqrt{3}k_{x}\\\
-5k_{-}&4\sqrt{3}k_{z}&3\sqrt{3}k_{x}&-4k_{z}\\\
4k_{z}&3\sqrt{3}k_{x}&-4\sqrt{3}k_{z}&-5k_{+}\\\
\sqrt{3}k_{x}&-4k_{z}&-5k_{+}&0\end{array}\\!\\!\right)\\!.\\!\\!\\!$ (165)

Finally, the pairing potential for the third component of $T_{2u}$ becomes

$\displaystyle\\!\\!\hat{\Delta}_{{\bf k}}$
$\displaystyle=|\mathbf{k}|\,\hat{T}_{2u}^{(3)}\,\hat{\mathcal{R}}$
$\displaystyle=\Delta\\!\left(\\!\\!\begin{array}[]{cccc}0&\\!\\!-5k_{-}&4k_{z}&\\!\\!i\sqrt{3}k_{y}\\\
-5k_{-}&\\!\\!4\sqrt{3}k_{z}&3i\sqrt{3}k_{y}&\\!\\!4k_{z}\\\
4k_{z}&\\!\\!3i\sqrt{3}k_{y}&4\sqrt{3}k_{z}&\\!\\!5k_{+}\\\
i\sqrt{3}k_{y}&\\!\\!4k_{z}&5k_{+}&\\!\\!0\end{array}\\!\\!\right).\\!\\!\\!\\!$
(170)

To sum up, the relation between cubic point group symmetry and full SO(3)
symmetry enabled us to explicitly obtain the cubic pairing matrices. This
method can also be applied to other point group symmetries.