Zero-Energy Bound State Trapped in Line-Shaped Vortex in Topological Superconductor

1. Introduction

Majorana bound states are unconventional zero-energy quasi-particles with non-Abelian statistics which are their own antiparticles. [1] These states have been proposed for storing quantum information [2] and for fault-tolerant quantum computing. The formation of Majorana bound states at the surface of a strong topological insulator (TI) with superconducting properties in the proximity of an s-wave superconductor (S) has been explored by Fu and Kane [3,4] and Sau et al. [5] Following Refs. [3] and [4] the surface states of a heterostructure in a TI with proximity-induced S were studied in different geometries. [5,6,7,8] The magnitude of the induced superconducting gap depends on the transparency of the TI/S junction and is in general smaller than the parent S gap. For low-lying excitations, the induced gap is a weak function of energy so that $\mathsf{\Delta }$ can be considered constant. Accordingly the problem is reduced to a 2D topological electron gas with superconducting BCS gap. An alternative platform to generate Majorana zero-energy modes in heterostructures consists of a semiconducting thin film sandwiched between a magnetic insulator and an s-wave superconductor. [9,10] In summary, a Majorana state can arise from a superconducting flux quantum trapped by a defect in a superconducting layer (Abrikosov vortex) [1,11,12,13,14], a layered heterostructure of a magnetic insulator, a semiconductor and a superconductor [9,10], in a tri-junction pair geometry [3,5] or a cylindrical cavity in a 2D superconductor. [7,8,15]

In this paper we investigate a flux quantum trapped in a line-shaped vortex core in a 2D topological superconductor. As a consequence of the geometry, all previously studied cases require polar coordinates, while in the present situation elliptic coordinates are more appropriate.

The electronic structure of in-gap states of a vortex in a 2D topological superconductor has been studied by numerous authors [5,7,8,15,16,17,18,19,20] by solving the corresponding Bogoliubov-de Gennes (BdG) equations. Key is the strong spin-orbit coupling leading to spin-momentum locking and a Berry phase of 1/2. This converts half-integer quantum numbers into integer ones and opens the possibility to create a zero-energy Majorana fermion. For a finite Fermi energy the low-energy excitations are equally spaced by an amount ${\mathsf{\Delta }}_{\infty }^2/E_F$, where ${\mathsf{\Delta }}_{\infty }$ is the superconducting gap far away from the vortex. The minigap from the Majorana mode to the first excitation can however be tuned by varying the Fermi energy and for $E_F\to 0$ the gap is of the order of ${\mathsf{\Delta }}_{\infty }$. [5,7,12]

Vortex states in superconducting graphene [21] are closely related to vortices with pointlike core in the present problem. The pseudospin parametrizing the two sublattices of the honeycomb lattice plays the role of the spin in the topological insulator. The BdG equations in superconducting graphene for energies close to the Dirac points (considering the spin, pseudospin and particle-holes) reduce to eight equations which decouple into two equivalent subsets of four equations each. Each subset is then equivalent to the present problem for a vortex trapped by a point defect with Dirac interaction. [22]

Analytical expressions for the low-energy bound state eigenvalues, the eigenfunctions and the local density of states (LDOS) for an isolated point vortex with strong spin-orbit coupling have been recently obtained [12,13,15] by solving the BdG equations using the method employed by Caroli, de Gennes and Matricon (CdeGM) [23,24] for a topological type II s-wave superconductor in terms of Bessel functions. The method consists of solving the BdG equations for small distances (compared to the correlation length $\xi$) from the core of the vortex, as well as for larger distances (still smaller than $\xi$). If the matching condition of these two solutions at an intermediate distance is independent of the distance from the vortex core, then we have a solution for the entire region of the vortex.

Experimentally, systems with large superconducting transition temperature are desirable. A zero-energy mode was detected via spin-selective Andreev reflections in the heterostructure Bi₂Te₃/NbSe₂[25] and with tunneling scanning spectroscopy in monolayers of the high-temperature superconductor FeTe_0.55Se_0.45. [26,27] To prove that the measured peak corresponds to a Majorana state requires the exchange of two such modes (qubit).

In the present paper we extend the calculation to a line-shaped vortex core (corresponding to a cut in the 2D superconductor confined between two foci) using elliptical coordinates rather than polar coordinates. The superconducting gap then vanishes along the entire line defect. In Sect. II we introduce the Rashba model and the BdG equations in elliptical coordinates, as well as the solution for the zero-energy Majorana state. The linear vortex traps one flux quantum. The general solution of the BdG equations is presented in the Appendix. Parametrizing the elliptical coordinates the problem is very similar to the solution of the BdG equations in spherical coordinates. Results are presented in Sect. III and Conclusions follow in Sect. IV.

Atomic line defects with zero-energy end-states have been observed in monolayer FeTe_0.5Se_0.5 high-$T_c$ superconductors. [28] This situation is different from the one treated here, where the cut in the superconducting layer joins the two end-points. The superconducting phase field winds around the position of each vortex, i.e. once around the linear cut in the present case, while once around the positions of each of the vortex cores in [28].

2. Rashba Interaction

2.1. The Model

The 2D electron gas represents the surface states of a 3D topological insulator. Superconductivity is induced via proximity by a type II s-wave superconductor. [5,29,30] The Rashba spin-orbit interaction locks the spin perpendicular to the momentum and the magnetic field is oriented perpendicular to the plane of the surface.

The wave function is a 4-component spinor, $\mathsf{\Psi }(\mathbf{r})={\left[{\psi }_{\uparrow }(\mathbf{r}){\psi }_{\downarrow }(\mathbf{r}){\psi }_{\uparrow }^{†}(\mathbf{r}){\psi }_{\downarrow }^{†}(\mathbf{r})\right]}^T$, and the Hamiltonian is given by $\mathcal{H}=\frac{1}{2}\int d^{2}r{\mathsf{\Psi }}^{†}(\mathbf{r}){\check{\mathcal{H}}}_B^{\perp }(\mathbf{r})\mathsf{\Psi }(\mathbf{r})$, where

$$\begin{array}{ccc}& & {\check{\mathcal{H}}}_B^{\perp }(\mathbf{r})=\left[\begin{array}{cc}\widehat{h}(\mathbf{r})& \widehat{\mathsf{\Delta }}(\mathbf{r})\\ -{\widehat{\mathsf{\Delta }}}^{\ast }(\mathbf{r})& -{\widehat{h}}^{\ast }(\mathbf{r})\end{array}\right],\hfill \end{array}$$

(1)

$$\begin{array}{ccc}& & \widehat{h}(\mathbf{r})=v_F\widehat{\sigma }·\left[\left(\mathbf{p}-{\displaystyle \frac{e}{c}}\mathbf{A}(\mathbf{r})\right)\times {\mathbf{e}}_z\right]-E_F\widehat{I},\hfill \end{array}$$

(2)

$$\begin{array}{ccc}& & \widehat{\mathsf{\Delta }}(\mathbf{r})=\mathsf{\Delta }(\mathbf{r})i{\widehat{\sigma }}_y,\hfill \end{array}$$

(3)

where ${\mathbf{e}}_z$ is the vector normal to the plane, $v_F$ the velocity in the Dirac cone, $\widehat{I}$ the unit matrix and $\mathbf{A}(\mathbf{r})$ the vector potential. Using similar arguments as in references [23,24] the vector potential and the magnetic field can be neglected in Equation (2).

2.2. Elliptic Coordinates

Since in this paper we consider a linearly shaped vortex, i.e. a linear cut in the 2D superconductor, it is appropriate to employ elliptical coordinates, [31] rather than the usual polar ones. Elliptic coordinates are a 2D orthogonal system in which the coordinate lines are confocal ellipses and hyperbolae. The two foci $F_1$ and $F_2$ are taken at $\pm a$ on the x-axis of the Cartesian coordinate system.

We define the elliptic coordinates $(\mu ,\nu )$ as

$$x=acosh(\mu )cos(\nu ),y=asinh(\mu )sin(\nu ),$$

(4)

where $\mu$ is a nonnegative real number and $\nu$ belongs to the interval $[0,2\pi )$. This can be rewritten on the complex plane as

$$X_{\pm }=x\pm iy=acosh(\mu \pm i\nu ).$$

(5)

Figure 1. (color online) Ellipses and hyperbolae with foci at $x=\pm a$ indicated by the blue crosses. The red straight line joining the two foci corresponds to $\mu =0$. The ellipses are parametrized by $\mu =0.2,0.4,0.6,0.8,1.0$ (solid curves), while the hyperbolae by $\nu =0,$$\pi /6,\pi /3,\pi /2,2\pi /3,5\pi /6,\pi ,7\pi /6,4\pi /3,3\pi /2,5\pi /3,11\pi /6$ (dashed curves).

Figure 1. (color online) Ellipses and hyperbolae with foci at $x=\pm a$ indicated by the blue crosses. The red straight line joining the two foci corresponds to $\mu =0$. The ellipses are parametrized by $\mu =0.2,0.4,0.6,0.8,1.0$ (solid curves), while the hyperbolae by $\nu =0,$$\pi /6,\pi /3,\pi /2,2\pi /3,5\pi /6,\pi ,7\pi /6,4\pi /3,3\pi /2,5\pi /3,11\pi /6$ (dashed curves).

These definitions correspond to ellipses and hyperbolae for constant $\mu$ and $\nu$, respectively, i.e.

$$\begin{array}{ccc}\hfill {\displaystyle \frac{x^2}{a^2{cosh}^2(\mu )}}+{\displaystyle \frac{y^2}{a^2{sinh}^2(\mu )}}& =& 1,\hfill \\ \hfill {\displaystyle \frac{x^2}{a^2{cos}^2(\nu )}}-{\displaystyle \frac{y^2}{a^2{sin}^2(\nu )}}& =& 1.\hfill \end{array}$$

(6)

An alternative definition of elliptical coordinates is $(\sigma ,\tau )$, where $\sigma =cosh(\mu )$ and $\tau =cos(\nu )$. The curves of constant $\sigma$ are then ellipses and those of constant $\tau$ hyperbolae. Hence, the coordinate $\tau$ belongs to the interval $[-1,1]$, whereas $\sigma$ must be greater than or equal to one. There is a simple geometrical relation of the distances to the foci $F_1$ and $F_2$. For any point in the plane, the sum $r_1+r_2$ of its distances to the foci equals $2a\sigma$, whereas their difference $r_1-r_2$ equals $2a\tau$. Thus, the distance to $F_1$ is $a(\sigma +\tau$), and the distance to $F_2$ is $a(\sigma -\tau )$.

According to Abrikosov [32] the radial dependence of an isolated superconducting flux point vortex is approximately $\mathsf{\Delta }(r)={\mathsf{\Delta }}_{\infty }tanh(r/\xi )$. For a distance r much shorter than the coherence length $\xi$, $\mathsf{\Delta }(r)$ is linear in r, if the vorticity is one. The sum of the distances to both foci is then $r_1+r_2=2a\sigma$ and due to the linearity in r, the sum of the pair potential for both vortices is $\mathsf{\Delta }(r_1)+\mathsf{\Delta }(r_2)-\mathsf{\Delta }(2a)=\mathsf{\Delta }\left[2a(\sigma -1)\right]$, i.e. it does not depend on $\tau$ and vanishes for $\mu =0$.

2.3. BdG Equations in Elliptic Coordinates

The momentum can be expressed in terms of the derivatives with respect to $X_{\pm }$, i.e.

$${\displaystyle \frac{\partial }{\partial X_{\pm }}}=-{\displaystyle \frac{i}{asinh(\mu \pm i\nu )}}\left[i{\displaystyle \frac{\partial }{\partial \mu }}\pm {\displaystyle \frac{\partial }{\partial \nu }}\right],$$

(7)

so that $\widehat{h}(\mathbf{r})=$

$$\begin{array}{c}\hfill \left[\begin{array}{cc}-E_F\hfill & -{\displaystyle \frac{i{v}_F}{asinh(\mu +i\nu )}}\left(i{\displaystyle \frac{\partial }{\partial \mu }}+{\displaystyle \frac{\partial }{\partial \nu }}\right)\hfill \\ {\displaystyle \frac{i{v}_F}{asinh(\mu -i\nu )}}\left(i{\displaystyle \frac{\partial }{\partial \mu }}-{\displaystyle \frac{\partial }{\partial \nu }}\right)\hfill & -E_F\hfill \end{array}\right].\end{array}$$

(8)

If only one flux quantum is contained in the vortex, we may write $\widehat{\mathsf{\Delta }}(\mathbf{r})=\mathsf{\Delta }\left[2a(\sigma -1)\right]i{\widehat{\sigma }}_y{e}^{-i\nu }$. The $e^{-i\nu }$ phase of $\widehat{\mathsf{\Delta }}(\mathbf{r})$ can be eliminated via a gauge transformation.

The field operators are expanded as

$$\psi (\mu ,\nu )={\displaystyle \frac{1}{\sqrt{2\pi }}}\sum _m{\psi }_m(\mu )e^{im\nu },$$

(9)

where m is an integer for vorticity one and similarly

$${\displaystyle \frac{1}{sinh(\mu \pm i\nu )}}=2\sum _{n=0}^{\infty }exp[-(2n+1)\mu \mp i(2n+1)\nu ].$$

(10)

The spinor $\mathsf{\Psi }$ is required to be a single-valued function of $\nu$, so that

$$\begin{array}{c}\hfill {\mathsf{\Psi }}_m(\mu ,\nu )=\left[\begin{array}{c}{\psi }_{\uparrow m}(\mu )e^{-i\nu }\\ {\psi }_{\downarrow m}(\mu )\\ {\psi }_{\uparrow m}^{†}(\mu )e^{i\nu }\\ {\psi }_{\downarrow m}^{†}(\mu )\end{array}\right].\end{array}$$

(11)

Applying the spinor to $\widehat{h}(\mathbf{r})$ we obtain

$$\begin{array}{c}\hfill \left[\begin{array}{cc}-E_F& {\displaystyle \frac{v_F}{asinh(\mu +i\nu )}}\left({\displaystyle \frac{\partial }{\partial \mu }}+m\right)\\ -{\displaystyle \frac{v_F}{asinh(\mu -i\nu )}}\left({\displaystyle \frac{\partial }{\partial \mu }}-m+1\right)& -E_F\end{array}\right]\left[\begin{array}{c}{\psi }_{\uparrow m}(\mu )e^{i(m-1)\nu }\\ {\psi }_{\downarrow m}(\mu )e^{im\nu }\end{array}\right].\end{array}$$

(12)

and

$$\begin{array}{c}\hfill \left[\begin{array}{cc}{E}_F& -{\displaystyle \frac{v_F}{asinh(\mu -i\nu )}}\left({\displaystyle \frac{\partial }{\partial \mu }}+m\right)\\ {\displaystyle \frac{v_F}{asinh(\mu +i\nu )}}\left({\displaystyle \frac{\partial }{\partial \mu }}-m-1\right)& E_F\end{array}\right]\left[\begin{array}{c}{\psi }_{\uparrow m}^{†}(\mu )e^{i(m+1)\nu }\\ {\psi }_{\downarrow m}^{†}(\mu )e^{im\nu }\end{array}\right].\end{array}$$

(13)

Inserting $\widehat{\mathsf{\Delta }}(\mathbf{r})$ and the Taylor expansion of $1/sinh(\mu \pm i\nu )$ [Equation (10)] the differential equation yields

$$\begin{array}{ccc}\hfill \sum _{n=0}^{\infty }& & \left[\begin{array}{cccc}-(E_F+E){\delta }_{n,0}& Z_n\left({\displaystyle \frac{\partial }{\partial \mu }}+m\right)& 0& \overline{\mathsf{\Delta }}{e}^{-i\nu }{\delta }_{n=0}\\ -Z_n^{\ast }\left({\displaystyle \frac{\partial }{\partial \mu }}-m+1\right)& -(E_F+E){\delta }_{n,0}& -\overline{\mathsf{\Delta }}{e}^{-i\nu }{\delta }_{n=0}& 0\\ 0& \overline{\mathsf{\Delta }}{e}^{i\nu }{\delta }_{n=0}& (E_F-E){\delta }_{n,0}& -Z_n^{\ast }\left({\displaystyle \frac{\partial }{\partial \mu }}+m\right)\\ -\overline{\mathsf{\Delta }}{e}^{i\nu }{\delta }_{n=0}& 0& Z_n\left({\displaystyle \frac{\partial }{\partial \mu }}-m-1\right)& (E_F-E){\delta }_{n,0}\end{array}\right]\hfill \end{array}$$

$$\begin{array}{c}\hfill \times \left[\begin{array}{c}{\psi }_{\uparrow mn}(\mu )e^{i(m-1)\nu }\\ {\psi }_{\downarrow mn}(\mu )e^{im\nu }\\ {\psi }_{\uparrow mn}^{†}(\mu )e^{i(m+1)\nu }\\ {\psi }_{\downarrow mn}^{†}(\mu )e^{im\nu }\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right],\end{array}$$

(14)

where $Z_n$ denotes ${\displaystyle \frac{2v_F}{a}}{e}^{-(2n+1)(\mu +i\nu )}$ and $\overline{\mathsf{\Delta }}=\mathsf{\Delta }\left[2a(\sigma -1)\right]$, with $\sigma$ being $cosh(\mu )$. Here we have incorporated the subindex n into the field operator $\psi$. It is easy to verify that ${\psi }_{\uparrow mn}={\psi }_{\downarrow mn}={\psi }_{\uparrow mn}^{†}={\psi }_{\downarrow mn}^{†}=0$ for $n\ne 0$. Only $n=0$ plays a relevant role and the subindex n can be suppressed. Below we denote the components of the spinor $\mathsf{\Psi }$ as $f_j^m(\mu )$ with $j=1,\cdots ,4$, where the index j contains the spin (up or down) and annihilation and creation symbols. The $\nu$-dependence of $f_j^m$ is $exp[-i\nu {\widehat{\tau }}_z(1+{\widehat{\sigma }}_z/2+im\nu ]$, where ${\widehat{\sigma }}_z$ and ${\widehat{\tau }}_z$ are Pauli matrices for the spin and particle/hole sectors, respectively.

The matrix equation reduces to four linearly coupled differential equations, in which the $\nu$ dependence cancels out and the solution is single valued,

$$\begin{array}{ccc}& & {\displaystyle \frac{2v_F}{a}}{e}^{-\mu }\left({\displaystyle \frac{\partial }{\partial \mu }}-m+1\right)f_1^m(\mu )+(E+E_F)f_2^m(\mu )+\overline{\mathsf{\Delta }}{f}_3^m(\mu )=0,\hfill \end{array}$$

(15)

$$\begin{array}{ccc}& & {\displaystyle \frac{2v_F}{a}}{e}^{-\mu }\left({\displaystyle \frac{\partial }{\partial \mu }}+m\right)f_2^m(\mu )-(E+E_F)f_1^m(\mu )+\overline{\mathsf{\Delta }}{f}_4^m(\mu )=0,\hfill \end{array}$$

(16)

$$\begin{array}{ccc}& & {\displaystyle \frac{2v_F}{a}}{e}^{-\mu }\left({\displaystyle \frac{\partial }{\partial \mu }}+m+1\right)f_3^m(\mu )-(E-E_F)f_4^m(\mu )+\overline{\mathsf{\Delta }}{f}_1^m(\mu )=0,\hfill \end{array}$$

(17)

$$\begin{array}{ccc}& & {\displaystyle \frac{2v_F}{a}}{e}^{-\mu }\left({\displaystyle \frac{\partial }{\partial \mu }}-m\right)f_4^m(\mu )+(E-E_F)f_3^m(\mu )+\overline{\mathsf{\Delta }}{f}_2^m(\mu )=0.\hfill \end{array}$$

(18)

Defining $\rho =e^{\mu }a/2$, $E_F=v_F{k}_F$, $\tilde{\rho }=\rho k_F$, $\tilde{E}=E/E_F$ and $\tilde{\mathsf{\Delta }}=\overline{\mathsf{\Delta }}/E_F$, the equations are dimensionless and have the following form

$$\begin{array}{ccc}& & \left({\displaystyle \frac{\partial }{\partial \tilde{\rho }}}-{\displaystyle \frac{m-1}{\tilde{\rho }}}\right)f_1^m(\tilde{\rho })+(1+\tilde{E})f_2^m(\tilde{\rho })+\tilde{\mathsf{\Delta }}{f}_3^m(\tilde{\rho })=0,\hfill \end{array}$$

(19)

$$\begin{array}{ccc}& & \left({\displaystyle \frac{\partial }{\partial \tilde{\rho }}}+{\displaystyle \frac{m}{\tilde{\rho }}}\right)f_2^m(\tilde{\rho })-(1+\tilde{E})f_1^m(\tilde{\rho })+\tilde{\mathsf{\Delta }}{f}_4^m(\tilde{\rho })=0,\hfill \end{array}$$

(20)

$$\begin{array}{ccc}& & \left({\displaystyle \frac{\partial }{\partial \tilde{\rho }}}+{\displaystyle \frac{m+1}{\tilde{\rho }}}\right)f_3^m(\tilde{\rho })+(1-\tilde{E})f_4^m(\tilde{\rho })+\tilde{\mathsf{\Delta }}{f}_1^m(\tilde{\rho })=0,\hfill \end{array}$$

(21)

$$\begin{array}{ccc}& & \left({\displaystyle \frac{\partial }{\partial \tilde{\rho }}}-{\displaystyle \frac{m}{\tilde{\rho }}}\right)f_4^m(\tilde{\rho })-(1-\tilde{E})f_3^m(\tilde{\rho })+\tilde{\mathsf{\Delta }}{f}_2^m(\tilde{\rho })=0.\hfill \end{array}$$

(22)

These equations are formally similar to those of a single pinned point vortex [7,8,12] described by polar coordinates. However $\rho (=e^{\mu }a/2)$ is not the radial coordinate so that the meaning of these equations in Cartesian coordinates is substantially different. $\nu$ and $\theta$ are both angular coordinates.

2.4. Majorana State

The Majorana bound state has zero energy and corresponds to $m=0$. It is easily verified that the solution of Equations (19)–(22) is given by

$$\begin{array}{ccc}& & f_1^0(\tilde{\rho })=C{J}_1(\tilde{\rho })e^{-K(\tilde{\rho })},\hfill \\ & & f_2^0(\tilde{\rho })=-C{J}_0(\tilde{\rho })e^{-K(\tilde{\rho })},\hfill \\ & & f_3^0(\tilde{\rho })=C{J}_1(\tilde{\rho })e^{-K(\tilde{\rho })},\hfill \\ & & f_4^0(\tilde{\rho })=-C{J}_0(\tilde{\rho })e^{-K(\tilde{\rho })},\hfill \end{array}$$

(23)

where $J_n$ are Bessel functions of integer order, C is a normalization constant and $K(\tilde{\rho })$ is given by

$$\begin{array}{ccc}& & {\displaystyle \frac{dK}{d\tilde{\rho }}}={\displaystyle \frac{\mathsf{\Delta }[2a(\sigma -1)]}{E_F}},\hfill \\ & & K(\tilde{\rho })={\int }_{a{k}_F/2}^{\tilde{\rho }}d{\tilde{\rho }}^{\prime }\mathsf{\Delta }\left[{\displaystyle \frac{2{\tilde{\rho }}^{\prime }}{k_F}}+{\displaystyle \frac{a^2{k}_F}{2{\tilde{\rho }}^{\prime }}}-2a\right]/E_F.\hfill \end{array}$$

(24)

Note that the definition of $K(\tilde{\rho })$ is different from the corresponding expression for a point vortex in polar coordinates. [12,13] The Majorana wave function is

$$\begin{array}{ccc}\hfill {\mathsf{\Psi }}_M& =& C\int d^{2}r{e}^{-K(\tilde{\rho })}[J_1(\tilde{\rho })e^{-i\nu }{\psi }_{\uparrow }(\mathbf{r})-J_0(\tilde{\rho }){\psi }_{\downarrow }(\mathbf{r})\hfill \\ & +& J_1(\tilde{\rho })e^{i\nu }{\psi }_{\uparrow }^{†}(\mathbf{r})-J_0(\tilde{\rho }){\psi }_{\downarrow }^{†}(\mathbf{r})].\hfill \end{array}$$

(25)

A Majorana state is self-adjoint and should satisfy ${\mathsf{\Psi }}_M={\mathsf{\Psi }}_M^{†}$. The adjoint of the first term in (25) yields the third term and the adjoint of the second term yields the fourth term.

3. Results

3.1. Excitation Energies

The solution of the BdG equations for the general case is outlined in the Appendix, where we have shown that the wave functions of the fermion excited states consist of two factors,

one being a Bessel function and the second one determines the exponential decay of the wave function. The latter, $exp(-K(\tilde{\rho }))$, is a function of the superconducting order parameter. The energies are equally spaced and depend on the function $K(\tilde{\rho })$, Equation (24). Since the superconducting gap is a linear function of distance with slope ${\mathsf{\Delta }}_{\infty }/\xi$, where $\xi$ is the coherence length, $K(\tilde{\rho })$ is easily evaluated yielding

$$\begin{array}{ccc}\hfill K(\tilde{\rho })& =& {\left({\displaystyle \frac{{\mathsf{\Delta }}_{\infty }}{E_F}}\right)}^2[{\tilde{\rho }}^2-2a{k}_F\tilde{\rho }\hfill \\ & & +{\displaystyle \frac{{(a{k}_F)}^2}{2}}\left(ln(2\tilde{\rho }/a{k}_F)+{\displaystyle \frac{3}{2}}\right)].\hfill \end{array}$$

(26)

As shown in Figure 2 the function $K(\tilde{\rho })$ is approximately parabolic and decreases with increasing the distance between the foci $2a$.

The energies are given by Equation (A34), where the matching point ${\tilde{\rho }}_c$ at the lower integration limit can be approximated by $a{k}_F/2$. It needs to be verified that as a function of the upper integration limit the value of the $E_m$ converges. This convergence is shown in Figure 3 for the same values of $a{k}_F$ as in Figure 2. The energy saturates into plateaus with increasing $\tilde{\rho }$. The gap between consecutive excited states decreases with increasing $a{k}_F$. For given $a{k}_F$ the excited states are equidistant, i.e. $E_m=m{E}_1$, where $m=0$ corresponds to the Majorana state.

3.2. Local Density of Bound States

The one-particle Green’s function for spin-component $\sigma$ is given by the wave functions and energy eigenvalues [37]

$$G_{\omega }^{\sigma }(\mathbf{r},{\mathbf{r}}^{\prime })=\sum _m\left({\displaystyle \frac{u_m^{\sigma }(\mathbf{r})u_m^{\sigma \ast }{(\mathbf{r})}^{\prime }}{\omega +i0-E_m}}+{\displaystyle \frac{v_m^{\sigma }(\mathbf{r})v_m^{\sigma \ast }{(\mathbf{r})}^{\prime }}{\omega +i0+E_m}}\right),$$

(27)

where u and v correspond to the functions $f_j^m(\mathbf{r})$, i.e. the wave functions for particles and holes. The local density of states (LDOS) for spin projection $\sigma$ is given by the imaginary part of the Green’s function for ${\mathbf{r}}^{\prime }\to \mathbf{r}$

$$N_{\sigma }(\omega ,\tilde{\rho })=\sum _m{N}_{m\sigma }(\omega ,\tilde{\rho })=\sum _m\left(|u_m^{\sigma }(\tilde{\rho }){|}^2\delta (\omega -E_m)+{\left|v_m^{\sigma }(\tilde{\rho })\right|}^2\delta (\omega +E_m)\right).$$

(28)

At finite temperature the $\delta$-function is substituted by minus the derivative of the Fermi function, and hence the peaks broaden with increasing temperature. [17]

The LDOS for our model with Rashba interaction is proportional to (up to a normalization constant for the spinor)

$$\begin{array}{ccc}\hfill N_{\uparrow }(\omega ,\tilde{\rho })& \propto & \sum _m\left[J_{m-1}{(\tilde{\rho })}^2\delta (\omega -E_m)+J_{m+1}{(\tilde{\rho })}^2\delta (\omega +E_m)\right]e^{-2K(\tilde{\rho })},\hfill \\ \hfill N_{\downarrow }(\omega ,\tilde{\rho })& \propto & \sum _m\left[J_m{(\tilde{\rho })}^2\delta (\omega -E_m)+J_m{(\tilde{\rho })}^2\delta (\omega +E_m)\right]e^{-2K(\tilde{\rho })}.\hfill \end{array}$$

(29)

The above expressions are formally identical to the ones for a point vortex with two important differences, namely, the function $K(\tilde{\rho })$ is not the same and in the present case $\tilde{\rho }$ is not the distance to the origin of the vortex, but a function of elliptical coordinates. The zero energy Majorana state corresponds to $m=0$ and $E_0=0$. The total LDOS ($N=N_{\uparrow }+N_{\downarrow }$) is invariant under the transformation $m\to -m$, i.e. $E_m\to -E_m$, and $u\leftrightarrow v$, i.e. the simultaneous interchange of the particle and hole amplitudes. It is interesting to notice that as for a point vortex the LDOS is not particle-hole symmetric.

To understand the lack of particle-hole symmetry we consider Equations (12) and (13) and arguments similar to Ref. [12]. For down-spin electrons the effective angular momentum is m, while for up-spin electrons it is $(m-1)$. On the other hand, for down-spin holes it is still m, but for up-spin holes it is $(m+1)$. As a consequence of the difference between electrons and holes with up-spin, the particle-hole symmetry is broken. The origin of this asymmetry is the Rashba spin-orbit interaction. Hence, the LDOS is asymmetric as a function of $\omega$.

The LDOS is measurable via scanning tunneling microscopy (STM) by fine-tuning the energy $\omega$ at a distance from the core of the linear vortex, parametrized by $\tilde{\rho }$.

The LDOS as a function of energy in units of the gap $E_1$ is shown in Figure 4 for three values of $\tilde{\rho }$. The intensity of the peaks depends on (i) the Bessel functions, (ii) the exponential decay due to $exp[-K(\tilde{\rho })]$, and (iii) on the separation of the two foci, $2a$. The LDOS also decreases rapidly as a function of temperature. The Majorana state ($m=0$) has always the same intensity as the $m=1$ peak. In general, the LDOS is symmetric about $\omega =0.5$. For $\tilde{\rho }=0$ and zero distance between the foci, $a=0$, i.e. for a simple point vortex, there are only two peaks, namely $m=0$ and 1, as a consequence of the Bessel functions, $J_m(\tilde{\rho })$, which are zero as $\tilde{\rho }\to 0$ for $m\ne 0$. More peaks appear at finite $\tilde{\rho }$ and/or finite $a{k}_F$, but for larger $\tilde{\rho }$ or $a{k}_F$ the intensity does not necessarily decrease monotonically with m as a consequence of the oscillations of the Bessel functions.

For example in Figure 4c for $a{k}_F=0.5$ and $\tilde{\rho }=2$ the $m=2$ peak has higher intensity than the $m=1$ state. This can be understood with the plot of the intensities as a function of $\tilde{\rho }$ displayed in Figure 5. For $\tilde{\rho }=2$ the $m=0$ and 1 peaks have smaller intensity than the $m=2$ and 3 peaks. But for larger energies the intensity of the peaks decreases rapidly. For $a=0$ the main maxima of the LDOS form concentric circles about the center of the core with their radius increasing with m. For $a\ne 0$, on the other hand, the main maxima of the LDOS follow ellipses.

The intensities of the peaks are shown in Figure 5 as a function of $\tilde{\rho }$ and for two values of $a{k}_F$. Note that $\tilde{\rho }\ge a{k}_F/2$. The peaks are denoted by two consecutive m values indicating, e.g. $[J_m{(\tilde{\rho })}^2+J_{m+1}{(\tilde{\rho })}^2]exp(-2K\tilde{\rho })$; the $m=0,1$ curve is black, the $m=1,2$ curve is red, the $m=2,3$ curve is green and the $m=3,4$ curve is blue. As already mentioned, the oscillations arise from the Bessel functions.

3.3. Spin Polarization

Of interest is also the spin-polarization of the peaks in the LDOS. For $a\to 0$ (point vortex) and $\tilde{\rho }=0$ the Majorana state is purely a down-spin peak. With increasing a and $\tilde{\rho }$ a growing up-spin component arises. On the other hand, the $m=1$ peak is predominantly up-spin, with growing down-spin components as a and $\tilde{\rho }$ increase. This pattern is shown in Figure 6 and can be expressed as $N_{up}(\omega ,\tilde{\rho })=N_{down}(\omega -1,\tilde{\rho })$ with $\omega$ in units of the excitation energy $E_1$. Note that $N_{down}(\omega ,\tilde{\rho })$ is a symmetric function of $\omega$, while $N_{up}(\omega ,\tilde{\rho })$ is even about $\omega =1$, and the sum of $N_{up}(\omega ,\tilde{\rho })$ and $N_{down}(\omega ,\tilde{\rho })$ yields the total LDOS, $N(\omega ,\tilde{\rho })$, which is symmetric about $\omega =0.5$.

3.4. Orbits in Cartesian Coordinates

So far we investigated the solution of the bound states in terms of the natural elliptic coordinates, $\tilde{\rho }$ and $\nu$. These variables do not correspond to the real space Cartesian coordinates x and y, defined by the elliptic coordinates Equation (4). We now consider the zero-energy Majorana state, which is given by the zero-order (down-spin density) and first-order (up-spin density) Bessel functions. The oscillations of the Bessel functions give rise to maxima and minima, which are shown in Figure 7a as a function of $\tilde{\rho }$ for $a{k}_F=5.0$. Here down-spin maxima are displayed in red and up-spin maxima in black. Note that solutions only exist for $\tilde{\rho }\ge a{k}_F/2$. As a function of $\tilde{\rho }$ the down-spin and up-spin maxima alternate and are approximately equally spaced.

Equation (4) transforms circular solution (natural coordinated for this problem) into elliptical real space coordinates. Note that for $\tilde{\rho }=a{k}_F/2$ the real space pattern is a horizontal segment between $-a$ and a, which corresponds to a vanishing superconducting gap. The order parameter grows linearly away from this segment. The down- and up-spin maxima ellipses alternate and grow with increasing $\tilde{\rho }$ (see Figure 7b).

4. Concluding Remarks

We investigated the bound states in the vortex core of a line-shaped cut of length $2a$ in a two-dimensional topological superconductor with Rashba interaction by solving the BdG equations following the analytic method outlined by Caroli et al. [23] The electron gas corresponds to the surface states of a 3D TI with proximity induced superconductivity from a nearby s wave superconductor. [3,6,7,9] This way the model is reduced to an effective 2D superconductor with a linear pinning defect between the two foci. The superconductor is gapless along the linear defect. The strong spin-orbit interaction locks the momentum and the spin perpendicularly. The characteristic energy scale for the spacing of the energy levels in the vortex is proportional to ${\mathsf{\Delta }}_{\infty }^2/E_F$.

The calculation yields a string of fermion bound states with energy $E_m$, $m\ne 0$, and a bound state with Majorana statistics for $m=0$ and $E=0$. In the present case the pinning defect has a linear shape of length $2a$. The results are similar to those of a point defect studied previously, [12,13,15] except that instead of concentric circles the curves of constant energy now are ellipses. The gapless region in the present case is the segment of length $2a$, while for the point defect the gap only vanishes at the singular point. The analytical expressions for the wave functions, consist of products of a Bessel function and an exponential decay, $exp(-2K(\tilde{\rho }))$, as function of the parameter $\tilde{\rho }$ (which for an ellipse is different than the radius of the circle). Given the wave functions we obtained an analytic expression of the LDOS for the bound states. This quantity is experimentally accessible via STM. [26]

We investigated the LDOS as a function of external frequency, temperature, the distance from the center of the vortex and the length $2a$ between the foci. The intensity of the peaks in the LDOS rapidly decreases with temperature due to the smearing of the Fermi function. It is interesting to notice that the particle-hole symmetry is broken in the LDOS as a consequence of the spin-orbit coupling. This can be traced to Equation (2). For instance, the Majorana peak is a pure down-spin state at the core of a point vortex, but becomes partially polarized away from the center and for finite a.

The main difference between the ordinary superconductor and the topological superconducting gas is the spin-locking. In the latter in a closed path the spin is forced to follow the momentum giving rise to a non-trivial Berry phase of 1/2. This converts the half-integer quantum numbers into integer ones and opens the possibility to the existence of a Majorana fermion.

The experimental search for zero-energy modes was successful in two systems, namely, heterostructures of Bi₂Te₃/NbSe₂[25,33,34] and monolayers of the high-temperature superconductor FeTe_0.55Se_0.45 on SrTiO₃(0011).[26,27,28,35] Both systems have the advantage of relatively large superconducting transition temperatures. A zero-energy mode, however does not imply a Majorana state, since it may as well originate from a magnetic impurity or an Andreev reflection. The verification that a zero-energy mode corresponds to a Majorana state requiresbraiding, i.e. exchanging the positions of two Majorana zero modes.

Waterfall-like STM spectra have been observed in the iron superconductors [26,28] and were theoretically studied in [36]. Zero-energy bound states simultaneously appearing at both ends of a 1D atomic line defect were discovered in [28]. These line defects resemble the ones studied in this paper.

Appendix A.1. Second Order Differential Equations

We first convert the first order differential equations, Equations (19)–(22), into second order differential equations.[12,13,23,24] For instance, we express $f_2^m$ from Equation (19) and insert it into Equation (20) and similar substitutions for $f_1^m$, $f_3^m$ and $f_4^m$. Defining $q_p=1+\tilde{E}$ and $q_h=1-\tilde{E}$ for particles and holes, respectively, we obtain

$$\begin{array}{ccc}& & \left[{\displaystyle \frac{{\partial }^2}{\partial {\tilde{\rho }}^2}}+{\displaystyle \frac{1}{\tilde{\rho }}}{\displaystyle \frac{\partial }{\partial \tilde{\rho }}}-{\displaystyle \frac{{(m-1)}^2}{{\tilde{\rho }}^2}}+q_p^2\right]f_1^m(\tilde{\rho })=q_p\tilde{\mathsf{\Delta }}{f}_4^m(\tilde{\rho })-\left({\displaystyle \frac{\partial }{\partial \rho }}+{\displaystyle \frac{m}{\tilde{\rho }}}\right)\tilde{\mathsf{\Delta }}{f}_3^m(\tilde{\rho }),\hfill \end{array}$$

(A1)

$$\begin{array}{ccc}& & \left[{\displaystyle \frac{{\partial }^2}{\partial {\tilde{\rho }}^2}}+{\displaystyle \frac{1}{\tilde{\rho }}}{\displaystyle \frac{\partial }{\partial \tilde{\rho }}}-{\displaystyle \frac{m^2}{{\tilde{\rho }}^2}}+q_p^2\right]f_2^m(\tilde{\rho })=-q_p\tilde{\mathsf{\Delta }}{f}_3^m(\tilde{\rho })-\left({\displaystyle \frac{\partial }{\partial \rho }}-{\displaystyle \frac{m-1}{\tilde{\rho }}}\right)\tilde{\mathsf{\Delta }}{f}_4^m(\tilde{\rho }),\hfill \end{array}$$

(A2)

$$\begin{array}{ccc}& & \left[{\displaystyle \frac{{\partial }^2}{\partial {\tilde{\rho }}^2}}+{\displaystyle \frac{1}{\tilde{\rho }}}{\displaystyle \frac{\partial }{\partial \tilde{\rho }}}-{\displaystyle \frac{{(m+1)}^2}{{\tilde{\rho }}^2}}+q_h^2\right]f_3^m(\tilde{\rho })=q_h\tilde{\mathsf{\Delta }}{f}_2^m(\tilde{\rho })-\left({\displaystyle \frac{\partial }{\partial \rho }}-{\displaystyle \frac{m}{\tilde{\rho }}}\right)\tilde{\mathsf{\Delta }}{f}_1^m(\tilde{\rho }),\hfill \end{array}$$

(A3)

$$\begin{array}{ccc}& & \left[{\displaystyle \frac{{\partial }^2}{\partial {\tilde{\rho }}^2}}+{\displaystyle \frac{1}{\tilde{\rho }}}{\displaystyle \frac{\partial }{\partial \tilde{\rho }}}-{\displaystyle \frac{m^2}{{\tilde{\rho }}^2}}+q_h^2\right]f_4^m(\tilde{\rho })=-q_h\tilde{\mathsf{\Delta }}{f}_1^m(\tilde{\rho })-\left({\displaystyle \frac{\partial }{\partial \rho }}+{\displaystyle \frac{m+1}{\tilde{\rho }}}\right)\tilde{\mathsf{\Delta }}{f}_2^m(\tilde{\rho }).\hfill \end{array}$$

(A4)

Appendix A.2. Solution for $\tilde{\rho }$ < $\tilde{\rho }c$

Since $2a\ll \xi$, where $\xi$ is the coherence length, for $\tilde{\rho }<{\tilde{\rho }}_c$ we may neglect $\tilde{\mathsf{\Delta }}$ in the core of the vortex. The equations for $\tilde{\rho }<{\tilde{\rho }}_c$ are then of the integer Bessel function type

$$\begin{array}{ccc}& & \left[{\displaystyle \frac{{\partial }^2}{\partial {\tilde{\rho }}^2}}+{\displaystyle \frac{1}{\tilde{\rho }}}{\displaystyle \frac{\partial }{\partial \tilde{\rho }}}-{\displaystyle \frac{{(m-1)}^2}{{\tilde{\rho }}^2}}+q_p^2\right]f_1^m(\tilde{\rho })=0,\hfill \end{array}$$

(A5)

$$\begin{array}{ccc}& & \left[{\displaystyle \frac{{\partial }^2}{\partial {\tilde{\rho }}^2}}+{\displaystyle \frac{1}{\tilde{\rho }}}{\displaystyle \frac{\partial }{\partial \tilde{\rho }}}-{\displaystyle \frac{m^2}{{\tilde{\rho }}^2}}+q_p^2\right]f_2^m(\tilde{\rho })=0,\hfill \end{array}$$

(A6)

$$\begin{array}{ccc}& & \left[{\displaystyle \frac{{\partial }^2}{\partial {\tilde{\rho }}^2}}+{\displaystyle \frac{1}{\tilde{\rho }}}{\displaystyle \frac{\partial }{\partial \tilde{\rho }}}-{\displaystyle \frac{{(m+1)}^2}{{\tilde{\rho }}^2}}+q_h^2\right]f_3^m(\tilde{\rho })=0,\hfill \end{array}$$

(A7)

$$\begin{array}{ccc}& & \left[{\displaystyle \frac{{\partial }^2}{\partial {\tilde{\rho }}^2}}+{\displaystyle \frac{1}{\tilde{\rho }}}{\displaystyle \frac{\partial }{\partial \tilde{\rho }}}-{\displaystyle \frac{m^2}{{\tilde{\rho }}^2}}+q_h^2\right]f_4^m(\tilde{\rho })=0.\hfill \end{array}$$

(A8)

For $m\ge 1$, the solutions for $\tilde{\rho }<{\tilde{\rho }}_c$ are then

$$\begin{array}{ccc}& & f_1^m(q_p\tilde{\rho })=A_1^m{J}_{m-1}(q_p\tilde{\rho }),\hfill \end{array}$$

(A9)

$$\begin{array}{ccc}& & f_2^m(q_p\tilde{\rho })=A_2^m{J}_m(q_p\tilde{\rho }),\hfill \end{array}$$

(A10)

$$\begin{array}{ccc}& & f_3^m(q_h\tilde{\rho })=A_3^m{J}_{m+1}(q_h\tilde{\rho }),\hfill \end{array}$$

(A11)

$$\begin{array}{ccc}& & f_4^m(q_h\tilde{\rho })=A_4^m{J}_m(q_h\tilde{\rho }).\hfill \end{array}$$

(A12)

The constants $A_j^m$ are not all independent. Inserting the solutions into the first order differential equations (still with $\mathsf{\Delta }=0$) we obtain $A_1^m=A_2^m$ and $A_3^m=-A_4^m$. $A_1^m$ and $A_4^m$ are independent for $\mathsf{\Delta }=0$, but are coupled for $\mathsf{\Delta }\ne 0$, so that $A_1=\sqrt{1+\tilde{E}}$ and $A_4=\sqrt{1-\tilde{E}}$ (see Appendix 4).

Appendix A.3. Solution for $\tilde{\rho }$ > $\tilde{\rho }c$

In this case $\mathsf{\Delta }$ plays an important role. An Ansatz for a solution of the second order differential equations is in terms of a Hankel function of the first kind times an envelope function $g_j$ [37]

$$f_j^m(\tilde{\rho })=B_j[H_m^{(1)}(\tilde{\rho })g_j(\tilde{\rho })+c.c.],$$

(A13)

where $\rho \ll \xi$ and the $B_j$ are constants. We postulate $B_1=B_2=-B_3=B_4=\frac{1}{2}B$, where the normalization constant B is equated to one for simplicity. The relative phases of the $B_j$ are the same as in Equations (A9)-(A12). This choice of $B_j$ is verified below where the wave function for $\tilde{\rho }<{\tilde{\rho }}_c$ is matched to the one for $\tilde{\rho }>{\tilde{\rho }}_c$.

It is convenient to rewrite Equation (A1) as

$$\begin{array}{ccc}& & \left[{\displaystyle \frac{{\partial }^2}{\partial {\tilde{\rho }}^2}}+{\displaystyle \frac{1}{\tilde{\rho }}}{\displaystyle \frac{\partial }{\partial \tilde{\rho }}}-{\displaystyle \frac{m^2}{{\tilde{\rho }}^2}}+1\right]f_1^m(\tilde{\rho })+\left({\tilde{E}}^2+2\tilde{E}+{\displaystyle \frac{2m-1}{{\tilde{\rho }}^2}}\right)f_1^m(\tilde{\rho })\hfill \\ & & =(1+\tilde{E})\tilde{\mathsf{\Delta }}{f}_4^m(\tilde{\rho })-\left({\displaystyle \frac{\partial }{\partial \rho }}+{\displaystyle \frac{m}{\tilde{\rho }}}\right)\tilde{\mathsf{\Delta }}{f}_3^m(\tilde{\rho }),\hfill \end{array}$$

(A14)

and similarly for the remaining three equations, Equations (A2)-(A4). The first line corresponds to the differential equation of the Hankel function. Inserting Equation (A13) into Equation (A14), dividing by $H_m^{(1)}(\tilde{\rho })$ and using the asymptotic expansion of the Hankel function for large argument,

$${\displaystyle \frac{1}{H_m^{(1)}(\tilde{\rho })}}{\displaystyle \frac{d{H}_m^{(1)}(\tilde{\rho })}{d\tilde{\rho }}}\approx -{\displaystyle \frac{1}{2\tilde{\rho }}}+i.$$

(A15)

we obtain the following differential equations for the functions $g_j(\tilde{\rho })$:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d^2{g}_1}{d{\tilde{\rho }}^2}}& +& 2i{\displaystyle \frac{d{g}_1}{d\tilde{\rho }}}+\left[{\tilde{E}}^2+2\tilde{E}+{\displaystyle \frac{2m-1}{{\tilde{\rho }}^2}}\right]g_1\hfill \end{array}$$

$$\begin{array}{ccc}& =& (1+\tilde{E})\tilde{\mathsf{\Delta }}{g}_4+\tilde{\mathsf{\Delta }}{\displaystyle \frac{d{g}_3}{d\tilde{\rho }}}+{\displaystyle \frac{d\tilde{\mathsf{\Delta }}}{d\tilde{\rho }}}{g}_3+\left[{\displaystyle \frac{2m-1}{2\tilde{\rho }}}+i\right]\tilde{\mathsf{\Delta }}{g}_3,\hfill \\ \hfill {\displaystyle \frac{d^2{g}_2}{d{\tilde{\rho }}^2}}& +& 2i{\displaystyle \frac{d{g}_2}{d\tilde{\rho }}}+\left[{\tilde{E}}^2+2\tilde{E}\right]g_2=(1+\tilde{E})\tilde{\mathsf{\Delta }}{g}_3\hfill \end{array}$$

(A16)

$$\begin{array}{ccc}& -& \tilde{\mathsf{\Delta }}{\displaystyle \frac{d{g}_4}{d\tilde{\rho }}}+\left[{\displaystyle \frac{2m-1}{2\tilde{\rho }}}-i\right]\tilde{\mathsf{\Delta }}{g}_4-{\displaystyle \frac{d\tilde{\mathsf{\Delta }}}{d\tilde{\rho }}}{g}_4,\hfill \\ \hfill {\displaystyle \frac{d^2{g}_3}{d{\tilde{\rho }}^2}}& +& 2i{\displaystyle \frac{d{g}_3}{d\tilde{\rho }}}+\left[{\tilde{E}}^2-2\tilde{E}-{\displaystyle \frac{2m+1}{{\tilde{\rho }}^2}}\right]g_3\hfill \end{array}$$

(A17)

$$\begin{array}{ccc}& =& -(1-\tilde{E})\tilde{\mathsf{\Delta }}{g}_2+\tilde{\mathsf{\Delta }}{\displaystyle \frac{d{g}_1}{d\tilde{\rho }}}+{\displaystyle \frac{d\tilde{\mathsf{\Delta }}}{d\tilde{\rho }}}{g}_1-\left[{\displaystyle \frac{2m+1}{2\tilde{\rho }}}-i\right]\tilde{\mathsf{\Delta }}{g}_1,\hfill \\ \hfill {\displaystyle \frac{d^2{g}_4}{d{\tilde{\rho }}^2}}& +& 2i{\displaystyle \frac{d{g}_4}{d\tilde{\rho }}}+\left[{\tilde{E}}^2-2\tilde{E}\right]g_4=-(1-\tilde{E})\tilde{\mathsf{\Delta }}{g}_1\hfill \end{array}$$

(A18)

$$\begin{array}{ccc}& -& \tilde{\mathsf{\Delta }}{\displaystyle \frac{d{g}_2}{d\tilde{\rho }}}-\left[{\displaystyle \frac{2m+1}{2\tilde{\rho }}}+i\right]\tilde{\mathsf{\Delta }}{g}_2-{\displaystyle \frac{d\tilde{\mathsf{\Delta }}}{d\tilde{\rho }}}{g}_2.\hfill \end{array}$$

(A19)

These equations are solved iteratively. To zeroth order we keep the dominant terms,

$$\begin{array}{c}\hfill 2i{\displaystyle \frac{d}{d\tilde{\rho }}}\left[\begin{array}{c}{g}_1^{(0)}\\ g_2^{(0)}\\ g_3^{(0)}\\ g_4^{(0)}\end{array}\right]=\tilde{\mathsf{\Delta }}\left[\begin{array}{c}{g}_4^{(0)}\\ g_3^{(0)}\\ -g_2^{(0)}\\ -g_1^{(0)}\end{array}\right]+i\tilde{\mathsf{\Delta }}\left[\begin{array}{c}{g}_3^{(0)}\\ -g_4^{(0)}\\ g_1^{(0)}\\ -g_2^{(0)}\end{array}\right].\end{array}$$

(A20)

The remaining terms are first order terms and will be treated perturbatively. The solution of Equation (A20) is of the form

$$\begin{array}{ccc}& & g_1^{(0)}(\tilde{\rho })=C{e}^{-K(\tilde{\rho })},g_2^{(0)}(\tilde{\rho })=-iC{e}^{-K(\tilde{\rho })},\hfill \\ & & g_3^{(0)}(\tilde{\rho })=-C{e}^{-K(\tilde{\rho })},g_4^{(0)}(\tilde{\rho })=-iC{e}^{-K(\tilde{\rho })},\hfill \end{array}$$

(A21)

where $K(\tilde{\rho })$ has been defined previously in Equation (24) and C is a unitary constant $e^{i\gamma }$ to be determined later.

The equations for the first order terms in Equations (A16)-(A19) are

$$\begin{array}{ccc}& & 2i{\displaystyle \frac{d}{d\tilde{\rho }}}\left[\begin{array}{c}{g}_1^{(1)}\\ g_2^{(1)}\\ g_3^{(1)}\\ g_4^{(1)}\end{array}\right]-\tilde{\mathsf{\Delta }}\left[\begin{array}{c}{g}_4^{(1)}\\ g_3^{(1)}\\ -g_2^{(1)}\\ -g_1^{(1)}\end{array}\right]-i\tilde{\mathsf{\Delta }}\left[\begin{array}{c}{g}_3^{(1)}\\ -g_4^{(1)}\\ g_1^{(1)}\\ -g_2^{(1)}\end{array}\right]=-{\displaystyle \frac{d^2}{d{\tilde{\rho }}^2}}\left[\begin{array}{c}{g}_1^{(0)}\\ g_2^{(0)}\\ g_3^{(0)}\\ g_4^{(0)}\end{array}\right]\hfill \\ & & -\left[\begin{array}{c}({\tilde{E}}^2+2\tilde{E}+{\displaystyle \frac{2m-1}{{\tilde{\rho }}^2}})g_1^{(0)}\\ ({\tilde{E}}^2+2\tilde{E})g_2^{(0)}\\ ({\tilde{E}}^2-2\tilde{E}-{\displaystyle \frac{2m+1}{{\tilde{\rho }}^2}})g_3^{(0)}\\ ({\tilde{E}}^2-2\tilde{E})g_4^{(0)}\end{array}\right]+\tilde{\mathsf{\Delta }}\tilde{E}\left[\begin{array}{c}{g}_4^{(0)}\\ g_3^{(0)}\\ g_2^{(0)}\\ g_1^{(0)}\end{array}\right]+\tilde{\mathsf{\Delta }}{\displaystyle \frac{d}{d\tilde{\rho }}}\left[\begin{array}{c}{g}_3^{(0)}\\ -g_4^{(0)}\\ g_1^{(0)}\\ -g_2^{(0)}\end{array}\right]\hfill \\ & & +{\displaystyle \frac{d\tilde{\mathsf{\Delta }}}{d\tilde{\rho }}}\left[\begin{array}{c}{g}_3^{(0)}\\ -g_4^{(0)}\\ g_1^{(0)}\\ -g_2^{(0)}\end{array}\right]+{\displaystyle \frac{\tilde{\mathsf{\Delta }}}{2\tilde{\rho }}}\left[\begin{array}{c}(2m-1)g_3^{(0)}\\ (2m-1)g_4^{(0)}\\ -(2m+1)g_1^{(0)}\\ -(2m+1)g_2^{(0)}\end{array}\right].\hfill \end{array}$$

(A22)

Inserting the zero-order solutions, Equation (A21), into Equation (A22) and considering the Ansatz $g_1^{(1)}=C{a}_1{e}^{-K}$, $g_2^{(1)}=-iC{a}_2{e}^{-K}$, $g_3^{(1)}=-C{a}_3{e}^{-K}$ and $g_4^{(1)}=-iC{a}_4{e}^{-K}$, we obtain

$$\begin{array}{ccc}& & 2i{\displaystyle \frac{d}{d\tilde{\rho }}}\left[\begin{array}{c}{a}_1\\ -i{a}_2\\ -a_3\\ -i{a}_4\end{array}\right]-2i\tilde{\mathsf{\Delta }}\left[\begin{array}{c}{a}_1\\ -i{a}_2\\ -a_3\\ -i{a}_4\end{array}\right]-\tilde{\mathsf{\Delta }}\left[\begin{array}{c}-i{a}_4\\ -a_3\\ i{a}_2\\ -a_1\end{array}\right]-i\tilde{\mathsf{\Delta }}\left[\begin{array}{c}-a_3\\ i{a}_4\\ a_1\\ i{a}_2\end{array}\right]\hfill \\ & & =-\left[\begin{array}{c}({\tilde{E}}^2+2\tilde{E}+{\displaystyle \frac{2m-1}{{\tilde{\rho }}^2}})\\ -i({\tilde{E}}^2+2\tilde{E})\\ -({\tilde{E}}^2-2\tilde{E}-{\displaystyle \frac{2m+1}{{\tilde{\rho }}^2}})\\ -i({\tilde{E}}^2-2\tilde{E})\end{array}\right]+\tilde{\mathsf{\Delta }}\tilde{E}\left[\begin{array}{c}-i\\ -1\\ -i\\ 1\end{array}\right]+{\displaystyle \frac{\tilde{\mathsf{\Delta }}}{2\tilde{\rho }}}\left[\begin{array}{c}-(2m-1)\\ -i(2m-1)\\ -(2m+1)\\ i(2m+1)\end{array}\right],\hfill \end{array}$$

(A23)

where terms of the order ${\displaystyle \frac{d\tilde{\mathsf{\Delta }}}{d\tilde{\rho }}}$ and ${\tilde{\mathsf{\Delta }}}^2$ vanish identically. All terms are proportional to $e^{-K\tilde{\rho }}$, so that this factor cancels out. The differential equations then take the following form

$$\begin{array}{ccc}& & 2{\displaystyle \frac{d}{d\tilde{\rho }}}\left[\begin{array}{c}{a}_1\\ a_2\\ a_3\\ a_4\end{array}\right]+\tilde{\mathsf{\Delta }}\left[\begin{array}{c}(a_3+a_4)-2a_1\\ (a_3+a_4)-2a_2\\ (a_1+a_2)-2a_3\\ (a_1+a_2)-2a_4\end{array}\right]=-\tilde{\mathsf{\Delta }}\tilde{E}\left[\begin{array}{c}1\\ 1\\ -1\\ -1\end{array}\right]\hfill \\ & & +i\left[\begin{array}{c}{\tilde{E}}^2+2\tilde{E}+{\displaystyle \frac{2m-1}{{\tilde{\rho }}^2}}\\ {\tilde{E}}^2+2\tilde{E}\\ {\tilde{E}}^2-2\tilde{E}-{\displaystyle \frac{2m+1}{{\tilde{\rho }}^2}}\\ {\tilde{E}}^2-2\tilde{E}\end{array}\right]+i{\displaystyle \frac{\tilde{\mathsf{\Delta }}}{2\tilde{\rho }}}\left[\begin{array}{c}2m-1\\ -2m-1\\ -2m+1\\ 2m+1\end{array}\right],\hfill \end{array}$$

(A24)

which can be decoupled by taking linear combinations,

$$\begin{array}{ccc}& & 2{\displaystyle \frac{d}{d\tilde{\rho }}}(a_1-a_2)-2\tilde{\mathsf{\Delta }}(a_1-a_2)=i(2m-1)\left[{\displaystyle \frac{1}{{\tilde{\rho }}^2}}+{\displaystyle \frac{\tilde{\mathsf{\Delta }}}{\tilde{\rho }}}\right]\hfill \\ & & 2{\displaystyle \frac{d}{d\tilde{\rho }}}(a_3-a_4)-2\tilde{\mathsf{\Delta }}(a_3-a_4)=-i(2m+1)\left[{\displaystyle \frac{1}{{\tilde{\rho }}^2}}+{\displaystyle \frac{\tilde{\mathsf{\Delta }}}{\tilde{\rho }}}\right]\hfill \\ & & {\displaystyle \frac{d}{d\tilde{\rho }}}(a_1+a_2+a_3+a_4)=2i{\tilde{E}}^2-i{\displaystyle \frac{1}{{\tilde{\rho }}^2}}\hfill \\ & & {\displaystyle \frac{d}{d\tilde{\rho }}}(a_1+a_2-a_3-a_4)-2\tilde{\mathsf{\Delta }}(a_1+a_2-a_3-a_4)=4i\tilde{E}+2i{\displaystyle \frac{m}{{\tilde{\rho }}^2}}-2\tilde{\mathsf{\Delta }}\tilde{E}.\hfill \end{array}$$

(A25)

The above equations can be integrated

$$\begin{array}{ccc}\hfill a_1-a_2& =& -i(m-\frac{1}{2}){\int }_{\tilde{\rho }}^{\infty }dxexp[K(\tilde{\rho })-K(x)]\left[{\displaystyle \frac{1}{x^2}}+{\displaystyle \frac{\tilde{\mathsf{\Delta }}}{x}}\right]\hfill \end{array}$$

(A26)

$$\begin{array}{ccc}\hfill a_3-a_4& =& i(m+\frac{1}{2}){\int }_{\tilde{\rho }}^{\infty }dxexp[K(\tilde{\rho })-K(x)]\left[{\displaystyle \frac{1}{x^2}}+{\displaystyle \frac{\tilde{\mathsf{\Delta }}}{x}}\right]\hfill \end{array}$$

(A27)

$$\begin{array}{ccc}\hfill a_1+a_2& +& a_3+a_4=2i{\tilde{E}}^2\tilde{\rho }+{\displaystyle \frac{i}{\tilde{\rho }}}\hfill \\ \hfill a_1+a_2& -& a_3-a_4=-i{\int }_{\tilde{\rho }}^{\infty }dxexp[2K(\tilde{\rho })-2K(x)]\left[4\tilde{E}+{\displaystyle \frac{2m}{x^2}}\right]\hfill \end{array}$$

(A28)

$$\begin{array}{ccc}& -& 2{\int }_0^{\tilde{\rho }}dxexp[2K(\tilde{\rho })-2K(x)]\tilde{E}\tilde{\mathsf{\Delta }}.\hfill \end{array}$$

(A29)

We are only interested in leading and next-leading terms and will disregard higher order terms. Small parameters are $\tilde{\mathsf{\Delta }}(\tilde{\rho })$, $\tilde{E}$ and $\tilde{\rho }$, so that the term ${\tilde{E}}^2\tilde{\rho }$ of Equation (A28) and the last term of Equation (A29) can be neglected.

The remaining two integrals can be integrated by parts:

$${\int }_{\tilde{\rho }}^{\infty }dxexp[K(\tilde{\rho })-K(x)]\left({\displaystyle \frac{1}{x^2}}+{\displaystyle \frac{\tilde{\mathsf{\Delta }}}{x}}\right)={\displaystyle \frac{1}{\tilde{\rho }}},$$

(A30)

and

$$\begin{array}{ccc}& & {\int }_{\tilde{\rho }}^{\infty }dxexp[2K(\tilde{\rho })-2K(x)]\left(\tilde{E}+{\displaystyle \frac{m}{2x^2}}\right)=-\tilde{E}\tilde{\rho }+{\displaystyle \frac{m}{2\tilde{\rho }}}+G(\tilde{\rho })\hfill \end{array}$$

(A31)

$$\begin{array}{ccc}& & G(\tilde{\rho })={\int }_{\tilde{\rho }}^{\infty }dxexp[2K(\tilde{\rho })-2K(x)]2\tilde{\mathsf{\Delta }}(x)\left(\tilde{E}x-{\displaystyle \frac{m}{2x}}\right).\hfill \end{array}$$

(A32)

The function $G(\tilde{\rho })$ equals zero at the matching point (see Appendix 4) and defines the bound state energies.

It is now straightforward to solve Equations (A26)-(A29) for $a_j$

$$\begin{array}{ccc}& & a_1=i\tilde{E}\tilde{\rho }-i{\displaystyle \frac{2m-1}{2\tilde{\rho }}}-iG(\tilde{\rho })\hfill \\ & & a_2=i\tilde{E}\tilde{\rho }-iG(\tilde{\rho })\hfill \\ & & a_3=-i\tilde{E}\tilde{\rho }+i{\displaystyle \frac{2m+1}{2\tilde{\rho }}}+iG(\tilde{\rho })\hfill \\ & & a_4=-i\tilde{E}\tilde{\rho }+iG(\tilde{\rho }).\hfill \end{array}$$

(A33)

This corresponds to the leading and next-leading solution for the wave function for large argument.

Appendix A.4. Matching of Wave Functions

The next step consists in matching the solution for the wave function for small $\tilde{\rho }$ and large $\tilde{\rho }$ at an intermediate value ${\tilde{\rho }}_c$ from the core of the vortex. The wave function is determined uniquely if this matching is independent of ${\tilde{\rho }}_c$. This condition on ${\tilde{\rho }}_c$ also implies that $G({\tilde{\rho }}_c)\approx 0$ and consequently the energy for the $m^{th}$ bound state is given by

$$E_m=m{\int }_{{\tilde{\rho }}_c}^{\infty }dxexp[-2K(x)]{\displaystyle \frac{\tilde{\mathsf{\Delta }}}{x}}/{\int }_{{\tilde{\rho }}_c}^{\infty }dxexp[-2K(x)].$$

(A34)

The excited states are then equally spaced energy levels.

The solutions for $\tilde{\rho }>{\tilde{\rho }}_c$ are given by $f_j^m(\tilde{\rho })=B_j[H_m^{(1)}(\tilde{\rho })g_j(\tilde{\rho })+c.c.]$, where the coefficients $B_j$ are defined just after (Equation (A13)). The functions $g_j$ were calculated consistently to first order of perturbation and can be written as an exponential, i.e. $g_j(\tilde{\rho })\propto exp[-K(\tilde{\rho })+a_j(\tilde{\rho })]$. This expression remains correct to first order. For $\tilde{\rho }<{\tilde{\rho }}_c$ the solutions are Equations (A9)–(A12), which have to be matched at ${\tilde{\rho }}_c$ following a procedure similar to Refs. [12] and [23]. To match the wave functions we use the asymptotic expansions for the Bessel and Hankel functions. To simplify we only explicitly work out the matching for the function $f_1^m$; the remaining three functions follow similarly.

For $\tilde{\rho }<{\tilde{\rho }}_c$ we have

$$\begin{array}{ccc}\hfill f_1^m({\tilde{\rho }}_c)& & =A_1{J}_{m-1}\left((1+\tilde{E}){\tilde{\rho }}_c\right)exp(-K({\tilde{\rho }}_c))\hfill \\ & & =\sqrt{{\displaystyle \frac{2}{\pi {\tilde{\rho }}_c}}}cos[(1+\tilde{E}){\tilde{\rho }}_c-{\displaystyle \frac{\pi }{2}}(m-1)-{\displaystyle \frac{\pi }{4}}\hfill \\ & & +{\displaystyle \frac{{(m-1)}^2-\frac{1}{4}}{2(1+\tilde{E}){\tilde{\rho }}_c}}]exp(-K({\tilde{\rho }}_c)),\hfill \end{array}$$

(A35)

and similarly for $\tilde{\rho }>{\tilde{\rho }}_c$ we obtain (with $B_1={\displaystyle \frac{1}{2}}$)

$$\begin{array}{ccc}\hfill f_1^m({\tilde{\rho }}_c)& & ={\displaystyle \frac{1}{2}}\left[H_m^{(1)}({\tilde{\rho }}_c)g_1({\tilde{\rho }}_c)+c.c.\right]B_1\hfill \\ & & =\sqrt{{\displaystyle \frac{2}{\pi {\tilde{\rho }}_c}}}\{exp\left[i\left({\tilde{\rho }}_c-{\displaystyle \frac{\pi m}{2}}-{\displaystyle \frac{\pi }{4}}+{\displaystyle \frac{m^2-\frac{1}{4}}{2{\tilde{\rho }}_c}}\right)\right]\hfill \\ & & \times exp\left[i\left(\gamma +\tilde{E}{\tilde{\rho }}_c-{\displaystyle \frac{2m-1}{2{\tilde{\rho }}_c}}\right)\right]+c.c.\}{\displaystyle \frac{1}{2}}exp(-K({\tilde{\rho }}_c)).\hfill \end{array}$$

(A36)

The phase $\gamma$ arises from the constant C for the solution for $\tilde{\rho }>{\tilde{\rho }}_c$. For $\gamma ={\displaystyle \frac{\pi }{2}}$ the two expressions, (A35) and (A36), are identical to next-leading order for all four functions, $j=1,\cdots ,4$.

Appendix A.5. Wave Functions for Excited States

Approximate expressions to leading order for the amplitudes of the wave functions for the excited bound states, $m\ge 1$, are given by

$$\begin{array}{ccc}& & f_1^m(\mathbf{r})=C^{\prime }{J}_{m-1}(\tilde{\rho })e^{-K(\tilde{\rho })}{e}^{i(m-1)\nu },\hfill \\ & & f_2^m(\mathbf{r})=C^{\prime }{J}_m(\tilde{\rho })e^{-K(\tilde{\rho })}{e}^{im\nu },\hfill \\ & & f_3^m(\mathbf{r})=-C^{\prime }{J}_{m+1}(\tilde{\rho })e^{-K(\tilde{\rho })}{e}^{i(m+1)\nu },\hfill \\ & & f_4^m(\mathbf{r})=C^{\prime }{J}_m(\tilde{\rho })e^{-K(\tilde{\rho })}{e}^{im\nu },\hfill \end{array}$$

(A37)

and the energy wave function with energy $E_m$ is then

$$\begin{array}{ccc}\hfill {\widehat{\psi }}_E& =& C^{\prime }\int d^{2}r{e}^{-K(\tilde{\rho })}[J_{m-1}(\tilde{\rho })e^{-i\nu }{\psi }_{\uparrow }(\mathbf{r})+J_m(\tilde{\rho }){\psi }_{\downarrow }(\mathbf{r})\hfill \\ & & -J_{m+1}(\tilde{\rho })e^{i\nu }{\psi }_{\uparrow }^{†}(\mathbf{r})+J_m(\tilde{\rho })e^{i\nu }{\psi }_{\downarrow }^{†}(\mathbf{r})]e^{im\nu }.\hfill \end{array}$$

(A38)

While for $m=0$ the wave function corresponds to a Majorana state (${\mathsf{\Psi }}_M={\mathsf{\Psi }}_M^{†}$), for $m\ne 0$ it represents an ordinary fermion wave function (not self-adjoint). Note that for $m=0$, $J_{-1}(-\tilde{\rho })=-J_1(\tilde{\rho })$, so that there is an additional minus sign, making the difference between a Majorana and ordinary fermion.