Beyond the effective length: How to analyze magnetic interference patterns of thin-film planar Josephson junctions with finite lateral dimensions

The magnetic field dependent critical current $I_{\text{c}}(B)$ of a Josephson junction is determined by the screening currents in its electrodes. In macroscopic junctions, a local vector potential drives the currents, however, in thin film planar junctions, with electrodes of finite size and various shapes, they are governed by non-local electrodynamics. This complicates the extraction of parameters such as the geometry of the effective junction area, the effective junction length and, the critical current density distribution from the $I_{\text{c}}(B)$ interference patterns. Here we provide a method to tackle this problem by simulating the phase differences that drive the shielding currents and use those to find $I_{\text{c}}(B)$. To this end, we extend the technique proposed by John Clem [Phys. Rev. B, \textbf{81}, 144515 (2010)] to find $I_{\text{c}}(B)$ for Josephson junctions separating a superconducting strip of length $L$ and width $W$ with rectangular, ellipsoid and rhomboid geometries. We find the periodicity of the interference pattern ($\Delta B$) to have geometry independent limits for $L \gg W$ and $L \ll W$. By fabricating elliptically shaped S$-$N$-$S junctions with various aspect ratios, we experimentally verify the $L/W$ dependence of $\Delta B$. Finally, we incorporate these results to correctly extract the distribution of critical currents in the junction by the Fourier analysis of $I_{\text{c}}(B)$, which makes these results essential for the correct analysis of topological channels in thin film planar Josephson junctions.


I. INTRODUCTION
Planar Josephson junctions are ubiquitous in modern solid state physics research, with examples ranging from topological junctions [1][2][3], high T c (grain boundary) junctions [4,5], gated-junctions that control supercurrent flow [6,7], graphene-based junctions [8,9], magnetic field sensors [10][11][12] and, junctions with a ferromagnetic weak link [13][14][15]. A major tool in analysing these junctions experimentally is the magnetic interference pattern observed in the critical current (I c (B)), the shape and periodicity of which can reveal, using Fourier transform, information about the underlying distribution of critical current in the weak link [16]. Often this Fourier analysis is carried out in terms of an effective junction length, given, for macroscopic junctions, by 2λ + d, where λ is the London penetration depth and d the thickness of the weak link. This effective length originates from the Meissner effect. However, when the junction is formed between two superconducting thin films, with a thickness below λ, the shielding currents running along the junction, responsible for the shape and periodicity of the magnetic interference of the critical current I c (B), are no longer determined by the Meissner effect in its macroscopic form (i.e., by the local vector potential). Rather they are determined by non-local electrodynamic effects [17][18][19][20].
In numerous theoretical and experimental studies, it was found that in thin film planar junctions, I c (B) becomes completely independent of λ and is solely determined by the geometry of the sample [20][21][22][23][24]. Moreover, John Clem provided a method to calculate I c (B) for planar junctions that are also restricted in their lateral size (i.e., a Josephson junction separating a rectangular superconducting strip of width W and length L in two halves) [23]. As experimental studies often deal with finite-size geometries, his theory is highly topical at the moment.
This paper bridges the gap between predicting the I c (B) of thin film planar junctions featuring finite lateral geometry, and the correct analysis of the experimental interference patterns used to extract the current density distribution. First we review the technique proposed by Clem and extend on his work by covering two more geometries: the ellipse and the rhomboid. We calculate I c (B) for these geometries, extract the periodicity of the interference pattern (∆B) for different ratios of L/W , and find ∆B to have two geometry independent limits for L W and L W . By fabricating elliptically shaped S−N−S junctions with different ratios of L/W , we experimentally verify the geometry dependence of ∆B. Finally, we adapt the well-known Fourier relation between I c (B) and the critical current density distribution for use on laterally finite thin film planar junctions. We find that altering the Fourier transform is crucial for predicting the location of possible current channels in thin film planar junctions.

II. REVIEW OF THE CLEM MODEL
We consider a normal metal Josephson junction (dimensions W JJ and d) that divides a symmetric superconducting thin film, having dimensions L and W , into two halves. Figure 1 shows a schematic of three of such films, having different geometries. The junction, colored red in Figure 1, is running along the y-direction from −W/2 to W/2 (i.e., W JJ = W ). Since we examine the thin film limit, the screening current density is assumed uniform along the thickness of the film, which effectively reduces the problem to a 2D one. We specifically con- sider the junction to be in the short junction limit, as the model by Clem treats an infinitesimally thin insulating tunnel junction. Furthermore, it is assumed that the electrode the electrode dimensions are smaller than the Pearl length, given by: Where t film the thickness of the superconducting films. This implies that the self fields originating from the screening currents are far smaller than the applied external field. Additionally we assume that the junction is in the narrow limit, meaning that the junction is less wide than the Josephson penetration length, which for planar junctions in the thin film limit is the given by [20,21,23]: Here t junc is the thickness of the junction (not necessarily equal to the thickness of the film), I c (0) its critical current at zero magnetic field, µ 0 is the vacuum permeability, and Φ 0 is the magnetic flux quantum. In order to calculate I c (B), we assume a sinusoidal current-phase relation J x = J c sin ϕ(y), where ϕ(y) is the gauge-invariant phase difference over the junction,  Figure 1.
which depends on the location along the junction. It can be evaluated within the framework of Ginzburg-Landau theory by considering the second Ginzburg-Landau equation, which is given as: Here A A A is the vector potential corresponding to the applied magnetic field (B B B = ∇ ∇ ∇ × A A A), and γ is the gauge covariant phase of the wavefunction describing the superconducting order parameter (given by Ψ = Ψ 0 e iγ [25]). Finally, θ is the gauge-invariant phase gradient (required by the fact that J J J is a gauge-invariant property). ϕ(y) is then given by integrating θ across the junction: In Figure 1d, we sketch a zoom of a junction, where we specify an integration contour under a magnetic induction of B B B = Bẑ. By integrating ∇ ∇ ∇γ along this contour and realizing that C ∇ ∇ ∇γ dl l l = 2πn, where n is an integer and sin (ϕ + 2πn) = sin (ϕ), we find: Here we have used Stokes theorem to evaluate the flux entering the contour and used the fact that the electrodes are mirror symmetric (J y ( d 2 , y) = -J y (− d 2 , y)). For macroscopic junctions J y,R ( d 2 , y ) = Bµ0 λL resulting from the the Meissner effect, leading to ϕ(y) = ϕ(0) + 2π(2λ+d)B Φ0 y, where we recognize the effective junction length. Since the junctions considered here are in the thin film limit, we take a different approach in evaluating (a) Gauge-covariant phase simulated in the right electrode for a disk-shaped planar Josephson junction, normalized to the applied magnetic field and width of the junction γΦ0/BW 2 . The junction is shown as a green line. This result allows for extracting the gauge-covariant phase along the junction. It follows the scaling of Eq. 14, and it is determined by a dimensionless function, which is plotted in (b). (c) Shows the interference pattern calculated using the result in (a) by numerically evaluating Equation 12 for different values of B. The typical interference pattern looks like a Fraunhofer pattern at first sight. However, the peak height decreases less strongly than 1/B, and the width of the side lobes is larger than half of the middle lobe, which is 10.76 mT wide. Furthermore, the width of the nth side lobe increases and reaches an asymptotic value for large values of n, which is evident from the inset of (c), where we plot the width of the nth side lobe. The width of the fifth side lobe is used for comparisons between simulations and experiments. J y ( d 2 , y ). First note that the supercurrent is conserved and therefore ∇ ∇ ∇ · J J J = 0. By choosing the convenient gauge A A A = −yBx, we find ∇ × A A A = Bẑ and ∇ · A A A = 0. Therefore, the divergence of the second Ginzburg-Landau equation (Eq. 3) reduces to: Therefore, we mapped the second Ginzburg-Landau equation onto the Laplace equation. With sufficient boundary conditions, it can be solved for a unique solution, which allows us to calculate J y ( d 2 , y). The boundary conditions arise from the prerequisite that no supercur-rent can exit the sample at its outer boundaries. Furthermore, we assume a weak Josephson coupling, meaning that the shielding currents in the electrodes are far larger than the Josephson currents between the electrodes, which we approximate as J x ( d 2 , y) = 0. Therefore, we can write: Wheren n n R is the unit vector, normal to the outer edges of the right electrode. Combined with the second Ginzburg-Landau equation, this leads to a set of Neumann boundary conditions: Which is sufficient to solve for γ(x, y). Next, Eq. 5 allows us to find the gauge-invariant phase difference over the junction ϕ(y). Note that we have conveniently chosen A y = 0. We then find: Therefore, ϕ(y) is given by the simple expression: Next, the current across the junction is given by J J J dS S S, yielding: We assume that the critical current density at zero field is distributed uniformly over the junction, yielding J c = tjuncW . Also, note that ϕ(0) is independent of y and therefore merely is a phase factor. The critical current is reached if we current-bias the junction by setting ϕ(0) = π/2, from which follows: We see that finding I c (B) becomes equal to a boundary condition problem of solving the Laplace equation in the geometry of the electrodes. Indeed, the solution is completely determined by the geometry of the sample and is independent of λ.

III. COMPARING DIFFERENT GEOMETRIES
As it is not trivial to find a general analytical solution to the boundary problem of Eq. 6 for the ellipsoid and rhomboid geometries, we solve the Laplace equation numerically using COMSOL Multiphysics 5.4. We define the right electrode geometry in 2D, divided into a triangular grid. Crucial for correctly solving Eq. 6, is a grid size that is small enough to capture small changes in γ and, on the edges,n n n R . We found a maximum element size (i.e., the grid edge size) of 0.01 ln (1 + L/W ) nanometer to be a good compromise between computation time and precision. Using trigonometry we evaluate A A A·n n n R for each geometry and list the corresponding boundary conditions in Table I (here the numbering corresponds to the numbers in Figure 1). In the Appendix, we provide a full derivation of each of the boundary conditions.

A. Simulation results
Clem showed that the analytical solution for the rectangular geometry is an infinite series of sines and hyperbolic tangents [23]. For the rectangle, this leads to the maximum in γ d 2 , y to occur at W/2, which can be approximated as: FIG. 3. Dimensionless measure of the period ∆B (the width of the fifth side lobe) of the calculated interference pattern Ic(B) for the three geometries. In (a) we plot this value on log-log scale versus the aspect ratio L/W , in (b) it is plotted versus the total electrode area A (i.e, combined area of left and right electrode), scaled by the W 2 . Figure (b) reveals two limits for ∆B for L W and L W . The first corresponds to the limit of an infinite superconducting strip ∆B = 1.842Φ0/W 2 , whereas in the latter we find ∆B = 2Φ0/A. Contrary to ∆B, Ic(B) itself is not geometry independent in this limit. from Eq. 14, for the limit L W in (a) and L W in (b). The maximum of these functions is located at y = |W/2| and equals unity. Therefore, ∆B (large n limit of the nth side lobe of Ic(B)) is universal for these limits. However, for the limit L W , f y W is not geometry independent, which entails that Ic(B) is not geometry independent as well, in this limit.
Here ζ is the Riemann zeta function. Now we generalize this approximation to include the other geometries. We find that the simulated γ d 2 , y universally follows: Where f y W is a dimensionless function defined by the specific geometry and A is the total surface area of the electrodes (i.e, combined area of left and right electrode). Note that we have substituted L W in the argument of the hyperbolic tangent for A W 2 ; the reason for this choice will become apparent below when discussing the period of the I c (B)-pattern. Figure 2a shows the calculated γ(x, y) for a disk geometry, normalized to the applied magnetic field and width of the electrodes γΦ 0 /BW 2 . We plot f y W for this disk in Figure 2b. By evaluating the integral of Eq. 12 numerically for different values of B, we calculate the interference pattern of a disk-shaped junction ( Figure  2c). the pattern resembles a Fraunhofer pattern at first sight. However, the peak height decreases less strongly than 1/B, and the width of the middle lobe is not twice the width of the side lobes. In the inset of Figure 2c, we plot the width of the nth side lobe (∆B n ); the width increases and reaches an asymptotic value for large n.
In order to compare the interference patterns of junctions of different geometry, we define the period of the oscillations to be the width of the fifth side lobe (∆B = ∆B 5 ). In the inset of Figure 2c, this is shown by the vertical reference line. The width of the fifth side lobe is not only sufficiently close to the asymptotic value but also experimentally accessible without the need for large magnetic fields. We now compare the periodicity of the interference patterns for different geometries by plotting the dimensionless value ∆BW 2 /Φ 0 as a function of the aspect ratio L/W in Figure 3a on a log-log scale. First, we find the results obtained on the rectangular junction to match the analytical results obtained by Clem [23]. Furthermore, the periodicity of the pattern increases as the sample dimensions are diminished. Finally, we evaluated the width of the junction (d) to be irrelevant in determining ∆B. Specifically, its contribution to the period is in the µT range for realistic sizes of d. The consequence is that ∆B is determined by the maximum of γ, i.e., γ( d 2 , W 2 ). ∆B reaches asymptotic values for the limits L W and L W for all three geometries. The value of ∆B becomes geometry independent in these limits, as revealed by rescaling the results from Figure 3a to a A W 2 dependence, displayed in Figure 3b. In the first limit, L W , all three geometries become an infinite superconducting strip. Here we retrieve ∆B = 1.842Φ 0 /W 2 , which matches literature [22,23]. In this limit, we find γ d 2 , y to follow: is a dimensionless function running from -1 to 1, plotted in Figure 4a. In the other limit, L W , Eq. 14 reduces to:  Figure  2c, the middle peak is twice as wide as the neighboring ones. (c) depicts a top-view false colored electron micrograph of an ellipse-shaped junction. Again we indicate the notches with white arrows; the scale bar represents 1 µm. In (d), we plot the corresponding interference pattern as a dV /dI color map, which is used to extract the periodicity of the oscillations [26].
independent of the underlying geometry and equal to unity, we find a geometry independent period, where ∆B = 2Φ 0 /A. We can generalize this concept to find a general expression for ∆B: Note that max(f y W ) ≈ 1 for all ratios L/W , and thus Eq. 18 can serve as a good approximation for ∆B. Therefore, we justify the relation of Eq. 14 as it demonstrates the emerging universal limits where ∆B = 2Φ 0 /A and ∆B = 1.842Φ 0 /W 2 , as well as provides a good approximation of ∆B between the limiting cases.
Although ∆B is geometry independent in the limit L W , I c (B) itself is not universal in this limit. This is caused by the fact that f y W differs between geometries for y = |W/2| (see Figure 4b). For the rectangular geometry, for example, this function is linear in y: f y W = 2y W . Therefore, we retrieve the Fraunhofer pattern, where L eff = L/2 + d. The effective length equals the length of a single superconducting electrode plus the junction length. This can be understood by considering that the screening currents trace loops in the electrodes, that reduce to two parallel and opposite current tracks, when L W . γ d 2 , y in the rhomboid geometry is radically different; it is well approximated by a sine function: f y W = sin πy W . This leads to an interference pattern that is far closer to the pattern shown in Figure 2c, and not a Fraunhofer pattern. In conclusion: the shape and periodicity of the I c (B)-pattern for low magnetic fields is independent of ∆B, which is universal for L W .

B. Comparison to experiments
In order to verify the dependence on the geometry, we fabricate five ellipse-shaped planar S−N−S junctions for different ratios of L/W . Besides, we make a rectangularshaped junction with dimensions well in the L W limit.
First, a four-probe contact geometry is patterned on Si substrates using electron-beam lithography. Next, an Ag (20 nm), MoGe (55 nm) bilayer is deposited by sputter deposition. Subsequently, we use Focused Ion beam (FIB) milling to structure elliptical devices in the bilayer. By applying an ultra-low beam current of 1.5 pA, the weak link is formed by a line cut in the MoGe layer at the center of the device. This completely removes the superconductor on top, but leaves a normal metal connection. The resulting trench separates the MoGe electrodes by a roughly 20 nm weak link, allowing Josephson coupling in this S−N−S system. Similar junctions, featuring a ferromagnetic layer, were fabricated in this manner, to study the interplay between supercurrents and ferromagnetic spin textures [13,14,27]. Figures 5a and 5c show false colored electron micrographs of two of such devices, for L = W and L = 4W respectively.
Two corresponding interference patterns obtained on the samples in 5a and 5c are shown in Figure 5b and 5d. Clearly, the period of the interference patterns scales with L/W . However, we find that the middle peak is twice the width of the neighboring ones and the amplitude of the side lobes of the I c (B)-pattern feature a similar width, instead of the asymptotic behavior predicted by our theory (see Figure 2c). This can be explained by considering that l ≈ 100 nm (Eq. 2; based on λ = 535 nm [28]), which is small with respect to W . Our samples are therefore not in the narrow junction limit and allow Josephson vortices to stabilize in the junction. The width of the middle lobe can therefore not be predicted by our theory. However, Boris et al. have shown that ∆B n for large n follows the predictions of non-local electrodynamics [20]. Therefore, we can compare the measured ∆B = ∆B 5 to our theoretical model.
To compare the period of the I c (B)-pattern to our theory, we plot ∆B for all measured samples along with the calculated values in Figure 6. By the blue star symbol we also mark the periodicity of the Co-based S−F−S disk-junctions discussed elsewhere [14]. Although there is a constant offset between the measured periodicity and the calculated values, the overall trend is well predicted.
This constant offset is due to a trivial side effect of the FIB structuring method: some parts of the bi-layer (i.e., the edges of the device) mill faster than the bulk of the material. Consequently, notches develop on the side of the device when fabricating the trench. These notches make the width of the weak link (W JJ ) shorter than the width of the electrodes (W ), which can result in a constant offset between experiments and the simulations, where it is assumed that W JJ = W . In order to show that we reach the geometrical independent limit for L W , we have fabricated a bar-shaped sample with L/W > 10. In the Supplemental Material we present a scanning electron micrograph of this device accompanied by the interference pattern obtained on this sample [29]. versus the aspect ratio L/W . The blue star indicates the periodicity of the cobaltbased disk junctions discussed in reference [14]. Although we can predicted the L/W -dependence, we find a constant offset between the experimental values and the simulations. This is due to the notches visible in Figure 5a and 5c, which makes the actual junction width (WJJ) shorter than the width of the electrodes (W ). To further illustrate this, we plot using the open green star.
In this limit we expect ∆BW 2 JJ /Φ 0 = 1.842 [10,12]. By inspection of the scanning electron micrograph we have extracted W JJ for the bar-shaped sample, which leads to ∆BW 2 JJ /Φ 0 = 1.70. By the green open star symbol, we plot ∆BW 2 JJ /Φ 0 for the bar-shaped sample in Figure 6. The error bars correspond to a 20 nm uncertainty in the junction length. [30] Another method of accounting for the influence of the notches is modifying the Fourier relation between the critical current density distribution J(y) and the magnetic interference pattern I c (B), which will be discussed in the next section.

IV. FOURIER ANALYSIS OF THIN FILM PLANAR JUNCTIONS.
In their 1971 paper, Dynes and Fulton found a Fourier relation between the current density distribution of a Josephson junction and its magnetic interference pattern [16]. This method has been used widely the last years in analysing supercurrents planar Josephson junctions [1, 2, 6-9, 13, 14, 31-33]. However, the original Fourier relation is developed for macroscopic junctions where the screening currents are Meissner-based. This section will give a brief review of the Dynes and Fulton method and will adapt the Fourier relation for the use in thin film planar junctions, which is essential for correctly interpreting interference patterns obtained on such junctions.
First we write the current phase relation in Eq. 11 as a complex expression and extend the integration bounds to infinity, since J c (y) = 0, for y > |W JJ /2|: (19) Here ϕ B is the gauge-invariant phase difference over the junction due to the magnetic induction. The critical current is given by the absolute value of the complex expression. Note that this equal to setting ϕ(0) = π/2 in Eq. 12: From this equation a general expression for a Fourier transform can be recognized. For a junction with macroscopic leads discussed above, we have ϕ B (B, y) = 2π(2λ+d)B Φ0 y and therefore: Here we have defined the reduced field β = (2λ+d)B

Φ0
, such that the position along the junctions y and β form conjugate variables. For the mesoscopic devices discussed here, this quantity needs to be replaced by Eq. 10, yielding: Where we omitted the contribution from the weak link, as its magnitude is negligible. Specifying γ d 2 , y using Eq. 14, we can define a new pair of conjugate variables: the lengthỹ = W f y W and the reduced field [34], to arrive at: Where we made a change of coordinates andJ c is defined as:J Here the function g ỹ W is the inverse of f y W , or g ỹ Figure 5b, carried out using three different methods. In (a) we use the formalism for macroscopic junctions (following Eq. 21, where L eff = 2λ + d), whereas in (b), we make use of the simulation data shown in Figure 2b (following Eq. 22). We indicate the boundaries of the electrodes (−W/2 and W/2) by solid reference lines and the boundaries of the actual weak link (−WJJ/2 and WJJ/2) by dotted reference lines. Only the method based on the simulations of the shielding currents correctly predicts the uniform current density distribution, which is limited to the actual junction only. Finally, in (c), we carried out the Fourier analysis using a linear approximation of f disk y W , circumventing the need for rescaling the axes, yet retaining the correct Jc(y).

Equation 24
is a Fourier transform that includes a rescaling of the axes to retrieve the actual current density distribution J c (y).
In Figure 7 we compare three different methods of obtaining the current density distribution extracted by the Fourier analysis from the data obtained on the diskshaped sample shown in Figure 5b. Specifically, Figure  7a shows the current density distribution obtained using the method for macroscopic junctions (i.e., following Eq. 21, using L eff = 2λ + d) and Figure 7b shows the Fourier transform based on our phase difference calculations (Eq. 22). The solid reference lines indicate the width of the electrodes (i.e., the disk diameter W ) and the dotted reference lines indicate the width of the actual junction as measured from the SEM micrograph (W JJ ). We only observe a constant distribution of critical current throughout the full width of the junction (expected for uniform S−N−S junctions) when we incorporate the calculations presented in this paper. Contrarily, the analysis based on the L eff = 2λ + d yields an unphysical concentration of critical current in the middle of the junction. Finally, note that the current is confined to the actual junction (W JJ ), not the full width of the superconducting film (W ). This explains the constant offset in Figure  6a.
Alternatively, we can use a linear approximation of γ d 2 , y to mitigate the need for rescaling the axes. Figure 7c shows the same Fourier analysis based on a linear approximation of f y W . Since the linear approximation of f disk y W breaks down near the edges, it yields less precise results at the junction boundaries. However, in the middle of the junction, the linear approximation of f y W is well suited for correctly analysing J c (y). For the technical details of carrying out the Fourier transform, the reader is referred to the Supplemental Material [29].

V. CONCLUSION
In conclusion, we analyzed the periodicity ∆B of the interference pattern I c (B) for thin film planar S−N−S Josephson junctions, both theoretically and experimentally. Specifically, we examine junctions separating rectangular, ellipsoid, and rhomboid films of width W and length L. By mapping the second Ginzburg-Landau equation to the two-dimensional Laplace equation, we solve I c (B) for different ratios of L/W . We show that ∆B has two universal limits for L W and L W , independent of the sample geometry. The first corresponds to an infinite superconducting strip, and the latter is caused by an emerging universal dependence of the phase difference on the junction electrode surface area. By fabricating elliptically-shaped S−N−S junctions, having different ratios for L/W , we experimentally verify the geometry dependence of ∆B. Lastly, we adapt the Fourier relation between I c (B) and the critical current density distribution to suit planar junctions in the thin film limit. This proves to be vital in correctly predicting the location of current channels in topological planar Josephson junctions.

VI. ACKNOWLEDGEMENTS
This work was supported by the Dutch Research Council (NWO) as part of the Frontiers of Nanoscience (NanoFront) program and through NWO projectruimte Grant No. 680.91.128. The work was also supported by EU Cost Action CA16218 (NANOCOHYBRI) and benefited from access to the Netherlands Centre for Electron Nanoscopy (NeCEN) at Leiden University.

Appendix: Derivation of the boundary conditions
As discussed in the main text, the Neumann boundary conditions for the gauge-invariant phase are given by: In this appendix, we will derive the results presented in Table I.
Combining the choice of the gauge A A A = −yBx witĥ n n n R =x for boundary 1, we find: For boundary 2 we obtain the same result, yet with a minus sign sincen n n R = −x. For boundary 3n n n R = ±ŷ, which yields (∇ ∇ ∇γ) ·n n n R ∼x ·ŷ = 0. Next, for boundary 4, parameterize the ellipse as L 2 cos tx + W 2 sin tŷ. The tangent is then given by the derivative to t, which is − L 2 sin tx + W 2 cos tŷ = − Ly Wx W x Lŷ . Here, in the second step, we transformed back to Cartesian coordinates lying on the ellipse. A vector perpendicular to the tangent, pointing inwards to the ellipse, is then given by: − W x Lx − Ly Wŷ . Normalizing yieldsn n n R : Taking the inner product with A A A, as in Eq. A.2, yields the boundary condition in Table I. Finally, for boundary 5, define the angle α as arctan(W/L). In that case, for y > 0, we findn n n R = − sin αx + − cos αŷ, such that Again yielding the boundary condition in table I. Note that the boundary condition is unchanged for y < 0, even though the y-component ofn n n R acquires a minus sign. This results from the choice of gauge (A y = 0).

S1. DATA OBTAINED ON A BAR-SHAPED SAMPLE
In the main text we present results obtained on elliptical samples, in this Supplemental Material we show the results obtained on a rectangular junction. In Figure S1 a scanning electron micrograph of this junctions accompanies a color plot of the I c (B)-pattern obtained on this sample. We have used this sample to examine the geometrical independent limit where L W . Since the oscillation amplitude of the bar-shaped sample decreases more strongly with increasing B than the ellipse-shaped samples, we cannot determine ∆B using the fifth side lobe of the pattern. Instead, we use the width of the fourth side lobe to establish ∆B. Following the discussion in the main text this yields: ∆BW 2 JJ /Φ 0 = 1.70. Interestingly, if we calculate ∆B by taking the average of the first four side lobes in the interference pattern, we find ∆BW 2 JJ /Φ 0 = 1.83.
Supporting Figure S1. The Fourier transform is defined in the main text as: In the right hand side, the transform I c is complex valued; its real and imaginary parts encode for the even and odd components in I c (β) respectively. Since the experimental interference patterns are mainly symmetric (i.e., an even function of the applied magnetic field), we can assume I c to be dominantly real: The real component of I c is an oscillating function that flips sign at each zero crossing. The imaginary part is significantly smaller than the real part, except at the zero-crossing where the even part vanishes. Therefore, the imaginary part of I c (I c,odd (β)) can be approximated by the critical current at the minima in the experimental interference pattern. Also I c,odd (β) is flipping its sign between each minimum and between the minima we approximate I c,odd (β) by linear interpolation. The inverse transform yieldingJ c (ỹ) from I c (β) is then given by: Figure S2 gives an overview of the subsequent steps of the Fourier analysis. First, I c is extracted from the experimental data by a defining a voltage threshold; this is depicted in Figure S2a for a disk-shaped junction (data in Figure 5b of the main text). We vertically translate the extracted I c values such that the global minimum equals zero current. This step in the data analysis prevents the overestimation of I c,odd (β), which would result in an overly anti-symmetric current density distribution. I c,even (β) is found by multiplying the translated I c by a flipping function that changes the sign of each subsequent lobe of the interference pattern, as can be observed in Figure S2b.
Supporting Figure S2. Overview of the subsequent steps of the Fourier transform analysis. Here we also rescale the field axis toβ. We follow the above procedure for finding I c,odd (β), which is depicted in Figure S2c. The corresponding critical current density distribution is found by a numerical Fourier transform carried out in Python using the Numpy package, yielding the distributionJ c (ỹ), depicted in Figure S2d. Finally both axes are rescaled using f disk y W to retrieve J c (y), which is shown in Figure 7a of the main text. For illustrative purposes, we have chosen to absorb part of the prefactor inỹ, as this this yields a larger contrast between Figure S2d and Figure 7a of the main text. However, as discussed in the main text, any choice ofỹ andβ is allowed, as long as it is consistent with γ. We indicate −W/2 and W/2 by solid reference lines and the boundaries of the actual weak link by dotted reference lines. For the Fourier transform using a linear approximation of f disk y W we definẽ y = y and absorb the fit of f disk y W intoβ.