Theory of surface spectroscopy for noncentrosymmetric superconductors

TL;DR

This paper develops a theoretical framework for understanding surface states and topological phases in noncentrosymmetric superconductors with various crystal symmetries, revealing complex surface phenomena and phase transitions.

Contribution

It provides a self-consistent calculation of the order parameter near surfaces and explores the topological phase diagram for different Fermi surface configurations.

Findings

Second phase transition below T_c for certain parameters

Intricate Andreev bound state structures at surfaces

Sharp topological transition between open and closed Fermi surfaces

Abstract

We study noncentrosymmetric superconductors with the tetrahedral $T_d$, tetragonal $C_{4v}$, and cubic point group $O$. The order parameter is computed self-consistently in the bulk and near a surface for several different singlet to triplet order parameter ratios. It is shown that a second phase transition below $T_c$ is possible for certain parameter values. In order to determine the surface orientation's effect on the order parameter suppression, the latter is calculated for a range of different surface orientations. For selected self-consistent order parameter profiles the surface density of states is calculated showing intricate structure of the Andreev bound states (ABS) as well as spin polarization. The topology's effect on the surface states and the tunnel conductance is thoroughly investigated, and a topological phase diagram is constructed for open and closed Fermi surfaces showing a sharp transition between the two for the cubic point group $O$.

Figures

Tables

Table 1. Table 1: The scaled bulk singlet to triplet ratios, r Δ bulk ≡ Δ s / ( Δ t ​ max ⁡ | 𝐥 ​ ( 𝐤 F ) | ) subscript superscript 𝑟 bulk Δ subscript Δ 𝑠 subscript Δ 𝑡 𝐥 subscript 𝐤 𝐹 r^{\text{bulk}}_{\Delta}\equiv\Delta_{s}/(\Delta_{t}\max|{\bf l}({\bf k}_{F})|) , chosen for the different point groups, O 𝑂 O , C 4 ​ v subscript 𝐶 4 𝑣 C_{4v} , and T d subscript 𝑇 𝑑 T_{d} , and their respective g 2 subscript 𝑔 2 g_{2} values used in this work.

	$g_2$					$r_{\Delta }^{\text{bulk}}$
$O$	$0.0$	0	0.26	0.38	0.50	0.62	0.74	0.86	0.98	1.1
$O$	$0.7$	0	0.20	0.33	0.46	0.59	0.72	0.84	0.97	1.1
$O$	$1.03$	0	0.09	0.23	0.38	0.52	0.67	0.81	0.96	1.1
$O$	$2.5$	0	0.17	0.30	0.44	0.57	0.70	0.83	0.97	1.1
$C_{4v}$	$0.0$	0	0.14	0.28	0.41	0.55	0.69	0.83	0.96	1.1
$C_{4v}$	$4.0$	0	0.14	0.28	0.41	0.55	0.69	0.83	0.96	1.1
$T_d$	N/A	0	0.14	0.28	0.41	0.55	0.69	0.83	0.96	1.1

Equations

$${\mathbf{l}}_{\mathbf{k}}\equiv \mathbf{l}(\mathbf{k})=g\mathbf{l}(g^{-1}\mathbf{k}),$$

$${\hat{\mathcal{H}}}_{\mathbf{k}}=\sum _{\mathbf{k}\alpha \beta }{c}_{\mathbf{k}\alpha }^{†}{\left({\xi }_{\mathbf{k}}{\sigma }_0+\alpha {\mathbf{l}}_{\mathbf{k}}\cdot \mathbf{\sigma }\right)}_{\alpha \beta }{c}_{\mathbf{k}\beta }$$

$${ϵ}_{\mathbf{k}}^{\text{CUB}}=t_1\left[\cos \left(k_x\right)+\cos \left(k_y\right)+\cos \left(k_z\right)\right]$$

$${ϵ}_{\mathbf{k}}^{\text{BCC}}=8t_1\cos \left(k_x/2\right)\cos \left(k_y/2\right)\cos \left(k_z/2\right),$$

$$U_{\mathbf{k}}=\left(\begin{array}{cc}\cos \left(\frac{{\theta }_l}{2}\right)& e^{-i{\varphi }_l}\sin \left(\frac{{\theta }_l}{2}\right)\\ -e^{i{\varphi }_l}\sin \left(\frac{{\theta }_l}{2}\right)& \cos \left(\frac{{\theta }_l}{2}\right)\end{array}\right),$$

$${\hat{\mathcal{H}}}_{\mathbf{k}}=\sum _{\mathbf{k}\lambda }{\xi }_{\mathbf{k}}^{\lambda }{b}_{\mathbf{k}\lambda }^{†}{b}_{\mathbf{k}\lambda }$$

$$i{\mathbf{v}}_F\cdot {\nabla }_{\mathbf{R}}{\hat{g}}^{\text{R,A,M}}+[z{\hat{\tau }}_3-\hat{\Delta }-{\hat{v}}_{\text{SO}},\hat{g}{]}^{\text{R,A,M}}=\hat{0}$$

$$\hat{\Delta }=\left(\begin{array}{cc}0& \Delta \\ \tilde{\Delta }& 0\end{array}\right)$$

$$\Delta (\mathbf{k})={\mathcal{Y}}_{\mathbf{k}}({\Delta }_s+{\Delta }_t{\mathbf{l}}_{\mathbf{k}}\cdot \mathbf{\sigma })i{\sigma }_2,$$

$${\Delta }_{\pm }(\mathbf{k})={\Delta }_s\pm {\Delta }_t|{\mathbf{l}}_{\mathbf{k}}|,$$

$$t_{\pm }(\mathbf{k})=-e^{\mp i{\varphi }_l(\mathbf{k})},{\varphi }_l(\mathbf{k})={\tan }^{-1}(l_y/l_x).$$

$${\hat{g}}^{\text{R,A,M}}$$

$$\mathcal{N}$$

$$\mathcal{G}$$

$$(i{\mathbf{v}}_F\cdot {\nabla }_{\mathbf{R}}+2z)\gamma$$

$$(i{\mathbf{v}}_F\cdot {\nabla }_{\mathbf{R}}-2z)\tilde{\gamma }$$

$$\gamma (\rho )={\gamma }_h+e^{i{\Omega }_1\rho }{\delta }_0{\left(e^{i{\Omega }_2\rho }+C(\rho ){\delta }_0\right)}^{-1}$$

$$U_{\underline{\mathbf{k}}}^{†}\Gamma (\underline{\mathbf{k}},\epsilon )U_{-\underline{\mathbf{k}}}^{\ast }={\Gamma }^s(\underline{\mathbf{k}},\epsilon )={\gamma }^s(\mathbf{k},\epsilon )=U_{\mathbf{k}}^{†}\gamma (\mathbf{k},\epsilon )U_{-\mathbf{k}}^{\ast }$$

$$U_{-\mathbf{k}}^T\tilde{\Gamma }(\mathbf{k},\epsilon )U_{\mathbf{k}}={\tilde{\Gamma }}^s(\mathbf{k},\epsilon )={\tilde{\gamma }}^s(\underline{\mathbf{k}},\epsilon )=U_{-\underline{\mathbf{k}}}^T\tilde{\gamma }(\underline{\mathbf{k}},\epsilon )U_{\underline{\mathbf{k}}},$$

$$V_{s_1{s}_2{s}_3{s}_4}(\mathbf{k},{\mathbf{k}}^{\prime })$$

$$v_s^2+v_t^2+v_m^2=1$$

$$(v_s,v_t,v_m)$$

$$V(\mathbf{k},{\mathbf{k}}^{\prime })=\frac{V}{2}{\mathcal{Y}}_{\mathbf{k}}{\mathcal{Y}}_{{\mathbf{k}}^{\prime }}^{\ast }\left(\begin{array}{cc}{v}_s+v_t|{\mathbf{l}}_{\mathbf{k}}||{\mathbf{l}}_{{\mathbf{k}}^{\prime }}|-v_m{l}_{+}& v_s-v_t|{\mathbf{l}}_{\mathbf{k}}||{\mathbf{l}}_{{\mathbf{k}}^{\prime }}|-v_m{l}_{-}\\ v_s-v_t|{\mathbf{l}}_{\mathbf{k}}||{\mathbf{l}}_{{\mathbf{k}}^{\prime }}|+v_m{l}_{-}& v_s+v_t|{\mathbf{l}}_{\mathbf{k}}||{\mathbf{l}}_{{\mathbf{k}}^{\prime }}|+v_m{l}_{+}\end{array}\right)$$

$$⟨\bullet ⟩$$

$$\left(\begin{array}{c}{\Delta }_{+}(\mathbf{k})\\ {\Delta }_{-}(\mathbf{k})\end{array}\right)$$

$$\mathfrak{f}(\mathbf{k},{ϵ}_n)$$

$$\left(\begin{array}{c}{\Delta }_s\\ {\Delta }_t\end{array}\right)$$

$$A_{\mathbf{k}}$$

$$\left(\begin{array}{c}{\Delta }_s\\ {\Delta }_t\end{array}\right)$$

$${\left[T^{⟨L_{\mathbf{k}}⟩}\right]}_{ij}\equiv T^{⟨{\left[L_{\mathbf{k}}\right]}_{ij}⟩},$$

$$L_{\mathbf{k}}$$

$$\ln \left(\frac{2e^{\gamma }{ϵ}_c}{\pi T_c}\right)⟨L_{\mathbf{k}}⟩\left(\begin{array}{c}{\Delta }_s\\ {\Delta }_t\end{array}\right)=\frac{1}{N_{F}V}\left(\begin{array}{c}{\Delta }_s\\ {\Delta }_t\end{array}\right)$$

Full text

I Introduction

Non-centrosymmetric materials lack a center of inversion in their crystal lattice. They have attracted increasing attention in recent years due to the fact that spin-orbit interaction has a strong effect on their physical properties.Bauer and Sigrist (2012); Qi and Zhang (2011); Yip (2014); Ando and Fu (2015); Nakajima et al. (2015); Rößler et al. (2011); Puggioni and Rondinelli (2014) In crystals with a center of inversion the band-diagonal elements of the spin-orbit interaction in a Bloch basis, ${\bf L}_{nn}({\bf k})$, vanish by symmetry. This is not the case for non-centrosymmetric materials, where these diagonal elements can be non-zero and indeed large (30-300 meV).Samokhin (2009) Anderson, in discussing heavy fermion materials, used group classification to study the possibilities for spin-triplet superconductivity in spin-orbit coupled materials.Anderson (1984) Experimental signatures of spin-triplet (as well as spin-singlet) pairing were found in the non-centrosymmetric heavy-fermion superconductor CePt3Si, discovered in 2004.Bauer et al. (2004) Since then many more non-centrosymmetric superconductors (NCSs) have been identified, including Y2C3,Singh and Mazin (2004) Li2(Pd1-xPtx)3B,Togano et al. (2004) CeIrSi3,Sugitani et al. (2006) UIr,Akazawa et al. (2004), BiPd, Joshi et al. (2011), and PbTaSe2,Ali et al. (2014) amongst others.Kimura et al. (2005); Settai et al. (2007); Honda et al. (2010); Kneidinger et al. (2014); Shibayama et al. (2007); Klimczuk et al. (2006); Zuev et al. (2007); Kneidinger et al. (2015) These materials show signs of both spin-singlet and spin-triplet supconductivity to a varying degree. The system Li2(PdxPt1-x)3B has been studied in more detail,Badica et al. (2007); Takeya et al. (2009) indicating that the difference between the two end compounds, $x=0$ and $x=1$, can at least in part be explained by a dominating triplet component for $x=0$, i.e. Li2Pt3B, whereas Li2Pd3B seems to have a dominating $s$-wave singlet component, indicated by the rather low value of the upper critical magnetic field extrapolated to zero temperature. Some systems, like LaNiC2 and LaNiGa2, are candidates for a non-unitary spin-triplet pairing state.Hillier et al. (2009)

Furthermore, it has been shown that, as spin-orbit interaction is time-reversal invariant, these superconductors can be topologically non-trivial.Sato and Fujimoto (2009); Tanaka et al. (2009); Kitaev (2009); Schnyder et al. (2009, 2008); Qi and Zhang (2011); Ando and Fu (2015); Samokhin (2015); Scheurer and Schmalian (2015); Scheurer (2016) The topology and the singlet-triplet admixture are both a consequence of the spin-orbit coupling (SOC) term in the Hamiltonian of these materials, which is derived from the non-relativistic limit of the Dirac equation and is proportional to the gradient of the crystal lattice potential. The lack of a center of inversion in the unit cell allows the gradient of the potential to be large throughout the Brillouin zone (BZ), and thus the SOC cannot be neglected. The above-mentioned property, that diagonal elements of the SOC in a Bloch basis are in general non-vanishing in non-centrosymmetric materials, allows to study the effects of the SOC in a minimal one-band model, Samokhin (2009) which is not possible in centrosymmetric materials.

In this paper we theoretically study NCSs with the emphasis on self-consistent superconducting order parameters for various surface orientations, as well as for all the topological phases of the crystal point groups $T_{d}$, $C_{4v}$, and $O$ with a closed Fermi surface. The SOC vector is expanded in terms of harmonic functions, constrained by the symmetries of the point group, to second order. The relative weight of the first and second order terms is parameterized by $g_{2}$. Second order terms are investigated for the point groups $C_{4v}$ and $O$: one non-zero value of $g_{2}$ for $C_{4v}$ and three for $O$. Besides the gapped topologically trivial phases all point groups have one non-trivial gapless phase, and $O$ has, for a closed Fermi surface, four non-trivial gapped phases, and we have chosen values of $g_{2}$ to correspond to these phases. In the literature the point group $C_{4v}$ with $g_{2}=0$ has been studied extensively,Vorontsov et al. (2008); P. A. Frigeri and Sigrist (2006); Frigeri et al. (2004a, b); Eschrig et al. (2012) as well as $O$ with values of $g_{2}$ equivalent to our choice of $g_{2}=0.7$.Yuan et al. (2006); Klam et al. (2012); Brydon et al. (2011) All results we present in this paper are self-consistent, and for parameter combinations not discussed so far in the literature. Non-self-consistent results for surface spectra for various point groups and surface orientations were obtained in Ref. Eschrig et al., 2012, and subsequently in Ref. Brydon et al., 2011. Topological aspects were in the focus of attention in Ref. Schnyder et al., 2012, whereas in Ref. Timm et al., 2015 the possibility of a surface instability was suggested.

II Theory

II.1 Normal state band dispersion

Within an effective one-band model, the SOC term in the Hamiltonian is given by $H^{\text{SO}}_{\bf k}=\alpha{\bf l}_{\bf k}\cdot\bm{\sigma}$, where $\alpha$ is the SOC strength, $\bm{\sigma}=(\sigma_{1},\sigma_{2},\sigma_{3})$ is a vector of spin Pauli matrices, and ${\bf l}_{\bf k}$ is the SOC vector which is real, invariant under crystal point group operations $g$,

[TABLE]

and odd in ${\bf k}$, ${\bf l}_{-{\bf k}}=-{\bf l}_{\bf k}$. We normalize the SOC vector such that its maximum magnitude within the BZ is unity, $\max_{\bf k\in BZ}|{\bf l}_{\bf k}|=1$.

The kinetic part of the normal-state Hamiltonian can thus be written as

[TABLE]

with $\xi_{\bf k}=\epsilon_{\bf k}-\mu$, where $\epsilon_{\bf k}$ is the band dispersion in the absence of SOC (we will use for simplicity a nearest-neighbor tight-binding dispersion), $\mu$ is the chemical potential, and $c_{{\bf k}\alpha}$ ($c^{\dagger}_{{\bf k}\alpha}$) are fermion annihilation (creation) operators for a quasiparticle with spin $\alpha\in\left\{\uparrow,\downarrow\right\}$. We will study simple cubic (CUB) and body centered cubic (BCC) lattices. The corresponding nearest-neighbor tight binding dispersions are

[TABLE]

and

[TABLE]

where $t_{1}$ is the hopping integral.

The point groups considered here are the cubic point group $O$, relevant for e.g. Li2PdxPt3-x;Togano et al. (2004); Chandra et al. (2005); Lee and Pickett (2005); Badica et al. (2007); Takeya et al. (2009) the tetragonal point group $C_{4v}$, relevant for e.g. CePt3Si;Bauer et al. (2004) and the tetrahedral point group $T_{d}$, relevant for e.g. Y2C3.Yuan et al. (2011) We use dispersion (3) for the cubic point group, $O$, and for sake of simplicity also for the tetragonal point group, $C_{4v}$, whereas dispersion (4) will be used for the tetrahedral point group $T_{d}$. The SOC vectors are obtained as lattice Fourier series, ${\bf l}_{\bf k}=\sum_{n}{\bf l}_{n}\sin({\bf k}\cdot{\bf R}_{n})$, where ${\bf R}_{n}$ are Bravais lattice vectors, and where the invariance under point group operations, Eq. (1), leads to restrictions on the ${\bf l}_{n}$.Samokhin (2009)

The Hamiltonian, Eq. (2), is diagonalized and brought to the so-called helicity basis by the canonical transformation $U_{\bf k}\left({\bf l}_{\bf k}\cdot\bm{\sigma}\right)U^{\dagger}_{\bf k}=|{\bf l}_{\bf k}|\sigma_{3}$, where

[TABLE]

with $\phi_{l}=\tan^{-1}(l_{y}/l_{x})$ and $\theta_{l}=\tan^{-1}(\sqrt{l^{2}_{x}+l^{2}_{y}}/l_{z})$ being the spherical angles of the SOC vector, ${\bf l}_{\bf k}=(l_{x},l_{y},l_{z})^{T}$, yielding

[TABLE]

where the helical index takes the values $\lambda=\{+,-\}$, and the helical band dispersion is given by $\xi^{\pm}_{\bf k}=\xi_{\bf k}\pm\alpha|{\bf l}_{\bf k}|$. Note that $\xi^{\lambda}_{\bf k}=\xi^{\lambda}_{-{\bf k}}$ even though the SOC vector is antisymmetric. This is a consequence of Eq. (2) being time-reversal invariant. Furthermore, the quasiparticle spin is fixed with respect to its momentum on each band, being parallel ($\lambda=+$) or antiparallel ($\lambda=-$) to ${\bf l}_{\bf k}$.

II.2 Superconducting state

Superconductivity is modeled within the Nambu-Gor’kov formalism. Under the canonical transformation defined above the Nambu spinor $\hat{C}_{\bf k}=(c_{{\bf k}\uparrow},c_{{\bf k}\downarrow},c^{\dagger}_{-{\bf k}\uparrow},c^{\dagger}_{-{\bf k}\downarrow})^{T}$ transforms into its helical equivalent $\hat{B}_{\bf k}=(b_{{\bf k}+},b_{{\bf k}-},b^{\dagger}_{-{\bf k}+},b^{\dagger}_{-{\bf k}-})^{T}\equiv\hat{U}_{\bf k}\hat{C}_{\bf k}$ with $\hat{U}_{\bf k}\equiv\text{diag}(U_{\bf k},U^{*}_{-\bf k})$, and the ”hat” denotes Nambu structure. It is straightforward to construct $4\times 4$ helical Green functions, e.g. the retarded $\hat{G}^{\text{R}}_{{\bf k}_{1}{\bf k}_{2}}(t_{1},t_{2})=-i\Theta(t_{1}-t_{2})\langle\{\hat{B}_{{\bf k}_{1}}(t_{1}),\hat{B}^{\dagger}_{{\bf k}_{2}}(t_{2})\}\rangle_{\mathcal{H}}$, where $\Theta$ is the Heaviside step function, $\langle\bullet\rangle_{\mathcal{H}}$ denotes a grand canonical average, $\{\bullet,\bullet\}$ is an anticommutator, and $\hat{B}_{\bf k}(t)$ a Heisenberg operator. The quasiclassical propagator is obtained by integrating out fast oscillations from the full Green functions. In the case when the magnitude of the SOC is much smaller than the Fermi energy, $\alpha\ll E_{F}$, it suffices to integrate over $\xi_{\bf k}$ and treat the SOC term perturbatively. For this case, in Wigner coordinates the quasiclassical propagator is given by $\check{g}({\bf k}_{F},{\bf R},\epsilon,t)=\int\text{d}\xi_{{\bf k}}\hat{\tau}_{3}\check{G}({\bf k},{\bf R},\epsilon,t)$, with ${\bf k}$ parameterized by $(\xi_{\bf k},{\bf k}_{F})$, $\xi_{{\bf k}}={\bf v}_{F}\cdot({\bf k}-{\bf k}_{F})$, $\hat{\bm{\tau}}=(\hat{\tau}_{1},\hat{\tau}_{2},\hat{\tau}_{3})$ are Pauli matrices in particle-hole space, and the ”check” denotes Keldysh matrix structure.

The SOC term enters the transport equations as a source term. Within this approximation the Eilenberger equation Eilenberger (1968) for the quasiclassical Green function takes the following form in the helicity basis

[TABLE]

with $z=i\epsilon_{n}=i\pi T(2n+1)$ for Matsubara, and $z=\epsilon\pm i0^{+}$ for retarded (advanced), quantities. $[\bullet,\bullet]$ is a commutator, the SOC term is $\hat{v}_{\text{SO}}=\alpha|{\bf l}_{{\bf k}_{F}}|\sigma_{3}\hat{\tau}_{0}$, and the gap has the form

[TABLE]

where the ”tilde operation” is defined as the particle-hole conjugate, $\tilde{Q}({\bf k}_{F},{\bf R},z,t)\equiv Q^{*}(-{\bf k}_{F},{\bf R},-z^{*},t)$. Eq. (7) is supplemented by the normalization condition $(\hat{g}^{\text{R,A,M}})^{2}=-\pi^{2}\hat{1}$. In order to simplify notation, we will henceforth drop the subscript $F$ at the Fermi momentum; all momenta in the quasiclassical theory are Fermi momenta. The subscript will be written out only when it is necessary to avoid confusion. We consider time-independent situations, such that the time variable $t$ will be dropped from here on.

The lack of a center of inversion allows for an admixture of spin-singlet and spin-triplet pairing.P. A. Frigeri and Sigrist (2006) The spin-triplet vector is set to be parallel to the SOC vector in order to maximize $T_{c}$.Frigeri et al. (2004a) In spin basis the order parameter is written

[TABLE]

where $\mathcal{Y}_{\bf k}$ is a crystal basis function corresponding to irreducible representation of the dominant pairing channel, and $\Delta_{s}$ and $\Delta_{t}$ are referred to as the singlet and triplet component, respectively. In the helicity basis the order parameter takes the form $\Delta({\bf k})=\mathcal{Y}_{\bf k}\cdot\text{diag}(\Delta_{+}({\bf k})t_{+}({\bf k}),\Delta_{-}({\bf k})t_{-}({\bf k}))$, where

[TABLE]

and the phase factors are given by

[TABLE]

Note that $t_{\pm}(-{\bf k})=-t_{\pm}({\bf k})$.

Eq. (7) can be parameterized in terms of coherence functions, $\gamma({\bf k},{\bf R},z)$ and $\tilde{\gamma}({\bf k},{\bf R},z)$,Eschrig (2000) in such a way as to automatically fulfill the normalization condition,

[TABLE]

where the top (bottom) sign corresponds to $\hat{g}^{\text{R}}$ ($\hat{g}^{\text{A}}$), and in the case of $\hat{g}^{\text{M}}$, to positive (negative) Matsubara frequencies. With this, Eq. (7) transforms into two decoupled Riccati differential equations,

[TABLE]

In the homogeneous case, i.e. in the bulk, the solution is $\gamma_{h}=\mathcal{Y}_{\bf k}\cdot\text{diag}(\gamma_{+}({\bf k})t_{+}({\bf k}),\gamma_{-}({\bf k})t_{-}({\bf k}))$ with the abbreviations $\gamma_{\pm}=-\Delta_{\pm}/(z+i\sqrt{|\mathcal{Y}_{\bf k}\Delta_{\pm}|^{2}-z^{2}})$. For this case the SOC term drops out.

The surface problem is treated by solving Eqs. (13)-(14) along classical trajectories parallel to ${\bf v}_{F}$, using the homogeneous solutions as initial conditions at a sufficient distance from the surface. This is done by discretizing the path and treating the order parameter as a series of step functions in the middle between the desired grid points. Each step is solved analytically.Eschrig (2009) Parameterizing the path as ${\bf R}={\bf R}_{0}+\rho{\bf v}_{F}$ and writing the order parameter $\Delta(\rho)=\Delta_{0}+\Theta(\rho)(\Delta_{1}-\Delta_{0})$ at one of these steps, $\gamma(\rho)$ with $\rho>0$ is given by

[TABLE]

with $\delta_{0}=[\gamma_{0}-\gamma_{h}]$, where $\gamma_{0}\equiv\gamma(0)$ is the initial value and $\gamma_{h}$ is the homogeneous solution for $\rho>0$, $\Omega_{1}=z-\gamma_{h}\tilde{\Delta}$ and $\Omega_{2}=-z+\tilde{\Delta}\gamma_{h}$, and $C(\rho)=C_{0}e^{i\Omega_{1}}-e^{i\Omega_{2}}C_{0}$, where $C_{0}$ is the solution to $C_{0}\Omega_{1}-\Omega_{2}C_{0}=\tilde{\Delta}$. The solution for $\tilde{\gamma}(\rho)$ is completely analogous.

The reflection at the surface is in leading approximation (as $\alpha\ll E_{F}$) considered to be specular in spin space, with the momentum component parallel to the surface, ${\bf k}_{\parallel}$, conserved. Writing the momentum for incoming trajectories ${\bf k}=(k_{\perp},{\bf k}_{\parallel})$ this gives the momentum for outgoing trajectories as $\underline{\bf k}=(-k_{\perp},{\bf k}_{\parallel})$. Following Ref. Eschrig, 2000, incoming (outgoing) quantities are written with lowercase (uppercase) symbols and the surface boundary conditions become

[TABLE]

and

[TABLE]

where the $s$ superscript indicates that the coherence functions are expressed in the spin basis.

II.3 Gap equation

The pairing potential in spin space can be written as a sum of singlet, triplet, and a mixture term P. A. Frigeri and Sigrist (2006)

[TABLE]

where $v_{s}$, $v_{t}$, and $v_{m}$ are free parameters that describe the relative coupling strength of each term, respectively, $V$ is the overall pairing potential strength, and $\mathcal{Y}_{\bf k}$ is the basis function of the irreducible representation with the highest $T_{\text{c}}$. To avoid ambiguity, we normalize the relative pairing strengths according to

[TABLE]

and for later reference introduce spherical coordinates

[TABLE]

In helicity space the pairing potential takes the form

[TABLE]

with $l_{\pm}=|{\bf l}_{\bf k}|\pm|{\bf l}_{{\bf k}^{\prime}}|$. The self-consistency equation in the Matsubara formalism is expressed in terms of Fermi surface averages $\langle\bullet\rangle$, defined as

[TABLE]

With this, the self-consistency equation takes the form

[TABLE]

where $\mathfrak{f}_{\pm}$ are defined by

[TABLE]

the phase factors are defined in Eq. (11), and $\epsilon_{c}$ is the BCS technical cutoff. Using the relations $\Delta_{s}=\frac{1}{2}(\Delta_{+}+\Delta_{-})$ and $\Delta_{t}|{\bf l}_{\bf k}|=\frac{1}{2}(\Delta_{+}-\Delta_{-})$, the implicit form of the self-consistency equation for the singlet and triplet components of the order parameter reads

[TABLE]

where

[TABLE]

After elimination of the cutoff and the pairing strength $V$ in favor of the superconducting transition temperature, one obtains

[TABLE]

where the matrix exponent in the logarithm, $T^{\left\langle L_{\bf k}\right\rangle}$, is taken element-wise, i.e.

[TABLE]

and with

[TABLE]

Furthermore, $\lambda_{\text{max}}\equiv\max\{\lambda_{1},\lambda_{2}\}$, and $\lambda_{1,2}$ are the eigenvalues of the matrix $\left\langle L_{\bf k}\right\rangle$. We follow Ref. Grein et al., 2013 in eliminating the cut-off dependence in close vicinity to the surface as well. For details on the numerical procedure to achieve self-consistency see appendix A.

II.4 Bulk superconducting phase

At $T=T_{c}$ the self-consistency equation reduces to

[TABLE]

where $\gamma=0.5772...$ is the Euler-Mascheroni constant.

The number of positive eigenvalues of $\langle L_{\bf k}\rangle$ determines the number of nucleation channels, with $T_{c}$ determined by the largest eigenvalue $\lambda_{\rm max}$. Using Eq. (LABEL:phiv), the eigenvalues can be mapped onto the unit sphere. How the number of nucleation channels depends on the spherical angles $\phi_{v}=\tan^{-1}(v_{t}/v_{s})$ and $\theta_{v}=\tan^{-1}(\sqrt{v^{2}_{s}+v^{2}_{t}}/v_{m})$ can be seen in Fig. 1. When both eigenvalues are positive there are two possible nucleation channels, the dominant and the subdominant one. The dominant channel is responsible for the transition to superconductivity due to its larger critical temperature. The dominant channel also determines the singlet to triplet order parameter ratio, $\Delta_{s}/\Delta_{t}$, and their relative sign. The subdominant channel nucleates at a lower temperature $T^{\text{sub.}}_{c}\leq T_{c}$. With a finite mixing, $v_{m}\not=0$, an admixture of singlet and triplet components is obtained. For certain choices of the parameters $(v_{s},v_{t},v_{m})$ it is possible to achieve a cross-over from dominating singlet component at $T=T_{c}$ to a dominating triplet component at $T=0$. An example in the single-channel region is $(v_{s},v_{t},v_{m})=(1,0,a/(\langle|{\bf l}_{{\bf k}}|^{2}\rangle-a^{2}))$ (ignoring normalization) with the parameter $a$ being slightly larger than the maximum value of the SOC vector on the Fermi surface, e.g. $a=1.01\max|{\bf l}_{{\bf k}}|$. This means that the topology of the system can be sensitive to its sub-critical temperature.

In certain parameter ranges for $(v_{s},v_{t},v_{m})$ it is possible to construct a configuration with two active channels in which the subdominant channel has a lower free energy at $T=0$, thus inducing a second phase transition below $T_{c}$. The simplest way to get a second phase transition is to choose $(v_{s},v_{t},v_{m})$ in such a way as to get a dominant channel with a large triplet component, as well as a rather large subdominant critical temperature. A example for such a choice is $(v_{s},v_{t},v_{m})=(0.999\langle|{\bf l}_{{\bf k}}|^{2}\rangle,1,0)$ (ignoring normalization) giving a dominant pure triplet channel, and a subdominant pure singlet channel. The subdominant critical temperature is $T^{\text{sub.}}_{c}=0.996T_{c}$ for the point groups and SOC vectors in table 1 (assuming $\mathcal{Y}_{\bf k}=1$). The condensation energy at zero temperature, assuming the same density of states on both Fermi surface sheets (which is the approximation employed here as the splitting is small), is given by

[TABLE]

and yields a lower value for the subdominant channel, for all three point groups and SOC vectors considered in this work, with this choice of $(v_{s},v_{t},v_{m})$.

II.5 Angle-resolved density of states

The angle-resolved surface density of states (DOS) is given by $N({\bf k},\epsilon)=-(2\pi)^{-1}N_{F}\text{Im}\text{Tr}_{\lambda}[\mathfrak{g}^{\text{R}}({\bf k},\epsilon)]$, or explicitly in terms of coherence functions

[TABLE]

The spin-resolved DOS along the quantization axis $j\in\{x,y,z\}$ is given by $N^{(j)}_{\pm}({\bf k},\epsilon)=N({\bf k},\epsilon)\pm N^{(j)}({\bf k},\epsilon)$, where

[TABLE]

Note that all quantities in Eq. (33) are expressed in the spin basis. Using a non-self-consistent order parameter, with the bulk solution all the way to the surface, it is straightforward to show that there are two classes of trajectories giving rise to Andreev bound states (ABS) at zero energy (see appendix B for details). Introducing the notation $\Upsilon_{\bf k}\equiv\text{sign}[\mathcal{Y}_{\bf k}\Delta_{-}({\bf k})]$ the first class of trajectories is simply given by $\Upsilon_{\underline{\bf k}}=-\Upsilon_{\bf k}\not=0$. With the spherical angles $(\phi_{l},\theta_{l})$ and $(\phi_{\underline{l}},\theta_{\underline{l}})$ corresponding to ${\bf l}_{\bf k}$ and ${\bf l}_{\underline{\bf k}}$ respectively, the second class is given by solutions to

[TABLE]

with the definition $F(\phi_{l},\theta_{l},\phi_{\underline{l}},\theta_{\underline{l}})=\cos(\theta_{l})\cos(\theta_{\underline{l}})+\cos(\phi_{l}-\phi_{\underline{l}})\sin(\theta_{l})\sin(\theta_{\underline{l}})$, provided that $(\Upsilon_{\bf k},\Upsilon_{\underline{\bf k}})=(0,-1)$, $(\Upsilon_{\bf k},\Upsilon_{\underline{\bf k}})=(-1,0)$, or $(\Upsilon_{\bf k},\Upsilon_{\underline{\bf k}})=(-1,-1)$. This second class of bound states arises due to the phase factors $t_{\pm}({\bf k})$ defined in Eq. (11), which can yield an extra phase shift of $\pi$. These results remain true for self-consistent order parameters as long as the gap does not completely close at some distance from the surface.

II.6 Point contact spectra

The point contact conductance between a normal metal and an NCS is computed using the following assumptions: the size of the point contact is much smaller than the coherence length but much larger than the Fermi wavelength, the Fermi surfaces on both sides of the interface are considered to be equal, and the proximity effect is ignored. The normal metal having index 1, and the NCS index 2, the scattering matrix of the interface in the spin/helicity basis is given by

[TABLE]

with the transmission amplitude

[TABLE]

where $t_{0}$ is the tunneling parameter and $\alpha_{\bf k}$ is the angle between the surface normal and the Fermi velocity of the outgoing trajectories in the normal metal. The reflection amplitude is given by $r=\sqrt{1-t^{2}}$. The zero-temperature tunnel conductance is given by Eschrig et al. (2012)

[TABLE]

where the expression is evaluated at $\epsilon=eV$, ${\bf v}_{{F}_{1}}$ is the Fermi velocity in the normal metal, $\langle\bullet\rangle_{\text{out}}$ indicates that the average is only for outgoing trajectories in the normal metal, $B(\epsilon)=S_{12}\left(\sigma_{0}+A_{2}(\epsilon)S_{22}\right)$,

[TABLE]

and $||\bullet||^{2}\equiv\frac{1}{2}\text{Tr}\left[(\bullet)(\bullet)^{\dagger}\right]$. The normal state conductance, $G_{N}$, is simply obtained by setting the coherence functions to zero.

II.7 Topology

We characterize the topology of a system by computing three topological invariants. The starting point is the Bogolioubov-de Gennes (BdG) Hamiltonian

[TABLE]

obeying time-reversal symmetry, $\mathcal{T}$, particle-hole symmetry, $\mathcal{C}$ , as well as the combined ’chiral’ symmetry $\mathcal{S}=i\mathcal{TC}$. The BdG Hamiltonian is thus of the symmetry class DIII.Schnyder and Brydon (2015) It anticommutes with $\mathcal{S}$ and in the basis where $\mathcal{S}$ is block diagonal $H$ becomes block off-diagonal, $\bar{H}=VHV^{\dagger}$. The flat-band block off-diagonal Hamiltonian $Q({\bf k})$ is constructed by projecting all bands above (below) the gap to $+1$ ($-1$)

[TABLE]

where $q({\bf k})$ is a $2\times 2$ matrix in the one-band model (we set for simplicity $\mathcal{Y}_{\bf k}=1$)

[TABLE]

with $A=\alpha+i\Delta_{t}$, $B_{\bf k}=\xi_{\bf k}+i\Delta_{s}$, $\lambda_{1}=\lambda^{-1}_{+}-\lambda^{-1}_{-}$, $\lambda_{2}=\lambda^{-1}_{+}+\lambda^{-1}_{-}$, where $\lambda_{\pm}=\left|A|{\bf l}_{\bf k}|\pm B_{\bf k}\right|$. Note that $Q({\bf k})$, and thus $q({\bf k})$, is ill-defined for nodal order parameters.

Fully gapped systems are classified by calculating the 3D winding number which is defined as

[TABLE]

where Einstein summation is implied, $\varepsilon^{abc}$ is the Levi-Civita pseudo-tensor, $a,b,c\in\{k_{x},k_{y},k_{z}\}$, and the integral is over the entire first BZ. From the definition of $q$ it is clear that $\nu$ is only well-defined if the order parameter on the negative helical Fermi surface does not have nodes, i.e. $\Delta_{-}({\bf k}^{-}_{F})\not=0$. There are two ways this can be true; either $\text{sign}[\Delta_{-}({\bf k}^{-}_{F})]=+1\;\forall{\bf k}^{-}_{F}\Longrightarrow\Delta_{s}/\Delta_{t}>\max|{\bf l}({\bf k}^{-}_{F})|$, or $\text{sign}[\Delta_{-}({\bf k}^{-}_{F})]=-1\;\forall{\bf k}^{-}_{F}\Longrightarrow\Delta_{s}/\Delta_{t}<\min|{\bf l}({\bf k}^{-}_{F})|$. We calculate $\nu$ numerically using the procedure in appendix C.

Nodal systems are classified by calculating the 1D winding number which is defined as

[TABLE]

where $l$ parameterizes the loop $\mathcal{L}$ in the BZ, and $\nabla_{l}$ is the directional gradient along this loop. The loop $\mathcal{L}$ cannot pass through nodes of the order parameter, but is other than that arbitrary. The 1D Hamiltonian for this loop is in general not time-reversal invariant and is thus of symmetry class AIII.Schnyder and Brydon (2015) In order to characterize a nodal phase the loop needs to be constructed in such a way as to always encircle a line node of $\Delta_{-}({\bf k}^{-}_{F})$ for any Fermi surface geometry.

With increasing singlet to triplet ratio the first nodes appear at the points where $\Delta_{s}/\Delta_{t}=\min|{\bf l}({\bf k}^{-}_{F})|$. Increasing $\Delta_{s}/\Delta_{t}$ further the nodal rings continue to be positioned around these points until they connect with one another. At this stage the nodal rings become positioned around the points where they eventually disappear $\Delta_{s}/\Delta_{t}=\max|{\bf l}({\bf k}^{-}_{F})|$. Thus a general loop should pass through the points where the nodal rings appear and disappear. This is accomplished by the loop

[TABLE]

where $\partial\text{BZ}$ is the BZ boundary, and the arrows do not necessarily imply straight lines.

In order to study the topology’s effect on the surface states the 1D winding number is also computed for straight noncontractible loops, i.e. loops traversing one or several of the three circles making up the BZ torus ${\bf T}^{3}=S^{1}\times S^{1}\times S^{1}$, that are perpendicular to the surface. Writing the momentum ${\bf k}=({\bf k}_{\parallel},k_{\perp})$ and the surface normal ${\bf n}=(l,m,n)$ the 1D winding number is written

[TABLE]

Restricting ourselves to time-reversal invariant noncontractible loops another topological invariant can be defined. Namely the $\mathbb{Z}_{2}$ invariant

[TABLE]

where ${\bf K}$ are time-reversal invariant momenta on the loop, and $\text{Pf}[\bullet]$ denotes the Pfaffian of an antisymmetric matrix $\bullet$. The 1D Hamiltonian for this loop is of the symmetry class DIII.Schnyder and Brydon (2015)

The singlet (triplet) component is said to be dominant if the inequality $\Delta_{s}/\Delta_{t}>\max|{\bf l}({\bf k}^{-}_{F})|$ is true (false). With a dominant singlet component the material is fully gapped. Increasing $\Delta_{s}$ and/or decreasing $\Delta_{t}$ the material becomes nodal and eventually fully gapped again if $\min|{\bf l}({\bf k}^{-}_{F})|>0$. As is shown below the dominance of either component is temperature dependent.

II.8 Surface band structure

The surface band structure is computed by first Fourier transforming the BdG Hamiltonian in the relative momentum coordinate $k_{\perp}$ in the direction of the surface normal ${\bf n}$,

[TABLE]

The helical dispersion, $\xi^{\lambda}_{\bf k}$, contains for the tight-binding approximation we use trigonometric functions whose Fourier transform give rise to a series of delta functions

[TABLE]

where $H_{-j}({\bf k}_{\parallel},{\bf R}-\frac{1}{2}\rho{\bf n})=H^{\dagger}_{j}({\bf k}_{\parallel},{\bf R}+\frac{1}{2}\rho{\bf n})$, $j$ is a layer index, and $\rho_{0}$ is the length one needs to move along the direction of the surface normal in order to return to a translation-equivalent point in the lattice unit cell. The sum has a finite number of terms, i.e. there exist a number $j_{c}$ such that $H_{j}=0\;:\;|j|>j_{c}$. The terms $H_{j}$ with $j\neq 0$ can be interpreted in terms of hopping across the layers. Discretizing the center-of-mass coordinate ${\bf R}$ in steps of $\rho_{0}$, the Schrödinger equation for $L$ layers can be written

[TABLE]

where $l=0,1,\dots,L-1$ and $j(l)=\min\{j_{c},l\}$ which takes care of the boundary conditions, i.e. no hopping across the boundary. Eq. (49) can be written more compactly as a matrix equation $H_{\text{eff}}({\bf k}_{\parallel})\psi({\bf k}_{\parallel})=E({\bf k}_{\parallel})\psi({\bf k}_{\parallel})$, and the band structure is given by the eigenvalues of $H_{\text{eff}}$, Non-trivial topology gives rise to zero-energy ABS. We are therefore mainly interested in the band structure close to zero energy. This allows us to avoid diagonalizing $H_{\text{eff}}$, and instead only compute the smallest magnitude eigenvalues using the Lanczos method. Note that the order parameter is suppressed at both surfaces.

III Numerical Results

In this work the SOC strength entering the quasiclassical calculations is considered to be much smaller than the Fermi energy, $\alpha\ll E_{F}$. In this case the Fermi surface is only weakly split. Ignoring this splitting, and the Fermi velocity renormalisation, the quasiparticles with opposite helicity are assigned to a single, common Fermi surface, and move coherently along classical trajectories. In addition, for the quasiclassical part of the numerical calculations, the Fermi surface is approximated to be spherical, with $|{\bf k}_{F}|$ being equal to the average of the Fermi surface defined by $\xi\left({\bf k}_{F}\right)=0$ with $(t_{1},\mu)=(-40\alpha,-50\alpha)$. Here, $t_{1}$ determines the bandwidth, which must be large compared to the Fermi-surface splitting in order for the approximation of equal Fermi surfaces for both helicities to be valid, and $\mu<0$ must smaller than $t_{1}-\alpha$ in order for the Fermi surface to be closed. The chosen values are consistent with the approximation of an approximately spherical Fermi surface. The SOC term enters the transport equations as a source term. In the following, we restrict our discussion to the maximally symmetric basis function corresponding to the irreducible representation $A_{1}$, i.e. $\mathcal{Y}_{\bf k}=1$.

III.1 The Cubic Point Group $O$

To next-nearest neighbors in the sum over Bravais lattice sites Samokhin (2009) the SOC vector corresponding to the cubic point group $O$, takes the form

[TABLE]

where $g_{2}$ is a free parameter which determines the relative weight between the first and second order contributions. Its magnitude and direction is illustrated in Fig. 2.

An important property of the SOC vector corresponding to the cubic point group is its lack of line nodes in the BZ, it only vanishes at specific points. With $g_{2}=0$ these points are simply $\Gamma$, X, M, and R [for the notation see Fig. 3(a)]. A finite value of $g_{2}$ brings about two more points. With $g_{2}>0$ they are positioned somewhere on the paths $\Gamma\rightarrow\text{R}$, and $\Gamma\rightarrow\text{M}$, and with $g_{2}<0$ on $\Gamma\rightarrow\text{R}$, and $\text{X}\rightarrow\text{R}$, in Fig. 3 (a). The exact positions, ${\bf k}^{*}$, of these points depend on the value of $g_{2}$, and are given by

[TABLE]

The lack of line nodes means that it is easy to construct a Fermi surface for which the minimum value of the SOC on the negative helical FS, $\min|{\bf l}({\bf k}^{-}_{F})|$, is not zero. The dependence of $\min|{\bf l}({\bf k}^{-}_{F})|$ on the chemical potential and the SOC parameter $g_{2}$ is shown in Fig. 3 (b). The SOC minimum is zero along certain lines in this parameter space. The line at $\mu=t_{1}$ marks the transition between open and closed FS, i.e. the FS is tangent to the X-point in the BZ. These lines in Fig. 3 (b) also mark the boundaries of fully gapped regions with different values of the 3D winding number $\nu$. This is demonstrated in figs. 3 (c) and (d) in which the topological phase diagram is shown for an open, $\mu=-20\alpha$, and closed, $\mu=-50\alpha$, Fermi surface respectively. White indicates that the system is fully gapped and topologically trivial, $\nu=0$, whereas the colored regions (excluding grey) indicate that the system is fully gapped and topologically non-trivial, $\nu\not=0$. Grey indicates a topologically non-trivial nodal phase, $N_{\mathcal{L}}=1$, with loops defined by Eq. (44).

The self-consistent order parameter is calculated for four different values of $g_{2}$, namely $g_{2}\in\{0,0.7,1.03,2.5\}$, one for each distinct gapped topologically non-trivial phase with a closed Fermi surface, i.e. the colored regions in Fig. 3 (d). This is done for nine values of the scaled bulk singlet to triplet ratio, $r^{\text{bulk}}_{\Delta}\in[0,1.1]$, with one active channel. These values are shown in table 1.

In order to investigate how the order parameter suppression depends on the surface orientation the order parameter is computed for a range of different surface normals, tracing out the path ${\bf n}=(1,0,0)\rightarrow(1,1,0)\rightarrow(1,1,1)\rightarrow(0,1,2)\rightarrow(1,0,0)$. As a measure of the suppression the ratio $r^{\text{bulk}}_{\Delta}/r^{\text{surf.}}_{\Delta}$, where $r^{\text{surf.}}_{\Delta}\equiv\Delta^{\text{surf.}}_{s}/(\Delta^{\text{surf.}}_{t}\max|{\bf l}_{{\bf k}_{F}}|)$ is the scaled surface singlet to triplet ratio, is plotted in Fig. 4 (a) - (d). The suppression is seen to be the largest for the surface normal ${\bf n}=(1,1,1)$

Along the same path of surface orientations the zero-bias conductance, computed with Eq. (37), is plotted in Fig. 4 (e) - (h). For singlet to triplet ratios in the interval $\min|{\bf l}_{{\bf k}_{F}}|<\Delta_{s}/\Delta_{t}<\max|{\bf l}_{{\bf k}_{F}}|$ very large zero-bias conductance is seen for all surface orientations except the two high symmetry axes, ${\bf n}=(1,0,0)$ and ${\bf n}=(1,1,0)$. This is due to there being no trajectories for which $\Delta_{-}$ changes sign upon reflection for these surface orientations. For all other surface orientations this is not the case, including the high symmetry axis ${\bf n}=(1,1,1)$. Note that all lines for which $\Delta_{s}/\Delta_{t}<\min|{\bf l}_{{\bf k}_{F}}|$ are degenerate, and the zero-bias conductance is zero for $\Delta_{s}/\Delta_{t}>\max|{\bf l}_{{\bf k}_{F}}|$. Furthermore, the surface suppression due to self-consistency does not affect the zero-bias conductance. This reflects the fact that the gap does not go to zero at some distance inwards from the surface for the obtained gap profiles.

The Andreev bound states (ABS) of NCSs have intricate structures and are spin polarized.Vorontsov et al. (2008) This is a consequence of the SOC being antisymmetric, $l_{\bf k}=-l_{-{\bf k}}$. States corresponding to different Andreev bound state branches have opposite spin polarization, and this spin polarization changes sign for reversed trajectories. As a result, the Andreev states carry spin current along the interface.Vorontsov et al. (2008) The existence of a surface spin current is a direct consequence of the spin-orbit coupling in the system.

As an example, the momentum angle-resolved and spin-resolved local density of states, $N^{(z)}(\phi,\epsilon)$, computed with Eq. (33), is plotted in Fig. 5(a) for momenta in the $xy$-plane (parameterized by the azimuthal angle $\phi$, the polar angle is $\theta=\pi/2$), at the surface with surface normal ${\bf n}=(1,0,0)$, for $g_{2}=1.03$ and a self-consistent pure triplet order parameter. An energy broadening $\epsilon\to\epsilon+i\delta$ with $\delta=10^{-2}$ was used, and the self-consistent order parameter was computed at $T=0.2T_{c}$. Red (blue) indicate relative polarization for spin up (down) quasiparticles. The spin polarization axis is along the $z$-axis and $N^{(x)}=N^{(y)}=0$. This is true for all values of $g_{2}$ with a pure triplet order parameter. However, the ABS structure is very different for the four $g_{2}$ values. Furthermore the spin polarization axis is found to be dependent on the singlet to triplet ratio, in addition to surface orientation.

The momentum-resolved zero-energy ABS for ${\bf n}=(1,1,1)$ are shown in Fig. 5 (b), computed with the bulk value of the order parameter all the way to the surface, assuming $r^{\text{bulk}}_{\Delta}=0.67$. The tunneling parameter was set to $t_{0}=10^{-\frac{1}{2}}$ (or $t_{0}^{2}=0.1$, making sure to be in the tunneling regime), and the broadening of the energies, $\epsilon\leftarrow\epsilon+i\delta$, with $\delta=10^{-3}$. The disk is the projection of the spherical Fermi surface onto the slab surface. Black indicates that there are no ABS for those momenta, green indicates ABS for which $(\Upsilon_{\bf k},\Upsilon_{\underline{\bf k}})=(+1,-1)$, and yellow $(\Upsilon_{\bf k},\Upsilon_{\underline{\bf k}})=(-1,+1)$. For this choice of surface orientation and singlet to triplet ratio these two types of trajectories are the only ones yielding ABS. This is not the case for lower singlet to triplet ratios, other $g_{2}$ values, and/or other surface orientations. Then there can exist solutions to Eq. (34). Indeed, for $\Delta_{s}/\Delta_{t}<\min|{\bf l}_{{\bf k}_{F}}|$ they are the only solutions yielding ABS. For $\Delta_{s}/\Delta_{t}>\max|{\bf l}_{{\bf k}_{F}}|$ no zero-energy ABS are seen.

Point contact conductance spectra for $g_{2}=1.03$, $t_{0}=10^{-\frac{1}{2}}$ and ${\bf n}=(1,1,1)$ are shown in Fig. 5 (c). A small energy broadening $\epsilon\to\epsilon+i\delta$ with $\delta=10^{-2}$ was used for the plot, except close to zero energy, where $\delta=10^{-5}$ was used (and 2.5 times as many momentum directions in the momentum average) in order to show the sharp zero-bias conductance peak. The transmission parameter is set to $t_{0}=10^{-\frac{1}{2}}$. Furthermore, $\max(\Delta)\equiv\Delta^{\text{bulk}}_{s}+\Delta^{\text{bulk}}_{t}\max(|{\bf l}_{{\bf k}_{\text{F}}}|)$, and the plots are shifted $0.2$ upwards from each other for the sake of visibility. The order parameters used are computed self-consistently at $T=0.2T_{\text{c}}$, and with only one active channel. The point contact conductance spectra differ widely between surface orientations and the values of $g_{2}$, in addition to the less pronounced difference between singlet to triplet ratios. The most striking difference is the appearance of zero-bias conductance peaks (ZBCPs) which are present for all singlet to triplet ratios in the interval $\min|{\bf l}_{{\bf k}_{F}}|<\Delta_{s}/\Delta_{t}<\max|{\bf l}_{{\bf k}_{F}}|$ provided there are trajectories with $\text{sign}[\Delta_{-}({\bf k})]=-\text{sign}[\Delta_{-}(\underline{\bf k})]$.

In Fig. 5 (d) the topological invariants $N_{(111)}$ and $W_{(111)}$ are plotted for $r^{\text{bulk}}_{\Delta}=0.67$. However, $W_{(111)}=1$, i.e. trivial, for this choice of parameters, and trivial topology is colored white. Light green/blue corresponds to $N_{(111)}=\pm 1$. The dashed circle indicates the projection of the spherical Fermi surface used in the quasiclassical calculations, i.e. Fig. 5 (a) - (c). Even though the Fermi surface is not spherical, it is clear that the zero-energy ABS are directly related to the topology. As is shown for the tetragonal point group $C_{4v}$ below, the ABS given by solutions to Eq. (34), for the relevant values of $(\Upsilon_{\bf k},\Upsilon_{\underline{\bf k}})$, is directly related to the $\mathbb{Z}_{2}$ invariant being non-trivial (i.e. $W_{(111)}=-1$).

Zero-energy states are present in the band structure whenever the aforementioned topological invariants have non-trivial values. In Fig. 5 (e) the surface band structure is shown for $r^{\text{bulk}}_{\Delta}=0.67$ along the $k^{\parallel}_{2}$-axis with $k^{\parallel}_{1}=0$, and $L=1.3\cdot 10^{4}$ layers. $N_{(111)}\not=0$ gives rise to singly degenerate zero-energy flat bands, one on each surface, with the corresponding wavefunctions decaying exponentially into the bulk. The surface momenta of the zero-energy flat bands are given by $N_{(111)}({\bf k}_{\parallel})\not=0$, which can be seen in Fig. 5 (f) where the lowest positive eigenvalue of $H_{\rm eff}$ [see Eq. (49)] is plotted for self-consistent order parameter. Note that the zero-energy flat-bands are given by the projection of non-trivial values of the 1D winding number.

III.2 The Tetragonal Point Group $C_{4v}$

To next-nearest neighbors in the sum over Bravais lattice sites Samokhin (2009) the SOC vector corresponding to the tetragonal point group $C_{4v}$ takes the form

[TABLE]

where $g_{2}$ determines the relative weight between first and second order contributions, just like for the cubic point group $O$. Its magnitude and direction on the Fermi surface is illustrated in Fig. 6.

But unlike $O$, this point group has line nodes of the SOC in the BZ. For all values of $g_{2}$ the SOC is identically zero along the three paths parallel to the $z$-axis, $\Gamma\rightarrow\text{Z}$, $\text{X}\rightarrow\text{R}$, and $\text{M}\rightarrow\text{A}$ in Fig. 7 (a). Given the simple cubic first order tight-binding dispersion, for the range of $\mu$ we study the line node $\Gamma\rightarrow\text{Z}$ intersects all closed, and the line node $\text{X}\rightarrow\text{R}$ intersects all open Fermi surfaces. Thus $\min|{\bf l}({\bf k}^{-}_{F})|=0$ for both cases. The transition between the two is therefore seamless, and there are no fully gapped phases with $\text{sign}[\Delta_{-}({\bf k}^{-}_{F})]=-1$. The only two distinct phases is a topologically trivial, $\nu=0$, and a nodal non-trivial phase, $N_{\mathcal{L}}=1$, shown in white and grey respectively in Fig. 7 for a closed Fermi surface, $\mu=-50\alpha$.

Despite there only being a single topologically non-trivial phase the order parameter is calculated self-consistently for the two values $g_{2}\in\{0,4\}$ in order to study the effect of second order contributions to the SOC vector. This is done for nine values of the scaled bulk singlet to triplet ratio, $r^{\text{bulk}}_{\Delta}\in[0,1.1]$, with one active channel. The exact values are shown in table 1. The order parameter is calculated with the same surface normals as for the cubic point group. How the order parameter suppression depends on the surface orientation can be seen in Fig. 8 (a) - (b). Here the greatest suppression is not for the surface normal ${\bf n}=(1,1,1)$, but rather ${\bf n}=(1,1,0)$, and ${\bf n}=(0,1,2)$ shows very little suppression.

The zero-bias conductances for the two $g_{2}$ values, are very dissimilar for surface normals in the $xy$-plane. With $g_{2}=0$, Fig. 8 (c), rather large conductances are seen for $0.69\leq r^{\text{bulk}}_{\Delta}\leq 0.96$ in between the high symmetry axes ${\bf n}=(1,0,0)$ and ${\bf n}=(1,1,0)$, with the largest for $r^{\text{bulk}}_{\Delta}=0.83$ and ${\bf n}\approx(1,0.44,0)$. The lines corresponding to $0<r^{\text{bulk}}_{\Delta}\leq 0.55$ are (almost) degenerate due to all of them having smaller singlet to triplet ratios than the rather small difference between the maximum and minimum value of the SOC in the $xy$-plane of the Fermi surface. Only a few trajectories around the poles contribute to the ZBCPs. There are no ZBCPs for ${\bf n}=(1,1,0)$ but rather a large dome-like feature which is interestingly higher than the peaks for the (almost) degenerate lines. With $g_{2}=4$, Fig. 8 (d), the lines corresponding to $0.69\leq r^{\text{bulk}}_{\Delta}\leq 0.96$ show a dip for ${\bf n}\approx(1,0.4,0)$ due to the higher order contributions in the SOC changing the shape of the nodal rings, causing their projection onto the surface to largely overlap for these singlet to triplet ratios.

The ABS are heavily affected by self-consistency and are, just like for $O$, spin polarized. In Fig. 9 (a) the quantity $N^{(z)}({\bf k})$, Eq. (33), is plotted in the $xy$-plane for a pure triplet order parameter and ${\bf n}=(1,0,0)$. The largest effect of self-consistency is seen for glancing trajectories and energies between approximately $|\epsilon|/\max(\Delta)\in[0.5,1]$. The ABS in this range are not present in the non-self-consistent case. For a pure triplet order parameter $N^{(x)}=N^{(y)}=0$.

The momentum-resolved zero-energy ABS, for $g_{2}=0$, $r^{\text{bulk}}_{\Delta}=0.69$ and ${\bf n}=(1,1,1)$, is shown in Fig. 9 (b). Here, not only the trajectories for which $\Upsilon_{\bf k}=-\Upsilon_{\underline{\bf k}}\not=0$, colored green and yellow, give rise to ABS, but also trajectories for which $\Upsilon_{\bf k}=\Upsilon_{\underline{\bf k}}=-1$ and Eq. (34) holds, colored magenta. This magenta line is there due to the SOC vanishing along the high symmetry axis $\Gamma\rightarrow\text{Z}$, see Fig. 7 (a), combined with the SOC vector being perpendicular to this axis. This line is present for $g_{2}=4$ as well, but only for $0\leq r^{\text{bulk}}_{\Delta}<0.69$, whereas it is present for $0\leq r^{\text{bulk}}_{\Delta}<0.96$ with $g_{2}=0$, amongst the ratios investigated.

In the tunnel conductance spectra, plotted for ${\bf n}=(1,1,1)$ and $t_{0}=10^{-\frac{1}{2}}$ in Fig. 9 (c), ZBCPs are seen for all scaled singlet to triplet ratios in the interval $r^{\text{bulk}}_{\Delta}\in(0,1)$, due to $\min|{\bf l}({\bf k})|=0$. Unlike the ZBCPs in the tunnel conductance spectra for $O$ with $g_{2}=1.03$, Fig. 5 (c), which emanate from valleys around $\epsilon=0$, most of the ZBCPs here emanate from a large dome. The domes are a consequence of the magenta colored ABS together with the ’flatness’ of the coherence functions in the denominator of the expression for the tunnel conductance when varying the momentum, such that trajectories with momenta in close vicinity to the ABS condition give rise to a large number of states that contribute considerably to the tunnel conductance. With increasing $r^{\text{bulk}}_{\Delta}$ these states decrease in number. For $r^{\text{bulk}}_{\Delta}>0.83$ they have completely disappeared and thus the dome is gone and the ZBCP emanates from a valley.

The topological invariants $N_{(111)}$ and $W_{(111)}$ are plotted in Fig. 9 (d). Light blue/green corresponds to $N_{(111)}=\pm 1$ and white to trivial values of both invariants. The dashed circle is the projection of the spherical Fermi surface used in the quasiclassical calculations. The red line is given by $W_{(111)}=-1$. Thus states corresponding to solutions of Eq. (34) are directly related to the $\mathbb{Z}_{2}$ invariant being non-trivial, and are topologically protected as well.

In Fig. 9 (e) the band structure is shown for $g_{2}=0$ and $r^{\text{bulk}}_{\Delta}=0.69$ along the $k^{\parallel}_{2}$-axis with $k^{\parallel}_{1}=0$, and $L=1.3\cdot 10^{4}$ layers. Note that the states corresponding to $W_{(111)}=-1$ are doubly degenerate on each surface. Just like for $O$ the zero-energy bands are given by the projection of the non-trivial values of the topological invariants, which can be seen in Fig. 9 (f).

III.3 The Tetrahedral Point Group $T_{d}$

To next-nearest neighbors in the sum over Bravais lattice sites Samokhin (2009) the SOC vector corresponding to the tetrahedral point group $T_{d}$ takes the form

[TABLE]

with no free parameter $g_{2}$ in contrast with $O$ and $C_{4v}$. It is illustrated in Fig. 10. This SOC exhibits line nodes in the BZ along the paths $\Gamma\rightarrow\text{P}\rightarrow\text{H}\rightarrow\Gamma$ and $\text{P}\rightarrow\text{N}$ in Fig. 11 (a). Just like for $C_{4v}$ the line nodes intersect the negative helical Fermi surface for all values of $\mu$, i.e. all Fermi surface geometries, given a BCC first order tight-binding dispersion. Thus $\min|{\bf l}({\bf k}^{-}_{F})|=0$ and there are no gapped phases with $\text{sign}[\Delta_{-}({\bf k}^{-}_{F})]=-1$, which can be seen in the topological phase diagram in Fig. 11 (b). Just like for $C_{4v}$ there are only two distinct topological phases; one gapped trivial, $(N_{\mathcal{L}},\nu)=(0,0)$, and a nodal non-trivial, $N_{\mathcal{L}}=1$, phase.

Due to there being no free parameter to vary the self-consistent order parameter is only calculated with this single SOC vector for this point group. This is done for nine values of the scaled bulk singlet to triplet ratio, $r^{\text{bulk}}_{\Delta}\in[0,1.1]$, with one active channel. The exact values are shown in table 1. The suppression of these order parameters is shown in Fig. 12 (a) for the same range of surface normals as for the previous examined point groups. Here the largest suppression is for ${\bf n}=(1,1,0)$ and ${\bf n}=(0,1,2)$ and barely any suppression at all for ${\bf n}=(1,1,1)$.

In Fig. 12 (b) the zero-bias conductance for these order parameters and surface normals are shown. Unsurprisingly the ZBC is very small for the high-symmetry axes ${\bf n}=(1,0,0)$ and ${\bf n}=(1,1,0)$, it is quite large, but still a local minima, for the high-symmetry axis ${\bf n}=(1,1,1)$, and larger still in between these surface normals.

The ABS are spin polarized for this point group as well. In Fig. 13 (a) the quantity $N^{(z)}({\bf k})$, Eq. (33), is plotted for a pure triplet order parameter and ${\bf n}=(1,0,0)$. For a pure triplet $N^{(x)}=N^{(y)}=0$. Self-consistency does not drastically alter the ABS for this surface normal due to the states being predominantly located at small energies for glancing trajectories.

The momentum-resolved zero-energy ABS for $r^{\text{bulk}}_{\Delta}=0.69$ are shown in Fig. 13 (b). The states in middle are from the non-overlapping parts of the projection of the nodal rings around the $\Gamma\rightarrow\text{P}$ high symmetry axis, and the ones around the edges of the disk from the projection of the nodal rings around $\Gamma\rightarrow\text{H}$ and $\Gamma\rightarrow\text{N}$.

Just like for the other point groups ZBCPs are seen in the tunnel conductance spectra, with $t_{0}=10^{-\frac{1}{2}}$, for ${\bf n}=(1,1,1)$ and singlet to triplet ratios in the interval $\min|{\bf l}_{{\bf k}_{F}}|<\Delta_{s}/\Delta_{t}<\max|{\bf l}_{{\bf k}_{F}}|$, i.e. $0<r^{\text{bulk}}_{\Delta}<1$, see Fig. 13 (c). ABS given by Eq. (34) only appear for $0\leq r^{\text{bulk}}_{\Delta}\leq 0.28$, and then not for ${\bf k}_{\parallel}=0$ which is the most important momentum when calculating the tunnel conductance 111The angle-resolved tunnel conductance is proportional to the DOS weighted by $\cos(\alpha_{\bf k})$ as well as the transmission amplitude squared.. Hence the ZBCPs emanating from valleys in the spectra.

The non-trivial values of the topological invariant $N_{(111)}$ ($W_{(111)}$ being trivial for this singlet to triplet ratio) is shown in Fig. 13 (d). The dashed circle is projection of the spherical Fermi surface used in the quasiclassical calculations. Compared to the other point groups considered the spherical Fermi surface approximation does not work as well due to the actual Fermi surface bulging out in the ${\bf k}=(1,1,1)$ direction. Furthermore, the BZ is not cubic and thus the line integral defining $N_{(lmn)}$ potentially goes through Fermi surfaces from adjacent BZs, which is precisely what happens for this surface normal. The slightly complicated structure near the circle thus stems from the partial overlap of the nodal rings of the adjacent Fermi surface combined with $N_{(lmn)}$ being additive.

Zero-energy flat bands are seen in the band structure for this point group and surface normal as well, shown for $r^{\text{bulk}}_{\Delta}=0.69$ along the $k^{\parallel}_{2}$-axis, with $k^{\parallel}_{1}=0$ and $L=1.3\cdot 10^{4}$ layers, in Fig. 13 (e).

As seen in the plot of the lowest positive band, Fig. 13 (f), the zero-energy states around the origin are somewhat patchy at this resolution, hence the slight gap at $k^{\parallel}_{2}\approx\pm 0.28|{\bf k}_{F}|$ in Fig. 13 (e). Furthermore, there is a small gap at $k^{\parallel}_{2}\approx\pm 0.8|{\bf k}_{F}|$ in Fig. 13 (e), but this is not an artifact of the lower resolution of the band structure compared to the topological invariant plot, Fig. 13 (d), as there is a small region separating $N_{(111)}=\pm 1$ at these momenta.

IV Conclusions

We have theoretically studied noncentrosymmetric superconductors self-consistently for the point groups $O$, $C_{4v}$, and $T_{d}$, with a closed Fermi surface. Four values of $g_{2}$, parameterizing the relative weight of first and second order contributions in the spin-orbit coupling (SOC), given by the Bravais lattice sum up to next-nearest neighbors, were chosen for $O$ in order to investigate all its gapped topological phases for a closed Fermi surface. The point groups $C_{4v}$ and $T_{d}$ were shown to have no gapped topological phases, yet two values of $g_{2}$ were chosen for $C_{4v}$ in order to study the effect of second order contributions in the SOC vector. For $T_{d}$ no higher order terms in the SOC were seen up next-nearest neighbors.

The reason for the existence of gapped topological phases for $O$ was shown to be due to the fact that the SOC only vanishes in the Brillouin zone at high symmetry points, whereas the SOC vanishes at certain high symmetry axes for $C_{4v}$ and $T_{d}$. It was shown for $O$ that the topology changes at the Lifshitz transition, i.e. at the transition point between an open and closed Fermi surface. This does not happen for $C_{4v}$ and $T_{d}$ and the Lifshitz transition is topologically seamless.

In the bulk it was shown that there are two distinct mixed states; with one or two nucleation channels. In both cases it was demonstrated that there is a possibility of a cross-over from dominating singlet to dominating triplet, or vice versa, with decreased temperature. Depending on the material this could be important if experiments are done at different temperatures. With two nucleation channels there is a possibility of a second phase transition at the subdominant critical temperature and it was shown by explicit construction that the subdominant channel for certain parameter values indeed has lower free energy. If this can be extended to more complicated Fermi surface geometries and parameter values remains to be seen.

The order parameter suppression’s dependence on surface orientation and singlet to triplet order parameter ration was studied for a range of different surface normals. The suppression was seen to be highly dependent on surface orientation.

The Andreev bound states (ABS) are found to be spin polarized with different polarization axes for different singlet to triplet ratios. The order parameter suppression affects the ABS heavily for glancing trajectories and sub-gap energies close to the gap, and less for smaller energies. Zero-energy states are not affected by the calculated suppression. Thus the zero-bias conductance peaks are present in the non-self-consistent tunnel conductance as well.

We showed that the zero-energy surface states are topological in nature. Thus it is clear that the calculated suppression should not affect the zero-energy states due to the gap not vanishing at any distance from the surface. If this can happen for other parameters and/or surface orientations is an open question.

Acknowledgements.

We appreciate the highly stimulating atmosphere within the Hubbard Theory Consortium and during the annual “Condensed Matter Physics in the City” events in London. N.W. would like to thank Roland Grein for discussions in the early stages of the work, and Patric Holmvall for discussions in the later stages. M.E. acknowledges support by EPSRC (Grant No. EP/J010618/1 and EP/N017242/1). N.W. acknowledges financial support by the Southeast Physics Network (SEPnet) for his Ph.D. study at Royal Holloway, University of London.

Appendix A Temperature Dependence of the Gap

The self-consistency equation for the order parameter Eq. (27), can be written symbolically in the form of a fixed point equation

[TABLE]

where $\bm{\Delta}=(\Delta_{s},\Delta_{t})^{T}$, and the function ${\bf F}(\bm{\Delta})$ is simply a short-hand notation for the right hand side of Eq. (27). Any $\bm{\Delta}$ that obeys eq. (56) is called a fixed point. Then a iteration scheme is employed to find a convergence to a fixed point. This yields a series of points $\bm{\Delta}_{1}$, $\bm{\Delta}_{2}$, $\ldots$, which hopefully converges to a solution. The procedure is said to have converged when the difference between iterations is sufficiently small

[TABLE]

where the number $c$ is the convergence criterion. In the bulk the fixed points can be obtained by computing ${\bf F}(\bm{\Delta})$ for a vast number of points.

We illustrate the method for the case of two attractive channels. Because the number of possible independent attractive fixed points is equal to the number of positive eigenvalues to the matrix $L$, one has for values of $(v_{s},v_{t},v_{m})$ in the yellow oval in Fig. 1 two nucleation channels. However, the subdominant channel does in general not nucleate at $T_{\text{c}}$, but at a lower temperature, $T^{\text{sub.}}_{\text{c}}<T_{\text{c}}$. Thus, if one follows the procedure in the previous paragraph for the initial guesses $\bm{\Delta}_{0}$ one will not see the possible transition to the subdominant channel. What is needed in this case is to calculate the order parameter with increasing temperature instead of decreasing. By computing a few iterations, $n\sim 20$, at a sufficiently low temperature, say $T=0.1T_{\text{c}}$ (which must be smaller than $T^{\text{sub.}}_{\text{c}}$ obviously), for a number of random initial guesses, an attractive fixed point corresponding to the subdominant channel is obtained, and is denoted $\bm{\Delta}^{\text{sub.}}$. For the lowest temperature the initial guess will thus be $\bm{\Delta}^{\text{sub.}}$, and subsequent guesses $\bm{\Delta}_{0}(T+\delta T)=\bm{\Delta}_{n}(T)$. By calculating the order parameter this way, it will converge to the subdominant channel value until $T\leq T^{\text{sub.}}_{\text{c}}$. At $T=T^{\text{sub.}}_{\text{c}}$ the order parameter transitions to the dominant channel value due to it being the only attractive fixed point at these temperatures, unless the subdominant channel value at $T=T^{\text{sub.}}_{\text{c}}$ is zero, $\bm{\Delta}_{n}(T^{\text{sub.}}_{\text{c}})={\bf 0}$, in which case it will stay zero. In this manner, the subdominant critical temperatures are obtained.

A choice of parameter values yielding an admixture of singlet and triplet with an attractive subdominant channel is e.g. $(v_{s},v_{t},v_{m})=(1,1/\langle|{\bf l}_{\bf k}|^{2}\rangle,-0.1)$ (ignoring normalization). In Fig. 14 examples for the fixed point iteration are shown for the point group $O$ (the plots look qualitatively similar for all other point groups). The grey blobs correspond to the function $f(\bm{\Delta})=|\bm{\Delta}-{\bf F}(\bm{\Delta})|$. Darker indicates smaller values of $f(\bm{\Delta})$, and pure black indicates the existence of a fixed point. The colored circles connected by lines show the convergence of 25 random initial guesses, $\bm{\Delta}_{0}$, progressing a number of iteration steps. The colors of the circles indicate the iteration number $n$. Starting with dark blue for $n=0$, transitioning through cyan, green, yellow, and ending with red for $n=n_{\text{max}}$. Any fixed point $\bm{\Delta}_{n_{\text{max}}}$ converges to is an attractive fixed point, however there are repulsive fixed points present in the subdominant channel. Concentrating on the subdominant channel, Fig. 14, one notices in addition to the attractive fixed points also two repulsive fixed points. From numerical investigations this seems to be a general feature, and the fixed points roughly fall on a parallelogram with the attractive fixed points at the vertices.

A criterion if a second phase transition exists can be obtained from the condensation energy, Eq. (31). Thus there is a second phase transition if it holds that

[TABLE]

at zero temperature. For certain parameters this is indeed the case. In general, how small $T^{\text{sub.}}_{\text{c}}$ can be without losing the second phase transition depends on the point group. The key to get a second phase transition is to choose $(v_{s},v_{t},v_{m})$ in such a way as to get a dominant channel with a large triplet component, as well as a rather large subdominant critical temperature.

Appendix B Zero-bias Andreev bound states

In the bulk, for zero energy and for real order parameter, the coherence functions take a particularly simple form,

[TABLE]

The values for the parameters $(v_{s},v_{t},v_{m})$ were chosen to yield positive singlet and triplet components, and we can therefore simplify the expressions further to

[TABLE]

where $\Upsilon_{\bf k}\equiv\text{sgn}\left[\Delta_{s}/\Delta_{t}-|{\bf l}_{\bf k}|\right]$. We are interested in zero-bias states protected by topology, for which it suffices to discuss the non-self-consistent order parameter, i.e. bulk values all the way to the surface. Because $\Upsilon_{\bf k}$ and $\Upsilon_{\underline{\bf k}}$ can take three values each there are $3^{2}=9$ different cases to consider. They are listed below together with the equations for surface ABS, using the (real) reflection amplitude $0\leq r\leq 1$. The solutions separate naturally into three groups. The first group

[TABLE]

has no solutions. Therefore, there can be no ZBCPs for $r^{\text{bulk}}_{\Delta}>\max|{\bf l}_{{\bf k}_{\text{F}}}|$, because $\Upsilon_{\bf k}=+1\;\forall{\bf k}$. The second group,

[TABLE]

requires $r=1$, as well as a sign change of $\Upsilon$ when reflected at the surface. The third group is

[TABLE]

where $F\left(\phi_{l},\theta_{l},\phi_{\underline{l}},\theta_{\underline{l}}\right)\in[-1,1]$ is defined after Eq. (34). Eq. (68) has no solution for real $r$, and the remaining equations have solutions only for $r=1$. Thus, there are only two classes of trajectories giving rise to ABS at zero energy, the first class given by Eqs. (66)-(67) and the second class given by Eqs. (69)-(B).

Fig. 15 shows solutions to Eqs. (69) - (B) and Eqs. (69) - (B) for two surface normal directions.

Appendix C 3D winding number

The 3D winding number is given by Schnyder et al. (2008)

[TABLE]

where $M_{a}=q^{-1}\partial_{a}q$. Introducing the notation $q({\bf k})=C_{0}({\bf k})\sigma_{0}+{\bf C}({\bf k})\cdot\bm{\sigma}$ with ${\bf C}=(C_{1},C_{2},C_{3})$, and $q^{-1}=\left[C^{2}_{0}-|{\bf C}|^{2}\right]^{-1}\left(C_{0}\sigma_{0}-{\bf C}\cdot\bm{\sigma}\right)$, and $R=\left[C^{2}_{0}-|{\bf C}|^{2}\right]$, and the $4\times 4$ matrix

[TABLE]

we find the exact formula for the trace

[TABLE]