Majorana Kramers pair in a nematic vortex

TL;DR

This paper demonstrates that a nematic vortex in an odd-parity topological superconductor can host Majorana Kramers pairs without breaking time-reversal symmetry, with potential realization in doped Bi2Se3.

Contribution

It introduces the concept of a TR-invariant nematic vortex hosting Majorana Kramers pairs and provides analytic solutions for zero modes in such vortices.

Findings

Nematic vortices can host Majorana Kramers pairs without magnetic fields.

Analytic solutions for zero modes in continuous nematic vortices are provided.

Discrete nematic vortices formed by domain meeting also support Majorana modes.

Abstract

A time-reversal (TR) invariant topological superconductor is characterized by a Kramers pair of Majorana zero-energy modes on boundaries and in a core of a TR invariant vortex. A vortex defect that preserves TR symmetry has remained primarily of theoretical interest, since typically a magnetic field, which explicitly breaks TR, needs to be applied to create vortices in superconductors. In this work, we show that an odd-parity topological superconductor with a nematic pairing order parameter can host a nematic vortex that preserves TR symmetry and binds a Majorana Kramers pair. Such a nematic superconductor could be realized in metal-doped Bi$_2$Se$_3$, as suggested by recent experiments. We provide an analytic solution for the zero modes in a continuous nematic vortex. In lattice, crystalline anisotropy can pin the two-component order parameter along high-symmetry directions. We show that a discrete nematic vortex, which forms when three nematic domains meet, also supports a TR pair of Majorana modes. Finally, we discuss possible experiments to probe the zero modes.

Figures

Equations

$$H$$

$${\Psi }_{\mathbf{r}}^{†}$$

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

$$\Delta (\mathbf{k})=\mathbf{\eta }\cdot \mathbf{F}(\mathbf{k}),$$

$$F_{\alpha }$$

$$F_{\alpha }^{(1)}$$

$$F_x^{(2)}$$

$$\Delta ({\mathcal{R}}_3\mathbf{\eta },{\mathcal{R}}_3{\mathbf{k}}_{\perp },k_z)={\mathcal{C}}_3\Delta (\mathbf{\eta },{\mathbf{k}}_{\perp },k_z){\mathcal{C}}_3^{†},$$

$$\Delta (\mathbf{\eta }={\eta }_x\hat{x},-k_x,k_y,k_z)$$

$$=$$

$$\Delta (\mathbf{k})=\mathbf{d}(\mathbf{k})\cdot \mathbf{\sigma }.$$

$$\mathbf{d}(\mathbf{k})={\eta }_x({\lambda }_2{k}_y+{\lambda }_3{k}_z,{\lambda }_2{k}_x,{\lambda }_1{k}_x).$$

$$\mathbf{k}=\pm (0,{\lambda }_3,-{\lambda }_2)k_F/\sqrt{{\lambda }_3^2+{\lambda }_2^2},$$

$$\mathbf{d}(\mathbf{k})={\eta }_y({\lambda }_2{k}_x,-{\lambda }_2{k}_y+{\lambda }_3{k}_z,{\lambda }_1{k}_y).$$

$$\mathbf{\eta }(\mathbf{r})=\eta (\rho )e^{in\varphi }(\cos m\varphi ,\sin m\varphi ),$$

$$\tilde{\mathcal{H}}$$

$$h_0$$

$$h_1$$

$$h_2$$

$$=-i\eta e^{i2\varphi }[-i({\partial }_{\rho }-\frac{1}{2\rho })+\frac{{\partial }_{\varphi }}{\rho }],$$

$$e^{i\ell \varphi }[e^{i\varphi }{f}_1(\rho ),e^{-i\varphi }{f}_2(\rho ),e^{i\varphi }{f}_3(\rho ),e^{-i\varphi }{f}_4(\rho ){]}^{\text{T}},$$

$${\gamma }_{\uparrow }^{†}$$

$${\gamma }_{\downarrow }^{†}$$

$$\displaystyle\big{[}\xi_{1},\xi_{2}\big{]}$$

$$\xi (\rho )$$

$$H_1=\int d\mathbf{r}(c_{\mathbf{r}\uparrow }^{†},c_{\mathbf{r}\downarrow }){\mathcal{H}}_1\left(\begin{array}{c}{c}_{\mathbf{r}\uparrow }\\ c_{\mathbf{r}\downarrow }^{†}\end{array}\right),{\mathcal{H}}_1=h_0{\tau }_z+{\lambda }_1{h}_1{\tau }_x,$$

$$h_0$$

$$h_1$$

$$\exp (iq\varphi )[f_1(\rho ),f_2(\rho ){]}^{\text{T}},$$

$$\frac{R}{\sin \theta }-L<\rho <\frac{R}{\sin \theta },$$

$$d_q^{†}=\int d\mathbf{r}{e}^{iq\varphi }{\psi }_q(\rho )[i{c}_{\mathbf{r}\uparrow }^{†}+c_{\mathbf{r}\downarrow }],$$

$${\tilde{H}}_2={\lambda }_2\sum _{p>0}{\epsilon }_p{d}_{1+p}^{†}{d}_{1-p}^{†}+{\epsilon }_p^{\ast }{d}_{1-p}{d}_{1+p},$$

$${\gamma }_{\uparrow }^{†}=d_1^{†}+d_1,{\gamma }_{\downarrow }^{†}=i(d_1^{†}-d_1),$$

$$\mathcal{F}=$$

$$+b_2({\eta }_x^2+{\eta }_y^2{)}^2+b_3[({\eta }_{+}^{\ast }{\eta }_{-}{)}^3+({\eta }_{-}^{\ast }{\eta }_{+}{)}^3].$$

$${\tilde{h}}_0{\xi }_1+{\lambda }_1{h}_1{\xi }_2^{\ast }-{\lambda }_2{\tilde{h}}_2{\xi }_1^{\ast }=0,$$

$${\tilde{h}}_0{\xi }_2^{\ast }-{\lambda }_1{h}_1{\xi }_1-{\lambda }_2{\tilde{h}}_2{\xi }_2=0,$$

Full text

I Introduction

Topological superconductors represent a paradigmatic system where topological effects and many-body interactions can have an interesting interplay. Topological superconductors with broken time reversal (TR) symmetry support Majorana zero-energy modes localized at the boundariesKitaev (2001) and in the vortex coresRead and Green (2000); Ivanov (2001). The non-Abelian statistics of Majorana modes may be utilized for fault-tolerant quantum computationKitaev (2003). The search for Majorana modes has been greatly influenced by the theoretical proposals based on topological insulator surface statesFu and Kane (2008) or Rashba wiresSau et al. (2010); Lutchyn et al. (2010) in proximity to ordinary $s$-wave superconductors in presence of magnetic field. Exciting experimental progressAlbrecht et al. (2016); Sun et al. (2016) has been made along the lines of these proposals. When TR symmetry is preserved, the Kramers pairs of Majorana modes have been predicted to exist at the boundaries of some families of topological superconductors Qi et al. (2009). TR symmetric topological superconductors have been theoretically proposed to be realized, for example, in CuxBi2Se3(Fu and Berg, 2010) and in materials with strong spin-orbit coupling in proximity to $\pi$ Josephson junctionKeselman et al. (2013); Schrade et al. (2015) or $s_{\pm}$-wave superconductorsZhang et al. (2013a). Among the expected exotic properties of Majorana Kramers pairs are non-Abelian statisticsLiu et al. (2014) and topological two-channel Kondo effectBao and Zhang (2016).

The Kramers pair of Majorana modes has also been predicted to appear in TR-invariant vortex core of a TR-invariant topological superconductorQi et al. (2009). TR transformation changes the parity of fermion number associated with the zero modes, and therefore acts as a supersymmetry in the vortexQi et al. (2009). However, it remains unclear how a vortex defect that preserves TR symmetry can be physically realized, because a vortex in a superconductor is typically induced by a magnetic field that explicitly breaks TR symmetry. Here we propose that a TR-symmetric vortex naturally appears in a topological superconductor with nematic order. The nematic superconductor has a two-component order parameterFu (2014), which can be viewed as a planar vector. We demonstrate analytically and numerically that a nematic vortex, such that only the vector part of the order parameter rotates, without $U(1)$ phase winding, is TR symmetric and binds a Majorana Kramers pair.

Our work is motivated by the recent experimental studies of superconducting $M$-doped Bi2Se3, where the metal element $M$ can be Cu, Nb and Sr. Shortly after the experimental discovery of superconductivity in CuxBi2Se3Hor et al. (2010), Fu and Berg first proposed that a TR-invariant topological superconductor could be realized in this systemFu and Berg (2010). Majorana surface states associated with the topological superconductor were expected to be observable by point contact measurementSasaki et al. (2011); however, scanning tunneling microscopy measurementsLevy et al. (2013) performed so far have not detected them. Despite the lack of observed surface states, measurements of bulk properties have revealed compelling evidence that the superconducting state in $M_{x}$Bi2Se3 spontaneously breaks the lattice discrete rotational symmetry and therefore is nematic, and potentially topological. The main experimental evidence comes from the two-fold anisotropy of various bulk responses to an external magnetic field, including Knight shiftMatano et al. (2016), specific heatYonezawa et al. (2017), upper critical fieldYonezawa et al. (2017); Pan et al. (2016) and magnetic torqueAsaba et al. (2017). Corresponding theories of these bulk properties have also been developedHashimoto et al. (2013); Nagai and Ota (2016); Venderbos et al. (2016a); Kozii et al. (2016).

Our paper is organized as follows. In Sec. II, we present the effective HamiltonianVenderbos et al. (2016b) for the nematic odd-parity superconductor, and review its symmetry properties and gap structure. In Sec. III, we study electronic states in the continuous nematic vortex and provide an analytic solution to the Majorana Kramers pair. The zero modes can be understood in two different limits. In the limit of only helical $p$-wave pairing, spin up (down) electrons form $p_{x}+ip_{y}$ ($p_{x}-ip_{y}$) Cooper pair condensates. A nematic vortex thus exactly realizes a TR-symmetric vortex discussed in Ref. Qi et al., 2009, where spin up and down electrons respectively experience a vortex and antivortex in the pairing order parameters. In the other limit of only polar $p$-wave pairing, the bulk of the superconductor is gapless. Surprisingly, we find that a single vortex can induce a flat band of zero energy modes. A finite helical $p$-wave pairing potential lifts all the zero modes to finite energy except for one pair of Majorana modes. In Sec. IV, we study a discrete nematic vortex, which can be experimentally realized at the core where three nematic domains meet. In agreement with the analysis of continuous vortex, the discrete nematic vortex also binds a Kramers pair of Majorana modes. When the out-of-plane dispersion is included, the Majorana modes evolve into helical modes. Finally in Sec. V, we discuss possible ways to experimentally probe the zero modes in nematic vortex.

II Hamiltonian and symmetry

Undoped Bi2Se3 is a strong topological insulator with both TR symmetry and inversion symmetry. Electronic bands in Bi2Se3 are therefore doubly degenerate at each $k$ point. For $M_{x}$Bi2Se3, we assume that the chemical potential lies in the conduction bands. When only the conduction bands near $\Gamma$ point are retained in a low-energy theory, the superconducting state in $M_{x}$Bi2Se3 can be described by a one-orbital modelFu and Berg (2010); Yip (2013); Venderbos et al. (2016b):

[TABLE]

Here $\boldsymbol{k}$ is the three-dimensional momentum operator, and $\uparrow$ and $\downarrow$ are the pseudo-spin label of the two degenerate states in the conduction band. $\tau_{\pm}=(\tau_{x}\pm i\tau_{y})/2$, and $\tau_{x,y,z}$ are Pauli matrices in Nambu space. $\mathcal{H}_{0}$ is the kinetic energy arising from band structure. We assume a parabolic dispersion: $\mathcal{H}_{0}(\boldsymbol{k})=\beta\boldsymbol{k}^{2}-\mu$, where $\mu>0$ is the chemical potential. The nematic superconductor has an odd-parity pairing in the two-component $E_{u}$ representation of the $D_{3d}$ point group. The pairing potential can be expressed as:

[TABLE]

where both $\boldsymbol{\eta}$ and $\boldsymbol{F}$ represent a two-component vector. $\boldsymbol{\eta}=(\eta_{x},\eta_{y})$ acts as the nematic order parameter. $\boldsymbol{F}$ can be decomposed into three parts when only linear $\boldsymbol{k}$ terms are retainedVenderbos et al. (2016b):

[TABLE]

where $\alpha=x,y$ and $F_{x,y}$ are the two components of $\boldsymbol{F}$. $\sigma_{x,y,z}$ are Pauli matrices in spin space. $\lambda_{1,2,3}$ are real parameters that characterize the relative strength of different pairing terms. When $\boldsymbol{\eta}$ is real, $\boldsymbol{F}^{(1)}$ and $\boldsymbol{F}^{(3)}$ represent polar $p$-wave pairing, and $\boldsymbol{F}^{(2)}$ is helical $p$-wave pairing where spin up (down) electrons respectively form $p_{x}+ip_{y}$ ($p_{x}-ip_{y}$) Cooper pairs.

We now elaborate on the symmetry properties when $\boldsymbol{\eta}$ is spatially uniform. The pairing potential has odd parity: $\Delta(-\boldsymbol{k})=-\Delta(\boldsymbol{k})$. TR symmetry is preserved if $\boldsymbol{\eta}$ is real: $T\mathcal{H}(\boldsymbol{k})T^{-1}=\mathcal{H}(-\boldsymbol{k})$, where $T=-i\sigma_{y}K$ is the TR operator and $K$ is the complex conjugate. The particle-hole (PH) symmetry can be expressed as: $\Xi\mathcal{H}(\boldsymbol{k})\Xi^{-1}=-\mathcal{H}(-\boldsymbol{k})$, where $\Xi=\sigma_{y}\tau_{y}K$.

The pairing potential $\Delta$ spontaneously breaks the three-fold rotational symmetry of Bi2Se3. Rotating $\boldsymbol{\eta}$ by $2\pi/3$ leads to a different but energetically degenerate pairing potential. This is manifested by the following identity:

[TABLE]

where $\boldsymbol{k}_{\perp}=(k_{x},k_{y})$ is the in-plane momentum. $\mathcal{R}_{3}$ is the real-space rotation matrix for a $2\pi/3$ rotation around $z$ axis, and $\mathcal{C}_{3}=\exp(-i\pi\sigma_{z}/3)$ is the rotation matrix in spin space.

In the $D_{3d}$ point group there are three mirror planes that are located at $x=0$ or its three-fold rotation counterparts. A mirror symmetry ($x\leftrightarrow-x$) is preserved if $\boldsymbol{\eta}$ is perpendicular to the corresponding mirror plane $x=0$:

[TABLE]

where $\mathcal{M}=\exp(-i\pi\sigma_{x}/2)$ is the mirror operator in spin space. The superconductor has point nodes protected by the mirror symmetry when $\boldsymbol{\eta}$ is along $\hat{x}$Zhang et al. (2013b); Yang et al. (2014); Fu (2014), as we discuss below.

As an odd-parity pairing potential, $\Delta$ can be expressed in terms of a vector $\boldsymbol{d}$:

[TABLE]

When $\boldsymbol{\eta}$ is real, the superconductivity gap is given by $|\boldsymbol{d}(\boldsymbol{k})|$ on the Fermi surfaceSigrist and Ueda (1991). We present the explicit form of $\boldsymbol{d}(\boldsymbol{k})$ for two representative cases. In the case of $\boldsymbol{\eta}=\eta_{x}\hat{x}$,

[TABLE]

Because of the mirror symmetry in (5), $d_{y,z}$ vanishes on the mirror plane $k_{x}=0$. Since $\boldsymbol{d}$ is odd under inversion, $|\boldsymbol{d}(\boldsymbol{k})|$ must have nodes on the mirror plane. The point nodes are located at:

[TABLE]

where $k_{F}$ is the Fermi wave vector $\sqrt{\mu/\beta}$. Note that the point nodes appear at non-zero $k_{z}$.

In another case $\boldsymbol{\eta}=\eta_{y}\hat{y}$,

[TABLE]

It is then impossible to make all three components of $\boldsymbol{d}(\boldsymbol{k})$ to vanish simultaneously on the Fermi surface, so the superconductor is fully gapped. The gap structures for the two different cases are schematically shown in Fig. 1. In summary, Hamiltonian (1) describes a topological Dirac superconductorYang et al. (2014) with surface Majorana arcs when $\boldsymbol{\eta}$ is perpendicular to one of the mirror planes, and otherwise a fully gapped odd-parity topological superconductor(Fu and Berg, 2010).

III Continuous nematic vortex

Since the nematic superconductor has a two-component order parameter, a continuous vortex defect along $\hat{z}$ is characterized by two winding numbers $(n,m)$:

[TABLE]

where $(\rho,\phi)$ are the polar coordinates of the in-plane position $\boldsymbol{r}_{\perp}$. The usual superconducting vortex has an integer phase winding number $n$ and a zero value of $m$. $n$ and $m$ can also take half integer values: $(n,m)=(\pm 1/2,\pm 1/2)$, which correspond to half vorticesIvanov (2001). Vortices with finite $n$ break TR symmetry. In this paper, we focus on a TR-invariant vortex with $(n,m)=(0,1)$, which we dub nematic vortex.

We study electronic states induced by the nematic vortex. As $k_{z}$ remains a good quantum number, we start by considering $k_{z}=0$. The dispersion of vortex states as a function of $k_{z}$ will be discussed in Sec. IV. The Hamiltonian (1) at $k_{z}=0$ and in the presence of nematic vortex can be expressed in polar coordinates as:

[TABLE]

where $\eta=|\boldsymbol{\eta}(\boldsymbol{r})|$, $k_{\pm}=k_{x}\pm ik_{y}$ and $\eta_{\pm}=\eta_{x}\pm i\eta_{y}$. In (11), the product between $\boldsymbol{\eta}$ and $\boldsymbol{k}$ is symmetrized as they no longer commute when $\boldsymbol{\eta}$ has spatial variation. In the nematic vortex, $\eta_{+}$ and $\eta_{-}$ respectively has a vortex $\exp(i\phi)$ and antivortex $\exp(-i\phi)$ structure. For the TR invariant vortex that we study, $\eta_{+}$ and $\eta_{-}$ vanish at the same location in the vortex core. In the following we will neglect the spatial variation in $|\boldsymbol{\eta}(\boldsymbol{r})|$ to simplify calculation.

$\tilde{\mathcal{H}}$ preserves both TR and PH symmetries. Therefore, the nematic vortex described by (11) is in class DIII and can be characterized by a $Z_{2}$ invariant that signals the presence or absence of Majorana Kramers doublet in the vortex coreTeo and Kane (2010). Here the existence of the zero modes is hinted by the limit in which $h_{1}$ term is absent in (11). In this limit, Hamiltonian (11) effectively describes a TR symmetric superconductor in which spin up and down electrons respectively form $p_{x}+ip_{y}$ and $p_{x}-ip_{y}$ Cooper pairs, and the superconducting phase has a vortex (antivortex) defect in the spin up (down) sector. Therefore, there must be a Kramers pair of Majorana modes when $\lambda_{1}=0$. By explicitly solving the eigenenergy problem of (11), we show that the zero modes are robust regardless of the exact values of $\lambda_{1,2}$.

Hamiltonian (11) has eigenvectors of the form:

[TABLE]

where $\ell$ takes integer values. For numerical diagonalization, $f_{i}(\rho)$ is expanded in terms of orthonormal Bessel functions on a finite-radius diskGygi and Schlüter (1991); Kawakami and Hu (2015). The energy spectrum as a function of $\ell$ is illustrated in Fig. 2(b). The numerical results indicate that a pair of zero-energy modes with $\ell=0$ are localized in the vortex core. Because of PH and TR symmetries, they must be Majorana modes and TR partners:

[TABLE]

which satisfy the Majorana condition ($\gamma_{\uparrow,\downarrow}^{\dagger}=\gamma_{\uparrow,\downarrow}$) and form Kramers doublet under TR transformation ($\gamma_{\uparrow}\rightarrow\gamma_{\downarrow}$, $\gamma_{\downarrow}\rightarrow-\gamma_{\uparrow}$). As shown in Appendix A,wave functions $\xi_{1,2}$ can be obtained analytically:

[TABLE]

where $\mathcal{N}$ represents a normalization factor and $J_{1}$ is Bessel function of the first order. In (14), $\lambda_{1,2}\geq 0$ is assumed, and $\tilde{\eta}$ is the effective total pairing strength $\eta\sqrt{\lambda_{1}^{2}+\lambda_{2}^{2}}$.

We analyze (14) in two different limits. In the limit $\lambda_{1}=0$, $\xi_{2}$ vanishes and $\gamma_{\uparrow,\downarrow}^{\dagger}$ is the Majorana mode in corresponding spin sector, as expected. In the other limit $\lambda_{2}=0$, the zero modes are also exponentially confined to the vortex center. This result is surprising, since the Hamiltonian (1) with zero $\lambda_{2}$ describes a gapless polar $p$-wave superconductor. Next we reveal another unexpected effect – that the polar $p$-wave superconductor supports not just one pair, but a flat band of zero modes localized around the nematic vortex core, despite the bulk being gapless.

Hamiltonian (1) with $k_{z}=0$ and $\lambda_{2}=0$ reduces to

[TABLE]

where $h_{0,1}$ is provided in (11). We reiterate their expressions below but introduce a new parameter $\theta$ that can be used to tune the system geometry,

[TABLE]

$\theta=\pi/2$ corresponds to a flat plane geometry. When $\theta<\pi/2$, Eq. (16) generalizes the problem on a plane to a cone with an opening angle $2\theta$, as illustrated in Fig. 3. This generalization is derived in Appendix B and is used to shed light on relationship between the vortex core states and edge modes in cylindrical geometry.

Similar to $\tilde{\mathcal{H}}$ in (11), the operator $\mathcal{H}_{1}$ in (16) has eigenstates of the form

[TABLE]

where $q$ is an integer. We solve (15) and (16) in the following range of $\rho$,

[TABLE]

with open boundary conditions. This specifies an annulus in the plane when $\theta=\pi/2$ and a truncated cone with a base radius $R$ when $\theta<\pi/2$.

An important special case is $\theta\to 0$. In this limit, $\rho\to\infty$, while $\rho\sin\theta$ stays at a finite value $R$. Geometrically, the truncated cone becomes a cylinder with a radius $R$ and a finite length $L$. Eq. (16) then describes exactly a uniform polar $p$-wave superconductor on the cylinder, as depicted in Fig. 3(c). It is known that a polar $p$-wave superconductor supports a flat band of zero-energy modes confined to edgesMartin and Morpurgo (2012).

Equations (16) and (18) state that the polar $p$-wave superconductor with a nematic vortex in the annular geometry can be related to a uniform one on the cylinder via a single parameter $\theta$. We solved the eigenvalue problem for different $\theta$ values using the ansatz (17). The energy spectrum as a function of $q$ is shown in Fig. 3. The zero-energy modes with small $|q|$ index on cylinder ($\theta=0$) remain at zero as $\theta$ increases, all the way to $\theta=\pi/2$.

Because of a PH symmetry of $\mathcal{H}_{1}$ ($\tau_{y}\mathcal{H}_{1}=-\mathcal{H}_{1}\tau_{y}$), the zero modes are eigenstates of $\tau_{y}$. Therefore the $q$th zero mode localized near the inner edge can be parametrized as:

[TABLE]

where $\psi_{q}(\rho)$ is the radial part of the wave function. The degree of confinement of $\psi_{q}(\rho)$ to the inner edge depends on the index $q$, as illustrated in Fig. 3(d). $\psi_{q}(\rho)$ moves away from the inner edge as $|q|$ increases due to a “centrifugal force” in (16). In the annulus geometry ($\theta=\pi/2$), the zero modes are robust when the inner radius becomes zero ($L\rightarrow R$). As shown in Appendix C, an analytic solution of the wave function $\psi_{q}(\rho)$ can be obtained at every integer $q$ for a disk with an infinite radius. The number of zero modes becomes finite on a finite-radius disk due to hybridization between modes deriving from the core and the outer edge. The stability of the flat band with respect to disorder and interactions is an interesting topic, beyond the scope of this paper. Previous work on related problems can be found, for example, in Ref Li et al., 2013; Potter and Lee, 2014; Baum et al., 2015.

Starting from the polar phase, we can examine the effect of the helical part of the pairing potential $\lambda_{2}(\boldsymbol{\eta}\cdot\boldsymbol{F}^{(2)}+\boldsymbol{F}^{(2)}\cdot\boldsymbol{\eta})/2$ by including it as a perturbation, and projecting it onto the zero modes $d_{q}^{\dagger}$. The projection leads to the following effective low energy Hamiltonian:

[TABLE]

where $p=0$ term is absent because of fermion anticommutation relation. Therefore, all zero modes except $d_{1}^{\dagger}$ are lifted to finite energy when $\lambda_{2}\neq 0$. $d_{1}^{\dagger}$ can be decomposed into two Majorana modes:

[TABLE]

which are exactly those presented in (13). As a summary, the Majorana Kramers doublet $\gamma_{\uparrow,\downarrow}^{\dagger}$ can be understood either from the helical limit $\lambda_{1}=0$ or the polar limit $\lambda_{2}=0$.

The degeneracy between $\gamma_{\uparrow}^{\dagger}$ and $\gamma_{\downarrow}^{\dagger}$ is protected by TR symmetry. TR changes the fermion parity associated with fermion mode $d_{1}$ and acts as: $d_{1}^{\dagger}\rightarrow-id_{1}$ and $d_{1}\rightarrow id_{1}^{\dagger}$, which is a manifestation of supersymmetry. A mass term $d_{1}^{\dagger}d_{1}$ is prohibited by TR symmetry. When TR symmetry is explicitly broken, for example by a Zeeman term of the form $\epsilon_{\text{Z}}\sigma_{z}$, generally there can be coupling between $\gamma_{\uparrow}^{\dagger}$ and $\gamma_{\downarrow}^{\dagger}$ that lifts the degeneracy.

IV Discrete nematic vortex

Because Bi2Se3 only has a discrete three-fold rotational symmetry, the nematic vector $\boldsymbol{\eta}$ is expected to be pinned to one of the high symmetry directions. The Ginzburg-Landau free energy for a uniform state isFu (2014):

[TABLE]

Here the sign of $b_{2}$ distinguishes between nematic and chiral states, which have respectively real and complex order parameters $\boldsymbol{\eta}$. We assume $b_{2}<0$ so that a nematic state is energetically favorable. The sixth-order terms in $\mathcal{F}$ capture the crystalline anisotropy. The sign of $b_{3}$ determines whether $\boldsymbol{\eta}$ is pinned to be perpendicular or parallel to one of the three mirror planes, for example, along $x$ or $y$ direction. Because of the lattice anisotropy, a continuous nematic vortex described by (10) is energetically unfavorable. Nevertheless, a discrete vortex can be realized around the core where three nematic domains meet, as illustrated in Fig. 4. Indeed, a recent experiment on the specific heat of CuxBi2Se3 under magnetic field has indicated the existence of multiple nematic domain walls within one sampleYonezawa et al. (2017). Therefore, a discrete nematic vortex may naturally appear due to local disorder or strain effects.

To numerically study the discrete vortex, we regularize the Hamiltonian (1) on a triangular lattice. Two representative types of discrete vortex are studied. In type I, $\boldsymbol{\eta}$ within each domain is parallel to a mirror plane. The bulk of each domain is therefore fully gapped, as illustrated in Fig. 1(b). The numerical calculation for $k_{z}=0$ indicates that there are four zero modes. Two of them are localized around the core center and the other two are bound to edges[see Fig. 4(a)], which are expected from the analysis of the continuous vortex. At finite $k_{z}$, the two zero modes localized around the core evolve into two branches that disperse linearly with $k_{z}$ as shown in Fig. 4(b) and therefore can be identified as helical Majorana modes.

In type II, $\boldsymbol{\eta}$ within each domain is perpendicular to a mirror plane. The bulk of each domain in this case is gapless, and the point nodes are at $k_{z}\neq 0$, as shown in Fig. 1(a). There are still two zero modes localized around the core for $k_{z}=0$, at which the bulk is fully gapped. At finite $k_{z}$, the helical modes merge into the bulk excitation continuum. The bulk excitations of the discrete vortex structure become gapless at a finite $k_{z}$ (Fig. 4(d)). The exact values of $k_{z}$ where the bulk excitation minimum are located are slightly different for the domain structure and the uniform state due to the presence of domain core and finite size effects. As $k_{z}$ may not be conserved in the presence of disorder, the zero modes in type II vortex are not as robust as those in type I.

The energy cost of a domain wall scales linearly with the area of the interface between two domains. This energy cost could be remedied by studying a sample with mesoscopic geometriesChung et al. (2007). Another possibility that discrete vortices could be stabilized is due to Imry-Ma scenarioImry and Ma (1975), as the nematic order parameter can couple to local disorder or strain fieldFu (2014).

V Discussion

The presence of the Majorana Kramers pair will lead to an enhanced density of states (DOS) near the vortex core at zero energy. To probe the enhanced DOS, it would be favorable to study materials in a thin film geometry so that the motion along $z$ is quantized. Scanning tunneling microscope (STM) is particularly suitable to search for the zero modes. STM has been successfully used to detect Majorana modes in different schemesNadj-Perge et al. (2014); Sun et al. (2016), and to identify nematic states in the quantum Hall regimeFeldman et al. (2016). Within each nematic domain of $M_{x}$Bi2Se3, the nematic order parameter will lead to a two-fold anisotropy in electronic band structure, which should be detectable in quasiparticle interference pattern in STMHofmann et al. (2013); Farrell et al. (2015); Lambert et al. (2017). Different nematic domains will have different interference pattern, and thus can also be distinguished.

In summary, a nematic vortex can be realized in odd-parity topological superconductors. It is an unusual type of topological defect that preserves time reversal symmetry and binds a Kramers pair of Majorana modes. Other interesting properties of nematic odd-parity topological superconductors are likely to await discovery.

VI Acknowledgment

We acknowledge support from Department of Energy, Office of Basic Energy Science, Materials Science and Engineering Division.

Appendix A Analytic solution to zero modes in nematic vortex

In this appendix, we present the details of obtaining analytic solution to Majorana zero modes. Equation (13) is an ansatz to zero modes of Hamiltonian (11). Wave functions $\xi_{1,2}$ satisfy the following equations:

[TABLE]

where $\tilde{h}_{0}$ and $\tilde{h}_{2}$ are given by:

[TABLE]

It is important to note that $h_{1}=i\tilde{h}_{2}$, which makes analytic solution possible. Equation (23) reduces to:

[TABLE]

To proceed we make the following ansatz:

[TABLE]

where $u$, $v$ and $\Lambda$ are parameters that need to be determined. $\mathcal{N}$ is just a normalization factor. Ansatz in (26) is proposed based on the following identity:

[TABLE]

Equation (25) then simplifies to an algebraic equation:

[TABLE]

When $u$ and $v$ are decomposed into their real and imaginary parts, equation (28) can be viewed as an eigenproblem with $\Lambda$ as the eigenvalue. There are two eigenvalues $\pm\sqrt{\lambda_{1}^{2}+\lambda_{2}^{2}}$. We choose $\Lambda=\sqrt{\lambda_{1}^{2}+\lambda_{2}^{2}}$ so that the zero modes are exponentially confined to the center. Given $\Lambda=\sqrt{\lambda_{1}^{2}+\lambda_{2}^{2}}$, there are two solutions for $(u,v)$:

[TABLE]

where $\lambda_{1,2}\geq 0$ is assumed without loss of generality. $(u_{1},v_{1})$ leads to the solution in (14). $(u_{2},v_{2})$ does not give rise to an independent solution but just switches the definition of $\gamma_{\uparrow}^{\dagger}$ and $\gamma_{\downarrow}^{\dagger}$.

Appendix B Cone geometry

In this appendix, we present some mathematical properties on the cone geometry. The coordinate on a cone is $(\rho,\phi)$, where $\rho$ is the radial coordinate measured from the tip of the cone, and $\phi$ is the azimuthal angle. The gradient and Laplace operators on the cone are

[TABLE]

where $\hat{\rho}$ and $\hat{\phi}$ are respectively the unit vectors along the radial and the azimuthal direction. $\theta$ is half of the opening angle of the cone, as shown in Fig. 3(b). The following identities are used to derive the Laplace operator from the gradient operator:

[TABLE]

We now show how the Hamiltonian in (15) and (16) are derived on the cone geometry. The kinetic energy has a simple connection to the Laplace operator: $h_{0}=\beta k_{\perp}^{2}-\mu=-\beta\nabla^{2}-\mu$. The pairing order parameter $\boldsymbol{\eta}$ varies in space in the following way:

[TABLE]

where $\phi_{0}$ is the angle between the order parameter $\boldsymbol{\eta}$ and the radial direction. In the main text, we have set $\phi_{0}=0$. $h_{1}$ has the expression:

[TABLE]

When $\phi_{0}=0$, the nematic order parameter $\boldsymbol{\eta}$ is parallel to the radial direction, providing a confinement potential. As $\phi_{0}$ increases, the confinement along radial direction becomes weaker. In the special case of $\phi_{0}=\pi/2$, the confinement is completely lost and there are no zero-energy modes. Therefore, the presence of zero modes in the nematic vortex of a polar $p$-wave superconductor requires that $\phi_{0}$ is not $\pi/2$.

Finally, the orthonormal basis functions on the truncated cone are:

[TABLE]

where $q$ and $j$ are integer numbers. $f_{q,j}$ vanishes at the two edges of the truncated cone. The energy spectrum on the cone geometry is obtained by expanding wave functions in terms of $f_{q,j}$.

Appendix C Analytic solution to the flat band of zero modes

In this appendix, we present analytic solutions to the flat band of zero modes in the nematic vortex of a polar $p$-wave superconductor. Here we consider a disk geometry with no hole ($L=R$) and an infinite radius($R\rightarrow+\infty$). The $q$th zero modes centered around the origin is presented in (19). The wave function $\psi_{q}(\rho)$ satisfies the following equation:

[TABLE]

which has the analytic solution

[TABLE]

where $J_{q}$ is the $q$th Bessel function. It is interesting to note that there is a zero-energy solution for every integer $q$ on the disk with infinite radius. As $|q|$ increases, $\psi_{q}(\rho)$ moves away from the origin. On a finite-radius disk, there are also zero modes bound to the outer edge, which can hybridize with $\psi_{q}(\rho)$ if the latter has a significant probability near the outer edge. Therefore, the number of zero modes on a finite-radius disk is finite, but it increases as $R$ increases.

To gain a deeper insight, we revisit the problem of a superconducting vortex with an odd winding number in a polar $p$-wave superconductorVolovik (1999). The effective Hamiltonian is:

[TABLE]

where $\kappa=\eta\exp[i(2q+1)\phi]$ represents the vortex with winding number $2q+1$, where $q$ is an integer. In polar coordinates, we have:

[TABLE]

It is known that $H_{q}$ supports one zero energy modeVolovik (1999). We look for zero mode of the following form:

[TABLE]

Here the wave function $\psi_{q}(\rho)$ should satisfy

[TABLE]

which shares exactly the same form as (35). This provides another point of view on the reason why a polar $p$-wave superconductor supports a flat band of zero modes in the presence of a nematic vortex.