Van der Waals ferromagnetic Josephson junctions

Superconductor-ferromagnet interfaces in two-dimensional heterostructures present a unique opportunity to study the interplay between superconductivity and ferromagnetism. The realization of such nanoscale heterostructures in van der Waals (vdW) crystals remains largely unexplored due to the challenge of making atomically-sharp interfaces from their layered structures. Here, we build a vdW ferromagnetic Josephson junction (JJ) by inserting a few-layer ferromagnetic insulator Cr2Ge2Te6 into two layers of superconductor NbSe2. The critical current and corresponding junction resistance exhibit a hysteretic and oscillatory behavior against in-plane magnetic fields, manifesting itself as a strong Josephson coupling state. Also, we observe a central minimum of critical current in some JJ devices as well as a nontrivial phase shift in SQUID structures, evidencing the coexistence of 0 and π phase in the junction region. Our study paves the way to exploring sensitive probes of weak magnetism and multifunctional building-blocks for phase-related superconducting circuits using vdW heterostructures.


1.
Additional transport results of device #01

. Josephson characterizations and Fiske steps
Additional transport properties of device #01 are shown here. Supplementary Fig.  1a is the temperature dependence of junction resistance (black) from 0.05 K to 300 K, and the calculated residual resistance ratio (RRR) at 10 K compared to 300 K is 0.15, which is of a similar magnitude to pure NbSe 2 (~ 0.12). A previous study in NbN/GdN/NbN JJ 1 showed that its junction behavior at high temperature is typical of a semiconducting barrier and then it turns into a continuous decrease in resistance below T Curie of GdN, due to the exchange splitting of the tunnel barrier with reduction of effective barrier height for one spin sign. However, the decrease of resistance in our device is dominated by NbSe 2 themselves. We have summarized the data of all devices, for each performs a metallic behavior before its superconducting transition temperature. In line with the small hysteresis of the IV curves as well as small R n , we have to address that the Cr 2 Ge 2 Te 6 has a small gap that can be doped easily by the environmental charge impurities. There's a high possibility that at the base temperature there is still some residual spin-polarized density of states (DOS) in the tunneling barrier itself. Supplementary Fig. 1b presents the IV curves measured at the different magnetic fields as labeled in the inset, and the arrows indicate two current steps at bias voltage because of the Fiske resonance. In device #01, we can see a first-order Fiske step in which the positions of the 2nd maximum in I c correspond to the 1st minimum. Besides, a second-order Fiske step occasionally appears in device #02, as shown in the main text Fig. 4a. The voltage position of the first Fiske step V 1 is 0.57 mV both indicated in the 4 mT (yellow) and 9.5 mT (blue) curves. The two-gap edge state of JJ device #01 is obtained from the integral of the dV/dI curve, and the value V 2 of ~ 2 mV is the same as the SIS tunnel junction results in the following sections.

Supplementary Fig. 1. Additional JJ characteristics of device #01. (a)
Temperature dependence of resistance of the junction (black) and NbSe 2 (red) from 0.05 K to 300 K, which shows a metallic barrier behavior above T c . (b) Differential conductance (blue) along with charge current (red) as a function of the bias voltage. The two-gap edge state value is ~ 2 mV, twice the gap of the superconductor electrodes (NbSe 2 ). Indeed the phenomenon of Fiske resonance is more essential driven by the voltage rather than the current. According to the Josephson equations, the supercurrent across a JJ oscillates at a constant frequency ω J = 2πV/Φ 0 at an applied voltage V 2 . When the oscillation frequency matches the nth harmonic of the JJ cavity mode as the spatial period of the Josephson current distribution matches the spatial period of the nth resonant electromagnetic mode of the junction, achieved by applying a magnetic field, it will result in the excitation of the cavity modes as the self-induced current steps at specific voltages V n denoted as Fiske steps [3][4][5] . Regrettably, we have used the current bias to drive the JJ in our experiments and mainly focused on its critical current, therefore we did not try the voltage-driven operations to reproduce the Fiske resonance for a comparison. Nevertheless, we try to integrate the current into the voltage coordinate by the differential resistance dV/dI to reveal the complementary details of the observed Fiske steps results.
The comparison between the Fiske steps' configuration driven by the current and the voltage is presented in Supplementary Fig. 3a. We can see that in the lower panel, the shaded areas indicated by the black and white dashed circles represent the specific excitation voltages V 1 and V 2 corresponding to the first and second-order Fiske resonances respectively. Although the original current bias acting as the coordinate has relatively higher accuracy in determining the peak positions of dV/dI, we can only roughly define the resonant modes in its integration plots because of the highly damping environment in our JJs as Supplementary Fig. 3b showed, where the resonances are not sharp enough at the almost constant voltages. The value of V 1 is defined as the valley position of dV/dI vs V at 0.57 mV ( Supplementary Fig. 4b), and the consequent specific capacitance C s and the parameter β c mentioned in the main text can be calculated for estimations. The first and the second-order Fiske steps are labeled by the black (V 1 ) and white (V 2 ) arrows with slight broadening due to the errors of integration. (b) dV/dI as a function of current (black) and its integrating voltage (red) measured at zero magnetic field.
Next, more characterizations are presented to confirm the flux dependence of the resonant modes, in which we select five line cuts of the voltage-driven plots at the different magnetic fields in the case of successive periods of its diffraction pattern, indicated as Supplementary Fig. 4a. The valley positions of the curves measured at 4, 9.5 and 15 mT shown in Supplementary Fig. 4c are nearly the same regardless of slight shifts, accounting for the clear first-order Fiske resonance at around 0.57 mV. The second-order Fiske resonance is observed in the curve of 12.5 mT at around 0.97 mV while the first-order step disappears. Compared to the single resonant mode results, they simultaneously occur in the curve of 18.5 mT shown in Supplementary  Fig. 4d, and the resonance voltages are very close to the values above. Therefore we can find it easier to identify the Fiske resonances in the voltage-driven plots as the reviewer suggested. To precisely determine the resonant voltages, we expect further improvements in the devices' qualities for negligible damping circumstance and higher-order Fiske steps, under which the internal losses in our JJs can be revealed more accurately.

Part 2. Magnetic response at low temperature
A detailed response to the parallel magnetic field in device #01 is shown in the next parts. Supplementary Fig. 5a shows the magnetoresistance results measured at 1.9 K. For the junction, its resistance increases below 60 mT and eventually decreases up to 9 T when spins in the intermediate Cr 2 Ge 2 Te 6 were aligned to the field direction. In the meanwhile, the Ising superconductors NbSe 2 with a large parallel critical field (at the top and bottom) were still in the superconducting state with zero resistance, suggesting that the entire contribution of junction resistance comes from the tunneling barrier. In Supplementary Fig. 5b we show that the illustrations of barrier moment as a sequence of applied field split into three stages from the junction resistance, which corresponds to a complete hysteresis magnetization loop of magnetic barriers.
During the experiments, we have used two measurement systems to carry out the transport data in different temperature ranges, therefore a comparison between them is needed. From the arrows in Supplementary Fig. 5c we can see the central positions of two sweep branches are slightly shifted, meanwhile, the oscillations at 0.1 K are evident, compared to 1.9 K due to fewer thermal perturbations. Nevertheless, the magnitude of hysteresis is almost the same at these two temperatures (for one starts from 30 mT and the other from 27 mT), which also proves the validity that we use ΔB max as an index for estimation from 7 K to 0.1 K.

Supplementary Fig. 5. Additional magnetic response of device #01. (a)
Magnetoresistance of the junction (black) and pure superconductor layers (red for the top and blue for the bottom) versus parallel magnetic field up to 9 T. (b) Oscillatory resistance patterns of the junction at 1.9 K. The black line is measured after zero-field cooling to induce an initial magnetization of the barrier, and the blue and red curves are responses in downward and upward directions of field sweeping, respectively. (c) Comparison of hysteretic resistance patterns between the two measurement systems in different temperature ranges (0.1 K in a dilution refrigerator and 1.9 K in a PPMS).
Since the hysteresis introduced by the magnetic barrier in JJs can be controlled by B max , we show a stacked plot of hysteretic junction resistance with B max ranging from 20 mT to 100 mT in Supplementary Fig. 6a. The coincidence of the opposite direction sweeping curves over 75 mT shown in the bottom indicates a saturated magnetized state, in contrast to smaller B max situations that the shift between the two curves arises immediately as the field direction turns back. Consequently, with a step-increase of B max , the central lobe (equivalent to a zero-flux state in the junction) moves as guided by arrows in Supplementary Fig. 6a, and the variation of its peak positions at downward or upward sweeping is presented in Supplementary Fig. 6b (1L represents the 1st maximum resistance at the left side of the central lobe while 1R at right side). Supplementary Fig. 6c depicts the separations between the two sweeping branches at different B max , from which we can clearly see a rapid increase in measured hysteresis since the initial magnetization, and then it gradually converges above 75 mT, the same as the conclusion reached from Supplementary Fig. 6a.

(a-c)
Junction resistance (blue for downward and red for upward) measured at 6 K, 6.2 K, and 6.5 K. The discrepancies in the center positions of the curves tend to be illegible as the temperature approaches T c , which causes a sudden drop of ΔB max value as the data shown in the main text Fig. 3(d). When the temperature is above T c ( ~ 6.4K), measured hysteresis is almost zero within the error of field resolution.

Supplementary Note 2. Discontinuity and irreversibility in I c patterns Part 1. Discontinuity in device #01
The breakpoints that appeared in I c patterns arise with the progressive magnetization of a highly polarized barrier, manifested as the abrupt jumps of I c evidently at higher field. In Supplementary Fig. 9, we first show the breakpoints measured in device #01 and the dV/dI curves on both sides of the breakpoint near 25 mT. The conspicuous discontinuities above 23 mT are in sharp contrast to the patterns at the lower field region (the single step of the applied field is 0.5 mT in the whole process). Nevertheless, the results in the higher field become less recognizable because of the non-zero junction residual resistance at zero bias, thus we cannot get a valid I c value there.  Supplementary Fig. 9a). The I c measured at 24.5 and 25 mT are almost identical (10.5 μA), but then encounter a sudden drop to almost zero as measured at 25.5 and 26 mT. The step of field resolution is 0.5 mT.

Part 2. Discontinuity and irreversibility in device #02
In device #02, the zero resistance state persists till -140 mT, in case that we can observe its first appearance of breakpoints at around -90 mT ( Supplementary Fig.  10a). When the field reverses back to zero, the breakpoints still maintain their dominance even in the range where I c has previously varied continuously ( Supplementary Fig. 10b). Apart from the interlayer coupling in the Cr 2 Ge 2 Te 6 barrier which prevents the flip back of spins as we have discussed in the main text, the assumptions on multi-domains state possessing discrepant coercivity in the junction region may also be suggested to explain the step-like behaviors in I c .

Supplementary Note 3. Excess Josephson current in thicker barrier device
For a thinner barrier like device #01 (~ 6.5 nm), the tunneling behavior shows a good coincidence to the AB relation, whose general form implies an invariance of I c and the R n depends on the measured superconducting energy gap. Note that the R n is ~ 10.6 Ω obtained from the linear part in its IV curves and the fitting results are shown in Supplementary Fig. 11a. The calculated I c R n product (V c ) at zero temperature is ~ 0.51 mV, which the magnitude is much smaller than the theoretical value of πΔ 0 /2e = 1.57 mV (bulk NbSe 2 energy gap: Δ 0 ~ 1 meV). The strong suppression shows that the tunneling probability of pairs in a magnetic barrier is lower than that for single electrons.
However, in a thicker barrier device #03 (~ 11nm), we find an apparent excess Josephson current substantially deviated from the AB relation at low temperature. The data and fitting results (with a data range of 4.6 K to 6 K) are shown in Supplementary Fig. 11b. The I c R n product obtained at the base temperature of 0.5 K is 1.43 mV, exceeding the diffusive limit from the AB relation fit at zero temperature (0.84 mV). We tend to attribute it to a possible multi-domain structure in few-layer Cr 2 Ge 2 Te 6 especially for thicker layers, and the introduced magnetic inhomogeneity by the relative size of domain walls implies a reduction of exchange energy in the ferromagnets, which accounts for the enhancement. Supplementary Fig. 11. Multi-domains induced excess Josephson current in the thick barrier device. (a) Temperature dependence of I c R n (V c ) product of device #01 (black open circles) with a 6.5 nm barrier, and the red solid line is the fit to AB relation for an insulating barrier which gives V c0 = 0.51 mV. The theoretical value is ~ 1.57 mV for NbSe 2 JJ. (b) Temperature dependence of I c R n product of device #03 with an 11 nm barrier and the largest V c acquired at 0.5 K is 1.43 mV. The red solid line is an AB relation fit from 4.6 K to 6 K which gives V c0 = 0.84 mV. The deviation between measured data and the fit is originated from the multi-domain structure in few-layer Cr 2 Ge 2 Te 6 especially for thicker samples, where the magnetic inhomogeneity can induce an equivalent reduction of exchange energy E ex .

Supplementary Note 4. Thickness-tuned crossover from Josephson coupling to quantum tunneling
In addition to the JJ devices with relatively thin Cr 2 Ge 2 Te 6 as barriers, we have also investigated the crossover from the Josephson coupling state to the quantum tunneling dominant situation tuned by the magnetic barrier thickness. For a thick insulating layer, the Cooper pair wavefunction Ψ penetrating the JJ has damped to zero in the middle of the barrier that cannot provide a stable phase coherence, thus the transport conductance of the junction here is dominant by the quantum tunneling effect through single-electron channels 6 (Supplementary Fig. 12a).
In the left inset of Supplementary Fig. 12b, we first display the normalized differential tunneling conductance (G s /G n ) at 2 K for a 14.5±0.4 nm thick Cr 2 Ge 2 Te 6 junction device #05. The separation between the peaks corresponds to the superconducting gap of bi-lateral few-layer NbSe 2 (20 ~ 40 nm), whose properties are similar to the bulk NbSe 2 . For the quantitative analysis of the tunneling spectra, we   Fig. 12b) is described by the BCS theory, and the extracted zero-temperature gap value is 2Δ 0 = 1.98 meV quantitatively consistent with the value 3.53 k B T c (k B is the Boltzmann constant and T c ~ 6.7 K is consistent with the resistance measurements). Dimensionless effective barrier strength Z is 0.87 at 2 K, which is comparable to other vdW tunneling materials 7,8 . The magnetic response of the tunnel junction versus the in-plane field is shown in Supplementary Fig. 12c, at zero bias voltage. The resistance is normalized to the value at different temperatures in the absence of field. The hysteresis shrinks in resistance amplitude and magnetic field range with increasing temperature, albeit it fully disappears above T c , similar to what we have observed in JJ devices. It's also noteworthy that the junction resistance near zero field is smaller than that in the higher field, and such a dip possibly accounts for a lower effective barrier height corresponding to zero magnetization compared to the fully-polarized state, which may also be influenced by the two spin-mixed S/F heterointerfaces to generate an anisotropic magnetoresistance 9,10 .

Supplementary Fig. 12. Quantum tunneling results of thicker barrier devices. (a)
Schematic views of the energy diagram of the SIS tunnel junction with the same superconductors, showing that the conductance as a function of bias voltage measures the BCS superconducting gap 2Δ. (b) Normalized differential tunneling conductance (G s /G n ) at 2K under zero field for a thick barrier device #05 (~ 14.5nm), and the BTK fitting gives the zero temperature BCS gap 2Δ 0 for NbSe 2 of 1.98 meV (black opened circles for experiment data and red solid lines for fit curves). (c) In-plane magnetoresistance of device #05 at 2 K, 4 K, and 8 K compared among the superconducting state and normal state.
Our results provide a barrier thickness-tuned crossover from quantum tunneling to Josephson coupling in NbSe 2 /Cr 2 Ge 2 Te 6 devices, nevertheless, we cannot confirm that it is a continuous process. Both of the effective penetration length and surface impurities are required to accomplish an entire zero resistance state in a highly transparent and defect-free junction, therefore we plot a summary of normalized residual resistance (R 2K /R 10K ) varied with barrier thickness in the range of 6 -15 nm. The results are divided into three regions: R 2K /R 10K > 1, R 2K /R 10K = 0 and intermediate step 0 < R 2K /R 10K < 1. All the devices' resistance in the Josephson coupling regime (6 ~ 11 nm) continuously decreases when cooling down, in sharp contrast to the tunneling devices (13 ~ 15 nm), however, for the others (11 ~ 13 nm) the junction resistance saturates to a finite non-zero value at low temperatures. A careful comparison of their temperature dependences is shown in Supplementary Fig. 13a. We suppose that the junction resistance is formed by a series connection of two NbSe 2 electrodes and the non-superconducting Cr 2 Ge 2 Te 6 part, while the former is zero below T c and the latter is the main cause of residual resistance. It is suggested that the effective penetration length of Ψ in the dirty limit should be comparable to the critical thickness for Cooper pair tunneling in these magnetic JJs. Interfacial scattering and impurities will lead to a finite residual resistance value even in JJ samples of thinner barriers, as suppressions of penetration length which is equivalent to thicker barriers.
Supplementary Fig. 13. Summary of thickness-dependent tunneling process. (a) Temperature dependence of different devices with a thickness of 6.5, 10, 10.5, 12.5, and 14.5 nm, respectively. The residual resistance measured at 2 K increases from 0 to over 1 suggesting the crossover from the Josephson coupling state to the quantum tunneling dominant state. (b) Normalized residual resistance at 2 K when tuning the Cr 2 Ge 2 Te 6 barrier thickness, in considerations of interface impurities and defects.

Supplementary Note 5. Magnetic characterizations of Cr 2 Ge 2 Te 6
The birth of microscopic magnetic domain structure in Cr 2 Ge 2 Te 6 is related to the competition between the magnetocrystalline anisotropy and magnetic dipolar interaction, where the latter one always prefers an in-plane magnetization and becomes significant in low-dimensional materials 11 . A multi-domain state is favored in thick ferromagnetic films due to the dipolar energy winning over the exchange energy 12 , while in some cases for ultrathin flakes whose easy axis is perpendicular to the film, dipolar interaction is expected to be very small, probably leading to a very  large domain size as shown in a six-layer Cr 2 Ge 2 Te 6 from Gong et al's results 13 . The dependence of the domain size (D) should both take the variety of domain wall energy (U w ) and also the film thickness (L) comparable with the range of exchange interactions (r e ) into account. The model put forward by Kooy and Enz 14 predicted that the domain size is proportional to the film thickness when L is much larger than the characteristic 'dipolar length' (D 0 ), however as L is reduced below D 0 while still larger than r e , Kaplan and Gehring 15 found that the domain size will dramatically increase as film thickness decreases in the ultrathin regime. Taking the magnetic parameters of Cr 2 Ge 2 Te 6 as shown above, its domain wall energy density (σ w ) and the domain wall size can be estimated 16 . Also in the bulk Cr 2 Ge 2 Te 6 , Guo et al 25 have found various multiple domain structures and symmetries stem from the competition between thermal fluctuation and anisotropy energy induced by magnetic dipolar interaction through magnetic force microscope (MFM) measurements. From the known observations of multi-domains in Cr 2 Ge 2 Te 6 nanoflakes or bulks, the characteristic domain size is around 70~120 nm, which is much larger than the few-layer Cr 2 Ge 2 Te 6 nanoflake thickness (<10 nm). To verify the magnetic properties of our Cr 2 Ge 2 Te 6 flakes with those reported, we have further performed some characterizations to give a comprehensive understanding.
We first characterized the intrinsic magnetic properties of Cr 2 Ge 2 Te 6 bulk crystals, where the magnetization measurements are performed through a SQUID magnetometer. The Curie temperature T c ~ 63 K was determined from the drop in the temperature dependence of its magnetization when a 100 mT magnetic field is applied to the in-plane (along a-b plane) and the out-of-plane (along c axis) of the crystal, as shown in Supplementary Fig. 14a. The bulk magnetizations below T c with an external magnetic field reveal a ferromagnetic behavior with magnetic anisotropy by the comparison of their saturation fields measured at 2 K, as shown in Supplementary Fig.  14b. With a larger saturation field along the in-plane direction, the Cr 2 Ge 2 Te 6 crystal shows an easy axis perpendicular to the a-b plane, in agreement with previous reports 17,26 . Supplementary Fig. 14. Magnetic characterizations of bulk Cr 2 Ge 2 Te 6 by SQUID magnetometer. (a) Temperature dependence of the bulk Cr 2 Ge 2 Te 6 sample's magnetization with a 100 mT magnetic field applied along the out-of-plane (black) and in-plane (red) directions. (b) Field dependence of its magnetizations at 2 K (black curves for out-of-plane direction and red curves for in-plane direction).
To further confirm the magnetic behavior in exfoliated Cr 2 Ge 2 Te 6 nanoflakes, we performed the Kerr rotation measurements to explore their properties changing with layer thickness. The measured temperature is at 5 K, below the critical temperature of NbSe 2 which is used as a superconducting electrode in our van der Waals (vdW) Josephson junction (JJ) devices. The original Kerr rotation signals under the out-of-plane magnetic field are present in Supplementary Fig. 15a, where we can obtain that the saturation magnetizations of Cr 2 Ge 2 Te 6 nanoflakes increase as the film thicknesses increase. In Supplementary Fig. 15b we present the normalized curves extracted from Supplementary Fig. 15a for clarity in comparison, and we can clearly see that for thinner Cr 2 Ge 2 Te 6 samples like 4 nm and 12 nm flakes, they only show monotonic transitions in the Kerr rotation hysteresis loops. Besides, the magnetic anisotropy is more obvious in the thinnest 4 nm Cr 2 Ge 2 Te 6 . Yet a slanted behavior, as well as multiple transitions from the Kerr rotation loops under lower magnetic field, appears in the 20 nm Cr 2 Ge 2 Te 6 sample, suggesting the formation of nonuniform magnetizations in this thickest Cr 2 Ge 2 Te 6 flake, probably arisen from such as multi-domains. Combined the results reported before with our characterizations on Cr 2 Ge 2 Te 6 samples, both the larger domain size up to 10 μm as a single-domain remanence in ultrathin nanoflakes and the smaller domain size around 100 nm along with a multi-domain structure in much thicker Cr 2 Ge 2 Te 6 samples or even bulk crystals have been confirmed. Although we haven't found the counterexamples in these two thickness regions for Cr 2 Ge 2 Te 6 , we cannot conclude a solid correlation between the equilibrium domain size of Cr 2 Ge 2 Te 6 with its thickness, nor could we prove the tendency of forming multi-domains as the film thickness increases approaching a 3D ferromagnet up to our considerations. The basic purpose to give such an argument is trying to explain the enhancement of Josephson current in the thick-barrier device, therefore we will only limit our discussions to the range of Cr 2 Ge 2 Te 6 thickness below 20 nm according to our experiments.

Supplementary Note 6. Barrier thickness dependence of Josephson junction parameters
The transport property of the tunnel junctions can be obtained from the thickness dependence of junction resistance area product R n A, and the result of our experiments is shown in Supplementary Fig. 16a containing all the CGT devices from the Josephson coupling state to the tunneling region as defined in the supplementary information. Samples within the thinner barrier region but still with residual resistance are also included. The normal state resistances R n of all the junctions were measured at 8 K above T c . We find an exponential increase of R n A fitted to an exp(d F /t) with the barrier thickness d F ranging from 6-15 nm (guided by the red dashed line), where the characteristic quasiparticle tunnel length t is calculated to be 2 nm, slightly larger than the value obtained by Idzuchi et al 27 of 1.3 nm. However, the normalized barrier resistance constant is around 3-4 Ω μm 2 , which is much lower than their result of 340 Ω·μm 2 . The relatively small value of our junctions may partially be attributed to the decrease of the semiconducting gap of thicker CGT than the ultrathin flakes as 1-6L they used, yet, the leakage of charge carriers across the barrier is more likely to be the main cause. Furthermore, we summarize the thickness dependence of the R n A product, the critical current density J c and the I c R n product of the only JJs whose resistances have exactly dropped to zero at 2 K in our PPMS instruments. As in Supplementary  Fig. 16b, except for the thickest barrier sample with an enhanced supercurrent mentioned before in the main text, we find that the J c decreases with increasing d F in the rest devices, but the I c R n product does not show a monotonic dependence. Besides, neither the R n A product nor the J c shows an evidently exponential change as expected for an effective insulating barrier, leading to the moderate variations of the I c R n product. The relatively slow decrease of the J c can also support the assumption of charge leakage to explain why our CGT barrier is not so insulating as it should be. Finally, we have to admit the individual differences in our samples reflected in their transport properties, as it is hard to control all the variables in the fabrication process of these vdW JJs, and we are also inclined to ascribe some observations above to an unphysical excuse. At the same time, the data based on these 5 JJs are not adequate to entirely reveal the superconducting nature beneath, therefore more detailed experiments are needed especially for the ultrathin CGT barriers in further studies.
Supplementary Fig. 16. Barrier thickness dependence of the parameters. (a) The resistance area R n A product of the junction (black rectangles) with different barrier thicknesses. The fitting result of the data by the exponential relation is indicated by the red dashed line. (b) The R n A product (black rectangles), critical current density J c (red circles), and characteristic voltage I c R n product (blue triangles) of the junction with different barrier thickness for absolute JJs whose resistances drop to exactly zero at 2 K.

Supplementary Note 7. Characterizations of single critical current due to the high damping environment in our Josephson junctions
The Josephson energy profile of a φ-JJ as a double-well potential should give the ground states phases of ±φ. Upon applying a bias current to the JJ corresponding to a tilt of the potential, the phase could be trapped either in the +φ or -φ state following a specific bias sweeping sequence 28,29 . When switching the JJ from the zero-voltage state to the voltage-state, the required tilt for the escape of the phase out of the potential well depends on the initial state (+φ or -φ) of the junction, at a higher I c+ or a lower I c-alternatingly as we could get in positive bias currents. Goldobin et al [29][30][31] have demonstrated the possibility to use I c± for the readout of the ±φ state. However, the I c-can only be observed together with the I c+ for low damping 30,32 (underdamped JJ) allowing a stationary phase motion, because the probability to find the phase trapped in a certain state is nearly constant for ±φ. In contrast, for the regime of high damping (overdamped JJ), the retrapping process is deterministic to predict the particular destination well of the Josephson phase, therefore only the higher I c+ can be measured. Due to non-negligible electronic noises in the external circuits of our instruments, as we have suggested in the main text, we were not able to evaluate the damping in our JJs accurately 33,34 , since their IVCs show tiny hystereses in the process of switching and retrapping like high damping JJs. For device #02 in which the unconventional Fraunhofer pattern was observed, we used two schemes to measure the two switching currents at zero magnetic field as predicted in a φ-JJ. The initial state of the JJ is prepared by sweeping from a large positive bias to zero. According to the comparison in Supplementary Fig. 17a, a switching current I c0 got by sweeping from zero to positive bias is almost identical to the other one I c-, got by sweeping from large negative bias to positive bias. In other words, we failed to observe the evidence of a φ-JJ following this method because of the high damping circumstance.
Supplementary Fig. 17. Sweeping schemes dependence of switching currents of device #02 at zero magnetic field. (a) A switching current I c0 is obtained by sweeping from large positive bias to zero bias and then back to positive bias (added 1.5 mV offset for clarity), while another switching current I c-is obtained by sweeping from large negative bias to positive currents. The comparison of these two currents labeled by the black dashed line gives the almost identical value, representing the failure in determining the φ-JJ due to the high damping circumstance. where m and n are integers.

Supplementary
As the magnetic homogeneity of the interlayer is satisfied, the magnetic flux Φ through the junction is the sum of the flux Φ M due to the ferromagnetic interlayer's magnetization M and the external field flux Φ H = M + H = 4 F + H where W is the junction width perpendicular to the magnetic field, d F is the thickness of the ferromagnetic barrier and Φ H is proportional to the magnetic field. When we get the R(B) curves as shown in the upper panel in Fig. 3a, we first define the central lobe as the zero-order index (n = 0), and the successive peaks and valleys are determined as the positions of the index equal to ±1, ±3/2, ±2…… Using these relations we can transform R(B) into Φ(B) dependence as the result presented in the lower panel of Fig. 3a.
The Φ(B) curves can further be transformed into M(B) curves by subtracting the external field flux term and the results are compared in Supplementary Fig. 18a. The magnitudes of the interlayer's remanences and the saturation magnetizations can be read from the M(B) curves. However due to the nonuniformity of periods of this device especially at high fields, rough approximations of them are expected rather than the exact estimations based on our data, therefore the direct derivations as the Φ(B) curves are preferred in the main text. Supplementary Fig. 18. Extracted Φ(B) and M(B) from R(B)