Supercurrent Interference in Semiconductor Nanowire Josephson Junctions

TL;DR

This paper models quantum transport in semiconductor nanowire Josephson junctions, revealing that critical current oscillations are interference effects of transverse sub-bands, influenced by disorder and phase-breaking processes, aligning with experimental observations.

Contribution

It introduces a detailed theoretical model of III-V semiconductor nanowire Josephson junctions, incorporating disorder and dephasing to explain critical current oscillations observed experimentally.

Findings

Critical current oscillations are interference effects of transverse sub-bands.

Disorder and gate voltage significantly affect oscillation patterns.

Phase-breaking processes are essential to match experimental data.

Abstract

Semiconductor-superconductor hybrid systems provide a promising platform for hosting unpaired Majorana fermions towards the realisation of fault-tolerant topological quantum computing. In this study, we employ the Keldysh Non-Equilibrium Green's function formalism to model quantum transport in normal-superconductor junctions. We analyze III-V semiconductor nanowire Josephson junctions (InAs/Nb) using a three-dimensional discrete lattice model described by the Bogolubov-de Gennes Hamiltonian in the tight-binding approximation, and compute the Andreev bound state spectrum and current-phase relations. Recent experiments [Zuo et al., Phys. Rev. Lett. 119,187704 (2017)] and [Gharavi et al., arXiv:1405.7455v2 (2014)] reveal critical current oscillations in these devices, and our simulations confirm these to be an interference effect of the transverse sub-bands in the nanowire. We add disorder to model coherent scattering and study its effect on the critical current oscillations, with an aim to gain a thorough understanding of the experiments. The oscillations in the disordered junction are highly sensitive to the particular realisation of the random disorder potential, and to the gate voltage. A macroscopic current measurement thus gives us information about the microscopic profile of the junction. Finally, we study dephasing in the channel by including elastic phase-breaking interactions. The oscillations thus obtained are in good qualitative agreement with the experimental data, and this signifies the essential role of phase-breaking processes in III-V semiconductor nanowire Josephson junctions.

Tables

Table 1. Table 1: Phase coherence length as a function of dephasing strength D 0 subscript 𝐷 0 D_{0}

$D_0$ (eV2)	$l_{\phi }$ (nm)
$1\times {10}^{-4}$	247
$5\times {10}^{-4}$	157
$1\times {10}^{-3}$	105

Table 2. Table 2: Phase coherence length as a function of dephasing strengths D 0 subscript 𝐷 0 D_{0}

$D_0$ (eV2)	rms($G$) $\left(2e^2/h\right)$	$B_c$ $\left(T\right)$	$l_{\phi }$ (nm)
$1\times {10}^{-4}$	0.038	0.56	247
$5\times {10}^{-4}$	0.0048	0.88	157
$1\times {10}^{-3}$	0.00098	1.32	105

Equations

$${\mathcal{H}}_0$$

$${\mathcal{H}}_p$$

$$\Delta (z)=\vartheta (-z){\Delta }_0{e}^{i{\chi }_L}+\vartheta (z-L){\Delta }_0{e}^{i{\chi }_R}$$

$$\left[\begin{array}{cc}{\mathcal{H}}_0& \Delta (z)\\ {\Delta }^{\ast }(z)& -{\mathcal{H}}_0^{\ast }\end{array}\right]\left[\begin{array}{c}u(z)\\ v(z)\end{array}\right]=E\left[\begin{array}{c}u(z)\\ v(z)\end{array}\right]$$

$${\alpha }_N=\left[\begin{array}{cc}2t-\mu & 0\\ 0& -2t+\mu \end{array}\right]$$

$${\alpha }_S=\left[\begin{array}{cc}2t-\mu & \Delta \\ {\Delta }^{\ast }& -2t+\mu \end{array}\right]$$

$$\beta =\left[\begin{array}{cc}-t& 0\\ 0& t\end{array}\right]$$

$$t\gg \mu \gg {\Delta }_0$$

$$H=\sum _i^n{c}_i^{†}{\alpha }_{N/S}{c}_i+\sum _{|i-j|=1}^n{c}_i^{†}\beta c_j$$

$${H}=\matrixquantity(\alpha_{N}&\beta&0&\dots&0\\ \beta^{\dagger}&\alpha_{N}&\beta&0&0\\ 0&\beta^{\dagger}&\alpha_{N}&\beta&\vdots\\ \vdots&0&\ddots&\ddots&\beta\\ 0&\dots&\dots&\beta^{\dagger}&\alpha_{N})$$

$$G^r(E)={\left(E\mathbb{I}+i\eta -\mathcal{H}-{\Sigma }_1^r-{\Sigma }_2^r\right)}^{-1}$$

$$d(E)=\frac{1}{2\pi }\text{Tr}\left[A(E)\right]=\frac{1}{2\pi }\text{Tr}\left[i\left(G^r(E)-G^a(E)\right)\right]$$

$$J_i(E)=\frac{2e}{h}f(E)\text{Tr}\left[\Re \left(G^a(E){\Sigma }_i^a(E)-G^r(E){\Sigma }_i^r(E)\right){\tau }_z\right]$$

$$I(\varphi )={\int }_{-\infty }^{\infty }{J}_i(E)dE$$

$${\mathcal{H}}_0=-\mu +\frac{-{\hslash }^2}{2m^{\ast }}\frac{{\partial }^2}{\partial z^2}+{\mathcal{H}}_T+U$$

$$[{\mathcal{H}}_T,{\hat{L}}_z]=0$$

$${\hat{L}}_z|\ell ⟩=\hslash \ell |\ell ⟩$$

$$\mathbf{A}=B\cdot x\hat{\mathbf{y}}$$

$${\alpha }_{N/S}=\left[\begin{array}{cc}h& {\Delta }_{N/S}\\ {\Delta }_{N/S}^{\ast }& -h^{\ast }\end{array}\right]$$

$${\beta }_{x,z}=\left[\begin{array}{cc}-t_{x,z}& 0\\ 0& t_{x,z}\end{array}\right]$$

$${\beta }_y=\left[\begin{array}{cc}-t_y{e}^{i2\pi n_x{\Phi }_a}& 0\\ 0& t_y{e}^{-i2\pi n_x{\Phi }_a}\end{array}\right]$$

$${\xi }_{\ell }=\frac{\hslash v_{F,\ell }}{2{\Delta }_0}$$

$$v_{F,\ell }=\sqrt{2\left(\mu -\frac{{\hslash }^2}{2m^{\ast }{R}^2}\left({\ell }^2+{\Phi }^2\right)\right)/m^{\ast }}$$

$$⟨{\hat{L}}_z⟩=\frac{\text{Tr}\left[G^n(E)\cdot {\hat{L}}_z\right]}{{\max }_E\{\text{Tr}\left[G^n(E)\right]\}}$$

$${\hat{L}}_z=\hat{x}\cdot {\hat{p}}_y-\hat{y}\cdot {\hat{p}}_x$$

$${\mathbb{I}}_x\otimes \hat{y}=a\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 2& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 2\end{array}\right)\hat{x}\otimes {\mathbb{I}}_y=a\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 2& 0\\ 0& 0& 0& 2\end{array}\right)$$

$${\hat{p}}_i=-i\frac{m}{\hslash }[{\hat{x}}_i,{\mathcal{H}}_0]$$

$$\left[\begin{array}{cc}{\mathcal{H}}_0& \Delta \\ {\Delta }^{\ast }& -{\mathcal{H}}_0^{\ast }\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]=E\left[\begin{array}{c}u\\ v\end{array}\right]$$

$${\mathcal{H}}_0=-\frac{{\hslash }^2}{2m^{\ast }}\frac{{\partial }^2}{\partial z^2}+\frac{1}{2m^{\ast }}{\left(-i\hslash \frac{1}{R}\frac{\partial }{\partial \theta }-e{A}_{\theta }\right)}^2-\mu$$

$$h_{\ell }=\frac{{\hslash }^2{k}_z^2}{2m^{\ast }}+\frac{{\hslash }^2{\ell }^2}{2m^{\ast }{R}^2}-\mu =\frac{{\hslash }^2{k}_z^2}{2m^{\ast }}-{\mu }_{\ell }$$

$${\mu }_{\ell }=\mu -\frac{{\hslash }^2{\ell }^2}{2m^{\ast }{R}^2}$$

Full text

I Introduction

Semiconductor-superconductor hybrid junctions have generated significant interest over the last decade. In particular, III-V semiconductor (InAs/InSb) nanowires in proximity to an s-wave superconductor have been extensively studied as a platform for topological superconductivityZhang et al. (2018); Mourik et al. (2012); Deng et al. (2016); Chen et al. (2017); Das et al. (2012); Rokhinson et al. (2012). Majorana bound states (MBSs) emerge as zero energy edge excitations in a gapped bulk spectrum of the topological superconducting nanowire Nayak et al. (2008); Majorana and Maiani (2006); Kitaev (2001); Beenakker (2013); Sau et al. (2010a, b); Alicea (2010, 2012); Lutchyn et al. (2010). Signatures of MBS have been reported as a zero-bias conductance peak in tunnelling experiments Zhang et al. (2018); Mourik et al. (2012); Deng et al. (2016); Chen et al. (2017); Das et al. (2012); Wimmer et al. (2011). The $4\pi$ Majorana-Josephson effect has been predicted and observed in nanowire Josephson junctions tuned to the topologically non-trivial regime Kitaev (2001); Rokhinson et al. (2012). With the massive progress being made with nanowire setups, it is anticipated that the focus will shift from the detection to the demonstration of non-Abelian statistics and finally to topological quantum information processingAlicea et al. (2011); Sarma et al. (2015); Aasen et al. (2016); Hyart et al. (2013); Karzig et al. (2017); Plugge et al. (2017). The 4$\pi$ Majorana-Josephson effect forms the basis of braiding and readout schemes of a recent topological qubit proposalStenger et al. (2019).

Josephson junctions based on semiconductor-superconductor hybrids form the basis for microwave quantum circuityLarsen et al. (2015), and superconducting qubitsde Lange et al. (2015); Hassler et al. (2011). They afford an attractive alternative for a scalable computing architecture with the possibility of an all-electric qubit controlKringhøj et al. (2018); Larsen et al. (2015); Casparis et al. (2018).

Several studies have focused on the structure of transverse subbands and magnetoconductance due to radial confinement in semiconductor nanowires Estévez Hernández et al. (2010a); Blömers et al. (2011a); Lahiri et al. (2018); Cayao et al. (2015) and carbon nanotubesLim et al. (2011). Recent experiments study the critical current as a function of the magnetic field and gate voltage in nanowire Josephson junctions tuned to the few-subband regimeGharavi et al. (2014); Zuo et al. (2017); Szombati (2017). For a magnetic field oriented along the nanowire axis, Zuo et al. measured a strong suppression of the critical current at fields on the order of 100 mT in InSb weaklinks with NbTiN contacts. At higher fields, the critical current exhibited local minima (nodes). Similar results were obtained by Gharavi et al. for InAs-Nb Josephson junctions. Unlike the Fraunhofer diffraction in wide planar junctions, the critical current nodes were aperiodic in the magnetic field, and highly sensitive to local fluctuations in the gate voltage. Motivated by these experiments, the object of this paper is to theoretically analyze few-mode nanowire Josephson junctions in a magnetic field oriented along the nanowire axis. We thus employ the Keldysh Non-Equilibrium Green’s Function formalism (NEGF)Keldysh (1964); Datta (1995, 2005); DuBois (1967); Rammer and Smith (1986); Zeng et al. (2003) to model quasiparticle transport in the junction, and compute the evolution of the critical current as a function of the axial field and chemical potential. Based on the simulations, we attribute the observed oscillations to the interference of the transverse subbands in the nanowire. These results are crucial in the design of Majorana setupsAasen et al. (2016); Hyart et al. (2013); Karzig et al. (2017); Plugge et al. (2017) and in interpreting experiments, particularly for those based on critical current measurementsCayao et al. (2017); San-Jose et al. (2013); Cayao et al. (2018); San-Jose et al. (2014).

Quantum transport traditionally involves excited states and the use of a variant of the Landauer-Büttiker’s scattering theoryDatta (1995); Beenakker and van Houten (1991); Singha et al. (2017) for performing transport calculations. This essentially involves solving the Schrödinger equation and an appropriate treatment of the boundary conditions. In a superconductor, however, zero-bias transport is essentially a ground state phenomenon supported by Cooper pairs condensed at the fermi levelDe Gennes (1999); Tinkham (2004). Blonder et al. generalised the scattering theory approach to hybrid semiconductor-superconductor junctions by solving the Bogoliubov-de Gennes equation across the N-S interfaceBlonder et al. (1982). Beenakker applied this formalism for mesoscopic N-S junctions, thus providing a multichannel generalization of Blonder’s resultsBeenakker (1992a). This technique has been prevalent in the literatureGharavi and Baugh (2015); Bagwell (1992) ever since, and it forms the basis for numerous simulation packages such as Kwant Groth et al. (2014). Despite its benefits, the scattering theory approach is not very convenient in dealing with disordered junctions. While phase-coherent scattering processes can be included via random on-site potentials, it is difficult to model phase-relaxing interactions. Moreover, this formalism becomes intractable whenever a self-consistent determination of the order parameter becomes necessary. This self-consistent computation can be performed using the correlation Green’s functionLevy Yeyati et al. (1995); Martín-Rodero et al. (1994), and various dephasing mechanisms such as electron-electron and electron-phonon interactions can be included through suitable self-energy operators in the NEGF formalism. The compatibility with phase-breaking processes is one of the main advantages of NEGF over scattering theory. Our results indicate dephasing to be essential in achieving qualitative agreement with the experiment, and this is one of the key takeaways of this paper.

This paper is organised as follows. We start with the Bogoliubov-de Gennes mean-field description of a one-dimensional nanowire Josephson junction (SNS). In Sec. II, we describe the junction in a tight-binding model and outline the key aspects of the NEGF formalism. The details of this formalism have been relegated to Appendix A. Employing the NEGF formalism, we compute the Andreev bound state (ABS) spectrum and current phase relationship (CPR) for this junction. Previous work almost exclusively focused on the Andreev approximation regime, which assumes the chemical potential of the nanowire ($\mu$) to be much larger than the superconducting order parameter ($\Delta_{0}$), i.e., $\mu\gg\Delta_{0}$Beenakker (1992b); Bagwell (1992); Kulik (1969); Blonder et al. (1982). We go beyond this Andreev approximation limit and investigate the bound states which anti-cross at a superconducting phase difference of $\pi$ between the leads. This anti-crossing in the ABS is further analyzed in Appendix B. In Sec. III, we model three-dimensional Josephson junctions in a discrete lattice tight-binding model (Fig. 1). The radial confinement gives rise to transverse angular momentum subbands which pick up characteristic phases in a magnetic field. Section III.2 details the procedure we follow to label these angular momentum subbands. In Sec. IV, we reproduce the critical supercurrent oscillations in the presence of an axial magnetic field. Our results confirm these observed oscillations to be arising from the interference between orbital channels of the junction. With the aim of gaining a thorough understanding of the experiments, we consider scattering processes in the nanowire and study the effect of disorder, gate voltage fluctuations, and phase-breaking processes on the critical current oscillations.

II Formalism

Superconducting correlations are induced in a proximitised semiconductor by electron-hole conversions at the interface, a process known as Andreev reflectionAndreev (1964, 1967). Low bias transport in normal (N)-superconductor (S) junctions involves Andreev reflections at the interface.

We first consider a one-dimensional SNS junction consisting of a semiconductor nanowire with supercondcuting contacts. We model this system using the Bogoliubov-de Gennes (BdG) mean-field Hamiltonian within the tight-binding approximation, $H=\mathcal{H}_{0}+\mathcal{H}_{\text{p}}$, where

[TABLE]

$\mathcal{H}_{0}$ is the single-particle effective Hamiltonian, $\psi_{\sigma}$ is the field operator with spin index $\sigma\in\{\uparrow,\downarrow\}$, $m^{*}$ is the electron effective mass, and $V$ models a potential energy induced in the junction. The chemical potential is defined as the energy difference between the lowest occupied subband and the Fermi energy, and is denoted by $\mu$. We assume an identical effective mass in the N and S regions thus neglecting the Fermi wave-vector mismatch at the interface. $\Delta(z)$ is the superconducting order parameter along the junction, which we assume to be constant with jump-discontinuities at the N/S interfaces

[TABLE]

where $\vartheta(x)$ is the unit step function at $x=0$, $\chi_{L,R}$ is the superconducting phase of the left and right leads respectively, and $\phi=\chi_{L}-\chi_{R}$ is the phase difference. In the $\left[\psi_{\uparrow}^{\dagger}(r),\psi_{\downarrow}(r)\right]$ Nambu basis, we have the BdG equation

[TABLE]

The device is divided into three parts – a normal semiconductor section with a length $L$ extended over $z\in[0,L]$, and semi-infinite superconducting contacts extending to $z=\pm\infty$ on either side of the semiconductor (Fig. 2).

We discretise the continuum model of Eqs. 1, 2 into a lattice model with a spacing of $a$. This is shown in Fig. 2. The superconductors are modelled as semi-infinite leads, while the number of lattice points in the normal region controls the length of the nanowire. The on-site tight-binding parameters in the normal and superconducting regions, in the Nambu representation are

[TABLE]

[TABLE]

where, $t={\hbar^{2}}/\left({2m^{*}a^{2}}\right)$ is the nearest neighbour tight binding hopping parameter. The hopping matrix is given by

[TABLE]

This is not an atomic model, but an “effective” discrete lattice description of the junction. Parabolic dispersion relations correspond to the parameter space $\mu\ll t$. Typical experiments are setup in the $\mu\gg\Delta_{0}$ (Andreev approximation) regime. Hence, the choice of the effective lattice parameter $a$ is bound by the inequalities

[TABLE]

The device Hamiltonian is subsequently written as

[TABLE]

where $c^{\dagger}_{i}$ is the creation operator of the Nambu spinor $\left[\psi_{\uparrow}^{\dagger}(r),\psi_{\downarrow}(r)\right]$ at site $i$, and $n=L/a$ is the number of sites in the device. The Hamiltonian of the normal region can be written in the general form

[TABLE]

II.1 Andreev bound states in SNS junctions

Andreev reflections at the N/S interfaces give rise to Andreev bound states in the semiconductor. We use the NEGF formalism to compute these bound state energies as a function of the superconducting phase difference ($\phi$) of the leads. The retarded Green’s function in the energy domain is given by

[TABLE]

where $E$ denotes the energy, $\mathbb{I}$ is the identity matrix and $\eta$ is an infinitesimal real constant. The Hamiltonian $\mathcal{H}$ is given by Eq. 10. The self-energy terms $\Sigma^{r}_{1,2}$ model the coupling of the device to the semi-infinite leads. The self-energy is not hermitian, and its anti-hermitian part is responsible for the finite lifetime of the electron in the device. This subsequently contributes to broadening the energy levels in the device. The self-energies are computed using the surface-Green’s function, which requires an iterative procedure as outlined in Appendix A.

We compute the density of states ($d$) in the nanowire as the trace of the spectral Green’s function

[TABLE]

The real-valued singularities of the density of states are the Andreev bound state (ABS) energies. This is computed as a function of the phase difference ($\phi$) of the order parameter of the contacts and is shown in Fig. 3. The parameters for this computation are consistent with the Andreev approximationAndreev (1964, 1967); Beenakker (1992b); Ashida et al. (1989) ($\mu\gg\Delta_{0}$). As discussed in Appendix B, the breakdown of this approximation is manifested as an avoided level crossing in the ABS spectrum.

II.2 Current phase relationship

The current-phase relationship (CPR) links the macroscopic current flow in the junction to the phase gradient of the superconducting order parameterJosephson (1962); Spanton et al. (2017). The traditional approach to computing the CPR involves a demarcation of the bound state and continuum currents. The bound state current involves transport in the sub-gap energy range ($\absolutevalue{E}<\Delta_{0}$) while the continuum current is supported by the continuous energy spectrum outside the gap. Once the ABS spectrum is computed from scattering theory, a thermodynamic relation is used to calculate the bound state current, and the transmission formalism is used for the continuum current. The total current is the sum of the bound state and continuum currentsKulik (1969).

By contrast, when using the NEGF formalism the current-energy density can be computed at contact $i$, as a function of the phase difference $\phi$ using the current operatorDatta (1995, 2005)

[TABLE]

where $f(E)=1/\left(\exp\left(E/k_{B}T\right)+1\right)$ is the Fermi-Dirac occupation probability for a given energy level and $k_{B}$ is the Boltzmann constant. $G^{r(a)}$ and $\Sigma_{i}^{r(a)}$ are the retarded (advanced) Green’s function and contact $i$ self-energy respectively. To incorporate the opposite charge of electrons and holes we use the Pauli-z operator ($\tau_{z}$) in the particle-hole Nambu space. This current operator is reviewed in the Appendix A. The total current at a phase difference $\phi$ is then given by

[TABLE]

Figures 4 and 4 compares the current phase-relations for a short ($L<\xi_{0}$) and long junction ($L>\xi_{0}$) respectively, as calculated from ideal scattering theory and NEGF. Here, $\xi_{0}$ is the healing length as defined in the following section (Eq. 22). By the term “ideal scattering theory” we refer to a scattering approach which explicitly neglects normal reflections at the N/S interface in a clean junctionKulik (1969); Bagwell (1992); Gharavi and Baugh (2015); Beenakker (1992b). With this assumption of no normal reflections, the bound states in a clean junction cross at $\phi=\pi$. Hence, there’s a discontinuity at $\phi=\pi$ in the CPR calculated using this method. The NEGF result is expected to match scattering theory exactly in the $\mu\gg\Delta_{0}$ limit.

III Transverse subbands in Josephson Junctions

We now consider a more realistic three-dimensional model of the junction with a magnetic field along the nanowire axis, parallel to the direction of current flow. Figure 1 illustrates a discrete lattice model of the three-dimensional nanowire. The junction is along the $z$-direction and the transverse subbands are on the $x$-$y$ plane. In III-V semiconductors (InAs, InN) the charge carriers are typically confined close to the surface due to a positive surface potential, forming a surface accumulation layer. In accordance with this, we use a shell conduction model by including a large surface confining potential $U$ at the core of the nanowire (yellow sites in Fig. 1). The radial confinement and azimuthal periodicity of the nanowire gives rise to transverse subbands.

The single-electron Hamiltonian of the nanowire is

[TABLE]

where $z$ is the longitudinal direction, $U$ is the surface confinement potential, and $\mathcal{H}_{T}$ is the Hamiltonian of the transverse modes.

For a cylindrical nanowire, the rotational symmetry about the longitudinal axis results in angular momentum $\left(\ell\right)$ subbands. This is because

[TABLE]

where $\hat{L}_{z}$ is the angular momentum operator in the z-direction. Hence, $\ell$ is a good quantum number. The $\ell$ subbands are eigenstates of the $\hat{L}_{z}$ operator, labelled by their eigenvalue

[TABLE]

We can consider the square cross-section in Fig. 1 as a perturbation to an ideal cylindrical geometry. Each of the transverse subbands in a square cross-section can be written as a superposition of angular momentum eigenstates. In Sec. III.2, we compute the average angular momentum of each transverse subband, and observe only a small difference (see Fig. 5) from the unperturbed quantized eigenvalues ($\hbar\ell$). We will thus work within the zeroth order of this perturbation and use the language of angular momentum subbands in our analysis.

The azimuthal motion of the Andreev quasiparticles couples with the applied magnetic field, resulting in a quasiparticle phase pickup. Oscillations in the maximal supercurrent (critical current) with field have been measured by Gharavi et al., and Zuo et al.. Unlike the Fraunhofer interference in wide planar junctions, the field is aligned with the current and the oscillations do not show any periodicity.

Using Peierls substitution, we include the orbital effect of the vector potential in the phase of the transverse hopping. For a constant magnetic field along the $z$-direction, the vector potential can be written as

[TABLE]

Within the tight-binding approximation, the on-site and hopping matrices in the particle-hole Nambu space are given by

[TABLE]

where ${h}=2t_{x}+2t_{z}+\absolutevalue{t_{y}}\left[2+(2\pi n_{x}\Phi_{a})^{2}\right]-\mu$

[TABLE]

[TABLE]

where $\Phi_{a}$ is the flux quanta per unit cell of the nanowire cross-section, and $n_{x}$ is the lattice site index in the $x-$direction. This factor alters the on-site energy ($\alpha_{N/S}$) and contributes a phase to the hopping term corresponding to the gauge chosen for the vector potential (Eq. 18).

III.1 Andreev bound states in a magnetic field

Figures 5(a), and 5(b) plot the subgap density of states as obtained from the spectral Green’s function (Eq. 12) for a nanowire with an InAs effective mass $m^{*}$ = 0.023$m_{e}$Nakwaski (1995) ($m_{e}$ is the bare electron mass), radius R = 30 nm and chemical potential $\mu=5\Delta_{0}$. As described in Appendix B, the bound states anti-cross at $\phi=\pi$ due to normal reflections at the N/S interfaces. A normalised flux of $\Phi=0.01$ ($0.03$) is applied in Fig. 5(a) (Fig. 5(b)), which breaks symmetry in the transverse direction (Eq. 18) and lifts the degeneracy of the $\pm\ell$ subbands. Here, $\Phi=B\cdot S/\Phi_{0}$, where $S$ is the cross-sectional area, and $\Phi_{0}=h/e$ is the flux quantum.

The chemical potential is adjusted to populate three subbands. In Fig. 5 we see that the $\ell=0$ subband remains unaffected, while the $\ell\neq 0$ subbands split in the presence of an axial flux. The splitting is proportional to the flux and the subband angular momentum. The process used to label the subbands is described in Sec. III.2.

The characteristic length-scale associated with an occupied subband is called the healing length $(\xi_{\ell})$Kulik (1969); Bagwell (1992), and is given by

[TABLE]

$v_{F,\ell}$ is the fermi-velocity and is given by

[TABLE]

We have a “short junction” when the nanowire is shorter than the healing length ($L<\xi_{\ell}$). Note that the healing length depends on the angular momentum quantum number, and the classification of the junction as long/short is subband dependent.

We now describe a procedure to label the angular-momentum subbands using the correlation Green’s function ($G^{n}$).

III.2 Average angular momentum of the transverse subbands

The angular momentum of the subbands can be computed as the expectation of the $\hat{L}_{z}$ operator (Eq. 17), however, we do not have access to the wavefunctions in a numerical simulation. We do have the correlation Green’s function $-iG^{<}=G^{n}$ which gives the particle-hole density per-unit energy. Using this, we find the expectation of the $\hat{L}_{z}$ operator as a function of energy

[TABLE]

where $\max_{E}\{\Tr\left[G^{n}(E)\right]\}$ returns the peak subgap quasiparticle concentration. In the Cartesian coordinate system, the $\hat{L}_{z}$ operator is written as

[TABLE]

The position operators $\hat{x},\hat{y}$ are diagonal in the tight-binding basis, with each entry a multiple of the lattice constant ($a$). For example, if we consider two points along the $x$ and $y$ axis, and one along $z$, we have the following position operators

[TABLE]

The linear momentum operators can be written as

[TABLE]

where the subscript $i$ is used to denote $x,y$ basis of the position and momentum operators. Using Eq. 25, 26, 27, we write the $\hat{L}_{z}$ operator and using Eq. 24 we find the expectation of $\hat{L}_{z}$ as a function of energy.

We employ the above procedure to compute the average angular momentum of the Andreev bound states in an SNS junction. The chemical potential is fixed to give us three occupied subbands ($\ell=0,\pm 1$). The average angular momentum of the subbands at $\Phi=0.01$, $\phi=\pi$ is shown in Fig. 5(c). Next, we increase the axial magnetic flux through the nanowire to $\Phi=0.03$. The $\ell=\pm 1$ subband states further split (Fig. 5(b)), and the angular momentum of the subbands is plotted in Fig. 5(d) for $\phi=\pi$. The zero angular momentum subband ($\ell=0$) has no azimuthal motion, and hence is unaffected by the axial field. A non-zero $\langle L_{z}\rangle$ for the $\ell=0$ subband (Figs. 5(c), 5(d)) results from its hybridisation with the nearby $\ell=-1$ subband. This consequently decreases the $\langle L_{z}\rangle$ magnitude for the $\ell=-1$ subband w.r.t. $\ell=+1$.

From Figs. 5(c),5(d) we infer that each transverse subband is primarily composed of a single angular momentum eigenstate. This resemblance to the subband structure of an ideal cylindrical nanowire is the basis for Sec. III.4, where we build an effective subband model by including of a fixed number of angular momentum subbands. In the next section we explain the importance of including only the zero angular momentum subband in the superconducting contacts.

III.3 Zero angular momentum subband in the Superconductor

The BdG Hamiltonian in a superconductor is given by

[TABLE]

For a cylindrical geometry with an azimuthal vector potential, $\mathcal{H}_{0}$ is given by

[TABLE]

As discussed, the radial confinement due to the nanowire’s cylindrical geometry gives rise to angular momentum subbands labelled by $\ell$.

We will analyze the eigenenergies of this superconductor in the presence and absence of a magnetic field.

III.3.1 Zero magnetic field, $\mathbf{B=0,A=0}$

Using the ansatz wavefunction $\exp(ik_{z}z)\exp(i\ell\theta)$, the diagonal elements of the BdG Hamiltonian can be written as

[TABLE]

where $\mu_{\ell}$ is the effective chemical potential

[TABLE]

The BdG Hamiltonian simplifies to

[TABLE]

and its eigenvalues $E$ are given by

[TABLE]

This is the well-known superconductor dispersion relation, with a gap of $\Delta_{0}$ on either side of the fermi level.

III.3.2 Constant axial magnetic field, $\mathbf{B}=B_{z}\hat{z}$

In the Coulomb gauge we can write the vector potential for this magnetic field as

[TABLE]

From Stoke’s law

[TABLE]

Exploiting the symmetry of the cylindrical geometry, the above equation can be simplified to

[TABLE]

Using the same ansatz $\exp(ik_{z}z)\exp(i\ell\theta)$, the diagonal elements can be written as

[TABLE]

[TABLE]

for the electron and hole parts respectively. The effective chemical potential $\mu_{\ell}$ is defined in Eq. 30, and the field-coupling term $\mathcal{E}_{\ell}=\frac{\hbar^{2}\left(2\ell\Phi\right)}{2m^{*}R^{2}}$. We note that

[TABLE]

The BdG Hamiltonian can then be written as

[TABLE]

and the eigenvalues $E$ are given by

[TABLE]

Thus, we see that a magnetic field “shifts” the superconducting gap. It is no longer centred at the fermi level.

While this may be a good model for a superconducting “nanowire”, experimental setups usually involve a superconductor sputtered into quasi-planar contactsGharavi et al. (2014); Zuo et al. (2017) which naturally support only the $\ell=0$ subband. The geometry of the superconducting contacts implies a large $\ell\neq 0$ subband energy, and can safely be assumed to remain unpopulated. This difference in geometry prompts the inclusion of $\ell\neq 0$ subbands in the nanowire, and their exclusion in the superconductor.

For $\ell=0$, we have $\mathcal{E}_{\ell}=0$, and thus the superconducting gap stays centered at the fermi level. This is illustrated in Fig. 6 – when $\ell\neq 0$ subbands are included in the superconductor (Fig. 6(a)), the ABS move vertically as a result of the shift in the superconducting gap. However, with only the zero angular momentum subband $\ell=0$ in the superconductor (Fig. 6(b)), the ABS curves are horizontally phase shifted whilst the superconducting gap remains fixed. This shift is proportional to the applied flux, and the angular momentum of quasiparticles occupying the subband.

From Eqs. 36, 37, the effective chemical potential for electron-like (hole-like) particles in the $\ell$ subband in the N-section is given by

[TABLE]

The electron and hole wavenumbers can then be written as a function of energy $(E)$

[TABLE]

III.4 The 1-dimensional effective subband model

As outlined above, it is important to ensure that we only have the $\ell=0$ subbands in the contacts. We also note from Eqs. 30, 41 that we can incorporate the effect of the angular momentum subbands via an effective potential $\mu_{\ell}$, and a field-coupling term $\mathcal{E}_{\ell}$.

The tight-binding Hamiltonian of the nanowire can be written as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $N_{\ell}$ is the number of subbands. (For example, $N_{\ell}=3\implies$ ($\ell={-1,0,1}$)) and $\mathbb{I}_{n}$ is the $n\times n$ identity matrix.

Meanwhile, the Hamiltonian of the contacts takes a similar form with

[TABLE]

IV Supercurrent Oscillations

IV.1 Clean junction

We compute the CPR of an SNS junction at finite axial magnetic fields, assuming a shell conduction model, and a nanowire diameter of 60 nm. Temperature is set to $T=100$ mK in all the simulations. In Fig. 7(b), we show the CPR as a function of the magnetic flux for a single occupied subband. With only the $\ell=0$ subband populated, there is no phase shift in the ABS, and the CPR retains its saw-tooth shape with a maximum near $\phi=\pi$. The critical current as a function of the flux is plotted in Fig. 7(a). The gradual fall in the critical current can be attributed to the decrease in average quasiparticle momentum with increasing flux, as shown in Eq. 23. Eventually, at $\Phi=4.04$ the band depopulates $\left[\min\left(\mu_{\ell}^{e},\mu_{\ell}^{h}\right)=0\right]$ and the current falls to zero. We observe in Fig. 7(a) that the critical current does not monotonically decrease to zero, particularly for $\Phi\in[3,4]$. The appearance of these small oscillations is due to the interference with the quasiparticles normally reflected from the N/S interfaces. As discussed in Appendix B, the discontinuity in the density of states gives rise to normal reflections. These reflected quasiparticles interfere and result in the non-monotonic decrease of the single subband critical current.

Next, we consider the case when three subbands are occupied ($|\ell|\leq 1$). The magnetic field evolution of the CPR is plotted in Fig. 8. Since the ABS for the $|\ell|=1$ subbands are phase shifted in presence of a flux (Fig. 6(b)), the total current is no longer maximum near $\phi=\pi$. The current in the junction is the sum over the individual subband currents, and consequently, the flux-dependent phase shift results in an interference pattern for the field evolution of the critical current. The phase shift in a subband CPR is proportional to the difference in the electron-hole wavenumbers $(k_{\ell}^{e}-k_{\ell}^{h})$ (Eq. 42), and the length $L$ of the junction. Hence, the fluxes at which the subband currents constructively interfere need not occur at integer multiples of the flux quantum $\Phi_{0}=h/e$.

In Fig. 8(a) we plot the critical current for three occupied subbands as a function of the axial flux. We see several oscillations of the critical current before the $|\ell|=1$ subbands depopulate at $\Phi=3.04$. At zero flux, the CPR of each subband is maximum near $\phi=\pi$ and hence they all add up constructively. As illustrated in Fig. 8(b), each subband contributes equally to the critical current.

As the flux is increased, the electron-hole pairs in the $|\ell|=1$ subbands pickup a phase and the subband CPRs no longer interfere constructively. Consequently, the critical current decreases with flux. At $\Phi=0.72$, the $\ell=0$ and $\ell=\pm 1$ subband CPRs are maximally out-of-phase, resulting in a local minima (node). The subband and total CPR at this node is plotted in Fig. 8(c). Beyond $\Phi=0.72$, the critical current switches phase from $\phi<\pi$ to $\phi>\pi$ and the current increases again. This increase persists till $\Phi=1.08$ at which point the current is maximum near $\phi=2\pi$. At this flux, the $|\ell|=1$ subband current peaks near $\phi=2\pi$ while it is negligible near $\phi=\pi$. Hence, this secondary peak – which only involves contribution from $|\ell|=1$ subbands – is approximately a two-third of the primary peak and corresponds to a phase pickup of $\pi$ in the aforementioned subbands. The subband and total CPR for the secondary peak is shown in Fig. 8(d). As noted earlier, the magnitudes of the primary and secondary peaks progressively diminish due to the decrease in average quasiparticle velocity.

Finally, we consider the situation when five subbands are occupied ($|\ell|\leq 2$). The critical current is plotted as a function of the magnetic flux in Fig. 9. Once again, at $\Phi=0$ the subband currents are all in-phase and constructively interfere to give a maximum. In presence of a magnetic field, the quasiparticles in the $|\ell|=1$ and $|\ell|=2$ subbands pick up different phases and hence, they do not appear to constructively interfere again in presence of a magnetic field to recover the zero field critical current.

The absence of such oscillations with a single occupied subband (Fig. 7(a)) confirms the subband supercurrent interference as the causal agent.

IV.2 Effect of disorder

In order to simulate experimentally relevant conditions, we include a random uncorrelated onsite disorder potential $u\in[-W,W]$ in the semiconductor. This models phase-coherent scattering events in the junction. We parameterise the disorder by the mean free path ($\lambda_{mf}$), which is estimated from the disorder-averaged normal state conductance ($g$) using the following relation

[TABLE]

$N_{\ell}$ is the number of subbands and $L$ is the length of the junction.

In Fig. 10 we plot the critical current oscillations in a nanowire for particular realisations of the disorder. While the initial decay and the oscillations are still present, the secondary maxima are suppressed. In a clean nanowire with a saw-tooth CPR (which peaks near $\phi=\pi$ at zero field), at a magnetic flux $\Phi^{*}$ the $|\ell|=1$ subbands pickup a phase of $\pi$ and their CPRs peak near $\phi=2\pi$. The $\ell=0$ subband retains its sawtooth CPR with a negligible current near $\phi=2\pi$. As described in the previous subsection, this results in the secondary maximum. Upon adding disorder to the nanowire, we depart from this saw-tooth CPR, tending towards a sinusoidal CPR which peaks further away from $\phi=\pi$ at zero field. Thus, there exists no $\Phi^{*}$ at which the $|\ell|=1$ subband current peaks while the $\ell=0$ subband current is negligible. As a consequence of the sinusoidal CPR, on picking up a phase of $\pi$ the $|\ell|=1$ subbands destructively interfere with the $\ell=0$ subband, and this causes the suppression of the secondary maxima.

In the presence of scatterers the effective path traversed by the quasiparticles increases and hence, the subbands destructively interfere at a lower flux. This is shown in Fig. 11, where the first crticial current node in a disordered junction occurs at a lower field as compared to the clean nanowire.

Furthermore, as shown by Zuo et al., the essential effect of disorder can be observed by the dependence of the critical current oscillations on the gate voltage. As shown in Fig. 12(a) for the clean nanowire, small variations in the gate voltage hardly cause any fluctuations in the oscillations. This is because small changes in the chemical potential do not change the number of occupied subbands and only weakly affects the quasiparticle transmission through the junction. However, in a disordered nanowire with a small mean free path, the quasiparticles traverse a longer path in the nanowire and hence, the critical current oscillations are significantly affected by the gate voltage. This is shown in Figs. 12(b),12(c) for two disorder realisations.

From Figs. 10,12 we infer that the critical current oscillations are highly sensitive to the gate voltage and the particular realisation of the disorder. Thus, a macroscopic current measurement indirectly gives us information about the microscopic specifics of the junction. However, while our model provides a qualitative understanding of the oscillations, the high sensitivity w.r.t. the microscopic parameters renders a quantitative description of the experiment highly challenging.

IV.3 Dephasing in the nanowire

The analysis presented in the previous sections described the phase-coherent flow of quasiparticles in the junction. We now include phase-breaking processes Golizadeh-Mojarad and Datta (2007); Datta (2005) that may arise from any time-dependent potential, such as the interaction of an electron with the surrounding bath of phonons, other electrons, or fast-fluctuating charge noise caused by traps in defects. Although it is non-trivial to identify the dominant source of dephasing, we can adopt a phenomenological model by introducing a Büttiker probe for the lattice backgroundDatta (2005).

We subsume these processes within the NEGF formalism by including a self-energy term for the lattice background $\Sigma^{r}_{s}$, proportional to the Green’s function and the emission-absorption coefficients. This calls for a self-consistent computation of the Green’s function and the bath self-energy Golizadeh-Mojarad and Datta (2007)

[TABLE]

where $\times$ denotes element by element multiplication. The elements of the matrix $D$ represent the correlation of the time-dependent interaction potential between pairs of lattice points. Adopting a homogeneous model, we assume uniform, elastic, and spatially uncorrelated interactions resulting in a diagonal $D$,

[TABLE]

for every pair of coordinates $z_{i}$ and $z_{j}$ in the nanowire. This model discards the off-diagonal elements of the Green’s function, and hence relaxes both the phase and momentum of the quasiparticles in the nanowire.

The magnitude of $D_{0}$ limits the phase relaxation length of the junction. Molecular beam epitaxy(MBE) and metalorganic vapour phase epitaxy grown InAs nanowires typically have a phase relaxation length on the order of a few-hundred nanometers $l_{\varphi}\sim 100-300$ nmBlömers et al. (2011b); Estévez Hernández et al. (2010b); Doh et al. (2007). In our model, $l_{\varphi}$ can be estimated from the statistical properties of universal conductance fluctuations (UCF)Blömers et al. (2011b); Estévez Hernández et al. (2010b); Doh et al. (2007) The details of this calculation is presented in Appendix D, and the results are tabulated here in Table 1

The critical current oscillations in the presence of elastic dephasing interactions, as listed in Table 1, are shown in Fig. 13. The green curve in Fig. 13 includes a random disorder potential profile in addition to the phase-breaking processes. With dephasing in the nanowire, the excess path traversed due to the disorder potentials does not result in a proportionate phase pick-up. Consequently, in contrast to Fig. 11, the inclusion of disorder does not result in a significant shift in the critical current nodes.

One of the main effects of dephasing is in pinning down the critical current nodes so that they’re less sensitive to disorder. A phase-coherent simulation would overestimate the disorder-induced quasiparticle phase pickup and consequently, the computed critical current oscillations are sensitive and strongly dependent on the microscopic disorder profile (see Fig. 10). However, measured oscillations from Zuo et al. exhibit only a gradual variation in the oscillations with changing disorder realizations – well modeled by the inclusion of dephasing in the nanowire. Further, dephasing is necessary to reproduce the observed reduction of the relative peak height for subsequent primary maxima. Phase-breaking processes restrict the coherent life time of the quasiparticles in the nanowire, which hinders a constructive interference at higher fields. This is well illustrated in Fig. 13, where the relative peak height of the second primary maximum (third peak) decreases with an increase in dephasing strength.

Magnetoconductance calculations by Lahiri et al. reveal the suppression of higher harmonics of the transmission characteristics in a nanowire with dephasing. This is manifest in experiments as an increasing inter-node spacingZuo et al. (2017); Gharavi et al. (2014) with respect to the axial field. As illustrated in Fig. 13, this feature is captured by the inclusion of phase-breaking processes in our model.

The oscillations in Fig. 13 are in good qualitative agreement with the experiments by Gharavi et al. and Zuo et al.. Thus, we infer that phase-breaking processes play a non-negligible role in III-V semiconductor nanowire Josephson junctions.

V Conclusion

In this paper we employed the Keldysh Non-Equilibrium Green’s Function formalism to model quantum transport in semiconductor nanowire Josephson junctions. In our analysis we used a three-dimensional discrete lattice model described by the Bogoliubov-de Gennes Hamiltonian in the tight-binding approximation, and computed the Andreev bound state spectrum and current-phase relations. We went beyond the Andreev approximation limit and investigated the avoided level crossing in the ABS spectrum. Our results confirm the measured critical current oscillations to arise from the subband supercurrent interference in presence of an axial magnetic field. The phase picked up by the quasiparticles depend on the difference of their wavenumbers, the length of the junction and the angular momentum quantum number. Thus, the oscillations do not show any periodicity in the flux quantum through the nanowire cross-section. We included phase-coherent scattering to model a disordered junction and investigated its effect on the critical current oscillations. We observed that the oscillations in the disordered junction are highly sensitive to the realisation of the random disorder potential, and on small fluctuations of the gate voltage. This high sensitivity makes a quantitative description of the experiment a challenging task. Nevertheless, a macroscopic current measurement conveys valuable information about the microscopic profile of the junction. We include elastic dephasing in the nanowire by modelling weak phase-breaking interactions. A good qualitative match of our results with the experiment is observed, and this underscores the role played by phase-breaking processes in III-V nanowire Josephson junctions.

The relevance of these results is emphasized by outlining the points of comparison with experimental data from Gharavi et al. and Zuo et al.. The data exhibits a strong suppression of the switching current at magnetic fields on the order of $100-500\text{ mT}$. Subsequently, oscillations with aperiodic nodes are observed in the field dependence of the switching current. As shown in Fig. 13, our simulations capture the characteristic features of this evolution. Furthermore, as observed in the experiments, the simulated oscillations display a strong gate-tunability and are not uniquely determined by the junction geometry. Finally, a phase-coherence length in the range $100-300$ nm of MBEBlömers et al. (2011b) and metalorganic vapor phase epitaxyEstévez Hernández et al. (2010b); Doh et al. (2007) grown InAs nanowire samples corroborates the dephasing in our model, as introduced by a uniform, spatially uncorrelated dephasing parameter $(D_{0})$.

*Acknowledgements: * The authors BM and PS would like to thank Prof. Supriyo Datta, Prof. Kantimay Das Gupta, Abhishek Sharma and Aritra Lahiri for useful discussions throughout this work. This work is an outcome of the Research and Development work undertaken in the project under the Visvesvaraya PhD Scheme of Ministry of Electronics and Information Technology, Government of India, being implemented by Digital India Corporation (formerly Media Lab Asia). This work was also supported by the Science and Engineering Research Board (SERB) of the Government of India under Grant number EMR/2017/002853. KG and JB acknowledge funding from the Natural Sciences and Engineering Research Council of Canada (NSERC).

Appendix A Non Equilibrium Green’s Function Formalism (NEGF)

The retarded Green’s function in the energy domain is given by

[TABLE]

where $\mathcal{H}$ is the Hamiltonian, $\Sigma_{1,2}^{r}$ are the retarded self-energies of the semi-infinite contacts, and $\eta$ is an infinitesimal real constant. The advanced Green’s function is the Hermitian conjugate of the retarded Green’s function ($G^{a}=G^{r\dagger}$). The Hamiltonian $\mathcal{H}$ is written in the tight-binding approximation Eq. 10. The surface Green’s functions ($g_{s}$) at each contact are recursively evaluated

[TABLE]

where the subscript $L,R$ labels the contact. Using this we compute the self-energy

[TABLE]

where $\sigma_{1}=\beta^{\dagger}g_{sL}\beta$, and $\sigma_{2}=\beta g_{sR}\beta^{\dagger}$

The anti-Hermitian part of the self-energy is responsible for the finite life-time of the quasiparticles in the junction and broadens the energy levels. This broadening matrix is denoted by $\Gamma_{i}$.

The Fermi functions in the particle-hole Nambu space is given by

[TABLE]

where $f(E,\mu)=1/\left(\exp\left[(E-\mu)/kT\right]+1\right)$ is the fermi function, and $V$ is the bias applied to the contact.

The lesser self-energy, or the inscattering matrix can be computed from the broadening matrix and fermi function as

[TABLE]

The lesser Green’s Function is then computed

[TABLE]

Next, we construct the current operator

[TABLE]

Electrons and holes travelling in the same direction carry opposite currents and hence, the current is given by the difference of the partial trace of the current operator over the electron and hole sub-spaces.

[TABLE]

This can be also be written as

[TABLE]

where $\tau_{z}$ is the Pauli operator in the particle-hole Nambu space. The total current is then evaluated by integrating the current-energy density

[TABLE]

There’s a small technical caveat to keep in mind when using the NEGF current operator.

Using the equations for the retarded and advanced Green’s functions,

[TABLE]

where $\Gamma=i(\Sigma^{r}-\Sigma^{a})$

Pre-multiplying by $G^{r}$ and post-multiplying by $G^{a}$,

[TABLE]

Multiplying the fermi-function,

[TABLE]

Thus,

[TABLE]

However, if we don’t consider the term proportional to the infinitesimal $\eta$ we end up with

[TABLE]

$G^{<}$ from Eq. 69 misses a term in the current proportional to ${\left(G^{r}G^{a}\right)}$, the trace of which increases with the number of bound-states. This leads to erroneous results for longer nanowires, which have a larger number of Andreev Bound States. The NEGF current-operator for contact $i$ is given by –

[TABLE]

Substituting $G^{n}=Af$ and $\Sigma^{in}=\Gamma f$, the current operator can be simplified to

[TABLE]

Appendix B Andreev bound state spectrum : Beyond Andreev approximation

Andreev reflections across an N/S interface were first analyzed by Blonder et al. in 1982, and have been prevalent in the literature ever since. Almost always, these results are derived under the Andreev approximationAndreev (1964, 1967); Beenakker (1992b); Ashida et al. (1989). This approximation deals with a regime where the chemical potential of the nanowire is much larger than the superconducting order parameter of the leads ($\mu\gg\Delta_{0}$). In this appendix, we analyze the implications of working in a regime where the Andreev approximation is not valid. Specifically, we consider the implications of being outside the Andreev approximation regime on the ABS spectrum in a clean 1-dimensional SNS junction.

The density of states in a superconductor is gapped by an energy $\Delta_{0}$ on either side of the fermi level. There is no gap in the normal state spectrum. Thus, quasiparticles in the subgap region face an energy barrier at the interface. This is an energy barrier between states with the same momentum, arising due to a difference in the order parameter ($0,\Delta_{0}$) across the interface. This has nothing do with an impurity or any non-ideality of the junction. However, it plays a role very similar to any impurity-induced barrier $U$ at the interface, i.e. it gives rise to normal reflections at the interfaces, which cannot be neglected when $\mu\gg\Delta_{0}$ is not valid, even for a clean junction.

We now examine the spectrum for the chemical potential $\mu$ comparable to $\Delta_{0}$. Figure. 14 plots the ABS spectrum for $\mu=0.5\Delta_{0}$.

We observe an avoided level crossing when $\mu\not\gg\Delta_{0}$, as shown in Fig. 14. This anti-crossing can be attributed to the normal reflections which become significant in this parameter space. When these normal reflections are negligible, we have independent rightward and leftward moving excitations in the nanowire, resulting in a crossing at $\phi=\pi$. However, once normal reflections become important, as is the case outside the Andreev approximation, the excitations moving in opposite directions get coupled to each other. This interaction between them, brought about by the normal reflections at the N/S interfaces, results in the anti-crossing.

In Fig. 15, we plot the ABS energy gap ($\delta_{\pi}$) at $\phi=\pi$ as a function of $\mu/\Delta_{0}$. It is evident from this plot that $\delta_{\pi}$ decreases as $\mu/\Delta_{0}$ increases. Thus, in the Andreev approximation regime ($\mu/\Delta_{0}\rightarrow\infty$) the bound states cross ($\delta_{\pi}\rightarrow 0$) at $\phi=\pi$. When $\Delta_{0}$ is kept constant and $\mu$ is varied, the gap varies as $(\mu/\Delta_{0})^{-1}$ for $\mu>\Delta_{0}$. In Appendix C we verify this dependence analytically by taking a scattering theory approach. Figure 15 plots the variation of the gap at $\phi=\pi$ with the nanowire length. These oscillations result from the interference of the waves reflected at the two S/N interfaces.

Appendix C The Andreev bound state quantisation condition outside the Andreev Approximation Regime

In this appendix we derive an expression for the Andreev bound state energy in a SNS junction, when the phase difference between the two superconductors is tuned to $\pi$. We show this to be non-zero for a finite $\mu/\Delta$ while approaching zero in the limit $\mu/\Delta\rightarrow\infty$, where $\mu$ is the chemical potential of the entire device and $\Delta$ is the induced superconducting gap.

C.1 Setting up the problem

We follow the notation from Kulik (1969). Consider first the transport of an electron across an N/S interface. It is well known that for a clean interface, an incident electron is perfectly Andreev reflected as a hole. However, there is also a component of the order $\Delta/\mu$ which is reflected as a normal electron. This is mentioned in Kulik’s work and can also be derived from first principles scattering theory at the N/S interface. We usually encounter cases when perfect Andreev reflection is assumed with no normal reflectionAndreev (1964, 1967); Beenakker (1992b); Ashida et al. (1989), but this is only true to the zeroth order in $\Delta/\mu$.

This normal reflection component will modify the Andreev bound state quantization conditionGharavi and Baugh (2015); Kulik (1969)

[TABLE]

where $\frac{1}{\gamma}=\frac{E}{\Delta}+i\sqrt{1-\frac{E^{2}}{\Delta^{2}}}$, $k_{0}$ and $k_{1}$ are the electron and hole wave vectors and $\chi=\chi_{2}-\chi_{1}$ is the phase difference between the superconducting order parameter across the junction.

Consider an SNS junction with a normal region length $d$ extended over $|z|<d/2$, and superconducting contacts defined over $|z|>d/2$. Writing the two-component wave functions

[TABLE]

The terms proportional to $B^{\prime}$, $C^{\prime}$ and $D^{\prime}$ are ignored when perfect Andreev reflection is assumed. We have included a $\Delta/\mu$ coefficient to emphasize a first order expansion beyond the perfect Andreev reflection scenario.

Equating the coefficients at $z=\pm d/2$, we get:

[TABLE]

Solving for $C$ and $D$ (to be compared with Eq 2.13 from Kulik), we get:

[TABLE]

[TABLE]

We now divide the above equations and keep terms to first order in $\Delta/\mu$,

[TABLE]

We will now focus on the qualitative behaviour of the solutions of the above equation. In order to do so, we simplify the above equation into a more tractable form:

[TABLE]

where the exact form of $\epsilon$ and $\eta$ can be derived from equation 80.

C.2 Expression for the energy at $\chi=\pi$

Let us tune the phase difference $\chi=\pi$. Defining $\phi=\cos^{-1}({E/\Delta})$, we can write $\gamma=e^{-i\phi}$. Using this relation, the quantization condition can be simplified to,

[TABLE]

Equating the real and imaginary parts,

[TABLE]

Since $\Delta/\mu\ll 1$, under the small angle approximation Eq. 84 can be simplified to

[TABLE]

For $E/\Delta\ll 1$

[TABLE]

This can be further simplified using \left(\text{\frac{\hbar^{2}k_{0/1}^{2}}{2m}=\mu\pm E}\right),

[TABLE]

[TABLE]

This is a finite energy for non-zero $\eta$ and goes to zero in the limit $\Delta/\mu\rightarrow 0$.

C.3 Comparison with numerics

Ignoring the term in the bracket from the final expression for E, we can simplify it to $E\approx\frac{\Delta^{2}}{\mu}$. We decided to observe the power law dependence of the gap using the numerical simulations. We considered two cases:

Fix $\mu$ and vary $\Delta$: we expect to see a behavior

[TABLE] 2. 2.

Fix $\Delta$ and vary $\mu$: we expect to see a behavior

[TABLE]

Figure 16 confirms this dependence of the ABS energy on $\mu/\Delta$ for the aforementioned cases.

Appendix D Phase relaxation length estimation

The dephasing in the nanowire can be parameterised by the phase relaxation length $l_{\varphi}$, which is a length scale over which the phase of the quasiparticles randomise. As explained in Sec. IV.3, phase-breaking processes are included via a self-energy for the lattice background $\Sigma_{s}=D\times G^{r}$; $D=D_{0}\mathbb{I}$. In this appendix, we estimate $l_{\varphi}$ as a function of the dephasing strength $D_{0}$.

D.1 Estimation from phase coherence lifetime

The phase relaxation length can be computed from the phase coherence lifetime $\tau_{\varphi}$

[TABLE]

where $v_{F}$ is the Fermi-velocity, and $\mathcal{D}$ is the diffusion constant. The anti-hermitian part of the lattice background self-energy limits the phase-coherent lifetime of the quasiparticles, and sets an energy scale for the problem. The phase coherence lifetime can thus be estimated as

[TABLE]

where $\Gamma_{s,d}$ is a diagonal element of $\Gamma_{s}=i\left(\Sigma_{s}-\Sigma_{s}^{\dagger}\right)$ is the broadening function corresponding to the lattice background. Note that $\Gamma_{s}$, and hence the $l_{\varphi}$ so estimated from Eq. 91 is a function of energy. The phase relaxation length is then reported as an average over the energy grid, $l_{\varphi}=\langle l_{\varphi}(E)\rangle_{E}$. Under the assumption of ballistic transport, we estimate an upper bound on the phase relaxation length $l_{\varphi}\underset{\sim}{<}150$ nm for $D_{0}=0.001$ eV2, and for $D_{0}=5\times 10^{-4}$ eV2, $l_{\varphi}\underset{\sim}{<}300$ nm. The Fermi-velocity decreases with field, and hence we observe a gradual monotonic degradation in the phase relaxation length.

D.2 Estimation from statistical properties of UCF

For a nanowire length comparable to $l_{\varphi}$, the normal-state conductance fluctuates with an amplitude of the order of $e^{2}/h$ in presence of a magnetic field. These aperiodic universal conductance fluctuations (UCF) measured in a magnetic field perpendicular to the nanowire axis can be analysed to extract information on phase coherent transport. The UCF originates from electron phase shifts resulting from the penetration of magnetic flux through closed electron trajectories. The conductance shows strong fluctuations for low dephasing strengths, while they are smeared out at higher coupling strengths.

The magnetoconductance fluctuation is denoted by $\delta G$

[TABLE]

where the average $\langle.\rangle$ is taken over the magnetic field $B$. The average fluctuation amplitude about the mean conductance is quantified by the root-mean-square $\text{rms}(\delta G)_{B}=\sqrt{\text{var}{(\delta G)}_{B}}$. The rms$(\delta G)_{B}$ decreases monotonically with $D_{0}$. The phase relaxation length $l_{\varphi}$ can be estimated from the analysis of the autocorrelation function $F$ of $\delta G$. The half-width half-maximum (HWHM) of $F$ corresponds to the correlation field $B_{c}$, which is a measure of a field range over which the phases of the interference path become uncorrelated.

[TABLE]

Assuming the phase relaxation length $(l_{\varphi})$ to be greater than the nanowire diameter $d$, we can extract $l_{\varphi}$ directly from the correlation fieldBlömers et al. (2011b); Estévez Hernández et al. (2010b); Beenakker and van Houten (1988)

[TABLE]

where $\gamma$ is a dimensionless prefactor depending on the transport regime. We work in the the dirty metal limit with $\gamma=0.95$Blömers et al. (2011b).

D.2.1 Results

This simulation involves normal-state low-bias transport in presence of a magnetic field oriented in a direction perpendicular to the nanowire axis. The nanowire length $L=200$ nm, and diameter $d=30$ nm. To model diffusive transport, an onsite random potential in the range $W\in[-1.5t,1.5t]$ is introduced at each point in the nanowire, where $t$ is the tight-binding hopping parameter. This corresponds to a mean-free path $\lambda_{mf}\approx 18$ nm. The magnetoconductance fluctuations are plotted in Fig. 18(a). The normalised autocorrelation $F/F(0)$ of $\delta G$ is shown in Fig. 18(b), and the extracted parameters are listed in Table 2.