Anomalous Fraunhofer patterns in gated Josephson junctions based on the bulk-insulating topological insulator BiSbTeSe2
Subhamoy Ghatak Oliver Breunig Fan Yang Zhiwei Wang A. A. Taskin Yoichi Ando







yang@ph2.uni-koeln.de



ando@ph2.uni-koeln.de
Physics Institute II, University of Cologne, Zülpicher Str. 77, 
50937 Köln, Germany
 
(These authors contributed equally to this work.)


One-dimensional Majorana modes are predicated to form in Josephson junctions based on three-dimensional topological insulators (TIs). While observations of supercurrents in Josephson junctions made on bulk-insulating TI samples are recently reported, the Fraunhofer patters observed in such TI-based Josephson junctions, which sometimes present anomalous features, are still not well understood. Here we report our study of highly gate-tunable TI-based Josephson junctions made of one of the most bulk-insulating TI materials, BiSbTeSe_2, and Al. The Fermi level can be tuned by gating across the Dirac point, and the high transparency of the Al/BiSbTeSe_2 interface is evinced by a high characteristic voltage and multiple Andreev reflections with peak indices reaching 12. Anomalous Fraunhofer patterns with missing lobes were observed in the entire range of gate voltage. We found that, by employing an advanced fitting procedure to use the maximum entropy method in a Monte Carlo algorithm, the anomalous Fraunhofer patterns are explained as a result of inhomogeneous supercurrent distributions on the TI surface in the junction. Besides establishing a highly promising fabrication technology, this work clarifies one of the important open issues regarding TI-based Josephson junctions.


Anomalous Fraunhofer patterns in gated Josephson junctions based on the bulk-insulating topological insulator BiSbTeSe_2
    Yoichi Ando
    December 30, 2023
========================================================================================================================




The superconducting proximity effect in the topological surface states of a three-dimensional (3D) topological insulator (TI) is of significant current interest, because it provides a promising means to engineer a two-dimensional (2D) topological superconducting state <cit.>. Because of the spin-momentum locking <cit.>, the spin degrees of freedom is frozen in the topological surface states and “spinless" superconductivity is established when pairing is induced in surface-state electrons  proximitized by a conventional s-wave superconductor. The spinless nature is the key to realizing topological superconductivity associated with Majorana fermions <cit.>, of which particle and hole excitations are identical. When Majorana fermions are spatially localized, they become Majorana zero-modes which obey non-Abelian statistics and are useful for topological quantum computation <cit.>. Hence, clarifying the nature of the proximity-induced superconductivity in the topological surface states of a 3D TI is not only interesting for finding topological superconductivity <cit.> but also crucial for high-impact applications.

Josephson junctions (JJs) are the basic devices to investigate the superconducting proximity effect in 3D TIs. In a JJ, Andreev bound states formed in the weak link mediate dissipationless currents <cit.>. For a JJ based on a 3D TI, a pair of Majorana bound states, whose energy changes with the phase difference ϕ between the two superconductor (SC) electrodes with a 4π periodicity, are predicted to form  <cit.>. In the past, a variety of SC-TI-SC JJs <cit.> have been experimentally investigated to address these Majorana bound states. Although many of these experiments were carried out in JJs made of bulk-conducting materials <cit.>, where a large part of the supercurrent is carried by conventional Andreev bound states, there have been several studies which were based on JJs made on bulk-insulating TIs <cit.>. However, the transparency of the interface between the SC and the TI surface states have not been so high compared to the case of InAs nanowires proximitized by epitaxially-grown Al <cit.>. Also, the experiments which employed electrostatic gating to approach the Dirac point have been relatively limited <cit.>.




Here, we present our study of highly gate-tunable Josephson junctions based on one of the most bulk-insulating TI materials, BiSbTeSe_2 <cit.>. A high characteristic voltage (expressed in the so-called “eI_cR_N/Δ ratio" of ∼0.4, among the largest ever reported <cit.>) demonstrates a high transparency of our SC/TI contact; this is also manifested in our observation of multiple Andreev reflections <cit.> with the peak index reaching 12. By varying the gate voltage V_g, the Fermi level E_F of BiSbTeSe_2 can be tuned in a broad range from inside the conduction band to below the Dirac point, establishing that proximity-induced surface superconductivity can be realized even at the Dirac point. Intriguingly, anomalous Fraunhofer patterns with missing or significantly suppressed lobes are observed in the entire range of V_g. Similar irregular Fraunhofer patterns have been reported in several cases <cit.>, but the origin has not been clarified. In this regard, it was claimed that anomalous Fraunhofer patterns observed in TI-based JJs could be a signature of Majorana bound states <cit.>. In the present work, we have been able to fit the anomalous Fraunhofer patterns by using a Monte Carlo algorithm involving the maximum entropy method and conclude that they result from an inhomogeneous supercurrent distribution.
In fact, it is understandable that the effect of current inhomogeneity becomes manifest in the surface-dominated transport regime where there is no bulk shunt.
This clarifies one of the important open issues regarding TI-based JJs, and our JJ fabrication technology paves the way for future Majorana devices in SC-TI hybrid systems.



The TI-based JJs reported here are fabricated on bulk-insulating BiSbTeSe_2 flakes.
A scanning electron microscope (SEM) image of a typical device is shown in Figure 1a as an example, and the device structure is schematically illustrated in Figure 1b. The flakes of BiSbTeSe_2 were exfoliated from a bulk single crystal and then transferred onto degenerately-doped Si substrates having a 290-nm-thick SiO_2 coating layer, which serves as the dielectric of the back gate. Flakes with a thickness t of ∼20 nm or thinner were selected for device fabrication, so that the E_F of the top surface can be controlled with the back gate. The JJ is defined by patterning two parallel Pt/Al/Au electrodes with a spacing L to cover the full width W of a selected flake; we call L and W the length and the width of the junction, respectively. Here the thin Pt layer works as a buffer layer to improve the transparency of Al/BiSbTeSe_2 contacts, and the Au capping layer helps to prevent Al electrodes from oxidization in air. We have investigated four JJ devices (denoted devices 1 – 4) with essentially the same design (see Supporting Figure S1 for their pictures). The dimensions of all the devices are shown in Table S1 of the Supporting Information. The measurements were performed in a pseudo-four-wire configuration as illustrated in Figure 1b, by using both ac and dc techniques (see detailed descriptions at the end).

The temperature dependence of the junction resistance R in device 1 is shown in Figure 1c. Superconducting transitions occur successively upon cooling: First, a  sharp drop in R is observed at T∼1.4 K, and then a broad transition continues down to 750 mK, below which the junction reaches a zero-resistance state. The first transition corresponds to the superconducting transition of Al electrodes, and the second transition is associated with the onset of Josephson coupling between the two Al electrodes.

A typical I-V curve of device 1 measured at V_g = 0 V is plotted in Figure 1d. For an applied dc current I_dc smaller than the critical current I_c=0.39 μA, the junction is in the superconducting state and the measured voltage V_dc remains zero. When I_dc exceeds I_c, a finite V_dc starts to develop, and the junction is driven into a resistive state. However, due to quasiparticle processes such as multiple Andreev reflection, the I-V curve is nonlinear at low V_dc. The linearity in the I-V curve gradually recovers when V_dc goes beyond 2Δ_Al/e; the inverse of the slope in this linear region gives a normal-state resistance R_N^I-V=186 Ω. An excess current I_e=0.94 μA is obtained by linearly extrapolating the I-V curve from the linear region to V_dc=0V, as illustrated in Figure 1c. Based on the Octavio-Tinkham-Blonder-Klapwijk theory <cit.>, the transparency T_J of the Al/Pt/BiSbTeSe_2 contact can be estimated from the values of R_N^I-V and I_e <cit.>, and T_J=0.77 is obtained for device 1. (see Supporting Information for details).

The I-V curves of our JJs show no hysteretic behaviour; the values of I_c obtained in up- and down-sweeps are identical. Such non-hysteretic I-Vs are usually seen in overdamped JJs. The Stewart-McCumber parameter of device 1 is estimated to be β_c=2eI_c(R_N^I-V)^2C/ħ≈ 1, where C is the shunt capacitance between the two Al electrodes. Similar values of β_c are obtained for other devices (Supporting Information). Such intermediate values of β_c indicate that our JJs lie in the region between the overdamped and underdamped limits. According to numerical calculations, JJs in such a region show non-hysteretic I-V characteristics similar to those of underdamped JJs <cit.>, which is consistent with our observation.


Sub-harmonic peaks at V_n = 2Δ / ne are often observed in the dI/dV vs V_dc curves of voltage-biased JJs which satisfy L<l_ϕ. Here, 2Δ is the superconducting energy gap, n=1,2... is the peak index and l_ϕ is the inelastic scattering length of the normal metal. Such sub-harmonic features are associated with multiple Andreev reflection (MAR) processes <cit.> taking place in the channel of a JJ. In the past, the MAR peaks have been observed in a variety of TI-based JJs <cit.>.

In the present experiment, clear MAR peaks were seen in the dI/dV vs V_dc curves of all four devices. A typical curve measured in device 2 is plotted in Figure 1e with assigned peak indices. The assignments of n are confirmed by the linear relation between  V_n and 1/n illustrated in the inset of Figure 1e.  A linear fit to the V_n(1/n) data gives Δ_Al=160 μeV, consistent with the superconducting gap Δ_Al reported for Al in the literature <cit.>.

A key feature in Figure 1e is the presence of high-index conductance peaks reaching n = 12, indicating that the inelastic scattering length l_ϕ in BiSbTeSe_2 exceeds 12L (= 840 nm), which is longer than what was reported for previous TI-based JJs <cit.>. Nevertheless, such a long l_ϕ is consistent with previous results of Aharonov-Bohm oscillations observed in square-ring interferometers based on  Bi_2Se_3 flakes <cit.>. In addition, the presence of high-index peaks also indicates a relatively high transparency T_J of the Al/BiSbTeSe_2 contacts. Numerical calculations <cit.> show that high-index features of MAR become the most visible when T_J takes a relatively high value of around 0.7-0.8, while they disappear in both transparent (T_J=1) and tunneling (T_J=0) limits, which is consistent with our estimate of T_J = 0.83 for device 2.

The amplitude of the conductance peak at V_n is expected to decay with increasing n, because the probability of MAR processes generally decreases with n. However, in our experiment, the n=12 peak shows an amplitude much larger than those of other peaks. The origin of the larger-than-expected amplitude of n=12 peak is not clear and the reason for the absence of peaks for n=7-11 is unknown. Nevertheless, similar MAR phenomena have been seen in previous experiments on TI-based JJs by Jauregui et al. <cit.> and on Ge/Si-nanowire-based JJs by Xiang et al. <cit.>.



As-grown BiSbTeSe_2 single crystals are bulk-insulating, with E_F located slightly below the Dirac point <cit.>. However, due to the unintentional doping introduced during the exfoliation and device fabrication processes, E_F in our samples is shifted into the n-type region. Nevertheless, the top-surface Fermi level of our devices can be tuned by the back-gate voltage V_g to the Dirac point, thanks to the small thickness (∼18 nm) of the TI flakes used here; in the following, V_DP denotes the V_g value needed to bring the Fermi level to the Dirac point. For a detailed analysis of the single-gate control of the Fermi levels of two surfaces involving a capacitive coupling, see Ref. .

The normal-state resistance R_N vs V_g curve of device 1 is plotted in Figure 1g. Here, R_N was measured by suppressing the superconductivity of Al with a magnetic field; note that R_N and R_N^I-V are not exactly the same, which probably indicates that the contact resistance is very sensitive to the magnetic field. The position of the Dirac point V_DP=-88V is determined by the peak in this curve, while the position of the conduction band bottom V_CB is estimated by using V_CB=en_2d/C_g+V_DP≈ 33V, where n_2d≈ 9.1 × 10^12 cm^-2 is the two-dimensional carrier density in BiSbTeSe_2 when E_F reaches the bottom of the conduction band, and C_g≈12 nF/cm^2 is the capacitance per unit area of the back gate (see Supporting Information for details). 

Figure 1f shows the colour mapping of dV/dI in device 1 scanned in the I_ dc vs V_ g plane. The critical current I_c is depicted by the white band at the edge of the superconducting region in deep blue colour. A finite I_c of 92 nA is observed at the Dirac point. In the n-type region where V_g > V_DP, I_c increases sharply with increasing V_g. However, in the p-type region where V_g<V_DP, I_c saturates when V_g decreases. In device 2 where E_F can be tuned deeper into the p-type region, we found that I_c actually increases slowly with decreasing V_g in the p-type region, although the rate of increase is very low (see Supporting Figure S2c). In the past, similar asymmetric I_c(V_g) behavior was also reported in TI-based JJs <cit.> and graphene <cit.>, and it was speculated that this behaviour is caused by the p-n junctions formed between the gated junction channel and the areas under the SC contacts.

The I_cR_N product of device 1 is calculated as a function of V_g using the data in Figures 1f and 1g, and shown in Figure 1h. The I_cR_N product is almost independent of V_g except for the region near the Dirac point. When E_F approaches the Dirac point from the n-type region, I_cR_N presents a minimum at around V_g-V_DP≈35V.
Similar behavior in I_cR_N was observed with less clarity in device 2
(see Supporting Figure S2d). The origin of this minimum is not clear. Most importantly, device 1 shows the eI_cR_N/Δ ratio of up to 0.37, which is among the largest values reported in the literature <cit.>, again proving the high transparency of the Al/BiSbTeSe_2 contacts realized in the present work.






The temperature dependence of I_c was measured at various V_g values, as plotted in Figure 2a. With decreasing T, I_c first increases exponentially and then saturates. The step-like increase in I_c at low temperature reported in a recent work on TI-based JJs using Nb electrodes <cit.> was not observed in our devices. A fit to I_c∝√(T) exp(-2π k_BT /E_Th) in the high-temperature range allows us to extract the Thouless energy E_Th (≡ħ D/L^2) at various V_g <cit.>; here, D is the average diffusion constant of supercurrent-carrying electrons and L is the length of the JJ.

The E_Th values obtained from the fittings are plotted in Figure 2b as a function of V_g. It is noteworthy that E_Th remains unchanged when E_F moves through the bottom of the conduction band at V_g=V_CB. In other words, the presence of bulk carriers does not change the average D of the supercurrent-carrying electrons. Since D of bulk carriers is expected to be much smaller than that of the surface electrons, the unchanged D suggests that the supercurrent is mainly carried by the surfaces states even when E_F is located inside the conduction band. This is in agreement with the conclusions in previous studies on planar TI-based JJs <cit.>, but contradicts the results obtained in Pb-Bi_2Se_3-Pb JJs fabricated along the c-direction of Bi_2Se_3 flakes <cit.>. As E_F approaches the Dirac point from the n-type region, E_Th starts to increase, indicating an increase in D of the surface electrons. Such an enhancement of D around the Dirac point is probably due to the smaller phase space available for scattering. In addition, the extracted E_Th is comparable to Δ_Al in the entire range of V_g, indicating that our junctions lie on the boundary between the short- and long-junction limit, similar to what was previously reported in Al-Bi_2Se_3-Al junctions<cit.>.


To determine the regime to which our JJs belong, it is important to estimate several characteristic lengths.
As discussed above, the lower boundary of the inelastic scattering length l_ϕ is estimated from the MAR results to be l_ϕ > 12L = 840 nm.
The superconducting coherence length ξ_N in BiSbTeSe_2 is given by ξ_N=√(ħ D/Δ_Al)=√(E_ThL^2/Δ_Al). With L = 70 nm measured in the SEM image, Δ_Al=160 μeV obtained from the MAR data, and E_Th = 187 – 236 μeV extracted from the I_c(T) curves, we obtain ξ_N=76-85nm.
The mean free path l_e of the supercurrent-carrying surface electrons can be calculated by using E_Th=ħ D/L^2 and D=v_Fl_e/2, where v_F=5.5×10^5m/s is the Fermi velocity of the surface states <cit.>, and one obtains l_e = 5.1 – 6.4 nm, which is much shorter than the previously reported value l_e∼40nm in Bi_1.5Sb_0.5Te_1.7Se_1.3 bulk crystals <cit.>. We tentatively attribute the short l_e in the BiSbTeSe_2 flakes to the surface disorder introduced in the device fabrication process. While a short l_e might seem a disadvantage for quantum computing, the readout of Majorana qubits can be performed without relying on ballistic transport <cit.>. 
The characteristic lengths estimated above point to the relation l_e≪ L ∼ξ_N≪ l_ϕ, which indicates that our JJs are diffusive junctions (L≫ l_e) in the intermediate region between the short- and long-junction limit (L ∼ξ_N).




For JJs satisfying L ≤ξ_N, the eI_cR_N/Δ ratio is expected to be of order 1 at T=0K <cit.>. However, in reality, the measured values of eI_cR_N/Δ are usually much smaller than the theoretical expectation due to the existence of resistive shunt channels or limited SC contact transparency <cit.>. The eI_cR_N/Δ ratio found in previous studies on TI-based JJs ranges between 0.01 and 0.31, all smaller than the ratio of 0.37 observed in device 1, indicating that the quality of our TI-based JJs is among the best in terms of bulk-insulating property and contact transparency.


For a trivial JJ with a homogeneous current distribution, the I_c vs B relation is expected to exhibit a so-called Fraunhofer pattern described by

    I_c(B)=I_c0|sin(πΦ_J/Φ_0)/πΦ_J/Φ_0|,

where I_c0 is the critical current at B=0, Φ_J is the flux threading through the junction channel and Φ_0=h/(2e) is the flux quantum.
To elucidate the I_c(B) behaviour of our JJs, the dV/dI data were measured as a function of I_dc and B at various V_g corresponding to the E_F positions from inside the conduction band to near the Dirac point. The I_c(B) patterns of device 1 are shown in Figures 3a–3d. All the observed patterns are essentially Fraunhofer-like, but bear anomalous features that deviate from a standard Fraunhofer pattern. They include: (i) an average period Δ B_avg≈ 1.44mT, which is about three times smaller than the expected period Δ B_0≈ 5.29mT estimated from the junction area S=(L+2λ_Al)W, where λ_Al (≈ 50nm) is the London penetration depth of bulk Al, and (ii) missing or significantly suppressed lobes which cannot be described by the formula for the Fraunhofer pattern.

Anomalous features in the Fraunhofer patterns have been reported for TI-based JJs in some of the previous works <cit.>. There was an attempt to explain the smaller-than-expected Δ B_avg by a phenomenological model involving quantized Majorana modes <cit.>. However, in the present case, it may partly be explained by a small 2D superfluid density (and hence longer effective λ_Al) in the 30-nm-thick Al electrodes which are thinner than the bulk λ_Al.
Also, it has been noted that the small Δ B_avg can be explained by the flux focusing effect <cit.>, which is commonly seen in planar JJs <cit.>. Due to the diamagnetism of a SC, the magnetic field near the SC electrodes is distorted, and half of the magnetic flux expelled from the SC electrodes is compressed into the junction; this effect may be quantified in terms of an extended effective junction area S_eff=WL_eff, where L_eff=(L+d_1/2+d_2/2) is the effective junction length considering the flux focusing (d_1 and d_2 are the widths of the left and right SC electrodes). Taking d_1 = d_2 =400nm and L=70nm for device 1, we obtain Δ B=Φ_0/(S_eff)=1.52mT, which is close to the observed average period Δ B_avg≈ 1.44mT. Therefore, in the following, we use S_eff instead of S for the reconstruction of the I_c(B) patterns.

Regarding the origin of the other feature, i.e. missing lobes in the Fraunhofer pattern, this has been left unexplained in the literature <cit.>.
Theoretically, an irregular I_c(B) pattern could be generated by (i) a non-sinusoidal current phase relation or (ii) an inhomogeneous supercurrent distribution in the JJ. In the past, much attention has been paid to the first possibility. For instance, it was predicted that the formation of 1D Majorana modes in a TI-based JJ gives rise to an I_c(B) pattern featuring a non-zero first minimum, deviating from a standard Fraunhofer pattern <cit.>. It was also predicted that the 4π periodicity due to Majorana zero-modes leads to a suppression of the odd-numbered (1st, 3rd, 5th, etc.) lobes in the Fraunhofer patttern <cit.>.

However, these theoretical possibilities do not really explain the disappearance of some particular lobes, and one had better consider the second possibility to understand the irregular I_c(B) patterns observed in our TI-based JJs. We set out to reconstruct the anomalous Fraunhofer patterns by considering an inhomogeneous supercurrent distribution j(x) along the width of the junction. The experimental data of I_c(B), extracted from the dV/dI(I_dc) data as the onset of rapid increase, are fitted by using a Monte Carlo algorithm to optimize the form of j(x) based on the maximum entropy method (see detailed descriptions at the end), which has not been applied to the present problem; this method removes the constraint of the standard Dynes-Fulton-type analysis <cit.> and allows to identify any form of inhomogeneous supercurrent distribution.
In Figure 3e we plot for three V_ g values the obtained forms of j(x), which reproduce all the main features in the I_c(B) patterns as shown with red dashed lines in Figures 3a–3c. In the normalized plot of j(x)/j_ avg shown in Figure 3f, one can see that the relative magnitude of the fluctuations is enhanced at V_ g = -60 V, which is close to the Dirac point. The enhancement of j(x) at the edges seen in Figure 3f is probably due to the contribution of the side surfaces <cit.> which also contain topological surface states.
We did not perform a fitting of the I_c(B) data for V_ g = -80 V, because the I_c values at this V_ g could not be reliably extracted due to the absence of true supercurrents in magnetic fields of only a few mT.

We note that the j(x) curves presented in Figure 3e are not the unique solutions which give good fittings to the data. This is partly because of the lack of the phase information in this inverse problem, but also because of the trade-off between the degree of fitting and the smoothness of the solution. Nevertheless, it is clear from the present analysis that the existence of an inhomogeneous j(x) is sufficient for causing the observed anomalous Fraunhofer patterns, most notably missing lobes.

An inhomogeneous j(x) can be caused by various reasons. One possible reason is the fluctuations in the junction length L(x) due to the roughness of the edges of Al electrodes. Since the critical current depends exponentially on L, for junctions with a short L, even small fluctuations in L(x) could result in a notable variation in j(x). Another possible reason is a position-dependent transparency of the Al/TI interface, which may arise from a granularity of the Al film. In addition, a possible variation of the surface carrier density due to charge puddles <cit.> may also be partially responsible for the inhomogeneous j(x).


To conclude, we have established a highly promising fabrication technology for TI-based JJs, which allowed us to elucidate that the anomalous Fraunhofer patterns are essentially due to an inhomogeneous supercurrent distribution in the proximitized surface states. This technology and the quantitative understanding gained in this work would make it possible to fabricate an optimized device to nail down the signature of Majorana fermions in the proximitized TI surface in near future.





Growth of bulk BiSbTeSe_2 crystals:
The bulk BiSbTeSe_2 single crystals were grown by melting high-purity (99.9999%) Bi, Sb, Te, and Se shots with the
molar ratio of 1:1:1:2 in an evacuated quartz tube via a modified Bridgman method, following our established recipe reported in Ref. . The grown crystals typically have a domain size of several mm and are easily cleaved along the (111) plane.

Device fabrication:
Ultrathin BiSbTeSe_2 flakes are obtained by mechanical exfoliation of bulk single crystals on degenerately-doped Si/SiO_2 wafers by a Scotch-tape technique. Promising flakes are initially identified under a laser microscope (Keyence VK-X200), and further characterized by atomic force microscopy. The flakes with an extremely smooth surface (without atomic terraces) are chosen for the fabrication of JJ devices using electron beam lithography. Previous studies on a similar bulk-insulating TI material have shown <cit.> that a very short junction channel (∼50 nm) is essential for realizing a sufficiently strong coupling between two SC electrodes to observe a Josephson supercurrent. To reliably fabricate short-channel junctions, we have developed our own recipe. A 160-nm-thick PMMA  layer was spin coated and baked at 120^∘C for 10 minutes. The reduced baking temperature ensures minimal unintentional doping due to creation of Se vacancies at elevated temperatures. The resist layer was exposed by a 25-kV electron beam. The width of the Al electrodes was kept around 400 nm to avoid vortex movements at finite magnetic fields in the superconducting state. To remove any resist residues in the contact areas after the development, the wafers were cleaned with oxygen plasma in a reactive-ion etching machine at 20 W for 7 s. Afterwards, the wafers were further cleaned in a HCl solution to remove any oxide layer in the contact areas. The metallization was carried out by sputtering 5-nm-thick Pt as buffer layer, followed by thermal evaporation of 30-nm-thick Al at the pressure of 5×10^-7 mbar. A 5-nm-thick Au layer was successively evaporated to protect the Al from oxidation, and the lift-off was done in warmed acetone.

Measurements:
The devices were cooled in a dry dilution refrigerator (Oxford Instruments TRITON 200) with a base temperature of less than 10 mK. To reduce the electromagnetic noise, the electrical lines are equipped with RC-filters at room temperature and at the 4-K plate, along with additional RC and copper-powder filters at the mixing chamber plate. The junctions were electrically characterized by dc current-voltage (I-V) and differential resistance (dV/ dI) measurements in a pseudo-four-probe configuration. For the dV/ dI measurement, an 18.6-Hz ac current with the amplitude of 1 nA was superimposed to a dc bias current. The magnetic-field dependence of the critical current (i.e. Fraunhofer pattern) was measured by performing the dV/ dI measurements while driving the superconducting magnet with a Keithley 2450 source meter.

Current density retrieval:
For a JJ with a sinusoidal current-phase relation, the experimentally observed critical current I_c is related to the modulus of the phase-dependent integration of the supercurrent density <cit.>,

    I_c^FT(B)= | ∫_-W/2^W/2 j(x) e^i φ(x,B) dx |  ,

which allows to deduce information about the distribution of the current density j and the phase difference φ within the junction.
In the present case of a narrow JJ (i.e. W ≪λ_J, where λ_J is the Josephson penetration length), φ can be assumed to depend linearly on the spatial position x along the junction and on the applied magnetic field B; namely, φ(x,B)=2π x L_eff B / Φ_0  + φ_0, where Φ_0=h/2e is the magnetic flux quantum and L_eff=L+2 ·d/2 is the effective length of the JJ to take into account the flux focusing effect of the superconducting electrode of width d <cit.>.
This phase factor makes the integration to be effectively a Fourier transform from the x coordinate to the B coordinate. In general, the supercurrent density J(x,y) may depend on both the spatial coordinates x and y within the junction area; however, we restrict ourselves to the dependence on x defined along the junction's width W, by considering the one-dimensional current density j(x)=∫_-L/2^L/2 J(x,y) dy.

As the critical current does not contain information on the complex phase of the Fourier transform, but only on its modulus, the current density j(x) cannot be obtained unambiguously from the experimental data I_c^exp(B). In fact, there is an infinitely large number of current profiles j(x) that yield a similar I_c^FT(B).
Most of previous reconstructions of the current distribution were essentially based on the Dynes-Fulton analysis <cit.>. Although the original Dynes-Fulton analysis utilized a relatively strong assumption that j(x) is an approximately even function <cit.>, some of the later attempts have relaxed this constraint <cit.>. 

In this work, we perform a novel type of analysis based on the ideas of Dynes and Fulton, but without requiring any presumptions; specifically, we identify the most physical solution by using a maximum entropy method <cit.>.
The outline of our fitting procedure is the following:
First, we discretize j(x) into N elements j_i and initialize them to a value of I_c^exp(0)/N. To ensure the agreement with the zero-field critical current, we keep ∑_i j_i=I_c^exp(0) to be fixed. Then we fit the data by maximizing the functional ℱ[j]= -λχ_j^2-∑_i j̃_i lnj̃_i using a Monte Carlo algorithm.
Here, χ_j^2=∑_k=0^M( I_c^FT(B_k)-I_c^exp(B_k))^2 is the square of the deviation from the experimental data taken at points B_k, and the second term is the entropy of the normalized current distribution j̃_i=j_i/∑_i j_i.
The weighting factor λ tunes between a pure least-square fit for λ→∞ and an increasing contribution of the maximum entropy method for λ→ 0. A larger λ gives a better fitting, but the obtained current distribution tends to become unphysically peaky; more weighting of the entropy term makes the resulting current distribution to be smoother and more physical, even though χ_j^2 becomes larger.
In our analysis, the value of λ is initially set to 1/(2M), and we increase (decrease) it during the Monte Carlo steps when the iteration goes to the wrong direction, i.e., when χ_j^2 becomes larger (smaller) than a given target value ξ; such a dynamical tuning of λ allows the square term χ_j^2 to converge to ξ while ensuring the resulting current distribution to be smooth.
Still, the algorithm does not provide a single and unique solution for the current distribution, and the smoothness of the distribution depends on the choice of ξ. To avoid over-fitting of the data, for each data set the maximal ξ was chosen so that the calculated  I_c^FT(B) reproduces the absent lobes and the overall shape of the Fraunhofer pattern. Typically, the fits converged after about 50,000 steps. Each data set was processed about 100 times to check for the robustness of the fitting result, and we found that the obtained current distributions – save for mirror symmetry – resembled each other within a few percent of deviation. From these we pick the solution having the largest entropy and show it in Figure 3.


Author contributions:
S.G., O.B. and F.Y contributed equally to this work and should be considered as co-first authors. 

Acknowledgement:
This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 741121) and was also supported by DFG (CRC1238 “Control and Dynamics of Quantum Materials", Projects A04 and B01). We thank D. Goldhaber-Gordon, A. Rosch, S. Trebst, and Y. Vinkler for helpful discussions.
O.B. acknowledges the support from Quantum Matter and Materials Program at the University of Cologne funded by the German Excellence Initiative.










10

TI-JJ-theory-Fu Fu, L.; Kane, C. L. Superconducting Proximity Effect and Majorana Fermions at the Surface of a Topological Insulator. Phys. Rev. Lett. 2008, 100, 96407.

Sato-Ando
Sato, M.; Ando, Y. Topological Superconductors: A Review. Rep. Prog. Phys. 2017, 80, 076501.

TI-review-Hasan Hasan, M. Z.; Kane, C. L. Colloquium: Topological insulators. Rev. Mod. Phys. 2010, 82, 3045-3067.


TI-review-Qi Qi, X.-L.; Zhang, S.-C. Topological Insulators and Superconductors. Rev. Mod. Phys. 2011, 83, 1057-1110.


TI-review-Ando Ando, Y. Topological Insulator Materials. J. Phys. Soc. Jpn. 2013, 82, 102001.


MF-review-Elliott Elliott, S. R.; Franz., M. Colloquium: Majorana Fermions in Nuclear, Particle, and Solid-State Physics. Rev. Mod. Phys. 2015, 87, 137-163.

QuantumCpt-review-Sarma Sarma, S. D.; Freedman, M.; Nayak, C. Majorana Zero Modes and Topological Quantum Computation npj Quantum Inf. 2015, 1, 15001.

ABS-Klapwijk Klapwijk, T. M. Proximity Effect from an Andreev Perspective. J. Supercond. Novel Magn. 2004, 17, 593-611.


Snelder
Snelder, M.; Veldhorst, M.; Golubov, A. A.; Brinkman, A. Andreev Bound States and Current-Phase Relations in Three-Dimensional Topological Insulators. Phys. Rev. B 2013, 87, 104507.

TI-JJ-Samarth Zhang, D.; Wang, J.; DaSilva, A. M.; Lee, J. S.; Gutierrez, H. R.; Chan, M. H. W.; Jain, J.; Samarth, N. Superconducting Proximity Effect and Possible Evidence for Pearl Vortices in a Candidate Topological Insulator. Phys. Rev. B 2011, 84, 165120.

TI-JJ-exp-Brinkman Veldhorst, M.; Snelder, M.; Hoek, M.; Gang, T.; Guduru, V. K.; Wang, X. L.; Zeitler, U.; van der Wiel, W. G.; Golubov, A. A.; Hilgenkamp, H.; et al. Josephson Supercurrent through a Topological Insulator Surface State. Nat. Mater. 2012, 11, 417-421.

TI-JJ-exp-LiLu Qu, F.; Yang, F.; Shen, J.; Ding, Y.; Chen, J.; Ji, Z.; Liu, G.; Fan, J.; Jing, X.; Yang, C.; et al. Strong Superconducting Proximity Effect in Pb-Bi_2Te_3 Hybrid Structures. Sci. Rep. 2012, 2, 339.

TI-JJ-exp-GoldhaberGordon Williams, J. R.; Bestwick, A. J.; Gallagher, P.; Hong, S. S.; Cui, Y.; Bleich, A. S.; Analytis, J. G.; Fisher, I. R.; Goldhaber-Gordon, D. Unconventional Josephson Effect in Hybrid Superconductor-Topological Insulator Devices. Phys. Rev. Lett. 2012, 109, 56803.

TI-JJ-exp-Lombardi Galletti, L.; Charpentier, S.; Iavarone, M.; Lucignano, P.; Massarotti, D.; Arpaia, R.; Suzuki, Y.; Kadowaki, K.; Bauch, T.; Tagliacozzo, A.; et al. Influence of Topological Edge States on the Properties of Al/Bi_2Se_3/Al Hybrid Josephson Devices. Phys. Rev. B 2014, 89, 134512.

TI-JJ-exp-Juelich Schüffelgen, P.; Rosenbach, D.; Li, C.; Schmitt, T.; Schleenvoigt, M.; Jalil, A. R.; Kölzer, J.; Wang, M.; Bennemann, B.; Parlak, U.; et al. Boosting Transparency in Topological Josephson Junctions via Stencil Lithography. 2017, arXiv:1711.01665 [cond-mat.supr-con].

TI-JJ-exp-Molenkamp Oostinga, Jeroen B.; Maier, Luis; Schüffelgen, Peter; Knott, Daniel; Ames, Christopher; Brüne, Christoph; Tkachov, Grigory; Buhmann, Hartmut; Molenkamp, Laurens W. Josephson Supercurrent through the Topological Surface States of Strained Bulk HgTe. Phys. Rev. X 2013, 3, 21007.

TI-JJ-exp-Mason Cho, S.; Dellabetta, B.; Yang, A.; Schneeloch, J.; Xu, Z.; Valla, T.; Gu, G.; Gilbert, M. J.; Mason, N. Symmetry Protected Josephson Supercurrents in Three-Dimensional Topological Insulators. Nat. Commun. 2013, 4, 1689.

TI-JJ-exp-Lee Lee, J.H.; Lee, G.H.; Park, J.; Lee, J.; Nam, S.G.; Shin, Y.S.; Kim, J.S.; Lee, H.J. Local and Nonlocal Fraunhofer-like Pattern from an Edge-Stepped Topological Surface Josephson Current Distribution. Nano Lett. 2014, 14, 5029-5034.

TI-JJ-exp-Brinkman2 Snelder, M.; Molenaar, C. G.; Pan, Y.; Wu, D.; Huang, Y. K.; de Visser, A.; Golubov, A. A.; van der Wiel, W. G.; Hilgenkamp, H.;  Golden, M. S.; et al. Josephson Supercurrent in a Topological Insulator without a Bulk Shunt. Supercond. Sci. Technol. 2014, 27, 104001.

Wiedenmann
Wiedenmann, J.; Bocquillon, E.; Deacon, R. S.; Hartinger, S.; Herrmann, O.;Klapwijk, T. M.; Maier, L.; Ames, C.; Brüne, C.; Gould, C.; et al. 4π-periodic Josephson Supercurrent in HgTe-based Topological Josephson Junctions. Nat. Commun. 2016, 7, 10303.

TI-JJ-exp-Chen Jauregui, L. A.; Kayyalha, M.; Kazakov, A.; Miotkowski, I.; Rokhinson, L. P.; Chen, Y. P.; Gate-tunable Supercurrent and Multiple Andreev Reflections in a Superconductor-Topological Insulator Nanoribbon-Superconductor Hybrid Device. Appl. Phys. Lett. 2018, 112, 093105.

Bocquillon
Bocquillon, E.; Deacon, R. S.; Wiedenmann, J.; Leubner, P.; Klapwijk, T. M.; Brüne, C.; Ishibashi, K.; Buhmann, H.; Molenkamp, L. W. Gapless Andreev Bound States in the Quantum Spin Hall Insulator HgTe. Nat. Nanotechnol. 2017, 12, 137-143.

Krogstrup
Krogstrup, P.; Ziino, N. L. B.; Chang, W.; Albrecht, S. M.; Madsen, M. H.; Johnson, E.; Nygård, J.,;Marcus, C. M.; Jespersen, T. S. Epitaxy of Semiconductor-Superconductor Nanowires. Nat. Mater. 2015, 14, 400-406.

RenBSTS
Ren, Z.; Taskin, A. A.; Sasaki, S.; Segawa, K.; Ando, Y. Optimizing Bi_2-xSb_xTe_3-ySe_y Solid Solutions to Approach the Intrinsic Topological Insulator Regime. Phys. Rev. B 2011, 84, 165311.

MAR-theory-Klapwijk Klapwijk, T. M.; Blonder, G. E.; Tinkham, M. Explanation of Subharmonic Energy Gap Structure in Superconducting Contacts. Physica B+C 1982, 109, 1657-1664.

Suominen
Suominen, H. J.; Danon, J.; Kjaergaard, M.; Flensberg, K.; Shabani, J.; Palmstrøm, C. J.; Nichele, F.; Marcus, C. M. Anomalous Fraunhofer Interference in Epitaxial Superconductor-Semiconductor Josephson Junctions. Phys. Rev. B 2017, 95, 035307.

OTBK-theory-Octavio Octavio, M.; Tinkham, M.;  Blonder, G. E.; Klapwijk, T. M. Subharmonic Energy-Gap Structure in Superconducting Constrictions. Phys. Rev. B 1983, 27, 6739.

OTBK-theory-Flensberg Flensberg, K.; Hansen, J. B.; Octavio, M. Subharmonic Energy-Gap Structure in Superconducting Weak Links. Phys. Rev. B 1988, 38, 8707.

JJ-book-Likharev Likharev, K. K. Dynamics of Josephson Junctions and Circuits; Gordon and Breach Science Publishers, 1986.

AB-effect-LiLu  Qu, F.;Yang, F.;Chen, J.;Shen, J.;Ding, Y.;Lu, J.;Song, Y.;Yang, H.;Liu, G.;Fan, J.; et al. Aharonov-Casher Effect in Bi_2Se_3 Square-Ring Interferometers. Phys. Rev. Lett. 2011, 107, 16802.

MAR-peak-index-Yeyati Cuevas, J. C.; Martin-Rodero, A.; Yeyati, A. L. Hamiltonian Approach to the Transport Properties of Superconducting Quantum Point Contacts. Phys. Rev. B 1996, 54, 7366.

MAR-peak-index-Bardas Averin, D.; Bardas, A. ac Josephson Effect in a Single Quantum Channel. Phys. Rev. Lett. 1995, 75, 1831.

Xiang2006 Xiang, J.; Vidan, A.; Tinkham, M.; Westervelt, R. M.; Lieber, C. M., Ge/Si nanowire mesoscopic Josephson junctions. Nature Nanotech. 2006, 1, 208.

BSTS2-sato Arakane, T.; Sato, T.; Souma, S.; Kosaka, K.; Nakayama, K. ; Komatsu, M.; Takahashi, T.; Ren, Z.; Segawa, K.; Ando, Y. Tunable Dirac Cone in the Topological Insulator Bi_2-xSb_xTe_3-ySe_y. Nat. Commun. 2012, 3, 636.

Yang2015 Yang, F.; Taskin, A. A.; Sasaki, S.; Segawa, K.; Ohno, Y.; Matsumoto, K.; Ando, Y., Dual-Gated Topological Insulator Thin-Film Device for Efficient Fermi-Level Tuning. ACS Nano 2015, 9, 4050-4055.

Graphine-JJ-exp-Geim Ben Shalom, M.; Zhu, M. J.; Fal'ko, V. I.; Mishchenko, A.; Kretinin, A. V.; Novoselov, K. S.; Woods, C. R.; Watanabe, K.; Taniguchi, T.; Geim, A. K.; et al. Quantum Oscillations of the Critical Current and High-Field Superconducting Proximity in Ballistic Graphene. Nat. Phys. 2016, 12, 318-322.



flux-focusing-Rosenthal Rosenthal, P. A. et al. Flux Focusing Effects in Planar Thin-Film Grain-Boundary Josephson Junctions. Appl. Phys. Lett. 1991, 59, 3482-3484.

LiangFu-2013 Potter, A. C.; Fu, L. Anomalous Supercurrent from Majorana States in Topological Insulator Josephson Junctions. Phys. Rev. B 2013, 88, 121109.

Baxevanis Baxevanis, B., Ostroukh, V. P.; Beenakker, C. W. J. Even-odd Flux Quanta Effect in the Fraunhofer Oscillations of an Edge-Channel Josephson Junction. Phys. Rev. B 2015, 91, 041409.

puddles-Breunig2017 Breunig, O.; Wang, Z.; Taskin, A. A.; Lux, J.; Rosch, A.; Ando, Y. Gigantic Negative Magnetoresistance in the Bulk of a Disordered Topological Insulator. Nat. Commun. 2017, 8, 15545.

CPR-review-Golubov Golubov, A. A.; Kupriyanov, M. Y.; Il Ichev, E. The Current-Phase Relation in Josephson Junctions. Rev. Mod. Phys. 2004, 76, 411.

BSTS-transport-Ando Taskin, A. A.; Ren, Z.; Sasaki, S.; Segawa, K.; Ando, Y. Observation of Dirac Holes and Electrons in a Topological Insulator. Phys. Rev. Lett. 2011, 107, 016801.

Schrade Schrade, C.; Fu, L. Majorana Superconducting Qubit. 2018, arXiv:1803.01002 [cond-mat.mes-hall].

JJ-review-Likharev Likharev, K. K. Superconducting Weak Links. Rev. Mod. Phys. 1979, 51, 101.

Dynes1971 Dynes, R. C.; Fulton, T. A. Supercurrent Density Distribution in Josephson Junctions. Phys. Rev. B 1971, 3, 3015.

Allen2016 Allen, M. T.; Shtanko, O.; Fulga, I. C.; Akhmerov, A. R.; Watanabe, K.; Taniguchi, T.; Jarillo-Herrero, P.; Levitov, L. S.; Yacoby, A. Spatially Resolved Edge Currents and Guided-Wave Electronic States in Graphene. Nat. Phys. 2016, 12, 128.

Zappe1975 Zappe, H. H. Determination of the Current Density Distribution in Josephson Tunnel Junctions. Phys. Rev. B 1975, 11, 2535.

Nesher1997 Nesher, O.; Ribak, E. N. Retrieval of Critical Current Distribution in Small Josephson Junctions. Appl. Phys. Lett. 1997, 71, 1249-1251.

Hui2014 Hui, H.-Y.; Lobos, A. M.; Sau, J. D.; Das Sarma, S. Proximity-induced Superconductivity and Josephson Critical Current in Quantum Spin Hall Systems. Phys. Rev. B 2014, 90, 224517.

Jaynes Jaynes, E.T. Information Theory and Statistical Mechanics. Phys. Rev. 1957, 106, 620-630.













 









 
Supporting Information

 






 
1. SEM images of all devices



The SEM images of all four devices are shown in Figure S1.





 
2. Characteristic parameters of all devices



Characteristic parameters of all four devices are listed in Table S1. Here, R_N denotes the normal-state resistance measured by applying a magnetic field B>B^Al_c or a large excitation current I_ac≫ I_c, while R_N^I-V is the normal-state resistance extracted from I-V curves. Since R_N is obtained with Al in the normal-conducting state one may naively expect that R_N>R_N^I-V. Yet, the opposite is observed, which may be explained by a change of the density of states in the TI due to the proximity effect <cit.>. 





 
3. Characterization of device 2





Several Al-BiSbTeSe_2-Al Josephson junctions were measured down to low temperatures, and all investigated devices showed consistent behaviour. In Figure S2, we present a full characterization of device 2. Essentially all the results of device 1 are reproduced in device 2, except that the E_F of device 2 is much closer to the Dirac point at V_g=0.

Similar to device 1, a two-step transition was seen in the cooling curve of device 2. The sharp transition at T ≈ 1.4K  is associated with the superconducting T_c of Al, and the broad transition occurring subsequently reflects the development of Josephson coupling. Due to the low carrier density at V_g=0, device 2 reaches zero resistance at a lower T compared to device 1 (presented in the main text).
The Dirac point of device 2 is located at V_g≈ -20V (which we define as V_ DP), and the E_F can be tuned deeply beneath the Dirac point by applying a large negative V_g, which enables us to investigate the I_c vs V_g relation in a wide range of the p-type region. The measured I_c(V_g) curve is illustrated by the white band in Figure S2(c). One can see that in the p-type region where V_g<V_DP≈-20 V, the I_c of device 2 actually increases slowly with decreasing V_g.
The I_cR_N product of device 2 plotted in Figure S2(d) shows a behavior similar to that of device 1, although the minimum (at V_g≈ -10 V) is less clear.


 
4. Fraunhofer patterns of device 2 at various V_g



In Figure S3, we present the 2D colour maps of the dV/dI data in the I_dc vs B plane measured in devices 2. Irregular Fraunhofer patterns with suppressed or missing lobes were observed in the entire range of the gate voltage. As E_F approaches the Dirac point, the Fraunhofer pattern becomes blurred, and the oscillations of I_c are only barely visible. This is presumably because the Josephson coupling is weakened at the Dirac point due to a low carrier density, and consequently the zero-resistance state is more prone to the thermal perturbation caused by electromagnetic noise.





 
5. Estimation of Stewart-McCumber parameter



Figure S4 shows the I-V curves of device 2 and device 3 measured in both forward and backward sweep directions. It is clear that the dc I-V characteristics are non-hysteretic. In the RCSJ model, such behaviour is expected for overdamped Josephson junctions in which the Stewart-McCumber parameter β_c (≡ 2π I_c (R_N^I-V)^2 C / Φ_0) is much smaller than unity (C is the shunt capacitance of the junction).
To calculate the β_c of our junctions, we estimated the shunt capacitance C by considering the two capacitors forming between the two Pt/Al/Au contact pads and the conductive Si substrate via the 290-nm-thick SiO_2 dielectric layer. It should be noted that these two capacitors are connected in series. Using the dielectric constant ϵ_r = 3.9 for SiO_2 and the areas of the contact pads of the devices, we  estimate C ≈ 10-30 fF. The shunt capacitance originating from the direct capacitive coupling between the two closely spaced Al electrodes at the junction is estimated to be roughly three orders of magnitude smaller than this C and hence is safely neglected.

Using the estimated values of C, we obtain β_c ≈ 0.2-2.5 in different
devices (see Table S1). In the framework of the RCSJ-model, these junctions are at the boundary between the underdamped (β_c ≫ 1) and overdamped (β_c ≪ 1) regimes. In the literature, a previous study by Goldhaber-Gordon's group reported TI-based Josephson junctions in this regime, and a non-hysteretic behaviour similar to ours was observed <cit.>. The numerical calculation for junctions in this regime also gave similar non-hysteretic results <cit.>. Therefore, the absence of
hysteresis in the dc I-V characteristics of our devices is consistent with their estimated β_c values.



 
6. Estimation of the transparency of Al-BiSbTeSe_2 contacts



 The transparency T_J of our Al-BiSbTeSe_2 superconducting contacts were estimated from the values of I_e and R_N^I-V by using the Octavio-Tinkham-Blonder-Klapwijk theory <cit.>. Here, I_e and R_N^I-V are the excess current and the normal-state resistance extracted from the I-V characteristics of the Josephson junctions, as defined in the main text. To identify I_e, the dc I-V curves were measured to a high bias beyond V_dc=2Δ_Al/e, where 2 Δ_Al≈ 320 μeV is the superconducting gap of Al. The linear part of the high-bias I-V curve was extrapolated to V_dc=0 and the intercept on the I_dc axis gives I_e. The R_N^I-V value was calculated as the inverse slope at high bias (i.e. V_dc>2Δ_Al/e). The value of Δ_Al was obtained from the multiple Andreev reflection peaks as discussed in the main text. In Ref. , the normalized excess current I_eR_N^I-V/Δ was numerically calculated as a function of the barrier strength Z, which is directly related to the transparency via T_J=1/(1+Z^2). Using the given relation between I_eR_N^I-V/Δ and Z, we estimated T_J≈ 0.65-0.86 in our devices, as listed in Table S1.






 
7. Estimation of V_CB



Considering the contributions from both top and bottom surfaces, the total two-dimensional carrier density n_2d coming from the topological surface states in a ultra-thin sample is approximately given by

    n_2d(k)=2 ·1/(2π)^2·π k^2=k^2/2π.

According to the ARPES results <cit.>, the conduction band bottom of BiSbTeSe_2 is located at E_CB≈0.2 eV (measured from the Dirac point), and the Fermi wave vector of the surface states at this energy, k_CB, is about 0.076 Å^-1. This gives n_2d(k_CB) = 9.1 × 10^12 cm^-2.

The gate voltage required to tune E_F from the Dirac point to the conduction band bottom is given by (V_CB-V_DP)=en_2d/C_g. Here, C_g=ϵ_0ϵ_r/d ≈ 12 nF/cm^2 is the capacitance per unit area of the back gate, with ϵ_r=3.9 and d=290 nm being the dielectric constant and the thickness of the SiO_2 dielectric layer, respectively. With these values, we obtain V_CB-V_DP≈ 121 V. For V_DP=-88 V in device 1 shown in the main text, one obtains V_CB≈33 V.


 
8. Applications of the maximum entropy fitting method



The presented anomalous Fraunhofer patterns have been analyzed using an advanced fitting scheme which is based on the maximum entropy principle. As a sanity check we applied the algorithm to an ideal pattern of a junction with a homogeneous current distribution given by I_c(B)=I_c0| sin(π x)/π x|, where x=Φ_J/Φ_0 and Φ_J is the flux through the junction. First, we initialized the current distribution to random numbers and then applied the same algorithm used to describe the anomalous patterns in the main text. Both the good fit to I_c(B) as well as the flatness of the reconstructed current distribution (Figure S5(a)) show that the fluctuating current distribution extracted in the main text is neither due to over-fitting of the data, nor an artifact, but the most physical distribution explaining the experimental data.

We also applied the method to published data of SQUID oscillations measured from topological insulator Josephson junctions <cit.> which have been argued by the authors to contain signatures of majorana bound states in the form of lifted SQUID nodes. To model the SQUID oscillations using our single junction fitting scheme, we forced the current distribution j(0<x<W) to be zero in the interval [w, W-w]. Here, w is the width of the individual junctions and (W-2w)· L_eff=A_eff represents the effective SQUID area with the effective length L_eff of the junctions. From the fit shown in Figure S5(b), we obtain A_eff=8.5 μm^2, W=13.3 μ m, L_eff=765 nm and w=1.1 μ m. These values are consistent with the observed oscillation period and with the device geometry, taking into account the effect of magnetic flux focusing. The reconstructed current density j(x) (inset of Figure S5(b)) shows moderate, yet finite spatial fluctuations whose amplitude and frequency are comparable to those discussed for our device in the main text. While the overall shape of the SQUID oscillations is well captured by our modeling, the lifting of the critical current close to the envelope's first node (marked by an arrow) cannot be fully accounted for. Even the tuning of the algorithm from maximum entropy towards least-square fitting does not improve the fit in this data range, which opens a room for discussing more involved mechanisms for the observed node lifting.



10
Gueron-2012 Guéron, S.; Pothier, H.; Birge, N. O.; Esteve, D.; Devoret, M. H. Unconventional Josephson effect in hybrid superconductor-topological insulator devices. Phys. Rev. Lett. 2012, 77, 3025.

Williams-PRL-2012 Williams, J. R.; et al. Unconventional Josephson effect in hybrid superconductor-topological insulator devices. Phys. Rev. Lett. 2012, 109, 56803.

JJ-book-Likharev Likharev, K. K. Dynamics of Josephson junctions and circuits (Gordon and Breach Science Publishers, 1986).

OTBK-1983 Octavio, M.; et al. Subharmonic energy-gap structure in superconducting constrictions. Phys. Rev. B 1983, 27, 6739.

OTBK-1988 Flensberg, K.; Hansen, J. B.; Octavio, M. Subharmonic energy-gap structure in superconducting weak links. Phys. Rev. B 1988, 38, 8707.

BSTS2-2012 Arakane, T.; et al. Tunable Dirac cone in the topological insulator Bi_2-xSb_xTe_3-ySe_y. Nat. Comm. 2012, 3, 636.

Kurter-2015 Kurter, C.; Finck, A. D. K.; Hor, Y. S.; van Harlingen, D. J. Evidence for an anomalous current–phase relation in topological insulator Josephson junctions. Nat. Comm. 2015, 6, 7130.






