Transverse profile and 3D spin canting of a Majorana state in carbon nanotubes

TL;DR

This paper provides a detailed 3D spatial and spin profile of Majorana bound states in carbon nanotubes, revealing their helical structure and symmetry properties through analytical and numerical methods.

Contribution

It introduces an analytical reconstruction of the Majorana wave function in CNTs considering spin and valley degrees of freedom, highlighting the lattice's chiral imprint.

Findings

Majorana wave function has a helical, anisotropic transverse structure.

Spin canting angle exhibits a spiral pattern around the CNT circumference.

Three contributions to the MBS are identified from Fourier analysis.

Abstract

The full spatial 3D profile of Majorana bound states (MBS) in a nanowire-like setup featuring a semiconducting carbon nanotube (CNT) as the central element is discussed. By accurate tight-binding calculations we show that the chiral nature of the CNT lattice is imprinted in the MBS wave function which has a helical structure, anisotropic in the transverse direction. The local spin canting angle displays a similar spiral pattern, varying around the CNT circumference. We reconstruct the intricate 3D profile of the MBS wave function analytically, using an effective low energy Hamiltonian accounting both for the electronic spin and valley degrees of freedom of the CNT. We find that the four components of the Majorana spinor are related by the three symmetries of our Bogoliubov-de Gennes (BdG) Hamiltonian, reducing the number of independent components to one. A Fourier transform analysis uncovers the presence of three contributions to the MBS, one from the $\Gamma$-point and one from each of the Fermi points, with further complexity added by the presence of two valley states in each contribution.

Tables

Table 1. Table 1: Values of k F subscript 𝑘 𝐹 k_{F} , κ Γ subscript 𝜅 Γ \kappa_{\Gamma} and κ F subscript 𝜅 𝐹 \kappa_{F} from the analytical calculation compared with values fitted from the numerics.

	Analytics ($\frac{1}{\mu m}$)	Fits ($\frac{1}{\mu m}$)
${\kappa }_{\Gamma }$	-7.94	-8.93
${\kappa }_F$	-6.56	-8.01
$k_F$	118.92	115.25

Equations

$$⟨{\vec{u}}_M\left(\vec{r}\right)⟩s_x{\vec{u}}_M\left(\vec{r}\right)$$

$$⟨{\vec{u}}_M\left(\vec{r}\right)⟩s_y{\vec{u}}_M\left(\vec{r}\right)$$

$$⟨{\vec{u}}_M\left(\vec{r}\right)⟩s_z{\vec{u}}_M\left(\vec{r}\right)$$

$${\theta }_{xy}(\vec{r})=\arctan \left(\frac{⟨s_y(\vec{r})⟩}{⟨s_x(\vec{r})⟩}\right)=\pi /2-2\arg (u_{M\uparrow }(\vec{r})).$$

$$H\left(k\right)=\left(\begin{array}{cccc}{\xi }_{K\uparrow }\left(k\right)& {\mu }_B{B}_{\perp }& {\Delta }_{K{K}^{\prime }}& 0\\ {\mu }_B{B}_{\perp }& {\xi }_{K\downarrow }\left(k\right)& 0& {\Delta }_{K{K}^{\prime }}\\ {\Delta }_{K{K}^{\prime }}& 0& {\xi }_{K^{\prime }\uparrow }\left(k\right)& {\mu }_B{B}_{\perp }\\ 0& {\Delta }_{K{K}^{\prime }}& {\mu }_B{B}_{\perp }& {\xi }_{K^{\prime }\downarrow }\left(k\right)\end{array}\right)\text{,}$$

$${\tilde{\mathcal{H}}}_{\text{BdG}}=\left(\begin{array}{cccc}{\tilde{E}}_1\left(k\right)& 0& {\tilde{\Delta }}_p\left(k\right)& -{\tilde{\Delta }}_s\left(k\right)\\ 0& {\tilde{E}}_2\left(k\right)& {\tilde{\Delta }}_s\left(k\right)& {\tilde{\Delta }}_p\left(k\right)\\ {\tilde{\Delta }}_p\left(k\right)& {\tilde{\Delta }}_s\left(k\right)& -{\tilde{E}}_1\left(k\right)& 0\\ -{\tilde{\Delta }}_s\left(k\right)& {\tilde{\Delta }}_p\left(k\right)& 0& -{\tilde{E}}_2\left(k\right)\end{array}\right)\text{.}$$

$${\hat{\mathcal{H}}}_{\text{BdG}}^{\pm }=\left(\begin{array}{cc}{\tilde{\xi }}_{\pm }\left(k\right)& {\tilde{\Delta }}_p\left(k\right)\\ {\tilde{\Delta }}_p\left(k\right)& -{\tilde{\xi }}_{\pm }\left(k\right)\end{array}\right)\text{.}$$

$${\tilde{\xi }}_{\pm }\left(k\right)=$$

$$\pm \frac{1}{2}\sqrt{{\left({\tilde{E}}_1\left(k\right)+{\tilde{E}}_2\left(k\right)\right)}^2+4{\tilde{\Delta }}_s^2\left(k\right)}.$$

$$|{\Psi }_M⟩=\frac{A_{\Gamma }}{\sqrt{2}}|{\Psi }_{\Gamma }⟩+\frac{A_R}{\sqrt{2}}|{\Psi }_{k_F}⟩+\frac{A_L}{\sqrt{2}}|{\Psi }_{-k_F}⟩\text{.}$$

$$\mathcal{P}{\Psi }_M\left(z\right)\stackrel{!}{=}{\Psi }_M\left(z\right)\text{,}$$

$${\Psi }_M\left(z=0\right)\stackrel{!}{=}0\text{,}$$

$${\int }_0^{\infty }dz{\left|{\Psi }_M\left(z\right)\right|}^2\stackrel{!}{=}1\text{.}$$

$${\Psi }_{M,1D}\left(z\right)=\left[\frac{A_{\Gamma }}{\sqrt{2}}{e}^{{\kappa }_{\Gamma }z}+\frac{A_R}{\sqrt{2}}{e}^{{\kappa }_{R}z}+\frac{A_L}{\sqrt{2}}{e}^{{\kappa }_R^{\star }z}\right]\left(\begin{array}{c}\mp i\\ 1\end{array}\right)\text{.}$$

$${\Psi }_{M,1D}\left(z\right)=\frac{\mathcal{N}}{2}\left({\psi }_{\parallel }\left(z\right)+{\psi }_{\parallel }^{\star }\left(z\right)\right)\frac{1}{\sqrt{2}}\left(\begin{array}{c}1-i\\ 1+i\end{array}\right)\text{,}$$

$${\psi }_{\parallel }\left(z\right)=\left(e^{{\kappa }_{F}z+i{k}_{F}z}-e^{{\kappa }_{\Gamma }z}\right)$$

$${\hat{\gamma }}_M$$

$${\hat{\gamma }}_M=\sum _{k,\tau ,s}\left(u_{\tau s}\left(k\right)c_{k\tau s}+v_{\tau s}\left(-k\right)c_{-k\tau s}^{†}\right)\text{,}$$

$$u_{\tau s}\left(k\right)=\tau s{\Lambda }_{\tau s}\left(k\right){\lambda }_s\left(k\right)\text{,}$$

$$v_{\tau s}\left(k\right)=\tau s{\Lambda }_{\tau s}\left(k\right){\lambda }_s^{\star }\left(k\right)\text{,}$$

$${\Lambda }_{\tau s}\left(k\right)=\{\begin{array}{ll}{a}_s\left(k\right)& \text{for }\tau =+1\text{,}\\ b_s\left(k\right)& \text{for }\tau =-1\text{,}\end{array}$$

$${\lambda }_s\left(k\right)=um\left(k\right)g\left(sk\right)-svn\left(k\right)h\left(sk\right)\text{.}$$

$$u_{p\uparrow }\left(\mathbf{r}\right)$$

$$+\frac{A_R}{\sqrt{2}}\left[e^{ip{\eta }_{K\uparrow k_F}}{u}_{K\uparrow }\left(k_F\right)e^{i{k}_{\perp }{r}_{\perp }}+e^{ip{\eta }_{K^{\prime }\uparrow k_F}}{u}_{K^{\prime }\uparrow }\left(k_F\right)e^{-i{k}_{\perp }{r}_{\perp }}\right]e^{{\kappa }_{F}z+i{k}_{F}z}$$

$$+\frac{A_L}{\sqrt{2}}\left[e^{ip{\eta }_{K\uparrow -k_F}}{u}_{K\uparrow }\left(-k_F\right)e^{i{k}_{\perp }{r}_{\perp }}+e^{ip{\eta }_{K^{\prime }\uparrow -k_F}}{u}_{K^{\prime }\uparrow }\left(-k_F\right)e^{-i{k}_{\perp }{r}_{\perp }}\right]e^{{\kappa }_{F}z-i{k}_{F}z}\text{.}$$

$$u_{p\uparrow }\left(\vec{r}\right)=$$

$$+$$

$$H_0=\sum _{\mathbf{k},s}{\gamma }_s\left(\mathbf{k}\right)a_{\mathbf{k}s}^{†}{b}_{\mathbf{k}s}+\text{h.c.},$$

$${\gamma }_s\left(\mathbf{k}\right)=t_{s,1}{e}^{i\mathbf{k}\cdot {\mathbf{a}}_1}+t_{s,2}{e}^{i\mathbf{k}\cdot {\mathbf{a}}_2}+t_{s,3},$$

$$H_{\text{CNT}}-\mu N=\sum _{k,\tau ,s}{\xi }_{\tau s}\left(k\right)c_{k\tau s}^{†}{c}_{k\tau s}\text{,}$$

$${\epsilon }_{\tau s}\left(k_{\perp },k\right)$$

$$H_{K{K}^{\prime }}=\sum _{k,s}{\Delta }_{K{K}^{\prime }}\left(c_{kKs}^{†}{c}_{k{K}^{\prime }s}+c_{k{K}^{\prime }s}^{†}{c}_{kKs}\right)\text{,}$$

$$H_Z=\sum _{k,\tau }{\mu }_B{B}_{\perp }\left(c_{k\tau \uparrow }^{†}{c}_{k\tau \downarrow }+c_{k\tau \downarrow }^{†}{c}_{k\tau \uparrow }\right)\text{.}$$

$$\left(\begin{array}{c}{c}_{kKs}\\ c_{k{K}^{\prime }s}\end{array}\right)=\left(\begin{array}{cc}{a}_s\left(k\right)& b_s\left(k\right)\\ -b_s\left(k\right)& a_s\left(k\right)\end{array}\right)\left(\begin{array}{c}{\alpha }_{ks}\\ {\beta }_{ks}\end{array}\right)\text{,}$$

$$a_s^2\left(k\right)$$

Full text

Over the past decade Majorana fermions have been of great interest in condensed matter physics. Under special conditions they arise as quasiparticles in superconductors,Aguado (2017) where they are zero energy eigenstates of the Bogoliubov-de Gennes (BdG) Hamiltonian and of the particle-hole symmetry operator. Theoretically such quasiparticles were predicted to appear in the elusive one-dimensional $p$-wave superconductors Kitaev (2001); but it is also possible to engineer $s$-wave systems in such a way that they mimic $p$-wave superconductivity Sato and Ando (2017). The most popular setup is based on semiconducting nanowires with large spin-orbit interaction and large $g$-factor in contact with a superconductor, which induces superconducting proximity correlations in the wire Lutchyn et al. (2010); Oreg et al. (2010). Although the experiments are by now very advanced Lutchyn et al. (2018), a definite proof that the reported signatures Mourik et al. (2012); Churchill et al. (2013); Deng et al. (2016); Zhang et al. (2017) are really due to the topologically non trivial Majorana bound states (MBS) is still missing. Thus, recent proposals have suggested to use local probes to infer exclusive properties of a MBS, such as its nonlocality and its peculiar spin canting structureLiu et al. (2017); Prada et al. (2017); Clarke (2017); Spanton et al. (2017); Hoffman et al. (2017); Schuray et al. (2018), or the maximal electron-hole content of the Majorana spinor Sticlet et al. (2012); Sedlmayr and Bena (2015). However, in order to exclude spurious effects, local experiments can be truly useful only if the spatial profile of the MBS is known with sufficient accuracy. This is very difficult to achieve for the case of the semiconducting nanowires, since their diameter of a few tens of nanometers and their length of several hundreds of nanometers do not allow for a microscopic calculation of the MBS wavefunction. Typically, the spatial profile is obtained with simple one-dimensional models Klinovaja and Loss (2012). The transverse profile has so far been obtained numerically for effective models: of core-shell nanowires in cylindrical Lim et al. (2013); Osca et al. (2014) and prismatic Manolescu et al. (2017); Stanescu et al. (2018), and of full nanowires in hexagonal Woods et al. (2018) geometries.

In this work we show that the spatial profile of MBS can be derived analytically with good accuracy in a setup which uses a carbon nanotube (CNT) in proximity with an $s$-wave superconductor. Similar to the nanowires, such CNTs can host MBS at their ends Egger and Flensberg (2012); Klinovaja et al. (2012); Sau and Tewari (2013); Hsu et al. (2015); Marganska et al. (2018). Due to their hollow character and small diameter, CNTs of several micrometers can be simulated numerically based on tight-binding models of carbon atoms on a rolled graphene lattice Izumida et al. (2009); Klinovaja et al. (2011). Such simulations allow one to accurately evaluate the excitation spectrum and local observables. Effective single-particle low energy models can be derived which well reproduce microscopic simulations Ando (2000).

In a recent paper Marganska et al. (2018) we have used a four-band and an effective one-band model to calculate the topological phase diagram and the energy spectrum of proximitized semiconducting CNTs in perpendicular magnetic field, see Fig. 1(a), with parameters obtained from a fit to the numerical spectra111In our tight binding model we consider one $p_{z}$ orbital per atom..

In this work we use the same models to analytically obtain the full 3D spatial profile of the Majorana wave function. First, we exploit our knowledge of the three symmetries of the effective BdG Hamiltonian in order to derive the relations between the four components of the Majorana spinor (see Fig. 1(e,f)), thus reducing the number of independent components to one. Second, we find that the presence of two angular momentum contributions (valleys) and the spin degree of freedom results in the formation of a composite, six-piece MBS whose 3D wave function has a distinctive spiral pattern with a $C_{2}$ symmetry, impossible to factorize into separate transverse and longitudinal profiles. Equally non-isotropic is the spin canting angle, a quantity encoding the relative phase of the spin up and spin down particle components of the Majorana wave function. A comparison with the numerical results for the MBS of a (12,4) CNT gives us confidence in the reliability of the effective model. Our results show that while simple 1D models can capture the important low energy properties of the BdG spectrum, they might miss crucial features present in the full 3D wave function. This can have profound implications in various setups, where the shape and local spin composition of an MBS are relevantPrada et al. (2017); Hoffman et al. (2017); Schuray et al. (2018).

The paper has the following structure. In Sec. I we discuss our microscopic model of the carbon nanotube, the symmetries of the BdG Hamiltonian in our setup and the resulting relations between the components of the Majorana spinor. In Sec. II we show and discuss the numerical results of the spin canting of the full 3D MBS. We proceed to reconstruct the MBS analytically. First we introduce in Sec. III the effective low energy model of the carbon nanotube, including the superconducting correlations. We also derive the form of the Majorana state in a continuum 1D approximation. In Sec. IV we calculate the 3D Majorana solution and determine its full spatial profile. Finally we compare the numerical results from the real-space tight-binding calculation with those of the analytical model.

I Model and its symmetries

Geometrically, a single wall carbon nanotube is equivalent to a rolled-up strip taken from the two-dimensional honeycomb of carbon atoms that makes up a graphene sheet Saito et al. (1998). The band structure of the CNT can be obtained from that of graphene by imposing periodic boundary conditions in the transverse direction, which quantize the transverse momentum, turning the two-dimensional dispersion of graphene into a series of 1D cuts, which are the CNTs one-dimensional subbands, shown schematically in Fig. 1(b). Effective low-energy Hamiltonians can be derived from the microscopic modelAndo (2000). Thus, like in graphene, the low-energy band structure in nanotubes consists of two distinct and time-conjugate valleys $K$ and $K^{\prime}$ which are indexed by the quantum number $\tau$ ($\tau=+1$ for $K$ valley and $\tau=-1$ for $K^{\prime}$ valley) (cf. Fig. 1(b)). However, the simple fact of being rolled up drastically modifies the band structure, leading to effects that are not present in graphene. These are a curvature induced band gap and an enhanced spin-orbit coupling Ando (2000); Kuemmeth et al. (2008); Izumida et al. (2009); Klinovaja et al. (2011). The spin-orbit coupling in the nanotubes results in an effective spin-orbit field directed along the tube axis, with the sign of the field given by $\tau s$, with $s$ the spin quantum number along the CNT. The CNT’s tiny diameter reduces the number of relevant transverse modes to exactly four in the low-energy regime, one for each spin and valley. In order to keep the low energy physics close to the $\Gamma$ point, we consider nanotubes of the zigzag class Marganska et al. (2015); Izumida et al. (2016), where the Dirac points are only slightly shifted from $k=0$. In order to open the gap at the $\Gamma$ point, we need to remove the Kramers degeneracy between the ($\tau,s$) and ($-\tau,-s$) states. The spin degeneracy can be removed by a transverse magnetic field, but only if the valleys are also mixed. Fortuitously, this happens automatically when the nanotube is in contact with the bulk superconductor, i.e. the source of the proximity effect. Its presence breaks the rotational symmetry of the tube, introducing mixing between the $K$ and $K^{\prime}$ valley. The resulting spectrum in a normal CNT is shown in Fig. 1(c).

The proximity to a superconducting substrate induces Cooper pairing in the CNT. The excitation spectrum of the system can be determined from the BdG Hamiltonian, where the superconducting correlations are treated in a mean-field approximation. In the microscopic model this corresponds to an on-site pairing term Uchoa and Castro Neto (2007), see Fig. 1(a), and using the Nambu spinor we can construct the microscopic BdG Hamiltonian of our system. To anticipate the discussion in Sec. III, in the reciprocal space this pairing yields both an inter-band ($\tilde{\Delta}_{s}$, with $s$-wave symmetry) and an in-band ($\tilde{\Delta}_{p}$) pairing, with $p$-wave symmetry, required for topological superconductivity. The two pairings are shown in Fig. 1(d).

The CNT alone has a crystalline symmetry of rotation by $\pi$ around an axis perpendicular to the CNT ($C_{2}^{\prime}$ axis in Fig. 1(a)). In consequence, the CNT on superconducting substrate is a topological crystalline superconductor Shiozaki and Sato (2014); Ando and Fu (2015) with $C_{2}^{\prime}$ axis oriented as shown in Fig. 1(a). In our setup, however, the $C_{2}^{\prime}$ symmetry is broken by the magnetic field parallel to the substrate and only the local symmetries remain.

The true time reversal symmetry is broken by the magnetic field. Nevertheless, the inspection of the single-particle Hamiltonian of our CNT setup in the real space Ando (2000); del Valle et al. (2011); Marganska et al. (2018) shows that all its dominant terms possess a local antiunitary symmetry, which commutes with the Hamiltonian. Its action on the basis states is defined by $\tilde{\mathcal{T}}c\,|is\rangle=-ic^{*}\,|i,-s\rangle$. Contrary to the true time reversal, $\tilde{\mathcal{T}}$ has bosonic nature $\tilde{\mathcal{T}}^{2}=1$. The $\tilde{\mathcal{T}}$ is discussed further in the Appendix A.2.

The second local symmetry is the particle-hole symmetry $\mathcal{P}$, inherent in all BdG systems. With the $\mathcal{P}$ and $\tilde{\mathcal{T}}$ symmetries combined, the BdG Hamiltonian of the nanotube is also chiral symmetric under $\mathcal{C}=\tilde{\mathcal{T}}\mathcal{P}$. When acting on the eigenstates of the finite system, expressed in the Nambu space as $\hat{\Psi}(\vec{r})=\sum_{s}[u_{s}(\vec{r})c_{s}(\vec{r})+v_{s}(\vec{r})c_{s}^{\dagger}(\vec{r})]$, these operators convert between the $u_{s}$ and $v_{s}$ components of the different states in the way shown schematically in Fig. 1(e). (The $\tilde{\mathcal{T}}$ relation has been noticed in Ref. Hoffman et al., 2017, although without attributing it to the presence of a pseudo-time-reversal symmetry.) The complementary relations holding in the reciprocal space, calculated in Sec. IV, are shown in Fig. 1(f). The presence of these three symmetries has a profound impact on the Majorana state.

The wave function of the Majorana bound state is given by $\Braket{\vec{r}}{\Psi_{M}}=\Psi_{M}\left(\vec{r}\right)$, where $\Ket{\Psi_{M}}=\hat{\gamma}_{M}\Ket{0}$ and $\hat{\gamma}_{M}^{\dagger}=\hat{\gamma}_{M}$ is the Majorana creation operator. Here $\vec{r}=\left(z,r_{\perp}\right)$, where $z$ and $r_{\perp}$ denote the longitudinal and the transverse components, respectively. The MBS is described by a spinor, $\Psi_{M}\left(\vec{r}\right)=\left(u_{M\uparrow}\left(\vec{r}\right)\text{, }u_{M\downarrow}\left(\vec{r}\right)\text{, }v_{M\uparrow}\left(\vec{r}\right)\text{, }v_{M\downarrow}\left(\vec{r}\right)\right)^{T}$, with $u_{Ms}\left(\vec{r}\right)$ and $v_{Ms}\left(\vec{r}\right)$ the electron and hole components, respectively, and $s$ indicating the spin degree of freedom. As detailed below, it is enough to find the $u_{M\uparrow}\left(\vec{r}\right)$ components and use the symmetries of the underlying Bogoliubov-de Gennes (BdG) Hamiltonian to determine the rest.

The first relation is a consequence of the fundamental property $\mathcal{P}\Psi_{M}\left(\vec{r}\right)\overset{!}{=}\Psi_{M}\left(\vec{r}\right)$ of a Majorana state. Thus the relation $\mathcal{P}u_{s}\left(\vec{r}\right)=v_{s}^{\star}\left(\vec{r}\right)$ becomes $u_{Ms}\left(\vec{r}\right)=v_{Ms}^{\star}\left(\vec{r}\right)$. As we will show in Section III, the MBS are also eigenstates of the chiral symmetry $\mathcal{C}$, implying $v_{Ms}\left(\vec{r}\right)=iu_{M,-s}\left(\vec{r}\right)$. Finally, since $\mathcal{C}=\tilde{\mathcal{T}}\mathcal{P}$, the Majorana state must be an eigenstate of $\tilde{\mathcal{T}}$ as well, yielding the last relation $u_{Ms}\left(\vec{r}\right)=-iu_{M,-s}^{\star}\left(\vec{r}\right)$. The relations illustrated in Fig. 1(e,f) become equalities within the Majorana spinor.

II Spin canting of the Majorana state

In the nanowire/quantum dot setups where the character of the potential MBS is determined by analyzing its coupling to the discrete levels of a quantum dot, the spin canting of the MBS turns out to play an important role. Prada et al. (2017); Hoffman et al. (2017); Schuray et al. (2018) If there is a mismatch between the spin of the MBS and that of the electron on the quantum dot, the coupling is suppressed. Thus we turn next to examine the local spin canting angle in our Majorana nanotube.

We first notice that the total spin of the Majorana particle, summed over both particle and hole contributions, is zero. Thus, we focus on the relative spin composition of the particle components, $(u_{M\uparrow},u_{M\downarrow})$. These are complex quantities for the considered CNT setup. The local expectation value for each spin direction in the particle sector is given by $\langle\vec{u}_{M}(\vec{r})|s_{\alpha}|\vec{u}_{M}(\vec{r})\rangle$, where $s_{\alpha}$ are the Pauli matrices, $\alpha=x,y,z$, and $\vec{u}_{M}\left(\vec{r}\right)=\left(u_{M\uparrow}\left(\vec{r}\right)\text{, }u_{M\downarrow}\left(\vec{r}\right)\right)^{T}$ is the electron component of the wave function.

Due to the symmetry relations, see Fig. 1(e) and Ref. Hoffman et al., 2017, for the Majorana state it holds

[TABLE]

The expectation value $\Braket{s_{z}}$ is zero because of the pseudo time-reversal symmetry. Knowing the values of $\langle s_{x}(\vec{r})\rangle$ and $\langle s_{y}(\vec{r})\rangle$ we can define a local spin direction in the plane perpendicular to the nanotube,

[TABLE]

The full 3D spatial profile of the wave function together with the local $\theta_{xy}(\vec{r})$ for our numerically obtained Majorana state is shown in Fig. 2(a). The distance from the CNT surface encodes the local amplitude of the MBS wave function, $|u_{M\uparrow}(\vec{r})|$, and the color scale maps $\theta_{xy}(\vec{r})$. The oscillation of $\theta_{xy}$ along $z$ with the same period as the MBS wave function is clearly visible. Further, Fig. 2(b) shows a zoom of the left end of the tube for the first peak of $|u_{M\uparrow}(\vec{r})|$ along $z$, polar angle $\varphi$ resolved and displaying the helical pattern of $\theta_{xy}$. Finally, Fig. 2(c) visualizes the local spin canting at the very left end of the nanotube, where the electron tunneling would occur. The spin canting angle distribution takes several different values at the edge atoms, with visible $C_{2}$ symmetry. Thus the tunneling from a putative quantum dot coupled to the left end is definitely different than in a nanowire, assumed to be isotropic. Whether this effect is helpful or detrimental for the experiment is not yet clear.

III Effective four and one-band model

The low energy Hamiltonian of a non-superconducting CNT in the basis $\{\ket{kK\uparrow}\text{,}\ket{kK\downarrow}\text{,}\ket{kK^{\prime}\uparrow}\text{,}\ket{kK^{\prime}\downarrow}\}$ is given by

[TABLE]

where $\xi_{\tau s}\left(k\right)=\varepsilon_{\tau s}\left(k\right)-\mu$ is the single-particle energy measured with respect to the chemical potential $\mu$, $\varepsilon_{\tau s}\left(k\right)$ is the single-particle energy of the electrons (see Eq. (A.1)), $\Delta_{KK^{\prime}}$ is the energy scale associated with the valley mixing and $\mu_{B}B_{\perp}$ is the Zeeman energy due to the perpendicular magnetic field $B_{\perp}$. Diagonalization of this Hamiltonian results in four spin- and valley-mixed bands shown in Fig. 1(b). We can safely neglect any contributions from disorder, because CNTs can be grown with ultraclean lattices. Cao et al. (2005); Deshpande et al. (2009); Jung et al. (2013) The Bloch Hamiltonian can be solved analytically with the assumption that the correlation induced by the magnetic field between lower (➀,➁) and the upper (➂,➃) pairs of bands is negligible Marganska et al. (2018). When the chemical potential is set in the lower gap at the $\Gamma$-point, this approximation allows us to consider only the lower bands $\tilde{E}_{1}\left(k\right)$ and $\tilde{E}_{2}\left(k\right)$; it holds for $\mu_{B}B_{\perp}$ smaller than both of the spin-orbit coupling and the valley mixing energy scales, which in our case are $\sim 2$ meV. The details of the calculation and a short discussion of the CNT properties is presented in the Appendix A.1.

In the eigenbasis of (2) with the two-band approximation the corresponding BdG Hamiltonian for our system is given by

[TABLE]

Out of the two superconducting pairing terms, $\tilde{\Delta}_{s}\left(k\right)=\tilde{\Delta}_{s}\left(-k\right)$ is an even function of $k$, while $\tilde{\Delta}_{p}\left(k\right)=-\tilde{\Delta}_{p}\left(-k\right)$ is an odd function of $k$, see Fig. 1(d). The pairing term $\tilde{\Delta}_{p}\left(k\right)$ can be viewed as a $p$-wave like gap. The BdG Hamiltonian (3) can be partly diagonalized, taking into account the blocks with the single particle energies $\tilde{E}_{1}\left(k\right)$, $\tilde{E}_{2}\left(k\right)$ and the superconducting gap $\tilde{\Delta}_{s}\left(k\right)$. Details of this calculation are given in the Appendix A.3. Then, the rotated BdG Hamiltonian is block-diagonal and the blocks are given by

[TABLE]

The quasiparticle energies $\tilde{\xi}_{\pm}\left(k\right)$ are

[TABLE]

The functions $\tilde{\xi}_{+}\left(k\right)$ and $\tilde{\Delta}_{p}\left(k\right)$ are sketched in Fig. 3(a). The low energy physics, relevant for the Majorana states, is described by the block $\hat{\mathcal{H}}_{\text{BdG}}^{+}$. The particle-hole symmetry operator for the $\hat{\mathcal{H}}_{\text{BdG}}^{+}$ block is $\mathcal{P}=\tau_{x}\mathcal{K}$, and the chiral symmetry operator is $\mathcal{C}=\tau_{y}$, where $\tau_{x,y,z}$ are the Pauli matrices acting in the two-dimensional subspace of each block.

IV Analytical reconstruction of the 3D MBS wave function

IV.1 1D Majorana profile

Majorana bound states are zero energy eigenstates of the BdG Hamiltonian and of the particle-hole symmetry operator. From the behavior of $\tilde{\xi}_{+}\left(k\right)$ we infer that the low-energy physics has three contributions: one from the $\Gamma$-point and one from each of the Fermi points. This ansatz is confirmed by the Fourier transforms for several azimuthal cuts ($\varphi=r_{\perp}/R=$ const) of the numerically obtained MBS wave function, shown in Fig. 3(b). One clearly sees one peak at the $\Gamma$-point and two peaks at opposite momenta. The peak locations are independent of $\varphi$ but their height is not. Furthermore, the peak at negative $k$ is larger. This is caused by the helical spin structure of the single-particle spectrum, shown in Fig. 1(b). The solution at $\pm k_{F}$ is generated mostly by the band ➀, and spin $\uparrow$ for this band is associated with $k<0$.

Thus, similar to some 1D models for nanowires Klinovaja and Loss (2012), the generic form of a Majorana state can be defined as

[TABLE]

We will later take into account the 3D nature of each of these three contributions and reconstruct the 3D spatial profile of the Majorana wave function. For now we approximate $\hat{\mathcal{H}}_{\text{BdG}}^{+}\approx\hat{\mathcal{H}}_{\text{BdG}}^{\Gamma}+\hat{\mathcal{H}}_{\text{BdG}}^{R}+\hat{\mathcal{H}}_{\text{BdG}}^{L}$, where we make Taylor expansions around the momenta $k=0$ and $k=\pm k_{F}$, with $k_{F}$ determined by the constraint $\tilde{\xi}_{+}\left(k_{F}\right)=0$. The details of the calculation are presented in Appendix B.

Crucially, the spinorial components of the solutions at each of the three $k$ points are the same, which allows us to combine them into a single state which is also an eigenstate of both $\mathcal{P}$ and $\mathcal{C}$. With the three contributions we can construct the 1D solution from the generic solution (5). It is characterized by an exponential decay governed by the imaginary wave vectors $\kappa_{i}$ ($i=\Gamma,L,R$). The coefficients can be determined by the three constraints

[TABLE]

From previous findings Marganska et al. (2018) we know that in the topological regime $\kappa_{\Gamma}\in\mathds{R}$ and $\kappa_{R}\text{, }\kappa_{L}\in\mathds{C}$. Moreover, it holds that $\text{Re}\left(\kappa_{R}\right)=\text{Re}\left(\kappa_{L}\right)$ and $\text{Im}\left(\kappa_{R}\right)=-\text{Im}\left(\kappa_{L}\right)$ $\Leftrightarrow\kappa_{R}=\kappa_{L}^{\star}$. Therefore, the wave function can be written as

[TABLE]

These eigenvectors are not eigenstates of the particle-hole operator $\mathcal{P}=\tau_{x}\mathcal{K}$, but we can multiply them by a complex number $c_{\pm}=\pm 1+i$ , such that they satify the Majorana constraint. Then, by applying the Majorana  (6a) and the boundary (6b) conditions we get the 1D solution, which is given by

[TABLE]

where

[TABLE]

encodes the dependence of the wave function on the longitudinal coordinate. The sum $\psi_{\parallel}(z)+\psi^{*}_{\parallel}(z)$ satisfies the boundary condition (6b), and $\mathcal{N}$ is the normalization constant determined from (6c). The contribution from the $\Gamma$-point is a pure evanescent state and from the contribution from the Fermi points we get a decaying oscillation with the wavevector $k_{F}$.

IV.2 Reconstructing the 3D profile

In the remaining part of this work we will provide the analytical form only for $u_{\uparrow}\left(\vec{r}\right)$ (dropping the $M$ subscript for compactness of notation), since the remaining Majorana spinor components can be obtained by the application of $\mathcal{P}$, $\tilde{\mathcal{T}}$ and $\mathcal{C}$ symmetries.

The Majorana operator to create the state (5) is defined as

[TABLE]

where $k\in\{\Gamma\text{, }k_{F}\text{, }-k_{F}\}$ and $u=v^{\star}=\frac{1-i}{\sqrt{2}}$. In order to find the analytical wave function we need to transform the wave function from the one-band back to the four-band model; this procedure is discussed in Appendix C. To express the Majorana state in the sublattice- and spin-resolved basis we need the transformations reversing (20), (25) and (31). At the end we obtain

[TABLE]

for $k\in\{\Gamma,\pm k_{F}\}$, where the coefficients $u_{\tau s}\left(k\right)$ correspond to the electron and $v_{\tau s}\left(k\right)$ to the hole contribution, respectively. We find a compact form for the coefficients

[TABLE]

with

[TABLE]

and (see Eq. (21) for $a_{s}(k)$ and $b_{s}(k)$)

[TABLE]

The coefficients $g(k),h(k)$ and $n(k),m(k)$ are found below, in Eqs. (26) and (32), respectively. By using the relations $a_{s}\left(-k\right)=b_{-s}\left(k\right)$, $g\left(-k\right)=h\left(k\right)$, $m\left(-k\right)=m\left(k\right)$, $n\left(-k\right)=n\left(k\right)$, we obtain $\Lambda_{\tau s}\left(k\right)=\Lambda_{-\tau-s}\left(-k\right)$ and $\lambda_{s}\left(k\right)=-i\lambda_{-s}^{\star}\left(-k\right)$. Finally, we arrive at the symmetry relations of the electron and hole coefficients $u_{\tau s}\left(k\right)$ and $v_{\tau s}\left(k\right)$ illustrated in Fig. 1(f).

We have now the expression of the wave function in conduction basis. In order to apply the boundary condition it must however be recast in the sublattice-resolved basis. In general for the transformation into the sublattice basis one needs also the valence band contribution. Here we can use the fact that, due to the high chemical potential, we are far away from the charge neutrality point and therefore the contribution from the valence band is negligible. With this the components in the sublattice basis are defined as $u_{p\tau s}\left(k\right)=e^{ip\eta_{\tau sk}}u_{\tau s}\left(k\right)$, where $\eta_{\tau s}\left(k\right)=\text{arg}\left(\gamma_{\tau s}\left(k\right)\right)$ is the phase of (15) in the low-energy regime, and $p=+1$ for $A$ sublattice and $p=-1$ for $B$ sublattice.

Since our nanotube is chiral, the open boundary conditions imply that the wave function must vanish on one end at the missing $A$ atoms and on the other end at the missing $B$ atomsdel Valle et al. (2011). We use therefore the open boundary condition $\Psi_{A}\left(z=0\text{, }r_{\perp}\right)\overset{!}{=}0$ $\forall r_{\perp}$. The wave function $u_{p\uparrow}\left(\vec{r}\right)$ is given by the superposition of the three contributions $k\in\{\Gamma\text{, }k_{F}\text{, }-k_{F}\}$ and the two valleys $K$ and $K^{\prime}$, each with its specific transverse profile $e^{i\tau k_{\perp}x_{\perp}}$:

[TABLE]

The amplitudes can be fixed by observing that the Majorana condition requires $A_{\Gamma}\in\mathds{R}$ and $A_{R}=A_{L}^{\star}$. From the open boundary condition in longitudinal direction we obtain a relation between $A_{R}$ and $A_{\Gamma}$; hence the particle component of the wave function can be written as

[TABLE]

The expressions for $u_{\tau s}(k)$ are given in Eq. (10), and for $\psi_{\parallel}(z)$ in Eq. (8). The spatial profile of the wave function is not trivial, in the sense that it cannot be factorized into separate longitudinal and transverse profiles, $u_{p\uparrow}\left(\vec{r}\right)\neq f\left(r_{\perp}\right)g\left(z\right)$. The absolute value $|A_{R}|$ is fixed by the normalization and its phase by the Majorana condition. Note that the transverse momentum $k_{\perp}$ is quantized by the periodic boundary condition. The Fermi wavevector $k_{F}$ is given by the position of the chemical potential $\mu$, and the characteristic decay lengths at $\Gamma$ and $\pm k_{F}$ by the parameters of the Hamiltonian at this $\mu$. Thus all factors in the wave function are in principle known from the analytics.

IV.3 Comparison between analytical and numerical results

In order to test the accuracy of our formula Eq. (13), we have performed a comparison between the analytical and numerical solutions for several 1D cuts of the full MBS profile, at varying values of the azimuthal angle $\varphi$. We fitted the numerical solutions with (13), finding for each cut the parameters $\kappa_{\Gamma},\kappa_{F},k_{F}$ and $A_{R}$.

The results for three values of the polar angle, $\varphi=0^{\circ},24.23^{\circ},114.23^{\circ}$ are shown in Fig. 4. The analytical model clearly reproduces very well the numerically obtained wave functions. However, due to the simplifications inherent in the effective one-band model, there are three aspects where we have to adjust for the lost information.

(i) In the microscopic model the $\mathcal{P}$ symmetry holds exactly (by construction), but $\tilde{\mathcal{T}}$ is minimally broken by two small effects. One is the presence of the weak spin-flip terms in the Hamiltonian, due to the enhanced spin-orbit coupling Ando (2000); Izumida et al. (2009); del Valle et al. (2011). The other is the small Peierls phase for the nearest neighbor hopping, due to the magnetic field Peierls (1933). Thus in the numerical solution the $\tilde{\mathcal{T}}$- and $\mathcal{C}$-related components of the Nambu spinor differ by about $\pm 3\%$. Removing the spin-flip and the Peierls phase restores the $\tilde{\mathcal{T}}$ and consequently also the $\mathcal{C}$ symmetries, see Appendix A.2 for details.

(ii) In the analytics we neglected some correlations due to the magnetic field. Further, we performed Taylor expansions around the three momenta $k=0,\pm k_{F}$. Thus, the values $\kappa_{\Gamma}$, $\kappa_{F}$ and $k_{F}$ from the analytics are slightly different from those which are obtained by fitting the numerical data using (13), see Tab. 1.

(iii) We implemented the valley mixing through a continuous potential ridge along the CNT/superconductor interface. This results in the coupling between the two valleys, but also in their coupling to higher transverse momentum bands which therefore also contribute, albeit very weakly, to the final Majorana state. In consequence, although we expect $A_{R}$ to be independent of $\varphi$, we obtain from the fitting procedure different $A_{R}$ for different $\varphi$ cuts, with the resulting values of $|A_{R}(\varphi)|$ shown in Fig. 5. We see that, although not constant, the amplitude $A_{R}$ is a weakly varying function of $\varphi$. Moreover, the data resolved for atoms at the same $z$ position show that $A_{R}$ is close to $\pi$-periodic. This is a consequence of the $C_{4}$ symmetry of our (12,4) CNT where the $K^{\prime}/K$ valley states carry the angular momentum $\ell=\pm 1$. Since the Majorana state is constructed predominantly from electron (and hole) states with $\ell=\pm 1$, the amplitude of its wave function, to which $A_{R}$ was fitted, has an approximate $C_{2}$ symmetry. This is also visible in Fig. 2(c), where the $C_{2}$ (instead of $C_{1}$) symmetry of spin texture arises from the factor of 2 in Eq. (1).

In Fig. 6 we show a comparison between the analytical and numerical results for $\mathrm{Re}(u_{A\uparrow}),\mathrm{Im}(u_{A\uparrow})$ and the resulting canting angle $\theta_{xy}(z)$ for $\varphi=0$. The slight discrepancy between the numerical and analytical values of the real and imaginary part of $u_{\uparrow}(\vec{r})$, shown in Fig. 6(a-b), is amplified in the spin canting angle behavior shown in Fig. 6(c). In particular, additional phase jumps are visible at positions where the real value in numerics is small and positive, while the analytical result is also small but negative. Nevertheless, the overall agreement is again good.

Conclusion

In this work we have shown in a combination of numerical modelling and analytical calculations how to determine the full spatial profile of the Majorana bound state in a proximitized semiconducting carbon nanotube. The wave function has three contributions: one from the $\Gamma$-point and one from each Fermi point, which is also supported by an analysis of the numerical data via a Fourier transformation. We find the symmetry relations which must be fulfilled by the components of the Majorana spinor. The excellent agreement between the analytically obtained and the numerically calculated spin and sublattice resolved spinor gives us confidence in the accuracy of the local observables further derived in this work. Despite being obtained for a CNT, our results might serve as a reference also for other systems where a microscopic calculation of the MBS spinor is not possible. The features which our model captures very well are: the three main momentum contributions to the MBS, the decaying behaviour of the wave function combined with its spiral pattern, its oscillation and the symmetries linking the different components of the Nambu spinor. We show that our analytical model fits very well the numerical data of the wave function obtained by a tight-binding calculation. Our results will be useful for modeling and interpreting the experimental results in a realistic quantum transport setup where the properties of the Majorana states are probed locally.

Acknowledgements.

The authors thank the Deutsche Forschungsgemeinschaft for financial support via GRK 1570 and IGK “Topological insulators” grants, as well as the JSPS for the KAKENHI Grants (Nos. JP15K05118, JP15KK0147, JP18H04282).

Appendix A CNT Spectrum

A.1 Single-particle spectrum

The simplest way of obtaining the Hamiltonian of a CNT in the momentum space representation is to use the zone folding approximationSaito et al. (1998). The Hamiltonian of a CNT can be written in the sublattice basis for $A$ and $B$ sublattice

[TABLE]

where $a_{\mathbf{k}s}^{\dagger}$ ($b_{\mathbf{k}s}^{\dagger}$) creates an electron on $A$ ($B$) sublattice with momentum $\mathbf{k}$ and spin $s$. The kinetic energy $\gamma\left(\mathbf{k}\right)$ is defined as

[TABLE]

where $t_{s,i}$ is the spin-dependent hopping parameter between an $A$ atom and its $i$-th neighbor, and $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ are the Bravais lattice vectors of the graphene lattice, see Fig. 1(a). The low-energy unperturbed CNT Hamiltonian $H_{\text{CNT}}$ can be obtained by an expansion of (14) around the Dirac points $\mathbf{k}=\bm{\kappa}+\tau\mathbf{K}$del Valle et al. (2011) and a rotation from sublattice into conduction/valence band basis. In the following we will assume that the chemical potential is in the conduction band, obtaining

[TABLE]

where $\xi_{\tau s}\left(k\right)=\varepsilon_{\tau s}\left(k\right)-\mu$, $\varepsilon_{\tau s}\left(k\right)=\left|\gamma_{\tau s}\left(k\right)\right|$ is the CNT single-particle energy in the conduction band, $\mu$ the chemical potential and $c_{k\tau s}^{\dagger}\Ket{0}=\Ket{k\tau s}$ define the basis of (2). The curvature of the CNT’s lattice results in both spin-dependent and spin-independent modifications, i.e. shifts in both transverse and longitudinal momentum. Thus, the single-particle energies of a CNT (2) for given transverse momentum $k_{\perp}$ and longitudinal momentum $k$ at low energies are given by

[TABLE]

where $K_{\perp},K_{\parallel}$ are the transverse and longitudinal components of momentum at the Dirac point $K$. The quantum numbers $\tau$ and $s$ are are defined in the main text. In the case of the (12,4) semiconducting nanotube, the numerical values of those momentum shifts in our calculations are $\Delta k_{\perp}^{c}=-22.83\,\mu\mathrm{m}^{-1}$, $\Delta k_{\parallel}^{c}=66.62\,\mu\mathrm{m}^{-1}$, $\Delta k_{\mathrm{SO}}=-2.917\,\mu\mathrm{m}^{-1}$, $K_{\parallel}=0$, and the lowest energy subbands shown in Fig. 1(b) have $k_{\perp}-\tau K_{\perp}=\tau/3R$. The value of $k_{\perp}$ for the $K$ valley subband in our nanotube is $-35/R$, where $R$ is the CNT radius. Note that the single-particle energies satisfy the time-reversal conjugation, $\varepsilon_{\tau s}\left(k\right)=\varepsilon_{-\tau-s}\left(-k\right)$.

The low-energy Bloch Hamiltonian (2) contains also the valley mixing and Zeeman field contributions, $H=H_{\text{CNT}}+H_{\Delta_{KK^{\prime}}}+H_{Z}$. Since the nanotube we are studying is of the zigzag class Marganska et al. (2015); Izumida et al. (2016), in order to mix the valleys it is enough to break only the rotational symmetry. In our setup we consider the valley mixing introduced by the presence of the substrate, modelling it as an electrostatic potential with a Gaussian distribution in the polar coordinate, $V(\varphi)=V_{0}\exp(-(\varphi-\pi/2)/\Delta\varphi)$. This corresponds to the substrate extending in the $xz$ plane. In the reciprocal space the valley-mixing term is given by

[TABLE]

and couples states with the same spin and $k$ but opposite valley. Following Ref. Marganska et al., 2018, we set $\Delta_{KK^{\prime}}=2.5$ meV.

The Zeeman effect with the field aplied along the $x$ axis couples opposite spins in the same valley,

[TABLE]

The CNT Hamiltonian (2) can be brought to a diagonal form by employing two unitary transformations. More details about the transformations can be found in Appendix D.1 of Ref. [Marganska et al., 2018]. The first transformation diagonalizes the Hamiltonian without Zeeman energy ($B_{\perp}=0$) and is defined as

[TABLE]

with $a_{s}\left(k\right)^{2}+b_{s}\left(k\right)^{2}=1$ and the following values of $a_{s}\left(k\right)$ and $b_{s}\left(k\right)$,

[TABLE]

where the energy eigenvalues are

[TABLE]

Due to the time-reversal conjugation of $\xi_{\tau s}\left(k\right)=\xi_{-\tau-s}\left(-k\right)$, it can be shown that $a_{s}\left(k\right)=b_{-s}\left(-k\right)$ and $E_{\pm s}\left(k\right)=E_{\pm-s}\left(-k\right)$.

Using equations (20) the Zeeman term can be expressed as

[TABLE]

The magnetic field $\tilde{B}_{\perp}$ couples the spins within the lower and upper band pair, while $B_{\perp}^{\star}$ couples the spins between band pairs. Both are symmetric in $k$, i.e. $\tilde{B}_{\perp}\left(k\right)=\tilde{B}_{\perp}\left(-k\right)$ and $B_{\perp}^{\star}\left(k\right)=B_{\perp}^{\star}\left(-k\right)$. This is a consequence of the pseudo-time reversal symmetry.

In the regime of small Zeeman energy, i.e. $\Delta E=\left|E_{+s}-E_{-s}\right|>\mu_{B}B_{\perp}$, the terms with $B_{\perp}^{\star}$ can be omitted. This allows us to treat the upper and lower pair of bands separately. We shall proceed to find the solutions for the lower band pair only, assuming that the chemical potential $\mu$ is tuned into the gap between the two energy bands $\tilde{E}_{1}$ and $\tilde{E}_{2}$. Therefore, we will neglect the influence of the bands $\tilde{E}_{3}$ and $\tilde{E}_{4}$ because those bands are not occupied. Then, the second transformation diagonalizing the Hamiltonian with magnetic field is defined as

[TABLE]

where the coefficients must satify $g^{2}\left(h\right)+h^{2}\left(k\right)=1$. The new quantum number in (25) $i\in\{1\text{,}2\}$ just reflects the ordering of the energy bands $E_{1}<E_{2}$. The coefficients $g$ and $h$ are defined as

[TABLE]

The coefficients satisfy the pseudo-time-reversal conjugation $g\left(k\right)=h\left(-k\right)$. Then, the single-particle energies of the full Hamiltonian with decoupled band pairs are

[TABLE]

The renormalized magnetic field opens a band gap at the $\Gamma$-point. The single-particle energies have the property $\tilde{E}_{i}\left(k\right)=\tilde{E}_{i}\left(-k\right)$ with $i\in\{1\text{,}2\}$ because $\tilde{B}_{\perp}\left(k\right)=\tilde{B}_{\perp}\left(-k\right)$. This pseudo-time reversal symmetry for conduction band states results in the relation depicted in Fig. 1(f). Since the single-particle states of a finite CNT in our setup contain both $\vec{k},s$ and $-\vec{k},-s$ contributions with equal weights, their spin components in the real space must also obey the relation shown in Fig. 1(e).

A.2 Pseudo-time reversal symmetry

The pseudo-time reversal invariance holds exactly for our effective model Hamiltonian  (2). For the real space Hamiltonian it is however broken by two effects, both absent in our four-band model. The first and smaller one is the presence of the Peierls phase.Peierls (1933) This phase can be safely neglected - a magnetic field of 120 T would result in only $10^{-3}$ of a flux quantum per each hexagonal plaquette.

The second and more important effect is the presence of nearest-neighbor hoppings with spin flip, which couple neighboring angular momentum subbands.Ando (2000); Izumida et al. (2009); del Valle et al. (2011) Including it would require bringing the number of subbands up to twelve. Figure 7 shows the strength of $\tilde{\mathcal{T}}$ breaking, quantified as the difference between the band ➀ (cf. Fig. 1) minima at $k>0$ and at $k<0$, as a function of $B$ and of $V_{0}$, as shown in the inset of Fig 7(a). Although the Peierls phase does contribute to the breaking of $\tilde{\mathcal{T}}$ when the spin flips are included, we see that neglecting the spin-fliping hoppings restores $\tilde{\mathcal{T}}$ completely. Nevertheless, even when all effects are present, for our parameters $B=14$ T and $V_{0}=0.4$ eV the pseudo-time reversal still holds down to $\mu$eV energy scales.

A.3 Superconducting spectrum

Including superconducting correlations on a mean-field level we add to the Hamiltonian a superconducting pairing term Uchoa and Castro Neto (2007), which is given by

[TABLE]

where $\Delta_{0}$ is the superconducting order parameter, which we take to be 0.4 meV. We can express the pairing Hamiltonian (28) in the eigenbasis of the CNT (2) and, after applying the approximations and transformations described in Appendix A.1, we obtain the BdG Hamiltonian (3) with the pairing terms

[TABLE]

We see that the pairing term $\tilde{\Delta}_{s}\left(k\right)$ has an even and $\tilde{\Delta}_{p}\left(k\right)$ an odd parity, as shown in Fig. 1(c).

The basis change which transforms (3) into (4) is given by

[TABLE]

with the normalization condition $m^{2}\left(k\right)+n^{2}\left(k\right)=1$ and the coefficients defined in the following way:

[TABLE]

Appendix B 1D Majorana bound state solutions

Majorana bound states are zero energy eigenstates of the BdG Hamiltonian and the particle-hole symmetry operator. The low-energy physics of the BdG Hamiltonian (4) is described by the block $\hat{\mathcal{H}}_{\text{BdG}}^{+}$. For the Majorana bound states we will approximate the BdG Hamiltonian by $\hat{\mathcal{H}}_{\text{BdG}}^{+}\approx\hat{\mathcal{H}}_{\text{BdG}}^{\Gamma}+\hat{\mathcal{H}}_{\text{BdG}}^{R}+\hat{\mathcal{H}}_{\text{BdG}}^{L}$ because the low-energy physics of $\hat{\mathcal{H}}_{\text{BdG}}^{+}$ has three contributions, as illustrated in Fig. 3(a). This is also supported by the numerics, see Fig. 3(b).

B.1 $\Gamma$-point contribution

The first contribution is coming from the $\Gamma$-point. Therefore, we obtain from a Taylor expansion around the $\Gamma$-point

[TABLE]

where $\frac{1}{m^{\star}}=\left.\frac{\partial^{2}\tilde{\xi}_{+}\left(k\right)}{\hbar^{2}\partial k^{2}}\right|_{k=0}$ and $\lambda=\left.\frac{\partial\tilde{\Delta}_{p}\left(k\right)}{\hbar\partial k}\right|_{k=0}$. Then the BdG Hamiltonian for $k\approx\Gamma$ becomes

[TABLE]

and the corresponding BdG equation reads

[TABLE]

Now, we interpret $k$ as the momentum operator $k\to\hat{k}=-i\partial_{z}$ and make the ansatz

[TABLE]

For the momentum $\kappa_{\Gamma}$ we need to solve the secular equation $\det\left(\hat{\mathcal{H}}_{\text{BdG}}^{\Gamma}-E\mathds{1}\right)\overset{!}{=}0$ for any energy $E$ and we obtain

[TABLE]

For zero energy modes the equation can be simplified

[TABLE]

The corresponding zero energy eigenvectors are given by

[TABLE]

B.2 Fermi point contribution

For the Fermi point contribution we need to linearize $\tilde{\xi}_{+}\left(k\right)$ around $k_{F}$ and $-k_{F}$, see Fig. 3. Then, we can define the following two Nambu spinors $\Psi_{R}=\left(d_{k+,R}\text{, }d_{-k+,L}^{\dagger}\right)$ and $\Psi_{L}=\left(d_{k+,L}\text{, }d_{-k+,R}^{\dagger}\right)$. The subscripts $R,L$ denote the right- and left-movers. The corresponding BdG Hamiltonians are given by

[TABLE]

[TABLE]

where for $\hat{\mathcal{H}}_{\text{BdG}}^{R}$ we have $k>0$ and for $\hat{\mathcal{H}}_{\text{BdG}}^{L}$ we have $k<0$. The corresponding BdG equation reads

[TABLE]

[TABLE]

With $k\to\hat{k}=-i\partial_{z}$ and making the ansatz

[TABLE]

we get the decay lengths $\kappa_{R}$ and $\kappa_{L}$ from the secular equations $\det\left(\hat{\mathcal{H}}_{\text{BdG}}^{\text{R/L}}-E\mathds{1}\right)\overset{!}{=}0$. The decay lengths for the zero energy modes become $\kappa_{R}=ik_{F}\mp\frac{\left|\tilde{\Delta}_{p}\left(k_{F}\right)\right|}{v_{F}\hbar}$ and $\kappa_{L}=-ik_{F}\mp\frac{\left|\tilde{\Delta}_{p}\left(k_{F}\right)\right|}{v_{F}\hbar}$. Furthermore, we get the two eigenvectors

[TABLE]

where we used $\text{sgn}\left(\tilde{\Delta}_{p}\left(k_{F}\right)\right)=-1$ and $\text{sgn}\left(\tilde{\Delta}_{p}\left(-k_{F}\right)\right)=+1$, see Fig. 1(b).

Appendix C Construction of 3D Majorana wave function

Explicitly, the coefficients of the electron and holes are parts of the Majorana bound state  (9) given by

[TABLE]