Current Conservation in the Self-Consistent Josephson Junction

TL;DR

This work shows that current conservation in Josephson junctions necessitates a fully self-consistent order parameter within the Bogoliubov–de Gennes framework. By introducing a phase-gradient method that imposes a constant phase gradient in the superconducting leads and self-consistently matching lead and junction currents, the authors obtain a current-conserving solution for quasi-1D SNS junctions. The resulting BdG spectrum and current–phase relation depart significantly from fixed-phase treatments, exhibiting backward skewness and lead-induced Doppler effects, with behavior dependent on junction length and gate voltage. The approach explains experimental observations in nanowire systems and provides a versatile, extensible framework for studying current conservation in more complex geometries and finite-temperature regimes.

Abstract

Conventional treatments of Josephson junctions (JJs) are typically not current-conserving. In the mean-field BCS theory, current conservation is only guaranteed if the superconducting order parameter is treated self-consistently. We show that this requirement has significant consequences for the current-phase relation (CPR) in certain regimes, where the current density in the superconducting leads is non-negligible. To this end, we introduce a numerical method for the self-consistent treatment of the BdG equations with current conservation for quasi-1D superconductor-normal (metal)-superconductor (SNS) JJs. Our model incorporates a phase gradient of the order parameter in the leads, which is set to match the Josephson current through the weak link. We compare our method to standard, non-current-conserving approaches by calculating the CPR for SNS JJs while varying lengths and gate voltages controlling the normal metal. We show that current conservation has significant implications for the Josephson harmonics and can weaken or even reverse forward skewness of the CPR.

Figures (7)

• Figure 1: (a) Schematic of a Josephson Junction with $L$ and $d$ the lengths of the superconducting leads and normal (insulating or metallic) weak link, respectively. The transverse extent of leads and weak link are identical, which implies that the current density should match. (b)--(e) Solid blue lines represent the phase-gradient-matching solution; dashed orange lines represent the fixed-phase method, and the dash-dotted green lines the square-well model for $\Delta$. The phase drop over the junction is $\delta\varphi=\pi/4$. (b) Absolute value of the order parameter $|\Delta(x)|$. (c) Phase profile of the order parameter $\varphi(x)$. (d) Current $j(x)$ throughout the system. (e) Current source $S(x)$ (the rhs of Eq. \ref{['eq:current_continuity_equation']}).
• Figure 2: BdG spectrum, dispersion, CPR, and Fourier components of the CPR, as obtained with the (top row) fixed-phase and (bottom row) phase-gradient methods (see text for details). (a) The BdG spectrum for varying phase drops over the junction. The Andreev states are identified by their localization within the normal metal, and are colored gray--blue for low and high degree of localization respectively. The Andreev states are blue. (b) The dispersion relation (Eq. \ref{['eq:dispersion']}) with the Andreev bound states displayed as blue dots. The plots clearly show the raising and lowering of the Andreev states at $\pm k_F$. The dispersion is drawn for $\delta\varphi = \pi/4$ as indicated by the vertical line in (a). (c) The CPR with the dashed line showing a sinusoidal CPR as reference. (d) Fourier components of the CPRs. The phase gradient in the leads strongly affects the Josephson harmonics. The length of the junction is $d \approx 2.2 \, \xi$.
• Figure 3: CPRs for junctions transitioning from the short junction ($d<\xi$) to the long junction regime ($d>\xi$); $\xi \approx 10$ sites. (a) CPRs for the phase-gradient and fixed-phase methods. (b) The evolution of the CPR maximum $\delta\varphi_{\mathrm{max}}$ with junction length. The colored dots represent numerical results, and the solid and dashed black lines derive from Refs. sonin2025theory and thuneberg2024Squarewell respectively. (c) The maximal values reached by the CPRs plotted against the junction length. The phase-gradient and fixed-phase methods agree remarkably well despite their difference in shape.
• Figure 4: (a) CPRs of the phase-gradient method with varying gate voltages applied (offset according to $V_G$). The black line tracks the maxima of the CPRs. The CPRs in the shaded region display backward skewness. (b) CPRs of the fixed-phase method. All CPRs are forward-skewed. (c) First two Fourier components (imaginary parts), $-\mathrm{Im}[a_{1,2}]$, for the phase-gradient (solid) and fixed-phase (dashed) methods. The shaded area marks the $V_G$-interval where $a_1$ and $a_2$ are of the same sign, indicating the backward skewness obtained by the phase-gradient method. $eV_G = -1\mu$ corresponds to depletion and marks the crossover from the SNS to SIS regime. In calculating these CPRs, $d/\xi \approx 3$.
• Figure 5: Reproduction of Fig. \ref{['fig:SNS_CC']} with the boundary regions included (shaded orange).