Evidence of Distributed Robust Surface Current Flow in 3 D Topological Insulators

The topologically protected conducting state is expected to exist on the entire surface of threedimensional topological insulators (TIs). Using concurrent measurements of the local and nonlocal conduction, we provide experimental evidence for the topological robustness of the surface-conducting states of bulk-insulating Bi1.5Sb0.5Te1.7Se1.3 crystalline flakes. The detailed investigation of local and nonlocal charge conductance on the top surfaces, combining with the comprehensive numerical simulation, reveals that the charge current is widely distributed over the entire surface of a TI. Our findings show evidence of the presence of the topologically protected conducting state at the side wall with irregularly stacked edges between the top and bottom surfaces. This study provides a reliable means of accurately characterizing the topological surface states with inherent nonlocal surface-dominant conducting channels in a TI.


I. INTRODUCTION
Topological insulators (TIs) show new quantum states of matter whose surface supports metallic conducting channels, even with the Fermi level (E F ) lying in the bulk energy gap [1][2][3][4][5].The surface-conducting channels of a TI are topologically protected, as long as the time-reversal symmetry is preserved [6,7].This is in sharp contrast to ordinary band insulators whose "accidental" surface state is not protected by physics.In addition, the nondegenerate helical spin texture of the topological surface-conducting state (TSS) [8,9] has attracted much attention, motivated by theoretical predictions of the emergence of exotic phenomena in the TSS, such as the magnetic monopole [10], topological exciton condensation [11], Majorana fermions [12], the quantum anomalous Hall effect [13], etc.The conduction nature of the TSS has been studied by tuning the E F well inside the bulk energy gap, which is essential for reducing the bulk conduction that hampers the observation of the phenomena specific to the TSS.Despite many efforts, however, it is still challenging to confirm the topologically protected conduction on the surface of a three-dimensional (3D) TI, with negligible conduction in the bulk [14][15][16][17][18][19][20][21][22][23].In particular, the existence of the robust surface-conducting state on the rough sides of a TI crystal, which has not been directly confirmed by transport measurements, is of great interest.
In general, the surface-conducting channels of a TI are admixed with the conducting channels formed in the bulk via the thermal activation among impurity levels, causing difficulties with accurately characterizing the surface conduction.Even with sufficiently suppressed bulk conduction, the existence of the robust surface conduction itself and the consequent distributed current flow on a TI makes the accurate estimation of the surface conductance and its precise characterization complicated and challenging.For the dominant surface conduction with a reduced bulk contribution, the current bias between the source and drain, arranged on the top surface of a TI, is not confined to the region between the two electrodes but flows over the entire TI surface.Because of this partial leakage of the bias current out of the source-drain region, the surface conductance cannot simply be determined in terms of the local source-drain potential difference (V loc ) and the bias current (I bias ), as I bias =V loc .Thus, in particular, a significant error can arise in estimating the backgate-voltage dependence of the top surface conductance, which is hardly affected by the backgating, as the backgating can significantly affect the current-flow profile on the top surface.Despite the issue of this distributed current flow on a TI surface, it has never been taken into account in estimating the surface conductance of a TI.This may have been the cause of the inconsistency reported previously in transport properties of TIs, such as the wide sample dependence of the magnetoresistances (MR) [15,17,[22][23][24][25][26][27][28][29][30] and the lack of consensus on the value of the Lande g factor of electrons in the TSS [15,25,[31][32][33][34][35].
In this study, we resolved the issue by concurrently monitoring the local (V loc ) and nonlocal (V nloc1 , V nloc2 ) potential differences in Bi 1.5 Sb 0.5 Te 1.7 Se 1.3 (BSTS) single crystalline flakes and matching them to numerical simulation [36] based on the finite-element method.With E F residing inside the bulk energy gap, owing to the electronhole carrier compensation in the material [37], we observed appreciable nonlocal signals.It is possible only in the presence of a robust conducting state, i.e., a TSS, over the entire surface of the crystal, including the sides, which are often rugged or rough with irregularly stacked edges between the top and bottom surfaces (see Appendix A).This distributed current flow, with a significant amount of current leaking into the bottom surface, suggests that the electrical surface conductance of a TI cannot be determined from local measurements (V loc ) only.This study provides a direct confirmation of the presence of the TSS, even on the rough sides of a TI, and demonstrates that the surface conductance can be determined accurately only with concurrent measurements and a comprehensive analysis of both the local and nonlocal potential differences on the top surface.

II. EXPERIMENTAL METHOD
Single crystals of BSTS were grown using the self-flux method [25,38].Stoichiometric mixtures of high-purity starting materials [Bi(5N), Sb(5N), Te(5N), Se(5N)] were loaded in an evacuated quartz ampoule, which was then heated up to 850 °C.After annealing at 850 °C for two days to enhance the material homogeneity, the melt mixture was slowly cooled down to 500 °C at a rate of 2 °C=h for a week.Before complete furnace cooling, it was post annealed at 600 °C for one more week to further improve the crystallinity.The stoichiometry and the high crystallinity of the single crystals were confirmed by the energy dispersive spectroscopy and the x-ray diffraction, respectively.Recent scanning tunneling microscopy and spectroscopy on our BSTS crystals showed that their surface states have linear dispersion and chirality from spin texture, thus verifying the topological nature.This implies that topological surface states remain intact even with alloying of constituent atoms and significant surface disorder [39].
For transport measurements, a BSTS flake was mechanically exfoliated onto a heavily doped Si substrate, capped with a 300-nm-thick SiO 2 layer [see Figs.1(a) and 1(b)], which allowed gate tunability of the E F at the bottom surface of the BSTS flake.To prevent electrical shorting at the sides of the BSTS flake, we electrically insulated the side surface by covering the edge of the flake with crosslinked poly(methyl methacrylate) (PMMA) 950 k C4 [inset of Fig. 1(a)] before thermally depositing the electrodes on the top.Electrical contacts were made using a standard electron (e)-beam lithography process followed by e-gun evaporation of Ti/Au (10 nm=300 nm) electrodes.The transport measurement was conducted via low-frequency lock-in techniques with a r.m.s.bias current of 100 nA at 17.77 Hz.The details of the sample dimensions are summarized in Table I.

Au ele ctr od e
Cro sslink ed PM MA

Physical parameters Length
Thickness of flake 110 nm Width of flake 5.5 μm Width of electrodes 0.8 μm Thickness of electrodes 310 nm Spacing between electrodes (edge to edge) 1.9 μm

A. Observation of the nonlocal potential difference
Figure 1(a) shows the scanning electron microscopy (SEM) image of a device.We measured V loc , V nloc1 , and V nloc2 concurrently, using the configuration shown in Fig. 1(b).The (þ) and (−) signs in Fig. 1(b) indicate the polarities of the voltage probes.Figure 1(c) shows the temperature (T) dependencies of V loc , V nloc1 , and V nloc2 .Appreciable nonlocal signals (V nloc1 and V nloc2 ) were observed; V nloc1 was evident, even at room temperature [V nloc1 ð300 KÞ ≈ 0.25 μV].Similar behavior was also observed in another BSTS flake of this study.In a standard Hall-bar geometry, nonlocal resistance can appear because of the quantum spin Hall effect [40], the spin Hall effect [41], or the quasiballistic transport mechanism [42].However, such causes are not applicable to our device without Hall-bar contacts [see Fig. 1(b)].The most plausible cause of the observed nonlocal signals, V nloc1 and V nloc2 , is the presence of conducting channels over the entire surface, including both the side and bottom surfaces, of the BSTS flake, which is more prominent when the conduction through the bulk is negligible.This feature is unique to the 3D TI but has not been directly confirmed by transport measurements.The simulation results, compared with the measured V bg and magnetic-field dependencies of V loc , V nloc1 , and V nloc2 , indicate that the observed finite nonlocal signals are a direct consequence of the current distribution, which was not confined between the source and drain but was widespread over the entire surface of the BSTS flake, with negligible current flow in its bulk.In addition, our experimental and simulation results commonly showed that V nloc1 and V nloc2 had similar V bg and magnetic-field dependencies, although they exhibited slight differences in their magnitudes.

B. Numerical simulation method
We performed numerical simulations on the current distribution, based on the finite-element method, using the commercial software package COMSOL [36].Details of the validity of the numerical simulation are described in Appendix C. For the numerical simulation, we used the actual shape and dimensions of the device (see Table I), which were determined from atomic force microscopy and SEM.The sheet conductance of the top (G top □ ) and bottom (G bot □ ) surfaces depends on the extent of the band bending at the two surfaces [8,43,44].In general, the values of G top □ and G bot □ are different.Therefore, numerical simulations were performed for different combinations of G top □ and G bot □ , with each ranging from 0.3 e 2 =h to 20 e 2 =h.
We also examined the influence of the bulk conductivity σ bulk on the values of V loc , V nloc1 , and V nloc2 .Numerically simulated values of V nloc1 , in Fig. 2(b), almost vanish for σ bulk ≳ 10 S=m in our BSTS flake.Figure 2 also shows that V loc and V nloc1 are almost independent of σ bulk when σ bulk < 10 −2 S=m for all combinations of G top □ and G bot □ , indicating that, in this case, the surface conduction predominates over the bulk conduction.
In our simulation, σ bulk is adopted in the range of 10 −3 ∼ 10 S=m.See Appendix C for details of the numerical simulation method.

C. V bg dependence of local and nonlocal signals
Figure 3(a) illustrates our simulation model.A bias current (I bias ¼ 100 nA) was injected from the source to the drain.S loc and S nloc1 represent the cross sections where the total currents in the respective local and nonlocal (V nloc1 ) regions are calculated.We also defined the horizontal bulk cross section (S bulk ) in addition to the top (S top ) and bottom (S bot ) surfaces.Figure 3(b) shows how transport parameters are extracted by matching the experimental data to the simulation results.We calculated a set of (V loc , V nloc1 , V nloc2 ) for each input parameter set of (σ bulk , G top □ , G bot □ ).A few million sets of (V loc , V nloc1 , V nloc2 ) for different combinations of input parameters σ bulk , G top □ , and G bot □ were prepared for matching between the data and the simulation results within uncertainty of jΔV loc j ≤ 1 μV, jΔV nloc1 j ≤ 0.3 μV, and jΔV nloc2 j ≤ 0.05 μV, which gave rise to an optimum value set of (σ bulk , G top □ , G bot □ ) corresponding to each data set.In this manner, from the V bg dependence of V loc , V nloc1 , and V nloc2 [Fig.3 3(d)].The measured V bg dependence of V loc in Fig. 3(c) has an opposite response to that of V nloc1 (as well as V nloc2 ) at 4.2 K.As the sheet resistance of the bottom surface (R bot □ ) increases, the current (I top ) on the top surface at S loc is expected to increase, resulting in an increase (decrease) of V loc (V nloc1 and V nloc2 ).By matching the simulation to the measured results, we confirmed that this qualitative reasoning is correct and that the opposite variation of V loc and V nloc1 with V bg in Fig. 3(c) is due to the ambipolar nature of the conduction in the bottom surface.
Figure 3 , where L is the center-to-center distance of the two voltage probes (6.0,2.5) 2. Numerical simulation results of the effect of σ bulk on (a) V loc and (b) of V loc and W is the width of our flake.This V bg dependence of R bot □ clearly shows the ambipolar nature of the bottom-surface conductance arising from a Diracfermionic dispersion relation, a distinguishing characteristic of the TSS, with the Dirac point at V bg ¼ −10 V.It is reminiscent of the V bg dependence of the sheet resistance in graphene, showing a maximum at the Dirac point.
However, as shown in Fig. 3(d), R top □ remains almost insensitive to backgating.It is apparent that, with a thickness of about 110 nm for our BSTS flake, the effect of V bg does not reach the top surface because of the screening at the bottom [43,44].This is consistent with the results of our previous studies [45].The resistance peak of R bot □ in Fig. 3(d) has a value of around 2.2 e 2 =h.This value is about 1=4 the minimum conductivity of graphene in the diffusive regime, with relatively low mobility ð3000-8000 cm 2 =VsÞ [46], which may be consistent with the fact that the TSS has 1=4 the degeneracy of graphene [47].Thus, in summary, since R bot □ increases as V bg approaches the Dirac point, the currents at S loc on the bottom surface and at S nloc1 on the top surface decrease simultaneously, while the current at S loc on the top surface increases.Strikingly, the V bg dependence of V loc and V nloc1 in Fig. 3(c) is a consequence of variations in the current flow at S loc and S nloc1 on the top surface, rather than the change in resistance on the top surface.
R exp □ , shown in Fig. 3(d) (triangles), is obtained from the observed value of V loc only, which is the conventional way of estimating the sheet resistance of samples.The relative discrepancy, ðjR exp □ − R top □ Þj=R top □ , between the "apparent" sample resistance taken on the top surface (R exp □ ) and the true value of the top-surface resistance (R top □ ), reached 30%-40%.A large discrepancy between R exp □ and R top □ (as well as This discrepancy is mainly caused by the distributed current flow over the entire TI surface and the difference in the values of G top □ and G bot □ , together with their V bg dependencies.Thus, the apparent conductance determined from the observed V loc , without a comprehensive analysis including information on V nloc1 and V nloc2 , will result in an erroneous estimation of the surface conductance and its gate dependence. We also obtained bulk conductivity, which turned out to be highly sensitive to small fluctuations of V nloc2 .Thus, only the upper limit of σ bulk was estimated meaningfully, at around 0.6 S=m.This value of σ bulk is a few tens of times smaller than the ones reported previously for extremely thick BSTS crystals (thickness of around 100 μm) [25,38,48].Because the metallic surface-conducting channels always coexist with thermally activated bulkconducting channels in a 3D TI, the measured resistance represents the combined resistance of these conducting channels connected in parallel.The surface conductance amounts to at least a few S=m [see Eq. (B4) in Appendix B].Thus, if the measured bulk conductivity (σ exp bulk ) in thick crystals decreases below about 10 S=m, similar to the values obtained in previous reports on BSTS [25,35,38,[49][50][51][52][53], σ exp bulk can no longer represent the true bulk conductivity (σ bulk ) of a 3D TI [see Appendix B for a discussion on the discrepancy between the measured apparent value (σ exp bulk ) and the true value (σ bulk ) of bulk conductivity in a 3D TI].For this reason, we suspect that the previously reported values of σ exp bulk of BSTS may have been significantly overestimated.
Simulation model (the thickness of BSTS is exaggerated for visual clarity).I bias ¼ 100 nA is injected from the source to the drain.S loc (S nloc1 ) represents the hypothetical plane used to extract the total current through the cross-sectional area between the two voltage probes of V loc (V nloc1 ).(b) Schematic diagram for methodology of this study, illustrating how to determine an optimal set of values for (σ bulk , G top □ , G bot □ ) by matching an observed set of (V loc , V nloc1 , V nloc2 ) to the simulated results.(c) V bg dependence of V loc (black squares) and V nloc1 (magenta circles) taken at 4.2 K.The inset shows Þ is also plotted for comparison, where L is the center-to-center distance of the two voltage probes of V loc and W is the width of our flake.
We also show that the different T dependencies among V loc and V nloc1;2 in Fig. 1(c), in the temperature range below about 70 K was caused by the combined T dependencies of G top □ , G bot □ , and σ bulk (see Appendix D).I nloc1 (green squares) in Fig. 4(e) shows the nonlocal current on the top surface at S nloc1 .In contrast to I bulk , an appreciable amount of current appears on the top surface at S nloc1 , with its direction opposite to that in the local-flow region.The large value of V nloc1 in Fig. 1(c) results from this I nloc1 , which cannot be realized without the robust surface-conducting channels in our BSTS flake, in particular, at the side surfaces, a feature of the TI that is in sharp contrast to an ordinary band insulator.In fact, even for an ordinary band insulator, a surface-conducting state can arise from accidental causes (e.g., surface impurities, surface reconstruction, surface band bending, or band-gap narrowing near the surface).However, existence of these "accidental surface states" are not topologically robust.Thus, applying V bg or, equivalently, shifting the E F can switch the accidental surface states on or off.However, as shown in Fig. 3(c), large values of V nloc1 are observed for all V bg , even with the ambipolar transport at the bottom surface in Fig. 3(d), which indicates the robustness of the surface-conducting state in our thin BSTS flake.In addition, if the top and bottom surfaces have opposite carrier types in an ordinary band insulator, a p-n junction would form somewhere between the top and bottom surfaces and should suppress much of the nonlocal surface current, I nloc1 .However, the observed nonvanishing V nloc1 in Fig. 3(c), irrespective of V bg , indicates that no such large conductance suppression develops on the sides of our BSTS flake (see Appendix A).Our observation of the robust and large V nloc1 cannot be explained in the absence of a gapless surface-conducting state, i.e., a Diracfermionic TSS.Thus, the side surfaces act as a filter that only passes the topologically protected surface currents contributing to the nonlocal potential difference.The current density profile at each plane, which is defined as in Fig. 3(a), is plotted using the same color scale for ðG top □ ; G bot □ Þ ¼ ð6.0; 2.5Þðe 2 =hÞ and σ bulk ¼ 0.06 S=m, which correspond to the observed values of V loc ¼ 146 μV, V nloc1 ¼ 2.72 μV, and V nloc2 ¼ 0.58 μV at 4.2 K for V bg ¼ 0 V.The solid lines in (a), (b), and (c) are contour lines describing the current flow (i.e., these lines are perpendicular to the equipotential lines).(d) V bg dependence of I top (red squares) [I bot (blue circles)], the current on the top (bottom) surface at the cross section of S loc depicted in Fig. 3(a).(e) V bg dependence of I side (black circles) [I bulk (magenta triangles)], the current on the side surface (inside bulk) at the cross section S loc depicted in Fig. 3(a).I nloc1 is the nonlocal current on the top surface at S nloc1 depicted in Fig. 3(a).

E. Magnetic-field dependence of local and nonlocal signals
The response of V loc to an external magnetic field in Fig. 5(a) shows a typical weak-antilocalization (WAL) behavior, which has been extensively studied previously [16,18,19,21,23,[55][56][57].In contrast to V loc , however, V nloc1;2 shows both positive and negative MR, depending on V bg .V nloc1;2 is affected by two factors: the change in the total current flowing through S nloc1 shown in Fig. 3(a) and the change in the resistance of R top .With an increasing magnetic field, R top also increases according to the WAL effect.If the current profile is not modified, then V nloc1;2 would increase with the magnetic field.However, because R bot also increases along with R top , with increasing B, V nloc1;2 depends on the consequent redistribution of current with V bg , which is determined by the variation in the set of ðσ bulk ; G top □ ; G bot □ Þ.By matching the observed results in Figs.5(a)-5(c) to the numerical simulation, we find that the complex magnetic-field dependence of V nloc1;2 in Figs.5(b) and 5(c) could also be attributed to the combined WAL effect of the top and bottom surfaces with different zero-field conductance [see Figs.Because the TSS belongs to the symplectic class, i.e., τ e , τ so ≪ τ ϕ with 1=τ s ¼ 0 (τ e : elastic scattering time; τ so : spin-orbit scattering time; τ ϕ : phase relaxation time; τ s : spin scattering time), the magnetoconductance correction of the TSS is expressed as where ψ is the digamma function, e is the electron charge, ℏ is Planck's constant divided by 2π, and l ϕ is the phase relaxation length [58].The theoretical value of α is −0.5 for each two-dimensional (2D) conducting channel for the symplectic class.The solid lines, in Figs.5(d) and 5(e), represent the best fits to Eq. ( 1), in good agreement with WAL behavior.Figures 5(f) and 5(g) show the V bg dependence of α and l ϕ , obtained from the best fits in Figs.5(d) and 5(e).The value of α on the bottom surface (circles) is close to −0.5, which agrees well with the prediction for a single conducting channel.The relatively large deviation of α from −0.5 near the Dirac point on the bottom surface can be understood by the presence of spatially nonuniform carrier density and carrier types (e-h puddles), which have been observed on the surface of a TI [59], similar to graphene [60].
In contrast to the bottom surface, the value of α on the top surface (squares) is close to −1 for all V bg .At the surface of a TI, a 2D electron gas (2DEG) is formed in a quantum well set up by band bending.The depth and magnitude of the band bending strongly depend on the surface impurity, carrier and impurity density, temperature, and external potential [61,62].In particular, the 2DEG that forms near the surface of a TI exhibits a large Rashba spin splitting due to strong spin-orbit coupling [8,43,44,63,64].Because the Rashba spin-split band also has a spinmomentum-locked helical spin texture, it exhibits a WAL effect with −1 < α < 0, depending on the magnitude of spin splitting or the interband scattering strength [45,58,65,66].Therefore, on the top surface, not only the TSS but also the Rashba-split 2DEG contributes to the WAL effect, giving rise to the total value of α close to −1.It is believed that the difference in the values of α on The solid lines in (d) and (e) are the best fits to Eq. ( 1).Each set of data, together with the fitting curve, is shifted vertically for clarity.(f),(g) V bg dependence of α and l ϕ obtained from the best fits to Eq. ( 1).Squares (circles) correspond to the top (bottom) surface.
the top and bottom surfaces is caused by doping of the top introduced during electrode patterning and subsequent processes.This result is consistent with our previous observation in BSTS [45].In Fig. 5(g), l ϕ on the top surface is fairly insensitive to V bg , in contrast to l ϕ on the bottom surface, which follows the sensitivity of R top □ and R bot □ on V bg , as shown in Fig. 3(d).The well-behaving magnetoconductance on the top and bottom surfaces, with highly relevant WAL parameters, extracted out of seemingly contradicting V loc ðBÞ and V nloc1;2 ðBÞ on the top surface, conclusively justifies our analysis scheme for the surface conductance in our BSTS TI crystal flakes.

IV. DISCUSSION
Most of the transport measurements on TI thin films or flakes have been carried out using electrodes arranged on the top surface.The results shown in our study point out that, if the conduction in a 3D TI becomes dominated by the surface conduction, the MR obtained from the top surface should be closely analyzed to extract the correct conduction parameters.The magnetic-field-dependent oscillation of V loc can be caused not only by the MR oscillation of the top surface but also by the field-dependent oscillation of I top that is affected by the MR oscillation at the bottom surface.This can lead to multiple-field periodicity in the fast Fourier transform of R vs 1 B .In fact, this multiple-field periodicity in Shubnikov-de Haas oscillations has been reported previously [67][68][69].However, even if the field periodicity corresponding to the top and bottom surfaces is resolved, one still cannot extract various valuable physical parameters, such as the cyclotron mass, mean-free time, and mobility, without precise determination of ΔR top and ΔR bot .To that end, V nloc1;2 , as well as V loc , should be measured using a measurement configuration similar to the one used in this study.The above argument is valid even when the current-biasing electrodes are arranged at the two ends of an elongated TI crystal in contact with its sides, as long as the top and bottom surfaces have substantially different sheet conductances or magnetoconductances.
In transport measurements, a conductance estimation using I bias =V loc is valid only when all of the injected bias current (I bias ) flows locally and uniformly across the two voltage probes.In the presence of both multiple conducting channels and subsequent distributed current flow on the surface of a 3D TI, the conductance cannot simply be estimated from the relation I bias =V loc .This issue of the distributed current flow becomes more serious when the surface conduction predominates over the bulk conduction.In our previous studies [45], which included only local measurements on a BSTS flake, grown under conditions and thicknesses similar to the ones in this study, we obtained around 6% of bulk conductance at 4.2 K. Comparing this with σ bulk obtained at 4.2 K in this study, we believe that the previous estimation was highly overestimated.

V. CONCLUSION
In conclusion, we report on the investigation of the distributed current flow over the entire surface of a 3D TI, determined by measuring both V loc and V nloc1;2 concurrently.With the topologically robust conducting state on the surface of a TI, the current flow is bound to be distributed all over the surface of a TI, which makes the estimation of the surface conductance based on the local potential measurements readily erroneous.We showed that the surface conductance of the top and bottom surfaces can be estimated with high precision by matching both the local and nonlocal transport measurements on the top surface to the numerical simulation, based on the finite-element method.The backgate-voltage, temperature, and B-field dependencies of the conductance on the top and bottom surfaces revealed the high physical relevance.This study offers a highly reliable means of accurately determining the surface conductance of a TSS on a 3D TI, which enables the accurate characterization of the topological conducting state existing on the surface of a TI.
Figure 6 shows magnified SEM images for different sides of the Bi 1.5 Sb 0.5 Te 1.7 Se 1.3 (BSTS) flake used in this study.Because the BSTS is a layered compound, as described in the main text, the side surfaces are rugged or rough when a crystal is exfoliated into a flake.In ordinary band insulators, surface-conducting states that may exist in side surfaces cannot be robust.However, in a three-dimensional topological insulator (TI), a gapless boundary conducting state should exist on any surface because of topological protection.Our observation of large nonlocal potential differences in our TI flake a direct consequence of the existence of this robust gapless surfaceconducting state on the BSTS TI material.

APPENDIX B: ESTIMATING BULK CONDUCTIVITY OF TI
The value of the base-temperature (4.2 K) bulk conductivity we obtained in this study turned out to be a few tens of times smaller than the values reported previously for thick BSTS single crystals [25,38,48].In the main text, we stated that σ bulk in the previous studies could have been overestimated because of the coexistence of metallic surface-conducting channels and thermally activated bulk channels.Here, we provide more discussions on that argument.In an ordinary semiconductor or a metal, the bulk conductivity is estimated from the relation where I in is the bias current, L is the spacing between two voltage probes, V loc is the observed potential difference, and t and W are the thickness and width of a sample, respectively.In fact, σ exp bulk , obtained from Eq. (B1), is the value averaged over all local values of the bulk conductivity, σ bulk .Therefore, estimating σ bulk of a material using σ exp bulk is valid only when the local fluctuation of σ bulk or relative difference of σ bulk between the surface and the bulk is sufficiently small.Estimating σ bulk using relation (B1) is valid for most materials.However, if multiple conducting channels are present and combined in parallel as in 3D TIs, the validity of relation (B1) is no longer guaranteed.
Conducting channels in a bulk TI can be divided into bulk and surface contributions.Separating the latter, the apparent value of bulk conductivity in relation (B1) is modified as where G 2D is the surface sheet conductance and t is the thickness of the sample.If σ bulk is sufficiently large so as to neglect the surface contribution, σ exp bulk ≈ σ bulk .This corresponds to ordinary semiconductors and conductors, with relatively low resistivity of about mΩ cm.
For a BSTS bulk crystal of thickness t ¼ 100 μm, with the Fermi level (E F ) of the surface state located around the middle between the bottom of the conduction band and the Dirac point, the carrier density is about n s ≈ 2 × 10 12 cm −2 (note that the maximum carrier density, without touching the conduction band, is about n s ≈ 5 × 10 12 cm −2 ).Then, assuming the carrier mobility to be around 1000 cm 2 =ðV • sÞ (these assumptions are reasonable, and a variation of the values does not affect the result of our discussion), along with the relation σ ¼ neμ, the sheet conductance of each surface becomes around 300 μS (corresponds to about 8.3 e 2 =h).For the top and bottom surfaces together, G 2D ≈ 600 μS.Thus, Eq. (B2) leads to If σ exp bulk is comparable to around 10 S=m, as in previous reports on BSTS [25,35,38,[49][50][51][52][53], σ exp bulk no longer represents the true value, σ bulk .Thus, in the presence of the separate surface-conducting channels in a TI with appreciable surface conductance, the apparent value of the bulk conductivity, σ exp bulk , can be significantly larger than the true value, σ bulk .
In addition to the topological surface-conducting channels, the conduction by the 2D electron gas formed in the quantum well due to the surface band bending also enhances the surface conductance G 2D .The surface band bending is a common feature at the surface of a semiconductor, especially in a narrow-band-gap semiconductor like TIs.Thus, in a thick bulk TI crystal, multiple surfaceconducting channels, other than the bulk-conducting channels, contribute to G 2D and enhance the apparent bulk conductivity, σ exp bulk , in Eq. (B1).In thin flakes, as in our samples, the dominant surface-conducting channels, connected in parallel to the bulk channels, make the estimation of the bulk conductivity impossible.In this study, by resorting to measurements of the nonlocal potential difference and the numerical simulation, we were able to set the upper limit value of σ bulk to about 0.6 S=m at 4.2 K.For the reasons discussed above, we suspect that the values from previous studies were overestimated.

APPENDIX C: DETAILS AND VALIDITY OF NUMERICAL SIMULATION METHOD
We used the commercial software package (COMSOL) [36] to calculate the current distribution on the surface of a 3D TI, based on the finite-element method.Figure 7(a) shows the top view of our simulation model.For the numerical simulation, we used the actual sample dimensions (see Table I), which were obtained from the atomic force microscopy and SEM.To exclude the unintended electrical contact at the sides of our BSTS flake, we covered the side edges underneath the Ti/Au electrodes with crosslinked PMMA 950 K C4 as an insulating layer.All the edge-to-edge spacings between two neighboring electrodes are identically 1.9 μm. Figure 7(b) is the cross section of the simulation model, and Fig. 7(c) shows a magnified view.We set the thickness of the surface-conducting channel (t s ) to be 5 nm in the simulation [see Fig. 7(c)].In fact, the spatial distribution of jψj 2 [ψ is the wave function of the TSS in a TI] is at most 1-2 quintuple layers (QL), which corresponds to 1-2 nm, depending on the van der Waals gap, i.e., the spacing between adjacent QLs [70][71][72].In addition to the TSS, ordinary conducting channels of 2DEG form on the surface of a TI.The range of jψj 2 of electrons a 2DEG can be up to nm, depending on the position of E F and the impurity density [8,43,44].Because carriers in both the TSS and the 2DEG have a common 2D nature, the accurate value of t s in numerical simulations is not crucial.The sheet conductances of the top and bottom surfaces (G top □ and G bot □ ) matter the most.In fact, there is little difference between the simulation models for different t s .Figures 7(d 7(b), because the aspect ratio of the cross section in our flake is very large, t s is hardly perceived; then, the difference between simulation models with different t s is almost negligible.
Figure 8 shows the simulation models for different t s (20 nm, 10 nm, 5 nm, 4 nm, 3 nm) with the same G top □ and G bot □ .The effect of t s variation on the simulation result of V loc and V nloc1 is very small.Even the small difference between t s ¼ 20 nm and t s ¼ 10 nm almost disappears as t s decreases below 5 nm [see the inset in Fig. 8(b)].Therefore, we conclude that the simulation results for the model of t s ¼ 1-2 nm should be almost the same as our model of t s ¼ 5 nm.Lastly, in the numerical simulation, we set G side □ equal to G top □ .Because the area of the side surface is very small compared to that of the top and bottom surfaces, the effect of a small ambiguity in the value of G side □ hardly affects the results of a set of (σ bulk ,G top □ , G bot □ ) unless the G side □ is significantly smaller than G top □ .

APPENDIX D: T DEPENDENCE OF LOCAL AND NONLOCAL SIGNALS
The T dependence of V nloc1 and V nloc2 in Fig. 1(c) in the main text is in discordance with that of V loc .V loc continues to decrease below around 80 K, while V nloc1 and V nloc2 slightly increase and saturate, respectively.This seems to be somewhat inconsistent with the behavior of metallic surfaces.Figures 9(a □ and G bot □ for V bg ¼ 0 V, respectively [the inset of Fig. 9(b) shows σ bulk vs T], which were obtained from the values of V loc , V nloc1 , and V nloc2 for T below 70 K (where the simulation results were largely insensitive to the variation of σ bulk ), in the same manner as used to obtain the temperature dependence of the variables in Fig. 3(d).
As mentioned in Appendix C, we set the t s to be 5 nm, providing the justification.However, as the temperature increases, this simple model is no longer applicable.At high temperatures (in our sample, T > 70 K), the boundary between the bulk-conducting region and the surfaceconducting layer becomes obscure as thermally excited bulk carriers increase and the carrier-confining surface potential well becomes broadened.Consequently, our simulation model in Fig. 7 can be only for a sufficiently low temperature region, where the boundary between surface and bulk-conducting channels can be well defined.We found that our simulation failed to provide solutions for the observed values of V loc , V nloc1 , and V nloc2 , for any combination set of (σ bulk , G top □ , G bot □ ), for temperatures above 70 K.Thus, here, we provide only a qualitative explanation of the observed temperature dependence of V loc , V nloc1 , and V nloc2 in Fig. 1(c).
The total local current between the source and drain constitutes most of the bias current and thus does not vary significantly over the entire temperature range of measurements.Here, the local current through the bulk, which dominates the bias current at high temperatures, is simply redistributed with decreasing T into top and bottom surfaces between the source and drain.Thus, V loc ðTÞ reflects the T dependence of the effective sample conductance via the bulk and surface-conducting channels connected in parallel between the source and drain.In the high-T region, where the bulk conductance dominates the effective resistance, V loc increases quickly with decreasing T.However, in the low-temperature range, where the effective resistance between the source and drain is mostly governed by the surface conductance with metallic T dependence, V loc decreases with decreasing T. The broad crossover of high-to-low T behavior in V loc took place at 80-100 K in our flake.In contrast, the T dependence of V nloc1;2 , which is determined by I nloc1;2 ðTÞ × 1=G top □ ðTÞ, is more complicated.Here, the minute nonlocal current I nloc1;2 ðTÞ arises from part of the bias current leaking out of the local region, as the current flow through the bulk is blocked in the local region.The T dependence of I nloc1;2 ðTÞ is determined by the combined effect from σ bulk ðTÞ, G top □ ðTÞ, and G bot □ ðTÞ.At sufficiently low T, below 80-100 K, however, as most of the bias current is taken over by the surface current, the rate of increase of I nloc1;2 with decreasing T drops sharply.Combining this with a slow metallic increase of G top □ ðTÞ gives rise to almost temperature-independent V nloc1;2 ðTÞ below 80-100 K as in Fig. 1(c).But, with decreasing T above 100 K, a combination of a faster increase in I nloc1;2 ðTÞ and a metallic increase in G top □ leads to a maximum in V nloc1;2 in the temperature range of 100 < T < 300 K. Detailed temperature dependencies of V loc , V nloc1 , and V nloc2 in Fig. 1(c □ and G bot □ , leading to the matching condition for the experimental data of V loc , V nloc1 , and V nloc2 simultaneously for temperatures in the range T < 70 K.The inset in (b) shows the corresponding temperature variation of σ bulk .Grey shaded regions indicate the range of the variables giving rise to the observed data (with fluctuations) within an uncertainty of jΔV loc j ≤ 1 μV, jΔV nloc1 j ≤ 0.3 μV, and jΔV nloc2 j ≤ 0.05 μV.The upper limit of σ bulk is ∼0.6 S=m.Symbols indicate the weighted average values.

FIG. 1 .
FIG. 1.(a) Scanning-electron-microscopy image of our BSTS thin flake, with Ti/Au (10 nm=300 nm thick) electrodes.The dark-blue regions under the Au electrodes (see the magnified view in the inset) are cross-linked PMMA insulating layers, which prevent shorting between the Au electrodes and the BSTS flake.The scale bar is 2 μm in length.The sample dimensions are summarized in Table I.(b) Schematic diagram of the measurement configuration.(c) Temperature dependence of V loc (black line), V nloc1 (blue line), and V nloc2 (red line).

Figures 4 (
Figures 4(a)-4(c) show the surface current profile using a common color code for S top , S bulk , and S bot of the simulation model [see Fig. 3(a)] for ðσ bulk ; G top □ ; G bot □ Þ ¼ ð0.06 S=m; 6.0 e 2 =h; 2.5 e 2 =hÞ, which corresponds to the observed values of V loc ¼ 146 μV, V nloc1 ¼ 2.72 μV, and V nloc2 ¼ 0.58 μV at 4.2 K and V bg ¼ 0 V.The solid lines in Figs.4(a)-4(c) are contour lines that show the current flow, perpendicular to the equipotential lines.Here, we assume vanishing off-diagonal components of the conductance tensor [54].A significant current appears on the bottom surface [Fig.4(c)].However, the current density in the bulk almost completely vanishes [Fig.4(b)].Figures 4(d) and 4(e) show the V bg dependence of the simulated total current on the top surface (I top ; red squares), bottom surface (I bot ; blue circles), two side surfaces (I side ; black circles), and bulk (I bulk , magenta triangles) at S loc .As expected, I bot becomes minimal at the Dirac point of the bottom surface, while I top and I side show the opposite V bg dependence.Even with the large cross-sectional area of the bulk, I bulk ½¼ I bias − ðI top þ I bot þ I side Þ is at most 3 pA [Fig.4(e)].I nloc1 (green squares) in Fig.4(e) shows the nonlocal current on the top surface at S nloc1 .In contrast to I bulk , an appreciable amount of current appears on the top surface at

FIG. 4 .
FIG. 4. (a),(b),(c) Numerical simulation results on the current distribution of our BSTS flake in Fig.1(a).The current density profile at each plane, which is defined as in Fig.3(a), is plotted using the same color scale for ðG top □ ; G bot □ Þ ¼ ð6.0; 2.5Þðe 2 =hÞ and σ bulk ¼ 0.06 S=m, which correspond to the observed values of V loc ¼ 146 μV, V nloc1 ¼ 2.72 μV, and V nloc2 ¼ 0.58 μV at 4.2 K for V bg ¼ 0 V.The solid lines in (a), (b), and (c) are contour lines describing the current flow (i.e., these lines are perpendicular to the equipotential lines).(d) V bg dependence of I top (red squares) [I bot (blue circles)], the current on the top (bottom) surface at the cross section of S loc depicted in Fig.3(a).(e) V bg dependence of I side (black circles) [I bulk (magenta triangles)], the current on the side surface (inside bulk) at the cross section S loc depicted in Fig.3(a).I nloc1 is the nonlocal current on the top surface at S nloc1 depicted in Fig.3(a).

FIG. 5 .
FIG. 5. V bg dependence of (a) V loc vs B, (b) V nloc1 vs B, and (c) V nloc2 vs B. Each curve is shifted vertically for clarity.The universal conductance fluctuations are conspicuous in V nloc1;2 .(d),(e) Simulation results of ΔG top □¼ G top □ ðBÞ − G top □ ðB ¼ 0Þ and ΔG bot □ ¼ G bot □ ðBÞ − G bot □ ðB ¼ 0Þ for a given V bg .The solid lines in (d) and (e) are the best fits to Eq. (1).Each set of data, together with the fitting curve, is shifted vertically for clarity.(f),(g) V bg dependence of α and l ϕ obtained from the best fits to Eq. (1).Squares (circles) correspond to the top (bottom) surface.

FIG. 6 .
FIG.6.SEM image of side surfaces of a BSTS flake.Yellow boxes are magnified images for each side part of a BSTS flake.
) and 7(e) show the different simulation results for two models with different t s .If t s increases [from Fig. 7(d) to Fig. 7(e)], the electrical current takes a shorter path near the edge region [see Fig. 7(e)].Then, the current density profile should change according to the simulation results.However, as shown in Fig.

FIG. 7 .FIG. 8 .
FIG. 7. (a) Schematic top view of the BSTS flake.(b) The cross section of the BSTS flake with the actual aspect ratio.The surface channel regions cannot be shown in this schematic figure because of the small surface thickness (about 5 nm) compared to the thickness of the flake (about 110 nm).(c) An exaggerated view of the cross section of the BSTS flake corresponding to the region shown in (b).(d),(e) Influence of the surface thickness on the current path near the edge.
) and 9(b) show the T dependence of G top

FIG. 9 .
FIG. 9. Simulation results of G top□ and G bot □ , leading to the matching condition for the experimental data of V loc , V nloc1 , and V nloc2 simultaneously for temperatures in the range T < 70 K.The inset in (b) shows the corresponding temperature variation of σ bulk .Grey shaded regions indicate the range of the variables giving rise to the observed data (with fluctuations) within an uncertainty of jΔV loc j ≤ 1 μV, jΔV nloc1 j ≤ 0.3 μV, and jΔV nloc2 j ≤ 0.05 μV.The upper limit of σ bulk is ∼0.6 S=m.Symbols indicate the weighted average values.

TABLE I .
Sample dimensions.These parameters are used in numerical simulations.