Majorana Zero-Modes in a Two-Dimensional $p$-Wave Superconductor

TL;DR

This paper investigates the emergence and manipulation of Majorana zero-modes at the edges of a two-dimensional topological superconductor, focusing on a triplet phase in doped Bi$_2$Se$_3$ and their potential use in quantum computing.

Contribution

It demonstrates how in-plane magnetic fields can create and control Majorana modes in a 2D topological superconductor, proposing new methods for their braiding.

Findings

Magnetic field opens a gap in surface modes, inducing Majorana zero-modes.

Majorana modes can be manipulated by changing the magnetic field orientation.

Potential applications in topological quantum computation.

Abstract

We analyze the formation of Majorana zero-modes at the edge of a two-dimensional topological superconductor. In particular, we study a time-reversal-invariant triplet phase that is likely to exist in doped Bi$_2$Se$_3$. Upon the introduction of an in-plane magnetic field to the superconductor, a gap is opened in the surface modes, which induces localized Majorana modes. The position of these modes can be simply manipulated by changing the orientation of the applied field, yielding novel methods for braiding these states with possible applications to topological quantum computation.

Tables

Table 1. Table 1: Relevant parameters, analytical expressions, and numerical estimates of different quantities discussed in the text. Here, k B subscript 𝑘 𝐵 k_{B} is Boltzmann’s constant, and μ e subscript 𝜇 𝑒 \mu_{e} is the magnetic moment of the electron. We estimate these parameters based on experimental data available for Cu x Bi 2 Se 3 Hor  et al. ( 2010 ); Wray  et al. ( 2010 ); Kriener  et al. ( 2011 ); Zhang  et al. ( 2009b ); Sandilands  et al. ( 2014 ); Wray  et al. ( 2011 ) . Details can be found in the Appendix.

Input parameters:	Radius	$R$		$1\mu m$
	Magnetic field	$B_Z$		1 T
Length and angular scales:	Fermi wavelength	$k_F$		${10}^{-1}{\AA }^{-1}$
	Electron density per layer	${\rho }_L$	$k_F^2/2\pi$	$1.6\times {10}^{13}{\mathrm{cm}}^{-2}$
	Coherence length	$\xi$	$0.2\times \left(\hslash v_F\right)/\left(k_B{T}_c\right)$	$2\times {10}^3\AA$
	Angular width of Majorana states	$1/b$	$\sqrt{{\Delta }_{\text{exp}}/\left({\mu }_e{B}_Z{k}_{F}R\right)}$	0.32 radians
Energy scales:	Fermi energy	${\overline{ϵ}}_F$	$v_F^2{k}_F^2/2E_G$	$k_B\times 2300K=200\mathrm{meV}$
	Critical temperature	$T_c$		3.8 K
	Superconducting gap	${\Delta }_{\text{exp}}$	$1.76k_B{T}_c$	$k_B\times 6.7K=0.6\mathrm{meV}$
	Quasi-particle gap	${\Delta }_M$	$\sqrt{\left(2{\mu }_e{B}_Z{\Delta }_{\text{exp}}\right)/\left(k_{F}R\right)}$	$k_B\times 0.3K=0.026\mathrm{meV}$
	Gap between Majorana modes	$\delta {ϵ}_M$	${\Delta }_M\exp \left(-2b^2\right)/\sqrt{2}$	$k_B\times {10}^{-8}K={10}^{-9}\mathrm{meV}$

Equations

$${\mathcal{H}}^{\text{2D}}(\mathbf{k})=E_G{\sigma }_x+v_F{\sigma }_z\left(k_x{s}_y-k_y{s}_x\right),$$

$$\mathcal{H}(\mathbf{k})=[{\mathcal{H}}^{\text{2D}}(\mathbf{k})-{ϵ}_F]{\tau }_z+{\Delta }_{\text{sc}}{\sigma }_y{s}_z{\tau }_x,$$

$$\mathcal{H}_{0}(\mathbf{k})=\left(\begin{array}[]{ccc}\mathcal{H}_{+}(\mathbf{k})&0\\ 0&\mathcal{H}_{-}(\mathbf{k})\end{array}\right),$$

$$\mathcal{H}_{\pm}(\boldsymbol{\kappa})=\left(\begin{array}[]{ccc}\pm\left(|\boldsymbol{\kappa}|^{2}-1\right)&\sqrt{2\gamma}\left(-i\kappa_{x}+\kappa_{y}\right)\\ \sqrt{2\gamma}\left(i\kappa_{x}+\kappa_{y}\right)&\pm\left(-|\boldsymbol{\kappa}|^{2}+1\right)\end{array}\right),$$

$${\mathcal{E}}_{\pm }^{\text{edge}}\approx \mp \left(\ell +\frac{1}{2}\right)\frac{\sqrt{2\gamma }}{\lambda }.$$

$${\mathcal{H}}_M$$

$${\mathcal{H}}^{\text{3D}}(\mathbf{k})$$

$$E(\mathbf{k})$$

$$|\Psi ⟩$$

$$\mathcal{H}(\mathbf{k})=({\mathcal{H}}^{\text{3D}}(\mathbf{k})-{ϵ}_F){\tau }_z,$$

$${\mathcal{H}}^{\text{red}}(\mathbf{k})$$

$$-{\mathcal{P}}_{+}\left\{v_F\left[\left(k_x{s}_y-k_y{s}_x\right){\sigma }_x{\tau }_z+v_z{k}_z{\sigma }_y\right]{\left(E_G{\sigma }_z{\tau }_z\right)}^{-1}{\Delta }_{\text{sc}}{\sigma }_y{s}_z{\tau }_x\right\}{\mathcal{P}}_{+}$$

$$-{\mathcal{P}}_{+}\left\{{\Delta }_{\text{sc}}{\sigma }_y{s}_z{\tau }_x{\left(E_G{\sigma }_z{\tau }_z\right)}^{-1}\left[v_F\left(k_x{s}_y-k_y{s}_x\right){\sigma }_x{\tau }_z+v_z{k}_z{\sigma }_y\right]\right\}{\mathcal{P}}_{+}.$$

$${\mathcal{H}}^{\text{red}}(\mathbf{k})=\left[\left(\frac{v_F^2|\mathbf{k}{|}^2+v_z^2{k}_z^2}{2E_G}-{\overline{ϵ}}_F\right){\tau }_z-\frac{2v_F{\Delta }_{\text{sc}}}{E_G}{\tau }_x\left(k_x{s}_x+k_y{s}_y\right)-\frac{2v_z{\Delta }_{\text{sc}}}{E_G}{\tau }_x{k}_z{s}_z\right]{\mathcal{P}}_{+},$$

$$k_{Fz}$$

$${\xi }_z^{-1}$$

$$\delta {ϵ}_z^n$$

$${\mathcal{H}}^{\text{2D}}(\mathbf{k})$$

$$\mathcal{H}_{z}=\bar{\epsilon}_{F}\left(\begin{array}[]{ccc}\frac{k_{z}^{2}}{k_{Fz}^{2}}-1&-\frac{2k_{z}}{\xi_{z}k_{Fz}^{2}}\\ -\frac{2k_{z}}{\xi_{z}k_{Fz}^{2}}&-\frac{k_{z}^{2}}{k_{Fz}^{2}}+1\end{array}\right).$$

$${\left(\frac{{ϵ}_0}{{\overline{ϵ}}_F}\right)}^2={\left(\frac{k_z^2}{k_{Fz}^2}-1\right)}^2+{\left(\frac{2k_z}{{\xi }_z{k}_{Fz}^2}\right)}^2.$$

$${\mathcal{H}}_{\pm }^{\text{red}}(\mathbf{k})$$

$$E(\mathbf{k})$$

$$|\Psi\rangle=\left(\begin{array}[]{cccc}\psi^{(e)}_{\uparrow}&\psi^{(h)}_{\downarrow}&-\psi^{(h)}_{\uparrow}&\psi^{(e)}_{\downarrow}\\ \end{array}\right)^{T},$$

$$\mathcal{H}_{0}(\mathbf{k})=\left(\begin{array}[]{ccc}\mathcal{H}_{+}(\mathbf{k})&0\\ 0&\mathcal{H}_{-}(\mathbf{k})\end{array}\right),$$

$${\mathcal{H}}_{+}(\mathbf{k})$$

$$\Psi_{\ell}(r,\theta)=e^{i\ell\theta}\left(\begin{array}[]{ccc}e^{i\theta}\psi_{+,\ell}(r)\\ \psi_{-,\ell}(r)\end{array}\right).$$

$$L_{z}=-i\partial_{\theta}\mathcal{I}-\frac{1}{2}s_{z}=\left(\begin{array}[]{ccc}\frac{\partial_{\theta}}{i}-\frac{1}{2}&0\\ 0&\frac{\partial_{\theta}}{i}+\frac{1}{2}\end{array}\right),$$

$$\left(\begin{array}[]{ccc}-\frac{v_{F}^{2}}{2E_{G}}\left(\partial_{r}^{2}+\frac{\partial_{r}}{r}-\frac{(\ell+1)^{2}}{r^{2}}\right)-\bar{\epsilon}_{F}&-\frac{2v_{F}\Delta_{\text{sc}}}{E_{G}}\left(\partial_{r}-\frac{\ell}{r}\right)\\ \frac{2v_{F}\Delta_{\text{sc}}}{E_{G}}\left(\partial_{r}+\frac{(\ell+1)}{r}\right)&\frac{v_{F}^{2}}{2E_{G}}\left(\partial_{r}^{2}+\frac{\partial_{r}}{r}-\frac{\ell^{2}}{r^{2}}\right)+\bar{\epsilon}_{F}\end{array}\right)\left(\begin{array}[]{ccc}\psi_{+,\ell}(r)\\ \psi_{-,\ell}(r)\end{array}\right)=E\left(\begin{array}[]{ccc}\psi_{+,\ell}(r)\\ \psi_{-,\ell}(r)\end{array}\right).$$

$${\psi }_{+,\ell }(r)=a{J}_{\ell +1}(kr),\text{and}{\psi }_{-,\ell }(r)=b{J}_{\ell }(kr),$$

$$\left(\begin{array}[]{ccc}\frac{v_{F}^{2}}{2E_{G}}k^{2}-\bar{\epsilon}_{F}&\frac{2v_{F}\Delta_{\text{sc}}}{E_{G}}k\\ \frac{2v_{F}\Delta_{\text{sc}}}{E_{G}}k&-\frac{v_{F}^{2}}{2E_{G}}k^{2}+\bar{\epsilon}_{F}\end{array}\right)\left(\begin{array}[]{ccc}a\\ b\end{array}\right)=E\left(\begin{array}[]{ccc}a\\ b\end{array}\right).$$

$$\Psi_{\ell}(r,\theta)=e^{i\ell\theta}\left(\begin{array}[]{ccc}-\frac{2v_{F}\Delta_{\text{sc}}}{E_{G}}ke^{i\theta}J_{\ell+1}(kr)\\ \left(\frac{v_{F}^{2}}{2E_{G}}k^{2}-\bar{\epsilon}_{F}-E\right)J_{\ell}(kr)\end{array}\right),$$

$$E=\pm \sqrt{{\left(\frac{v_F^2}{2E_G}{k}^2-{\overline{ϵ}}_F\right)}^2+{\left(\frac{2v_F{\Delta }_{\text{sc}}}{E_G}k\right)}^2}.$$

Full text

Many approaches to topological quantum computation are based on the creation and manipulation of massless Majorana states Kitaev (2001). These arise as excitations in a two-dimensional system when a fermion is effectively split into two parts, with each part localized far away from the other in space. Since fermions are fundamental particles, Majorana states are always generated in pairs. Such states are known to occur in topological superconductors Hasan and Kane (2010); Qi and Zhang (2011); Alicea (2012); Leijnse and Flensberg (2012); Ando and Fu (2015); Aasen et al. (2016); Kallin and Berlinsky (2016); Sato and Ando (2016) and have been predicted to exist in the $\nu=5/2$ fractional quantum Hall effect Read and Green (2000); Fendley et al. (2007). Because of their nontrivial half-fermion statistics, braiding, or exchanging, Majorana states is a non-Abelian process which takes place within the space of degenerate ground states. Quantum gates can be built by simply braiding Majorana states Alicea (2012). Currently, realistic schemes for braiding Majorana states require tri-wire junctions Alicea et al. (2011), and either electrostatic gates Sau et al. (2011) or controllable magnetic fluxes van Heck et al. (2012).

In this Letter, we propose and study a solid-state platform for the creation and manipulation of Majorana states, motivated by current research on natural topological superconductors. Our proposed system is modeled on doped Bi2Se3, but the results we present here are more general, and can be extended to other superconductors and to the superfluid B phase of 3He. In what follows, we describe the general features of the our Majorana platform, and show by detailed calculations how these Majorana states arise. We then discuss a simple scheme for manipulating Majorana excitations on superconducting discs. Using a simple extension, we show that one can study our model in arbitrary geometries, providing further impetus for experimental realization. This is described in full detail in the Appendix.

We start with a thin disc of a topological superconductor material with thickness smaller than the superconducting coherence length, as sketched in Fig. 1. In this limit, the surface states at the top and bottom surfaces hybridize and develop a small gap. The only remaining sub-gap states are quasi-one-dimensional gapless Majorana bands localized at the edges. If time-reversal symmetry (TRS) is not broken in the superconducting phase, there are two bands related by TRS inside the gap. An in-plane magnetic field applied to the system, as shown in Fig. 1, breaks time-reversal symmetry. This field hybridizes the two sets of Andreev states, and opens a gap in the energy spectrum. The sign and magnitude of the gap are determined by the normal of the field to the edge and is opposite where the field enters the disc to where the field exits the system. This leads to the formation of Majorana modes located inside the bulk energy gap. These Majorana edge modes are localized near the boundary of the disc, at the points where the field is parallel to the edge.

Of the different order parameters that have been proposed to describe the superconducting phase of doped Bi2Se3 Fu and Berg (2010); Venderbos et al. (2016); Chirolli et al. (2016), we will consider a time-reversal-invariant, odd-parity, triplet phase. Our analysis applies to similar Majorana modes appear in neutral fermionic superfluids, such as 3He Levitin et al. (2013, 2014); Saitoh et al. (2016) (see also Tsutsumi et al. (2011)), as well as to gapped atomic Fermi superfluids Ku et al. (2012). In general, the two requirements for our proposed system are: (i) (effective) two-dimensionality, and (ii) existence of gapless counter-propagating Andreev edge modes. The first condition requires the thickness of the system to be smaller than the superconducting coherence length. The second condition excludes a two-dimensional, gapped chiral $p_{x}\pm ip_{y}$ superconductor, as the corresponding edge modes flow in one direction only. Our results can be generalized to other topological superconducting phases, provided that the gap does not vanish on the Fermi surface, in particular to two-component time-reversal-invariant nematic phases Venderbos et al. (2016); Chirolli et al. (2016). Related artificial topological superconductors can be created using the proximity effect Fu and Kane (2008); Zhang et al. (2013); Bernardo et al. (2017).

We describe the system using a two-orbital, $\mathbf{k}\cdot\mathbf{p}$ Hamiltonian proposed for the topological insulator Bi2Se3  Zhang et al. (2009a); Liu et al. (2010); Hsieh and Fu (2012a). We consider quasi-two-dimensional systems, of thickness $d$ such that $k_{Fz}^{-1}\ll d\lesssim\xi_{z}$, where $k_{Fz}$ and $\xi_{z}$ are the Fermi wave-vector and the penetration length along the $z$-axis. The extraction of an effective low-energy two-dimensional Hamiltonian is described in the Appendix. We obtain

[TABLE]

where we denote by $\sigma_{i}$ and $s_{i}$ Pauli matrices that act on orbital space $\{A,B\}$ and spin space $\{\uparrow,\downarrow\}$, respectively, with $i=x,y,z$. We can induce superconductivity in this system via chemical doping. In the mean-field limit, we represent this by a Bogoliubov-de Gennes (BdG) Hamiltonian

[TABLE]

where the $\tau_{i}$ are Pauli matrices that act on the Nambu particle-hole space, $\epsilon_{F}$ is the chemical potential, and the coupling $\Delta_{\text{sc}}\sigma_{y}s_{z}\tau_{x}$ describes a triplet $p$-wave, time-reversal-invariant superconductor.

For Bi2Se3, and other systems of interest, the chemical potential lies far away from the valence band and the magnitude of the superconducting gap is small, i.e. $\bar{\epsilon}_{F}=\epsilon_{F}-E_{G}\lesssim 2E_{G}$ and $\Delta_{\text{sc}}\ll E_{G}$. In this limit, we can neglect the valence band, and project the BdG Hamiltonian onto the conduction band only. To lowest order in $\Delta_{\text{sc}}$, the Hamiltonian reduces to two $2\times 2$ Hamiltonians, $\mathcal{H}_{\pm}$. Thus, we write the projected Hamiltonian as

[TABLE]

where $\mathcal{H}_{\pm}(\mathbf{k})=\pm\left[\frac{v_{F}^{2}|\mathbf{k}|^{2}}{2E_{G}}-\bar{\epsilon}_{F}\right]\tau_{z}+\frac{2v_{F}\Delta_{\text{sc}}}{E_{G}}\left[k_{y}\tau_{x}+k_{x}\tau_{y}\right]$ represents two chiral superconductors with direction-dependent superconducting gaps at the Fermi surface, $\Delta_{\text{sc}}^{\text{Fermi}}=2v_{F}k_{F}\Delta_{\text{sc}}e^{i\varphi}/E_{G},$ where $\varphi=\arctan(k_{y}/k_{x})$ and $k_{F}=\sqrt{2E_{G}\bar{\epsilon}_{F}}/v_{F}.$ The measured superconducting gap is $2\Delta_{\text{exp}}=2|\Delta_{\text{sc}}^{\text{Fermi}}|=4\Delta_{\text{sc}}v_{F}k_{F}/E_{G}$. Each chiral superconductor has a single edge mode with a preferred direction  Read and Green (2000); Qi et al. (2009). The two $p$-wave branches are related by time-reversal symmetry, making the total Hamiltonian time-reversal invariant. Therefore, when solving for the eigenstates of $\mathcal{H}_{0}$ we restrict our attention to finding the states of $\mathcal{H}_{+}$, and determine the eigenstates of $\mathcal{H}_{-}$ by reversing time.

We now study the edge modes on a disc of radius $R$ by diagonalizing $\mathcal{H}_{0}$ in that domain. In the following, we use dimensionless variables by expressing all energies in units of $\bar{\epsilon}_{F}$. Then, we can write

[TABLE]

where $\boldsymbol{\kappa}=\mathbf{k}/k_{F},$ and $\gamma=\Delta_{\text{exp}}^{2}/2\bar{\epsilon}_{F}^{2}\ll 1$ sets the relative scale of the superconducting gap. In our notation, the electron and hole components are ordered as |\Psi\rangle=\left(\begin{array}[]{cccc}\psi_{\uparrow}^{(e)}&\psi_{\downarrow}^{(h)}&-\psi_{\uparrow}^{(h)}&\psi_{\downarrow}^{(e)}\end{array}\right)^{T}. The eigenfunctions of the Hamiltonian in Eq. (4) on a disc are two-component spinors, where each component contains the product of radial Bessel functions with wave-vector $\kappa=|\boldsymbol{\kappa}|$ and angular momentum $\ell$, as discussed in detail in the appendix. The upper and lower components differ by one unit of angular momentum  Stone and Roy (2004). We are interested in the case where $\kappa$ is complex, corresponding to edge modes with an exponentially decaying wavefunction. This occurs for energies, $E$, below the superconducting gap, i.e. $E/\bar{\epsilon}_{F}=\mathcal{E}<\delta=\Delta_{\text{sc}}/\bar{\epsilon}_{F}.$ By imposing the boundary condition that the wavefunction vanishes on the edge, we can find the energy spectrum for the edge states. In the semiclassical limit where the system size is much larger than the Fermi wavelength $\lambda=k_{F}R\gg 1,$ and to first order in $\lambda^{-1}$, the spectra of edge modes for the $p^{+}$ and $p^{-}$ superconductors are

[TABLE]

We note that Eq. (5) is valid up to $|\ell|\approx\lambda\sqrt{1-\gamma/2}$; for larger angular momenta, the states lie outside the superconducting gap, and they are delocalized. Due to circular symmetry, the probability density of the edge states is uniform along the perimeter of the disc. The decay length towards the interior of the disc is $r_{\ast}\propto\xi=\bar{\epsilon}_{F}/(\Delta_{\text{exp}}k_{F})$. This length should be smaller than the disc’s radius. Table 1 list estimates of the coherence length and of other parameters. Scalar disorder leads to pair breaking in $p$-wave superconductors, although the spin-orbit coupling reduces this effect Michaeli and Fu (2012). The parameters in Table 1 use the critical temperatures taken from experiments Hor et al. (2010), which already include the effect of disorder. Finally, it is important to note that scalar disorder does not mix the mid-gap states.

We now consider a symmetry-breaking perturbation that induces zero-mode Majorana states. We apply an in-plane magnetic field that breaks time-reversal symmetry and induces a superconducting gap that vanishes at boundary points where the field is tangent to the disc. We expect to find localized zero-modes at these points. The perturbation is of the form $\mathcal{H}_{Z}=\tau_{Z}\left(\mathcal{E}_{Z}\mathbf{n}\cdot\mathbf{s}\right),$ where $\mathbf{n}$ points in the direction of the field, $\mathcal{E}_{Z}=\mu_{e}B_{Z}/\bar{\epsilon}_{F}$ is the Zeeman coupling scaled by $\bar{\epsilon}_{F},$ and $\mu_{e}$ is the magnetic moment of the electron, see Fig. 1. The Zeeman field couples electrons of opposite spins from the different branches, and likewise, holes from different branches. Note that a parallel magnetic field weakly perturbs the superconducting phase Clogston (1962); Chandrasekhar (1962). We work in the weak-field limit, where the Zeeman energy is much smaller than the superconducting gap. Then, we can neglect the bulk states, and truncate the Hilbert space to the (unperturbed) edge states only. We project the full Hamiltonian $\mathcal{H}=\mathcal{H}_{0}+\mathcal{H}_{Z}$ onto this basis and diagonalize the truncated Hamiltonian. We find two Majorana zero-modes separated in energy from the first excited quasi-particle states by an energy gap $\Delta_{M}$, which can be expressed as $\Delta_{M}/\bar{\epsilon}_{F}\approx\left(8\gamma/\lambda^{2}\times\mathcal{E}_{Z}^{2}\right)^{1/4}$ for large $\lambda$. This implies stability of the zero-mode states as they are well-separated from the other edge modes. A second important scale is the splitting between the two Majorana states. Since the wavefunction at either side of the disc is approximately Gaussian, with angular width is $\langle\theta^{2}\rangle\approx\sqrt{2\gamma}/({\cal E}_{Z}\lambda)$, this leads to a splitting which decays exponentially with $R$. As illustrated here, we can manipulate the gaps that determine the stability of the Majorana modes by simply adjusting the strength of the applied field and the size of the disc.

As the edge modes are localized at boundary points where the superconducting gap changes sign, we can rotate these modes by simply rotating the Zeeman field, as shown Fig. 3. This gives us a simple way to perform particle exchanges for braiding and other purposes. We now explore this possibility in a three-disc configuration. In the presence of tunneling junctions between discs, the Hamiltonian restricted to the Majorana states is

[TABLE]

where $\gamma_{i}$ and $\gamma_{j}$ are Majorana operators in neighboring discs, and $t_{i,j}(\phi)$ is the hopping between those Majorana states. This term depends on the orientation of the field, $\phi$. By rotating the magnetic field, states in different discs can be brought into contact, as illustrated in Fig. 4. For instance, the change of the exchange field which takes the left configuration into the center one exchanges states 3 and 1, while the change from the center configuration to the right one exchanges states 2 and 6. This scheme shows a simple way to manipulate the Majorana fermions. A complete braiding protocol is outside the scope of this work. An alternative proposal, based on a one-dimensonal ring of magnetic atoms on a superconductor has been recently discussed in Li et al. (2016) (note that the braiding scheme discussed there can be extended to our proposal, using overlapping dots). For other realizations, see also Nadj-Perge et al. (2013, 2014).

The appeal of this approach is that no tri-junctions are required, and only a rotation of the magnetic field is needed. These operations are carried out without the recourse to external gates or magnetic fluxes. If such fluxes and gates are added to this system, more ways to correlate the Majorana particles are induced, leading to new functionalities. Furthermore, more complex dot geometries can also be a platform for more ambitious engineering efforts. We explore some of these possibilities in the Appendix.

It may also be the case that the magnetic field is due to spontaneously polarized magnetic moments or to an additional chiral superconducting component of the order parameter Chirolli et al. (2016). Then, quantum fluctuations of the field will lead to additional interactions between the Majorana states. Finally, a large array of quantum dots can serve as a platform for a surface code for topological quantum computation Bravyi and Kitaev (1988); Vijay et al. (2015); Vijay and Fu (2016).

One question that needs to be addressed is whether we can perform the magnetic field rotation without substantially mixing in low-lying fermionic excitations (so-called quasiparticle poisoning Sarma et al. (2015)). As shown in detail in the Appendix, in the limit of weak coupling, we can adjust the separation in energy scales to be large enough to assure that this does not happen. There is a subtle interplay between effective tunneling-supported mixing of the Majorana particles and quasiparticle poisoning that deserves further investigation.

In summary, we have analyzed the emergence of localized Majorana modes at the edges of finite two-dimensional topological superconductors. Though we have based our model on the time-reversal-invariant triplet superconducting phase that is likely to exist in doped Bi2Se${}_{3},$ our results apply more generally to systems where the thickness is much smaller than the bulk coherence length. In particular, the triplet superconducting phase of doped Bi2Se3 bears qualitative resemblance to the B phase of superfluid 3He, and therefore, we expect our results to hold in that context as well. Furthermore, the analysis remains valid for other gapped superconducting phases with gapless edge modes, such as the nematic phase, also proposed for doped Bi2Se3. Similar Majorana states can also be expected in fermionic superfluids based on cold atoms.

Our system provides an appealing alternative Majorana platform precisely because of its simplicity. We can control the stability of the Majorana states by widening the gap to the next excited states with a simple adjustment of the field strength. The positions of the Majorana states can be modified by changing the orientation of the applied magnetic field. Already with a single disc, we can exchange two Majorana particles by rotating the field by $\pi.$ In arrays with many quantum dots, the field can be used to modulate the interaction between Majorana states in different dots, and to exchange them, without requiring the existence of tri-junctions, electrostatic gates, or magnetic fluxes. These properties make the proposed model an interesting candidate for the experimental realization of a Majorana-based qubit for topological quantum computation.

Appendix A Derivation of the Projected Hamiltonian

In this section, we derive the projected Hamiltonian describing the two-dimensional $p$-wave topological superconductor used in the main text. A simple four-band model for the valence and conduction bands of the topological insulator Bi2Se3 near $\Gamma$ is Zhang et al. (2009a); Liu et al. (2010); Hsieh and Fu (2012b)

[TABLE]

where parameters $E_{G},$ $v_{z},$ and $v$ are determined from experimental data. The four bands correspond to the two dominant Wannier orbitals $\{A,B\}$, each with spin $\{\uparrow,\downarrow\}$. The $2\times 2$ Pauli matrices $\sigma_{i}$ and $s_{i}$, with $i=x,y,z,$ act on orbital and spin spaces respectively. We use units where $\hbar=1$, so that all dimensions can be expressed in terms of energies and lengths.

The bands are doubly degenerate with dispersion

[TABLE]

The gap between the conduction band and the valence band is thus $2E_{G}$. In order to describe a superconductor, we need to extend the Hilbert space to include hole wavefunctions and a corresponding hole Hamiltonian. Using the convention given in Ref. Hsieh and Fu (2012b), the full wavefunction now contains eight components,

[TABLE]

In this basis, the $8\times 8$ Hamiltonian is simply

[TABLE]

where $\epsilon_{F}$ is the Fermi energy, and $\boldsymbol{\tau}$ is a third set of Pauli matrices that couples particles and holes of the same spin and in the same orbital. The particles and holes are now mixed by a triplet, time-reversal-invariant coupling of the form $\Delta=\Delta_{\text{sc}}\sigma_{y}s_{z}\tau_{x}$ that induces a superconducting gap proportional to $\Delta_{\text{sc}}$.

If the energy gap is large, i.e., $2E_{G}\gtrsim\epsilon_{F}-E_{G}$ and $E_{G}\gg\Delta_{\text{sc}}$, only the lower part of the conduction band is occupied and we can neglect the valence band, as the conduction and valence bands are separated by a large energy gap $2E_{G}$. To simplify the result, we first perform a rotation in orbital space by $\pi/2$ about the $y$ axis, so that $\sigma_{y}\rightarrow\sigma_{y},\sigma_{x}\rightarrow\sigma_{z},\sigma_{z}\rightarrow-\sigma_{x}$. In this new representation, the states at the edge of the conduction and valence bands are eigenstates of $\sigma_{z}$, with eigenvalues $\pm 1$. The operators ${\cal P}_{\pm}=({\cal I}\pm\sigma_{z})/2$, where ${\cal I}$ is the identity matrix, thus project the wavefunctions onto the conduction and valence bands respectively. The kinetic energy in Eq. (A1), which is now proportional to $\sigma_{x}$, and the gap, which remains proportional to $\sigma_{y}$, mix the conduction and valence bands. Using perturbation theory, we can expand the Hamiltonian to first order in $\Delta_{\text{sc}}$ and then project onto the conduction band to find the reduced Hamiltonian

[TABLE]

Using $\sigma_{z}^{-1}=\sigma_{z},\tau_{z}^{-1}=\tau_{z}$, and adding the chemical potential that was omitted in Eq. (A6), we obtain

[TABLE]

where $\bar{\epsilon}_{F}=\epsilon_{F}-E_{G}$.

This Hamiltonian defines a three-dimensional topological superconductor. The Fermi wave-vector and the coherence length along to the $z$-direction are given by

[TABLE]

We now analyze a quasi-two-dimensional system of thickness $d$, where $k_{Fz}^{-1}\ll d\lesssim\xi_{z}$. In the absence of superconductivity, the Hamiltonian can be separated as ${\cal H}^{\text{red}}={\cal H}_{\parallel}+{\cal H}_{z}$, where ${\cal H}_{\parallel}$ and ${\cal H}_{z}$ are the in-plane and out-of-plane components of the Hamiltonian, respectively. The electronic states can then be divided into sub-bands, $\epsilon_{n}({\bf k})=\epsilon_{z}^{n}+\epsilon_{\parallel}({\bf k})$. Here, $\epsilon_{z}^{n}$ and $\epsilon_{\parallel}({\bf k})$ are eigenvalues of ${\cal H}_{z}$ and ${\cal H}_{\parallel}$, respectively, and $\mathbf{k}$ is the two-dimensional wave-vector in the $x$-$y$ plane. The separation of levels due to the quantization in the $z$-direction is, approximately,

[TABLE]

where $\Delta_{\text{exp}}$ is the value of the measured superconducting gap at the Fermi energy, see Eq. (A7). In order to define a low-energy effective theory for the electronic structure, we study first the lowest energy state for $|{\bf k}|\approx 0$ and energy $|\epsilon|\sim\Delta_{\text{sc}}\ll\delta\epsilon_{z}^{n}$.

For $|{\bf k}|=0$, the Hamiltonian in Eq. (A7) can be split into two $2\times 2$ Hamiltonians, one which mixes electrons and holes with spin up, and another one which mixes electrons and holes with spin down. The two Hamiltonians are related by time-reversal symmetry. For $d\gtrsim\xi_{z}$ the lowest states of these Hamiltonians are localized at the top and bottom surfaces. The decay length of these states is given by the coherence length $\xi_{z}$. When $d\lesssim\xi_{z}$, these states hybridize, and give rise to two Andreev states at energies $\pm\epsilon_{0}\lesssim\Delta_{\text{sc}}$. Thus, the combination of the two $2\times 2$ Hamiltonians gives rise to four low-energy states, $|i\rangle,i=1,\cdots,4$, related by particle-hole and time-reversal symmetries.

A generic wavefunction in this low-energy sector can be written as $\sum_{i}\psi_{i}(\mathbf{r})|i\rangle$, where $\mathbf{r}$ is a two dimensional vector. Using the basis of states $|i\rangle$, we can write an effective Hamiltonian whose states are four-component spinors. Without loss of generality, we can label the two states at energy $+\epsilon_{0}$ as $|e\uparrow\rangle$ and $|e\downarrow\rangle$, and the two spinors with energy $-\epsilon_{0}$ as $|h\uparrow\rangle$ and $|h\downarrow\rangle$. We assume that the matrix elements of the spin Pauli matrices $s_{i}$ and the electron-hole Pauli matrices $\tau_{i}$ ($i=x,y,z$) have the same values in this basis as in the original basis of spins and electron-hole states. Hence, we can write the effective Hamiltonian for the in-plane wavefunction in terms of $x$ and $y$ coordinates as

[TABLE]

This Hamiltonian can be split into two independent Hamiltonians related by time-reversal symmetry. These Hamiltonians mix electrons with spin up and holes with spin down, and vice versa.

We make an approximate decomposition into a state with small $|\mathbf{k}|$, and a lowest energy $k_{z}$ state. To good approximation, the Hamiltonian depends only on $k_{z}$,

[TABLE]

A standard plane-wave solution to this Hamiltonian has eigenvalues

[TABLE]

This spectrum has an energy gap of about $4/\xi_{z}k_{Fz}$ due to superconductivity.

We can solve for the eigenenergies; values of the parameter $\epsilon_{0}$ in Eq. (A10) are shown in Fig. A1 for a realistic choice of parameters. Here we require that the superconducting gap is 0.6 meV. For energy eigenvalues inside the gap, $\epsilon_{0}/\bar{\epsilon}_{F}<2/\xi_{z}k_{Fz},$ we have states that are weakly localized at the top and bottom surfaces, with an exponential decay slower than the thickness of the disc. As illustrated in Fig. A1, there is a wide range of values for $d$ for which we find a state inside or near the gap.

We note that for $k_{Fz}^{-1}=20\leavevmode\nobreak\ \mathrm{\AA}$, we have $\xi_{z}=4000\leavevmode\nobreak\ \mathrm{\AA}$ (a smaller Fermi wavevector leads to an even larger coherence length for a fixed gap). This shows that our system needs not be atomistically thin, but can have a thickness on the order of a thousand angstroms, similar to the coherence length. Henceforth, without loss of generality, we assume we have chosen a particular thickness $d$ for the disc where no localization at the surfaces occurs. We also absorb $\epsilon_{0}$ into the definition of $\bar{\epsilon}_{F}$. This does not change the value of $\bar{\epsilon}_{F}$ appreciably because $\bar{\epsilon}_{F}$ is much larger than $\epsilon_{0}.$

The two reduced Hamiltonians can be interpreted as describing two metals with quadratic dispersion and superconducting gaps at the Fermi surface with angular dependence $\Delta_{\text{sc}}^{\text{Fermi}}=2v_{F}k_{F}\Delta_{\text{sc}}e^{i\varphi}/E_{G}$, where $\varphi=\arctan(k_{y}/k_{x})$ is the angular coordinate in momentum space, and $k_{F}=\sqrt{2E_{G}\bar{\epsilon}_{F}/v_{F}^{2}}$. In experiments, the more relevant energy scale is the superconducting energy gap measured at the Fermi surface $2\Delta_{\text{exp}}=2|\Delta_{\text{sc}}^{\text{Fermi}}|=4v_{F}k_{F}\Delta_{\text{sc}}/E_{G}.$ The two Hamiltonians are

[TABLE]

The two Hamiltonians are related by time inversion. Each Hamiltonian describes a chiral $p$-wave superconductor Read and Green (2000). The energy dispersion is

[TABLE]

Using the convention where

[TABLE]

the full Hamiltonian should be written as

[TABLE]

where

[TABLE]

and ${\cal H}_{-}(\mathbf{k})=-{\cal H}_{+}(-\mathbf{k})$.

Appendix B Edge States on a Disc

In this section, we solve for the edge modes on a disc. As noted above, the low-energy states in the conduction band of the two-dimensional superconductor doped Bi2Se3 can be described by two effective $p$-wave Hamiltonians that couple electrons and holes of opposite spins defined by Eqs. (A18) and (A21). The $p^{+}$ and $p^{-}$ superconducting branches are related by time-reversal symmetry. Thus, since eigenstates of $\mathcal{H}_{+}$ are also eigenstates of $\mathcal{H}_{-}$, we restrict our attention to $\mathcal{H}_{+}$ on a circle of radius $R$. In polar coordinates, we have $\pm ik_{x}+k_{y}=\pm e^{\pm i\theta}\left(\partial_{r}\mp i\partial_{\theta}/r\right)$, and $k_{x}^{2}+k_{y}^{2}=-\nabla^{2}=-\partial_{r}^{2}-\partial_{r}/r-\partial^{2}_{\theta}/r^{2}$. Using separation of variables, we can write our wavefunction in coordinate space as

[TABLE]

To insist that the wavefunction is single-valued, we have that $\ell$ must be an integer. In this case, $\ell$ takes the role of describing our angular momentum quantum number. Let us now show that indeed, $\ell$ can be used as a good quantum number. The angular momentum operator in this case is

[TABLE]

where $\mathcal{I}$ is the identity matrix. The spectrum of $L_{z}$ is $L_{z}\Psi_{\ell}=\left(\ell+\frac{1}{2}\right)\Psi_{\ell}$. By using the fact that partial derivatives commute, it immediately follows that $[\mathcal{H}_{+},L_{z}]=0$. Consequently, we can diagonalize both $\mathcal{H}_{+}$ and $L_{z}$ simultaneously, and use $\ell$ as a good quantum number.

Using the form of $\Psi$ in Eq. (A22), we find that the radial wavefunctions satisfy the following eigenvalue problem

[TABLE]

To solve for the radial wavefunction, we make an ansatz that the wavefunction components are Bessel functions based on the motivation that eigenfunctions of the Laplacian are Bessel functions

[TABLE]

where $a$ and $b$ are constants, $J_{\ell}(x)$ is the Bessel function of the first kind of order $\ell,$ and $k$ is the linear momentum in the radial direction. Substituting this ansatz into Eq. (A24), our problem is reduced to a simple eigenvalue problem

[TABLE]

This has unnormalized solutions

[TABLE]

where

[TABLE]

Having found the eigenstates of our topological superconductor, we can now study the effect of the edge. We impose hard-wall boundary condition whereby we require the wavefunction to vanish at the perimeter of the disc $\Psi_{\ell}(R,\theta)=0$. This condition is equivalent to continuity of the Dirac wavefunction at the boundary where we have Andreev reflections Stone and Roy (2004). In order to satisfy this condition, the wavefunction for a particular energy $E$ must be a linear combination of outgoing and reflected waves.

Let us now simplify by using dimensionless variables, where we express lengths in terms of the Fermi momentum $k_{F}=\sqrt{2E_{G}\bar{\epsilon}_{F}}/v_{F}$ as defined in the main text, and all energies in terms of $\bar{\epsilon}_{F}$,

[TABLE]

where $\gamma=\frac{4\Delta_{\text{sc}}^{2}}{E_{G}\bar{\epsilon}_{F}}=\frac{\Delta_{\text{exp}}^{2}}{2\bar{\epsilon}_{F}^{2}}\ll 1$ measures the superconducting gap relative to $E_{G}$ and $\bar{\epsilon}_{F}$ and $\mathcal{E}=E/\bar{\epsilon}_{F}$. In terms of these parameters, the full wavefunction is

[TABLE]

where $\rho=r/R$ is the fractional radius, and $\lambda=k_{F}R$ is the dimensionless system size. Setting both components of the above equation to zero at $\rho=1$, we obtain the following transcendental equation for the allowed energies

[TABLE]

We can numerically solve this transcendental equation for some typical parameter values, given that we stay within the regime of validity for our effective Hamiltonian. To obtain an edge state with energy below the gap, we can consider the limit of a very large disc, as done in Stone and Roy (2004). Suppose that $\lambda\gg 1$, then for a very small angle $\theta\ll 1$, the boundary resembles that of the infinite half plane for which we know that the edge state has energy $\mathcal{E}^{\text{edge}}=-\frac{2v_{F}k_{x}\Delta_{\text{sc}}}{E_{G}\bar{\epsilon}_{F}}$, where $k_{x}$ is the linear momentum along the edge Stone and Roy (2004); Fendley et al. (2007). Translating this result to our geometry in the semi-classical limit, we have $|\mathbf{L}|=|\mathbf{r}\times\mathbf{p}|=\ell\approx Rk_{x}.$ Thus, approximately, the energy of the edge state in the semi-classical limit is $\mathcal{E}_{\text{edge}}\approx-\frac{\ell}{\ell_{0}}\delta$, where $\delta=\Delta_{\text{sc}}/\bar{\epsilon}_{F}$ and $\ell_{0}=RE_{G}/2v_{F}$. We can tighten this estimate further using a general argument of particle-hole symmetry; that is, if $\mathcal{H}\Psi=E\Psi$, then $\mathcal{H}(\tau_{x}\mathcal{C}\Psi)=-E(\tau_{x}\mathcal{C}\Psi),$ where $\mathcal{C}$ is the complex-conjugation operator. With this, we see that a state of angular momentum $\ell$ has the opposite energy to a state of angular momentum $-(\ell+1)$. Therefore, we find that the edge energy is

[TABLE]

This can be further improved by a more detailed asymptotic analysis. We find

[TABLE]

with a numerical indication that the next order has a more complex dependence on $\ell$.

We compare in Fig. A2 the functional form of Eq. (A32) to the lowest energy found from the transcendental equation in Eq. (A31). To leading order in $\lambda^{-1}$, there is good agreement between the numerical result and the approximation in Eq. (A32). This energy spectrum lies inside the superconducting gap, and corresponds to edge states at the perimeter of the disc. One such edge state for $\ell=0$ is shown in Fig. A3. The edge mode we have found is the chiral edge state we have sought from the beginning. From Eq. (A32) we see that each of the Hamiltonians $\mathcal{H}_{\pm}$ describes a single chiral band with unidirectional group velocity—negative for the case described here. This reflects the fact that a single chiral Hamiltonian breaks time-reversal invariance.

To recover time-reversal invariance, we now consider the full Hamiltonian, Eq. (A18), of the two-dimensional topological superconductor. The spectrum of edge modes for $\mathcal{H}_{-}$ consists of the time-reversed states of $\mathcal{H}_{+}$. Suppose a state $\Psi_{+}^{\text{edge}}$ is an eigenstate of $\mathcal{H}_{+}$ of angular momentum $\ell$, then the time-reversed state with angular momentum $-(\ell+1)$ is an eigenstate of $\mathcal{H}_{-}$. Therefore, the edge spectrum of $\mathcal{H}_{-}$ is

[TABLE]

This is just the chiral spectrum for modes propagating in the other direction. Taken together, Eqs. A32 and A34 give the full spectrum of the two edge modes for our time-reversal-invariant topological superconductor. The band structure of this system is shown in Fig. A2.

Appendix C Majorana States in the Presence of a Zeeman Field

In this section, we consider the effect of a constant in-plane Zeeman field defined by

[TABLE]

where $\mathcal{E}_{Z}$ is the energy magnitude of the field, $\mathbf{n}$ is a unit normal along the direction of the field, and $\phi$ is the angle of that field. This Zeeman field couples electrons of opposite spins and holes of opposite spins, both with the same sign. By construction, this addition to the Hamiltonian explicitly breaks time-reversal symmetry. The full Hamiltonian now becomes

[TABLE]

where $\mathcal{H}_{0}$ is defined by Eq. (A18). We consider $\mathcal{H}_{Z}$ as a weak perturbation to $\mathcal{H}_{0}$, in comparison to the low-lying edge states with energies below the superconducting gap. That is, we assume $\mathcal{E}_{Z}\ll\delta$. So we can use lowest-order perturbation theory to assess how the energy spectrum of the low-lying edge states is changed. We first write explicitly the full edge wavefunctions. With $C_{+}$ and $C_{-}$ in Eq. (A30) related by the vanishing boundary condition

[TABLE]

we can write the wavefunction in the $\mathcal{H}_{+}$ subspace as

[TABLE]

where $\mathcal{N}$ is a normalization constant. We can simplify the wavefunction further by noting that

[TABLE]

where the last equality follows from inspection of Eq. (A29) that for $\mathcal{E}<\delta$, $\kappa_{\pm}$ is complex, and consequently, $\kappa_{+}^{*}=\kappa_{-}$. Using this, it is clear that the radial part of Eq. (A38) is purely imaginary

[TABLE]

The eigenstates of the $\mathcal{H}_{-}$ subspace are simply the time-reversed versions of the $\mathcal{H}_{+}$ states. We obtain the $\mathcal{H}_{-}$ eigenstates by taking the complex conjugate of Eq. (A40), mapping $\ell\mapsto-(\ell+1)$, inverting the order of the components, and adding an appropriate minus sign

[TABLE]

In this form, we have electron and hole components in the order defined in Eq. (A17). We can now study the coupled Hamiltonian by diagonalizing $\mathcal{H}_{Z}$ in this basis. There is a natural cut-off on the value of $\ell$ using the fact that the square-root in Eq. (A29) must be real

[TABLE]

Here, we see that $0<\gamma<2.$ However, as we will show in a later section, it is, in fact, the case that for a range of semi-realistic parameters, $\gamma\ll 1.$ The diagonal elements of the matrix then contain the unperturbed energies, with the upper-left block containing the $p^{+}$ spectrum, and the lower-right block containing the $p^{-}$ spectrum. The Zeeman interaction only mixes electrons and holes of opposite spins, and thus couples the $p^{+}$ and $p^{-}$ subspaces. The matrix elements are of the form

[TABLE]

where $f_{\text{rad}}$ and $g_{\text{rad}}$ are radial integrals. Therefore, we see that this perturbation only couples states of angular momenta $\ell$ and $\ell\pm 1$. This is shown schematically in Fig. A4. The Hamiltonian can now be diagonalized numerically to find the spectrum of edge modes. Since angular momentum is not conserved because the perturbation is not axis-symmetric, we label states of the new spectrum by $n$ instead of $\ell$. One can think of $n$ as the average angular momentum after hybridization from the perturbation. The result of a particular realization of this perturbed model is shown in Fig. A5. In this case, we use $\lambda=50$ and $\mathcal{E}_{Z}=1/20$. As shown in Fig. A5, adding this perturbation generates a gap, but there remain two localized modes at almost zero energy. This is due to particle-hole symmetry that is preserved by the perturbation. In other words, for every positive-energy state, there exists a negative-energy state of the same amplitude in the unperturbed spectrum. The addition of the Zeeman perturbation shifts the energies, but preserves the fact that there is a hole state for every particle state. Hence, the states with approximately zero energy come in pairs. Furthermore, we find numerically that the radial integrals in Eq. (A43) are almost exactly 1/2, with small deviations only when $\ell$ gets close to $\pm\ell_{\text{max}}$.

We now estimate the magnitude of resulting gap due to the perturbation. As we have established, the spectra of the unperturbed edge states can be well-approximated by a linear energy dispersion of $\ell+1/2$ with opposite chirality. Therefore, we can write the Hamiltonian corresponding to this two-state system as

[TABLE]

The Zeeman field couples states in one branch with angular momentum $\ell$ with states in the other branch with angular momentum $\ell\pm 1$. We can correspondingly decompose the wavefunctions as

[TABLE]

where the quantities $a_{\ell,\pm}$ are amplitudes. Next, we use ${\cal H}_{Z}|\ell,+\rangle=\frac{\mathcal{E}_{z}}{2}(|\ell+1,-\rangle+|\ell-1,-\rangle)$. We can make a gauge transformation which changes the sign of one half of the amplitudes $a_{\ell,\pm}$, which is nothing more than a $\pi/2$ rotation of the axes. Thus, if these new amplitudes are slowly varying, we find that for the matrix representation of ${\cal H}_{Z}$, we can write ${\cal H}_{Z}a_{\ell,+}\approx\mathcal{E}_{z}\frac{\partial a_{\ell,-}}{\partial\ell}$. The Hamiltonian can now be approximately expressed as

[TABLE]

where $\sigma_{x}$ and $\sigma_{z}$ are Pauli matrices, and we define $\hat{\ell}$ and $\hat{p}_{\ell}$ as the position and momentum operators. By squaring the Hamiltonian, we find

[TABLE]

This squared Hamiltonian describes a harmonic oscillator. The eigenenergies are, using $m$ as the oscillator quantum number,

[TABLE]

This shows that the zero-energy mode occurs once, and every other energy occurs twice. We can thus express the final spectrum of the edge Majorana modes in terms of $n=\pm m$ as

[TABLE]

where $n=...,-2,-1,0,1,2,.$.. The wavefunction of the zero-energy Majorana modes is a Gaussian in angular momentum space with harmonic oscillator length parameter

[TABLE]

This parameter is $\gtrsim 1$, which agrees with the fact that we have assumed a slow variation of the coefficients $a_{\ell,\pm}$. The implication is that the Gaussian width in the dual, and more physically relevant, variable $\theta$ is given by $1/b_{\ell}$ and is reasonably narrow. This implies that the wavefunction of the Majorana zero-modes is narrowly localized in angular space. We can see that this width decreases as we increase the radius of the disc or the coupling to the magnetic field.

The expression (A49) is valid when the width in $\theta$ is so small that the Majorana wavefunctions do not significantly overlap. As can be seen in Fig. A5, the harmonic spectrum fits well for realistic parameters, but works best for small $n$, since the approximation of a very narrow state is not completely satisfied. We can improve the agreement by increasing the value of $\lambda$. From these results, we can estimate that the gap between the Majorana zero-modes and the nearest quasiparticle states is given by

[TABLE]

which is larger than the spacing between the first quasiparticle states by a factor of $1/(\sqrt{2}-1)\approx 2.4$.

The derivation above assumes three simplifications, which are all irrelevant for the estimate of the quasi-particle gap $\Delta_{M}$, but have consequences for the energies of the mid-gap states that are slightly displaced from zero. We thus have a gap in the spectrum of the quasi-zero-modes, which we shall denote by $\delta\epsilon_{M}$. This effect that the Majorana states are actually only quasi-Majorana states is common to all approaches to topological computing.

The three approximations are, in order of the severity of their impact on $\delta\epsilon_{M}$:

The approximation that there is no upper limit for $\ell$, whereas we know that $\ell_{\text{max}}$ is finite (see Eq. (A42)), 2. 2.

The assumption that the radial matrix elements in Eq. (A43) are all $1/2$, and 3. 3.

The neglect of the periodicity of $\theta$ in Eq. (A46).

We now investigate all of these here. The first point can be dealt with exactly by solving the harmonic oscillator problem with a cut-off on $\ell$ by placing our harmonic oscillator into an infinite well in $\ell$ space with hard walls located at $\pm\ell_{\text{max}}.$ This will only make a significant contribution when the wavefunction in $\ell$ is wide compared to $\ell_{\text{max}}$, and is proportional to $\exp(-b_{\ell}^{2}/\ell_{\text{max}}^{2})$, which is extremely small for a very wide range of parameters. Therefore, this contribution to the splitting of the mid-gap energies is not relevant for the system considered here. Since there is so little sensitivity to the behavior at large values of $\ell$, the second point is also not important; we find that the dominant source of splitting comes from the periodicity in $\theta$.

It is straightforward to show that, after making the same change-of-sign transformation on the coefficient $a_{\ell}$ as before, but not ignoring the fact that $\theta$ is periodic, we find that the Hamiltonian (A46) including periodicity is

[TABLE]

Once again we can square this Hamiltonian, but the result is not as illuminating as before,

[TABLE]

which is a set of two Mathieu-like equations with shifted derivative with a coupling term.

We have not attempted to solve this problem in a general way, even though there are many indications that this may be possible. Instead, we use a combination of numerical calculations with some analytical insights to find an appropriate dependence of the splitting of the mid-gap energies on this periodicity in $\theta.$ We start from Eq. (A52), and perform a scaling inspired by the harmonic oscillator results (for $\Delta_{M}$ see Eq. (A51))

[TABLE]

Thus we know that all eigenenergies can be exactly be expressed in the form $\Delta_{M}f_{n}(b_{\ell}),$ where $f_{n}$ is some function that depends on the energy level $n$ and the parameter $b_{\ell}.$ Indeed, numerical diagonalizations with a finite cutoff $\ell_{\text{max}}$ give exactly this scaling behavior. Actually, in the case of small splitting, we find that a good approximation for $f_{0}$ is given by

[TABLE]

with no deviations until $f_{0}$ gets as large as $0.1$–well outside the regime of interest. More typical values of $b_{\ell}=2-3$ give $f_{0}=2\times 10^{-4}-10^{-8}$. This establishes that the splitting of the mid-gap modes is very small.

Appendix D Estimates of Physical Parameters

In this section, we give estimates of the physical parameters to hopefully guide future experimental attempts to realize our proposal. We note that although our estimates will based primarily on data that are available on CuxBi2Se${}_{3},$ our proposal is not limited to just this two-dimensional material, as we have already explained in the main text.

First, we estimate the physical length scale of the disc. The three-dimensional electron density of CuxBi2Se3 has been reported to be about $n\approx 2\times 10^{20}$ cm*-3* Hor et al. (2010). This value varies slightly depending on dopant concentration of Cu Kriener et al. (2011). Using the fact that each quintuple layer of CuxBi2Se3 is about $1$ nm Zhang et al. (2009b), we find that the electron density per layer is $\rho_{L}\approx 1.6\times 10^{13}$ cm${}^{-2}.$ The Fermi wavelength is then given by $k_{F}=\sqrt{2\pi\rho_{L}}=10^{-1}$ Å${}^{-1}.$ We want the Fermi wavelength to be much smaller than the radius of the disc. Therefore, we choose $\lambda$ to be large. For instance, with $\lambda=100,$ the corresponding radius is $R=1$ $\mu$m. For the simulations shown in this paper, we have chosen $\lambda=50$ to improve visibility as the edge modes become too localized to clearly see with increasing radius.

Next, we determine the various energy scales of the proposal. The critical temperature $T_{c}$ is found experimentally to be 3.8 K Hor et al. (2010). The superconducting gap is then given by $\Delta_{\text{exp}}=1.76k_{B}T_{c}=k_{B}\times 6.7\,\mathrm{K}=0.6\,\mathrm{meV}$ . We note that scalar disorder can break $p$-wave superconducting pairs, leading to a lower $T_{c}.$ This means that it is possible in future experiments that cleaner samples can host superconductivity with a higher $T_{c}$ and larger gap. It has been reported that the effective mass $m_{\text{eff}}$ in CuxBi2Se3 is about 0.3 $m_{e}$ electron masses Sandilands et al. (2014). Using this, we can estimate the Fermi energy $\epsilon_{F}=\hbar^{2}(3\pi^{2}n)^{2/3}/2m_{\text{eff}}\approx 350$ meV. Now, if we take $E_{G}=150$ meV Wray et al. (2011), this then gives the shifted Fermi energy $\bar{\epsilon}_{F}=\epsilon_{F}-E_{G}=v_{F}^{2}k_{F}^{2}/2E_{G}=200\,\mathrm{meV}=k_{B}\times 2300\,\mathrm{K}.$ The coherence length in the BCS theory is given by $\xi=\hbar v_{F}/\pi\Delta_{\text{exp}}=\bar{\epsilon}_{F}/\Delta_{\text{exp}}k_{F}\approx 3.4\times 10^{3}$ Å. This calculation justifies our assumption of quasi-two-dimensionality because the coherence length is much larger than the thickness of each quintuple layer of CuxBi2Se3.

Finally, we estimate the physical quantities relevant when a Zeeman field $B_{Z}=1$ T is applied. When a magnetic field is applied, Majorana zero-modes emerge, localized in angular space where the field is tangent to the disc’s surface. The angular width of the Majorana states is given by the inverse of Eq. (A50), $b_{\ell}^{-1}=\sqrt{\Delta_{\text{exp}}/\mu_{e}B_{Z}Rk_{F}}\approx 0.32,$ where $\mu_{e}$ is the Bohr magneton of the electron. The energy gap between the zero-modes and the first quasiparticle excitations is $\Delta_{M}=\sqrt{2\mu_{e}B_{Z}\Delta_{\text{exp}}/Rk_{F}}\approx 0.026$ meV $=k_{B}\times 0.3$ K. The gap between the quasi-zero-energy Majorana modes is $\delta\epsilon_{M}=\Delta_{M}\exp(-2b^{2})/\sqrt{2}\approx 10^{-9}$ meV $=k_{B}\times 10^{-8}$ K.

We note that the above parameters are semi-realistic estimates. Using these, the dimensionless paramters are $\lambda=100,$ $\gamma<10^{-5}$ and $\mathcal{E}_{Z}<10^{-3}.$ However, in numerical simulations, we have chosen these dimensionless parameters to be $\lambda=50,$ $\gamma=1/16,$ and $\mathcal{E}_{Z}=1/20$ for improved visibility.

Appendix E Tight-Binding Calculations

In this section, we present alternative calculations of our model using the tight-binding approach. This can either be thought of as a microscopic realization of the continuum model or a finite-difference discretization. We choose to work on a square lattice. However, the results obtained for system sizes much larger than the lattice spacing should be independent of this choice. The process of nearest-neighbor hopping is most easily described on a bipartite lattice. For a single chiral superconductor, corresponding to a single spin projection, the gap couples the conduction band to the valence band. Each occupied site has an energy $\pm\epsilon_{F}$ (plus for particles, minus for holes), and the hopping constant is $\pm t$ (where the sign indicates the direction of hopping). The superconducting gap $\Delta_{\text{sc}}$ couples particles to holes with the usual Dirac matrix structure, as shown in Fig. A6. The full model then combines the spin-up and spin-down copies of the chiral superconductor with conjugate Hamiltonians. The Zeeman term couples particles with spin up and spin down, and similarly for the holes. In summary, if $\mathbf{r}=(i,j)$ labels the sites, $k=1,2$ labels particles or holes respectively, $l=1,2$ labels spin up or down respectively, $\boldsymbol{\rho}$ links nearest-neighbor points, $a$ is the lattice spacing, $z$ is the coordination number, and the indices on the Pauli matrices denote a matrix element, then the Hamiltonian has the following form on any bipartite lattice:

[TABLE]

This has the continuum limit

[TABLE]

which makes it simple to read off the relation between the parameters of the continuum and tight-binding models. With this tight-binding model, we can study any geometry, and on any bipartite lattice contained within that geometry. For the results reported here, we use a square lattice, but we have checked that we can obtain the same answers for a hexagonal lattice. We first study a disc, or rather, an approximation to a disc built from squares, to compare with our analytic calculations on the same geometry. We find that the numerical results are identical to those of the continuum model, as shown in Fig. A7a. What cannot be seen from that figure is the angular dependence of the resulting wavefunctions; these are in agreement with the continuum model as well. If we add an in-plane magnetic field, we find the Majorana zero-modes as seen before, shown in Fig. A7b-c. We note that even if the field is not aligned along a symmetry axis of the lattice, where the lattice has more effect, our tight-binding model still produces consistent results, as illustrated in Fig. A7c.

Next, due to our choice of a square lattice, the simplest geometry to study is that of a 2D quantum wire, the de-facto standard for topological quantum computing. We look for states with energies below the superconducting gap. In the case of a disc, all of these states are boundary states localized near the edge. However, in a wire geometry, we no longer have this property as some of these mid-gap states behave more like bulk states. As can be seen in Fig. A8i-ii, in the absence of a Zeeman field, we find states in the $p_{x}+ip_{y}$ superconductor that are localized at the ends of the wire in the two lowest-energy states. However, for the next state higher in energy, the mid-gap state looks more like a bulk state, as shown in Fig. A8iii. When we turn on the in-plane magnetic field perpendicular to the long axis of the wire, we get two Majorana states, each located at one end of the wire, as shown in Fig. A8vi-vii. The numerical results violate reflection symmetry, due to an exact degeneracy in the spectrum. As we rotate the field, we see that this degeneracy is broken, and we move from states with the Majorana modes located at the ends to delocalized states as the field aligns with the long axis. This can be seen most clearly in the energy-level plot in Fig. A8, where it is shown how the degeneracy in zero-energy states is lifted as we get close to the case where the field is parallel to the wire direction.

Of more interest to quantum-information processing is entangling the Majorana states in two adjacent discs connected by thin wires by rotating the magnetic field. Again, we first look at the $p_{x}+ip_{y}$ model in this setup in the absence of a Zeeman field, as shown in Fig. A9. We see that inside the discs, there is no distinction between the behavior of electrons and holes, but in the connecting wire, there is an interesting difference. This has not yet been investigated in detail in this work, but will likely have an impact on how to construct analytic models of this situation. We now move onto assessing the model in the presence of a Zeeman field. For a field that is not parallel to the connecting wire, we find the expected Majorana states, as seen in Fig. A10b-c. As we turn on the field parallel to the wire, we find that more states come into play, as shown in Fig. A10d-f. We see very low energy states that are completely located only in the connecting wire, which is clearly a large perturbation on the situation of two unconnected discs.

Since we want to reduce the coupling between the two superconducting discs, we consider a long, thin connecting wire next, with results shown in Fig. A11. We see that in such a situation, we can instead realize a weak perturbation, but there is still non-zero probability of the zero-modes penetrating the wire, and thus the Majorana states still couple.

A most straightforward way to assess the coupling quantitatively is to look at the energy spectra of the two setups of either a thick or a thin connecting wire. Comparing Fig. A10a and Fig. A11a, we see that for a thicker wire, we get a strong mixing between quasiparticles and Majorana states, i.e. a clear example of quasiparticle poisoning. For a much thinner wire, we get almost no coupling between the Majoranas. The tunneling splitting is much smaller than the gap to the first quasiparticle state. The latter case of a thinner wire would thus make a better candidate for connecting the three-dot device discussed in the main text for the purposes of performing quantum logic.