Giant fractional Shapiro steps in anisotropic Josephson junction arrays

Giant fractional Shapiro steps have been observed in Josephson junction arrays as resulting from magnetic flux quantization in the two-dimensional array. We demonstrate experimentally the appearance of giant fractional Shapiro steps in anisotropic Josephson junction arrays as unambiguous evidence of a skewed current phase relationship. Introducing anisotropy in the array results in a giant collective high frequency response that reflects the properties of a single junction, as evidenced by the observation of a Fraunhofer like magnetic field dependence of the total critical current of the system. The observed phase dynamics can be perfectly captured within an extended resistively shunted Josephson junction model. These results directly indicate the potential of Josephson junction arrays to explore the current phase relation in a very broad frequency range (down to 50 MHz) and in a wide variety of novel link materials exhibiting non-conventional current phase relationships. Giant fractional Shapiro steps have been observed in Josephson junction arrays as resulting from magnetic flux quantization in the two-dimensional array. Here, the authors demonstrate the observation of giant fractional Shapiro steps in an anisotropic Josephson junction array, as unambiguous evidence of a skewed current phase relationship.

T he current-phase relationship (CΦR), describes the relation between the supercurrent, I s , through a Josephson junction (JJ) and the gauge invariant phase difference across the junction, Δθ. The exact relation is determined by the details of the Andreev bound state (ABS) spectrum that exists in the junction [1][2][3] . For a conventional superconductor-insulator-superconductor (SIS) junction, the current-phase relation is sinusoidal, I s = I c sin (Δθ), where I c is the junction critical current. With the advances in materials sciences and the realization of novel link materials, exhibiting non-conventional CΦRs, JJ with new functionalities arise having a tremendous potential for applications [4][5][6][7][8][9] .
To investigate the harmonic content of the current-phase relationship in these novel weak link materials, experiments relying on detecting the phase-locked response of a current driven Josephson junction when subjected to a high frequency radiation field of frequency, ν rf , have recently been used [10][11][12][13][14][15] . In case the weak link is characterized by a conventional sinusoidal current-phase relationship, this resonant response manifests itself as constant voltage plateaus, so called Shapiro steps, in the voltage versus current (VI)-characteristics of the junction at voltages V n = nΦ 0 ν rf , where n 2 Z and Φ 0 the fluxquantum 16 . However, for a weak link having a 2π periodic non-sinusoidal CΦR, the response exhibits, in addition to the conventional integer Shapiro steps, steps at fractional values V n/q = (n/q)Φ 0 ν rf , where the qth (q ≠ n) fractional step originates from the phase-locked response with the qth harmonic of the CΦR [17][18][19] . Similarly, the 4π periodicity of the Majorana bound state spectrum in topological Josephson junctions, is expected to leave a fingerprint in the Shapiro step response 10 .
Despite the experimental accessibility of these techniques, the detailed interpretation of such experiments is complicated by the non-linear nature of the Josephson response 10,15,20,21 . Moreover, the use of a single junction makes these experiments challenging mainly due to its low response (or radiation power) and small impedance 12,22 . Josephson junction arrays (JJAs) containing many junctions provide the natural alternative to provide the necessary enhancement of the coherent response and the design degree of freedom allows to tune their impedance for optimum integration in the measurement/excitation circuit [23][24][25] . Although JJAs hold a strong promise to study the CΦR, the array geometry has been predominantly studied as a model system for the Kosterlitz-Thouless (KT) transition and the fully frustrated XYmodel [26][27][28][29][30][31] . Whereas, the impact of the particular CΦR on their high frequency response remains until now unexplored.
In this paper, we first explore the dependence of the giant Shapiro response on the anisotropy of the array. We demonstrate that for highly anisotropic arrays, the giant collective high frequency response reflects the properties of a single junction. Further, we demonstrate experimentally the appearance of giant fractional Shapiro steps in anisotropic Josephson junction arrays as unambiguous evidence of a non-sinusoidal CΦR.

Results
We present transport measurements on 2D proximity-induced superconductor-normal metal-superconductor (SNS) Josephson junction arrays having very different lattice parameters (Fig. 1a). The JJAs consist of a rectangular lattice of square 25 nm-thick molybdenum germanium (Mo 79 Ge 21 ) superconducting islands, on top of a metallic gold (Au) film that couples the superconducting islands through the proximity effect 32 . We define the Fig. 1 Sample layout and DC transport measurements. a A schematic illustration, together with a false-coloured scanning electron microscopy image (of Device 1), of the superconductor-normal metal (S-N) hybrid devices studied in this work. The black scale bar corresponds with 200 nm. b The temperature, T, and DC current, I T , dependence of the Device 2 differential resistance, dV T /dI T , under a zero applied magnetic field condition. NS, VS and Z-VS indicate the normal state, voltage state and zero-voltage state of the Josephson junction array, respectively. Concerning the discontinuity of the VS-NS transitions around T = 1. 8  anisotropy, η = d y /d x , of the rectangular array as the ratio of the spacing between islands perpendicular, d y , and along the current direction, d x . For both devices the spacing along the current direction is quasi identical, d x = 250 nm, while the spacing perpendicular to the current direction, d y , varies. Device 1 has an anisotropy of η = 0.8. Whereas, Device 2 has a large anisotropy of η = 4. Further, Device 1 has N x = 200 columns of N y = 13 square islands, while for Device 2, N x × N y = 160 × 6. We define the xdirection to be colinear with the current bias and the magnetic field is always applied along the z-direction. Due to the similarities between the physics of JJAs and the XY-model, we express the applied magnetic field in reduced units, f = H/H 1 , known in the JJA literature as frustration. Where H 1 = Φ 0 /A, with A the area of the unit cell of the JJA. For Device 1, H 1 = 3.7 mT, while for Device 2 having the larger unit cell, H 1 = 1.8 mT.
Both devices show global superconductivity having a critical temperature of T c = 4.0 K and T c = 4.15 K, for Device 1 and 2, respectively 33,34 . To obtain the characteristics of an individual junction out of measurements reflecting the collective response of the whole array, we measured the VI-characteristics in a highly symmetric situation (zero frustration) where each junction probes the same current. This condition is achieved under zero-field, where the applied current, I T , is distributed equally between the rows, in a way that each junction parallel to the current is crossed by a current I = I T /N y 35 . In such state, the junctions perpendicular to the current are never activated, while the voltage across each parallel junction is V = V T /(N x − 1), with V T the corresponding total value measured along the array. The temperature dependence of the differential resistance, dV T /dI T (T), of the arrays under zero external magnetic field while sweeping the bias current for temperatures ranging from 0.3 to 5 K is shown in Fig. 1b (for Device 2).
It is clear from Fig. 1b that for currents below the critical current, I c (T), the JJA is in the zero-voltage state, where the current is a pure supercurrent. When the current bias exceeds the critical current, the extra current is carried by a resistive channel of quasiparticle excitations and the JJA resides in the voltage state, which exhibits a linear current-voltage relation (see Fig. 1c, d).
Finally, when the current bias reaches the de-pairing current, the sample will transit to the normal state. A clear fingerprint of the nature of the weak link is given by the temperature dependence of I c (T). As shown in Fig. 1e, I c (T) exhibits an exponential decay with T at high temperatures and saturates at lower temperatures, consistent with Likharev's conventional theory of long diffusive SNS junctions (see Supplementary Information Note 1.1) 34 . In this case, the CΦR can be written as: where I c = max(|I s |) is the junction's critical current and the temperature dependent parameter, ' T ð Þ ¼ I c dI s dΔθ À1 π þ 1, has a simple geometrical interpretation, as it determines the level of skewness of the CΦR 2,35 Whereas, at high temperatures the current-phase relationship is harmonic ' ¼ 0, at low temperatures the current-phase relationship becomes considerably forward-skewed 0 < ' < 1. As a result, the array of long diffusive SNS-junctions (d x /ξ N (T c )~5, with ξ N the normal metal coherence length, (see Supplementary Information Note 1.1)) provides a prototypical system to investigate the effect of the harmonic content in the CΦR on the high frequency response, as it can be easily tuned by temperature 3,36 .
Giant Shapiro response at integer frustration. We initially investigate the high frequency response at high temperatures, where the current-phase relation of an SNS junction is sinusoidal, ' ¼ 0.
Due to the smallness of the fluxquantum, Φ À1 0 ¼ 483:6 MHz=μV, typically microwave excitation frequencies (ν rf > 1 GHz) are necessary in order to observe measurable Shapiro steps, V n . By constructing a (N x − 1) × (N y − 1) Josephson junction network the Shapiro step response is multiplied by a factor N x − 1, resulting in a giant Shapiro response at V T n ¼ ðN x À 1ÞV n and becomes, for the array size of the devices considered in this work, measurable at radio frequencies (RF) as low as 50 MHz. The scaling in the transverse direction, N y , increases the overall critical current of the device.
A typical Shapiro response for Device 1 is shown in Fig. 2a at f = −1, which similarly to f = 0 results in a highly symmetric (nonfrustrated) situation where the whole JJA behaves as a single JJ. Both the time averaged VI-scan at T = 3 K and the differential resistance V T /dI T , are shown using a current bias, I T = I dc + I rf (t). Where, I dc is a DC current and I rf (t) = I rf cos(2πν rf t) is an RF excitation of amplitude I rf and frequency, ν rf = 200 MHz. Note, that since the RF excitation is experimentally applied by a function generator, the conversion between the applied AC voltage, V rf , and the transmitted AC current, I rf , is a priori unknown (see Supplementary Information Note 2). The pronounced Shapiro response at the voltage values V T n ¼ nðN x À 1ÞΦ 0 ν rf , with N x = 200, indicates that the whole Josephson junction network has a coherent response where Φ 0 is the fluxquantum, ν rf the RF frequency of the drive and N x = 200. b The experimentally measured frequency dependence of the Shapiro steps for a fixed RF excitation amplitude of V rf = 100 mV for Device 1. Shapiro steps are clearly reflected in the map, as the zones of low differential resistance (yellow). The black stripe lines indicate the expected theoretical nth Shapiro step response according to under RF excitation, in which all N x − 1 = 199 junctions contribute, despite the presence of an unavoidable disorder in the system due to sample fabrication 37 . Figure 2b, shows the experimentally measured frequency dependence for a fixed RF excitation amplitude for Device 1. Shapiro steps are clearly reflected in the map, as the zones of low differential resistance (yellow). The black stripe lines indicate the expected theoretical nth Shapiro step response. The linear evolution of integer Shapiro steps down to 50 MHz, directly show the gained resolution of the JJA compared to a single junction and its potential to explore the properties of a single junction using the collective response of the JJA 12 .
Recent works, exploring topological weak links, have shown that a Shapiro map is a powerful representation of the response at a fixed RF frequency to obtain information regarding the weak link CΦR 10,11,38 . A Shapiro map traces the widths of the observed plateaus as a function of the AC drive amplitude and DC current bias. For a current biased system, the width of the Shapiro current plateaus can deviate from the ideal Bessel dependence, expected for a JJ which is voltage biased. For the current biased case the current width of the nth Shapiro step depends in a non-trivial way on the reduced frequency 21,[37][38][39][40] : where R is the resistance of the Junction. For an harmonic CΦR, it has been shown in refs. 20,21 that the current step widths of the integer Shapiro steps follow qualitatively a Bessel like dependence ΔI n~Ic J n (I rf /Ω), with J n the nth order Bessel function of the first kind. Only in the case of low AC current amplitudes, low reduced frequencies and low Shapiro step numbers, the step widths will deviate qualitatively from the Bessel like dependence. Figure 3a shows the experimentally measured Shapiro map, as a function of the AC drive amplitude and DC current for a fixed radio frequency, ν rf = 150 MHz (corresponding to Ω~0.12), for Device 1. At the lowest RF powers, the influence of the irradiation is small and only the n = 0 plateau, corresponding with the zerovoltage state is observed. Upon increasing the RF power, the map starts to show higher order Shapiro plateaus. Only integer Shapiro steps are present in the map. It is clear that depending on the RF amplitude a different Shapiro step response will be observed. For example, while for V rf = 240 mV, we will encounter dominantly odd Shapiro steps. For V rf = 300 mV we will encounter dominantly even Shapiro steps in the VI-characteristic. Similar observations can be made when varying the RF frequency or temperature (See Fig. 1b and Supplementary Information Note 3). This even-odd effect 41 , results in the conspicuous absence of steps at particular RF amplitudes or frequencies and is attributed to the oscillatory dependence of the width of the Shapiro current plateaus on the reduced frequency 21 . Figure 3b shows the Shapiro response obtained within the RSJ model (see Methods section) for a single junction having a sinusoidal CΦR. The simulated response shows excellent quantitative agreement with experiment, apart from a small error (clearly observed at I dc = 0 μA) in the oscillation period of the Shapiro step widths with the AC drive amplitude. This error results from a small mismatch between the I rf current value used in the extended RSJ model and the V rf value, which is experimentally applied. This indicates that our JJA can be well described by the current driven RSJ model for a single junction. Figure 3c shows the evolution of the first n = 0, 1, 2, 3 Shapiro current step widths with applied RF amplitude, extracted from the experimental data (dots) in Fig. 3a. Although, the behavior of the step widths for the current biased case are known to deviate from the ideal Bessel dependence, for the high value of the reduced frequency, Ω = 0.12, a qualitatively similar behavior is found as indicated by the fit (solid lines) 21,39 .
Giant Shapiro response at non-zero frustrations. The applied magnetic field provides an interesting experimental knob to modify the response of the system as it has a direct impact on the critical current and, as a result, on the reduced frequency. In addition, a magnetic field induces frustration in the JJA and a detailed investigation is required to differentiate both effects. Figure 4a shows the measured and simulated contour map of the differential resistance as a function of frustration and reduced DC voltage, V T =V T 1 , for T = 3 K. For f ∈ [−2, 2] strong Shapiro features can be observed at integer values of the frustration, whereas shallow features can also be observed near half-integer values of the frustration (see Supplementary Information Note 4 and 5). The observation of giant Shapiro steps at fractional or integer multiples of the frustration is due to coherent dynamics of the vortex lattice under RF excitation when commensurability exists between the magnetic field induced vortex lattice and the underlying Josephson junction lattice [26][27][28][29][30] . Deviating from these high symmetry conditions (i.e. inducing frustration), results in a suppression of the coherent response of the JJA. The high symmetry conditions at particular filling factors are independently In particular the experimentally obtained dependence of the differential resistance, dV T /dI T , as a function of the AC drive amplitude, V rf , and DC current, I T dc , for a fixed radio frequency, ν rf = 150 MHz, for Device 1 at T = 3.3 K is shown in a. The Shapiro map generated from the extended resistively shunted Josephson junction (RSJ) model (see Methods section) for a single junction is shown in b. In the latter case, the differential resistance is multiplied by a factor 199/13, the current by a factor 13 and the voltage with a factor 199, in order to obtain the corresponding cumulative array response. c The width of the first four Shapiro steps, ΔI n for n = 0, 1, 2, 3, are plotted in function of the RF amplitude, I T rf . The lines show a fit using a Bessel dependence with parameters, A = 30.12 μA and β = 0.0767 μA −1 .
confirmed trough increased critical current values at the same magnetic field values for T = 2 K (Fig. 4b).
The high anisotropy of the array in Device 2 (and consequently the critical current, see Supplementary Information Note 1.4), is expected to result in a strong suppression of effects related to flux quantization in the array. Consequently, the anisotropic JJA is expected to reflect single JJ behaviour, while simultaneously keeping the advantages of high device critical currents and giant Shapiro responses, the array geometry provides. However, disorder in the array, inherent to the fabrication process, might induce broadening of the Shapiro steps 42 . Figure 4c shows the measured and simulated contour map of the differential resistance as a function of the applied magnetic field and reduced DC voltage, V T =V T 1 , for Device 2 at T = 1.8 K. The observed magnetic field dependence can be explained by the variation of the critical current with field at T = 2 K, shown in Fig. 4d, exhibiting the typical Fraunhofer-like dependence, illustrative for a single uniform weak link. We can estimate the junction area by fitting the I c (B, T) with the Fraunhofer function (see Supplementary Information Note 1.3). This provides us with an estimate of the effective area of the parallel junctions to be 0.173 μm 2 , which corresponds well with the geometric area of the parallel junction, 0.125 μm 2 . This single junction behaviour of the highly anisotropic JJA, Device 2, is further evidenced by performing simulations within the extended RSJ model on a single junction, showing excellent agreement with experimental observations (see Supplementary Information Note 6).
Giant fractional Shapiro response. The enhanced Shapiro response in anisotropic JJAs, free from spurious effects due to flux quantization, provides an excellent playground to investigate the CΦR by probing their RF response and potentially their enhanced RF emission spectra 12 . As discussed, at low temperatures the CΦR of the SNS junctions is expected to become considerably skewed. In this case it is expected that fractional Shapiro steps can appear due to phase locking between the incident RF radiation and a higher harmonic of the CΦR [17][18][19] . Figure 5a (negative field values), shows the experimentally obtained field dependence of the Shapiro response for the anisotropic JJA (Device 2) at T = 0.4 K. Some clear differences can be observed with the map at T = 1.8 K (Fig. 4c). The most pronounced difference is the observation of giant half-integer Shapiro steps at larger values of the frustration (see Fig. 5b). As the anisotropic array exhibits single junction behaviour, the half-integer steps have to originate from an intrinsic property of the single junctions within the array. To this end, we show in Fig. 5a (positive field values) the simulated field dependence for a 6 × 6 array using the anisotropic array junction parameters (for simulations on a single junction, see Supplementary Information Note 7), including the possibility of a single-valued, but skewed Fig. 4 Field dependence of the Shapiro response and field dependence of the critical current. a Field dependence of the Shapiro response for Device 1, presented as a contour map of the differential resistance, dV T /dI T , versus frustration, f, and the DC voltage, V T , at T = 3 K. The DC voltage is scaled by where Φ 0 is the fluxquantum and ν rf the radio frequency (RF) of the drive. The experimental measurement using an RF excitation of ν rf = 150 MHz and V rf = 100 mV, is shown for negative frustrations. The result obtained within the extended resistively shunted Josephson junction (RSJ)-model for a 13 × 13 array, using an RF excitation, I rf = 7.69 μA applied to each junction is shown for positive frustrations. The differential resistance obtained from the simulation for a 13 × 13 array is multiplied by a factor 13(N x − 1)/(12N y ), in order to match the value for the array size, (N x − 1) × N y of the experiment, with (N x , N y ) = (200, 13). b Frustration dependence of the voltage, V T , versus current, I T , characteristics of Device 1 at T = 2 K when no RF excitation is applied, the black contour line indicates the evolution of the critical current as determined by using a voltage criterion of 0.1 μV. c Field dependence of the Shapiro response for Device 2. The experimentally measured contour map of the differential resistance versus frustration, f and the DC voltage, V T , for T = 1.8 K is shown for negative frustrations. The RF excitation parameters are ν rf = 175 MHz and V rf = 800 mV. The DC voltage is scaled by V T 1 ¼ ðN x À 1ÞΦ 0 ν rf ¼ 57:5 μV. The result obtained within the extended RSJ-model for a single junction, for an RF current bias of I rf = 200 μA applied to the junction, is shown for positive frustrations. The differential resistance obtained from the simulation for a single junction is multiplied by a factor (N x − 1)/N y , in order to match the value for the array size, (N x − 1) × N y of the experiment, with (N x , N y ) = (160, 6). At low fields f = [−6, 0] the experimental data shows qualitatively similar modulations with the magnetic field, f, which cannot be resolved with the current colour scale. d Field dependence of the differential resistance of Device 2 at T = 2 K when no RF excitation is applied, the black contour line indicates the evolution of the critical current using a resistance criterion of 5 Ohm.
Further, Fig. 5a shows the evolution of dominantly even steps at low frustrations, |f| < 2, to dominantly odd steps at frustrations 2 < |f| < 4, and eventually to a combination of half-integer and integer steps at |f| > 4. The simulated maps feature a similar field dependence. This intriguing field dependence can be directly related to the evolution of the Shapiro step widths with the reduced frequency, Ω. For Device 2, the reduced frequency changes considerably from Ω = 0.0063 to Ω = 0.0307 when increasing the magnetic field from f = 0 to f = 10 due to the Fraunhofer-like field-dependence of the critical current.  (Fig. 5a) and indeed confirms that at f = 0 (Ω = 0.0063) no fractional Shapiro steps can be observed, while at f = 10 (Ω = 0.0307) half-integer steps show a pronounced fingerprint. Further, a myriad of fractional resonances, corresponding to phase locking to different harmonics of the CΦR, can be found, as shown in the inset. These, however, are considerably narrower than the half-integer ones and thereby more susceptible to being washed out by thermal fluctuations. The complete evolution of the Shapiro step widths with reduced frequency is shown in Fig. 6b as obtained for a reduced RF amplitude of I rf /I c = 6.04, corresponding to the f = 0 case. It is clear that for this reduced RF amplitude fractional steps become appreciable only above Ω = 0.01 (see Supplementary Note 8).

Discussion
In this work, we demonstrated that for a prototypical SNS-based JJA, the presence of fractional steps at the lowest explored temperature of T = 0.4 K, demonstrates unambiguously that the CΦR is considerably skewed, with ' ¼ 0:5. This value corresponds well with the value expected within the equilibrium quasi classical theory 3,35 . Note however, that non-equilibrium effects, inherent to the measurement technique, can modify the occupation of the ABS that persist up to high temperatures 14,36,[43][44][45][46][47] . The giant Shapiro response in a JJA, in contrast to a single junction, enables to measure the RF response using sub gigahertz frequencies. It allows, similarly as in Dassonneville et al. 46 , to explore the impact of these nonequilibrium processes on the ABS over a wide range of frequencies.
Further, we have shown that a complete picture of the response is necessary in order to obtain information about the CΦR from the complex non-linear response.

Methods
Device fabrication. The devices shown in Fig. 1a consist of a rectangular lattice of square superconducting islands, which are deposited on top of a metallic gold (Au) film, patterned for four-point transport measurements using standard e-beam lithography. The superconducting square islands are 25 nm-thick molybdenum germanium (Mo 79 Ge 21 ) thin films of size a = 500 nm deposited using pulsed layer deposition. To prevent oxidation they are in-situ covered by a 5 nm-thick Au capping layer. The Au film underneath the array is deposited on a Si/SiO 2 substrate using molecular beam epitaxy. A 5 nm-thick titanium sticking layer is used to assure good adhesion of the Au film onto the substrate. Device 1 has a 40 nm-thick Au base layer. Whereas, Device 2 has a 30 nm-thick Au base layer.
Electrical characterization. The devices are wired to a chip carrier that can be placed in two variable temperature cryostats. Measurements for temperatures in the range 0.3-1.6 K are performed in a He 3 system, while the high temperature range is explored in a He 4 flow cryostat. Due to the reduced cooling power of the He 3 system, the SN-transition in Fig. 1b shifts towards lower currents. For the fourterminal transport measurements, characterizing the sample (Fig. 1b), we measure the differential resistance using a standard lock-in technique.
where Φ 0 is the fluxquantum and ν rf the radio frequency of the drive. The result obtained using the extended resistively shunted Josephson junction (RSJ)model (see Methods section) for a 6 × 6 array, including the field dependence of the critical current, a skewed current-phase relationship (CΦR) with ' ¼ 0:5 (see Methods section) and using an RF excitation, I rf = 300 μA applied to each junction, is shown for positive frustrations. The value for the differential resistance is multiplied by a factor 6(N x − 1)/(5N y ), in order to match the value for the array size, (N x − 1) × N y of the experiment, with (N x , N y ) = (160, 6). b The average differential resistance, dV T dI T D E 7<f<9 , in the range 7 < f < 9 versus the reduced DC voltage, indicating clearly the presence of fractional Shapiro steps.
The Shapiro response presented in Figs. 4a, c and 5a are measured by simultaneously biasing the JJA with a DC current bias (Keithley 6221) and an RF AC voltage drive (Rhode&Schwarz SMY signal generator), while measuring the differential resistance using a standard lock-in technique. The corresponding voltage value is obtained by integrating numerically the differential resistance. The Shapiro responses presented in Figs. 2b and 3a are not performed with a lock-in amplifier in order to obtain directly dV T /dI T . Instead we measured V T I T -curves and calculated the derivative, dV T /dI T , using a 9-point central difference stencil method.
Extended RSJ-model. To understand our experimental findings, we model the square lattice of N x × N y superconducting islands on a thin metallic layer as an array comprising overdamped Josephson junctions, each shunted by a resistive load 26,34,39,48,49 . The junctions forming the array correspond to the SNS-junctions formed between nearest neighbouring islands as illustrated in Fig. 1a. Within this extended resistively shunted JJ model (RSJ-model), the current I ij from the ith to the jth island is given by: The first term on the right hand side in Eq. 3 describes the supercurrent, I s;ij , through the junction which is a function of the gauge invariant phase difference over the junction, θ ij − a ij . Where, θ ij = θ j − θ i and θ i is the phase on island i. Further, a ij ¼ 2π A Á dl and A is the vector potential. The relation between the supercurrent and the gauge invariant phase difference is determined by the nature and shape of the weak link and can be described, for the majority of junctions, roughly by two constants, namely, the critical current I c;ij and ' ij 35 .
Note that ' ij has a simple geometrical interpretation and determines the skewness of the CΦR. For a sinusoidal CΦR, ' ij ¼ 0; In case the curve is forwardskewed and single valued as expected for high transmission SNS-junctions, 0 < ' ij < 1. The second term on the right hand side in Eq.3 is the current through a resistive shunt having resistance R ij and it stems from a dissipative quasiparticle current. Finally, the third term on the right hand side in Eq.3 is a current contribution due to thermal fluctuations, taken into account as a Langevin noise in the junction, with covariance hγ t ð Þγ t 0 ð Þi ¼ 2k B T R δðt À t 0 Þ.
Since in our experiment, the effective penetration depth Λ, exceeds the width of the transport bridge in the whole temperature range, Λ > N y (a + d y ), the self-field induced by the current through the junctions can be neglected and the vector potential can be approximated by A ¼ Bxŷ, where B is the magnitude of the applied magnetic field 50 . The equations of motion for each junction are coupled to each other by the requirement that the current coming in-and out each superconducting island is conserved, that is, for each SC island i, X where j is restricted to the neighboring cells and I i is the current injected directly into the island due to an external current bias. Meaning I i = I for the islands at the left boundary, I i = −I for those at the right boundary, and 0 otherwise. As the impedance of the junctions are much smaller than the impedance of the external drive and system, describing the drive as a current biased system is justified 22 . The system of coupled stochastic equations is integrated following a predictor-corrector scheme with time step Δt = 0.01τ J , where τ J = ħ/2eRI c is the phase particle relaxation time, which sets the typical time scale for the phase to adapt to a change of the drive current. The junctions parameters R ij and the critical current I c;ij (B, T), can be determined from experiments (see Supplementary Information Note 1.2). The field and temperature dependence of I c;ij is explicitly included in the simulations. More details about the numerical method can be found in the Supplementary Information Note 8.

Data availability
The data that support the findings of this work are available from the corresponding author upon reasonable request.