Vortex bound states of charge and magnetic fluctuations-induced topological superconductors in heterostructures

TL;DR

This paper investigates vortex bound states in topological superconductor heterostructures induced by charge and magnetic fluctuations, revealing zero-energy states that could host Majorana modes without proximity effects.

Contribution

It introduces a minimal model for vortex states in fluctuation-induced topological superconductors, showing zero-energy bound states across different pairing symmetries.

Findings

Zero-energy vortex bound states exist for all pairing symmetries.

Majorana modes could be observed without proximity-induced superconductivity.

Vortex states are accessible at relatively high temperatures.

Abstract

The helical electron states on the surface of topological insulators or elemental Bismuth become unstable toward superconducting pairing formation when coupled to the charge or magnetic fluctuations. The latter gives rise to pairing instability in chiral channels $d_{xy}\pm i d_{x^2-y^2}$, as has been observed recently in epitaxial Bi/Ni bilayer system at relatively high temperature, while the former favors a pairing with zero total angular momentum. Motivated by this observation we study the vortex bound states in these superconducting states. We consider a minimal model describing the superconductivity in the presence of a vortex in the superconducting order parameter. We show that zero-energy states appear in the spectrum of the vortex core for all pairing symmetries. Our findings may facilitate the observation of Majorana modes bounded to the vortices in heterostructures with no need for a proximity-induced superconductivity and relatively large value of $\Delta/E_F$.

Equations

$$S_e=\int d\tau d\mathbf{r}{\Psi }^{†}\left({\partial }_{\tau }+{\tilde{H}}_0\right)\Psi .$$

$${\tilde{H}}_0=\int d\mathbf{r}{\Psi }^{†}(\mathbf{r})\left(v_F[\mathbf{\sigma }\times \mathbf{p}{]}_z-\mu )\Psi (\mathbf{r}\right).$$

$${\psi }_{\uparrow \mathbf{k}}=\frac{1}{\sqrt{2}}({\psi }_{\mathbf{k}+}+{\psi }_{\mathbf{k}-}),{\psi }_{\downarrow \mathbf{k}}=\frac{i{e}^{-i{\varphi }_{\mathbf{k}}}}{\sqrt{2}}({\psi }_{\mathbf{k}+}-{\psi }_{\mathbf{k}-}).$$

$${\psi }_{\uparrow }(\mathbf{r})=\frac{1}{\sqrt{2}}[{\psi }_{+}(\mathbf{r})+{\psi }_{-}(\mathbf{r})],$$

$${\psi }_{\downarrow }(\mathbf{r})=\frac{i}{\sqrt{2}}\sum _{\mathbf{k}}{e}^{-i{\varphi }_{\mathbf{k}}}({\psi }_{\mathbf{k}+}-{\psi }_{\mathbf{k}-})e^{i\mathbf{k}\cdot \mathbf{r}}.$$

$$e^{-i{\varphi }_{\mathbf{k}}}\approx \frac{k_x}{k_F}-i\frac{k_y}{k_F},$$

$${\psi }_{\downarrow }(\mathbf{r})=\frac{1}{\sqrt{2}{k}_F}({\partial }_x-i{\partial }_y)\left[{\psi }_{+}(\mathbf{r})-{\psi }_{-}(\mathbf{r})\right].$$

$${\psi }_{\uparrow }(\mathbf{r})\approx \frac{1}{\sqrt{2}}{\psi }_{+}(\mathbf{r}),{\psi }_{\downarrow }(\mathbf{r})\approx \frac{1}{\sqrt{2}{k}_F}({\partial }_x-i{\partial }_y){\psi }_{+}(\mathbf{r}).$$

$$H_0=-\int d\mathbf{r}{\psi }_{+}^{†}(\mathbf{r})\left(\frac{v_F}{k_F}{\nabla }^2+\mu \right){\psi }_{+}(\mathbf{r}),$$

$$S_M=\frac{{\rho }_s}{2}\int d\tau d\mathbf{r}\left[-i(\mathbf{l}\times {\partial }_{\tau }\mathbf{l})+\kappa (\nabla \mathbf{l}{)}^2\right],$$

$$S_{eM}=g\int d\tau d\mathbf{r}{\Psi }^{†}\left(\mathbf{l}\cdot \mathbf{\sigma }\right)\Psi ,$$

$$S_M=\frac{T}{2}\sum _q{D}^{-1}(\mathbf{q})b_q^{†}{b}_q,$$

$$S_{eM}=\frac{gT}{2}\sum _{k,q,\alpha \beta }\left(b_q{\psi }_{k\alpha }^{†}{\sigma }_{\alpha \beta }^z{\psi }_{k+q,\beta }+h.c.\right),$$

$$S_c=\frac{T}{2\mathcal{A}}\sum _{k,k^{\prime }}U(\mathbf{k},{\mathbf{k}}^{\prime })e^{-i({\varphi }_k-{\varphi }_{k^{\prime }})}{\psi }_k^{†}{\psi }_{-k}^{†}{\psi }_{-k^{\prime }}{\psi }_{k^{\prime }},$$

$$U(\mathbf{k},{\mathbf{k}}^{\prime })=\sum _{l=\text{even}}{U}_l{e}^{il({\varphi }_k-{\varphi }_{k^{\prime }})}.$$

$$H_{\Delta }=\sum _{\mathbf{k}}\Delta (|\mathbf{k}|)e^{i(l-1){\varphi }_k}{\psi }_k^{†}{\psi }_{-k}^{†}+h.c.$$

$$H_{\Delta }=\sum _{\mathbf{k}}\frac{\Delta (|\mathbf{k}|)}{k_F}(k_x+i{k}_y){\psi }_k^{†}{\psi }_{-k}^{†}+h.c.,$$

$$\begin{gathered} H=\int d\mathbf{r}[-{\psi }^{†}\left(\frac{v_F}{k_F}{\nabla }^2+\mu \right)\psi +i\frac{\Delta (r)}{2k_F}{\psi }^{†}\{e^{in\theta },{\partial }_x \\ +i{\partial }_y\}{\psi }^{†}+i\frac{\Delta (r)}{2k_F}\psi \{e^{in\theta },{\partial }_x-i{\partial }_y\}\psi ], \end{gathered}$$

$$[H_{MF},{\gamma }_i^{†}]=E_i{\gamma }_i^{†}.$$

$${\gamma }_i^{†}=\int d\mathbf{r}(u_i(\mathbf{r}){\psi }^{†}(\mathbf{r})+v_i(\mathbf{r})\psi (\mathbf{r})).$$

$$H_{BdG}=\left(\begin{array}{ccc}-\frac{v_F}{k_F}{\nabla }^2-\mu & & i\frac{\Delta (r)}{2k_F}\left\{e^{in\theta },{\partial }_x+i{\partial }_y\right\}\\ i\frac{\Delta (r)}{2k_F}\left\{e^{-in\theta },{\partial }_x-i{\partial }_y\right\}& & \frac{v_F}{k_F}{\nabla }^2+\mu \end{array}\right).$$

$$H_{BdG}=\left(\begin{array}{cc}-\frac{v_F}{k_F}\left[\frac{1}{r}{\partial }_r(r{\partial }_r)+\frac{1}{r^2}{\partial }_{\theta }^2\right]-\mu & i\frac{\Delta (r)}{k_F}{e}^{i{n}^{\prime }\theta }\left({\partial }_r+i\frac{1}{r}{\partial }_{\theta }-\frac{n}{2r}\right)\\ & \\ i\frac{\Delta (r)}{k_F}{e}^{-i{n}^{\prime }\theta }\left({\partial }_r-i\frac{1}{r}{\partial }_{\theta }-\frac{n}{2r}\right)& \frac{v_F}{k_F}\left[\frac{1}{r}{\partial }_r(r{\partial }_r)+\frac{1}{r^2}{\partial }_{\theta }^2\right]+\mu \end{array}\right),$$

$${\partial }_r^{2}u+\frac{1}{r}{\partial }_{r}u-\frac{m_{+}^2}{r^2}u+\frac{k_F\mu }{v_F}u-i\frac{\Delta (r)}{v_F}\left({\partial }_{r}v-\frac{2m-1}{2r}v\right)=-\frac{k_{F}E}{v_F}u,$$

$${\partial }_r^{2}v+\frac{1}{r}{\partial }_{r}v-\frac{m_{-}^2}{r^2}v+\frac{k_F\mu }{v_F}v+i\frac{\Delta (r)}{v_F}\left({\partial }_{r}u+\frac{2m+1}{2r}u\right)=\frac{k_{F}E}{v_F}v,$$

$$\tilde{E}=m\frac{{\int }_0^{\infty }(\frac{\tilde{\Delta }(x^{\prime })}{x^{\prime }}+\frac{n^{\prime }}{x^{\prime 2}})e^{-2\chi (x^{\prime })}d{x}^{\prime }}{{\int }_0^{\infty }{e}^{-2\chi (x^{\prime })}d{x}^{\prime }}$$

$$\begin{gathered} H_{BdG}=(h_0-\mu ){\tau }_3+i\frac{\Delta (r)}{2k_F^3}\left\{e^{in\theta },({\partial }_x+i{\partial }_y{)}^3\right\}{\tau }_{+} \\ +i\frac{\Delta (r)}{2k_F^3}\left\{e^{-in\theta },({\partial }_x-i{\partial }_y{)}^3\right\}{\tau }_{-}, \end{gathered}$$

$$H=-i\mathbf{v}\cdot \nabla {\tau }_3+\Delta (r)\cos ({\theta }^{\prime }){\tau }_1+\Delta (r)\sin ({\theta }^{\prime }){\tau }_2,$$

$$\begin{gathered} H=-iv{\partial }_x{\tau }_3+\Delta (r)\cos (n\theta +(3-n)\varphi ){\tau }_1 \\ +\Delta (r)\sin (n\theta +(3-n)\varphi ){\tau }_2 \end{gathered}$$

$$H=-iv{\partial }_x{\tau }_3+\Delta (r)\cos (n\theta ){\tau }_1+\Delta (r)\sin (n\theta ){\tau }_2.$$

$$E_m={\omega }_0\left(-n\pi +2\pi \left(m+\frac{1}{2}\right)\right),$$

$${\partial }_x^{2}u+\frac{1}{x}{\partial }_{x}u-\frac{m_{+}^2}{x^2}u+u-i\tilde{\Delta }(x)\left({\partial }_{x}v-\frac{2m-1}{2x}v\right)=-\tilde{E}u$$

Full text

I Introduction

For many years it was a common belief that quantum mechanical particles are either bosons or fermions whose wave functions take a plus or minus sign, respectively, upon the exchange of two identical particles. But in the past forty years it has been realized that this picture for point-like particles is correct only for a space with dimensions equal or greater than three. There are other possibilities in two-dimensional (2D) space Wilczek (1982a). In 2D systems, particles’ wave functions can in general acquire a complex phase upon turning one particle around another one, the so-called Abelian anyons Leinaas and Myrheim (1977); Wilczek (1982b). The first physical realization of Abelian anyons occurred in systems exhibiting the fractional quantum Hall effect Halperin (1984); Arovas et al. (1984). For the non-Abelian anyons, on the other hand, the multi-component wave functions live in a degenerate subspace Fröhlich (1988), and the interchange of particles amounts to a unitary evolution matrix within the degenerate subspace. The non-commutative structure of matrices promises a platform for fault-tolerant topological quantum computations Nayak et al. (2008).

Majorana fermions, a special class of non-Abelian anyons, discovered first by theoretical particle physicist Ettore Majorana Majorana (2008) in 1937, have the property that they are their own anti-particles. If $\gamma_{i}$ and $\gamma^{\dagger}_{i}$ are donated as the annihilation and creation operators for a Majorana fermion in a quantum state $|i\rangle$ then $\gamma_{i}=\gamma^{\dagger}_{i}$ and $\left\{\gamma_{i},\gamma_{j}\right\}=2\delta_{ij}$. Majorana fermions were not observed in elementary particle physics. But the developed concepts were traced in condensed matter physics years after the theoretical discovery Kopnin and Salomaa (1991); Moore and Read (1991); Read and Green (2000). In the seminal work by Reed and Green Read and Green (2000), it is shown that a boundary between a 2D topological $p+ip$ superconducto and a trivial one hosts a single Majorana mode, and in a vortex core a Majorana bound state (MBS) is formed. A prime example of a system hosting MBS was introduced by Kitaev Kitaev (2001). The model is a chain of spin-less fermions with superconducting $p$-wave order parameter. Under certain conditions, where the bulk of the system is topologically non-trivial, two unpaired MBSs appear at the ends of an open chain.

The discovery of topological insulators (TI) Hasan and Kane (2010); Hasan and Moore (2011) provided a boost to the realization of Majorana fermions in solid state systems. In a celebrated work by Fu and Kane in Ref. [Fu and Kane, 2008] it is shown that the MBSs appear in the vortex cores of a conventional $s$-wave superconductor proximitized to the surface of TIs arising from the spin-momentum locked structure of the surface states of the TI and the underlying Berry phase, which makes the $s$-wave pairing formally like a $p$-wave pairing upon projection onto the surface states. Viewed the surface states as massive Dirac electrons, the vortex bound states of the corresponding fermion-vortex problem have been studied Jackiw and Rossi (1981); Seradjeh (2008). The MBSs can also be realized in doped TIs Hosur et al. (2011); Chiu et al. (2012). Moreover, the surface of TIs can be replaced by more conventional semiconductors with strong Rashba spin-orbit coupling. The latter lifts the spin degeneracy and creates multiple spin-momentum locked Fermi surfaces. An external magnetic or Zeeman field is required to tune a quantum phase transition from an induced $s$-wave superconductor to a topological superconductor with Majorana fermions dispersing along the edges of 2D systems Sau et al. (2010a, b); Alicea (2010) or bounded to the end points of semiconductor nano-wires Lutchyn et al. (2010); Oreg et al. (2010).

At the heart of these settings for generating a topological phase with MBSs lies the existence of three conventional ingredients: a semiconductor quantum well with strong Rashba coupling, an $s$-wave superconductor, and a ferromagnetic insulator. The Zeeman field produced by the ferromagnetic layer should be strong enough to remove the extra Fermi surfaces, and simultaneously should be weak for the induced superconducting pairing amplitude to survive, a condition which severely restricts the choice of materials. Moreover, the induced superconductting gap is rather small to allow for resolving MBSs Wang et al. (2012). Therefore, it is highly demanding to look for heterostructures with less impeding ingredients.

The recently discovered superconductivity in epitaxial Bismuth/Nickel (Bi/Ni) bilayer heterostructure with relatively large transition temperature $T_{c}\approx 4.2$K may provide an example of an intrinsic topological superconductor with chiral $d_{xy}\pm id_{x^{2}-y^{2}}$ order parameter Gong et al. (2017). The nodeless structure of the proposed gap function is consistent with measurements of frequency-dependent optical conductivity in time-domain THz spectroscopy Chauhan et al. (2019). It is also argued that the superconductivity could result from the bulk alloys near the interface Chao (2019), but it turns out the surface superconductivity is consistent with thinckness dependent of transition temperature Gong et al. (2017). A schematic representation of the bilayer system in shown in Fig. 1. The main advantage of the latter system over the quantum well structures discussed above is that here the superconductivity is intrinsically driven by the ferromagnetic fluctuations circumventing the proximity to an extra $s$-wave superconductor layer. Now the important question is: what are the vortex bound states in an intrinsic $d_{xy}\pm id_{x^{2}-y^{2}}$ topological superconductor in Bi/Ni bilayer system? The aim of this paper is to theoretically answer this question. We show that the Berry phase effects change the angular momentum of the order parameter by one giving rise to odd parity with zero-energy states. Furthermore, in contrary to conventional superconductors, large ratios of superconducting gap to Fermi energy $\Delta/E_{F}\approx 10^{-2}-10^{-1}$ in Bi/Ni allow for MBS to remain well-separted from the low-lying excited states.

For non-chiral $d$-wave superconductor $d_{x^{2}-y^{2}}$ the existence of exteneded Franz and Tešanović (1998) and localized core states Knezevic and Radovic (1999); Kato and Maki (2001) have been discucussed. The analoge of discrete Caroli-de Gennes-Matricon core states Caroli et al. (1964) for spin-rotatinally symetric chiral $d$-wave superconductors has been studied Franz and Tešanović (1998); Rosenstein et al. (2013), but no zero modes are reported. The heterostrucure of topological insulator Bi2Se3 film on a nodal $d$-wave superconductor Bi2Sr2CaCu2O8+δ has been experimentally studied recently where a rather large induced and isoptropic superconducting gap was reported Wang et al. (2013). The isotropically nodeless gap was attributed to an emergent s-wave component on the surface of TI due to broken four-fold symmetry with possible MBSs in vortex cores Li et al. (2015). Note that our system is also distinct from the case of $d_{x^{2}-y^{2}}$ superconductor proximitized to an electron gas studied in Ref. [Ortiz et al., 2018], where an external Zeeman field is applied to the system while in our work the time-reversal symmetry is already broken by the nature of chiral pairing.

The paper is organized as follows. In Sec. II we derive the single-particle Hamiltonian describing the helical electron states near the Fermi surface, then in Sec. III we derive the pairing correlations driven by magnetic fluctuations. In Sec. IV we study the vortex bound states for various cases and the existence of zero-energy states, and finally Sec. V concludes.

II the real space projected non-interacting Hamiltonian

We begin with the electronic structure of the surface of Bi exposed to the vacuum as shown in Fig. 1. Most likely, the interface adjacent to Ni doesn’t contribute to superconductivity due to strong pair breaking effects, which is consistent with the thickness dependency of $T_{c}$ Gong et al. (2017, 2015). The strong Rashba coupling splits the surface electron states to spin-momentum locked states with the largest pocket centered around the center of surface of Brillouin zone. For simplicity we only consider this pocket and, hence, the Bi thin film in our setup in Fig. 1 can be replaced with a surface of TI as well. Therefore the helical electron states are described by the following action in Euclidean-time formalism:

[TABLE]

The non-interacting Hamiltonian $\tilde{H}_{0}$ reads as

[TABLE]

where $\Psi=(\psi_{\uparrow},\psi_{\downarrow})^{T}$ with $\psi_{\updownarrow}$ as electron annihilation operators, $v_{F}$ is the magnitude of spin-orbit interaction or equivalently the Fermi velocity of electrons at the surface of TI, $\mu$ is the chemical potential, $\bm{\sigma}$ is a vector of Pauli matrices, and $\mathbf{p}=-i\bm{\nabla}$. We set $\hbar=k_{B}=1$ throughout.

Since we are interested in electron states near the Fermi surface, we use helical eigenstates $|\mathbf{k},\pm\rangle=(1,\pm ie^{-i\phi_{k}})^{T}/\sqrt{2}$ in momentum space, where $\tilde{H}_{0}(\mathbf{k})|\mathbf{k},\pm\rangle=\pm\varepsilon_{\mathbf{k}}|\mathbf{k},\pm\rangle$ with $\varepsilon_{\mathbf{k}}=v_{F}|\mathbf{k}|$. We assume that the Fermi level crosses the “+” band and is located well above the node. We first write the annihilation operators in the spin basis in terms of annihilation operators in the helical “$\pm$” basis:

[TABLE]

Then in real space they become

[TABLE]

To perform the momentum sum in Eq. (5) we approximate the angular exponential term using

[TABLE]

which is justified so long as the electron states near the Fermi surface are involved in the physical processes of interest such as pairing instabilities. Here $k_{F}=\mu/v_{F}$ is the Fermi momentum. The field operator $\psi_{\downarrow}(\mathbf{r})$ then reads as

[TABLE]

By projection to the Fermi surface we obtain

[TABLE]

Rewritting Eq. (2) by use of these expressions, the projected form of the non-interacting Hamiltonian reads as follows:

[TABLE]

which is not dissimilar to the Hamiltonian of a 2D Fermi gas via the identification $v_{F}/k_{F}=1/2m$. Hereafter we drop the subindex in the field and write $\psi_{+}(\mathbf{r})\equiv\psi(\mathbf{r})$.

III Model of magnetic fluctuations and pairing hamiltonian

To make the structure of the paper self-contained, in this section we present the details of a minimal model describing the superconductivity in the bilayer structure shown in Fig. 1 which is largely based on the Ref. [Gong et al., 2017]. In the regime of interest the superconducting $T_{c}$ is much lower than the Curie temperature of ferromagnetic Ni layer, and hence the ferromagnet is deep inside the ordered phase. We assume that the in-plane moments are aligned along the $y$ direction (see Fig. 1) and the low-energy fluctuations of the magnetic moments, the spin waves, are described by the vector $\mathbf{l}(\tau,\mathbf{r})$ in which $\mathbf{l}\cdot\hat{y}=0$. The magnetic fluctuations and their coupling to electrons are described, respectively, by the following actionsKargarian et al. (2016):

[TABLE]

and

[TABLE]

where $\rho_{s}$ is the density of magnetic moments, $\kappa$ is the characteristic of spin waves, and $g$ is the strength of interaction between electron spins and magnetic moments. For simplicity and in the interest of formation of Cooper pairs with zero center-of-mass momentum we only consider the out-of-plane fluctuations denoted by $l_{z}\equiv b$.

By taking Fourier transform to momentum space, the Eqs. (10-11) become

[TABLE]

and

[TABLE]

where $T$ is the temperature, $q=(\mathbf{q},2n\pi T)$, $k=(\mathbf{k},(2n+1)\pi T)$ with $n$ as an integer, and $D(\mathbf{q})=1/(\kappa\rho_{s}|\mathbf{q}|^{2}+\zeta)$ is the magnon propagator with a small gap $\zeta$ due to anisotropyTserkovnyak and Loss (2012). Upon integrating out the bosonic field $b$ and subsequent projection of electron fields $\psi_{k\alpha}$ to the Fermi surface described by effective spinless fermion operators $\psi_{k}$ in Eq. (3), we obtain the following effective interaction between electrons in the Cooper channelGong et al. (2017)

[TABLE]

where $\mathcal{A}$ is the area of the system. Here $U(\mathbf{k},\mathbf{k}^{\prime})$ is an even function of its arguments and can be expanded in angular harmonics as

[TABLE]

Therefore the even angular momentum components of the interaction matrix contribute to the odd component of the condensate $f=\langle\psi_{-k}\psi_{k}\rangle$. This is a direct result of the non-trivial topology of Dirac Fermions. The effective angular momentum of the condensate $f$ is decreased by one due to the Berry phase. Inserting Eq. (15) in Eq. (14) and decoupling the interaction in the Cooper channels $f$, we obtain the mean-field BCS Hamiltonian as follows.

[TABLE]

In this work we only consider the channels with the lowest angular momenta $l=0,\pm 2$. The superconducting instability in channels $l=\pm 2$ is driven by the magnetic fluctuations being relevant to Bi/Ni bilayer system, while the instability with $l=0$ arises from phonons or charge fluctuationsGor’kov and Rashba (2001), i.e., $\sigma^{z}\rightarrow\mathbf{1}$ in Eq. (11). In our formalism below we study all cases.

IV spectrum of vortex bound states

The formulation and arguments presented in preceding sections provide a minimal superconducting Hamiltonian using Eqs. (9) and (16), i.e., $H=H_{0}+H_{\Delta}$. In the following subsections we first derive the corresponding Bogoliubov-de Gennes (BdG) equations for each channel $l$ and then study the spectrum of vortex bound states.

IV.1 The cases with $l=2$ and $l=0$

By inspection we see that for both cases the phases of the pairing in Eq.(16) are simply complex conjugate of each other, and thus, they can be treated within the same formalism. We present the details of calculations for $l=2$ and will briefly discuss the $l=0$ case at the end of this subsection. For the former case the pairing term $H_{\Delta}$ in Eq. (16) is written as

[TABLE]

where we use Eq. (6), assuming pairing occurs near the Fermi surface. To introduce a vortex in the order parameter, we assume that the space profile of pairing gap in the polar coordinate is $\Delta(\mathbf{r})=\Delta(r)e^{in\theta}$, where $r$ is measured from the center of vortex and $n$ denotes the winding of the vortex, the degree of vorticity. Thus, the full mean-field Hamiltonian of this system in real space can be recast as

[TABLE]

where $\{A,B\}=(AB+BA)/2$ is a symmetric operator. We define $\gamma_{i}^{\dagger}$ as a creation operator for the $i$-th eigenstate of the mean-field Hamiltonian satisfying

[TABLE]

The operator $\gamma^{\dagger}$ is used to diagonalize the Hamiltonian and hence can be written as a linear combinations of $\psi$’s:

[TABLE]

From now on we drop the index $i$ for simplicity. Using Eq. (18) and Eq. (20) in Eq. (19), we get a system of differential equations of the form $H_{BdG}\varphi(\mathbf{r})=E\varphi(\mathbf{r})$, where $\varphi(\mathbf{r})=(u(\mathbf{r}),v(\mathbf{r}))^{T}$ and $H_{BdG}$ is the BdG Hamiltonian

[TABLE]

Rewriting the differential operators in polar coordinates, Eq. (21) assumes the following form:

[TABLE]

where $n^{\prime}=n+1$. We use a pseudo-rotation operator defined by the unitary transformation $U(\theta)=e^{-i(m+\frac{n^{\prime}}{2}\tau^{z})\theta}$, where $\tau^{z}$ is the Pauli matrix acting in particle-hole space, to remove the phase dependency of the pairing gap Sau et al. (2010b). That is, we write the wave function $\varphi(\mathbf{r})$ in the form $\varphi(r,\theta)=e^{i(m+\frac{n^{\prime}}{2}\tau^{z})\theta}\varphi(r)$. The possible values for $m$ are determined by the single-valued condition of wave function implying that $m$ is an integer (half-integer) for even (odd) $n^{\prime}$. Using this transformation the eigenvalue problem turns into a set of differential equations for $v(r)$ and $u(r)$ as

[TABLE]

where $m_{\pm}=(2m\pm n^{\prime})/2$. We see that the equations are not symmetric under $n\rightarrow-n$. Note that the shift in $n^{\prime}$ relative to winding $n$ by one comes from the $\pi$-Berry phase of the electron states on the Fermi surface. The latter phase shifts the relative angular momentum of pairs by one Zeytinoglu et al. . Therefore the bound states of cores with opposite phase winding around the vortex would have different energy spectra. A set of equations similar to those quoted in Eqs.(23-24) is presented for a $p$-wave superconductorRead and Green (2000); Matsumoto and Heeb (2001), where the kinetic terms are usually in the long wavelength limit and it’s assumed that the chemical potential is negative in the vortex core (the strong coupling phase) and positive outside (the weak coupling phase) with $\Delta(r)$ as a constant. In our case however we assume that $\mu$ to be constant and take a space-varying order parameter like conventional superconductors.

Due to the $p$-wave structure of the Hamiltonian (22) the vortex core hosts a Majorana fermion. In Appendix A we employ an approach similar to Ref.[Sau et al., 2010b] and explicitly show that a zero-energy state exists in the vortex core. Furthermore, in order to get more insight into the spectrum of the bound states we use a long-wave approximationCaroli et al. (1964); Kopnin and Salomaa (1991) for the wave function (see Appendix B for details) and show that the spectrum reads as

[TABLE]

where $\chi(x)=\int_{0}^{x}\frac{\tilde{\Delta}(x^{\prime})}{2}\,dx^{\prime}$, and we used dimensionless parameters $x=k_{F}r$ and $\tilde{\Delta}=\Delta/\mu$ and $\tilde{E}=E/\mu$. One has to note that we are a little cavalier in using the semiclassical approach, since the value of $\Delta/E_{F}$ in our system, as we discuss more in Sec.V, is rather large compared to conventional superconductors. Within the approximations used it turns out that the vortices with $n^{\prime}\neq 0$ would have very large energy if $m\neq 0$ simultaneously. Let us consider a vortex with the lowest value of vorticity $n^{\prime}=0$ corresponding to $n=-1$ as shown schematically in Fig. 1. The condition $U(2\pi)=1$ implies that the $m$ has to be an integer number with $m=0$ corresponding to a zero-energy state.

As we mentioned at the beginning of this subsection the cases with $l=2$ and $l=0$ can be treated on equal footing, since the corresponding equations are the same. The latter case, $l=0$, yields $n^{\prime}=n-1$ and a vortex with lowest winding number will have $n=1$. Again the vortex can host a zero-energy state.

IV.2 The Case with $l=-2$

In this case the BdG equations become third order and an analytical solution for them is a formidable task if not impossible. To circumvent this problem, we use the semi-classical approximation used in the analysis of the Andreev bound states in superconductorsAndreev (1964) and closely follow Refs.[Stone, 1996] and [Volovik, 1999]. The BdG equation in this case is:

[TABLE]

where $h_{0}=-\frac{v_{F}}{k_{F}}\nabla^{2}$ is the kinetic energy. For solving BdG equations we use an ansatz for the wave function as $\Psi=\varphi(\mathbf{r})e^{i\mathbf{q}\cdot\mathbf{r}}$ with an approximation that the momentum $\mathbf{q}$ is restricted to the Fermi surface, i.e. $\mathbf{q}=k_{F}(\cos\phi,\sin\phi)$ known as the momentum of a quasi-particle in the Andreev approximation. Using $\Psi$ in Eq. (26), we obtain

[TABLE]

which acts only on $\varphi(\mathbf{r})$. Here we defined $\mathbf{v}=(2v_{F}/k_{F})\mathbf{q}$ and $\theta^{\prime}=n\theta+3\phi$. We rotate the coordinates such that the new $x$-axis becomes parallel to $\mathbf{q}$:

[TABLE]

Then the $\phi$ dependence in the Hamiltonian can be removed using the transformation $\varphi\rightarrow e^{i(3-n)\phi\tau_{3}/2}\varphi$:

[TABLE]

This is a quasi one-dimensional problem derived in Ref.[Stone, 1996] (see Eq. (3.10) in the latter reference with replacement $\theta\rightarrow-n\theta$). To proceed we define an impact parameter for quasi-particles as $b=r\sin\theta$ which measures the minimum distance of the quasi particle trajectory from the origin of the vortex core. For the pairing gap we use a profile as $\Delta(r)=\Delta\Theta(r-R)$, where $\Theta(x)$ is the usual step function and $R$ is the radius of the vortex. The latter is of order of $R\simeq v_{F}/\Delta$. With these assumptions the quasi one-dimensional model Eq.(29) can be solved to obtain the energy spectrum of the bound states. For small values of $bk_{F}\ll 1$, corresponding to trajectories passing near the origin, the spectrum reads as

[TABLE]

where $\omega_{0}=v_{F}/2R$ is the angular velocity of the superfluid Stone (1996). Now it is clearly seen that for $n=1$ the spectrum becomes $E_{m}=2\pi\omega_{0}m$ and a zero mode corresponds to $m=0$. Therefore the vortex bound states for the $l=-2$ case also contain a zero-energy mode.

V Conclusions

This work is mainly motivated by the efforts put forward in recent years to find Majorana bound states in the vortex core of superconducting states. We proposed the superconducting epitaxial Bi/Ni bilayer as a platform to create and manipulate the Majorana states. The system has an advantage over the heterostructues proposed in the literatures in that the chiral superconducting states are created intrinsically due to the magnetic fluctuations of the ferromagnetic layer circumventing the need for a proximixed superconducting layer. The heterostructure here can be replaced by other materials combinations, e. g., a thin film of topological insulator Bi2Te3 deposited on the magnetic insulator layer FeTe Manna et al. (2017); He et al. (2014), or superconducting states in oxide interfaces Scheurer and Schmalian (2015), making our proposal for creating and manipulating of zero modes experimentally feasible.

The chiral superconducting states in Bi/Ni bilayer are characterized by total angular momentum $l=\pm 2$ corresponding to $d_{xy}\pm d_{x^{2}-y^{2}}$, which break the time-reversal symmetry. We showed that the underlying strong spin-orbit coupling alter the bound state spectrum in the vortex core. In particular we demonstrated that a zero-energy state corresponding to Majorana bound state appears at the vortex core for both cases of the pairing wave functions. We also showed that the case with total angular momentum $l=0$, which intrinsically does not break the time-reversal symmetry and my be induced by charge fluctuations, can also support a zero-energy state. The set-up studied in our work, as shown in Fig. 1, should be contrasted with proposals in the literatures where the superconductivity is induced by proximity. Our finings may motivate the search for Majorana zero modes in vortices in heterostructures with no need for proximity to an extra superconducting layer.

Another peculiar aspect of vortex bound states in Bi/Ni is that the Majorana bound state remains well separated from the low-lying excited states due to a relatively large value of $\Delta/E_{F}$. The superconducting gap is estimated to be about $\Delta\approx 0.7$ meV Chauhan et al. (2019) and the Fermi energy $E_{F}$ is about 27 meV for hole pockets and 10 meV for electron pockets Hofmann (2006) yielding a ratio of about $10^{-2}-10^{-1}$. The ratio is by an order of magnitude larger than the corresponding values for conventional superconductors. A clear observation of discrete bound states in iron-based superconductor FeTe0.55Se0.45Chen et al. (2018) is reported due to large value of $\Delta/E_{F}$. The surface of latter compound was shown to be a topological superconductor Zhang et al. (2018) hosting a well-resolved Majorana bound state Wang et al. (2018). Therefore the same sort of well-resolved bound states and Majorana zero mode should be observed in epitaxial Bi/Ni bilayer.

In summary, our work offers the chiral superconductor in epitaxial Bi/Ni bilayer as a new platform to explore the vortex states with (i) robust Majorana bound state due to large ratio of $\Delta/E_{F}$, (ii) less complexity in the heterostructure, and (iii) relatively high transition temperature.

VI Acknowledgments

The authors would like to acknowledge the support from the Sharif University of Technology under Grant No. G690208.

Appendix A Explicit calculation of zero energy state

Using the dimensionless parameters $x=k_{F}r$ and $\tilde{\Delta}=\Delta/\mu$ and $\tilde{E}=E/\mu$, the Eqs.(23-24) can be rewritten as

[TABLE]

Here we consider (31) and (32) and follow Ref.[Sau et al., 2010b] for $\tilde{E}=0$. We assume that the vortex boundary is at $x=x_{0}$ and take $\tilde{\Delta(x)}$ to be zero for $x<x_{0}$ and a constant value for $x>x_{0}$. For the region inside the vortex the equations become a set of decoupled Bessel equations with the general analytical solution:

[TABLE]

Where $J_{m}(x)$ is the Bessel function of the first kind. A closed solution of equations for $x>x_{0}$ is not tractable. Instead we try to find the asymptotic solution of equations in the limit $x\gg 1$. In this limit (31) and (32) become:

[TABLE]

[TABLE]

Looking for a decaying solution of the form $\begin{pmatrix}u\\ v\end{pmatrix}=\begin{pmatrix}u_{0}\\ v_{0}\end{pmatrix}e^{\kappa x}$ one gets:

[TABLE]

Where $f_{\pm}(x)$ are decaying functions with the asymptotic form $f_{\pm}(x)\rightarrow e^{-\kappa_{\pm}x}$ and $\kappa_{\pm}=\frac{|\tilde{\Delta}|}{2}\pm\sqrt{(\frac{\tilde{\Delta}}{2})^{2}-1}$. We have to match the solutions for inside and outside of the vortex at the vortex boundary. The condition $\varphi(x_{0}^{-})=\varphi(x_{0}^{+})$ gives:

[TABLE]

[TABLE]

$\varphi^{\prime}(x_{0}^{-})=\varphi^{\prime}(x_{0}^{+})$ yields:

[TABLE]

[TABLE]

These equations alongside with the normalization condition $\int(u^{2}(r)+v^{2}(r))d\mathbf{r}=1$ should be satisfied in order to have a solution. In general, it is not possible to satisfy all of these conditions by adjusting only 4 unknowns ($A_{1}$-$A_{4}$) and the problem is over-constrained and a zero mode solution does not exist. Nevertheless, in the special case of $n^{\prime}=0$ where $m_{+}=m_{-}$ from (37) and (38) we have $A_{2}=-iA_{1}$ which makes (39) and (40) identical and therefore there are only 3 independent boundary conditions which alongside with the normalization condition assign a unique value to $A_{1}$-$A_{4}$. This is compatible with (25) in which a zero energy state exists only if $n^{\prime}=0$.

Appendix B Energy spectrum of bound states

Following Ref.[Kopnin and Salomaa, 1991] and Ref.[Caroli et al., 1964], we assume that the wave functions in (31) and (32) take the following form:

[TABLE]

where $H^{1}_{q}$ and $H^{2}_{q}$ are the Hankel functions of first and second kinds, respectively, and $f$ and $g$ are slowly varying functions. We insert the above ansatz in Eqs. (31-32) and neglect the second derivatives of $f_{\pm}$ and $g_{\pm}$. Using the asymptotic behavior of Hankel functions, we obtain the following differential equations governing $f_{\pm}$ and $g_{\pm}$:

[TABLE]

[TABLE]

The low-energy spectrum of bound states in the vortex lies within the superconducting bulk gap. Therefore a natural assumption is to assume $\tilde{E}<\tilde{\Delta}$ in the equations above, otherwise there would be no bound states at the vortex core. Physically the bound states result from the Andreev reflections of quasiparticles in the vortex core Stone (1996). We also assume that $x\gg 1$ which means that we are considering the long distance behavior of the system. Then by treating the expressions on the right-hand side as perturbations, we obtain the following expressions for $f$ and $g$ up to first order:

[TABLE]

where $\chi(x)=\int_{0}^{x}\frac{\tilde{\Delta}(x^{\prime})}{2}\,dx^{\prime}$. In order to avoid the singularity of Hankel functions at the origin we should take $A_{1}=A_{2}$ and the second terms should vanish as $x\rightarrow 0$. These boundary conditions eventually lead to the following expression for the energy spectrum of vortex bound states $\tilde{E}$ in (25).