Large Second-Order Josephson Effect in Planar Superconductor-Semiconductor Junctions

Growing interest in the integration of new materials into quantum devices offers opportunities to study proximity effects, such as between superconductors and semiconductors. Elements as basic as planar junctions, with two superconductors placed side-by-side, can host Majorana zero modes, and be used as nonlinear quantum circuit elements [1, 2, 3]. Depending on the junction material, planar junctions can also be used to search for exotic phenomena such as triplet superconductivity [4].

The primary characteristic of a Josephson junction (JJ) is the current phase relation (CPR) [5, 6]. A CPR connects the supercurrent to the phase difference across the junction and is calculated from the weak link energy spectrum and interface transparency. Among other properties, it predicts the electromagnetic response of a Josephson junction in a circuit.

The CPR of a superconductor-insulator-superconductor (SIS) Josephson junction, as originally derived, is sinusoidal, $I({\varphi})=I_{c}\sin\varphi$, where $I$ is the supercurrent, $\varphi$ is the phase difference between two superconducting leads, and $I_{c}$ is the critical current. More generally, a CPR can be decomposed into a Fourier series, $I({\varphi})=\sum_{n}I_{n}\sin(n\varphi)$ [6]. In real SIS devices, higher-order components may be present but the amplitude would be small [7]. In superconductor-normal metal-superconductor (SNS) JJs, the CPR can deviate significantly from the sinusoidal form due to contributions from Andreev bound states which are subgap quasiparticle states.

In the clean SNS limit, the CPR is predicted to be linear or a skewed sine function [8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. Even in junctions with skewed CPR, the first harmonic, $\sin(\varphi)$ is typically dominating. $4\pi$-periodic, or the half-integer harmonic in CPR is predicted in topological superconductors and explored in a variety of materials [18, 19].

Large second-harmonic CPR was searched for in a variety of junctions at the so-called 0-$\pi$ transition [20, 21, 22, 23, 24, 25] or in 45^∘-twisted high-temperature superconductor junctions [26, 27, 28, 29, 30]. In these cases the second harmonic can be observed because the first harmonic is cancelled.

In SNS junctions all harmonics are present, including the second harmonic. However because of the presence of even higher-order terms, we do not expect double modulation. If a SQUID consisting of such junctions is flux-biased to $\pi$, the first harmonic can be canceled, resulting in the domination of the second harmonic signatures  [31, 32, 33, 34, 35]. This effect is used to create the $\pi$-periodic energy-phase relation for so-called “$0-\pi$ qubit” [36, 37, 33, 38].

A possible experimental signature of the second-order Josephson effect is half-integer Shapiro steps [39], which can also have a variety of other origins due to phase locking of Josephson vortices or quasiparticle dynamics [40, 41]. Another type of evidence is in supercurrent interference patterns where additional minima or kinks are observed at values of half-integer magnetic flux quanta [25, 20].

In our experiments, we find Josephson junctions with an unusually large second harmonic in planar junctions based on an InAs quantum well. We find confirmation of this in several measurements. First, we observe that the superconducting quantum interference devices (SQUIDs) made of two planar InAs junctions exhibit extra kinks in the flux modulation. The shape of the SQUID characteristics evolves with junction asymmetry, which can be explained by a simple two-component CPR model. Second, single planar junction supercurrent diffraction patterns show kinks near half-quantum of applied flux. Finally, half-integer Shapiro steps are also observed. In our single junction measurements, the observation of the large $\sin(2\varphi)$ term is not related to a cancellation of the first-order term, e.g. at the 0-$\pi$ transition.

The amplitude of the second harmonic from several of these measurements is found to be around 0.4 of the first harmonic. The sign of the second harmonic is determined to be negative from the model that includes the loop inductance (see supplementary information). Negative sign is expected upon decomposition of a skewed sinusoidal function into Fourier components. Because of the long mean free path on the order of the junction length and high interface transparency, it is not surprising if the CPR has higher-order terms. However, a skewed function should contain all sinusoidal harmonics, including those higher than two. We do not observe an apparent third-order or higher-order terms, which is surprising and requires further studies.

Planar InAs quantum well junctions are prepared using the nanowire shadow method [42, 43]. Because the junction area is not subjected to chemical or mechanical etching, this approach preserves the quantum well in the junction area. Superconductors used are Al or Sn. We focus on Al devices in the main text. Data from Sn devices are available in supplementary materials [44]. Standard electron beam lithography and wet etching is used to pattern the mesa outside the junction. Measurements are performed in a dilution refrigerator at 50 mK unless otherwise stated.

We first present examples of dc-SQUID patterns with extra features that we study in the context of the $\sin(2\varphi)$ terms. Device SQUID-1 is depicted schematically in Fig. 1(a) and an optical microscope image is shown in Fig. 1(b). The device has two Josephson junctions in parallel, JJ_a and JJ_b. Two nanowire top-gates $V_{g,a}$ and $V_{g,b}$ are used for tuning critical currents. We set supercurrent to zero in one junction with its gate and measure the gate dependence of the switching current $I_{sw}$ in the other junction [Fig. 1(c)]. $I_{sw}$ is typically referred to as critical current, though the true Josephson critical current can be higher than the switching current. $I_{sw}$ quenches at gate voltages near $-100$ mV and saturates above 300 mV in both junctions. Both junctions have the same nominal geometry, so their magnetic flux modulation, or diffraction, patterns are expected to have similar periods. This can be confirmed in Fig. 1(d), where there is only one Fraunhofer-like low-frequency modulation while both junctions are in the superconducting state. With the Fraunhofer period of 0.29 mT and the junction width of 5 $\mu$m, we get an effective junction length of 1.4 $\mu$m (an order of magnitude larger than the typical physical length) for both JJs. The long effective length is likely due to large London penetration depth which is typical for thin film superconductors. We use the junction magnetic flux $\Phi_{JJ}$, given in the unit of $\Phi_{0}=h/2e$. The SQUID flux and the junction flux are applied from a large superconducting magnet and cannot be independently controlled. The high frequency oscillations within the Fraunhofer-like envelope are due to interference between the two junctions. The period 1.57 $\mu$T gives a SQUID area of $1.31\times 10^{3}$ $\mu$m² which is similar to the area of the inner loop ($1.28\times 10^{3}$ $\mu$m²). The oscillations are shown in detail in Fig. 2.

SQUID-1 modulation patterns are shown in Figs. 2(a)-2(c) for different junction flux $\Phi_{JJ}$. Two junctions are tuned to have the same zero-field switching current so that the SQUID is symmetric ($I_{sw,a}=I_{sw,b}=0.9$ $\mu$A). At $\Phi_{JJ}=0$, the pattern exhibits kinks near SQUID oscillation minima [Fig. 2(a) yellow arrows, zoomed-in data in Fig. S9]. We also observe that, despite tuning the SQUID to the symmetric point, the nodes of the modulation patterns are lifted from zero. Pattern at $\Phi_{JJ}=0.5$ [Fig. 2(b)] does not contain extra features. Here, $I_{sw}$ follows a $|\cos B|$ curve as expected for standard symmetric SQUID without high-order harmonics in its CPR. When $\Phi_{JJ}$ is 0.9, the pattern has a distinct double-modulation character. Extra minima in $I_{sw}$ are observed near half-period of SQUID modulation [Fig. 2(c) yellow arrow].

Next, we extract $I_{sw}(\Delta B)$ traces from data like in Fig. 2 and fit them to a basic model for a SQUID with two-harmonic CPR junctions $I(\varphi)=I_{1}\sin(\varphi)+I_{2}\sin(2\varphi)$ (see supplementary materials for model description and detailed parameters [44]). The numerical model uses $I_{2}/I_{1}$=0.4 the ratio between the second and the first Josephson harmonics in the CPR, for both junctions a and b.

In Fig. 3(a), the device is tuned to a symmetric state where the zero-field switching currents in JJ_a and JJ_b are similar ($I_{sw,a}=I_{sw,b}=0.9$ $\mu$A). As the junction flux $\Phi_{JJ}$ increases, the amplitude of $I_{sw}$ shifts towards lower values due to the Fraunhofer-like envelope. The model reproduces the basic features. Kinks, pointed out by single-ended arrows in panel (a), are most apparent near junction flux $\Phi_{JJ}=0$ (zoomed-in data in Fig. S9). Node lifting (double arrows) is also reproduced by the model. While it is typically associated with asymmetric SQUIDs, in this model it originates from the second-order Josephson effect. Additional features disappear both in the data and in the model near $\Phi_{JJ}=0.5$, where the pattern closely follows the standard $|\cos B|$ SQUID relation. When $\Phi_{JJ}$ approaches 0.8 the experimental curve becomes skewed while the simulated curve is more symmetric.

In Fig. 3(b) SQUID-1 is tuned by gates, and is asymmetric while $\Phi_{JJ}$ is 0. $I_{sw,b}$ is fixed at $1.47$ $\mu$A while $I_{sw,a}$ is changing. The amplitude of the oscillation decreases as $I_{sw,a}$ decreases. The oscillations are asymmetric. At $I_{sw,a}=1.2$ $\mu$A, two kinks (arrows) are at different $I_{sw}$ heights. As $I_{sw,a}$ decreases, one kink becomes difficult to resolve. In this case skewed patterns are captured by the model, e.g. at $I_{sw,a}=0.2$ $\mu$A. More examples of SQUID data are presented in supplementary materials from this device and additional SQUIDs [44].

We first discuss the situation where $\Phi_{JJ}=0$. Our SQUID model shows that for a symmetric SQUID at $\Phi_{JJ}=0$, the nodes in $I_{sw}$ are lifted from 0 and two additional kinks appear as $I_{2}$ increases [Fig. S3(a)]. These two signatures are highlighted in Fig. 3(a). Thus the presence of additional kinks, and the node lifting are consistent with a significant second harmonic.

As the SQUID is made more asymmetric in Fig. 3(b), the SQUID modulation appears more skewed, meaning that the y-axis position of one kink decreases while the position of the other increases. When $I_{sw,a}<<I_{sw,b}$ , $JJ_{b}$ passes a nearly fixed supercurrent (thus a constant phase difference) to maximize the total switching current $I_{sw}$ while the oscillation comes mostly from $JJ_{a}$. As a result, $I_{sw}(\Delta B)$ traces out the CPR curve of $JJ_{a}$ with a constant shift. This regime is reached for $I_{sw,a}=0.6$ $\mu$A and 0.4 $\mu$A. We notice that at $I_{sw,a}=0.2$ $\mu$A the kink is not resolved in the experimental curve but persists in the simulated curve. This may be due to a decrease in the junction transparency when $I_{sw,a}$ is small which decreases the amplitude of the second harmonic in the CPR, but can also be due to the decreased accuracy of extraction at lower currents.

Near $\Phi_{JJ}=0.5$, $I_{sw}$ oscillates like a $|\cos{B}|$ function (Eq. S8) and shows no kinks or extra minima. The device looks like an ordinary SQUID. This is because at this junction flux value the second harmonic within each junction is cancelled. While the first harmonic has diffraction minima at $\Phi_{JJ}=1$, the second harmonic has nodes at $\Phi_{JJ}=0.5$ where $I_{2}$=0.

The deviation from simulations at $\Phi_{JJ}=0.8$ in Fig. 3(a) and $\Phi_{JJ}=0.9$ in Fig. 2(c) can be explained by a small difference between two JJs’ Fraunhofer periods. The actual flux is closer to 1 in one junction than the other when $\Phi_{JJ}$ is near 1. This difference makes the SQUID more and more asymmetric near $\Phi_{JJ}=1$.

We also find second Josephson harmonic signatures in single junctions (Fig. 4). Similar results from other single JJs including those made from Sn/InAs 2DEG are available in Refs. [44, 42]. The overall envelope $I_{sw}(B)$ is Fraunhofer-like, which is standard for single-harmonic uniform and sinusoidal junctions. However, in deviation from this behavior, kinks are observed at half-integer magnetic flux values [Fig. 4(a) black arrows]. The kink is most clear within the first lobe of the Fraunhofer-like modulation. The second lobe is skewed which is consistent with a minor kink. The pattern is not symmetric in positive-negative field but is inversion symmetric in field-current four-quadrant view. This is consistent with the self-field effect in these extended junctions.

We extract $I_{sw}$ and fit it with three different models in Fig. 4(b), detailed in supplementary materials [44]. The solid blue trace shows the basic Fraunhofer diffraction pattern for a single-component CPR. The two-component CPR model reproduces the kinks as shown in solid orange trace. Another model shown with dashed green line is for non-uniform critical current density and is discussed in the Alternative Explanations block.

Both non-uniformly distributed critical current and second harmonic in CPR produce half-period kinks. We notice that measured $I_{sw}$ of the single junction is smaller than the simulated value near zero field. This may be due to a non-uniformly distributed supercurrent or the fact that the experimentally measured switching current is smaller than the theoretical critical current [45]. Because the fit does not reproduce the data simultaneously in the central lobe and in the side lobes, we do not rely on the values of $I_{2}/I_{1}$ extracted from this analysis (see supplementary materials for details [44]). Half-periodic kinks were observed in anodic-oxidation Al/2DEG junctions but the analysis of the second-order Josephson effect has not been performed [46].

Shapiro steps in JJ-1 are presented at a microwave frequency of 6.742 GHz. A Shapiro step is a minimum in $dV/dI$ [Fig. 5(a)]. Another common form of plotting Shapiro step data is to bin data points by measured voltage across the junction [47]. In this form steps become straight bright lines in the histogram at Josephson voltages [Fig. 5(b)]. Apart from the usual integer steps, half-integer Shapiro steps are also observed at 3/2, 5/2 and 7/2. The missing step at 1/2 is likely due to self-heating [48, 49, 50].

While half-integer steps are clear and observed at relatively high frequencies, we do not observe steps at higher denominators such as 1/3 or 1/4. This is another confirmation that higher order Josephson harmonics, such as 3rd and 4th order, are not apparent in the data. Steps at 1/3 have been reported before in InSb nanowire junctions [51].

In these junctions we also observe missing integer steps at lower microwave frequencies, discussed in a separate manuscript [50]. This has been reported as a signature of the 4$\pi$ Josephson effect characterized by the $\sin(\varphi/2)$ CPR. However, we observe the missing step pattern in the non-topological regime where fractional Josephson effect is not expected. We explain this through a combination of fine-tuning and low signal levels.

To illustrate that missing Shapiro steps are weak as evidence of unusual features in the CPR, we draw attention to the fact that the step at $n=1/2$ is missing from our data. We do not take this as evidence of an “integer” Josephson effect.

Asymmetries between negative and positive switching currents in SQUIDs. The large second-order Josephson effect can produce asymmetries between positive and negative bias switching currents in junctions and SQUIDs [52, 53, 14, 34, 54]. These phenomena are heavily studied under the name “superconducting diode” in recent literature. To understand how the asymmetry arises we provide a simulation of our 2-junction 2-component CPR model in positive and negative bias (Figs. 6(a)). As can be seen, a kink in the positive bias is aligned with a dip in negative bias, and vice versa, an effect entirely due to the two-component CPR. This behavior is also present in the experiment (Fig. 6(b)). At the same time, the experimental data differ from the model in that the switching current maxima are not aligned for positive and negative bias in the experiment, but they are in the model. This difference is due to the finite inductance of the SQUID loop, of order 160 $\mu$H, resulting in finite phase winding due to circulating currents, an inductive effect [55, 56]. This effect provides another mechanism for the asymmetry. The simulation including both the inductance and the second Josephson harmonic shows that increasing the inductance would suppress signatures due to the second harmonic, e.g., the kinks at quarter flux values (Fig. S6).

Sign of the second harmonic term. Simulation results in the main text do not include the inductive effect, therefore they do not reveal the sign of the second harmonic term. This is because $I_{sw}$ is unchanged if both signs of the second harmonic term and the external flux are flipped [44]. The presence of finite SQUID inductance provides a simple way to determine the sign of the second harmonic term. The model combining both second harmonic and inductive effects shows that maxima in negative and positive $I_{sw}$s move towards kinks near half-integer flux values if the second-harmonic term is negative, and vice versa (Figs. S6(b) and S7(b)). The movement of $I_{sw}$ maxima in the experiment (Fig. 6(b)) suggests a negative second harmonic term. This is expected if the second-order effect originates from high-transparency of the junction and from the skewed current-phase relation, which yields a negative second-order Fourier component. Note that the sign of the second-order term can also be determined by checking the directions of of the applied field and current flow.

Half-periodic kinks in a single JJ diffraction patterns can be explained by a non-uniformly distributed supercurrent. We reproduce these kinks with the model of a three-level supercurrent density distribution (Fig. 4(b) and supplementary materials [44]). The non-uniform supercurrent may be caused by resist residue introduced during fabrication. We observe stripes of residue in our first batch (Tab. S1, Al-chip-1) of junctions which includes JJ-1, caused by double electron beam exposure (Fig. S1). These residues are avoided and are not observed in newer devices without double exposure (Figs. 1(b),S23(b),S23(c)), which includes SQUID-1 and more devices in Fig. S23 and Ref. [42] - many of which do show additional modulation such as extra nodes and kinks.

It is worth noting that the non-uniform model produces a “lifted odd node” diffraction pattern that has been used to argue for the observation of the exotic fractional Josephson effect, characterized by 4$\pi$-periodic current phase relations associated with Majorana modes [57]. We obtain the diffraction pattern of the same overall shape without the need to have this physics present [Fig. 4(b)]. In planar junctions it is difficult to know whether a node corresponds to integer or half-integer flux. This is because the junction area for supercurrent diffraction purposes can be much larger than the lithographic area. In our work, for instance, the junction length is close to 1.5 microns while the lithographic length is $\sim$150 nm. This makes it hard to distinguish a half-integer kink from a lifted integer node.

Another aspect we consider is whether patterns such as those extracted in Fig. 3 or Fig. 4(b) are true representations of the switching current evolution. In some of the junctions the apparent extra modulation (blue area in colorplots, highlighted in Fig. S24(c) by changing the color) appears to be above the true switching current and may arise due to the evolution of finite-voltage resonances in magnetic field. One origin for these resonances is multiple Andreev reflections (MAR). See Figs. S15-S18 for examples.

This type of artefacts does not explain all of our data. For example, in Fig. 4(a) we do not see MAR or faint finite-voltage state regions. Another argument we can give in support of the $\sin(2\varphi)$ origins of the SQUID patterns is the detailed agreement between our numerical model and the data. See supplementary materials for an extended discussion [44].

Half-integer Shapiro steps have been attributed to effects not related to the second-order Josephson effect, such as a phase-locking of Josephson half-vortices to microwaves, and non-equilibrium quasiparticle dynamics [40, 41]. Hence their demonstration cannot by itself be used to claim the presence of $\sin(2\varphi)$ terms. We provide Shapiro data as additional rather than key evidence.

Taken together, we argue that the evidence from SQUIDs, single junctions and agreement with the model, are all consistent with a strong second-order Josephson effect. While alternative explanations can be used to support individual measurements, none of them can explain all of the data, though multiple factors may be present simultaneously with a lower likelihood.

In summary, we studied the current-phase relation in Josephson junctions and SQUIDs made with the nanowire shadow-mask method. Signatures due to a strong second-order harmonic in the CPR are observed. The simulation shows good agreement to experimental $I_{sw}$ data in the SQUID device suggesting a ratio of the second to first harmonics of $I_{2}/I_{1}=-0.4$.

In the clean limit and at zero temperature, the CPR is a skewed 2$\pi$-periodic function [6] with many higher order terms. After Fourier decomposition, one gets $|I_{2}/I_{1}|=0.4$ and $|I_{3}/I_{1}|\approx 0.26$. So, in an ideal ballistic junction, the magnitude of the second harmonic we obtain would not necessarily be surprising. The surprising is the combination of doubly-modulated junction characteristics together with the lack of manifested higher order terms such as the third term. This indicates an unusual situation in which either the second order term is unusually large, or the higher order terms are suppressed through an undetermined mechanism.

The origin of such a strong second-order harmonic, not accompanied by even higher order terms (skewed CPR) still needs to be understood. The sensitivity of observed signatures to surface treatment can be studied. Higher harmonics can be used to engineer nonlinearity in quantum circuits.

This study is divided into two periods, which correspond to period 2 and period 3 in Ref. [42].

The first period was between August 2018 to June 2019, including sample preparation, device fabrication and measurements. We explored the first-generation devices (Al, InSb nanowire, without gates). More than 7 devices on 1 chip are measured during 2 cooldowns in a dilution refrigerator, producing about 8900 datasets.

The second period was between March 2021 to February 2022. We explored the second-generation devices, (Al or Sn, InAs/HfO_x nanowire, with self-aligned nanowire gates). 62 devices on 6 chips are measured during 8 cooldowns in dilution refrigerators, producing about 5700 datasets.

Curated library of data extending beyond what is presented in the paper, as well as simulation and data processing code are available at [58].

We thank J. Stenger, D. Pekker, D. Van Harlingen for discussions. We thank E. Bakkers, G. Badawy and S. Gazibegovic for providing InSb nanowires. We acknowledge the use of shared facilities of the NSF Materials Research Science and Engineering Center (MRSEC) at the University of California Santa Barbara (DMR 1720256) and the Nanotech UCSB Nanofabrication Facility.

Work supported by the ANR-NSF PIRE:HYBRID OISE-1743717, NSF Quantum Foundry funded via the Q-AMASE-i program under award DMR-1906325, the Transatlantic Research Partnership and IRP-CNRS HYNATOQ, U.S. ONR and ARO.

A.-H.C, H.W., and M.H. grew InAs nanowires and the dielectric layer. M.P., J.S.L., C.P.D., A.P.M., S.D.H., J.T.D., and C.J.P. grew quantum wells and superconducting films. A.Z. and P.Z. fabricated devices. L.J., V.V.d.S., A.Z., and P.Z. performed measurements. L.J. and P.Z. did the simulation. A.Z. and L.J. prepared draft versions of this manuscript. P.Z. and S.M.F. wrote the manuscript with inputs from all authors.

The magnetic field induces an extra phase variation inside the junction (Eq. S3), leading to the supercurrent interference. The switching current ($I_{sw}$) of a single JJ in a magnetic field can be written as

where $\Phi_{JJ}$ is the magnetic flux in the junction normalized by the superconducting flux quantum ($\Phi_{0}=h/2e$), $\varphi_{0}$ is the phase difference which is a free parameter here for calculating $I_{sw}$, $W$ is the width (perpendicular to the direction of the current) of the junction, $I_{1}$ and $I_{2}$ are amplitudes of the first- and second-order harmonics in the CPR ($I_{n}$ can be regarded as the global current and $I_{n}/W$ as the current density), we ignore higher-order harmonics in our model, $y$ is the position in the $W$ direction. The integration in Eq. S2 is independent of $W$, so we can set $W=1$ for simplicity.

If $I_{1}$ is a constant and $I_{2}=0$, Eq. S2 can be simplified as

and $I_{sw}$ reduces to $\left|I_{1}\mathrm{sinc}\left(\pi\Phi_{JJ}\right)\right|$ which resembles the Fraunhofer diffraction.

Three models for JJs are used in Fig. 4(b):

1. 1.

   CPR has only the 1st harmonic, $I_{1}\neq 0$, $I_{2}=0$.

2. 2.

   CPR has 1st and 2nd harmonics, $I_{1}\neq 0$, $I_{2}\neq 0$.

3. 3.

   CPR has only the 1st harmonic and $I_{1}$ is non-uniformly distributed. In this model, $I_{1}(y)$ is a three-step piece-wise function. It equals to $j_{\text{center}}$ in the middle ($W/3<y<2W/3$) and $j_{\text{side}}$ on two sides.

As sketched in Fig. 1(a), a SQUID consists of two JJs, i.e., JJ_a and JJ_b. We denote magnetic fluxes (normalized by $h/2e$) in the enclosed area of the SQUID, JJ_a, and JJ_b by $\Phi_{SQUID}$, $\Phi_{JJ,a}=\Phi_{SQUID}/r_{a}$ and $\Phi_{JJ,b}=\Phi_{SQUID}/r_{b}$. Here $r_{a}$ and $r_{b}$ are ratios between the SQUID-enclosed area and junction areas. We denote amplitudes of harmonics in JJ_a and JJ_b by $I_{1,a}$, $I_{2,a}$, $I_{1,b}$, $I_{2,b}$.

Similar to Eq. S1-S3, the switching current of a SQUID can be calculated by the following equations

where $\delta=2\pi(\Phi_{JJ,a}+\Phi_{SQUID})=2\pi(1/r_{a}+1)\Phi_{SQUID}$. $W_{a}$ and $W_{b}$ are junction widths. The integration terms are independent of widths because $W_{a}$ and $W_{b}$ are absorbed under the transformation $y\rightarrow y/W_{i}$.

In the simplest situation, i.e, $r_{a},r_{b}>>1$, two junctions are identical, and there is only the first harmonic in the CPR, $I_{sw}$ reduces to

which is a high-frequency SQUID oscillation (the cosine term) modulated by a low-frequency Fraunhofer oscillation ($I_{sw,JJ}$).

In our simulation, we assume $I_{2,i}/I_{1,i}$ is a constant independent of the junction index $i$ and the gate voltage. Parameters used for SQUID-1 in Fig. 3 are as follows, $r_{i}=185$, $I_{1,i}=\alpha I_{sw,i}$, $I_{2,i}=0.4\alpha I_{sw,i}$, $i\in\{a,b\}$. $r_{i}$ is extracted from periods of the JJ oscillation and the SQUID oscillation. $I_{sw,i}$ is the measured zero-field switching current in junction $i$ [Fig. 1(c)]. $\alpha$ is a fitting parameter which is chosen to be $0.91$ in Fig. 3(a) and $0.82$ in Fig. 3(b). The difference in $\alpha$ may arise due to the uncertainty in tuning nominal switching currents by gates or higher-order harmonics that are not considered in the simulation.

We model the inductive SQUID following Ref. [56], but including higher order terms. The currents through junctions a and b are:

where $\varphi_{a}$ and $\varphi_{b}$ are phase differences in junctions a and b, respectively. $\gamma$ is the amplitude of the second harmonic. For simplicity, we use normalized currents, $i_{a}=I_{a}/I_{0}$ and $i_{b}=I_{b}/I_{0}$. The total current $i$ and circulating current $j$ are:

$\varphi_{a}$ and $\varphi_{b}$ are connected by the equation:

where $\phi_{ext}=\Phi_{ext}/\Phi_{0}$ is the normalized flux applied by the external field, $\beta=2LI_{0}/\Phi_{0}$ is the normalized inductance, $n$ is an arbitrary integer. Here we ignore the inductance difference between the two arms of the SQUID for simplicity. A more general case can be found in Ref. [56]. Note that when $\gamma\neq 0$, $I_{0}$ is different from the critical current by a factor $f(\gamma)$. We should use $\beta^{\prime}=(I_{c,a}+I_{c,b})L/\Phi_{0}=f(\gamma)\beta$ instead of $\beta$ to compare with the experimental parameters.

The maximized (minimized) $i$ is achieved when both $i_{a}$ and $i_{b}$ are maximized (minimized), which gives $\varphi_{a}=\varphi_{b}$, $i_{a}=\pm(1-\alpha)f(\gamma)$, $i_{b}=\pm(1+\alpha)f(\gamma)$, $+$ ($-$) for the maximum (minimum). The external fluxes where $i$ reaches maximum (minimum) can be calculated by substituting these into Eq. S13:

The dependence of critical current on $\phi_{ext}$ is calculated using the following procedure. First, at every $\phi_{ext}$, we search for $(\varphi_{a},\varphi_{b})$ pairs satisfying Eq. S13. Second, we find the minimum and maximum currents among valid $(\varphi_{a},\varphi_{b})$ pairs, which are the negative and positive critical currents. The results of inductive modeling can be found in Figs. S5 and S6.

In the JJ model and the non-inductive SQUID model, changing the sign of the second harmonic term is equivalent to changing the sign of the external flux. This is because Eqs. S2 and S6 are unchanged under the transformation $\{I_{2},\Phi_{JJ},\varphi_{0}\}\rightarrow\{-I_{2},-\Phi_{JJ},\pi-\varphi_{0}\}$ and $\{I_{2,a},I_{2,b},\Phi_{SQUID},\varphi_{0}\}\rightarrow\{-I_{2,a},-I_{2,b},-\Phi_{SQUID},\pi-\varphi_{0}\}$, respectively.

In the inductive SQUID model, changing the sign of the second harmonic term is equivalent to changing both signs of the external flux and the normalized inductance $\beta$. This is because Eqs. S9 and S10 are unchanged under the transformation $\{\gamma,\phi_{ext},\beta,\varphi_{a}\}\rightarrow\{-\gamma,-\phi_{ext},-\beta,\pi-\varphi_{a}\}$.

The results for flipping the signs of coefficients in above models can be found in Figs. S7 and S8.

.

.