Highly skewed current–phase relation in superconductor–topological insulator–superconductor Josephson junctions

Three-dimensional topological insulators (TIs) in proximity with superconductors are expected to exhibit exotic phenomena, such as topological superconductivity (TSC) and Majorana-bound states (MBS), which may have applications in topological quantum computation. In superconductor–TI–superconductor Josephson junctions, the supercurrent versus the phase difference between the superconductors, referred to as the current–phase relation (CPR), reveals important information including the nature of the superconducting transport. Here, we study the induced superconductivity in gate-tunable Josephson junctions (JJs) made from topological insulator BiSbTeSe2 with superconducting Nb electrodes. We observe highly skewed (non-sinusoidal) CPR in these junctions. The critical current, or the magnitude of the CPR, increases with decreasing temperature down to the lowest accessible temperature (T ~ 20 mK), revealing the existence of low-energy modes in our junctions. The gate dependence shows that close to the Dirac point the CPR becomes less skewed, indicating the transport is more diffusive, most likely due to the presence of electron/hole puddles and charge inhomogeneity. Our experiments provide strong evidence that superconductivity is induced in the highly ballistic topological surface states (TSS) in our gate-tunable TI-based JJs. Furthermore, the measured CPR is in good agreement with the prediction of a model which calculates the phase-dependent eigenstate energies in our system, considering the finite width of the electrodes, as well as the TSS wave functions extending over the entire circumference of the TI.


INTRODUCTION
Three-dimensional (3D) topological insulators (TIs) are a new class of quantum matters and are characterized by an insulating bulk and conducting topological surface states (TSS). These TSS are spin-helical with linear Dirac fermion-like energy-momentum dispersion. 1,2 The TSS of 3D TIs in proximity to s-wave superconductors are among the top candidates proposed to realize topological superconductors, 3 capable to support Majoranabound states (MBS) and promising for future applications in topological quantum computing. 4,5 A Josephson junction (JJ) made of a TI with two superconducting contacts is one of the most common platforms to study the nature of the induced superconductivity in TIs and possible topological superconductivity. One of the fundamental properties of a JJ is its supercurrent (I) as a function of the phase (φ) difference between the two superconducting contacts, referred to as the current-phase relation (CPR), where the maximum of I(φ) is the critical current (I C ) of the JJ. Given the topological protection or the prohibited backscattering from non-magnetic impurities in the TSS of 3D TIs, 1,2 superconductor-TI-superconductor (S-TI-S) junctions are expected to demonstrate novel features in their CPR. [6][7][8][9] While for conventional junctions the CPR is 2π-periodic, for TI-based JJs the CPR is predicted to have an additional 4π-periodic component. 3,10 This 4π-periodicity originates from the zero-energy crossing (at φ = π) of the Andreev bound states (ABS) and is protected by the fermion parity conservation. However, if the temporal variation of φ is slower than the quasiparticle poisoning time, the 2π-periodicity of the CPR is restored, which can mask the unique topological nature of the JJs. [10][11][12] Nonetheless, in this case the topologically protected modes can give rise to a highly nonsinusoidal ("skewed") 2π-periodic CPR similar to a perfectly ballistic (scattering free) JJ. 11,13,14 JJs have been experimentally studied in 3D TIs [15][16][17][18][19][20][21][22][23][24][25][26] including TI nanoribbons (TINRs). 27,28 These earlier studies of TI-based JJs have reported indirect signatures of skewed CPR using a variety of techniques, such as temperature-dependent I C and Fraunhofer patterns. 17,20,25 However, it has been challenging to observe significant skewness in direct phase-sensitive measurements of CPR. 18,28 In this work, we fabricate S-TI-S junctions based on the topological insulator BiSbTeSe 2 flakes, which have an insulating bulk and demonstrate TSS-dominated electrical properties at low temperatures. 29,30 We measure the CPR in the S-TI-S junctions using an asymmetric superconducting quantum interference device (SQUID). 31 Remarkably, the measured CPR in our S-TI-S junctions are highly skewed, revealing that the superconducting transport is carried by the ballistic TSS in our TI JJs. Furthermore, we observe that the skewness depends on the back-gate voltage (V g ) and is the smallest close to the charge neutrality point (CNP). We present a theoretical model based on the induced superconductivity in the ballistic TSS of the TI. This model takes into account the finite-size (of both Nb and TI) and proximity effects, and relates the induced supercurrent to the TSS that extend over the entire circumference of the TI. The calculated energy spectrum (energy vs. phase φ) of the junction reveals the existence of extremely low-energy modes that exist over the entire range of phases, i.e. 0 ≤ φ < 2π. The computed CPR from the theory is in good agreement with the experimental results.

RESULTS AND DISCUSSION
We adapt an asymmetric SQUID technique 31 to measure the CPR in our TI (BiSbTeSe 2 )-based JJ. Our BiSbTeSe 2 crystals are among the most bulk-insulating 3D TIs, where the Fermi energy lies within the bulk bandgap and inside the TSS, as verified by the angle-resolved photoemission spectroscopy (ARPES) and transport measurements. 29 Exfoliated thin films of this material exhibit ambipolar field effect, as well as several signatures of topological transport through the spin-helical Dirac fermion TSS, including the half-integer quantum Hall effect and π Berry phase. 29,30 Furthermore, we have recently observed an anomalous enhancement of the critical current in JJs based on BiSbTeSe 2 nanoribbons, demonstrating the induced superconductivity in the TSS. 28 Figure 1a shows a scanning electron microscope (SEM) image of an asymmetric SQUID with a BiSbTeSe 2 flake (sample A). The asymmetric SQUID consists of two JJs in parallel. The first JJ is the S-TI-S junction with an unknown CPR, I(φ), and is highlighted by the dashed white rectangle in Fig. 1a. The second JJ is a conventional S-S′-S junction (REF junction), where S and S′ are 300 and 80 nm wide Nb lines, respectively. The data presented here comes from two devices, sample A (width W~2 μm and thickness t~40 nm) and sample B (W~4 μm and t~13 nm). All our measurements are performed in a dilution refrigerator with a base temperature (T) of~20 mK.
Supercurrent (I REF ) is the superconducting magnetic flux quantum, and L s is the self-inductance of the SQUID loop. We can estimate L s~5 pH 32 for our SQUID assuming a constant current density, Nb London penetration depth λ~39 nm 33 and Nb film thickness t~40 nm. Although penetration depth may not be applied easily in our system, the calculated L s is acceptable within a factor of 2. Since L s I/Φ 0~0 .04 ≪ 1, we can ignore the contribution of the selfinductance L s in the phase difference, thus φ − φ R = 2πΦ B /Φ 0 and . In general, we can reconstruct the CPR of the TI-based junction using an analytical approach presented in ref. 34 However, in our asymmetric SQUID, the REF junction is designed such that I C REF ≫ I C , where I C is the critical current of the S-TI-S junction. Thus φ R ≈ π/2 when the SQUID reaches its critical current with I C SQUID~I C REF + I(2πΦ B /Φ 0 + π/2). Therefore, the modulation of the I C SQUID vs. B = Φ B /S (in period of B 0 = Φ 0 /S) will provide a very good approximation of the CPR, particularly the skewness of I(φ) in the TI-based JJ. Figure 1b depicts the set-up for the measurement of the CPR in our TI-based JJs. In order to reduce the uncertainty of the measured I C SQUID due to thermal and quantum fluctuations, we use a square wave pulsed current (frequency f~17 Hz) with 50% duty cycle to bias the SQUID. The voltage (V S ) across the SQUID is measured by a lock-in amplifier and the magnetic flux in the SQUID is varied by an external magnet. We also performed dc measurements, which resulted in qualitatively similar behavior to the pulsed measurements. However, the dc mode is more sensitive to thermal fluctuations and the Joule heating effect. Therefore, the CPR is noisier and its amplitude (i.e., the critical current I C ) is smaller compared with the pulsed measurement.
For a fixed Φ B , once the amplitude of the pulsed current is increased above I C SQUID , a non-zero voltage appears across the SQUID. Figure 1c depicts a color map of V S as functions of I SQUID and the external magnetic field (B) applied out of the plane of the SQUID. The solid white curve highlights I C SQUID vs B. We estimate I C REF~1 8 μA (the dashed red line in Fig. 1c) by taking average of Figure 1d shows the supercurrent I(φ) normalized by its amplitude (I C ) vs. φ measured in sample A for V g = 0 V at T = 20 mK (red symbols). The measured CPR does not exhibit any hysteresis vs. magnetic flux. However, since the absolute value of the flux in the SQUID is unknown (for instance, due to a remnant field), we shift the experimental curve in the horizontal direction such that φ = 2πΦ B /Φ 0 . The measured CPR is contrasted with a reference sin(φ) shown by the dashed blue curve in Fig.1d. The measured CPR in sample A is forward skewed, i.e. its maximum occurs at φ = 0.75π (instead of π/2 for sin(φ)).
It is also noted that the measured CPR in Fig. 1d (at T = 20 mK) is non-anti-symmetric (NAS) around zero, i.e., not symmetric under simultaneous sign changes in Φ B and I (note for a 2π periodic function such as the CPR, being antisymmetric around zero is equivalent to also being antisymmetric around π). A similar NAS current vs. flux curve can also be seen in the theoretical and experimental results of ref. 31 for atomic contacts. Although those authors do not explicitly discuss this NAS, they do discuss the difference between the observed current-flux relation (which is taken as a measure of the CPR) and a theoretically predicted CPR (based on conventional models which ignores stochastic switching processes) with an explanation that we find plausible in our case as well. In particular, like ref., 31 we observe NAS only for the highly skewed CPR and not for the more sinusoidal ones. This NAS is an artifact of the assumption that switching to the normal state occurs when the phase φ R across the REF junction is exactly π/2, used in reconstructing the CPR from experimental data. In practice, the superconducting to normal transition of the SQUID is governed by the stochastic switching processes, as well as thermal fluctuations, and the tilted washboard potential used to take these processes into account 31 is not symmetric under a sign change of the applied flux. As a result, even if the underlying CPR of the TI junction in our case is anti-symmetric (AS) around zero, i.e. symmetric under simultaneous sign change in φ and I, the one reconstructed from the experiment may be NAS. We further note that while the stochastic processes could smear the "real CPR" and introduce artifacts, such as NAS (not present in the real CPR) in the measured CPR, such NAS observed in the experimental CPR does not mask the qualitative fact that the underlying real CPR of the junction is skewed. In other words, as found in ref. 31 for junctions with the sinusoidal CPR, the current vs. flux curves are AS. However, for junctions with non-sinusoidal CPR, the current vs. flux curve becomes NAS. The above argument is further supported by the observation that when we use different TI flakes or use gate tuning (in either case the reference junction and SQUID arm are not affected), we can obtain less skewed or non-skewed (sinusoidal) CPR and at the same time the NAS goes away. Moreover, a similar measurement (on a TI nanoribbon JJ) using an asymmetric SQUID with the same geometry as the one used here resulted in a sinusoidal CPR. 28 Therefore, we conclude the skewness of the CPR observed in sample A originates from the TI flake and not the REF junction or the stochastic switching process. Figure 2a depicts the normalized supercurrent I/I C vs. φ measured at a few different temperatures in sample A. The amplitude of the CPR (i.e. I C ) as a function of T is plotted in Fig. 2b. We observe that the CPR remains highly non-sinusoidal up to T4 00 mK, but becomes nearly sinusoidal at higher T = 1.3 K. Furthermore, I C exhibits a strong T dependence and increases as we decrease the temperature down to the lowest accessible T = 20 mK. Such a temperature dependence is in contrast to that of conventional junctions, where I C is expected to saturate at low temperatures. 35 Figure 2c depicts the amplitude of the fast Fourier transform (FFT) normalized by the amplitude of the first harmonic as a function of 2π/φ = B 0 /B, where B 0 = Φ 0 /S~1.1 G and S1 6 μm 2 is the area of the SQUID. The FFT is calculated for the data taken at T = 20 mK in the range −10π ≤ φ ≤ 10π and reveals that the CPR can be described by a Fourier series with up to six harmonics. The blue and black curves are predictions of a general model for ballistic junction and our model for TI junction, respectively, and will be discussed later. In order to describe the shape of the CPR in our samples, we define the total harmonic distortion (THD) as where A j is the amplitude of the j th harmonic. Figure 2d depicts THD, A 2 /A 1 , and A 3 /A 1 vs. T in sample A at V g = 0 V. We observe that THD, A 2 /A 1 , and A 3 /A 1 are nearly temperature independent up to T~400 mK. Moreover, at T = 1.3 K, A 3 /A 1~0 and THD~A 2 /A 1 , indicating that at this temperature, only the first and second harmonics are present in the CPR. Thus, the CPR of the TI junction is less skewed compared to that at the base temperature. Figure 3a demonstrates the CPR measured at different V g 's for sample B at T = 30 mK. The inset of Fig. 3b depicts the twoterminal resistance R of the SQUID vs. V g measured at T = 8 K, above the critical temperature (T C Nb~7 K) of the Nb electrodes. Sample B exhibits a strong gate dependence and an ambipolar field-effect in its normal-state resistance with the CNP at V CNP− 15 V. We also observe that in sample B the skewness changes as a function of V g . Figure 3b plots the THD vs. V g for both sample A (red) and sample B (blue). We note that the CPR is most skewed in sample B at V g = 30 V, where the chemical potential is inside the bulk bandgap yet away from the CNP (see the inset of Fig. 3b). The reduced skewness at V g~0 V may be a result of the charge inhomogeneity and electron/hole puddles near the CNP.

THEORETICAL MODELING
In this section, we introduce a theoretical model based on the induced superconductivity in the spin-helical surface states of TIs. Our model builds upon the existing model proposed by Fu and Kane 3 and adapts it to systems of the type studied in the present experiment. Specifically, it considers the finite size of the superconducting contacts as well as the TI flakes. Since the superconducting (Nb) contacts in our case are only 300 nm wide (a value comparable to the expected coherence length ξ = ħv F /Δ 0~3 30 nm of the junction), we cannot assume existence of the ABS (confined inside the junction) but should suppose instead that the surface state wavefunction extends over the entire circumference of the sample (see Fig. 4a). We denote the circumference C x and define the longitudinal coordinate x to be in the range −C x /2 ≤ x ≤ C x /2. We adopt the Hamiltonian of Fu and Kane 3 and take the pairing amplitude to be a piecewise constant function of x as follows: Δ(x) = Δ 0 exp(iφ/2) for L/2 < x < L/2 + b, Δ 0 exp(−iφ/2) for −L/2−b < x < −L/2, and zero otherwise. Here L and b are the separation and width of the contacts, respectively. The wavefunction is subject to antiperiodic boundary conditions in x. 36 In this simple model, we assume that the system is translationally invariant in the y direction, so the wavefunction depends on y as exp(ik y y) for some k y . This renders the problem effectively onedimensional.
Our explanation of the CPR and temperature dependence of the critical current is based on an interplay between the finite-size and proximity effects. To compute the CPR, we first rewrite, following ref., 3    eigenvalue problem for fHg into a matrix problem, which we diagonalize numerically for various values of the phase difference φ. fHg has a particle-hole symmetry, which stems from using four fermionic components in place of two: at each φ, the energy levels come in ±E pairs. In terms of the nonnegative levels E n ≥ 0 (one from each pair), the total free energy at finite temperature T is 37,38 and the current is obtained as IðφÞ ¼ 2e h À Á dF dφ . As we increase the number of k x (or j) participating in the expansion, F(φ) suffers from an ultraviolet divergence, but the current does not. To calculate finite temperature properties, we replace Δ 0 above with the Tdependent superconducting gap Δ(T) modeled using the BCS selfconsistent equation. 39 The energy spectrum (±E n vs. φ) for sample A for the modes with k y = 0 and energies within the gap, |E n | ≤ Δ 0 is shown in Fig. 4b. Interestingly, we observe modes with energies much smaller than Δ 0 that extend over the entire range of φ, see red curves in Fig. 4b. These low-energy states lead to the non-saturation of the junction's critical current down to our lowest accessible temperature (T2 0 mK) as seen in Fig. 2b in the theoretical (blue) curve, consistent with the experimental data (symbols).
Because the wavefunction extends over the entire circumference C x , while the Nb contacts occupy only a small part of it, the energy scale of the low-energy modes is only a fraction of the full Δ 0 . Our results can be understood qualitatively using the perturbation theory. For Δ 0 = 0, the energies are 3 so there is a strictly zero energy state whenever k′, defined above, equals one of the quantized freefermion momenta where n ≥ 0 is an integer. When Δ 0 > 0, these states are gapped roughly by 2bΔ 0 /C x . For sample A with C x~6 μm and the contact width b~300 nm, this is about 0.1Δ 0 . Crucially, these low-energy states exist for the entire range of phases, 0 ≤ φ ≤ 2π, in contrast for instance to the case of a conventional ballistic junction (the Kulik-Omelyanchuk theory 40 ), where the minimal excitation energy remains on the order Δ 0 except for a narrow vicinity of φ = π.
For a given μ, the condition (3) will be satisfied better for some k y than for others. In practice, the translational invariance in the y direction is not precise, so k y is not an exact quantum number. Nevertheless, because of the large size of the TI flake in the transverse (y) direction, the quantization interval for k y is small, so unless μ is exceptionally close to zero, we expect there will be a significant number of modes for which Eq. (3) is satisfied to a good accuracy. We therefore adopt the simple model in which we calculate the supercurrent for k y = 0 only and multiply the result by an effective number of transverse channels N ch to account for the contribution of all the modes. We determine N ch by matching the overall magnitudes of the experimental and computed critical currents. We find N ch~1 9 and N ch~4 6 for sample A at V g = 0 V and sample B at V g = 30 V, respectively. The role of the chemical potential in this effectively one-dimensional model is played by μ' = ħv F k', which is now considered as a parameter. It is distinct from the true chemical potential μ, obtained from the gate voltage , where e is the electron charge and C g = 12 nF cm −2 is the parallel plate capacitance per unit area of a 300-nm SiO 2 .
We plot the computed CPR for sample A as solid curves in Fig. 2a, where an excellent agreement with the measured data is observed. The blue curve in Fig. 2c is the FFT calculated for the theoretical CPR (in the range −10π < φ < 10π and at T = 20 mK) of a perfectly ballistic short junction: 35,41 IðϕÞ where Δ(T) is the temperature-dependent superconducting gap of the junction obtained from the BCS theory. 39 Notably, the experimentally observed A 2 /A 1 in Fig. 2c is within 3% of that predicted for the fully ballistic junction (blue curve) and the THD = 0.46 extracted from our measured CPR at T = 20 mK is within 20% of the theoretical ballistic limit (THD = 0.55), indicating sample A is nearly ballistic. The black curve in Fig. 2c plots the FFT of the CPR calculated using our theoretical model (Fig. 4). We observe that the FFT of the CPR, calculated using our model, is in reasonable agreement with the FFT of the measured CPR. In contrast, the perfectly ballistic model (blue curve) notably overestimates the higher harmonics (A 3 and above). The computed CPR for sample B is plotted with dashed curves in Fig. 3a at two different V g 's. Theoretical calculations were done with μ′ = 0 and 50 meV, respectively. While the theoretical CPR agrees well with the experiment for V g = 30 V, we see a deviation between the theory and experiment for V g = 0 V. Sample B is much thinner (t~13 nm) than sample A (t~40 nm). When a TI becomes sufficiently thin, there may be a gap opening in the TSS close to the Dirac point due to hybridization of the top and bottom surface states. This gap causes the TI to transition into a trivial insulator. Moreover, there are electron-hole puddles and charge inhomogeneity near the CNP. Therefore, the transport may be more diffusive, i.e. the CPR is more sinusoidal, close to the CNP due to effects of disorder and hybridization. Such effects are not included in our theory and may be responsible for the discrepancies between the calculated and measured CPR at V g = 0 V in sample B. Lastly, we note that even though both samples are exfoliated from the same crystal and undergone similar fabrication processes, sample to sample variations may still be present and play a role for the observed differences between the CPR's of samples A and B. In our previous experiments on S-TINR-S JJs, 28 even though we have also observed evidence that the superconductivity is induced in the TSS, we only observe a sinusoidal CPR. A possible reason for this is a much smaller transverse size (C y ) of the TINR compared to the flakes used in the current work. As a consequence, k y is quantized in larger units (i.e. 2π/C y ), and the condition (3) is less readily satisfied. Effectively, the small transverse size generates a gap in the TSS spectrum, preventing the occurrence of low-energy states and rendering the CPR more sinusoidal at our experimental temperatures. 28 A similar explanation may be relevant also for sample B of the present paper near the CNP.
We lastly note that the non-sinusoidal CPR was previously reported in ref. 18 However, the observed CPR in our TI flakes are more skewed. Moreover, our TI flakes are more disordered compared to HgTe quantum wells of ref. 18 For instance, the normal-state transport in HgTe samples shows mobilities as large as 10,000-40,000 cm 2 V −1 s −1 , 17 while our TI flakes have normalstate mobility of~1000 cm 2 V −1 s −1 . 29 Therefore, the observed skewed CPR and the underlying ballistic superconducting transport in our disordered TI flakes provide strong evidence for the topological protection of the modes with k y near zero. This is further corroborated by the fact that in our previously studied junctions based on TINRs, 28 where such "topology" is changed (since Dirac point becomes gapped, removing the k y~0 modes and backscattering is now allowed for the remaining modes with finite k y ), the observed CPR becomes sinusoidal (skewness disappears).
In conclusion, we have measured the CPR, one of the fundamental properties of a JJ, in a topological insulator BiSbTeSe 2 -based JJ using an asymmetric SQUID technique. We observed highly forward-skewed CPR, indicating that the transport through the TSS of the TI junction was close to ballistic. Temperature and gate dependence of the CPR were also studied, where we observed that CPR became more sinusoidal at high temperatures (T~1.3 K) and close to the CNP. The reduced skewness near CNP was an indication of diffusive transport and was associated with the existence of electron-hole puddles and charge inhomogeneity in the very thin TI. Moreover, we developed a theoretical model that considered induced superconductivity in the spin-helical TSS of TIs. Our model assumed that the surface states can extend over the entire circumference of the TI. The predicted skewness of the CPR and the dependence on the temperature were consistent with our experimental observations. Overall, the experiment and the theory pointed toward robust features that make our TI system an excellent candidate to observe topological superconductivity and MBS.

Sample preparation
High-quality single crystals of BiSbTeSe 2 were grown using the Bridgman technique as described elsewhere. 29 We obtain BiSbTeSe 2 flakes using the standard scotch-tape exfoliation technique and transfer them onto a 300nm-thick SiO 2 /500-μm-thick highly doped Si substrates, which are used as back gates. We then locate the BiSbTeSe 2 flakes with different width and thickness using an optical microscope. Subsequently, electron beam lithography is performed to define a SQUID consisting of the TI-based JJ and the REF junction (based on a narrower Nb line). The electrode separation, L, in the TI-based JJ is~100 nm. Finally, a thin layer (t~40 nm) of superconducting Nb is deposited in a DC sputtering system. Prior to the Nb deposition, a brief (~3 s) in situ Ar ion milling is used to clean the interface between Nb and the TI flake.