Selfconsistent 3D model of SN-N-NS Josephson junctions

We develop a quantitative model describing the distribution of the supercurrent density and density of states in SN-N-NS type Josephson junctions in three dimensions (S is a superconductor and N is a normal metal). The model is based on the self-consistent solution of the quasiclassical Usadel equations using the finite element method. We investigate the influence of the proximity effect on the properties of the junction as a function of phase difference across the structure for various spatial dimensions and material parameters of S, N metals. The results are consistent with analytical solutions in the thin N layer limit and show consistent behavior for a large range of junction parameters. The results may serve to design nanoscale Josephson junctions for use in superconducting digital circuits.


Introduction
One of the important problems in the development of superconducting electronics is the design of Josephson junctions, which are nonlinear elements of superconducting circuits. These junctions should have high enough critical parameters such as I c and I C R N product (where I c is a critical current and R N is a normal resistance of a junction). In this respect, the reduction of the size of Josephson junctions is an important problem. Going to nanoscale, one naturally encounters the question about the local distribution of currents and electronic spectra in a junction. This problem was not properly addressed theoretically up to now in three-dimensional geometries. * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Planar SN-N-NS junctions (S stands for a superconductor, and N for a normal metal) with variable-thickness bridge geometry are currently considered as promising elements of superconducting circuits, which are scalable down to nanometer size and have high enough values of critical parameters. For a proper understanding of the properties of these structures, theoretical models are invaluable. These models should provide both an explanation for experimental results and good estimates for tunable parameters. The vast majority of the models developed in the past are based on the quasiclassical Green's functions method [1,2]. So far, these models are limited to effective one-dimensional models under assumptions on the junction geometry and interface parameters to enable analytical or numerical solutions. However, these quasi-onedimensional models are not suitable to describe situations that are realized in actual superconducting circuits.
For junctions consisting of multiple layers, such as SN or superconductor-ferromagnet (SF) type junctions the junction geometry can be quite complex and the interfaces between the layers cause complicated proximity coupling. These types of junctions can only be described by one-dimensional models in a few limiting cases. For a more general analysis of these types of junctions higher-dimensional models are required. There are also physical phenomena that cannot be studied in one dimension such as magnetic vortices and topological excitations.
Recently, the interest in higher-dimensional models for superconducting junctions has increased [3][4][5][6][7]. Numerical solutions of these models yield predictions for junction properties for non-trivial geometries and general values of interface parameters and allow for a better comparison with experiments than one-dimensional models.
In this paper, we present a three-dimensional model of the SN-N-NS Josephson junction in the dirty limit. Technically the problem is reduced to an effective 2D model, which we solve numerically on a 2D mesh using the finite element method to calculate the distribution of the current density and the density of states (DOS) inside the junction and compare our results with calculations done in the one-dimensional limit.

Model
The SN-N-NS Josephson junction consists of two superconducting electrodes of lengths (L − s)/2 and thickness d S , which are located a distance of s apart. The electrodes are connected by a thin normal metal strip of thickness d N and length L stacked on a dielectric substrate. The superconducting and normal materials both have a uniform width W. Far away from the weak link the superconducting electrodes go over into thicker bulk superconducting electrodes S ′ made of the same material. All materials are deposited on a dielectric substrate. The geometry is depicted in figure 1. In the coming model only the S and N layers are considered, the other layers enter only as effective boundary conditions. The axis origin lies at the top center of the normal metal strip. The center of the normal metal is positioned at y = −d N /2.
In experimental realizations of SN-N-NS junctions usually thin N films are used. There are two reasons for this choice. The first is to minimize the suppression of superconductivity induced into the N metal from its part located far from the SN interface. The film is in the normal state since superconducting correlations decay along the y-axis on the scale of the coherence length ξ N . The second reason is that the suppression of superconductivity in the S electrodes is smaller for thinner N layers. It is for this reason that the typical thickness of the N film is of the order of ξ N or smaller. The electron mean free path in thin films is limited by the film thickness.
Therefore, we can assume that all materials satisfy the dirty limit condition l j ≪ ξ j , where j = N, S denotes either the normal metal strip or the superconducting electrode, l j is the mean free path and ξ j the superconducting coherence length. Furthermore, we assume that the width W of the junction is much smaller than the Josephson penetration depth λ J so that magnetic effects inside the junction can be neglected.
Under these assumptions the junction can be described in the framework of the two-dimensional Usadel equations [8], which have the form: HereĜ j is the 2×2 matrix Usadel Green's function in layer j, i is the imaginary unit, andτ 3 and∆ are the third Pauli matrix and the pair potential matrix respectively: Furthermore E is the energy, D j the diffusion coefficient in layer j D j = v F l j /3, with v F the Fermi velocity. ∆ is the pair potential with ∆ * denoting it is complex conjugate. Inside the normal metal strip ∆ vanishes.

Stationary equations
In order to determine stationary properties of the junction, such as current densities, the Usadel equations need to be solved for specific imaginary energies E n = −iω n , determined by the Matsubara frequencies ω n = πT(2n + 1). For these energies it is convenient to parametrize the Green's functionĜ using the Φ-parametrization [9]. In the superconducting layers this turns (1) into: while in the normal metal strip it satisfies: Here G j = ω n / √ ω 2 n + |Φ j | 2 and Φ j are the normal and anomalous Green's functions in layer j, t = T/T c and ξ = ξ S /ξ N is the ratio between the coherence lengths in the superconductors and the normal metal, defined by ξ S = (D S /2πT c ) 1/2 , ξ N = (D N /2πT c ) 1/2 . Furthermore all energy quantities are normalized to πT c and all lengths are normalized to ξ N . Equation (4) is the self-consistency equation and needs to be solved together with (3) to determine the pair potential ∆.
The Usadel equations must be supplemented by proper boundary conditions at all boundaries of the domain to provide a unique solution. Far away from the weak link at the boundaries Γ + and Γ − the superconducting electrodes go over into the thicker bulk superconducting electrodes S ′ with a uniform phase. At this boundary the phase and magnitude of the Green's functions can thus be fixed: where ∆ 0 is the Bardeen-Cooper-Schrieffer (BCS) energy gap [10] and φ is the phase difference between the bulk electrodes. These boundary conditions are applicable especially whenever the junction is used to close a superconducting ring such as in SQUID structures. At the boundaries Γ SN between the superconductor and normal metal the Kupriyanov-Lukichev boundary conditions [11] can be applied to ensure continuity of the supercurrent flowing through this boundary: with the interface parameters: where R B is the specific boundary resistance of the SN boundary, ρ j is the resistivity of layer j and ⃗ n denotes the outward normal at the boundary. At all other boundaries the boundary conditions ensure zero current flow: From the solution of the Usadel equations the current density inside the junction can be calculated: This allows to express the current density in terms of the normal layer resistance: with R N = ρNs dNW .

Spectral equations
To determine the energy-dependent properties of the junctions the Usadel equations need to be solved for real energies E. In this case it is more convenient to parametrize the Green's functionĜ using the θ-parametrization [12,13], which transforms equation (1) for a superconductor into: The equations in the normal layer are equal to (12), but without the factor ξ and with zero pair potential. These equations in the S and N layers need to be solved simultaneously to determine the Green's functions θ j and χ j , which in turn determine the local DOS N(E,⃗ r) = Re{cos θ j }. The boundary conditions to the spectral equations can be determined in a similar fashion as for the stationary equations. Equation (6) turns into the conditions: The Kupriyanov-Lukichev boundary conditions take the form: Lastly, the no-current boundary conditions are given by:

Current density
Using the finite element method equations (3)-(9) were solved self-consistently on a 2D mesh representing the SN-N-NS junction. Details of the numerical method can be found in the appendix. The current density at each point was calculated via (10). An example of the distribution of the current density throughout the junction is shown in figures 2(a) and (b). It can be seen that the current density is symmetric and varies strongly throughout the junction. It is the highest inside the weak link region at the center of the normal metal. It peaks at the edges of the weak link region due to the sharp angle between the end of the superconductors and the normal metal strip. Inside the superconductors the current density is negligible compared with the current density inside the normal metal layer, even at high values of γ. By virtue of (7) the ycomponent of the current density is continuous over the SNinterface, ensuring continuity of the current perpendicular to the interface.
In the case d N ≪ 1, γ ≪ 1 the SN-N-NS junction was transformed into an effective 1D problem in the normal metal strip and was solved analytically in [14] in the limit s ≪ 1 for the cases of a small and large parameter γ BM .
We extended these results to two dimensions by including the effect of finite normal layer thickness d N and proximity effect parameter γ on the critical current. The effect of these parameters on the temperature dependence of the critical current is shown in figures 3(a) and (b). The current density decreases monotonically, both as a function of d N and of γ. The influence of d N is stronger than that of γ. Figures 4(a) and (b) show the dependence of the current density on d N and γ. The current density decreases approximately exponentially, both with d N and γ. In the limit d N − → 0, γ − → 0 the results in [14] are recovered.

Density of states
In analogy to the procedure in [14] for the Φ-parametrization, equations (12)- (15) can be reduced to effective 1D equations inside the normal metal layer in the case d N ≪ 1, γ ≪ 1: Subject to: Here γ BM = γ B d N . In the short junction limit s ≪ 1 the last term in (16) can be neglected and the general solution in the region 0 ≤ x ≤ s/2 can be obtained analytically [15], The constants χ 0 , j E and α are determined by imposing the boundary conditions at x = 0 and requiring continuity of the solution at x = s/2: At the center of the junction the DOS is then: In the case of zero phase difference the DOS is independent of position and this expression reduces to: It follows from (31) that the DOS contains two singularities.
One at E = ∆ 0 and another at E/∆ 0 = z 0 , with: Here β = γ BM ∆ 0 . The DOS is shown in figure 5. For nonzero values of γ BM the singularity in the DOS at E = ∆ 0 splits into two peaks and the energy gap is lowered. In agreement with the results obtained in [16] the singularity at E = ∆ 0 becomes increasingly sharp as γ BM is increased. For γ BM ≪ 1 the energy gap decreases approximately linearly with increasing γ BM . The energy gap is thus proportional to the transparency of the SN-boundary. This behavior of the DOS has also been derived for other types of SN-junctions [16,17]. When a phase difference over the junction is introduced the DOS inside the normal metal strip becomes positiondependent. The two peaks in the DOS for finite γ BM respond differently to an increase of phase as seen in figure 6. For a nonzero phase difference the singularity at E = ∆ 0 turns into a local maximum, whose height decreases with φ. The other singularity remains independent of φ, but the singularity will shift to smaller energies, eventually disappearing at φ = π when the DOS returns to its non-superconducting value. The magnitude of the energy gap thus decreases as φ goes to π.
For perfect transparency of the interfaces (γ B = 0) the dependence of the energy gap at the center of the normal metal on the phase difference is E g = ∆ 0 cos φ/2, which corresponds to the expression previously found in [15]. Even at small nonzero values of γ B the dependence of the singularity on φ seems to approximately satisfy this relationship.  The above results were extended using a 2D model based on the finite element method to take into account larger values of d N , γ and s. With this model the DOS can be calculated not only in the normal metal, but also in the superconductors. Figure 7 shows the distribution of the DOS inside the junction. The energy was chosen below the BCS energy gap to ensure that a bulk superconductor would have zero DOS. Far away from the weak link region the superconductor indeed has zero DOS, but due to the inverse proximity effect with the normal layer the DOS increases close to the weak link.
The dependence of the DOS on the phase difference and energy at the center of the normal metal is shown in figure 8. For small phase differences the two peaks in the DOS at E = ∆ 0 and E ≈ z 0 ∆ 0 are clearly visible. Both peaks become increasingly more narrow as the phase difference increases and vanish at a phase difference φ = π. The magnitude of the energy gap depends on the phase difference as well, having its maximum at zero phase difference and closing completely at φ = π. Even though γ = 1, and the short junction limit condition is not satisfied, the shape of the energy gap still resembles the limiting behavior of cos(φ/2) well.

Discussion and conclusion
In this paper, we studied the behavior of the current density and DOS in the SN-N-NS junction. We have shown that these quantities depend strongly on the junction dimensions, phase difference and interface parameters.
Using the finite element method we have extended a previous one-dimensional model of the SN-N-NS junction to two dimensions to obtain a physically more accurate model and extended it for DOS calculations.
In two dimensions, the critical current in the junction shows a strong resemblance to the one-dimensional case. Its maximal value is approximately inversely proportional to d N and even at large values of γ. The critical current density of the SN-N-NS junction is mostly determined by the properties of the normal metal and the transparency of the SN interfaces.
The DOS inside the junction varies strongly inside both the N, and the S layer as a function of position and energy. By virtue of the proximity effect a subgap in the DOS is induced in both layers comparable to that found in one-dimensional models. Even at small values of the parameter γ this subgap persists in the superconducting layer close to the normal metal.
Our model can be used for accurate design of nanoscale Josephson junctions based on the desired distribution of currents and electronic spectra in the junction. In particular, it can be used as a first step to address the problem of power dissipation in junctions with high current densities.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.

Acknowledgments
M Yu K thanks for support from Russian Science Foundation (Project No. 20-12-00130) and is grateful to the Interdisciplinary Scientific-Educational School of the Moscow State University 'Photonic and Quantum Technologies. Digital Medicine.'

Appendix. Numerical method
The Usadel equations describing the SN-N-NS junction were solved numerically in both 1D and 2D using a finite element method written in MatLab. To apply this method the Usadel equations and their boundary conditions are written in an equivalent weak formulation. For the stationary equations this is: Here S is the domain of one of the superconductors and η S is an arbitrary test function that vanishes on Γ ± . The formulation in the normal metal layer is similar. For the spectral equations the weak formulation is given by: The computational domain is represented using a mesh of triangular elements. The Usadel equations can be transformed into a nonlinear system of equations by means of a standard Galerkin method using Gaussian quadrature rules [18]. In the Galerkin method the Green's functions are approximated by an expansion into basis functions. The main difficulty for applying the Galerkin method is that in general the Green's functions are discontinuous over the internal boundary Γ SN .
To resolve this we approximate the solution separately in each layer by piecewise first order Lagrangian basis functions and match the solutions in the different layers at Γ SN using the Kupriyanov-Lukichev boundary conditions. For the solution of the nonlinear system a Newton-Raphson method with line search [19] was used, which provides both fast convergence and good stability even when the solution is close to singular.
In addition to the Usadel equations the self-consistency equation needs to be solved for the pair potential ∆. This equation can be solved using the Picard iteration method: where k is the iteration number. The convergence of this method is however quite slow. To improve convergence we applied Anderson's acceleration method with restarts to the self-consistency equation [20][21][22][23]. The above described iterative process was repeated until both the pair potential and the Green's function had converged to within a tolerance of 10 −5 . To check the converge of our method we determined the dependence of the numerically calculated current density on the mesh size as shown in figure A1. It can be seen that mesh refinement reduces the numerical error linearly as expected and that even for mesh element sizes around 0.1ξ 2 N the numerical error remains smaller than 2%.