Multiple Andreev reflection and critical current in topological superconducting nanowire junctions

Semiconducting nanowires (NWs) with a strong spin–orbit (SO) coupling in the proximity of s-wave superconductors and in the presence of an external Zeeman magnetic field B are a promising platform to study Majorana physics. Theory predicts that above a critical field $B_{\mathrm{ c}}\equiv \sqrt {\mu ^2+\Delta ^2}$, defined in terms of the Fermi energy μ and the induced s-wave pairing Δ, the wire undergoes a topological transition into a phase hosting zero-energy Majorana bound states (MBSs) at the ends of the wire [1, 2]. Recent experiments have reported measurements of differential conductance dI/dV that support the existence of such MBSs at normal–superconductor (NS) junctions in InSb [3, 4] and InAs [5] NWs. The main result of these experiments is an emergent zero-bias anomaly (ZBA) in dI/dV as B increases. In this context, the ZBA results from tunnelling into the MBS [6–8]. Although these experiments are partially consistent with the MBS interpretation [9–12], some important features such as the expected superconducting gap inversion were not observed. Moreover, other mechanisms that give rise to ZBAs, such as disorder [12–15], Kondo physics [16] or Andreev bound states (ABSs) [17, 18], cannot be completely ruled out. In particular, Zeeman-resolved ABSs in NWs with charging effects [18] can give magnetic field dependences essentially indistinguishable from some of the claimed Majorana experiments [5]. Furthermore, it has been recently pointed out that even ZBAs similar to the 0.7 anomaly in quantum point contacts may play a role in single-barrier structures [19].

Stronger evidence could be provided by the observation of non-Abelian interference (braiding) [20], or by transport in phase-sensitive superconductor–normal–superconductor (SNS) junctions. The latter approach, which typically involves the measurement of an anomalous 'fractional' 4π-periodic ac Josephson effect [21–23], is much less demanding than performing braiding. Realistically, however, the fractional effect, detected through, for instance, the absence of odd steps in Shapiro experiments [23–26], may be difficult to measure (dissipation is expected to destroy it in the steady state), or may even develop without relation to topology [27]. Although it has been shown that the 4π periodicity survives in the dynamics, such as noise and transients [28–30], simpler experimental probes of MBSs are extremely desirable.

Here we propose the multiple Andreev reflection (MAR) current in voltage-biased SNS junctions made of NWs [31, 32] as an alternative, remarkably powerful, yet simple tool to study the topological transition. This is made possible by the direct effect that gap inversion, MBS formation and fermion-parity conservation have on the MAR current I_dc(V ) at various junction transparencies T_N. For tunnel junctions, I_dc(V ) traces the closing and reopening of the superconducting gap at B_c, Δ_eff ∼ |B − B_c|. This gap inversion can be shown to be a true topological transition by tuning the junction to perfect transparency T_N = 1. In this regime, the limiting current I_dc(V →0) shows signatures of parity conservation effects that are responsible for the fractional Josephson current in the presence of MBSs, but which, in contrast to the latter, survive in the steady-state limit. Moreover, the detailed dependence of MAR as a function of T_N has a fundamental advantage over NS junctions in that it contains information about the peculiar dependence of MBS hybridization with superconductor phase difference ϕ, despite not requiring any external control on it. Similarly, we show that another important phase-insensitive quantity, the critical current I_c, remains unexpectedly finite for all B due to a significant continuum contribution, and exhibits an anomaly at the topological transition.

This paper is organized as follows. In section 2, we review the model for Rashba NWs in the presence of both s-wave superconducting pairing and an external Zeeman field and describe how a spinless p-wave superconductor regime can be achieved. In particular, we discuss how the problem can be understood in terms of two independent p-wave superconductors, originated from the Rashba helical bands, and weakly coupled by an interband pairing term. This two-band description is very useful in order to understand the main results of this paper. Such results are discussed in section 3 which is divided into two parts. The first part (subsection 3.1) is devoted to the ABSs which are confined in the junction. The detailed evolution of these ABSs as the system undergoes a topological transition has not been discussed in the literature, to the best of our knowledge, and becomes essential in order to gain a deep understanding of transport across the junction, which is discussed in section 3.2. In this part we study the ac Josephson effect in NW SNS junctions, with a focus on how the MAR currents reflect the topological transition in the NWs as the Zeeman field increases. In particular, we present a thorough analysis of MAR transport in topological SNS junctions for arbitrary transparency of the normal part. Finally, we discuss in section 4 how the critical current I_c does not vanish at the critical point where the system becomes gapless and, importantly, how I_c provides direct evidence of the topological transition.

A single one-dimensional NW in the normal state is described by the Hamiltonian

$$\begin{eqnarray*} &&H_0=\frac{p^2}{2m^*}-\mu+\alpha_{\mathrm{so}}\sigma_y p+B\sigma_x, \end{eqnarray*}$$

where m* is the effective mass, α_so is the SO coupling, μ is the Fermi energy and σ_i is the spin Pauli matrices. An external magnetic field $\mathcal {B}$ along the wire produces a Zeeman splitting $B=\frac {1}{2}g\mu _{\mathrm{ B}} \mathcal {B}$, where μ_B is the Bohr magneton and g is the wire g-factor. The Nambu Hamiltonian

$$\begin{equation} H=\left[\begin{array}{@{}cc@{}} H_0 & -\mathrm{i}\Delta\sigma_y\\ \mathrm{i}\Delta^*\sigma_y & -H^*_0\\ \end{array}\right] \end{equation} \tag{ 1 }$$

models the NW in the presence of an induced s-wave superconducting pairing Δ (here assumed real without loss of generality). The essential ingredient for a topological superconductor is an effective p-wave pairing acting on a single ('spinless') fermionic species [21]. SO coupling splits NW states into two subbands of opposite helicity at B = 0. At finite B, these two subbands, which we label + and − (black and orange lines in figure 1(a)), have spins canted away from the SO axis. The s-wave pairing Δ, expressed in the ± basis, takes the form of an intraband p-wave $\Delta ^{++/--}_{\mathrm{p}}(p)= \pm \mathrm{i} p\Delta \alpha _{\mathrm{so}}/\sqrt {B^2+(\alpha _{\mathrm{so}} p)^2}$, plus an interband s-wave pairing $\Delta ^{+-}_{\mathrm{s}}(p)=\Delta B/\sqrt {B^2+(\alpha _{\mathrm{so}} p)^2}$ [33]. Without the latter, the problem decouples into two independent p-wave superconductors, while Δ⁺⁻_s acts as a weak coupling between them. Each quasi-independent ± sector has a different (B-dependent) gap, which we call Δ₋ (at small p) and Δ₊ (large p), see figures 1(a) and (b). While Δ₊ remains roughly constant with B (for strong SO coupling [9, 11]), Δ₋ vanishes linearly as B approaches the critical field, Δ₋ ≈ |B − B_c|.⁴ This closing and reopening (gap inversion) signals a topological transition, induced by the effective removal of the 'minus' sector away from the low-energy problem. Below B_c the NW is composed of two spinless p-wave superconductors, and is therefore topologically trivial. Above B_c, Δ₋ is no longer a p-wave gap, but rather a normal (Zeeman) spectral gap already present in the normal state, transforming the wire into a single-species p-wave superconductor with non-trivial topology. This phase contains MBSs, protected by the effective gap Δ_eff = Min(Δ₊,Δ₋), at the wire ends. Above a certain field B⁽²⁾_c, the gap Δ_eff saturates at Δ₊ and the physics of superconducting helical edge states in spin-Hall insulators is recovered [22, 28].

Zoom In Zoom Out Reset image size

In the previous section, we described how a semiconducting NW with a strong SO coupling in the proximity of an s-wave superconductor and in the presence of an external Zeeman magnetic field B behaves as a topological superconductor above a critical field B_c. Here we are concerned with the effects of this topological transition on the MAR current I_dc(V ) across junctions formed with such NWs. In particular, we consider SNS junctions of different normal transparencies T_N. Experimentally, such geometry can be fabricated by partially covering a single NW with two superconducting leads and leaving an uncovered normal region in the middle. The coupling of the normal part of the NW to the superconducting leads can be tuned by local control of the electron density with a gated constriction. This can be realized by using, for instance, bottom-gates forming a quantum point contact, see figure 2. Such geometry has been successfully implemented experimentally in [19] for NS junctions, where control of the coupling between the superconducting and normal sections from near pinch-off (tunnelling limit) to the multichannel regime is demonstrated. For simplicity, here we focus on short SNS junctions⁵ with single channel NWs. For computation purposes, we consider a discretization of the continuum model equation (1) for the Rashba NW into a tight-binding lattice with a small lattice spacing a. This transforms terms containing the momentum operator p into nearest-neighbour hopping matrices v. Namely $H_0=\sum _{i}c^+_ihc_i+\sum _{\langle ij\rangle }c^+_i vc_j + \mathrm{h.c.}$, with

$$\begin{eqnarray*} h=\left(\begin{array}{@{}cc@{}} 2t-\mu & B \\ B & 2t-\mu \end{array}\right),\quad v=\left(\begin{array}{@{}cc@{}} -t & \frac{\hbar}{2a}\alpha_{\mathrm{so}} \\ -\frac{\hbar}{2a}\alpha_{\mathrm{so}} & -t \end{array}\right) \end{eqnarray*}$$

are matrices in spin space and t = ℏ²/2m*a². The pairing is incorporated like in equation (1). A short SNS junction is modelled by suppressing the hopping matrix v₀ = νv between two sites in the middle of the wire, which represent the junction. The dimensionless factor ν∈[0,1] controls the junction's normal transparency at B = 0, which we denote as T_N(ν). A phase difference ϕ across the junction is implemented by multiplying Δ to the left and right of the junction by e^∓iϕ/2, respectively. Despite the simplicity of this description, it contains the relevant physics of a short SNS junction. As has been shown for standard junctions [36], such physics essentially depend on the contact normal transmission as well as the voltage drop across it⁶. Thus, we expect that a more detailed modelling, including e.g. a spatially dependent voltage drop, would only modify the effective transmission T_N(ν) which defines the different regimes we shall explain in the following.

Zoom In Zoom Out Reset image size

3.1. Andreev bound states

In such short SNS junctions, an ABS should form for each of the two p-wave sectors described in section 2 for B < B_c, while only one, associated with Δ₊, should remain for B > B_c. To support this picture, we present calculations of the local density of states (LDOS) at the junction in the transparent limit (T_N = 1). This LDOS is peaked at the energy _± of the ABS, which is a function of the phase difference ϕ across the junction. For B = 0 (figure 3(a)) the two ABSs are degenerate and confined within the gap Δ.⁷ As the Zeeman field increases, figure 3(b), the two ABSs split and the system develops the two distinct gaps Δ₊ and Δ₋ described in section 2. Note that both ABSs are truly bound at energies below the lowest gap Δ₋, but only quasibound in the energy window Δ₋ <  < Δ₊. This is readily apparent in the plots as a broadening of the ABS resonances (see e.g. figure 3(c)). As B approaches the critical field B_c, Δ₋ gets reduced and becomes exactly zero at B = B_c. Note that at this point the upper ABS has reached zero energy at ϕ = π and is quasibound for all energies, figure 3(d). Upon entering the topological phase (B ⩾ B_c), Δ₋ reopens, but one of the ABSs of the problem has disappeared (figure 3(e)). The surviving ABS associated with Δ₊ arises due to the hybridization of the two emerging MBSs across the junction. Global fermion-parity conservation protects the ϕ = π level crossing. Owing to the residual Δ⁺⁻_s coupling between the two sectors, the Δ₊ Andreev state is once more quasibound in the energy window Δ₋ <  < Δ₊. At high enough magnetic fields, Δ₊ is the smallest gap of the problem and hence the Majorana ABS is truly bound, figure 3(f). In long junctions, more ABSs can be confined in the junction [9, 35, 38]. These extra ABSs coexist with the ones described here and may, for example, anti-cross with the Majorana-like Δ₊ Andreev level, affecting its character near zero energy [9].

Zoom In Zoom Out Reset image size

3.2. Ac Josephson effect and multiple Andreev reflection currents

Under a constant voltage bias V , the pairings Δ to the left and right of the junction acquire an opposite and time-dependent phase difference, ϕ(t) = 2eV t/ℏ. This induces Landau–Zener transitions between the ABSs and into the continuum, thereby developing a time-dependent Josephson current with both I_dc and I_ac components. Such is the point of view in e.g. [29, 39]. Alternatively, ϕ(t) can be gauged away into the hopping across the junction, $v_0(t)=\nu \,\mathrm{e}^{-\mathrm{i}\frac {{eV}}{\hbar } t\tau _z}\sum _{\sigma \sigma '} c^+_{r\sigma }v_{\sigma \sigma '}c_{l\sigma '} + \mathrm{h.c.}$, where τ_z is the z-Pauli matrix in Nambu space. By employing Keldysh–Floquet theory [36, 40], we obtain the stationary-state time-dependent ac Josephson current $I(t)=\sum _n \,\mathrm{e}^{\mathrm{i}n\frac {{eV}}{\hbar } t}I_n$ (note that only even harmonics survive, see appendix A for full details). Here, we concentrate on the dc-current I_dc = I₀. The results for I_dc(V ) at small, intermediate and full transparency are summarized in figures 4(a)–(c) for increasing values of B spanning the topological transition.

Zoom In Zoom Out Reset image size

3.2.1. Tunnelling regime

For non-topological tunnel junctions, dc-transport vanishes below an abrupt threshold voltage V_t = 2Δ_eff/e = 2Δ₋/e (figure 4(a), blue curves). This well-known result follows from the fact that there are no quasiparticle excitations in the decoupled wires for energy ∈(− Δ_eff,Δ_eff) if B < B_c. Indeed, to second order in perturbation theory in ν, the MAR current takes the form of a convolution between A₀(ω) and A₀(ω ± eV /ℏ), where A₀ is the decoupled (ν = 0) spectral density at each side of the junction (appendix B). (The trace of A₀(), proportional to the local density of states (LDOS), is shown in figure 5(a).) Hence, as B increases, the tunnelling current threshold follows the closing of the gap in the LDOS, until V_t vanishes and I_dc becomes linear in small V at B_c (black curve). As B > B_c, the gap reopens, but the threshold is now halved to V_t = Δ_eff/e (red curves).⁸ The change, easily detectable as a halving of the slope of the threshold dV_t/dB across B_c, is due to the emergence of an intra-gap zero-energy MBS in the topological phase (see zero-energy peak in figure 5(a)), which opens a tunnelling transport channel from or into the new zero-energy state. Moreover, when B = B⁽²⁾_c, Δ₋ surpasses Δ₊, and Δ_eff saturates at Δ₊. This is directly visible in V_t(B) as a kink at B⁽²⁾_c (see dashed and dotted lines in figure 4(a))⁹.

Zoom In Zoom Out Reset image size

3.2.2. Intermediate transparency regime

As transparency increases, subharmonic MAR steps develop at voltages V_t/n = 2Δ_eff/en ( $n=2,3,4,\ldots$) see figure 4(b). The specific profile of each step with V still contains information on the LDOS of the junction at energies around Δ_eff. At B = 0, the power-law LDOS for || > Δ results in a staircase-like curve I_dc(V ) (blue line in figure 4(d)). This shape is roughly preserved up to B = B_c. For B > B_c the MAR profile changes qualitatively, however. The subharmonic threshold voltages V_t/n are halved (since V_t is halved), and the MAR current profile becomes oscillatory instead of step-like. A blowup of the oscillations is presented in figure 4(d) (red curve), together with guidelines for the corresponding V_t/n in grey.

The emergence of oscillatory MAR steps, which here is connected to the formation of zero-energy peaks in the LDOS owing to the localized MBSs, is well known in Josephson junctions containing a resonant level [42–44]. Note, however, that the oscillations in a topologically trivial system, such as for instance a quantum dot between two superconductors, arise at odd fractions of 2Δ_eff, i.e. at voltages 2Δ_eff/(2n − 1)e, instead of the Δ_eff/en of the Majorana case. Interestingly, this difference is ultimately due to the fact that a resonant level spatially localized within the junction cannot carry current directly into the reservoirs, while a zero-energy MBS (essentially half a non local fermion) can. This same situation arises in d-wave Josephson junctions, which also exhibit oscillatory Δ_eff/en MAR subharmonics owing to the presence of mid gap states [45].

3.2.3. Transparent limit

In the limit T_N → 1, ABS energies _±(ϕ) (figure 3(c)) touch the continuum at ϕ = 0. This has an important consequence. From the Landau–Zener point of view of the ac Josephson effect [39], the time dependence of ϕ(t) = 2eV t/ℏ for an arbitrarily small V will induce the escape of any quasiparticle occupying an ABS into the continuum after a single ϕ(t) cycle. A given ABS becomes occupied with high probability in each cycle around ϕ = π if the rate ℏ dϕ(t)/dt = 2 eV exceeds its energy minimum (π) ≡ δ_π. (Recall that this energy is finite, since the Andreev approximation does not apply, see appendix C.) One quasiparticle is then injected into the continuum per cycle, and a finite I_dc(V ≳δ_π/e) arises. Below such voltage, however, the ABS remains empty, so that if δ_π is finite, as is the case of a realistic non-topological junction (see figure 3(a)–(c)), one obtains I_dc(V →0) = 0 (valid for any transparency at B < B_c). This is in contrast to the conventional B = 0, T_N = 1 result $I(V\to 0)=4\Delta \mathcal {G}_0/e$, predicted within the Andreev approximation ( $\mathcal {G}_0=e^2/h$).

After the topological transition, this picture changes dramatically. The two MBSs at each side of the junction hybridize for a given ϕ into a single ABS. This seemingly innocent change has a notable consequence. Since fermion parity in the superconducting wires is globally preserved, an anticrossing at ϕ = π, which would represent a mixing of a state with one and zero fermions in the lone ABS, is forbidden. Parity conservation therefore imposes δ_π = 0 in the presence of MBSs, irrespective of T_N or μ/Δ.¹⁰ This is a true topologically protected property of the junction, and gives rise to a finite $I_{\mathrm{ dc}}(V\to 0)=2\Delta _{\mathrm{eff}} \mathcal {G}_0/e$, i.e. half the value expected for the non-topological junction in the Andreev approximation. This abrupt change is shown in figures 4(c) and (e).

The I_dc(V →0) MAR current in transparent junctions, therefore, directly probes the emergence of parity protection.

In the transparent limit, a supercurrent peak [31, 32, 46] may hinder the experimental identification of the I_dc(V →0) limit, but itself holds valuable information about the transition. The critical current I_c may be computed in general by maximizing the V =0 (time-independent) current I(ϕ) with respect to ϕ (including the contribution from the continuum). For a short transparent junction at B = 0, I_c is maximum and is equal to I⁰_c ≡ eΔ/ℏ in the Andreev approximation. Figure 5(b) shows I_c for increasing values of B. Naively, one may expect that a junction without a superconducting gap should not carry a finite supercurrent, but this is not the case here. At B = B_c, I_c is finite¹¹, while Δ_eff = 0 (the junction LDOS at criticality is also gapless, see figure 5(a)). This gapless supercurrent comes from the ₊(ϕ) quasi-bound Andreev state in the continuum, which contributes almost as if it was a subgap ABS. It is thus a reasonable approximation to write I_c as the sum of the critical current from each ABS. For B < B_c, $I_{\mathrm{ c}}\approx \frac {1}{2}I_{\mathrm{ c}}^0(\Delta _++\Delta _-)/\Delta$. The Δ₋ contribution, however, should not be included for B > B_c, leading to a discontinuity in ∂I_c/∂B. This simple model gives a qualitative fit (dotted line in figure 5(b)) to the exact numerics (solid line), with deviations coming from corrections to the Andreev approximation, and contributions above Δ₊. Additional deviations in experiments, coming e.g. from the finite impedance of the electromagnetic environment, are not expected to alter the discontinuity in ∂I_c/∂B, which remains a signature of the topological transition.

In conclusion, we have shown that the dc current in voltage-biased Josephson junctions is a flexible experimental probe into the various aspects of the topological superconducting transition in semiconducting NWs. Tuning the junction transparency one may obtain evidence of MBS formation as conclusive as a fractional Josephson effect, without requiring control of the junction phase. Moreover, we have found that the critical current in the wire does not vanish at the transition due to above-gap contributions, although its derivative with B exhibits a discontinuity as a result of the disappearance of one ABS. This behaviour of I_c provides direct evidence of the topological transition. MAR spectroscopy and critical current measurements in NWs similar to the ones studied here have already been reported [31, 32].

Although we have focused here on the simplest case (single-band, short junction limit), we expect the main features of the topological transition to remain robust under more general conditions. Preliminary results in the quasi-one-dimensional multiband case show that I_c is a non-monotonic function for increasing magnetic fields. For weak interband SO mixing [48], the behaviour discussed in figure 5(b) can be generalized and I_c presents a series of minima at different fields corresponding to the topological transition of each subband.

Importantly, the alternative physical scenarios, such as, for instance, disorder [12–15] or ABSs [18], that produce ZBAs in NS junctions (and thus mimic Majorana physics), cannot give the distinct features associated with global parity that were discussed here for SNS junctions. We therefore believe that experiments along the lines discussed in this paper could provide the first unambiguous report of a topological transition in NWs, and the emergence of MBSs.

We acknowledge the support of the European Research Council, the Spanish Research Council CSIC through the JAE-Predoc Program (JC) and the Spanish Ministry of Economy and Innovation through grant numbers FIS2012-33521, FIS2011-23713, FIS2010-21883 and FIS2009-08744 and the Ramón y Cajal Program (EP).

Consider a mesoscopic system composed of two semi-infinite leads (labelled L and R), each in thermal equilibrium at the same temperature T and with the same chemical potential μ = 0. Each lead has a finite s-wave superconducting pairing Δ_α, where α = L,R. A central system (α = S), which may or may not be superconducting, is coupled to both leads through operator v. In its Nambu form, the Hamiltonian of the system reads

$$\begin{eqnarray*} \skew3\hat{H}=\frac{1}{2}\sum_{ij}(\begin{array}{@{}c|c@{}} c_j & c^+_j \end{array}) H_{ij}\left(\begin{array}{@{}c@{}} c_j\\ \hline c^+_j \end{array}\right), \end{eqnarray*}$$

where the Nambu Hamiltonian matrix takes the general form

$$\begin{eqnarray*} H=\left(\begin{array}{@{}ccc|ccc@{}} h_{\rm L} & v^+ & 0 & \Delta_{\rm L} & 0 & 0 \\ v & h_{\mathrm{S}} & v^+ & 0 & \Delta_{\mathrm{S}} & 0 \\ 0 & v & h_{\rm R} & 0 & 0 & \Delta _{\rm R}\\\hline \Delta_{\rm L}^+ & 0 & 0 & -h_{\rm L}^* & (-v^+)^* & 0 \\ 0 & \Delta_{\mathrm{S}}^+ & 0 & -v^* & -h_{\mathrm{S}}^* & (-v^+)^* \\ 0 & 0 & \Delta_{\rm R}^+ & 0 & -v^* & -h_{\rm R}^{*} \end{array}\right). \end{eqnarray*}$$

Here h_α is the normal Hamiltonian for each section of the system. The blocks delimited by lines denote the Nambu particle, hole and pairing sectors.

If we apply a left–right voltage bias V through the junction, the Bardeen–Cooper– Schrieffer (BCS) pairing of the leads will become time dependent, Δ_L/R → e^{±iV t}Δ_L/R, while h_L/R → h_L/R ± V/2 (we take e = ℏ = 1). Both these changes can be gauged away from the leads and into the system by properly redefining c⁺_i → c⁺_i(t) = e^{±iV t/2}c⁺_i. This transformation is done also inside the system S, thereby effectively dividing it into two, the portion with an e^{iV t/2} phase (denoted S_L), and the portion with the opposite phase (denoted S_R). This restores H to its unbiased form, save for a new time dependence in h_S → h_S(V t), which is constrained to the coupling between the S_L and S_R,

$$\begin{eqnarray*} h_{\mathrm{S}}(Vt)=\left(\begin{array}{@{}cc@{}} h_{\mathrm{S}_{\rm L}} & \mathrm{e}^{-\mathrm{i}Vt}v^+_0\\ \mathrm{e}^{\mathrm{i}Vt}v_0 & h_{\mathrm{S}_{\rm R}} \end{array}\right). \end{eqnarray*}$$

It is important to note that H(t) is periodic, with an angular frequency ω₀ = V . In the steady state limit (at long times t after switching on the potential V ) all response functions and observables will exhibit the same time periodicity (all transient effects are assumed to be completely damped away). In particular, the steady-state current I(t) = I(t + 2π/ω₀), so that

$$\begin{eqnarray*} &&I(t)=\sum_n \mathrm{e}^{\mathrm{i}n\omega_0 t}I_n \end{eqnarray*}$$

for some harmonic amplitudes I_n, in general complex, that satisfy I_n = I*_−n since I(t) is real. This current can be computed using the Keldysh Green's function formalism [49]. The standard expression for I(t) is computed starting from the definition of I(t) = ∂_tN_L, where N_L is the total number of fermions in the left lead. By using the Heisenberg equation and the Keldysh–Dyson equation, one arrives at

$$\begin{eqnarray*} &&I(t)=\mathrm{Re}[J(t)], \end{eqnarray*}$$

where

$$\begin{eqnarray*} &&J(t)=\frac{2e}{\hbar} \int \mathrm{d}t'~\mathrm{Tr}\{[ G^{\mathrm{r}}(t,t')\Sigma_{\rm L}^<(t',t) + G^<(t,t')\Sigma_{\rm L}^{\mathrm{a}}(t',t)]\tau_z \}. \end{eqnarray*}$$

The z-Pauli matrix τ_z above acts on the Nambu particle–hole sector,

$$\begin{eqnarray*} \tau_z=\left(\begin{array}{@{}c|c@{}} \mathbbm{1} & 0\\\hline 0 & -\mathbbm{1} \end{array}\right). \end{eqnarray*}$$

The self-energy from the left lead is defined as Σ^a,<_L(t',t) = v g^a,<_L(t',t)v⁺, where g_L(t',t) = g_L(t' − t) stands for the left lead's propagator, when decoupled from the system (this propagator depends only on the time difference since the decoupled lead is time independent in this gauge). We define the Fourier transform of g as

$$\begin{eqnarray*} &&g(\omega)=\int_{-\infty}^\infty \mathrm{d}t \,\mathrm{e}^{\mathrm{i}\omega t} g(t). \end{eqnarray*}$$

The retarded propagator in Fourier space is

$$\begin{eqnarray*} &&g_{\rm L}^{\mathrm{r}}(\omega)=\frac{1}{\omega-h_{\rm L}+\mathrm{i}\eta} \end{eqnarray*}$$

while the advanced $g^{\mathrm{a}}_{\mathrm{ L}}(\omega )=[g^{\mathrm{r}}_{\mathrm{ L}}(\omega)]^+$. One can compute g^<_L(ω) = if(ω)A_L(ω), where f(ω) = 1/(e^ω/k_BT + 1) is the Fermi distribution in the leads, and A_L(ω) = i(g^r_L(ω) − g^a_L(ω)) is the Nambu spectral function. The g^r_L/R (and in particular A_L/R(ω)) is assumed known, or at least easily obtainable from h_L/R and v. Finally, the Green functions G^r(t',t) and G^<(t',t) correspond to the propagator for the full system, including the coupling to the leads. (Note that, in practice, since G is inside a trace in J(t), only matrix elements of G inside the S portion of the full system are needed.) The retarded G^r satisfies the equation of motion

$$\begin{eqnarray*} &&\left[\mathrm{i}\partial_{t'}-H(t')\right]G(t',t)=\delta(t'-t) \end{eqnarray*}$$

while G^< (when projected onto the finite-dimensional system S) satisfies the Keldysh relation

$$\begin{eqnarray*} G^<(t',t)&=&\int \mathrm{d}t_1\,\mathrm{d}t_2 G^{\mathrm{r}}(t',t_1)\\ &&\times[\Sigma_{\rm L}^<(t_1-t_2)+\Sigma_{\rm R}^<(t_1-t_2)]G^{\mathrm{a}}(t_2,t). \end{eqnarray*}$$

Since h_S in H is time dependent, G propagators depend on two times; unlike Σ_L/R or g_L/R they are not Fourier diagonal. Instead, we can exploit the steady-state condition, which reads

$$\begin{eqnarray*} &&G(t',t)=G\left(t'+\frac{2\pi}{\omega_0},t+\frac{2\pi}{\omega_0}\right) \end{eqnarray*}$$

to expand the system's G as a Fourier transform in t' − t and a Fourier series in t. We define

$$\begin{eqnarray*} &&G(t',t)=\sum_n \mathrm{e}^{-\mathrm{i}n\omega_0 t}\int_{-\infty}^\infty \frac{\mathrm{d}\epsilon}{2\pi}\mathrm{e}^{-\mathrm{i} \epsilon(t'-t)}G_{n}(\epsilon). \end{eqnarray*}$$

The natural question is how the equation of motion is expressed in terms of the harmonics G_n(). It takes the most convenient form if we redefine G_n() (where is unbounded) in terms of the quasienergy $\tilde {\epsilon }\in [0,\hbar \omega _0]$, i.e. $\epsilon =\tilde \epsilon +m\omega _0$:

$$\begin{eqnarray*} &&G_{mn}(\tilde\epsilon)=G_{m-n}(\tilde\epsilon+m\omega_0). \end{eqnarray*}$$

This has the advantage that the equation of motion translates to a matrix equation analogous to that of a static system in Fourier space

$$\begin{eqnarray*} &&\sum_m (\tilde \epsilon+n'\omega_0-H_{n'm})G^{\mathrm{r}}_{mn}(\tilde\epsilon)=\delta_{n'n}, \end{eqnarray*}$$

where

$$\begin{eqnarray*} &&H_{n'n}=\int \mathrm{d}t \,\mathrm{e}^{\mathrm{i}(n'-n)t} H(t). \end{eqnarray*}$$

This is known as the Floquet description of the steady-state dynamics in terms of sidebands, which appear formally as a new quantum number n. Time-dependent portions of H(t) act as a coupling between different sidebands. The effective Hamiltonian for the nth sideband is the static portion of H(t), shifted by −nω₀. One therefore sometimes defines the Floquet 'Hamiltonian' of the system as

$$\begin{eqnarray*} &&\boldsymbol{h}_{{\mathrm{S}}_{nm}}=h_{{\mathrm{S}}_{nm}}-n\omega_0\delta_{nm}, \end{eqnarray*}$$

where, as before, ${h_{\mathrm{S}}}_{n'n}=\int \mathrm{d}t \,\mathrm{e}^{\mathrm{i}(n'-n)t} {h_{\mathrm{S}}}(t)$. Likewise, one may define the Floquet self-energies as

$$\begin{eqnarray*} &&{\boldsymbol{\Sigma}}_{{\rm L}_{nm}}(\tilde\epsilon)=\delta_{nm}\Sigma_{\rm L/R}(\tilde\epsilon+n\omega_0), \end{eqnarray*}$$

(since the leads are static, Σ is sideband-diagonal).

The Floquet equation of motion for $G^{\mathrm{r}}_{nm}(\tilde \epsilon )$ can be solved like in the case of a static system. Within the S portion of the system, we have

$$\begin{eqnarray*} &&\boldsymbol{G}^{\mathrm{r}}(\tilde\epsilon)=[\tilde\epsilon-{\boldsymbol{h}}_{\mathrm{S}}-{\boldsymbol{\Sigma}}_{\rm L}^{\mathrm{r}}(\tilde \epsilon)-{\boldsymbol{\Sigma}}_{\rm R}^{\mathrm{r}}(\tilde \epsilon)]^{-1}. \end{eqnarray*}$$

Boldface denotes the sideband structure implicit in all the above matrices. Similarly, the Keldysh relation takes the simple form

$$\begin{eqnarray*} &&\boldsymbol{G}^{<}(\tilde \epsilon)=\boldsymbol{G}^{\mathrm{r}}(\tilde \epsilon)[{\boldsymbol{\Sigma}}_{\rm L}^<(\tilde \epsilon)+{\boldsymbol{\Sigma}}_{\rm R}^<(\tilde \epsilon)]\boldsymbol{G}^{\mathrm{a}}(\tilde \epsilon). \end{eqnarray*}$$

Finally, the time averaged current I_dc ≡ I₀ takes the form

$$\begin{equation} I_{\rm dc}=\frac{2e}{h} \int_0^{\hbar\omega_0} \mathrm{d}\tilde\epsilon~\mathrm{Re}\,\mathrm{Tr}\{[\boldsymbol{G}^{\mathrm{r}}(\tilde\epsilon){\boldsymbol{\Sigma}}_{\rm L}^<(\tilde\epsilon)+{\boldsymbol{G}}^{<}(\tilde\epsilon){\boldsymbol{\Sigma}}_{\rm L}^{\mathrm{a}}(\tilde\epsilon)]\tau_z \}, \end{equation} \tag{ A.1 }$$

where the trace includes the sideband index. In a practical computation, the number of sidebands that must be employed is finite, and depends on the applied voltage bias V (the typical number scales as n_max ∼ v₀/V ). We employ an adaptive scheme that increases the number of sidebands recursively until convergence for each value of V .

It is possible to solve the I_dc current explicitly in the tunnelling limit. To a leading (second) order in the left–right coupling v₀, equation (A.1) reduces, after some algebra, to

$$\begin{equation} I_{\rm dc}\approx\frac{e}{\pi} \mathrm{Re}\int \mathrm{d}\omega \left[f(\omega-\omega_0)-f(\omega)\right] \mathrm{Tr}\left\{A_{11}^{\mathrm{L}}(\omega) v_0^+ A_{11}^{\mathrm{R}}(\omega-\omega_0)v_0\right\}, \end{equation} \tag{ B.1 }$$

where A^α₁₁ is the particle–particle Nambu 2 × 2 matrix block of the spectral function of the α = L,R decoupled wire,

$$\begin{eqnarray*} A^{\alpha}(\omega)=\left(\begin{array}{@{}cc@{}} A_{11}^\alpha(\omega) & A_{12}^\alpha(\omega)\\ \left[A_{12}^\alpha(\omega)\right]^+ & -\left[A_{11}^{\alpha}(-\omega)\right]^* \end{array}\right), \end{eqnarray*}$$

and the trace is taken over spin space. The trace of A^α(ω) is proportional to the LDOS. Figure B.1 shows results for the tunnel current using equation (B.1) for increasing Zeeman fields. Overall, the agreement with full numerics in figure 4(a) is very good and, importantly, all the relevant features such as, for instance, the closing of the gap, are captured by this tunnelling approximation.

Zoom In Zoom Out Reset image size

It is conventional, in the study of hybrid superconducting-normal junctions, to assume the limit in which the Fermi energy μ of the metal under consideration is much greater than the superconducting gap, and any other energy E involved in the problem, μ ≫ Δ,E. This is known as the Andreev approximation. In essence, it allows one to regard the normal system as featureless, with constant Fermi velocity and density of states. In this case, a number of simplifications can be carried out in the computation of equilibrium transport properties. One important consequence of the approximation in the context of our work is that, in a short SNS junction with phase difference ϕ, and symmetric under time-reversal symmetry (B = 0 in our case), two degenerate Andreev states will appear of energy $\epsilon (\phi )=|\Delta |\sqrt {1-T_{\mathrm{ N}}^2\sin ^2(\phi /2)}$ [37]. This immediately implies that at perfect normal transparency T_N = 1, the ABSs will reach zero energy at ϕ = π. In other words, in the Andreev approximation, the ABS energy minimum in the non-topological phase will be (π) = δ_π = 0.

While in the topological phase, a zero δ_π is a robust property, protected by parity conservation, a δ_π = 0 in the non-topological phase is accidental, a direct consequence of the Andreev approximation, and is not protected by any symmetry. In fact, any deviation from the Andreev approximation will induce a finite splitting δ_π. In semiconducting NWs such as the ones considered in this work, this correction is very relevant. Indeed, for the NW to undergo a superconducting topological transition at a reasonable Zeeman field B, μ/Δ must remain relatively small (the wire must be close to depletion, and far from the Andreev approximation), otherwise the critical field $B_{\mathrm{ c}}=\sqrt {\mu ^2+\Delta ^2}$ would be physically unreachable. The splitting δ_π, therefore, remains a relevant quantity in the formation and detection of MBSs.

The value of δ_π in our system may be computed numerically. The most efficient way is to consider a short SNS junction with finite length superconductors, T_N = 1, B = 0 and a phase difference ϕ = π. Since this system is closed, an exact diagonalization of the tight-binding Nambu Hamiltonian yields a minimum eigenvalue that is exactly δ_π if the S leads are long enough (longer than the coherence length). We find that this quantity is finite in the case of wires close to depletion, μ ≳ Δ, and that it vanishes as one approaches the Andreev approximation regime μ ≫ Δ, see figure C.1. More specifically, we have found that δ_π scales as δ_π = c₁Δ²/(μ + Δc₂) for some c_1,2 > 0, within very good precision.

Zoom In Zoom Out Reset image size