{'fraction_non_alphanumeric': 0.036288242150923924, 'fraction_numerical': 0.012210735510500469, 'mean_word_length': 4.500525541307547, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Supercurrent Interference in Semiconductor Nanowire Josephson Junctions

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.

Signatures of MBS have been reported as a zero-bias conductance peak in tunnelling experiments 1-5, 16 . The 4π Majorana-Josephson effect has been predicted and observed in nanowire Josephson junctions tuned to the topologically non-trivial regime 6,9 . 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 processing [17][18][19][20][21][22] . The 4π Majorana-Josephson effect forms the basis of braiding and readout schemes of a recent topological qubit proposal 23 .
Josephson junctions based on semiconductorsuperconductor hybrids form the basis for microwave quantum circuity 24 , and superconducting qubits 25,26 . They afford an attractive alternative for a scalable computing architecture with the possibility of an all-electric qubit control 24,27,28 .
Several studies have focused on the structure of transverse subbands and magnetoconductance due to radial confinement in semiconductor nanowires [29][30][31][32] and carbon nanotubes 33 . 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 regime [34][35][36] . 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) [37][38][39][40][41][42] 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 setups [19][20][21][22] and in interpreting experiments, particularly for those based on critical current measurements [43][44][45][46] .
Quantum transport traditionally involves excited states and the use of a variant of the Landauer-Büttiker's scattering theory 38,47,48 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 level 49,50 . Blonder et al. generalised the scattering theory approach to hybrid semiconductor-superconductor junctions by solving the Bogoliubov-de Gennes equation across the N-S interface 51 . Beenakker applied this formalism for mesoscopic N-S junctions, thus providing a multichannel generalization of Blonder's results 52 . This technique has been prevalent in the literature 53,54 ever since, and it forms the basis for numerous simulation packages such as Kwant 55 . 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 function 56,57 , 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 (µ) to be much larger than the superconducting order parameter (∆ 0 ), i.e., µ ∆ 0 51,54,58,59 . We go beyond this Andreev approximation limit and investigate the bound states which anti-cross at a superconducting phase difference of π 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 B 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 phasebreaking 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 reflection 60,61 . Low bias transport in normal (N)-superconductor (S) junctions involves Andreev reflections at the interface. We first consider a one-dimensional SNS junction con- The length L of the junction is defined by the superconducting contact separation. A set of bottom gates tune the local chemical potential. Our model treats the junction as a normal (N) cuboidal cross-section (green) with flat superconducting (S) leads (blue). (b) The discrete lattice structure of our model -highlighted in a section of the nanowire. The red spheres form the lattice sites in an effective tight-binding approximation. The length of the nanowire is controlled by the number of layers of the N-region. A potential U on the inner lattice points (yellow sites) confines the particles to surface of the nanowire. The transverse square cross-section is 60 nm wide.
sisting 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, (1) H 0 is the single-particle effective Hamiltonian, ψ σ is the field operator with spin index σ ∈ {↑, ↓}, 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 µ.
We assume an identical effective mass in the N and S regions thus neglecting the Fermi wave-vector mismatch at the interface. ∆(z) is the superconducting order parameter along the junction, which we assume to be constant with jump-discontinuities at the N/S interfaces where ϑ(x) is the unit step function at x = 0, χ L,R is the superconducting phase of the left and right leads respectively, and φ = χ L − χ R is the phase difference.
In the ψ † ↑ (r), ψ ↓ (r) Nambu basis, we have the BdG equation The device is divided into three parts -a normal semiconductor section with a length L extended over z ∈ [0, L], and semi-infinite superconducting contacts extending to z = ±∞ on either side of the semiconductor (Fig. 2).
FIG. 2. Schematic of the SNS junction with a semi-infinite superconducting leads and an N device region. The length L of the nanowire is given by the number of lattice points (n) in the N-region [L = n · a], where a is the effective lattice constant. α and β are the tight-binding onsite and nearest neighbour coupling terms respectively. ∆0 and χ is the magnitude and phase of the superconducting order parameter.
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 where, t = 2 / 2m * a 2 is the nearest neighbour tight binding hopping parameter. The hopping matrix is given by This is not an atomic model, but an "effective" discrete lattice description of the junction. Parabolic dispersion relations correspond to the parameter space µ t. Typical experiments are setup in the µ ∆ 0 (Andreev approximation) regime. Hence, the choice of the effective lattice parameter a is bound by the inequalities The device Hamiltonian is subsequently written as where c † i is the creation operator of the Nambu spinor ψ † ↑ (r), ψ ↓ (r) 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

A. 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 (φ) of the leads. The retarded Green's function in the energy domain is given by where E denotes the energy, I is the identity matrix and η is an infinitesimal real constant. The Hamiltonian H is given by Eq. 10. The self-energy terms Σ 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 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 (φ) of the order parameter of the contacts and is shown in Fig. 3. The parameters for this computation are consistent with the Andreev approximation 58,60-62 (µ ∆ 0 ). As discussed in Appendix B, the breakdown of this approximation is manifested as an avoided level crossing in the ABS spectrum.
FIG. 3. Andreev bound state spectrum as a function of the superconductor phase difference for a clean, short 1dimensional SNS junction. The junction is tuned into the Andreev approximation regime with µ = 30∆0.

B. Current phase relationship
The current-phase relationship (CPR) links the macroscopic current flow in the junction to the phase gradient of the superconducting order parameter 63,64 . 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 (|E| < ∆ 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 currents 59 .
By contrast, when using the NEGF formalism the current-energy density can be computed at contact i, as a function of the phase difference φ using the current operator 38,39 where f (E) = 1/ (exp (E/k B T ) + 1) is the Fermi-Dirac occupation probability for a given energy level and k B is the Boltzmann constant. G r(a) and Σ r(a) i are the retarded (advanced) Green's function and contact i selfenergy respectively. To incorporate the opposite charge of electrons and holes we use the Pauli-z operator (τ z ) in the particle-hole Nambu space. This current operator is reviewed in the Appendix A. The total current at a phase difference φ is then given by Figures 4(a) and 4(b) compares the current phaserelations for a short (L < ξ 0 ) and long junction (L > ξ 0 ) respectively, as calculated from ideal scattering theory and NEGF. Here, ξ 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 junction 53,54,58,59 . With this assumption of no normal reflections, the bound states in a clean junction cross at φ = π. Hence, there's a discontinuity at φ = π in the CPR calculated using this method. The NEGF result is expected to match scattering theory exactly in the µ ∆ 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 threedimensional nanowire. The junction is along the zdirection 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 FIG. 4. The total current-phase-relation in an SNS junction computed from the NEGF current operator (orange) is compared with the ideal scattering theory (purple) calculations (see text). The CPR is plotted in the (a) short channel (L=40 nm) and (b) long channel (L = 400 nm) limit, chemical potential µ = 30∆0 and healing length ξ0 = 222 nm. We recover the expected saw-tooth profile for the total current using the NEGF formalism.
where z is the longitudinal direction, U is the surface confinement potential, and H T is the Hamiltonian of the transverse modes. For a cylindrical nanowire, the rotational symmetry about the longitudinal axis results in angular momentum ( ) subbands. This is because whereL z is the angular momentum operator in the zdirection. Hence, is a good quantum number. The subbands are eigenstates of theL z operator, labelled by their eigenvalueL 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 B, we compute the average angular momentum of each transverse subband, and observe only a small difference (see Fig. 5) from the unperturbed quantized eigenvalues ( ). 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 Within the tight-binding approximation, the on-site and hopping matrices in the particle-hole Nambu space are given by where Φ 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 (α N/S ) and contributes a phase to the hopping term corresponding to the gauge chosen for the vector potential (Eq. 18).

A. 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.023m e 65 (m e is the bare electron mass), radius R = 30 nm and chemical potential µ = 5∆ 0 . As described in Appendix B, the bound states anti-cross at φ = π due to normal reflections at the N/S interfaces. A normalised flux of Φ = 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 ± subbands. Here, Φ = B·S/Φ 0 , where S is the cross-sectional area, and Φ 0 = h/e is the flux quantum.
The chemical potential is adjusted to populate three subbands. In Fig. 5 we see that the = 0 subband remains unaffected, while the = 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 B. The characteristic length-scale associated with an occupied subband is called the healing length (ξ ) 54,59 , and is given by v F, is the fermi-velocity and is given by We have a "short junction" when the nanowire is shorter than the healing length (L < ξ ). 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 angularmomentum subbands using the correlation Green's function (G n ).

B. Average angular momentum of the transverse subbands
The angular momentum of the subbands can be computed as the expectation of theL 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 particlehole density per-unit energy. Using this, we find the expectation of theL z operator as a function of energy where max E {Tr [G n (E)]} returns the peak subgap quasiparticle concentration. In the Cartesian coordinate system, theL z operator is written aŝ The position operatorsx,ŷ are diagonal in the tightbinding 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 The linear momentum operators can be written aŝ where the subscript i is used to denote x, y basis of the position and momentum operators. Using Eq. 25, 26, 27, we write theL z operator and using Eq. 24 we find the expectation ofL 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 ( = 0, ±1). The average angular momentum of the subbands at Φ = 0.01, φ = π is shown in Fig. 5(c). Next, we increase the axial magnetic flux through the nanowire to Φ = 0.03. The = ±1 subband states further split ( Fig. 5(b)), and the angular momentum of the subbands is plotted in Fig. 5(d) for φ = π. The zero angular momentum subband ( = 0) has no azimuthal motion, and hence is unaffected by the axial field. A non-zero L z for the = 0 subband (Figs. 5(c), 5(d)) results from its hybridisation with the nearby = −1 subband. This consequently decreases the L z magnitude for the = −1 subband w.r.t. = +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 D, 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.

C. Zero angular momentum subband in the Superconductor
The BdG Hamiltonian in a superconductor is given by For a cylindrical geometry with an azimuthal vector potential, H 0 is given by As discussed, the radial confinement due to the nanowire's cylindrical geometry gives rise to angular momentum subbands labelled by . We will analyze the eigenenergies of this superconductor in the presence and absence of a magnetic field. Using the ansatz wavefunction exp(ik z z) exp(i θ), the diagonal elements of the BdG Hamiltonian can be written as The BdG Hamiltonian simplifies to and its eigenvalues E are given by This is the well-known superconductor dispersion relation, with a gap of ∆ 0 on either side of the fermi level.

Constant axial magnetic field, B = Bzẑ
In the Coulomb gauge we can write the vector potential for this magnetic field as From Stoke's law Exploiting the symmetry of the cylindrical geometry, the above equation can be simplified to Using the same ansatz exp(ik z z) exp(i θ), the diagonal elements can be written as for the electron and hole parts respectively. The effective chemical potential µ is defined in Eq. 30, and the field-coupling term E = The BdG Hamiltonian can then be written as and the eigenvalues E are given by 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 contacts 34,35 which naturally support only the = 0 subband. The geometry of the superconducting contacts implies a large = 0 subband energy, and can safely be assumed to remain unpopulated. This difference in geometry prompts the inclusion of = 0 subbands in the nanowire, and their exclusion in the superconductor. For = 0, we have E = 0, and thus the superconducting gap stays centered at the fermi level. This is illustrated in Fig. 6 -when = 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 = 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 subband in the N-section is given by The electron and hole wavenumbers can then be written as a function of energy (E) D. The 1-dimensional effective subband model As outlined above, it is important to ensure that we only have the = 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 µ , and a field-coupling term E .
The tight-binding Hamiltonian of the nanowire can be written as where N is the number of subbands. (For example, N = 3 =⇒ ( = −1, 0, 1)) and I n is the n × n identity matrix.
Meanwhile, the Hamiltonian of the contacts takes a similar form with IV.

A. 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 = 0 subband populated, there is no phase shift in the ABS, and the CPR retains its saw-tooth shape with a maximum near φ = π. 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 Φ = 4.04 the band depopulates min µ e , µ h = 0 and the current falls to zero. We observe in Fig. 7(a) that the critical current does not monotonically decrease to zero, particularly for Φ ∈ [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 (| | ≤ 1). The magnetic field evolution of the CPR is plotted in Fig. 8. Since the ABS for the | | = 1 subbands are phase shifted in presence of a flux ( Fig. 6(b)), the total current is no longer maximum near φ = π. The current in the junction is the sum over the individual subband currents, and consequently, the fluxdependent 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 e −k 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 Φ 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 | | = 1 subbands depopulate at Φ = 3.04. At zero flux, the CPR of each subband is maximum near φ = π 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 | | = 1 subbands pickup a phase and the subband CPRs no longer interfere constructively. Consequently, the critical current decreases with flux. At Φ = 0.72, the = 0 and = ±1 subband CPRs are maximally out-ofphase, resulting in a local minima (node). The subband and total CPR at this node is plotted in Fig. 8(c). Beyond Φ = 0.72, the critical current switches phase from φ < π to φ > π and the current increases again. This increase persists till Φ = 1.08 at which point the cur- rent is maximum near φ = 2π. At this flux, the | | = 1 subband current peaks near φ = 2π while it is negligible near φ = π. Hence, this secondary peak -which only involves contribution from | | = 1 subbands -is approximately a two-third of the primary peak and corresponds to a phase pickup of π 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 (| | ≤ 2). The critical current is plotted as a function of the magnetic flux in Fig. 9. Once again, at Φ = 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 | | = 1 and | | = 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.

B. Effect of disorder
In order to simulate experimentally relevant conditions, we include a random uncorrelated onsite disorder potential u ∈ [−W, W ] in the semiconductor. This models phase-coherent scattering events in the junction. We parameterise the disorder by the mean free path (λ mf ), which is estimated from the disorder-averaged normal state conductance (g) using the following relation N 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 φ = π at zero field), at a magnetic flux Φ * the | | = 1 subbands pickup a phase of π and their CPRs peak near φ = 2π. The = 0 subband retains its sawtooth CPR with a negligible current near φ = 2π. 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 φ = π at zero field. Thus, there exists no Φ * at which the | | = 1 subband current peaks while the = 0 subband current is negligible. As a consequence of the sinusoidal CPR, on picking up a phase of π the | | = 1 subbands destructively interfere with the = 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   FIG. 11. The critical current as a function of the applied magnetic flux for a clean and disordered (λ mf = 30, 80 nm) junction. Each critical current curve is normalized to its respective zero field value. We observe that the first node in the disordered junction occurs at a lower field as compared to a clean junction. 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.

C. 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 39,66 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 background 39 .
We subsume these processes within the NEGF formalism by including a self-energy term for the lattice background Σ r s , proportional to the Green's function and the emission-absorption coefficients. This calls for a selfconsistent computation of the Green's function and the bath self-energy 66 where × 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, 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 ϕ ∼ 100 − 300 nm [67][68][69] . In our model, l ϕ can be estimated from the statistical properties of universal conductance fluctuations (UCF) [67][68][69] The details of this calculation is presented in Appendix D, and the results are tabulated here in Table I The critical current oscillations in the presence of elastic dephasing interactions, as listed in Table I, are shown in Fig. 13(a). The green curve in Fig. 13(b) 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 pickup. 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(a), 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 spacing 34,35 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.  Table I. (a) The nanowire is free of disorder and only phasebreaking processes are involved. (b) The green curve corresponds to a nanowire with a random onsite potential distribution resulting in a mean free path λ mf = 30 nm, in addition to the phase-breaking processes with D0 = 1 × 10 −4 eV 2 . This is plotted in comparison with the corresponding disorder-free nanowire in (a) (purple curve). Each critical current curve is normalized to its respective zero field value. The simulations were performed for L = 160 nm, µ = 30∆0.

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 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 MBE 67 and metalorganic vapor phase epitaxy 68,69 grown InAs nanowire samples corroborates the dephasing in our model, as introduced by a uniform, spatially uncorrelated dephasing parameter (D 0 ).  The retarded Green's function in the energy domain is given by where H is the Hamiltonian, Σ r 1,2 are the retarded selfenergies of the semi-infinite contacts, and η is an infinitesimal real constant. The advanced Green's function is the Hermitian conjugate of the retarded Green's function (G a = G r † ). The Hamiltonian H is written in the tight-binding approximation Eq. 10. The surface Green's functions (g s ) at each contact are recursively evaluated where the subscript L, R labels the contact. Using this we compute the self-energy where σ 1 = β † g sL β, and σ 2 = βg sR β † 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 Γ i .
The Fermi functions in the particle-hole Nambu space is given by where f (E, µ) = 1/ (exp [(E − µ)/kT ] + 1) 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 The lesser Green's Function is then computed Next, we construct the current operator 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.
This can be also be written as where τ z is the Pauli operator in the particle-hole Nambu space. The total current is then evaluated by integrating the current-energy density 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, Pre-multiplying by G r and post-multiplying by G a , Multiplying the fermi-function, Thus, However, if we don't consider the term proportional to the infinitesimal η we end up with G < from Eq. A18 misses a term in the current proportional to (G r G a ), 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 - Substituting G n = Af and Σ in = Γf , the current operator can be simplified to  58,[60][61][62] . This approximation deals with a regime where the chemical potential of the nanowire is much larger than the superconducting order parameter of the leads (µ ∆ 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 ∆ 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, ∆ 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 µ ∆ 0 is not valid, even for a clean junction.
We now examine the spectrum for the chemical potential µ comparable to ∆ 0 . Figure. 14 plots the ABS spectrum for µ = 0.5∆ 0 .
FIG. 14. Andreev Bound state spectrum as a function of the superconductor phase difference for a clean 1-dimensional SNS junction. The junction is tuned out of the Andreev approximation regime with µ = 3∆0. The avoided level crossing is due to the normal reflections at the interfaces which couple the leftward and rightward moving quasiparticles.
We observe an avoided level crossing when µ ∆ 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 cross-ing at φ = π. 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(a), we plot the ABS energy gap (δ π ) at φ = π as a function of µ/∆ 0 . It is evident from this plot that δ π decreases as µ/∆ 0 increases. Thus, in the Andreev approximation regime (µ/∆ 0 → ∞) the bound states cross (δ π → 0) at φ = π. When ∆ 0 is kept constant and µ is varied, the gap varies as (µ/∆ 0 ) −1 for µ > ∆ 0 . In Appendix C we verify this dependence analytically by taking a scattering theory approach. Figure 15(b) plots the variation of the gap at φ = π with the nanowire length. These oscillations result from the interference of the waves reflected at the two S/N interfaces.  15. (a) The Andreev bound state energy gap at a superconductor phase difference φ = π is plotted as a function of the ratio of the chemical potential and the order parameter µ/∆0. The states will cross at φ = π as µ/∆0 → ∞ (b) The variation in the Andreev bound state energy gap for a fixed ratio of chemical potential to order parameter µ/∆0 at the superconductor phase difference φ = π as a function of the nanowire length L. These oscillations arise from the interference of the normally reflected quasiparticles at the two N/S interfaces. 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 π. We show this to be non-zero for a finite µ/∆ while approaching zero in the limit µ/∆ → ∞, where µ is the chemical potential of the entire device and ∆ is the induced superconducting gap.

Setting up the problem
We follow the notation from Kulik 59 . 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 ∆/µ 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 reflection 58,60-62 , but this is only true to the zeroth order in ∆/µ. This normal reflection component will modify the Andreev bound state quantization condition 53,59 where 1 γ = E ∆ +i 1 − E 2 ∆ 2 , k 0 and k 1 are the electron and hole wave vectors and χ = χ 2 −χ 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 The terms proportional to B , C and D are ignored when perfect Andreev reflection is assumed. We have included a ∆/µ coefficient to emphasize a first order expansion beyond the perfect Andreev reflection scenario.
Equating the coefficients at z = ±d/2, we get: Solving for C and D (to be compared with Eq 2.13 from Kulik), we get: We now divide the above equations and keep terms to first order in ∆/µ, 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: where the exact form of and η can be derived from equation C9.

Expression for the energy at χ = π
Let us tune the phase difference χ = π. Defining φ = cos −1 (E/∆), we can write γ = e −iφ . Using this relation, the quantization condition can be simplified to, Equating the real and imaginary parts, Since ∆/µ 1, under the small angle approximation Eq. C13 can be simplified to This can be further simplified using This is a finite energy for non-zero η and goes to zero in the limit ∆/µ → 0.

Comparison with numerics
Ignoring the term in the bracket from the final expression for E, we can simplify it to E ≈ ∆ 2 µ . We decided to observe the power law dependence of the gap using the numerical simulations. We considered two cases: 1. Fix µ and vary ∆: we expect to see a behavior   16. The dependence of the ABS energy at χ = π as a function of µ/∆ is investigated. δπ is the gap in the ABS spectrum at χ = π, and is thus twice the absolute value of the ABS energy. (a) µ kept constant and the ∆ is varied. We have µ/∆ on the x-axis (in log scale) showing the expected 1/x 2 scaling behaviour. (b) µ is varied and ∆ is kept constant. The expected 1/x scaling is observed for µ > ∆. the quasiparticles, and sets an energy scale for the problem. The phase coherence lifetime can thus be estimated as where Γ s,d is a diagonal element of Γ s = i Σ s − Σ † s is the broadening function corresponding to the lattice background. Note that Γ s , and hence the l ϕ so estimated from Eq. D1 is a function of energy. The phase relaxation length is then reported as an average over the energy grid, l ϕ = l ϕ (E) E . Under the assumption of ballistic transport, we estimate an upper bound on the phase relaxation length l ϕ < ∼ 150 nm for D 0 = 0.001 eV 2 , and for D 0 = 5 × 10 −4 eV 2 , l ϕ < ∼ 300 nm. The Fermivelocity decreases with field, and hence we observe a gradual monotonic degradation in the phase relaxation length. The estimated lϕ ∼ 150 nm, and lϕ ∼ 300 nm for (a) and (b) respectively. The monotonic decrease in lϕ is due to the reduction in Fermi velocity with an applied field. The simulations were performed for L = 160 nm, µ = 30∆0.

Estimation from statistical properties of UCF
For a nanowire length comparable to l ϕ , the normalstate 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 δG where the average . is taken over the magnetic field B.
The average fluctuation amplitude about the mean conductance is quantified by the root-mean-square rms(δG) B = var(δG) B . The rms(δG) B decreases monotonically with D 0 . The phase relaxation length l ϕ can be estimated from the analysis of the autocorrelation function F of δ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. Assuming the phase relaxation length (l ϕ ) to be greater than the nanowire diameter d, we can extract l ϕ directly from the correlation field 67,68,70 where γ is a dimensionless prefactor depending on the transport regime. We work in the the dirty metal limit with γ = 0.95 67 .

a. 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 ∈ [−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 λ mf ≈ 18 nm. The magnetoconductance fluctuations are plotted in Fig. 18(a). The normalised autocorrelation F/F (0) of δG is shown in Fig. 18(b), and the extracted parameters are listed in Table II.