Current-Phase Relation of Ballistic Graphene Josephson Junctions.

The current-phase relation (CPR) of a Josephson junction (JJ) determines how the supercurrent evolves with the superconducting phase difference across the junction. Knowledge of the CPR is essential in order to understand the response of a JJ to various external parameters. Despite the rising interest in ultraclean encapsulated graphene JJs, the CPR of such junctions remains unknown. Here, we use a fully gate-tunable graphene superconducting quantum intereference device (SQUID) to determine the CPR of ballistic graphene JJs. Each of the two JJs in the SQUID is made with graphene encapsulated in hexagonal boron nitride. By independently controlling the critical current of the JJs, we can operate the SQUID either in a symmetric or asymmetric configuration. The highly asymmetric SQUID allows us to phase-bias one of the JJs and thereby directly obtain its CPR. The CPR is found to be skewed, deviating significantly from a sinusoidal form. The skewness can be tuned with the gate voltage and oscillates in antiphase with Fabry-Pérot resistance oscillations of the ballistic graphene cavity. We compare our experiments with tight-binding calculations that include realistic graphene-superconductor interfaces and find a good qualitative agreement.

The past few years have seen remarkable progress in the study of graphene-superconductor hybrids. This surge in interest has primarily been driven by the ability to combine high-quality graphene with superconductors via clean interfaces and has led to several experimental advances. These include the observation of specular Andreev reflection,[1] crossed Andreev reflections,[2] and superconducting proximity effects in ballistic graphene Josephson junctions (JJs).[3−7] In a majority of these studies the device comprises of graphene encapsulated in hexagonal boron nitride (BN) contacted along the edge by a superconductor. The encapsulation in BN keeps the graphene clean, while the edge contacting scheme provides transparent interfaces. In particular, ballistic JJs fabricated in this manner have been central to recent studies of novel Andreev bound states in perpendicular magnetic fields,[4] edge-mode superconductivity,[5] and supercurrents in the quantum Hall regime.[6] However, to date there have been no measurements of the Josephson current phase relation (CPR) in these systems.

The CPR is arguably one of the most basic properties of a JJ and provides information about the Andreev bound state (ABS) spectrum in the junction. While typical superconductor–insulator–superconductor (SIS) JJs exhibit a sinusoidal CPR, deviations from this behavior can be present in superconductor–normal–superconductor (SNS) junctions. Examples of this include JJs with high transmission such as nanowires[8,9] and atomic point contacts,[10,11] where the CPR contains significant higher frequency components. Furthermore, the periodicity of the CPR itself can be different from 2π for more exotic systems such as topological JJs.[12] For graphene JJs, there have been several numerical estimates of the CPR which take into account its linear dispersion relation.[13−17] More recently, ballistic graphene JJs operated in large magnetic fields have been predicted to undergo a topological transition,[18] which should be detectable via direct CPR measurements. However, the experimental determination of the CPR in graphene has been restricted to junctions that are either in the diffusive limit[19] or in a geometry that does not allow gate control of the junction properties.[20]

Here, we use a direct current (dc) superconducting quantum interference device (SQUID) to directly determine the CPR in encapsulated graphene JJs. These graphene SQUIDs stand out from previous studies[21,22] in two important ways. First, the superconducting contacts are made with molybdenum rhenium (MoRe), which allows us to operate the SQUID up to 4.2 K. More importantly, our SQUID consists of graphene JJs that are ballistic and independently tunable, thereby allowing full electrostatic control over the SQUID response. By applying appropriate gate voltages we can continuously tune from a symmetric to an asymmetric SQUID. We show that the asymmetric configuration allows us to directly extract the CPR from flux periodic oscillations in the critical current of the SQUID. The CPR is found to be nonsinusoidal, displaying a prominent forward skewing. This skewness can be tuned over a large range with the gate voltage and shows correlations with Fabry-Pérot (FP) resistance oscillations in the ballistic cavity. We complement our experiments with tight-binding simulations that go beyond the short junction limit and explicitly take into account realistic graphene–superconductor interfaces.

Figure a shows a scanning electron micrograph and cross-sectional schematic of a device. It consists of two encapsulated graphene JJs contacted with MoRe, incorporated in a SQUID loop. The fabrication strategy is similar to earlier work[3] and further details are provided in the Supporting Information (SI). The left (L-JJ)/right (R-JJ) JJs can be tuned independently by applying voltages (VL/VR) to local top gates. The junctions are intentionally designed to have a geometrical asymmetry, which ensures that the critical current of R-JJ (IcR) is larger than that of L-JJ (IcL) at the same carrier density. We report on two devices (Dev1 and Dev2) both of which have the same lithographic dimensions (L × W) for L-JJ (400 nm × 2 μm). The dimensions of R-JJ for Dev1 and Dev2 are 400 nm × 4 μm and 400 nm × 8 μm, respectively. All measurements were performed using a dc current bias applied across the SQUID in a dilution refrigerator with a base temperature of 40 mK.

Figure 1

(a) Scanning electron micrograph of the graphene dc-SQUID (Dev1) along with a cross-sectional schematic. Gate voltages VL and VR independently control the carrier density of the left and right junction, respectively. (b) Resistance R across the SQUID versus VL and VR, demonstrating independent control of carrier type and density in the JJs. The bias current (I) for these measurements was fixed at 500 nA. (c) Line trace taken along the dashed white line in (b) showing Fabry-Pérot oscillations in the hole-doped regime. (d) Differential resistance dV/dI as a function of dc current bias I and magnetic field B with the SQUID operated in a symmetric configuration (VL = +10 V and VR = +2.5 V). Flux-periodic oscillations are clearly visible with a slowly decaying envelope arising from the interference pattern of a single JJ. (e) V–I plots [extracted from (d)] for different values of magnetic flux Φ showing a nearly 100% modulation of the critical current. All measurements shown here are performed at T = 4.2 K.

Figure b shows the variation in the normal state resistance (R) of the SQUID with VL and VR at T = 4.2 K. The device was biased with a relatively large current of 500 nA, which is larger than the critical current of the SQUID for most of the gate range. Figure c shows a single trace taken along the white dashed line of Figure b, where R-JJ is held at the charge neutrality point (CNP). We see clear FP oscillations on the hole (p) doped region due to the formation of n–p junctions at the superconductor–graphene interfaces.[3,4] Furthermore, the criss-cross pattern seen in the lower left quadrant of Figure b indicates that both graphene junctions are in the ballistic limit and that there is no cross-talk between the individual gates. We note that when VR > 3 V the critical current of the SQUID (Ic) is larger than the applied current bias, and a zero-resistance state is thus visible even at 4.2 K. Having established the fact that our JJs are in the ballistic regime, we now look in more detail at the behavior of the SQUID. At T = 4.2 K, we first tune the gate voltages (VL = +10 V, VR = +2.5 V) such that the SQUID is in a symmetric configuration and IcR = IcL. Figure d shows the variation in differential resistance dV/dI with current bias I and magnetic field B, where we observe clear oscillations in Ic with magnetic flux. In this configuration, the modulation in Ic is nearly 100%, as seen by the individual V–I traces in Figure e. The slow decay in the maximum value of Ic arises due to the (Fraunhofer) magnetic field response of a single junction. The devices were designed such that this background was negligible around B = 0 (that is, the SQUID loop area was kept much larger than the JJ area). Minimizing this background is important for a reliable determination of the CPR, as we will see below.

We now turn our attention to the flux-dependent response of a highly asymmetric SQUID (IcR ≫ IcL), a condition that can be readily achieved by tuning the gate voltages appropriately. To outline the working principle of the device, we start with the assumption that both JJs have a sinusoidal CPR (a more general treatment can be found in the SI). So, the critical current of the SQUID can be written as Ic = IcL sin θ + IcR sin δ, where θ (δ) is the phase drop across L-JJ (R-JJ). When an external magnetic flux (Φ) threads through the SQUID loop, the flux and phase are related by δ – θ = 2πΦ/Φ0, assuming the loop inductance is negligible. Now, when IcR ≫ IcL the phase difference across R-JJ is very close to π/2. Thus, Ic(Φ) ≈ IcR + IcL sin(2πΦ/Φ0 + π/2) and the flux-dependence of Ic directly represents the CPR of L-JJ, that is, Ic(Φ) ≈ IcR + Is(ϕ), where Is is the supercurrent through L-JJ and ϕ is the phase drop across it. This principle of using an asymmetric SQUID to probe the CPR has been employed in the past for systems such as point contacts[10,11] and vertical graphene JJs,[20] where an SIS junction (with a well-known sinusoidal CPR) was used as the reference junction. In our case, the reference junction is also a graphene JJ, where the CPR is not known a priori. We show (see SI) that this does not affect our ability to probe the CPR, provided time reversal symmetry is not broken, meaning that the CPR satisfies the condition Is(ϕ) = −Is(−ϕ).[23] Throughout the remainder of the text we use R-JJ as the reference junction (larger critical current), and L-JJ is the junction under study (smaller critical current).

Figure a shows the typical magnetic response of the asymmetric SQUID at T = 40 mK with VL = −4 V (fixed) and different values of VR. For the most asymmetric configuration (VR = +10 V) Ic oscillates around a fixed value of roughly 6 μA (IcR) with an amplitude of about 500 nA (IcL). Using the arguments described above, this Ic(Φ) curve can be converted to Is(ϕ), as shown in Figure b. Here Is is the normalized supercurrent defined as (Ic – IcR)/IcL. We note that there is an uncertainty (less than one period) in the exact position of zero B. This, combined with the unknown CPR of the reference graphene JJ, makes it important to do this conversion carefully, and we describe the details in the SI. The CPR shows a clear deviation from a sinusoidal form, showing a prominent forward skewing (that is, Is peaks at ϕ > π/2). We define the skewness of the CPR as S = (2ϕmax/π) – 1,[19] where ϕmax is the phase for which the supercurrent is maximum.

Figure 2

(a) Variation of Ic with Φ for VL = −4 V and VR = +10, +5, and +3 V at 40 mK. Solid black lines are results from RCSJ simulations of the SQUID. (b) Variation of supercurrent Is = (Ic – IcR)/IcL with phase ϕ extracted from the top curve in (a). ϕmax indicates the phase at which Is reaches a maximum and is noticeably different from π/2, indicating a forward skewed CPR.

To be certain that we are indeed measuring the CPR of L-JJ, we perform some important checks. We keep IcL fixed and reduce IcR (by reducing VR). Figure a shows that reducing IcR merely shifts the Ic(Φ) downward and therefore does not affect the extracted CPR, as one would expect. Furthermore, we use the experimentally determined CPR (from Figure b), the junction asymmetry, and loop inductance as inputs for the resistively and capacitively shunted junction (RCSJ) model to compute the expected SQUID response (see SI for details of the simulations). These plots (solid lines) show an excellent agreement between simulations and experiment, thus confirming that the asymmetry of our SQUID is sufficient to reliably estimate the CPR of L-JJ. Furthermore, it shows that there are no significant effects of inductance in our measurements, which could potentially complicate the extraction of the CPR from Ic(Φ) in an asymmetric SQUID.[24]

To study the gate dependence of the CPR we fix VR at +10 V (to maximize IcR) and study the change in S with VL (Figure a) for Dev1 and Dev2. The SI shows the Is(ϕ) curves used to extract S. For both devices we find that S is larger on the n-side as compared to the p-side, showing a dip close to the CNP. We note that Dev2 allows us to probe the CPR up to a larger range on the n-side due to its larger geometric asymmetry (see SI for other measurements on Dev2). We expect the skewness to depend strongly on the total transmission through the graphene JJ, which should depend on (a) the number of conducting channels in the graphene, as well as (b) the transparency of the graphene-superconductor interface. The gate voltage VL obviously changes the Fermi wave vector but it also changes the contact resistance,[25] which plays a significant role in determining S. This can be seen most clearly for Dev2 for high n-doping, where S saturates, despite the fact that IcL continues to increase up to the largest measured density (see inset). At large p-doping, S also seems to saturate but a closer look (Figure b) shows that S oscillates in antiphase with the FP oscillations in resistance. This clearly indicates that in this regime the CPR is modulated by phase coherent interference effects similar to the FP oscillations that affect the total transmission.

Figure 3

(a) Variation of skewness S as a function of carrier density n for Dev1 and Dev2. The larger geometric asymmetry of Dev2 (see text) allows one to reliably probe the CPR up to higher n-doping. Inset shows the variation of IcL with density. (b) A finer scan for Dev2 shows that S oscillates with carrier density in the p-doped regime in antiphase with Fabry–Pérot oscillations in the resistance.

We complement our measurements with a minimal theoretical model by solving the corresponding Bogoliubov–de Gennes (BdG) equations to calculate the CPR in graphene JJs. To set the stage, we note that SNS junctions can be characterized by the quasiparticle mean free path lf in the normal (N) region and the coherence length ξ0 = ℏvF/Δ, where vF is the Fermi velocity in N. In our devices L ≪ lf, that is, they are in the ballistic regime, and therefore we neglect impurity scattering in our calculations. Taking vF ≈ 106 m/s for graphene and Δ ≈ 1.2 meV for MoRe, one finds ξ0 = 548 nm, which means that in our junctions L ≲ ξ0, that is, they are not in the strict short junction limit L ≪ ξ0. Consequently, the Josephson current is carried not only by discrete Andreev bound states (ABSs) but also by states in the continuum (CONT).[26−28] For this reason, we numerically solve the BdG equations using a tight-binding (TB) model (see Figure a) and a recently developed numerical approach[17,29] which handles the ABSs and states in the continuum on equal footing. The description of both the normal region and the superconducting terminals is based on the nearest-neighbor TB model of graphene.[30] The on-site complex pair-potential Δ is finite only in the superconducting terminals and changes as a step-function at the N–S interface. The results presented here are calculated using the top-contact geometry (Figure a), a model with perfect edge contacts is discussed in the SI. As observed experimentally, we take into account n-doping from the superconducting contacts (see Figure a inset). If the junction is gated into hole-doping, a FP cavity is formed by the two n–p junctions in the vicinity of the left and right superconducting terminals. The length L* of this FP cavity depends on the gate voltage,[4] for stronger hole-doping the n–p junctions shift closer to the contacts. For further details of the model see SI.

Figure 4

(a) The geometry of the system used in the calculations. The superconducting leads are attached in a top-contact geometry to the normal graphene sheet and overlap with the normal graphene sheet over NL unit cells. γ denotes the nearest-neighbor intralayer hopping in the leads and in the graphene sheet, while γ1 is the nearest-neighbor interlayer hopping. A periodic boundary condition is applied in the y-direction. (Inset) Top view of the system. Because of doping from the S contacts, the normal graphene region is assumed to be n-doped up to a distance x1 (x2) from the left (right) contact. The distance L* = x2 – x1 is the effective cavity length which depends on the gate voltage applied to the junction. (b) The contribution of the ABSs (red) and continuum CONT (blue) to the total supercurrent (black) as a function of the phase difference for an n-doped junction (n = 0.9 × 1011 cm–2 and L/ξ0 = 0.73). (c) The skewness S as a function of doping of the junction. The regimes i–iii indicated by the rectangles are further discussed in the text. Dashed lines show the average S in the p- and n-doped regime. (d) The skewness (red circles, left axis) and normal state resistance (blue, right axis) versus doping for strong p-doping of the junction.

Turning now to the CPR calculations, Figure b shows separately the contribution of the ABS and the continuum to the supercurrent. Because L ≲ ξ0, the latter contribution is not negligible and affects both the value of the critical current and the skewness of the CPR. We note that our calculations can qualitatively reproduce the doping dependence of Ic (see SI), however the obtained values of Ic are about 2.5 times larger than the measurements (Figure a inset). The exact reason for this discrepancy is not known, but a similar disagreement between theory and experiment was also found in ref (4). Focusing now on the skewness, in Figure c we show the calculated skewness S as a function of the doping of the junction at zero temperature. We consider three regimes: (i) strongly p-doped junction; (ii) large n-doping, (iii) the region around the CNP. We start with the discussion of (i). It is well established that in this case the p–n junctions lead to FP oscillations in the normal resistance as well as in the critical current[3,4] of graphene JJs. Our calculations, shown in Figure d, indicate that due to FP interference the skewness also displays oscillations as a function of doping around an average value of S ≈ 0.23. As already mentioned, similar oscillations are present in the normal state resistance R, however, we find that R oscillates in antiphase with the skewness. Compared to the measurements (Figure b), our calculations therefore reproduce the phase relation between the oscillations of the skewness and R and give a qualitatively good agreement with the measured values of the skewness. In the strong n-doped regime (case (ii) the calculated average skewness of S = 0.27 is larger than for p-doped junctions, and very close to the measured values. Small oscillations of S are still present in our results and they are due to the n–n′ interfaces, that is, the difference in the doping close to the contacts (for x < x1 and x > x2, see Figure a inset) and the junction region (x1 < x < x2), which enhances backscattering. Our calculations therefore predict smaller skewness for p-doped than for n-doped junctions. The enhancement of S in the n-doped regime can be clearly seen in the measurements of Figure a. We note that previous theoretical work,[16] which calculated the spatial dependence of the pairing amplitude self-consistently, predicted a skewness of S ≈ 0.15 for n-doped samples with L < ξ0, while a nonself-consistent calculation which took into account only the contribution of the ABS yielded S ≈ 0.42.[16] The comparison of these results to ours and to the measurements suggest that the skewness may depend quite sensitively on the S–N interface as well as on the L/ξ0 ratio and that our approach captures the most important effects in these junctions. Finally, we briefly discuss the case (iii), where the measurements show a suppression of the skewness as the CNP is approached. The measured values of S ∼ 0.1 are similar to those found in diffusive junctions[19] but are significantly lower than the theoretical prediction of S = 0.26 in the short junction limit[13] at the CNP. This suppression of S is not reproduced in our calculations, instead, we find rapid oscillations as the CNP is approached from the p-doped regime. This discrepancy is likely to be due to effects that are not included in our ballistic model, such as charged scatterers that are poorly screened in this regime, or scattering at the edges, which is more relevant at low densities when only a few open channels are present.

Finally, we study the effect of temperature on the CPR of these JJs. In Figure a, we compare the CPR in the n-doped regime (VL = +1 V; n = 0.9 × 1011 cm–2) at 40 mK and 4.2 K. One clearly sees that at 4.2 K the CPR is sinusoidal. This is consistent with our observation that the critical current modulation of the SQUID is nearly 100% at 4.2 K (Figure d), a condition which can only be achieved if the CPR is sinusoidal. Figure b shows the full temperature dependence of S for two representative values of electron and hole doping (see SI for the corresponding temperature dependence of the critical currents). The reduction in skewness with temperature is consistent with the thermal population of excited Andreev bound states and continuum states.[14,16,17,19] We compare our results with numerical estimates which take into account this effect, along with the temperature dependence of the pair potential (black curves in Figure b). While the qualitative behavior of both (experimental and numerical) curves are similar, the experimentally determined skewness reaches zero (sinusoidal CPR) faster than the numerics. At this point, it is difficult to ascertain the exact reason for this discrepancy but one possible explanation for this is that the induced superconducting gap in the graphene is somewhat smaller than the bulk MoRe gap, resulting in a faster decay.

Figure 5

(a) CPR for VL = +1 V (n = 0.9 × 1011 cm–2) at 40 mK (upper curve) and 4.2 K (lower curve). Solid line shows the calculated CPR. A forward skewness is clearly seen in the curve at 40 mK but is absent at 4.2 K. (b) Variation of S with temperature for electron and hole doping. Increasing the temperature suppresses higher harmonics in the CPR, thereby reducing S until it vanishes near 4.2 K and the curves become sinusoidal. Black lines show the results of tight binding simulations.

In conclusion, we have used a fully gate-tunable graphene based SQUID to provide measurements of the current-phase relation in ballistic Josephson junctions made with encapsulated graphene. We show that the CPR is nonsinusoidal and can be controlled by a gate voltage. We complement our experiments with tight binding simulations to show that the junction length and nature of the superconductor–graphene interface play an important role in determining the CPR. We believe that the simplicity of our device architecture and measurement scheme should make it possible to use such devices for studies of the CPR in topologically nontrivial graphene Josephson junctions.