Bimodal conductance distribution of Kitaev edge modes in topological superconductors

Gapless edge states are a striking manifestation of topological protection in two-dimensional systems. First studied in connection with the quantum Hall effect in a strong magnetic field [1, 2], they are now known to exist also in the presence of time-reversal symmetry (for topological insulators (TIs)) or particle-hole symmetry (for topological superconductors) [3, 4]. The edge current can carry charge or heat, it can be uni-directional (chiral) or bi-directional (helical), but in each manifestation there is no backscattering—so that the corresponding electrical or thermal conductance is quantized, independent of system size.

Isotropic two-dimensional superconductors with spin-triplet p-wave pairing belong to the class of topological superconductors, with charge-neutral gapless edge states. Depending on the absence or presence of time-reversal symmetry, the edge modes can be chiral (class D) or helical (class DIII). It is known that the quantization of the thermal conductance breaks down if the two-dimensional superconductor is strongly anisotropic; the edge states remain, but backscattering by disorder is no longer forbidden by a topological invariant [5–8]. One might surmise that the edge states will localize on length scales L larger than the mean free path ℓ, with an exponentially decaying conductance $\propto {\rm exp} \left( -L/\ell \right)$, but that is not what happens. Instead, in [6] an anomalously slow (super-Ohmic) scaling $\propto \sqrt{\ell /L}$ was found, unlike that of any known one-dimensional system. A statistical symmetry (translational invariance of the disorder distribution) was identified as the origin of the topological protection [6].

Here we study this remarkable delocalization of edge states in the framework of the Kitaev Hamiltonian [9] of randomly coupled Majorana fermions. We contrast the two realizations of the model illustrated in figure 1: at the edge of a two-dimensional superconductor (Kitaev edge, figure 1(a)) and as a one-dimensional chain of nanoparticles (Kitaev chain, figure 1(b)). While the Kitaev chain allows for delocalization, this requires a fine-tuning to the critical point of the topological phase diagram [10, 11]. Generically, the conductance of the Kitaev chain has a log-normal distribution peaked at an exponentially small value [12–14], because disorder drives the system away from the gapless critical point into the gapped phase. In contrast, for the Kitaev edge we find a bimodal conductance distribution, with a second peak near the quantized conductance of the clean system. (Compare dashed and solid curves in figure 1(c).) The $\sqrt{\ell /L}$ weight of this second peak produces the super-Ohmic conductance scaling of [6].

Zoom In Zoom Out Reset image size

We explain the difference in conductance distributions in terms of the different way in which translational invariance of the disorder distribution is realized in the two systems: in the Kitaev edge all nearest-neighbor coupling strengths of Majorana fermions are statistically equivalent, while in the Kitaev chain even and odd-numbered couplings are inequivalent. Finally, we show how the charge-neutral edge modes of the topological superconductor can be detected in an electrical—rather than thermal—measurement, by considering the shot noise of time-dependent current fluctuations.

The outline of this paper is as follows. In section 2 we introduce model Hamiltonians for p-wave superconductors with chiral or helical edge states and calculate the topological phase diagram in the presence of both anisotropy and disorder. The topological phase transitions are identified by considering the bulk conductance and the associated topological invariants. Edge conductance in the topologically non-trivial phases is studied in section 3. In section 4 we contrast the conductance distributions of the Kitaev edge and the Kitaev chain. Electrical, rather than thermal, detection of the edge modes is discussed in section 5. We conclude in section 6.

The topological phase diagram of clean chiral p-wave superconductors (or superfluids) was studied in [5, 7, 8]. Here we show how the topologically distinct phases evolve when we include disorder, for both chiral and helical pair potentials, which as we will see have a qualitatively different phase diagram. Numerical calculations on a disordered tight-binding model are compared with analytical calculations of the phase boundaries in self-consistent Born approximation (SCBA).

2.1. Model Hamiltonians

Superconductors with broken spin–rotation symmetry are in symmetry class D or DIII in the Altland–Zirnbauer classification [17], depending on whether time-reversal symmetry is broken or not. In both symmetry classes the Bogoliubov-De Gennes Hamiltonian $\mathcal{H}\left( {\boldsymbol{k}} \right)$ has electron-hole symmetry,

$$\begin{eqnarray}&&{{\tau }_{x}}{{\mathcal{H}}^{{}}}\left( {\boldsymbol{k}} \right){{\tau }_{x}}=-{{\mathcal{H}}^{*}}\left( -{\boldsymbol{k}} \right),\end{eqnarray} \tag{ 1 }$$

where the Pauli matrix ${{\tau }_{i}}$ acts on the electron-hole degree of freedom. In class DIII there is additionally the time-reversal symmetry

$$\begin{eqnarray}&&{{\sigma }_{y}}\mathcal{H}\left( {\boldsymbol{k}} \right){{\sigma }_{y}}={{\mathcal{H}}^{*}}\left( -{\boldsymbol{k}} \right),\end{eqnarray} \tag{ 2 }$$

with ${{\sigma }_{i}}$ acting on the spin degree of freedom and $\hbar {\boldsymbol{k}}$ the momentum.

The minimal class-D Hamiltonian, constrained by equation (1), has the form

$$\begin{eqnarray}&&{{\mathcal{H}}_{{\rm D}}}({\boldsymbol{k}} )=\epsilon ({\boldsymbol{k}} ){{\tau }_{z}}+{{\Delta }_{x}}{{\tau }_{x}}{\rm sin} {{k}_{x}}+{{\Delta }_{y}}{{\tau }_{y}}{\rm sin} {{k}_{y}},\end{eqnarray} \tag{ 3a }$$

$$\begin{eqnarray}&&\epsilon ({\boldsymbol{k}} )=-2{{t}_{x}}{\rm cos} {{k}_{x}}-2{{t}_{y}}{\rm cos} {{k}_{y}}-\mu ,\end{eqnarray} \tag{ 3b }$$

where $\left( {{\Delta }_{x}},{{\Delta }_{y}} \right)=\left( \Delta ,\alpha \Delta \right)$ is the anisotropic amplitude of the chiral p-wave pair potential (in a gauge where it is real), $\left( {{t}_{x}},{{t}_{y}} \right)=\left( t,\alpha t \right)$ is the anisotropic hopping amplitude, and μ is the chemical potential. The parameter $\alpha \in \left[ 0,1 \right]$ measures the degree of anisotropy, with $\alpha \to 1$ the isotropic limit. We consider equal-spin pairing, so the spin degree of freedom does not appear in ${{\mathcal{H}}_{{\rm D}}}$.

In class DIII the additional constraint (2) is satisfied by taking two time-reversed copies of the Hamiltonian (3),

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{H}}_{{\rm DIII}}}({\boldsymbol{k}} ) & = & \epsilon ({\boldsymbol{k}} )({{\sigma }_{0}}\otimes {{\tau }_{z}})+{{\Delta }_{x}}({{\sigma }_{z}}\otimes {{\tau }_{x}}){\rm sin} {{k}_{x}} \\ {} & {} & +{{\Delta }_{y}}({{\sigma }_{0}}\otimes {{\tau }_{y}}){\rm sin} {{k}_{y}}+K({{\sigma }_{y}}\otimes {{\tau }_{y}}), \\ \end{array}\end{eqnarray} \tag{ 4 }$$

coupled with strength K.

The Hamiltonians (3) and (4) are discretized on a two-dimensional square lattice of size ${{L}_{x}}\times {{L}_{y}}$ (lattice constant $a\equiv 1$). Electrostatic disorder (strength δ) is added by randomly varying μ, independently for each lattice site and uniformly in the interval $\left[ \mu -\delta ,\mu +\delta \right]$. We study thermal conduction by attaching disorder-free leads at two ends of the lattice, connected to reservoirs at temperature ${{T}_{0}}$ and ${{T}_{0}}+\delta T$ (see figure 2). The scattering matrix,

$$\begin{eqnarray}S=\left( \begin{array}{ccccccccccccccc} r & t \\ t^{\prime} & r^{\prime} \\ \end{array} \right),\end{eqnarray} \tag{ 5 }$$

evaluated at the Fermi level (E = 0) determines the thermal conductance

$$\begin{eqnarray}&&G={{G}_{0}}\;{\rm Tr}\;{{t}^{\dagger }}t,\ \ {{G}_{0}}={{\pi }^{2}}k_{{\rm B}}^{2}{{T}_{0}}/6\;h,\end{eqnarray} \tag{ 6 }$$

in the low-temperature, linear response regime. The numerical calculations are performed using the kwant tight-binding code [18].

Zoom In Zoom Out Reset image size

2.2. Class D phase diagram

2.2.1. Clean limit

We first discuss the phase diagram of the class-D Hamiltonian (3) in the clean limit of refs. [5, 7, 8], before including the effects of disorder. Without disorder the momentum is a good quantum number and one can search for gap closings in the Brillouin zone. These occur at the four high-symmetry points ${{k}_{x}},{{k}_{y}}\in \left\{ 0,\pi \right\}$, for chemical potentials $\mu =\pm 2t\left( 1\pm \alpha \right)$. In the μ–α plane the four gapless lines are boundaries separating five topologically distinct insulating phases, see figure 3.

Zoom In Zoom Out Reset image size

The number of chiral edge modes is given by the Chern number ν [3, 4, 19], being the winding number of eigenstates in the Brillouin zone. We compute this two-dimensional topological invariant from the scattering matrix rather than from an integral over the Brillouin zone, in a formulation that can be applied directly to disordered systems [20, 21],

$$\begin{eqnarray}&&\nu =\frac{1}{2\pi i}\int _{0}^{2\pi }{\rm d}\phi \;\frac{{\rm d}}{{\rm d}\phi }\;{\rm ln} \;{\rm det} \;r\left( \phi \right).\end{eqnarray} \tag{ 7 }$$

Here $r\left( \phi \right)$ is the reflection block of the scattering matrix (5) in the geometry of figure 2(a), with leads attached to x = 0 and $x={{L}_{x}}$ and twisted periodic boundary conditions ³ on the scattering state $\psi \left( x,y \right)$ in the y-direction: $\psi \left( x,0 \right)={{{\rm e}}^{i\phi }}\psi \left( x,{{L}_{y}} \right)$.

For $2t\left( 1-\alpha \right)<|\mu |<2t\left( 1+\alpha \right)$ the system is topologically non-trivial, with $\nu ={\rm sign}\;\mu$ and a chiral Majorana edge mode. (The sign of ν gives the direction of propagation.) The absence of backscattering leads to a quantized thermal edge conductance $G={{G}_{0}}$. This characterizes the strong TI.

When $|\mu |>2t\left( 1+\alpha \right)$ or $|\mu |<2t\left( 1-\alpha \right)$ the Chern number $\nu =0$, so there are no chiral edge modes. These regions in the phase diagram are distinguished by an alternative 'weak' topological invariant ${{\nu }_{y}}$ [22, 23]. Again, to prepare ourselves for disorder effects, we use a scattering matrix formulation rather than a Brillouin zone formulation [5, 7]. The two-dimensional weak topological invariant is the strong topological invariant in one dimension lower, which in class D is given by the determinant of the reflection matrix [11],

$$\begin{eqnarray}&&{{\nu }_{y}}={\rm det} \;r\left( \phi =0 \right).\end{eqnarray} \tag{ 8 }$$

The dimensional reduction is implemented by evaluating $r\left( \phi \right)$ at $\phi =0$, so for periodic boundary conditions in the y-direction. When the Chern number $\nu =0$ the weak invariant ${{\nu }_{y}}$ may be equal to $+1$ (trivial insulator) or $-1$ (weak topological insulator, WTI).

For $|\mu |>2t\left( 1+\alpha \right)$ we are in the topologically trivial phase, with ${{\nu }_{y}}=1$ and no edge modes at all. In contrast, when $|\mu |<2t\left( 1-\alpha \right)$ the system is a WTI, with ${{\nu }_{y}}=-1$ and non-chiral Majorana modes on the edges in the y-direction.

2.2.2. Disorder effects

Having described the phase diagram of the system in the clean limit, we now turn to the effects of disorder. Sufficiently strong disorder can convert a class-D superconductor that is insulating in the bulk into a thermal metal [24–26]. To search for this topological phase transition we take the two-terminal geometry of figure 2(a) with periodic boundary conditions in the y-direction, in order to focus on the metallic or insulating nature of the bulk. (We will consider edge conduction in section 3.)

Numerical results for the disorder-averaged thermal conductance $\left\langle G \right\rangle$ are shown in figure 4, as a function of chemical potential μ and disorder strength δ. One can see that both the TI and WTI phases are robust to disorder, up to about $\delta \approx t$. For stronger disorder there is a TI-to-thermal metal phase transition, followed by a transition to a topologically trivial Anderson insulator. The phase boundaries between TI, WTI, and thermal metal are in quite good agreement with those calculated in SCBA (dashed lines in figure 4, see the appendix for details of the calculation). The transition to an Anderson insulator at strong disorder is out of reach of that approximation.

Zoom In Zoom Out Reset image size

The distinct topological nature of the TI, WTI, and trivial phase is confirmed by a calculation of the topological invariants $\nu ,{{\nu }_{y}}$, see figure 5. In the bulk insulating phases these are quantized numbers: $\nu \in \left\{ -1,0,1 \right\}$ (a so-called $\mathbb{Z}$ invariant), while ${{\nu }_{y}}\in \left\{ -1,1 \right\}$ (a ${{\mathbb{Z}}_{2}}$ invariant). At the topological phase transitions, when the bulk gap closes, both ν and ${{\nu }_{y}}$ are free to vary between these integer values, resulting in the smooth transitions shown in figure 5.

Zoom In Zoom Out Reset image size

2.3. Class DIII phase diagram

2.3.1. Clean limit

We now turn to the phase diagram of the class DIII Hamiltonian (4), first without disorder. It is convenient to rotate the Hamiltonian to a block off-diagonal form,

$$\begin{eqnarray}U{{\mathcal{H}}_{{\rm DIII}}}{{U}^{\dagger }}=\left( \begin{array}{ccccccccccccccc} 0 & A \\ {{A}^{\dagger }} & 0 \\ \end{array} \right),\;\;U={\rm exp} (-\frac{1}{4}{\rm i}\pi \;{{\sigma }_{x}}\otimes {{\tau }_{x}}),\end{eqnarray} \tag{ 9a }$$

$$\begin{eqnarray}\begin{array}{rcl} A & = & {\rm i}\Delta {{\sigma }_{0}}{\rm sin} {{k}_{x}}+\alpha \Delta {{\sigma }_{z}}{\rm sin} {{k}_{y}} \\ {} & {} & +(\mu +2t{\rm cos} {{k}_{x}}+2\alpha t{\rm cos} {{k}_{y}}-{\rm i}k){{\sigma }_{y}}. \\ \end{array}\end{eqnarray} \tag{ 9b }$$

At the gap closings of ${{H}_{{\rm DIII}}}$ the determinant of A vanishes, which happens when

$$\begin{eqnarray}&&\mu /t=-2{\rm cos} {{k}_{x}}-2\alpha {\rm cos} {{k}_{y}},\end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray}&&{{K}^{2}}/{{\Delta }^{2}}={{{\rm sin} }^{2}}{{k}_{x}}+{{\alpha }^{2}}{{{\rm sin} }^{2}}{{k}_{y}}.\end{eqnarray} \tag{ 10b }$$

The gap closings identify the boundaries of insulating phases, as shown in figure 6. While in class D the gap closes along a line in the phase diagram, in class DIII there are extended gapless regions of a metallic phase separating the insulating phases. (This is a generic feature of helical p-wave superconductors [27].)

Zoom In Zoom Out Reset image size

We can distinguish five distinct insulating regions of phase space. For weak anisotropy, $\alpha >K/\Delta$, we find two $\pm \mu$ symmetric insulating phases bounded by

$$\begin{eqnarray}&&2-2\sqrt{{{\alpha }^{2}}-{{K}^{2}}/{{\Delta }^{2}}}<|\mu /t|<2+2\sqrt{{{\alpha }^{2}}-{{K}^{2}}/{{\Delta }^{2}}}.\end{eqnarray} \tag{ 11 }$$

When the anisotropy reaches the critical value

$$\begin{eqnarray}&&{{\alpha }_{c}}={{\left( 1-{{K}^{2}}/{{\Delta }^{2}} \right)}^{1/2}},\end{eqnarray} \tag{ 12 }$$

a third insulating phase appears centered around $\mu =0$, in the interval

$$\begin{eqnarray}&&|\mu /t|<2\left( {{\alpha }_{c}}-\alpha \right),\ \ {\rm for}\ \ \alpha <{{\alpha }_{c}}.\end{eqnarray} \tag{ 13 }$$

Additionally, for any amount of anisotropy there are insulating phases at large chemical potentials, with boundaries given by

$$\begin{eqnarray}|\mu /t|<\left\{ \begin{array}{ccccccccccccccc} 2\alpha +2\sqrt{1-{{K}^{2}}/{{\Delta }^{2}}}, & {\rm for}\;\;\alpha <{{\alpha }_{c}} \\ 2\sqrt{2}\sqrt{1+{{\alpha }^{2}}-{{K}^{2}}/{{\Delta }^{2}}}, & {\rm for}\;\;\alpha >{{\alpha }_{c}}. \\ \end{array} \right.\end{eqnarray} \tag{ 14 }$$

To test the topological properties of these phases, we compute the associated DIII topological invariants in a scattering formulation (so that we can directly apply it to disordered systems in the next subsection). The strong topological invariant [21],

$$\begin{eqnarray}&&Q={\rm Pf}\;\left[ {\rm i}{{\sigma }_{y}}r\left( \phi =0 \right) \right]\times {\rm Pf}\;\left[ {\rm i}{{\sigma }_{y}}r\left( \phi =\pi \right) \right],\end{eqnarray} \tag{ 15 }$$

is determined by the Pfaffians of the reflection matrix with periodic ($\phi =0$) and anti-periodic ($\phi =\pi$) boundary conditions in the y-direction. In view of the time-reversal symmetry condition (2), the matrix ${\rm i}{{\sigma }_{y}}r\left( \phi \right)$ is anti-symmetric for $\phi =0,\pi$, so the Pfaffian exists.

The insulating regions delimited by equation (11) are topologically non-trivial ($Q=-1$), with helical Majorana edge states and quantized thermal conductance $G=2{{G}_{0}}$. All other insulating phases have Q = 1. The ones appearing at large chemical potentials, bounded by equation (14), are topologically trivial, without edge states. However, the phase which develops at the critical anisotropy ${{\alpha }_{c}}$, bounded by equation (13), has Q = 1 but still supports gapless modes on edges oriented in the y-direction. The weak topological invariant ${{Q}_{y}}=-1$ of this phase is obtained by dimensional reduction to the one-dimensional class-DIII topological invariant [28],

$$\begin{eqnarray}&&{{Q}_{y}}={\rm Pf}\;\left[ {\rm i}{{\sigma }_{y}}r\left( \phi =0 \right) \right].\end{eqnarray} \tag{ 16 }$$

2.3.2. Disorder effects

Figure 7 shows the effect of disorder on the topological phases, probed by calculating the thermal conductance in the geometry of figure 2(a) with periodic boundary conditions in the y-direction. Comparison with the class-D phase diagram of figure 4 shows as a qualitative difference that the thermal metal phase extends down to zero disorder. This behavior is also captured by the SCBA (dashed curves), see the appendix for details of the calculation.

Zoom In Zoom Out Reset image size

The scattering matrix formulas (15) and (16) for the class-DIII strong and weak topological invariants are applied to a disordered system in figure 8. These are both ${{\mathbb{Z}}_{2}}$ invariants, meaning that they take on the values $\pm 1$ when the bulk is insulating: $Q={{Q}_{y}}=1$ in the trivial insulator, $Q=-1$, ${{Q}_{y}}=\pm 1$ in the TI, and Q = 1, ${{Q}_{y}}=-1$ in the WTI.

Zoom In Zoom Out Reset image size

So far we studied thermal conduction in the geometry of figure 2(a) with periodic boundary conditions in the transverse direction, in order to eliminate edge contributions and focus on bulk properties. To study edge conduction in the TI and WTI phases we now take the geometry of figure 2(b), with leads attached to $y=0,{{L}_{y}}$ and hard-wall boundary conditions at $x=0,{{L}_{x}}$. We again first consider the clean case and then add the effects of disorder, for both symmetry classes D and DIII.

3.1. Clean case

The TI and WTI phases both have gapless edge states, the difference being that the TI edge states appear on all edges while the WTI edge states exist only at two of the four edges (see figure 9). In the geometry of figure 2(b) we can probe the edge conductance in both phases. Without disorder the conductance is system-size independent, because there is no backscattering, and the difference between the TI and WTI phases is simply a factor of two.

Zoom In Zoom Out Reset image size

This conductance doubling at the TI-to-WTI transition is shown in figure 10. In class D it happens because the chiral edge state of the TI phase can propagate in both directions in the WTI phase. In class DIII we start out with helical edge states in the TI phase, with direction of propagation tied to the spin degree of freedom. In the WTI phase this helicity is lost, so now both spin-up and spin-down can propagate in both directions and the conductance is doubled.

Zoom In Zoom Out Reset image size

3.2. Disorder effects

The addition of disorder has no effect on the conductance in the TI phase, since backscattering of chiral or helical edge states is forbidden. The edge states in the WTI are neither chiral nor helical, so disorder does cause backscattering and reduces the edge conductance. However, as discovered in [6], the recovery of translational invariance upon ensemble averaging prevents localization of the WTI edge states. Instead of an exponential decay of the conductance with the length $L\equiv {{L}_{y}}$ of the edge, there is only an algebraic $1/\sqrt{L}$ decay. In figure 11 we show this super-Ohmic conductance scaling for the WTI phase of the class-D and class-DIII Hamiltonians (3) and (4).

Zoom In Zoom Out Reset image size

The absence of localization at the edge of the anisotropic p-wave superconductor is puzzling if one tries to understand it starting from the limit $\alpha \to 0$ of strong anisotropy. Then the system can be thought of as an array of weakly coupled nanowires with overlapping Majorana zero-modes at the end points, a so-called Kitaev ladder [29–31]. The effective edge Hamiltonian is the Kitaev Hamiltonian [9] in class D, or two time-reversed copies of it in class DIII. The disordered one-dimensional Kitaev model, called the Kitaev chain, is known to be an insulator [12–14]—so how do the Kitaev edge modes avoid localization?

To clarify the situation we contrast the two class-D systems. (Class DIII is similar.) The Kitaev Hamiltonian

$$\begin{eqnarray}&&{{H}_{{\rm K}}}=\mathop{\sum }\limits_{n=1}^{2N}{\rm i}{{\kappa }_{n}}{{\gamma }_{n}}{{\gamma }_{n+1}}\end{eqnarray} \tag{ 17 }$$

describes the nearest-neigbor coupling (coupling strength ${{\kappa }_{n}}$) of $2N$ Majorana fermion operators ${{\gamma }_{n}}$. These are Hermitian operators, ${{\gamma }_{n}}=\gamma _{n}^{\dagger }$, with anti-commutation relation ${{\gamma }_{n}}{{\gamma }_{m}}+{{\gamma }_{m}}{{\gamma }_{n}}=2{{\delta }_{nm}}$. To obtain a closed system the Majorana's are assumed to lie on a ring, so that ${{\kappa }_{2N}}$ couples ${{\gamma }_{2N}}$ to ${{\gamma }_{2N+1}}\equiv {{\gamma }_{1}}$.

This one-dimensional system in symmetry class D has a ${{\mathbb{Z}}_{2}}$ topological invariant [9],

$$\begin{eqnarray}&&{{Q}_{{\rm K}}}={\rm sign}\left( {\rm Pf}\;{{A}^{+}} \right)\left( {\rm Pf}\;{{A}^{-}} \right),\end{eqnarray} \tag{ 18 }$$

determined by the Pfaffians of a pair of real antisymmetric matrices ${{A}^{\pm }}$, having non-zero matrix elements

$$\begin{eqnarray}&&A_{n,n+1}^{\pm }=-A_{n+1,n}^{\pm }={{\kappa }_{n}},\ \ 1\leqslant n\leqslant 2N-1,\end{eqnarray} \tag{ 19a }$$

$$\begin{eqnarray}&&A_{2N,1}^{\pm }=-A_{1,2N}^{\pm }=\pm {{\kappa }_{2N}}.\end{eqnarray} \tag{ 19b }$$

Evaluation of the Pfaffians gives

$$\begin{eqnarray}&&{{Q}_{{\rm K}}}={\rm sign}\left( \mathop{\prod }\limits_{n=1}^{N}\kappa _{2n-1}^{2}-\mathop{\prod }\limits_{n=1}^{N}\kappa _{2n}^{2} \right).\end{eqnarray} \tag{ 20 }$$

Translational invariance of the disorder ensemble means two completely different things for the Kitaev edge and for the Kitaev chain. For the Kitaev edge it means that the coupling strengths ${{\kappa }_{n}}$ between adjacent Majorana's all have the same distribution. The disorder average $\left\langle {{Q}_{{\rm K}}} \right\rangle$ of the topological invariant then vanishes, which is why the Kitaev edge is called a critical WTI [6]. In contrast, as illustrated in figure 1, for the Kitaev chain translational invariance means that the ${{\kappa }_{n}}$ with n even or those with n odd have the same distribution, but the distributions of ${{\kappa }_{2n}}$ and ${{\kappa }_{2n-1}}$ are unrelated. The topological invariant is then non-zero on average, so the Kitaev chain is not critical [12, 14].

The implication for the transmission probability T follows if we remove the coupling between ${{\gamma }_{2N}}$ and ${{\gamma }_{1}}$, so that we can introduce transmission and reflection amplitudes $t=1/{\rm cosh} \alpha$, $r={\rm tanh} \alpha$. The Lyapunov exponent α determines both T and ${{Q}_{{\rm K}}}$,

$$\begin{eqnarray}&&{{Q}_{{\rm K}}}={\rm sign}\;\alpha ,\ \ T=1/{{{\rm cosh} }^{2}}\;\alpha ,\end{eqnarray} \tag{ 21 }$$

and has a Gaussian distribution $P\left( \alpha \right)$ [10, 14]. The variance ${\rm Var}\;\alpha =L/\ell$ is determined by the mean free path ℓ for backscattering along the edge, of length $L\gg \ell$. The mean $\left\langle \alpha \right\rangle =L/\xi$ defines the localization length ξ. A vanishing $\left\langle {{Q}_{{\rm K}}} \right\rangle$ implies that the median of $P\left( \alpha \right)$ is zero, and since it's Gaussian also $\left\langle \alpha \right\rangle =0\Rightarrow \xi =\infty$.

For the Kitaev chain $\left\langle {{Q}_{{\rm K}}} \right\rangle \ne 0$ and hence ξ is finite, so the transmission probability has a log-normal distribution peaked at $T={{{\rm e}}^{-2L/\xi }}$, with an exponentially decaying average transmission [12, 14]. In contrast, for the Kitaev edge $\left\langle {{Q}_{{\rm K}}} \right\rangle =0\Rightarrow \left\langle \alpha \right\rangle =0\Rightarrow \xi =\infty$. The Gaussian distribution of the Lyapunov exponent then produces a bimodal distribution of the transmission probability for $L\gg \ell$,

$$\begin{eqnarray}\begin{array}{rcl} P(T) & = & \int _{-\infty }^{\infty }{\rm d}\alpha \;\delta (T-1/{{{\rm cosh} }^{2}}\;\alpha ){{\left( 2\pi L/\ell \right)}^{-1/2}}{{{\rm e}}^{-{{\alpha }^{2}}\ell /2L}} \\ {} & = & {{\left( \ell /2\pi L \right)}^{1/2}}\;{{T}^{-1}}{{\left( 1-T \right)}^{-1/2}} \\ {} & {} & \times {\rm exp} \left[ -(\ell /2L){\rm arcos}{{{\rm h}}^{2}}\;({{T}^{-1/2}}) \right], \\ \end{array}\end{eqnarray} \tag{ 22 }$$

peaked near T = 0 and T = 1, with average decaying algebraically as $\left\langle T \right\rangle =\sqrt{2\ell /\pi L}$.

We have tested the result (22) in a computer simulation of the anisotropic p-wave superconductor, with class-D Hamiltonian (3). In the geometry of figure 2(b) both edges at x = 0 and $x={{L}_{x}}$ contribute to the thermal conductance in the WTI phase, but one edge can be removed by reducing the width of the contacts to the interval $0\leqslant x\leqslant {{L}_{x}}/2$. Results are shown in figure 12. With the mean free path ℓ as a single fit parameter, the transition from a uni-modal distribution to a bimodal distribution upon increasing $L/\ell$ is well-described by equation (22).

Zoom In Zoom Out Reset image size

So far we considered thermal conduction as the probe of edge state transport. Electrical detection would be more convenient experimentally, and this is possible by adapting the nanowire setup of [11]. All contacts are now at the same temperature ${{T}_{0}}$, the superconductor is grounded as well as one of the metal contacts (number 2), and the other contact (number 1) is biased at voltage ${{V}_{1}}$. The electrical current into the grounded contact 2 fluctuates in time with noise power ${{P}_{12}}$. This is dominated by shot noise, at low temperatures when thermal noise can be neglected (${{k}_{{\rm B}}}{{T}_{0}}\ll e{{V}_{1}}$).

The noise power is given in terms of the transmission matrix by [32]

$$\begin{eqnarray}\begin{array}{rcl} {{P}_{12}}/{{P}_{0}} & = & {\rm Tr}\;\left( t_{ee}^{\dagger }t_{ee}^{{}}+t_{he}^{\dagger }t_{he}^{{}} \right)-{\rm Tr}\;{{\left( t_{ee}^{\dagger }t_{ee}^{{}}-t_{he}^{\dagger }t_{he}^{{}} \right)}^{2}} \\ {} & = & \frac{1}{2}\;{\rm Tr}\;{{t}^{\dagger }}t-\frac{1}{2}\;{\rm Tr}{{({{\tau }_{z}}{{t}^{\dagger }}{{\tau }_{z}}t)}^{2}}, \\ \end{array}\end{eqnarray} \tag{ 23 }$$

with ${{P}_{0}}={{e}^{3}}{{V}_{1}}/h$. The subscripts $e,h$ indicate transmission from electron to electron (${{t}_{ee}}$) or from electron to hole (${{t}_{he}}$), and we used electron-hole symmetry in the second equation to rewrite the whole expression in terms of the full transmission matrix t.

As derived in [11], when the transmission is via an unpaired Majorana mode, the second trace in equation (23) vanishes identically so the electrical shot noise is directly related to the thermal conductance: ${{P}_{12}}/{{P}_{0}}=\frac{1}{2}G/{{G}_{0}}$. This applies to symmetry class D. More generally, in both symmetry classes D and DIII the two quantities ${{P}_{12}}$ and G have the same $1/\sqrt{L}$ scaling in the WTI phase ⁴ , compare figures 11 and 13.

Zoom In Zoom Out Reset image size

The Kitaev model [9] is paradigmatic for topological superconductivity and Majorana zero-modes, and for that reason has been studied extensively [33, 34]. Here we have shown that the realization of this model at the edge of a two-dimensional system (what we have called the 'Kitaev edge') is fundamentally different from its strictly one-dimensional counterpart, the Kitaev chain. The difference, summarized in figure 1, manifests itself in the different distribution of the thermal conductance, peaked at exponentially small value in the Kitaev chain [12–14] while the Kitaev edge has a second peak at the conductance quantum.

As a possible physical realization of Kitaev edge modes we have studied in some detail a model of an anisotropic two-dimensional chiral p-wave superconductor [5–8], as well as its time-reversally symmetric (helical) counterpart. Both can produce WTI with Kitaev edge modes, but while they appear at any amount of anisotropy for chiral p-wave pairing, the helical p-wave pairing requires a threshold anisotropy (compare figures 3 and 6). We have demonstrated the robustness of the WTI phase to disorder by numerical simulations, in good agreement with analytical calculations of the phase boundaries in self-consistent Born approximation (figures 4 and 7).

Experimentally the transition into the WTI phase can be detected, on length scales below the mean free path, via the doubling of the thermal conductance (figure 10), and on larger length scales via the super-Ohmic scaling (figure 11). Because of the complexity of thermal transport measurements at low temperatures, we have proposed an alternative fully electrical method of detection, using the electrical shot noise power (figure 13).

Discussions with A R Akhmerov are gratefully acknowledged. This research was supported by the Foundation for Fundamental Research on Matter (FOM), the Netherlands Organization for Scientific Research (NWO/OCW), and an ERC Synergy Grant.

We calculate the phase diagram in the presence of electrostatic disorder (strength δ as defined in the main text) using the SCBA. Below we provide details of the calculation for the class-DIII Hamiltonian (4). The corresponding results for the class-D Hamiltonian (3) are simply obtained by taking the vanishing coupling $K\to 0$ limit and are summarized at the end of this appendix.

We calculate the disorder-averaged density of states from the self-energy $\Sigma$, defined by

$$\begin{eqnarray}&&\begin{array}{l} \frac{1}{E+{\rm i}{{0}^{+}}-{{H}_{{\rm DIII}}}-\Sigma } \\ \qquad =\left\langle \frac{1}{E+{\rm i}{{0}^{+}}-{{H}_{{\rm DIII}}}-{{H}_{{\rm disorder}}}} \right\rangle . \\ \end{array}\end{eqnarray} \tag{ A.1 }$$

The SCBA self-energy at the Fermi level (E = 0) is given by

$$\begin{eqnarray}&&\Sigma =\frac{1}{3}{{\delta }^{2}}\mathop{\sum }\limits_{{\boldsymbol{k}} }{{\tau }_{z}}\frac{1}{{\rm i}{{0}^{+}}-{{H}_{{\rm DIII}}}({\boldsymbol{k}} )-\Sigma }{{\tau }_{z}}.\end{eqnarray} \tag{ A.2 }$$

The sum over ${\boldsymbol{k}}$ ranges over the first Brillouin zone and in the continuum limit

$$\begin{eqnarray}&&\mathop{\sum }\limits_{{\boldsymbol{k}} }\mapsto \frac{1}{4{{\pi }^{2}}}\int _{-\pi }^{\pi }{\rm d}{{k}_{x}}\int _{-\pi }^{\pi }{\rm d}{{k}_{y}}.\end{eqnarray} \tag{ A.3 }$$

The SCBA self-energy is a ${\boldsymbol{k}}$-independent 4×4 matrix with spin and electron-hole degrees of freedom

$$\begin{eqnarray}&&\Sigma =({{\sigma }_{0}}\otimes {{\tau }_{z}})\delta \mu -({{\sigma }_{y}}\otimes {{\tau }_{y}})\delta K-({{\sigma }_{0}}\otimes {{\tau }_{0}}){\rm i}\gamma .\end{eqnarray} \tag{ A.4 }$$

The terms $\delta \mu$ and $\delta K$ renormalize the chemical potential and coupling respectively. Both terms account for a disorder induced shift of the phase boundaries between the trivial insulator, TI, WTI, and thermal metal. The term γ produces a finite density of states, induced by the disorder. Such a finite density of states may indicate a thermal metal or a trivial Anderson insulator, but it cannot distinguish between the two.

We substitute the self-energy (A.4) into equation (A.2) and observe that the right-hand-side depends only on the renormalized chemical potential $\tilde{\mu }=\mu -\delta \mu$ and coupling $\tilde{K}=K-\delta K$. Denoting the renormalized Hamiltonian by ${{\tilde{H}}_{{\rm DIII}}}\left( {\boldsymbol{k}} \right)\equiv {{H}_{{\rm DIII}}}\left( {\boldsymbol{k}} \right)+\left( {{\sigma }_{0}}\otimes {{\tau }_{z}} \right)\delta \mu -\left( {{\sigma }_{y}}\otimes {{\tau }_{y}} \right)\delta K$ we can write the following identity:

$$\begin{eqnarray}&&\begin{array}{l} {{\left( {\rm i}\gamma -{{\tilde{H}}_{{\rm DIII}}}({\boldsymbol{k}} ) \right)}^{-1}}=\left( {\rm i}\gamma +{{\tilde{H}}_{{\rm DIII}}}({\boldsymbol{k}} ) \right) \\ \cdot \left( {{\gamma }^{2}}+{{f}_{\tilde{\mu },\tilde{K}}}({\boldsymbol{k}} )-2\tilde{K}({{\Delta }_{x}}{\rm sin} {{k}_{x}}\;{{\sigma }_{z}}{{\tau }_{z}}+{{\Delta }_{y}}{\rm sin} {{k}_{y}}\;{{\sigma }_{y}}) \right) \\ \cdot {{\left( 4{{\tilde{K}}^{2}}\;g({\boldsymbol{k}} )-{{\left( {{\gamma }^{2}}+{{f}_{\tilde{\mu },\tilde{K}}}({\boldsymbol{k}} ) \right)}^{2}} \right)}^{-1}}, \\ \end{array}\end{eqnarray} \tag{ A.5 }$$

where

$$\begin{eqnarray}&&g({\boldsymbol{k}} )=\Delta _{x}^{2}{{{\rm sin} }^{2}}\;{{k}_{x}}+\Delta _{y}^{2}{{{\rm sin} }^{2}}\;{{k}_{y}},\end{eqnarray} \tag{ A.6a }$$

$$\begin{eqnarray}&&{{f}_{\tilde{\mu },\tilde{K}}}({\boldsymbol{k}} )={{\epsilon }_{\tilde{\mu }}}{{({\boldsymbol{k}} )}^{2}}+{{\tilde{K}}^{2}}+g({\boldsymbol{k}} ).\end{eqnarray} \tag{ A.6b }$$

Because of symmetry only terms even in ${\boldsymbol{k}}$ contribute to ${{\sum }_{{\boldsymbol{k}} }}$. We find from equation (A.2) three coupled equations for the parameters $\delta \mu$, $\delta K$, and γ, which completely determine the SCBA self-energy:

$$\begin{eqnarray}&&\delta \mu =\frac{1}{3}{{\delta }^{2}}\mathop{\sum }\limits_{{\boldsymbol{k}} }\frac{{{\epsilon }_{\tilde{\mu }}}({\boldsymbol{k}} )\left( {{\gamma }^{2}}+{{f}_{\tilde{\mu },\tilde{K}}}({\boldsymbol{k}} ) \right)}{4{{\tilde{K}}^{2}}\;g({\boldsymbol{k}} )-{{\left( {{\gamma }^{2}}+{{f}_{\tilde{\mu },\tilde{K}}}({\boldsymbol{k}} ) \right)}^{2}}},\end{eqnarray} \tag{ A.7a }$$

$$\begin{eqnarray}&&\delta K=\frac{1}{3}{{\delta }^{2}}\mathop{\sum }\limits_{{\boldsymbol{k}} }\tilde{K}\frac{{{\gamma }^{2}}+{{f}_{\tilde{\mu },\tilde{K}}}({\boldsymbol{k}} )-2\;g({\boldsymbol{k}} )}{4{{\tilde{K}}^{2}}\;g({\boldsymbol{k}} )-{{\left( {{\gamma }^{2}}+{{f}_{\tilde{\mu },\tilde{K}}}({\boldsymbol{k}} ) \right)}^{2}}},\end{eqnarray} \tag{ A.7b }$$

$$\begin{eqnarray}&&\gamma =\frac{1}{3}{{\delta }^{2}}\mathop{\sum }\limits_{{\boldsymbol{k}} }\frac{-\gamma ({{\gamma }^{2}}+{{f}_{\tilde{\mu },\tilde{K}}}({\boldsymbol{k}} ))}{4{{\tilde{K}}^{2}}\;g({\boldsymbol{k}} )-{{\left( {{\gamma }^{2}}+{{f}_{\tilde{\mu },\tilde{K}}}({\boldsymbol{k}} ) \right)}^{2}}}.\end{eqnarray} \tag{ A.7c }$$

We address first the shift of phase boundaries in the weak disorder case. To this end we set $\gamma =0$ assuming that the disorder is too weak to induce a finite density of states. We are looking for solutions of the clean system gap closing conditions equations (11), (13), (14) expressed in terms of the renormalized parameters $\mu \to \tilde{\mu }$, $K\to \tilde{K}$. The gap closing condition defines the phase boundary and can be expressed as $\tilde{K}=\tilde{K}\left( \tilde{\mu } \right)$ (with a different function $\tilde{K}\left( \tilde{\mu } \right)$ for each boundary). We rewrite the SCBA equations (A.7a ) with $\gamma =0$ in the form

$$\begin{eqnarray}&&\mu =\frac{1}{3}{{\delta }^{2}}F(\tilde{\mu },\tilde{K})+\tilde{\mu },\end{eqnarray} \tag{ A.8a }$$

$$\begin{eqnarray}&&{{\delta }^{2}}=3(K-\tilde{K}){{\left( G(\tilde{\mu },\tilde{K}) \right)}^{-1}},\end{eqnarray} \tag{ A.8b }$$

where

$$\begin{eqnarray}&&F(\tilde{\mu },\tilde{K})=\mathop{\sum }\limits_{{\boldsymbol{k}} }\frac{{{\epsilon }_{\tilde{\mu }}}({\boldsymbol{k}} ){{f}_{\tilde{\mu },\tilde{K}}}({\boldsymbol{k}} )}{4{{\tilde{K}}^{2}}\;g({\boldsymbol{k}} )-f_{\tilde{\mu },\tilde{K}}^{2}({\boldsymbol{k}} )},\end{eqnarray} \tag{ A.9 }$$

$$\begin{eqnarray}&&G(\tilde{\mu },\tilde{K})=\mathop{\sum }\limits_{{\boldsymbol{k}} }\tilde{K}\frac{{{f}_{\tilde{\mu },\tilde{K}}}({\boldsymbol{k}} )-2\;g({\boldsymbol{k}} )}{4{{\tilde{K}}^{2}}\;g({\boldsymbol{k}} )-f_{\tilde{\mu },\tilde{K}}^{2}({\boldsymbol{k}} )}.\end{eqnarray} \tag{ A.10 }$$

In this way we obtain the parametric solution for the disorder strength ${{\delta }^{2}}\left[ \tilde{\mu },\tilde{K}\left( \tilde{\mu } \right) \right]$ and the gap parameter $\mu \left[ \tilde{\mu },{{\delta }^{2}}\left( \tilde{\mu } \right) \right]$ along the phase boundary. We vary the parameter $\tilde{\mu }$ away from the clean system solution μ at $\delta =0$. The sums over the Brillouin zone are computed numerically in the continuum limit. The resulting parametric phase boundaries $\left[ \mu \left( \tilde{\mu } \right),\delta \left( \tilde{\mu } \right) \right]$ separate insulating and gapless phases at low to moderate disorder.

For sufficiently large disorder $\delta >{{\delta }_{c}}$ the SCBA equations (A.7a ) may support solutions with non-zero γ indicating the onset of a finite density of states at zero energy. This marks the transitions from the strong and WTIs to the thermal metal at strong disorder. The ${{\delta }_{c}}$ dependence on μ and K follows from solutions of the SCBA equations at infinitesimal $\gamma \ne 0$

$$\begin{eqnarray}&&1=\frac{\delta _{c}^{2}}{3}\mathop{\sum }\limits_{{\boldsymbol{k}} }\frac{{{f}_{\tilde{\mu },\tilde{K}}}({\boldsymbol{k}} )}{f_{\tilde{\mu },\tilde{K}}^{2}({\boldsymbol{k}} )-4{{\tilde{K}}^{2}}\;g({\boldsymbol{k}} )}\equiv \frac{\delta _{c}^{2}}{3}H(\tilde{\mu },\tilde{K}),\end{eqnarray} \tag{ A.11a }$$

$$\begin{eqnarray}&&\tilde{\mu }=-\frac{1}{3}\delta _{c}^{2}F(\tilde{\mu },\tilde{K})+\mu ,\end{eqnarray} \tag{ A.11b }$$

$$\begin{eqnarray}&&\tilde{K}=-\frac{1}{3}\delta _{c}^{2}G(\tilde{\mu },\tilde{K})+K.\end{eqnarray} \tag{ A.11c }$$

To determine ${{\delta }_{c}}$ we first search numerically (using Steffensen iteration) for fixed point solutions of equations (A.11b ), (A.11c ),

$$\begin{eqnarray}&&\left( \begin{array}{ccccccccccccccc} \tilde{\mu } \\ \tilde{K} \\ \end{array} \right)=\frac{-1}{H(\tilde{\mu },\tilde{K})}\left( \begin{array}{ccccccccccccccc} F(\tilde{\mu },\tilde{K}) \\ G(\tilde{\mu },\tilde{K}) \\ \end{array} \right)+\left( \begin{array}{ccccccccccccccc} \mu \\ K \\ \end{array} \right),\end{eqnarray} \tag{ A.12 }$$

for a given value of the chemical potential μ and coupling K. Finally we compute ${{\delta }_{c}}$ from (A.11a ) for the obtained solutions $\left( \tilde{\mu },\tilde{K} \right)$.

Both the parametric solutions to the renormalized gap closing conditions and the computed μ-dependence (for fixed K) of the critical disorder ${{\delta }_{c}}$ are shown as white dashed lines in figure 7.

This was all for class DIII. The formulas for class D correspond to the $K\to 0$ limit. Denoting by ${{E}_{\mu }}\left( {\boldsymbol{k}} \right)$ the excitation spectrum of the class-D Hamiltonian (3) we find

$$\begin{eqnarray}&&E_{\mu }^{2}({\boldsymbol{k}} )={{f}_{\mu ,K\to 0}}({\boldsymbol{k}} )=\epsilon _{\mu }^{2}({\boldsymbol{k}} )+\Delta _{x}^{2}{{{\rm sin} }^{2}}\;{{k}_{x}}+\Delta _{y}^{2}{{{\rm sin} }^{2}}\;{{k}_{y}}.\end{eqnarray} \tag{ A.13 }$$

The shift of the phase boundaries at low disorder is obtained by imposing the gap closing conditions for the renormalized chemical potential $\tilde{\mu }=\pm 2t\left( 1\pm \alpha \right)$. From equation (A.8a ) we directly obtain the boundary position as a function of the disorder strength

$$\begin{eqnarray}&&\mu ={{\left. \left[ -\frac{1}{3}{{\delta }^{2}}\mathop{\sum }\limits_{{\boldsymbol{k}} }\frac{{{\epsilon }_{\tilde{\mu }}}({\boldsymbol{k}} )}{E_{\tilde{\mu }}^{2}({\boldsymbol{k}} )}+\tilde{\mu } \right] \right|}_{\tilde{\mu }=\pm 2t(1\pm \alpha )}}.\end{eqnarray} \tag{ A.14 }$$

For the isotropic case $\alpha =1$ the central $\mu =0$ transition separating the two TI phases is not renormalized by disorder. The transitions separating the TI phases from the trivial insulator phase are shifted according to $\mu =\pm \left( 4t+0.130\;{{\delta }^{2}} \right)$. For $\alpha =1/2$ the transition line between the WTI and TI phases is given by $\mu =\pm \left( 3t+0.039\;{{\delta }^{2}} \right)$ and the transition between TI and the trivial phase by $\mu =\pm \left( t+0.180\;{{\delta }^{2}} \right)$. These phase boundaries are shown as dashed lines in figure 4.

The critical disorder lines separating the TI and WTI phases from the thermal metal at large disorder can be found from the SCBA solutions equations (A.11a ), (A.11b ) at infinitesimal $\gamma \ne 0$. In class-D we can directly parametrize such solutions by the renormalized chemical potential $\tilde{\mu }$. We obtain

$$\begin{eqnarray}&&\delta _{c}^{2}(\tilde{\mu })=3{{\left( \mathop{\sum }\limits_{{\boldsymbol{k}} }\frac{1}{E_{\tilde{\mu }}^{2}({\boldsymbol{k}} )} \right)}^{-1}},\end{eqnarray} \tag{ A.15 }$$

$$\begin{eqnarray}&&\mu \left[ \tilde{\mu },{{\delta }_{c}}(\tilde{\mu }) \right]=\tilde{\mu }-\frac{1}{3}\delta _{c}^{2}(\tilde{\mu })\mathop{\sum }\limits_{{\boldsymbol{k}} }\frac{{{\epsilon }_{\tilde{\mu }}}({\boldsymbol{k}} )}{E_{\tilde{\mu }}^{2}({\boldsymbol{k}} )}.\end{eqnarray} \tag{ A.16 }$$

The parametric dependence for $\tilde{\mu }\in \left[ -4,4 \right]$ calculated from the above equations is also included in figure 4.