Signatures of topological phase transitions in mesoscopic superconducting rings

Our analysis starts with the Kitaev model of a 1D spinless p-wave superconductor [1, 16]

$$\begin{equation} H_{\rm TB}=-\mu_{\rm TB} \sum_{j=1}^N c^\dagger_j c_j-\sum_{j=1}^{N-1} (t c_j^\dagger c_{j+1}+\Delta_{\rm TB} c_jc_{j+1}+\mathrm{h.c.}), \end{equation} \noindent \tag{ 1 }$$

which describes a wire of N sites. Electrons on site j are annihilated by c_j, hop between neighboring sites with hopping amplitude t and have chemical potential μ_TB. For all numerical calculations in this paper we choose t = 1. The p-wave pairing strength is given by Δ_TB. Here, we label both the chemical potential and the pairing strength by the subscript TB to distinguish these parameters of the tight-binding model (1) from their analogues in the continuum model introduced below. The wire can be closed into a ring with a weak link by an additional hopping term between sites 1 and N,

$$\begin{equation} H_{\rm T}=-t' c_N^\dagger c_1 + \mathrm{h.c.}, \end{equation} \noindent \tag{ 2 }$$

with hopping amplitude t'. We assume that charging effects are weak and can be neglected (see [41, 42] for the consequences of charging in ring-like structures).

For an infinite and uniform wire, the Kitaev Hamiltonian (1) exhibits a phase transition when the chemical potential μ_TB crosses one of the band edges. The system is in a topological (nontopological) superconducting phase when the chemical potential is within (outside) the interval [ − 2t,2t], i.e. within (outside) the band at vanishing pairing Δ_TB = 0. The spectrum exhibits a superconducting gap on both sides of the topological phase transition, while the gap closes at the topological critical points |μ_TB| = 2t. It is thus natural to introduce the chemical potential measured from the lower band edge, i.e. μ = μ_TB + 2t.

In the vicinity of the band edges (say the lower band edge) and thus of the topological phase transition, we can make a continuum approximation to the tight-binding model (1). We will mostly employ the tight-binding model in the first part of the paper, while we partially find it more convenient to rely on the continuum approximation in dealing with effects of disorder in section 4. The continuum model is formulated in terms of the corresponding Bogoliubov–de Gennes Hamiltonian [1, 16]

$$\begin{equation} H=\left[\begin{array}{@{}cc@{}} \displaystyle \frac{p^2}{2m}+V(x)-\mu & \displaystyle \frac{1}{2}\left\lbrace\Delta'(x) ,p\right\rbrace\\ \displaystyle \frac{1}{2}\left\lbrace\Delta'(x),p\right\rbrace & -\left(\displaystyle \frac{p^2}{2m}+V(x)-\mu\right) \end{array}\right], \end{equation} \noindent \tag{ 3 }$$

where Δ'(x) is the p-wave pairing strength and the curly brackets denote the anticommutator. Here, we have included a disorder potential V (x) which we will return to in more detail in section 4. For V (x) = 0 the bulk spectrum of the continuum model is given by

$$\begin{equation} \epsilon_p=\pm\left[\! \left(\frac{p^2}{2m}-\mu\right)^2+|\Delta'|^2 p^2\!\right]^{1/2}, \end{equation} \noindent \tag{ 4 }$$

which becomes gapless for μ = 0. This indicates the above-mentioned topological phase transition between a topological phase with μ > 0 and a nontopological phase for μ < 0.

In a semi-infinite wire, the topological phase is characterized by an MBS localized near its end point. The MBS has zero energy and a wave function that decays exponentially into the wire on the scale of the superconducting coherence length ξ. In a finite wire, the MBS localized at the two ends of the wire hybridize and form a conventional Dirac fermion whose energy ₀ scales as the overlap of the two Majorana end states which is exponentially small in the length L of the wire. The wavefunction of the MBS depends on the parameter regime (see, e.g., [43]). This is easily seen by determining the allowed wavevectors at zero energy from equation (4), which yields

$$\begin{equation} p_0 =\pm \mathrm{i}\,m|\Delta'|\pm \sqrt{2m\mu-m^2|\Delta'|^2}. \end{equation} \tag{ 5 }$$

• (i)  
  $\mu\gg m\Delta'^2$: Deep in the topological phase, the bulk excitation spectrum equation (4) has two minima around $\pm p_{\mathrm {F}}=\pm \sqrt {2m\mu }$ with a gap Δ⁽ⁱ⁾_eff ≈ p_FΔ' (see figure 1(a)). According to equation (5), the Majorana wavefunctions decay on the scale ξ = 1/mΔ' and oscillate with a much shorter period 1/p_F. In a finite wire the hybridization energy is given by (see the appendix) $\epsilon _0= 2\Delta 'p_{\mathrm {F}}|\sin (p_{\mathrm {F}} L)|\exp (-L/\xi )$, which has accidental degeneracies at integer values of p_FL/π.
• (ii)  
  $\mu\ll m\Delta'^2$: Near the topological phase transition at μ = 0, the excitation spectrum has only a single minimum at p = 0 with a gap of order μ (see figure 1(b)). At low energies, we can neglect the kinetic energy in equation (3) and the spinless p-wave superconductor can be approximately described by the Dirac Hamiltonian

  $$\begin{equation} H\simeq -\mu\tau_z+\Delta' p\tau_x. \end{equation} \tag{ 6 }$$

  Equation (5) gives p₀ ≈ ± i μ/Δ', so that the spatial extent of the Majorana wavefunction is governed by the coherence length ξ = Δ'/μ, which diverges at the topological phase transition. In contrast to the previous regime, the end-state energy does not exhibit oscillations, $\epsilon _0 \propto \exp (-L/\xi )$.

Zoom In Zoom Out Reset image size

In the presence of a magnetic flux threading the ring, both the tunneling amplitude and the pairing strength become complex and acquire a phase. The precise nature of these phases depends on the physical realization of the Kitaev chain. We illustrate this point by discussing two possible setups based on the proposal to realize the Kitaev chain in a semiconductor wire proximity coupled to an s-wave superconductor [9, 10], as illustrated in figure 2:

• (a)  
  The s-wave superconductor is interrupted underneath the weak link in the quantum wire. Current can flow around the loop only through the semiconductor weak link.
• (b)  
  The s-wave superconductor forms a closed ring and a weak link exists only in the semiconductor wire. The current through the weak link of the semiconductor will, in general, be only a small perturbation of the current flowing through the superconductor.

We assume that the thickness of the superconducting ring is small compared to both its London penetration depth and its superconducting coherence length ξ_SC. The supercurrent flowing in the superconductor is given by [44]

$$\begin{equation} J_{\mathrm{s}}=\frac{2e}{m^*}|\psi|^2\left(\hbar\nabla\varphi-2 eA\right), \end{equation} \tag{ 7 }$$

where m* and |ψ|² are the effective mass and density of the superconducting electrons and φ denotes the phase of the s-wave order parameter. The p-wave pairing potential in the quantum wire inherits its phase φ from the s-wave superconductor underneath via the proximity effect. (The effective p-wave order parameter may have an additional phase shift that depends on geometric details such as the direction of the Zeeman field and the spin–orbit coupling; however, these contributions lead to constant offsets of the phase which are unaffected by the magnetic flux.) The vector potential A oriented along the wire is related to the Aharonov–Bohm flux ϕ through

$$\begin{equation} \phi=\oint\mathrm{d}x A(x), \end{equation} \tag{ 8 }$$

where the integral is taken around the ring of circumference L.

Zoom In Zoom Out Reset image size

The phase of the order parameter φ is different for the two setups illustrated in figure 2. In setup (a), no supercurrent is able to flow since the loop is interrupted, J_s = 0. If we choose a gauge in which the vector potential is uniform around the ring, A(x) = ϕ/L, the phase φ of the order parameter becomes φ(x) = 4π(ϕ/ϕ₀)(x/L) in terms of the normal-metal flux quantum ϕ₀ = h/e. In setup (b), the supercurrent around the ring is governed by fluxoid quantization φ(x + L) = φ(x) + 2πn, with the integer n labeling the fluxoid states. In a gauge in which A(x) = ϕ/L, this implies that ∇φ = 2πn/L, yielding a supercurrent of J_s = (2e/m*)|ψ|²[2πℏn/L − 2eA]. Here, [x] denotes the integer closest to x. In thermodynamic equilibrium, the system realizes the fluxoid state of lowest energy and thus of lowest supercurrent, i.e. n = [ϕ/(ϕ₀/2)].

Within the chosen gauge, in setup (a) the hopping amplitude and the pair potential in the tight-binding Hamiltonian in equation (1) take the form t → t e^{i 2πϕ/Nϕ₀} and Δ_TB → Δ_TB e^{i 4π(ϕ/ϕ₀)(j/N)}. Alternatively, one can perform the gauge transformation c_j → c_j e^{−i (j−1/2)2πϕ/Nϕ₀}, which eliminates the phase from the pair potential. In this new gauge, both the pair potential and the hopping amplitude t in the interior of the ring are real while the hopping amplitude across the weak link acquires a phase factor, t' → t' e^{i 2πϕ/ϕ₀}. Our numerical results will be obtained for this representation of the tight-binding model. In contrast, in setup (b), we find that Δ_TB → Δ_TB e^{i 2π[ϕ/(ϕ₀/2)](j/N)} for the pair potential (note the closest integer symbol [.] in the exponent), while t → t e^{i 2πϕ/Nϕ₀} as well as t' → t' e^{i 2πϕ/Nϕ₀}. As in the previous case (a), we can eliminate the phase of the pair potential by a gauge transformation. However, this no longer eliminates the phase of the hopping matrix element t. Instead, one finds that t → t e^{i (π/N){ϕ/(ϕ₀/2)−[ϕ/(ϕ₀/2)]}} and t' → t' e^{i (π/N){ϕ/(ϕ₀/2)+(N−1)[ϕ/(ϕ₀/2)]}}. The fact that we can no longer eliminate the magnetic flux from the bulk of the wire is a manifestation of the fact that supercurrents in the s-wave superconductor modify the spectrum of the quantum wire [8].

Clearly, the effective Kitaev chain is quite different for setups (a) and (b). In the remainder of this paper, we will focus on setup (a) where the flux enters only into the tunneling Hamiltonian representing the weak link. In this setting, the current in the semiconductor wire of interest here is experimentally more accessible since there is no background current in the bulk s-wave superconductor unlike in setup (b).

We first briefly review the Josephson effect of two semi-infinite wires connected at their ends through a weak link (or equivalently, a ring of infinite circumference), as originally considered by Kitaev [1]. The corresponding low-energy excitation spectrum as a function of flux is sketched in figure 3(a). Due to the Majorana end states, there are two subgap states whose energies are governed by the tunneling amplitude across the weak link. While each individual level is periodic in flux with period h/e, the overall spectrum is h/2e periodic. As a result, the thermodynamic ground state energy of the system—and thus the Josephson current in strict thermodynamic equilibrium—are h/2e periodic.

Zoom In Zoom Out Reset image size

At the same time, the h/e periodicity of the individual subgap states is a direct consequence of the Majorana nature of the endstates. This signature of MF can be brought out in measurements of the Josephson current if the fermion parity of the system is a good quantum number. The level crossing of the two Majorana subgap states in figure 3(a) is then protected by fermion-parity conservation. As a result, since there is only a single level crossing per superconducting flux quantum, the system necessarily goes from the ground state to an excited state (or vice versa) when changing the flux by h/2e. During this process, the excited state is unable to relax to the ground state since this would require a change in fermion parity. Thus, the system only returns to its initial state after a change in flux of h/e, which corresponds to the fractional Josephson effect.

For rings with finite circumference, the two MBS localized at the two banks of the weak link hybridize not only through the tunnel coupling across the weak link but also because of the overlap of their wavefunctions in the topological superconductor. In the previous subsection, we considered the situation in which this interior hybridization is vanishingly small compared to tunneling across the weak link. Conversely, when the tunneling across the weak link is negligible compared to the interior hybridization, the splitting of the Majoranas due to the interior overlap does not depend on flux. Weak tunneling across the junction will then cause a small h/e-periodic modulation of the split Majorana levels with flux. In this situation, even the thermodynamic ground state energy becomes h/e periodic, regardless of the presence or absence of fermion-parity violating processes. In fact, of the two h/e-periodic levels, the negative-energy level (which is occupied in equilibrium) corresponds to an even-parity ground state while the positive-energy level is occupied in the odd-parity first excited state. Weak fermion-parity violating processes will not destroy the h/e-periodic Josephson current as the two levels no longer cross as a function of flux.

The full crossover of the Bogoliubov–de Gennes spectrum as the interior overlap of the MBS increases is illustrated with numerical results in figure 3(b) (see section 3.4 for details of the numerical calculations). They confirm the above picture for the limit of strong overlap. But they also show that an h/e-periodic contribution to the equilibrium Josephson current exists even when the interior splitting is of the order of or smaller than the tunnel coupling across the weak link. Indeed, the interior overlap essentially pushes one of the two states (dashed line) up in energy, while it pushes its particle–hole conjugate state (solid line) down. At small interior overlaps, this shifts the two level crossings (initially at ϕ₀/4 and 3ϕ₀/4) outwards towards a flux of zero and one flux quantum ϕ₀. Note that the level crossings remain intact, protected by fermion-parity conservation. However, once the level crossings reach a flux of zero and ϕ₀, respectively, the levels merely touch at these points. Thus, fermion parity no longer protects the levels from splitting, and indeed one state remains at finite and negative energies at all values of flux while, symmetrically, its particle–hole conjugate state remains at finite and positive energies.

Consider now the Josephson current as a function of flux in the presence of weak but finite fermion-parity violating processes. Specifically, we assume that the flux is varied by h/e on a time scale that is large compared to the relaxation time of the fermion parity while, at the same time, the fermion-parity violating processes are weak compared to the hybridization of the MBS so that the Bogoliubov–de Gennes spectra in figure 3 are relevant. In this case, the Josephson current is essentially h/2e periodic deep in the topological phase, where L ≫ ξ. However, as the system approaches the topological phase transition, ξ grows and, hence, the hybridization of the MBS through the interior of the ring increases. As a result, the h/e-periodic contribution to the current increases. Conversely, the MBS disappear on the nontopological side of the phase transition where the Josephson current thus reverts to h/2e periodicity. As a result, we expect a peak in the h/e-periodic Josephson current near the topological phase transition, whose measurement would constitute a clear signature of the topological phase transition and the formation of MF.

This expectation is confirmed by the numerical results shown in figure 4(a), where the corresponding Fourier coefficient $A_{h/e}=(2e/h)\int _0^{h/e}\mathrm {d}\phi I(\phi )\sin (2\pi e\phi /h)$ of the equilibrium Josephson current I(ϕ) is plotted as a function of chemical potential. A_h/e(μ) exhibits a peak in the topological phase (μ > 0), which moves closer to the topological phase transition at μ = 0 as the ring circumference increases (see figure 4(b)).

Zoom In Zoom Out Reset image size

Deep in the topological phase the MBS are localized at the weak link. Approaching the phase transition at μ = 0, the MBS delocalize. On the one hand, this causes an increase in the overlap of the MBS in the interior of the topological superconductor. As discussed above, this leads to an increase of the h/e-periodic Josephson current. On the other hand, however, the probability density of the MBS near the weak link decreases, causing a suppression of the hybridization of the MBS across the weak link and hence of the h/e-periodic Josephson current. Thus, the peak occurs for the value of μ where the interior overlap splitting is equal to the tunnel coupling. Since the interior overlap is exponentially small in L/ξ while the hybridization across the weak link is roughly independent of the ring's circumference L, the peak position shifts towards the phase transition point at μ = 0 with increasing L (cf figure 4(c)). Since at the same time the h/e-periodic Josephson current becomes suppressed when the systems is approaching the phase transition, the peak is more pronounced in shorter rings.

A more quantitative description can be developed by restricting the Hamiltonian to the low-energy subspace spanned by the two MBS. The projection of the tunneling Hamiltonian across the weak link onto this subspace yields

$$\begin{equation} H_{\rm T}=-\Gamma\cos(2\pi \phi/\phi_0) (d_{\rm M}^\dagger d_{\rm M}-1/2), \end{equation} \noindent \tag{ 9 }$$

where d_M is the Dirac fermion constructed from the two MBS. The parameter Γ measures the tunnel coupling of the MBS across the weak link and is given by (cf equation (A.4) in the appendix)

$$\begin{equation} \Gamma=\frac{ t'\mu (4t-\mu) \Delta_{\mathrm{TB}} }{ t(t+\Delta_{\mathrm{TB}})^2} . \end{equation} \noindent \tag{ 10 }$$

Here the factor of μ accounts for the probability density of the Majorana wavefunction at the junction, which vanishes at the phase transition.

The overlap of the Majorana end-states in the interior of the wire leads to an additional coupling (cf the appendix)

$$\begin{eqnarray*} &&H_{\mathrm{overlap}}=\epsilon_0 \left(d_{\mathrm{M}}^{\dag} d_{\mathrm{M}}-1/2\right), \end{eqnarray*} \noindent$$

where

$$\begin{equation} \epsilon_0 =2 \mu \exp(-L/\xi) \end{equation} \noindent \tag{ 11 }$$

measures the strength of the overlap.

Combining these two contributions for a mesoscopic ring near the topological phase transition ($\mu\ll m{\Delta^\prime}^2$), the effective low-energy Hamiltonian reads as

$$\begin{equation} H_{\rm eff}=\left[\epsilon_0-\Gamma\cos\left(\frac{2\pi\phi}{\phi_0}\right)\right]\left(d_{\rm M}^\dagger d_{\rm M}-1/2\right). \end{equation} \tag{ 12 }$$

The Bogoliubov–de Gennes spectrum of this Hamiltonian reproduces the numerically calculated subgap spectra depicted in figure 3(b).

In principle, both the negative energy continuum states as well as the negative energy subgap state contribute to the equilibrium Josephson current. If we denote the sum over all negative excitation energies by E₀(ϕ), we can write the equilibrium Josephson current as I(ϕ) = −∂_ϕE₀(ϕ). However, it is natural to expect and will be corroborated by our numerical results that the Josephson current is dominated by the contribution of the subgap state I(ϕ) ≃ ∂_ϕ|₀ − Γcos(2πϕ/ϕ₀)|/2. Thus, it is straightforward to compute the h/e-periodic Fourier component of the Josephson current,

$$\begin{equation} A_{h/e} = \left\lbrace\begin{array}{@{}ll@{}} \displaystyle \frac{e\Gamma}{\pi \hbar} \left[\displaystyle \frac{\epsilon_0}{\Gamma}\sqrt{1-\displaystyle \frac{\epsilon_0^2}{\Gamma^2}} + \arcsin\left( \displaystyle \frac{\epsilon_0}{\Gamma} \right) \right], & \epsilon_0 <\Gamma,\\ \displaystyle \frac{e\Gamma}{2\hbar}, & \epsilon_0 >\Gamma. \end{array}\right. \end{equation} \tag{ 13 }$$

In the next section, we have compared this analytical result with numerics and found nice agreement.

To obtain numerical results for the Josephson current, we solve the Hamiltonian defined in equations (1) and (2) by exact diagonalization. Figure 4(a) compares the amplitude of the h/e-periodic component as a function of chemical potential with the analytical result in equation (13). The numerical results agree well with the behavior predicted by the low-energy model, except for small deviations in the immediate vicinity of the phase transition at μ = 0. In the inset of figure 4(a), we compare the analytical and numerical results for the quantities Γ and ₀ appearing in the low-energy Hamiltonian. While the model correctly captures ₀ in the regime of interest, there are deviations of Γ near μ = 0. These discrepancies are readily understood as a consequence of the finite circumference of the ring. Although the coherence length diverges at the phase transition, the MBS can delocalize at most throughout the entire length of the ring; there remains a finite probability density of the MBS wavefunction at the weak link.

Figure 4(b) shows that the left flank of the peak of A_h/e and eΓ/h deviates slightly. This deviation is a measure of the size of the bulk contribution to the h/e-periodic current. The latter can thus be seen to be small, justifying our focus on the low-energy Hamiltonian (12) describing the MBS only. In the same figure, the h/2e component is plotted, showing that the h/e-periodic Josephson current exceeds the h/2e component. This is a consequence of the tunneling regime that favors single-electron tunneling over the tunneling of Cooper pairs.

In figure 4(c), we show how the position of the peak in A_h/e depends on the circumference of the ring. We find that the value of μ where the peak occurs scales as 1/L. This result can be understood as follows. Γ is essentially independent of the length of the ring, while ₀ scales as ${\sim }\exp (-L/\xi )$. As we have seen above, the peak occurs at ₀ = Γ. For given t, Δ_TB and t', Γ is fixed and the peak occurs at a constant value of the ratio L/ξ. Since ξ ∼ 1/μ, the value of μ where the peak occurs scales as 1/L. Also note that the above-mentioned tail of the peak at μ ⩽ 0 originating from finite-size corrections is more pronounced in shorter rings.

In this section, we investigate the fate of the peak in the equilibrium h/e-periodic Josephson current in the presence of disorder. Our main results are the following.

• 1.  
  The typical peak height is not affected by disorder as long as the mean free path is longer than the circumference of the ring. Thus the signature persists in the presence of moderate disorder.
• 2.  
  For stronger disorder the peak height decreases and the peak position is shifted to lower chemical potentials.

To study the effect of disorder we add a random onsite potential $\sum _i V_i c_i^\dagger c_i$ to the tight-binding Hamiltonian (1) and (2), where the V_i are taken from a uniform distribution over the interval [ − W,+W]. The mean free path is then related to the disorder strength as l∝1/W².⁵ To obtain numerical results we compute the spectrum by exact diagonalization. Disorder affects the h/e-periodic Josephson current by introducing fluctuations in the quantities ₀ and Γ. While Γ is mainly affected by local fluctuations of the probability density of the Majorana wavefunction at the junction, ₀ fluctuates due to the disorder potential in the entire ring. The interior overlap in disordered wires has been investigated previously for the continuum model (3) in regime (i), i.e. $\mu\gg m\Delta'^2$ [21], where disorder leads to an increase of ₀ and subsequently to a disorder-induced phase transition to the nontopological phase.

Figure 5(a) shows numerical results for the h/e-periodic Josephson current for a few disorder configurations. The peaks in the presence of disorder (green dashed curves) are of comparable height as the peak in the clean ring (black solid curve). The peak shifts as a function of chemical potential, which indicates fluctuations of the coherence length due to disorder.

Zoom In Zoom Out Reset image size

Surprisingly, the peak shifts to lower chemical potentials, corresponding to a decrease in ₀ with disorder in stark contrast to the known case of large μ. This implies that the topological phase is stabilized by disorder if the system is close to the phase transition. To investigate this further we plot the height and position of the peak maxima of many disorder configurations as a color code histogram for l > L in figure 5(b) and l < L in figure 5(c). Indeed the average peak height is comparable to the one in the clean case for l > L. When l ≲ L the average peak height starts to decrease. The histogram in figure 5(c) confirms that the peak is shifted to lower chemical potentials on average.

To understand this behavior we analyze the probability distributions of ₀ and Γ over the disorder ensemble. In figure 6, we show numerical results for the histogram of Γ corresponding to the two ensembles in figures 5(b) and (c) at μ = 0.4. For weak disorder the distribution is symmetric with a mean near the zero-disorder tunnel coupling. For larger disorder when l < L the distribution becomes wider and asymmetric and the average decreases.

Zoom In Zoom Out Reset image size

In order to determine the probability distribution for ₀ we now turn to the continuum Hamiltonian (3) for a wire of length L without tunnel junction. To model short-range correlated disorder in the continuum model, we include a disorder potential with zero average 〈V (x)〉 = 0 and correlation function 〈V (x)V (x')〉 = γδ(x − x'). For this model, we employ a numerical method based on a scattering matrix approach [20, 45, 46]. From the scattering matrix S we obtain the lowest energy eigenstate ₀ by finding the roots of det(1 − S()). In this model, the probability distribution of the hybridization energy ₀ has been shown to be log-normal in [21]. Specifically, it was shown that the log-normal distribution is governed by

$$\begin{equation} \begin{array}{@{}ll@{}} \displaystyle \langle\ln\left(\epsilon_0/2\Delta_{\rm eff}\right)\rangle =-L\left(\frac{1}{\xi}-\frac{1}{2l}\right),\\ \displaystyle {\rm var} \ln\left(\epsilon_0/2\Delta_{\rm eff}\right) =\frac{L}{2l} \end{array} \end{equation} \tag{ 14 }$$

for regime (i). The distribution function reflects the disorder-induced phase transition to the nontopological state at ξ = 2l.

The numerical results are presented in figure 7. In the inset, we show that the mean of ln(₀/2Δ_eff) is indeed linear in L with the slope depending on disorder strength. This slope is plotted as a function of inverse mean free path in figure 7. The data for $\mu=300\,m\Delta'^2$ (blue dots) agree well with the prediction equation (14) with the definitions l = v²_F/γ and ξ = 1/mΔ'.

Zoom In Zoom Out Reset image size

The same plot also shows data corresponding to regime (ii), i.e. $\mu\ll m\Delta'^2$, marked by red crosses. Here, we have $l = \Delta'^2/\gamma$ and ξ = Δ'/μ. Clearly, the behavior is qualitatively different from regime (i), since disorder decreases ₀ rather than increasing it. This is consistent with the shift of the peak of A_h/e to lower μ.

In order to gain analytical insight we now derive the probability distribution of ₀ in regime (ii) extending the results of [21]. The relevant momenta at low energies in this regime are near p = 0 (cf figure 1(b)). Linearizing the dispersion around this point yields the Dirac Hamiltonian equation (6), where the disorder potential enters as a random mass term. Since the disorder potential is short-range correlated, it couples high- and low-momentum degrees of freedom in the original Hamiltonian. Thus a proper linearization of the Hamiltonian requires one to project out the high-momentum states, which renormalizes the gap.

For a strictly linear model with a random mass term, the overlap ₀ has a log-normal distribution [21],

$$\begin{equation} \begin{array}{@{}ll@{}} \displaystyle \langle\ln\left(\epsilon_0/2\Delta_{\rm eff}\right)\rangle =-\frac{L}{\xi},\\ \displaystyle {\rm var} \ln\left(\epsilon_0/2\Delta_{\rm eff}\right) =\frac{L}{l}. \end{array} \end{equation} \tag{ 15 }$$

Thus for the Dirac Hamiltonian the mean of $\ln \left (\epsilon _0/2\Delta _{\mathrm { eff}}\right )$ does not depend on disorder. A systematic linearization of the disordered spinless p-wave superconductor in the vicinity of the topological phase transition effectively renormalizes the chemical potential μ and hence the coherence length ξ = Δ'/μ.

We start by defining the projection operators $P=\sum _{|p|<p_1} |\psi _p\rangle \langle\psi _p|$ onto the low momentum subspace and Q = 1 − P, where {|ψ_p〉}_p is a complete set of momentum eigenstates of the clean Hamiltonian. The relevant momentum scale for this projection is given by p₁ = mΔ', since for p ≪ p₁, the term p²/2m constitutes the lowest energy scale of the Kitaev Hamiltonian. Furthermore, we assume that the disorder potential does not affect high momenta p₁ ≫ 1/l. We can now project the clean Kitaev Hamiltonian H to the low- and high-energy subspaces,

$$\begin{equation} PHP \simeq P\left[(-\mu+V(x))\tau_z+\Delta'p\tau_x\right]P, \end{equation} \tag{ 16 }$$

$$\begin{equation} QHQ \simeq Q\left(p^2/2m\right)\tau_zQ. \end{equation} \noindent \tag{ 17 }$$

Both subspaces are exclusively mixed by the disorder potential PHQ = PV (x)τ_zQ. To second order in V , the correction to the low-energy Hamiltonian is then given by

$$\begin{eqnarray} \delta H(p)&\simeq& \langle \psi_p|PHQ\left(\epsilon_p-QHQ\right)^{-1}QHP|\psi_p \rangle \nonumber\\\ms\ms &\simeq&\sum_{|p'|>p_1}V_{pp'}\frac{1}{\epsilon_p-p'^2/2m\tau_z}V_{p'p}. \end{eqnarray} \noindent \tag{ 18 }$$

Here we used the short notation V_pp' = 〈ψ_p | V (x) | ψ_p'〉. Averaging over disorder, we obtain

$$\begin{equation} \langle \delta H(p) \rangle \simeq-\sum_{|p'|>p_1}\frac{2m\gamma}{p'^2}\tau_z\sim-\frac{\gamma}{\Delta'}\tau_z. \end{equation} \noindent \tag{ 19 }$$

Thus the renormalization produces a contribution to the low-energy Hamiltonian which has the same structure as the chemical potential term. Hence we find a renormalized chemical potential μ' = μ + λγ/Δ' with a numerical factor λ > 0 that cannot be determined from this argument. Thus disorder enters the final result through the renormalized coherence length

$$\begin{equation} \frac{1}{\xi}\rightarrow \frac{1}{\xi}+\frac{\lambda}{l}. \end{equation} \noindent \tag{ 20 }$$

The data in figure 7 confirm equations (15) and (20) and determine the unknown numerical prefactor to be λ = 1/2. Thus for $\mu\ll m\Delta'^2$, ₀ has a log-normal distribution with the mean and variance given by

$$\begin{equation} \begin{array}{@{}ll@{}} \langle\ln\left(\epsilon_0/2\Delta_{\rm eff} \right)\rangle =-L\left(\displaystyle \frac{1}{\xi}+\displaystyle \frac{1}{2l}\right),\\ {\rm var} \ln\left(\epsilon_0/2\Delta_{\rm eff}\right) =\displaystyle \frac{L}{l}. \end{array} \end{equation}\noindent \tag{ 21 }$$

This result is very similar to equation (14) where, however, the disorder correction to the decay length enters with opposite sign. This underlines the contrast between the two regimes, i.e. that disorder drives the system further into the topological phase when it is close to the phase transition, but away from it for larger chemical potentials. Specifically, a spinless p-wave superconducting wire with negative chemical potential may exhibit edge states with an energy exponentially small in L as long as disorder is strong enough.

Combining the disorder-induced fluctuations of Γ and ₀ we can understand the suppression of the peak in the h/e-periodic Josephson current in figure 5(c) for l < L. While ₀ is decreased on average for a given μ with increasing disorder, Γ does not increase at the same time and thus the average peak height decreases. However, the fluctuations of Γ and ₀ become larger as disorder increases such that for single disorder configurations significant peaks are still possible even if the average peak height decreases.

Motivated by the contrasting probability distributions of ₀ in the regimes of large and small μ we numerically calculate the phase diagram of the continuum model (3) as a function of μ and γ, particularly paying attention to the region near the topological phase transition of the clean model. By means of the scattering matrix approach also used in the last section we compute the determinant of the reflection matrix of a wire of length L at  = 0 which approaches the values +1 and −1 as L → ∞ in the nontopological and topological phases, respectively [32, 47].

The resulting phase diagram is plotted in figure 8. From equation (14) we infer that the topological phase transition occurs for ξ = 2l in the regime $\mu\gg m\Delta'^2$. Using the definitions of l and ξ in this regime, we obtain the phase boundary γ⁽ⁱ⁾_c(μ) = 4Δ'μ. This is compared with the numerical results in the inset of figure 8. The numerically calculated phase boundary γ^num_c(μ) has only sublinear deviations from the predicted line, so that the ratio γ^num_c(μ)/2μΔ' = ξ/l approaches the value 2 for μ → ∞ as expected.

Zoom In Zoom Out Reset image size

However, near μ = 0 the behavior is qualitatively different. Here, disorder can induce a topological phase for μ < 0 as well as a re-entrant nontopological phase at larger disorder. From equation (21) we find the condition ξ = −2l for the phase boundary. This corresponds to

$$\begin{equation} \gamma^{(\mathrm{ii})}_{\mathrm{c}}(\mu)=-2\Delta'\mu, \end{equation} \noindent \tag{ 22 }$$

which we find to agree well with the numerical results for the phase diagram (see the dashed line in figure 8). Thus the phase diagram confirms that weak disorder leads to an enhancement of the chemical potential range of the topological phase, while for stronger disorder the range decreases again.

In order to compute the effective Josephson coupling we neglect the effect of the overlap of the two Majorana wavefunctions in the topological part of the ring. We therefore consider the limit of a junction between two half-infinite topological sectors of the wire, the right one on sites j∈[1,∞) and the left one on j∈(∞,−1). They are both described by the Hamiltonian (1) with pairing strength Δ_TB e^i ϕ_a with a = L,R for the left and right sectors, respectively. The hopping between the two sectors from equation (2) now simply reads H_T = −t'(c^†₋₁c₁ + c^†₁c₋₁)—cf figure A.1.

Zoom In Zoom Out Reset image size

We include the effect of the tunneling between the two wires perturbatively. The low-energy excitations of this model for t' = 0 are represented by the left and right zero-energy Majorana operators [1]

$$\begin{equation} b_{\mathrm{R}} =A \sum_{j=1}^\infty (x_+^j-x_-^j)\gamma_{B,j}^{(\mathrm{R})}, \end{equation} \tag{ A.1 }$$

$$\begin{equation} b_{\mathrm{L}}= A \sum_{j=-1}^{-\infty} (x_+^{-j}-x_-^{-j})\gamma_{A,j}^{(\mathrm{L})}, \end{equation} \tag{ A.2 }$$

where $x_{\pm } = (- \mu _{\mathrm {TB}} \pm \sqrt {\mu _{\mathrm {TB}}^2 -4t^2 + 4 \Delta _{\mathrm { TB}}^2 }) / (2t +2 \Delta _{\rm TB})$, $A=(\sum _{j} |x_+^j-x_-^j|^2)^{-1/2}$ is a normalization constant, and the operators γ_A(B),j are defined via c_J = (γ^(R)_B,j e^{−i ϕ_R/2} + i γ^(R)_A,j e^i ϕ_R/2)/2. The projection of H_T onto the subspace spanned by the operators b_R, b_L leads to the effective coupling between the Majorana states. The Hamiltonian can be rewritten in terms of ordinary fermion operators d_M = (b_R + i b_L)/2 to take the form presented in equation (9) with Γ = t'A²|x₊ − x₋|². For simplicity, we consider $\mu <2t(1-\sqrt {1-\Delta _{\mathrm { TB}}^2/t^2})$, which corresponds to the condition $\mu<m\Delta'^2$ in the continuum model (3) and ensures that x_± are real (see the discussion of the continuum model in section 2.1). Explicitly, one obtains

$$\begin{equation} \Gamma=\frac{ t'\mu (4t-\mu) \Delta_{\mathrm{TB}} }{ t(t+\Delta_{\mathrm{TB}})^2}. \end{equation} \tag{ A.3 }$$

In the continuum limit, when μ ≪ Δ_TB ≪ t the expression simplifies to

$$\begin{equation} \Gamma \approx 4 t' \mu \Delta_{\mathrm{TB}} / t^2 . \end{equation} \tag{ A.4 }$$

In order to compute the energy splitting ₀, we employ an alternative method working directly in the continuum limit. Similar to before, we neglect here the Majoranas interaction through the Josephson junction. The problem is then completely equivalent to calculating the energy splitting of two Majorana at the end of a wire of length L. We start considering the regime (ii) described by the Hamiltonian in equation (6), H = μ(x)τ_z + Δ'(x)pτ_x, with the specific choice of parameters (cf figure A.1)

$$\begin{eqnarray} &&\begin{array}{ll} \displaystyle \mu(x) = -V_0 [\Theta(-x)+\Theta(x-L)] +\mu \Theta(x)\Theta(L-x) , \\ \displaystyle \Delta'(x) = \Delta' , \end{array} \end{eqnarray} \tag{ A.5 }$$

where V₀ > 0 and μ > 0 guarantee that the wire is in a topological phase in [0,L] and in a nontopological phase otherwise.

We determine the Majorana wavefunctions by a perturbative approach. We first find the wavefunctions of the MF localized at one of the interfaces between the topological and the insulating region, thus fully ignoring the existence of the other Majorana state. We label the corresponding states at the interfaces at x = 0 and at L as |L〉 and |R〉, respectively. We then project the Hamiltonian onto the Majorana subspace, to obtain their effective interaction, $H_{\mathrm {overlap}}= \sum _{\alpha ,\beta =L,R} \langle \alpha |H |\beta \rangle \, |\alpha \rangle \langle \beta |$. Eventually, we will be interested in the limit V₀ → ∞ corresponding to a high insulating barrier outside the wire.

The Majorana state at x = 0 is the zero-energy eigenstate of the Hamiltonian H_L, defined again by the same Hamiltonian as in equation (6), but with the choice of parameters

$$\begin{eqnarray} &&\begin{array}{@{}ll@{}} \displaystyle \mu(x) = -V_0 \Theta(-x) +\mu \Theta(x) , \\ \displaystyle \Delta'(x) =\Delta' , \end{array} \end{eqnarray} \tag{ A.6 }$$

as depicted in figure A.1. Solving this equation separately for x > 0 and x < 0 leads to zero-energy states with imaginary momenta. Namely, we can write

$$\begin{equation} | L\rangle = \mathbf{v}_{-,L}\, \mathrm{e}^{\alpha_0 x} \Theta(-x) + \mathbf{v}_{+,L}\, \mathrm{e}^{-\alpha x} \Theta(x), \, \end{equation} \tag{ A.7 }$$

where α = μ/|Δ'| is the inverse coherence length in the wire and α₀ = V₀/|Δ'|. The 2D vectors v_±,L are determined by the continuity of the wavefunction and its derivative at the interface and by the wavefunction normalization. They are given by

$$\begin{equation} \mathbf{v}_{+,L}=\mathbf{v}_{-,L}= \sqrt{\frac{\alpha \alpha_0}{\alpha+ \alpha_0}} \left( \begin{array}{@{}c@{}} 1 \\ -\mathrm{i} \end{array} \right). \end{equation} \tag{ A.8 }$$

In the limit V₀ → ∞ that we are interested in, they reduce to

$$\begin{equation} \mathbf{v}_{+,L}=\mathbf{v}_{-,L} = \sqrt{\alpha} \left( \begin{array}{@{}c@{}} 1\\ -\mathrm{i} \end{array} \right). \end{equation} \tag{ A.9 }$$

In a completely analogous way, we can calculate the zero-energy eigenstate of the Hamiltonian H_R defined once more by the Hamiltonian in equation (6), now with (cf figure A.1)

$$\begin{equation} \mu(x) = -V_0 \Theta(x-L) +\mu \Theta(L-x) , \end{equation} \tag{ A.10 }$$

$$\begin{equation} \Delta'(x) =\Delta' . \end{equation} \tag{ A.11 }$$

In this case, the zero-energy eigenstate reads

$$\begin{equation} |R\rangle = \sqrt{\frac{\alpha \alpha_0}{\alpha+ \alpha_0}} \left( \begin{array}{@{}c@{}} 1 \\ \mathrm{i} \end{array} \right) \left[ \mathrm{e}^{-\alpha_0 (x-L)} \Theta(x-L) + \mathrm{e}^{\alpha (x-L)} \Theta(L-x) \right] . \end{equation} \tag{ A.12 }$$

In the limit V₀ → ∞ the prefactor reduces to $\sqrt {\alpha }$. Note that the particle–hole superposition has different phases in |L〉 and |R〉.

We can then project the full Hamiltonian onto the low-energy subspace spanned by the two Majorana states. In doing so, we conveniently rewrite H = H_L + V_R = H_R + V_L, where V_R = −(V₀ + μ)Θ(x − L)τ_z and V_L = −(V₀ + μ)Θ(−x)τ_z, and the spatial dependence of the various terms is presented in figure A.1. We can then compute, e.g., 〈L|H|L〉 = 〈L|H_L|L〉 + 〈L|V_R|L〉 = 〈L|V_R|L〉, and all the other matrix elements in a similar way, to obtain

$$\begin{equation} H_{\mathrm{overlap}} = \left(\begin{array}{@{}cc@{}} \langle L |V_{\mathrm{R}} | L \rangle & \langle L |V_{\mathrm{L}} | R \rangle \\ \langle R |V_{\mathrm{L}} | L \rangle & \langle R |V_{\mathrm{L}} | R \rangle \end{array}\right) = \left(\begin{array}{@{}cc@{}} 0 & -2 \mu \mathrm{e}^{- \alpha L} \\ -2 \mu \mathrm{e}^{- \alpha L} & 0 \end{array}\right) . \end{equation} \tag{ A.13 }$$

This leads to the energy splitting ₀ = 2μe^−αL, as reported in section 3.3.

In complete analogy, one can perform the calculation of the energy splitting in the regime $\mu \gg m \Delta'^2$. In this case the appropriate Hamiltonian, H, is that in equation (3), with the same choice of parameters as in equation (A.6). Again we start considering the left interface, looking for zero-energy solutions of the H_L, i.e. the Hamiltonian in equation (3) with the choice of parameters as in equation (A.7). We introduce $a_0 = V_0/(m \Delta'^2) >0$ and $a = \mu /(m \Delta'^2) > 0$. We can write the solutions as

$$\begin{equation} \fl | L\rangle =\sqrt{m \Delta'}\sqrt{\frac{a}{1+2 a}} \left(\begin{array}{@{}c@{}} 1\\ -\mathrm{i} \end{array}\right) \left[ (1+\eta) \,\mathrm{e}^{\mathrm{i}\,\kappa_-x} \Theta(-x) + \left( \,\mathrm{e}^{\mathrm{i}\, k_+x} +\eta \,\mathrm{e}^{\mathrm{i}\, k_-x}\right) \Theta(x) \right], \end{equation} \tag{ A.14 }$$

$$\begin{eqnarray} &&\fl |R\rangle =\sqrt{m \Delta'}\sqrt{\frac{a}{1+2 a}} \left(\begin{array}{@{}c@{}} 1\\ \mathrm{i} \end{array}\right) \left[ \left( \,\mathrm{e}^{-\mathrm{i}\, k_+(x-L)} +\eta \,\mathrm{e}^{-\mathrm{i}\, k_-(x-L)}\right) \Theta(L-x)+ (1+\eta)\,\mathrm{e}^{-\mathrm{i}\,\kappa_-x} \Theta(x-L) \right] ,\nonumber\\ \end{eqnarray} \tag{ A.15 }$$

where $\kappa _- = -\mathrm {i} m|\Delta '| (\sqrt {2 a_0})$, $k_{\pm }=m |\Delta '| (\pm \sqrt {2 a} +\mathrm {i})$ and η = (κ₋ − k₊)/(k₋ − κ₋). In all the expressions, we are neglecting $\mathcal {O}(\sqrt {a}/\sqrt {a_0}, 1/\sqrt {a_0})$, consistent with the condition a₀ ≫ a ≫ 1, reflecting the limit V₀ → ∞ and the regime under consideration. The matrix elements of the effective Hamiltonian are, in this case, 〈L|H|L〉 = 〈R|H|R〉 = 0 and 〈L|H|R〉 = 〈R|H|L〉 ≡ − ₀, with

$$\begin{equation} \fl \epsilon_0=\Delta' \sqrt{2m \mu} \left[ \sin ( \sqrt{2\,m \mu} L) - \sqrt{a/a_0} (\sqrt{2} -1) \cos ( \sqrt2 {m \mu}) L +\mathcal{O} (\sqrt{a} /a_0^2) \right]. \end{equation} \tag{ A.16 }$$