††thanks: These two authors contributed equally.††thanks: These two authors
contributed equally.

# Signatures of long-range-correlated disorder in the magnetotransport of
ultrathin topological insulators

D. Nandi  Department of Physics, Harvard University, Cambridge, MA 02138, USA
Francis Bitter Magnet Lab and Plasma Science and Fusion Centre, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA  B. Skinner  Department of
Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA  G.H.
Lee  Department of Physics, Harvard University, Cambridge, MA 02138, USA
Department of Physics, Pohang University of Science and Technology, Pohang
790-784, Republic of Korea  K.-F. Huang  Department of Physics, Harvard
University, Cambridge, MA 02138, USA  K. Shain  Department of Physics, Harvard
University, Cambridge, MA 02138, USA  Cui-Zu Chang  Francis Bitter Magnet Lab
and Plasma Science and Fusion Centre, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA  Department of Physics, The Pennsylvania State
University, University Park, PA 16802-6300, USA  Y. Ou  Francis Bitter Magnet
Lab and Plasma Science and Fusion Centre, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  S.-P. Lee  Department of Physics,
University of Alberta, Edmonton, Alberta T6G 2E1, Canada  Department of
Physics and Astronomy, The John Hopkins University, Baltimore, Maryland 21218,
USA  J. Ward  Department of Physics, Harvard University, Cambridge, MA 02138,
USA  J.S. Moodera  Francis Bitter Magnet Lab and Plasma Science and Fusion
Centre, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
Department of Physics, Massachussetts Institute of Technology, Cambridge, MA
02139, USA  P. Kim  Department of Physics, Harvard University, Cambridge, MA
02138, USA  B.I. Halperin  Department of Physics, Harvard University,
Cambridge, MA 02138, USA  A. Yacoby  Department of Physics, Harvard
University, Cambridge, MA 02138, USA<EMAIL_ADDRESS>
###### Abstract

In an ultrathin topological insulator (TI) film, a hybridization gap opens in
the TI surface states, and the system is expected to become either a trivial
insulator or a quantum spin Hall insulator when the chemical potential is
within the hybridization gap. Here we show, however, that these insulating
states are destroyed by the presence of a large and long-range-correlated
disorder potential, which converts the expected insulator into a metal. We
perform transport measurements in ultrathin, dual-gated topological insulator
films as a function of temperature, gate voltage, and magnetic field, and we
observe a metallic-like, non-quantized conductivity, which exhibits a weak
antilocalization-like cusp at low magnetic field and gives way to a
nonsaturating linear magnetoresistance at large field. We explain these
results by considering the disordered network of electron- and hole-type
puddles induced by charged impurities. We argue theoretically that such
disorder can produce an insulator-to-metal transition as a function of
increasing disorder strength, and we derive a condition on the band gap and
the impurity concentration necessary to observe the insulating state. We also
explain the linear magnetoresistance in terms of strong spatial fluctuations
of the local conductivity, using both numerical simulations and a theoretical
scaling argument.

##  I Introduction

Three-dimensional topological insulators (TIs) are an exotic state of matter
in which gapless electronic excitations exist at the surface of a bulk system
with gapped conduction and valence bands.Hsieh  _et al._ (2009) These surface
states exhibit a number of interesting phenomena associated with their linear
dispersion and spin-momentum locking, Hasan and Kane (2010) including magnetic
monopole responses to applied electric field Qi  _et al._ (2009) and a strong
magnetoelectric effect. Essin  _et al._ (2009) Angle-resolved photoemission
spectroscopy (ARPES) measurements have identified gapless Dirac surface states
in several materials, including Bi1-xSbx, Bi2Se3, Sb2Te3, Bi2Te3. Hasan and
Kane (2010)

When a TI crystal is made very thin, however, the nature of the surface states
undergoes a significant change. In such ultrathin TI films, electrons have a
finite amplitude for quantum tunneling between the top and bottom surfaces,
resulting in a hybridization gap for the surface states whose magnitude
depends on the film thickness $d$. Under the right conditions, this gap can
stabilize the quantum spin Hall state, which is characterized by one-
dimensional, helical edge states around the border of the TI surface (as has
been observed in the 2D TI HgTe/CdTe König  _et al._ (2007); Roth  _et al._
(2009)).

One can therefore expect a basic dichotomy of possibilities for an undoped,
ultrathin TI film. Either the system becomes a trivial insulator, with a
vanishing conductivity in the limit of zero temperature, or it becomes a
quantum spin Hall insulator, with a quantized conductance. The fate of the TI
film vis-a-vis these two possibilities is predicted to depend in a nontrivial
way on the value of the thickness $d$, with the system oscillating between a
quantum spin Hall and a trivial insulating state as a function of
thickness.Liu  _et al._ (2010)

In this paper, however, we find evidence for a third possibility, outside of
this dichotomy, in which a hybridization gap exists but the insulating state
is destroyed by the presence of long-range-correlated disorder. We measure the
resistivity of ultrathin films of pristine (Bi,Sb)2Te3 as a function of
temperature, chemical potential, and magnetic field, and we find a number of
features that suggest that a dominant role is played by long-range-correlated
disorder, which arises inevitably due to charged impurities in the film and in
the substrate. Our films are 4 quintuple layers (QLs) thick, which is
predicted to produce a quantum spin Hall insulator, but we find instead a
finite, nonquantized value of the resistance in the limit of zero temperature.
We provide an explanation for this observation in terms of an insulator-to-
metal transition produced by increasing long-ranged disorder. We also observe
a prominent weak anti-localization correction and a large linear
magnetoresistance at high magnetic fields, which we explain in terms of strong
spatial fluctuations of the local conductivity.

The remainder of this paper is organized as follows. In Sec. II we briefly
describe our sample preparation and measurement setup. Section III describes
our zero-field measurements, and presents theoretical arguments to explain a
resistivity that is both finite and non-quantized in the limit of zero
temperature. Section IV presents results for the resistivity as a function of
magnetic field and gate voltage, along with a theoretical discussion of both
the weak antilocalization corrections and the linear magnetoresistance that we
observe.

##  II Sample preparation and measurement setup

Our TI films are made from 4 QLs of (Bi,Sb)2Te3 grown on a SrTiO3(111)
substrate using molecular-beam epitaxy in an ultra-high vacuum. Each QL layer
is $\sim 1$ nm thick. The Bi, Sb and Te effusion cells as well as the SrTiO3
(111) substrate are held at high temperature in order to ensure precise
control of surface stoichiometry. The crystallinity of the films is monitored
by reflection high-energy electron diffraction (RHEED) pattern. To prevent
oxidation in ambient conditions, $2$ nm of Tellurium and $2$ nm of Alumina
capping layer is grown on top of the films. The samples tend to degrade at
high temperatures, and hence all the processing was done at no higher than 100
∘C. Contacts were made using e-beam lithography. In the dual gated devices,
the top gate dielectric was made from 20 nm of HfO2 grown by atomic layer
deposition.Nandi  _et al._ (2018) Magneto-transport measurements are performed
in a dilution refrigerator with a 8 T magnet and using low frequency ac lockin
techniques. Results in this paper are taken from five Hall bar devices, which
we denote H1, H2, H3, H4 and H5 .

##  III Transport at zero magnetic field

As mentioned in the Introduction, hybridization between the two parallel TI
surfaces leads to a gap opening at the Dirac point of the surface dispersion
relation (as illustrated in Fig. 1). The size $\Delta$ of this gap is
generally expected to be between 5 meV and 50 meV; this range encompasses
estimates from density functional theory and tight-binding models for the gaps
in 4-QL-thick Bi2Se3 and Bi2Te3.Liu  _et al._ (2010) For situations where the
chemical potential resides in the middle of the gap and there is no band
bending, one would expect the TI surface to become an insulator with an
activation energy $\Delta/2$ for the conductivity. Such a state would have
either zero conductivity in the limit of zero temperature (if the system is a
trivial insulator) or a quantized conductance $2e^{2}/h$ (if the system forms
a quantum spin Hall state).

Figure 1: Schematic picture of the dispersion of the Dirac surface states (red
lines). Hybridization between top and bottom surfaces opens a small gap
$\Delta$ at the $\Gamma$ point. The bulk conduction and valence band states
are denoted by the upper and lower blue shaded areas, respectively.

Using a dual-gated FET setup, we shift the chemical potential for both the top
and bottom surfaces of our samples across a wide range in order to search for
these insulating states. At large negative voltages the chemical potential
resides far below the energy of the Dirac point, while at large positive
voltages the chemical potential is high above. When the two gate voltages are
chosen such that both surfaces are at the charge-neutral point (CNP), the
system assumes its maximally insulating state. This behavior is shown in Fig.
2 for five different samples as a function of the back gate voltage. (The
behavior as a function of both back and top gate voltages is discussed in the
Supplementary Material.)

Figure 2: Sheet resistivity of 5 different Hall bar devices — denoted H1 (blue
curve), H2 (brown), H3 (green), H4 (yellow) and H5 (purple) — as a function of
back gate voltage. The top gate voltage is held fixed at $V_{tg}$ = 0. The
measurement temperature was $\sim$ 30 mK.

Contrary to the expectation for a clean system, our measurements reveal a
conductivity that is neither insulating-like nor quantized. Indeed, Fig. 2
shows that the resistance takes a value of order $h/e^{2}$ at the CNP, but
this value varies from one sample to the next. As shown in Fig. 3, the
resistivity depends very weakly on temperature, even though the temperature is
far below $\Delta/k_{\text{B}}\gtrsim 50$ K. Using the traditional description
of an undoped semiconductor, one would predict an activated dependence
$\rho_{xx}\propto\exp[\Delta/(2k_{\text{B}}T)]$.

Figure 3: (a) The longitudinal resistivity $\rho_{xx}$ of device H1 as a
function of back gate voltage for different values of the temperature. (b) The
maximal resistivity (corresponding to the CNP) of device H1 is plotted as a
function of temperature, showing only a weak, logarithmic dependence on
temperature.

Taken together, these two observations suggest that the system is not well
described by either the clean band insulator or quantum spin Hall insulator
states. The most trivial explanation for our results would be that the surface
bands simply do not have a band gap. For example, in 4-QL-thick films of
Bi2Te3 the $\Gamma$ point lies at a local minimum of the valence band, Yazyev
_et al._ (2010) and consequently there is no finite window of energy with zero
density of states. In Sb2Te3, tunneling experiments have shown that films of
the same thickness have a thermodynamic gap that is no larger than a few meV.
Jiang  _et al._ (2012) Still, it is worth considering whether there is
another, more interesting explanation for the lack of insulating behavior,
especially since similar results have been recorded for transport in thin
films of Bi2Se3, Kim  _et al._ (2011) despite a Dirac point that lies well
outside the bulk valence band states and a gap $\Delta$ on the order of tens
of meV. Liu  _et al._ (2010); Zhang  _et al._ (2010)

The apparent breakdown of the clean insulator picture can be rationalized by
considering the effects of long-ranged disorder induced by Coulomb impurities,
which exist both in the TI film and in the substrate. Such impurities are
known to provide long-wavelength fluctuations of the local Fermi energy, which
provide a finite density of states at zero energy due to band-bending. Skinner
_et al._ (2013) This random band-bending effect is illustrated in Fig. 4. When
the surface gap $\Delta$ is sufficiently large, the fluctuations of the Fermi
energy lead to the formation of isolated electron and hole puddles [red and
blue regions of Fig. 4(a), respectively], separated by insulating tunnel
barriers (white regions). When the chemical potential $\mu$ is precisely in
the middle of the band gap (which we define as $\mu=0$), electron and hole
puddles appear in equal number. Similar picture of fluctuating Fermi energy in
the presence of insulating gap has been also adopted in graphene nanoribbons
Stampfer  _et al._ (2009) and dual-gated bilayer graphene Lee  _et al._
(2015a). One can estimate the condition for maintaining a good insulating
state at $\mu=0$ by demanding that the typical tunneling action $S$ for
electron tunneling across such a barrier satisfy $S\gg\hbar$.

Figure 4: (a) Schematic illustration of the disorder potential landscape.
Charged impurities produce a slowly-varying Coulomb potential that locally
creates puddles of electrons (red) and holes (blue). The typical correlation
length of the potential, $r_{s}$, is labeled, and thin black lines show
contours of constant potential. White regions denote tunnel barriers between
electron and hole puddles, which exist at large enough band gap $\Delta$. (b)
A schematic illustration of band bending, showing the energy along some
particular direction on the surface. The surface band gap $\Delta$ is labeled,
along with the typical magnitude $\Gamma$ of the disorder potential and the
width $x_{t}$ of the tunnel barrier between electron and hole puddles.

To estimate the action $S$, we first assume that the typical magnitude
$\Gamma$ of the disorder potential is sufficiently large that
$\Gamma\gg\Delta$. We also assume that the typical correlation length $r_{s}$
of the potential is much larger than the film thickness. Both of these
assumptions are validated below. Under these assumptions one can use numerical
estimates for $\Gamma$ and $r_{s}$ based on a gapless TI surface. Skinner  _et
al._ (2014), 111The theoretical estimates of Ref. Skinner  _et al._ , 2014 are
slightly modified for our case, since the two opposite TI surfaces screen the
disorder in tandem and therefore the density of states should be doubled
relative to the case of a single surface. These estimates give

$\Gamma\simeq\left(\frac{\pi^{3}\alpha_{s}^{2}}{2}\right)^{1/6}\hbar
vN_{\text{i}}^{1/3}$ (1)

and

$r_{s}\simeq\frac{N_{\text{i}}^{-1/3}}{(2\alpha_{s})^{4/3}},$ (2)

where $N_{\text{i}}$ is the (three-dimensional) concentration of impurities in
the TI film and substrate, $\hbar$ is the reduced Planck constant, $v$ is the
Dirac velocity, and $\alpha_{s}=e^{2}/(4\pi\varepsilon_{0}\varepsilon\hbar v)$
is the effective fine structure constant, with $\varepsilon$ the effective
dielectric constant and $\varepsilon_{0}$ the vacuum permittivity. If one
assumes an impurity concentration of order $10^{19}$ cm-3, a Dirac velocity
$v$ of order $5\times 10^{5}$ m/s, and an effective dielectric constant
$\varepsilon$ as large as several hundred (due to the close proximity of the
highly-polarizable SrTiO3 substrate), then $\alpha_{s}$ is of order $0.1$, the
disorder potential $\Gamma$ is of order $20$ meV, and $r_{s}\sim 40$ nm. It is
worth noting that at chemical potentials far from the CNP, both $r_{s}$ and
$\Gamma$ will generally be smaller than their $\mu=0$ values.

In order to estimate the width $x_{t}$ of the typical spatial separation
between electron and hole puddles, one can notice that the typical in-plane
electric field is $F\sim\Gamma/(er_{s})$, so that $x_{t}$ is given by
$eFx_{t}\sim\Delta$. Solving for $x_{t}$ and substituting Eqs. (1) and (2)
gives $x_{t}\sim
4\pi\varepsilon_{0}\varepsilon\Delta/(\alpha_{s}^{2/3}e^{2}N_{\text{i}}^{2/3})$.
The typical tunneling action between electron and hole puddles can be
estimated as the product of the height $\sim\Delta$ of the tunnel barrier and
the time $\sim x_{t}/v$ needed to traverse it. So $S\sim\Delta
x_{t}/v\sim\hbar\Delta^{2}/[\alpha_{s}^{5/3}(\hbar vN_{\text{i}}^{1/3})^{2}]$.
For the system to be insulating, one must have $S\gg\hbar$, which is
equivalent to

$\Delta\gg\alpha_{s}^{5/6}\hbar vN_{\text{i}}^{1/3}.$ (3)

Equation (3) can be viewed as a generic requirement for the existence of an
insulating state in a gapped 2D system on a substrate with charged impurities.
That is, either the gap $\Delta$ must be large enough or the impurity
concentration $N_{\text{i}}$ must be small enough that Eq. (3) is satisfied.

In our samples, the right-hand side of Eq. (3) is of order $10$ meV. For much
smaller values of the gap, one can say that electron and hole puddles are
well-connected by quantum tunneling, and there is no meaningful “insulating
barrier” between them. Our samples apparently correspond to such a situation,
where Eq. (3) is violated, so that one can think of the surface as effectively
ungapped even though $\Delta$ is finite. Producing a well-insulating TI thin
film apparently requires either a larger hybridization gap $\Delta$ or a much
smaller impurity concentration $N_{\text{i}}$. For the remainder of this paper
we set $\Delta=0$ when discussing transport.

The zero-field DC resistivity is given by $\rho=(h/e^{2})/(k_{\text{F}}\ell)$,
where $k_{\text{F}}$ is the typical Fermi wave vector and $\ell$ is the
electron mean free path. For the “puddled” scenario depicted in Fig. 4, the
typical value of $k_{\text{F}}$ at zero chemical potential is
$\sim\alpha_{s}^{1/3}N_{\text{i}}^{1/3}$, while the mean free path is of the
same order as $r_{s}$. Thus, the resistivity at zero chemical potential is
given by Skinner  _et al._ (2014)

$\rho_{\text{max}}\simeq\frac{h}{e^{2}}\alpha_{s}\ln(1/\alpha_{s}).$ (4)

For our samples, this expression gives a value of order $\approx 0.3h/e^{2}$.

Our picture of conduction through a spatially-disordered landscape is also
consistent with measurements of the superconducting proximity effect, which we
present in the Supplementary Material.

##  IV Magnetotransport

We also study the electron transport under the application of a perpendicular
magnetic field. The measured longitudinal resistance $R_{xx}$ is plotted in
Fig. 5 for device H4 as a function of the field strength $B$ and the back gate
voltage $V_{\text{bg}}$. For a given magnetic field, the resistance is
maximized near the CNP, which for this device corresponds to
$V_{\text{bg}}\approx 16$ V. The resistance also increases monotonically as a
function of $B$. In general, we observe an asymmetry between positive and
negative voltages relative to the CNP, with negative voltages generally
corresponding to smaller resistance. This asymmetry suggests that negative
values of the chemical potential correspond to a larger density of states than
positive values of the chemical potential, which may arise either because of
curvature of the Dirac band or because of proximity of the Dirac point to the
bulk valence band states (as depicted in Fig. 1). In our measurement
conditions, the Hall resistance $R_{xy}$ is everywhere much smaller than the
longitudinal resistance $R_{xx}$; this is shown explicitly in the
Supplementary Material. Thus we can approximate the conductivity
$\sigma\simeq(L/w)/R_{xx}$, where $L/w\approx 2$ is the aspect ratio of the
sample H5. The leading order correction to this expression is of order
$(R_{xy}/R_{xx})^{2}$, which is smaller than $5\%$ throughout the regime of
our measurements.

Figure 5: Longitudinal resistance $R_{xx}$ of device H5 as a function of back
gate voltage $V_{\text{bg}}$ and magnetic field $B$. The measurement
temperature $T=50$ mK. $R_{xy}$ is everywhere much smaller than $R_{xx}$.
Vertical dashed lines demarcate the three regimes of voltage depicted in Fig.
6(a)–(c).

Our data show two notable features as a function of magnetic field. For any
given gate voltage there is a sharp cusp in $R_{xx}(B)$ near $B=0$, which
previous experimental studies have attributed to weak antilocalization
(WAL).Steinberg  _et al._ (2011); Liu  _et al._ (2012); Tian  _et al._ (2014);
Zhang  _et al._ (2013); Kim  _et al._ (2011) The correction
$\Delta\sigma_{\text{WAL}}$ to the conductivity associated with WAL is
described by the theory of Hikami, Larkin, and Nagaoka: Hikami  _et al._
(1980)

$\Delta\sigma_{\text{WAL}}=\frac{\alpha}{\pi}\frac{e^{2}}{h}\left[\psi\left(\frac{1}{2}+\frac{\hbar}{4eBL_{\phi}^{2}}\right)-\ln\left(\frac{\hbar}{4eBL_{\phi}^{2}}\right)\right].$
(5)

Here $\psi(z)$ is the digamma function, $L_{\phi}$ is the phase coherence
length, and $\alpha$ is a numerical coefficient defined so that $\alpha<0$
indicates WAL and $\alpha>0$ corresponds to weak localization. The quantity
$2|\alpha|$ is usually associated with the number of parallel conduction
channels.

In Fig. 6 we plot $\Delta\sigma_{\text{WAL}}=\sigma(B)-\sigma(0)$ as a
function of magnetic field for different values of the back gate voltage. The
sharp, logarithmic cusp of $\Delta\sigma_{\text{WAL}}(B)$ is consistent with
Eq. (5), and we can perform good fits in the range
$-2\,\textrm{T}<B<2\,\textrm{T}$ in order to extract the parameters $\alpha$
and $L_{\phi}$. Interestingly, for low enough voltages that the chemical
potential is far below the CNP, we find that all measured curves
$\Delta\sigma_{\text{WAL}}(B)$ are identical, irrespective of the gate voltage
[Fig. 6(a)]. Similarly, all measured values of $\Delta\sigma_{\text{WAL}}(B)$
for voltages far above the CNP also collapse onto a single curve [Fig. 6(c)],
although this curve is distinct from the one corresponding to low voltages. In
the intermediate voltage range $9\textrm{ V}\leq V_{\text{bg}}\leq
23.4\textrm{ V}$ (between the two dashed vertical lines in Fig. 6), the
behavior of $\Delta\sigma_{\text{WAL}}(B)$ transitions smoothly from one
limiting curve to the other [Fig. 6(b)]. We interpret these two limiting
curves as corresponding to chemical potentials either far below or far above
the Dirac point. In the latter case, the current is carried only by the linear
Dirac surface states above the gap. In the opposite limit of small
$V_{\text{bg}}$, the current moves through an admixture of Dirac surface
states and conduction band states. For intermediate voltages, the disorder
potential mixes these two behaviors spatially by random band bending.

Figure 6: (a)–(c) The measured WAL correction to the conductivity of device H5
as a function of magnetic field for different ranges of the back gate voltage
$V_{\text{bg}}$ (the range of $V_{\text{bg}}$ is indicated in the title of
each plot, and is equivalent to the three ranges demarcated in Fig. 5). For
each plot the different curves correspond to different values of
$V_{\text{bg}}$. (d) The value of the constant $\alpha$ extracted from a fit
to Eq. (5) as a function of $V_{\text{bg}}$. (e) The extracted phase coherence
length as a function of $V_{\text{bg}}$. The measurement temperature is $T=50$
mK.

By fitting our data to Eq. (5), we are able to extract estimates for the
constant $\alpha$ and the phase coherence length $L_{\phi}$ as a function of
the back gate voltage. These results are shown in Fig. 6(d) and (e),
respectively. It is worth remarking that the inferred value of $\alpha$ is
everywhere close to $-1$, as one might expect for a conduction process with
two parallel channels (arising from the two parallel surfaces). The estimated
phase coherence length is on the order of several hundred nm, consistent with
previous studies at low temperature.Steinberg  _et al._ (2011); Tian  _et al._
(2014)

At larger magnetic field, the WAL correction gives way to a resistance that
increases linearly with magnetic field strength, with no evidence of
saturation. As shown in Fig. 7, this linear magnetoresistance (LMR) effect is
most prominent near the CNP. As the chemical potential is moved away from the
CNP in either direction, the slope of the LMR is reduced.

Multiple explanations have been proposed during the last few decades for
nonsaturating LMR in 2D electron systems. For example, Wang and Lei have
proposed a mechanism for LMR on a TI surface based on Zeeman splitting.Wang
and Lei (2012) However, the magnitude of the LMR associated with this
mechanism is much smaller than the value we observe. Indeed, to explain our
largest observed LMR slope with the Zeeman splitting mechanism would
apparently require an electron $g$-factor of several hundred, which seems
inconsistent with transport experiments in a tilted magnetic field. Chen  _et
al._ (2011) Other authors have explored more generic, semiclassical
explanations for LMR, and have shown how it can arise from either persistent
gradients of electron density Ilan  _et al._ (2006) or mesoscopic spatial
fluctuations of the mobility. Parish and Littlewood (2005); Kozlova  _et al._
(2012) Such fluctuations are commonly treated using either resistor network
models Parish and Littlewood (2005) or effective medium approaches Knap  _et
al._ (2014); Ping  _et al._ (2014), which give largely equivalent results.
Ramakrishnan  _et al._ (2017)

Figure 7:  (a) – (c) The resistance $R_{xx}$ of Hall bar device H5 as a
function of magnetic field is shown for different values of the back gate
voltage, whose value for each plot is indicated by the arrow pointing
downward. The thick black line in each plot shows the experimental data, and
the thin red line is a straight-line guide to the eye that indicates the
linear slope of $R_{xx}$ versus $|B|$. (d) The value of the slope $dR_{xx}/dB$
of device H5 at large $B$ is plotted as a function of back gate voltage
$V_{\text{bg}}$. $V_{\text{CNP}}\approx 16$ V indicates the voltage at the
charge neutral point.

Similar to these latter approaches, we suggest a way to understand our results
based on a simple model in which the system is described by a local, Drude-
like conductivity tensor that varies as a function of position due to spatial
fluctuations in the electron density. In particular, we suppose that one can
define local longitudinal and Hall conductivities $\sigma_{xx}(\mathbf{r})$
and $\sigma_{xy}(\mathbf{r})$, respectively, which vary as a function of
position $\mathbf{r}$. Such a description is generally valid so long as
variations in the electron density occur over a length scale $r_{s}$ that is
much longer than the mean free path $\ell$ or the Fermi wavelength $\sim
k_{\text{F}}^{-1}$.

In the Drude model, the ratio $\sigma_{xy}/\sigma_{xx}=\omega_{c}\tau$, where
$\omega_{c}$ is the cyclotron frequency and $\tau$ is the transport scattering
time. For a gapless Dirac system

$\omega_{c}\tau=\frac{eB\ell}{\hbar k_{\text{F}}},$ (6)

where $k_{\text{F}}$ is the local value of the Fermi momentum and $\ell=v\tau$
is the local mean free path. (As declared above, we are again ignoring the
effects of any small band gap $\Delta$.) Note that the ratio
$\sigma_{xy}/\sigma_{xx}$ becomes large when $B$ is sufficiently large. The
value of the local Fermi momentum can be described by the Thomas-Fermi
equation,

$E_{\text{F}}[k_{\text{F}}(\mathbf{r})]-e\phi(\mathbf{r})=\mu,$ (7)

where $E_{\text{F}}=\hbar vk_{\text{F}}\times\mathrm{sign}(\mu+e\phi)$ is the
Fermi energy relative to the Dirac point and $\phi(\mathbf{r})$ is the
electrostatic potential.

In order to understand the appearance of LMR within this model, consider first
the case when the system is close enough to the CNP that electron and hole
puddles exist in almost equal number, $|\mu|\ll\Gamma$. This is the regime
where the LMR is observed to be most prominent experimentally, and one can
understand its appearance using the following scaling arguments. Near the CNP,
electron and hole puddles are nearly equally abundant, and the local
conductivity at the boundary between puddles is small because of the locally
vanishing value of the electron density. Consequently, the current across the
system is forced to pass through narrow “pinch points” of the random
potential, where the electric potential is close to zero and adjacent
electron- or hole-type puddles are narrowly separated. These pinch points
provide the bottleneck for conduction, and they become more prominent with
increasing magnetic field. Ruzin (1993) (See the Supplementary Material for
simulated images of current flow.) One can think that an order unity number of
such pinch points exist per square area $\sim r_{s}^{2}$, and consequently, if
$G$ is the typical conductance of the pinch point, then the longitudinal
resistivity of the system is $\rho_{xx}\sim 1/G$.

To estimate the typical conductance $G$ of the pinch point, one can exploit
the result for the (two-terminal) conductance of a square with nonzero Hall
conductivityHans and Friedrich (1958); Jensen and Smith (1972); Abanin and
Levitov (2008): $G_{\square}=\sqrt{\sigma_{xx}^{2}+\sigma_{xy}^{2}}$. If the
magnetic field $B$ is large enough, then $\sigma_{xy}\gg\sigma_{xx}$ at the
pinch point and we arrive at the relation $\rho_{xx}\sim
1/|\sigma_{xy}^{(0)}|$, where $\sigma_{xy}^{(0)}$ represents the Hall
conductivity at the pinch point. At $\omega_{c}\tau\gg 1$ the Hall
conductivity $\sigma_{xy}\simeq en/B=\pm ek_{\text{F}}^{2}/(2\pi B)$, and so
our result for the longitudinal resistivity is equivalent to

$\rho_{xx}\sim\frac{h}{e^{2}}\frac{eB}{\hbar[k_{\text{F}}^{(0)}]^{2}},$ (8)

where $k_{\text{F}}^{(0)}$ represents the typical Fermi momentum at the pinch
point.

In general, pinch points are locations where the random potential is close to
zero. Thus, if the chemical potential is not too close to zero, then Eq. (7)
implies $k_{\text{F}}^{(0)}=|\mu|/(\hbar v)$. As the chemical potential is
shifted away from the CNP, the corresponding value of $k_{\text{F}}^{(0)}$
increases and the slope of the magnetoresistance declines. This is consistent
with the experimental result in Fig. 7. On the other hand, as one approaches
the CNP very closely, both the typical spatial size of the pinch point and its
typical Fermi momentum are reduced, and the resistance increases. While the
relation $k_{\text{F}}^{(0)}=|\mu|/(\hbar v)$ implies a divergence of the
resistance at $\mu\rightarrow 0$, such a divergence may be truncated by the
finite mean free path. In other words, since the local conductivity is not
well-defined at scales shorter than the mean free path $\ell$, one can think
that the minimal size of the pinch point is $\sim\ell$, and consequently that
the minimal value of $k_{\text{F}}^{(0)}$ is
$\sim\overline{k_{\text{F}}}\ell/r_{s}$, where $\overline{k_{\text{F}}}$ is
the typical value of $k_{\text{F}}$ near the center of an electron or hole
puddle. Using the estimates for $\overline{k_{\text{F}}}$ and $\ell$ presented
below, and inserting the expression for $k_{\text{F}}^{(0)}$ into Eq. (8),
gives a maximum magnetoresistance slope of order $\rho_{xx}\sim 0.03(h/e^{2})$
per Tesla of field. This is consistent in order of magnitude with our measured
result.

In order to test our scaling arguments quantitatively, we implemented numeric,
finite-element simulations of current flow through a Hall bar geometry with a
correlated disorder potential and a Drude-like conductivity tensor having a
local value of $k_{\text{F}}(\mathbf{r})$ given by Eq. (7). For simplicity,
our simulations assume a transport scattering time $\tau$ that is independent
of energy or position. While there is no reason a priori to expect this
assumption to be accurate quantitatively, the scaling argument leading to Eq.
(8) suggests that the longitudinal resistivity $\rho_{xx}$ becomes independent
of $\tau$ at sufficiently large field. As shown in Fig. 8, the simulation
consistently reproduces the LMR trend, as well as the decline in the LMR slope
with increasing chemical potential. Details of the simulation method are
provided in the Supplementary Material (along with results for $\rho_{xy}$,
which are consistent with experiment). Within the assumption of a constant
scattering time, one can fit the experimental data at the CNP quantitatively
by setting the root-mean-square Fermi momentum $\overline{k_{\text{F}}}=0.88$
nm-1, and the mean free path $\ell=7.4$ nm. More details about the fitting are
provided in the Supplementary Material.

Figure 8: Magnetoresistance for a simulated Hall bar with random disorder. (a)
Two example curves are plotted for $\rho_{xx}$ as a function of the
dimensionless magnetic field $B^{*}=B\ell^{2}/(\hbar e)$. Both curves
correspond to $\overline{k_{\text{F}}}\ell=6$, where $\overline{k_{\text{F}}}$
is the root-mean-square deviation of the Fermi momentum from its mean value
and $\ell$ is the mean free path, which is taken to be a constant. The curves
are labeled by their corresponding value of the chemical potential $\mu$,
normalized to the root-mean-square amplitude $\Gamma$ of the disorder
potential. (b) Shows a fit to the experimental data (thick blue curve,
corresponding to device H5) for $\rho_{xx}(B)$ at the charge-neutral point
($\mu=0$). The simulation data (black line with error bars) corresponds to
$\overline{k_{\text{F}}}=0.88$ nm-1 and $\overline{k_{\text{F}}}\ell=6.5$. (c)
The magnetoresistance slope $d\rho_{xx}/dB^{*}$ is plotted as a function of
chemical potential $\mu$, normalized to the disorder potential amplitude
$\Gamma$. The slope is calculated by a linear fit to simulation data in the
interval $0.2<B^{*}<1$, and in this example $\overline{k_{\text{F}}}\ell=6$ is
held constant.

##  V Conclusion

In this paper we have presented experimental results for the resistivity of
thin TI films as a function of temperature, chemical potential, and magnetic
field. In the absence of disorder, these systems are predicted to form a
quantum spin Hall state. We find, however, that the transport in our system is
dominated by long-range fluctuations of the disorder potential, presumably
induced by charged impurities in the film and in the substrate. In particular,
such long-ranged disorder creates a random landscape of $p$\- and $n$-type
regions, and this landscape destroys the insulating state. Our theoretical
arguments suggest that one may reach an insulating state only when the gap is
large enough and the impurity concentration is low enough that Eq. (3) is
satisfied.

The magnetotransport shows signs of both weak antilocalization and a
nonsaturating linear magnetoresistance. The WAL correction
$\Delta\sigma_{\text{WAL}}(B)$ is described well by the usual Hikami-Larkin-
Nagaoka theory, Eq. (5), with two parallel conduction channels
($\alpha\approx-1$). We also find that $\Delta\sigma_{\text{WAL}}(B)$
collapses onto one of two curves when the chemical potential is far from the
CNP. We have shown that the linear magnetoresistance can be interpreted as the
result of spatial fluctuations in the local conductivity arising from strong
disorder fluctuations. Our estimate of the linear MR slope and its dependence
on chemical potential are both consistent with observations.

Taken together, our results provide new understanding of electron transport in
ultrathin topological insulators, and may bring us closer to realizing ideal
quantum spin Hall insulators. More broadly, our improved understanding of the
disorder-induced insulator-to-metal transition and LMR may be important for a
wide class of disordered 2D electron systems.

###### Acknowledgements.

It is a pleasure to thank I. Sodemann, A. Nahum, S. Simon, J. T. Chalker, A.
Kamenev, and Y. P. Chen for valuable discussions. We are indebted to Di Wei,
Lucas Orona, Tony Zhou, Charlotte Bøttcher, Michael Kosowsky and Andrew Pierce
for invaluable help with fabrication and measurements. AY, DN, KS and JW
acknowledge the support from Gordon and Betty Moore Foundation Grant No. 4531,
NSF grant No. DMR-1708688, ARO Grant No. W911NF16-1-0491, ARO Grant No.
W911NF-17-1-0023, and ARO grant No. W911NF-18-1-0316. Fabrication for this
work was supported by DOE Award Number DE-SC0001819. BS was supported as part
of the MIT Center for Excitonics, an Energy Frontier Research Center funded by
the U.S. Department of Energy, Office of Science, Basic Energy Sciences under
Award No. DE-SC0001088. BIH acknowledges support from the STC Center for
Integrated Quantum Materials under NSF grant DMR-1231319. PK, GHL and KH
acknowledge support from NSF Grant No. DMR-1420634. JSM, CZC and YO
acknowledge the support from NSF Grant No. DMR-1700137, ONR Grant No.
N00014-16-1-2657 and the Center for Integrated Quantum Materials under NSF
Grant No. DMR-1231319. CZC also thanks the support from Alfred P. Sloan
Research Fellowship and ARO Young Investigator Program Award (W911NF1810198).

##  VI Supplemental Material for ‘Signatures of long-range-correlated disorder
in the magnetotransport of ultrathin topological insulators’

###  VI.1 Hall measurements

Figure 9: (a) Longitudinal resistance $R_{xx}$ and (b) Transverse Hall
resistance $R_{xy}$ for Hall bar device H3.

In the main text we focused on the longitudinal resistance $R_{xx}$ and
ignored the transverse Hall resistance $R_{xy}$. Here we present results for
$R_{xy}$, and we show that it is everywhere much smaller than $R_{xx}$. In
Fig. 9 we plot $R_{xx}$ and $R_{xy}$ for device H3 as a function of back gate
voltage $V_{\text{bg}}$ and magnetic field strength $B$.

From the $R_{xy}$ data, we calculate the inverse Hall coefficient $1/R_{H}$,
as shown in Fig. 10(a). The strong asymmetry in $R_{H}$ on different sides of
the charge neutral point (CNP) reflects an asymmetry in the band structure of
the TI films. The minimum value of $|1/R_{H}|\sim 3\times 10^{12}$ cm-2
suggests a typical carrier density of electron and hole puddles. Similar
results for $R_{H}$ were seen in Ref. Lee  _et al._ , 2015b.

Figure 10: (a) Inverse Hall coefficient as a function of back gate voltage.
(b) Ratio of $(\rho_{xy}/\rho_{xx})^{2}$ plotted as a function of magnetic
field and backgate voltage. Data corresponds to device H3, which has an aspect
ratio $L/w\approx 1$.

In the presence of a perpendicular magnetic field, the conductivity is given
in general by

$\displaystyle\sigma_{xx}$
$\displaystyle=\frac{\rho_{xx}}{(\rho_{xx}^{2}+\rho_{xy}^{2})}$
$\displaystyle=\frac{1}{\rho_{xx}(1+\rho_{xy}^{2}/\rho_{xx}^{2})}.$ (9)

Figure 10(b) shows that for all relevant magnetic field and backgate voltages,
the quantity $(\rho_{xy}/\rho_{xx})^{2}\ll 1$. Hence, the global conductivity
can be well approximated by $\sigma_{xx}\approx 1/\rho_{xx}$.

###  VI.2 Dual Gating

Figure 11: Dependence of the longitudinal resistance $R_{xx}$ on top and
backgate voltages for device H4.

In some of our samples, both a back gate and a top gate were fabricated,
allowing the chemical potential to be modulated by two independent gate
voltages. In Fig. 11 we present the low-temperature longitudinal resistance
$R_{xx}$ for sample H4 as a function of the top and back gate voltages,
$V_{\text{tg}}$ and $V_{\text{bg}}$, respectively. For this sample, the
longitudinal resistance exhibits a maximum of $\approx 0.5h/e^{2}$ at the CNP,
and there is no sign of an insulating state with large resistance.

###  VI.3 Superconducting proximity effect

Figure 12: (a) Optical image and (b) Fraunhoffer pattern of device JJ1.

Coupling superconductivity to topological insulator (TI) surface states is
predicted to create a px+ipy superconductor with Majorana bound states.Fu and
Kane (2008, 2008); Tanaka  _et al._ (2009); Potter and Fu (2013) Such Majorana
states may enable topologically protected qubits for quantum computation,Moore
(2010) and consequently there has been an intensive search for Majorana modes
in exfoliated 3D TIs Williams  _et al._ (2012); Rokhinson  _et al._ ; Kurter
_et al._ (2015). Achieving Josephson coupling is an important step towards
engineering topological superconductivity.

We investigated superconducting Josephson coupling mediated by the TI thin
film, as depicted schematically in Fig. 12(a). An important consideration for
achieving Josephson coupling is a transparent superconducting contact with
minimum contact resistance. To make transparent contacts, we etched Tellurium
capping layer covering topological insulator film (Bi,Sb)2Te3, and in-situ
evaporated Ti (5 nm) / Nb (2.5 nm)/ NbN (50 nm) without breaking vacuum in the
evaporation chamber. The device dimensions were 6 $\mu$m $\times$ 0.1 $\mu$m,
as shown in Fig. 12(b).

The interference pattern of the critical current $I_{c}$($B$) of the Josephson
junction under perpendicular magnetic field $B$ gives valuable information
about the Josephson current density $J_{c}$($x$) along the junction width. In
the case of a spatially uniform Josephson current density, the interference
pattern corresponds to a single slit Fraunhofer pattern where the lobes of
$I_{c}$ decays with $1/B$. On the other hand, for TIs in which conduction
occurs dominantly at the edges of the system, the interference pattern is
expected to resemble the double slit-like diffraction pattern of a dc SQUID
with non-decaying lobes. In this case the critical current $I_{c}(B)$ is given
by

$I_{c}(B)=I_{c}\left|\cos\left(\frac{\pi AB}{\phi_{0}}\right)\right|,$ (10)

where $\phi_{0}$ is the magnetic flux quantum and $A$ is the area of the
junction. Similar analyses of the interference pattern of Josephson junctions
has shown edge dominant conduction in HgTe/HgCdTe 2D topological insulator
Hart  _et al._ (2014) and Bi1.5Sb0.5Te1.7Se1.3 3D TIs. Lee  _et al._ (2014)

Figure 13: Example numeric solutions for the electrochemical potential $V$ and
the current $j_{x}$ in the $x$-direction for a given realization of the random
potential $P(x,y)$. (a) and (b) correspond to zero magnetic field, and (c) and
(d) correspond to the same random potential with a field $B^{*}=2$. (b) and
(d) have the same color scale. All images correspond to zero chemical
potential, $\mu=0$.

In our devices, the normal junction resistance $R_{N}\approx 60$ $\Omega$ and
the maximum Josephson current $I_{c}\approx 7\mu$A at the temperature of 50 mK
give $I_{c}R_{N}\approx 420\mu$eV, which corresponds to about
$\Delta_{\text{sc}}/6e$, where $\Delta_{\text{sc}}$ is the superconducting
gap. The Josephson junction displays an interference pattern, as shown in Fig.
12(c), that cannot be explained by a uniform Josephson current distribution.
The extracted period of the quasi-periodic oscillation of $I_{c}$ is about 35
G, which matches well with the expected value of $\phi_{0}/A=34$ G. Imperfect
constructive and destructive interferences visible in the pattern are
indicative of the effect of random disorder in the junction. Proximity-induced
Josephson coupling and the similar interference pattern were also observed by
us in Josephson junctions made of Ti (5 nm)/Al (50 nm).

Figure 14: Simulated values of the longitudinal (a) and Hall (b) resistivities
as a function of magnetic field and chemical potential. The vertical axis
corresponds to the dimensionless magnetic field $B^{*}$, while the horizontal
axis corresponds to the chemical potential $\mu$ in units of the root-mean-
square disorder potential $\Gamma$. In this example
$\overline{k_{\text{F}}}\ell=6$ everywhere. Compare to Fig. 9, and note that
the scale of $\rho_{xy}$ is much smaller than $\rho_{xx}$.

###  VI.4 Numeric simulations of current flow

In the main text we presented results for the magnetoresistance based on a
model where the local conductivity tensor $\hat{\sigma}(\mathbf{r})$ varies
with the position $\mathbf{r}$ due to spatial variations in the electron/hole
concentration. Here we present more details about our numeric simulations of
the resistivity.

Within the Drude model, the conductivity tensor has the form

$\hat{\sigma}(\mathbf{r})=\left(\begin{array}[]{cc}\sigma_{xx}(\mathbf{r})&-\sigma_{xy}(\mathbf{r})\\\
\sigma_{xy}(\mathbf{r})&\sigma_{xx}(\mathbf{r})\end{array}\right),$ (11)

where the values of $\sigma_{xx}(\mathbf{r})$ and $\sigma_{xy}(\mathbf{r})$
are related to the value of the local electrostatic potential by

$\displaystyle\sigma_{xx}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{h}k_{\text{F}}(\mathbf{r})\ell\frac{1}{1+[\omega_{c}(\mathbf{r})\tau]^{2}},$
$\displaystyle\sigma_{xy}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{h}k_{\text{F}}(\mathbf{r})\ell\frac{\omega_{c}(\mathbf{r})\tau}{1+[\omega_{c}(\mathbf{r})\tau]^{2}},$
(12)

where $k_{\text{F}}(\mathbf{r})$ is the local value of the Fermi momentum,
$\omega_{c}(\mathbf{r})$ is the local value of the cyclotron frequency, $\tau$
is the transport scattering time and $\ell=v\tau$ is the mean free path. In
principle, both $\ell$ and $\tau$ may have a dependency on the local Fermi
energy, and therefore on position. However, for simplicity in our numeric
simulations we take $\ell$ and $\tau$ to be constants with no spatial
variation. The local Fermi momentum $k_{\text{F}}(\mathbf{r})$ is related to
the disorder potential $\phi(\mathbf{r})$ by the Thomas-Fermi equation, Eq.
(7) of the main text.

In order to address the problem computationally, we define the following
dimensionless units. First, we define a normalized electrochemical potential

$P(\mathbf{r})=\frac{e\phi(\mathbf{r})+\mu}{\sqrt{\langle(e\phi)^{2}\rangle_{\mathbf{r}}}},$
(13)

so that $P(\mathbf{r})$ has a standard deviation of unity. Here, $\mu$ is the
chemical potential and
$\sqrt{\langle(e\phi)^{2}\rangle_{\mathbf{r}}}\equiv\Gamma$ is the root-mean-
square value of the disorder potential. The local value of the Fermi momentum
is then given by

$k_{\text{F}}(\mathbf{r})=\overline{k_{\text{F}}}\left|P(\mathbf{r})\right|,$
(14)

where $\overline{k_{\text{F}}}=\Gamma/(\hbar v)$ represents the root-mean-
square value of the Fermi momentum at zero chemical potential.

We also define a dimensionless magnetic field strength

$B^{*}\equiv\frac{eB\ell^{2}}{\hbar}=\left(\frac{\ell}{\ell_{B}}\right)^{2},$
(15)

where $\ell_{B}=\sqrt{\hbar/(eB)}$ is the magnetic length. With these
definitions,
$\omega_{c}\tau=(\overline{k_{\text{F}}}\ell)^{-1}B^{*}/P(\mathbf{r})$, and we
can rewrite the elements of the conductivity tensor in dimensionless form as

$\displaystyle\frac{\sigma_{xx}(\mathbf{r})}{e^{2}/h}$ $\displaystyle=$
$\displaystyle\frac{\overline{k_{\text{F}}}\ell\left|P(\mathbf{r})\right|}{1+\left(\frac{1}{\overline{k_{\text{F}}}\ell}\right)^{2}\left[\frac{B^{*}}{P(\mathbf{r})}\right]^{2}},$
$\displaystyle\frac{\sigma_{xy}(\mathbf{r})}{e^{2}/h}$ $\displaystyle=$
$\displaystyle\frac{B^{*}\textrm{sign}\\{P(\mathbf{r})\\}}{1+\left(\frac{1}{\overline{k_{\text{F}}}\ell}\right)^{2}\left[\frac{B^{*}}{P(\mathbf{r})}\right]^{2}}.$
(16)

Written in this form, the system is characterized by three dimensionless
parameters: $B^{*}$, $\mu/\Gamma$, and $\overline{k_{\text{F}}}\ell$. For a
given choice of these parameters, and for a given realization of the random
potential $P(\mathbf{r})$, one can solve for the current density
$\mathbf{j}(\mathbf{r})$ through the system by solving the continuity equation

$\mathbf{\nabla}\cdot\mathbf{j}=0,$ (17)

where $\mathbf{j}=-\hat{\sigma}\mathbf{\nabla}V(\mathbf{r})$, and
$V(\mathbf{r})$ is the local deviation of the electrochemical potential away
from equilibrium. Below we present results based on a finite element solution
to Eq. (17). For definiteness, we choose our simulated system to be a Hall bar
with length $L=8\xi$ and width $w=4\xi$, where $\xi$ is the correlation length
of the potential (equivalent to $r_{s}$ in the main text). Applying a unit
voltage across the long end of the bar gives the boundary conditions
$V(x=0)=1$, $V(x=L)=0$, and $j_{y}(y=0)=j_{y}(y=w)=0$. The finite-element mesh
size was $\xi/20$.

The random potential was taken to be a Fourier series with random coefficients
whose magnitude decreases at high wave vector. Specifically, we write

$P(x,y)=\textrm{Re}\left\\{\sum_{m,n=-\infty}^{\infty}c_{m,n}e^{i(k_{m}x+k_{n}y)}\right\\}+\frac{\mu}{\Gamma},$
(18)

where $c_{m,n}$ are random coefficients and

$k_{m}=\frac{\pi m}{L},\hskip 14.22636ptk_{n}=\frac{\pi n}{w}.$ (19)

The coefficients $c_{m,n}$ are chosen to have a random phase in the complex
plane and a random magnitude bounded by

$|c_{m,n}|^{2}\leq e^{-(k_{m}^{2}+k_{n}^{2})\xi^{2}/2}.$ (20)

The normalization of the potential is such that $c_{0,0}=0$ and
$\frac{1}{2}\sum_{m,n}|c_{m,n}|^{2}=1$.

Figure 13 shows a typical numerical solution for the potential $V(x,y)$ and
the current $j_{x}$ in the $x$ direction at $\mu=0$. Figures 13(a)-(b) show
the potential and current at zero field, $B^{*}=0$, while Figs. 13(c)-(d) show
the same system at large field $B^{*}=2$. One can see from these images the
strong effect of “pinch points” in the random potential, where the current is
concentrated at narrow constrictions between puddles. This focusing of the
current becomes more exaggerated at large magnetic field, and correspondingly
the electrochemical potential $V(\mathbf{r})$ drops abruptly at pinch points
and becomes relatively constant far from the boundaries between n- and p-type
regions.

For a given realization of the random potential, the longitudinal resistivity
can be defined from the simulation as $\rho_{xx}=\Delta V/(L\langle
j_{x}\rangle)$, where $\langle j_{x}\rangle$ is the area-averaged current
density in the $z$ direction and $\Delta V\equiv 1$ is the voltage drop across
the system. One can estimate the dependence of $\rho_{xx}$ on magnetic field
$B^{*}$ by averaging the simulated resistivity $\rho_{xx}$ over many
realizations of the random potential for a given value of $B^{*}$. Our
numerical results, including those shown in Fig. 8 of the main text, are
averaged over 100 such realizations. One can also define the Hall resistivity
$\rho_{xy}$ by numerically extracting the transverse voltage $V_{H}$ across
the Hall bar at the midpoint, $x=4\xi$, for a given realization. The Hall
resistivity is defined by $\rho_{xy}=V_{H}/(w\langle j_{x}\rangle)$.

In Fig. 14 we show the values of $\rho_{xx}$ and $\rho_{xy}$ as produced by
our simulation for a range of magnetic fields $B^{*}$ and $\mu/\Gamma$. As in
the experimental data (see Fig. 9), $\rho_{xy}$ is much smaller than
$\rho_{xy}$ throughout the range of interest.

In order to fit the experimental data at $\mu=0$ to the numerical simulation,
we calculate numerically the resistivity $\rho_{xx}/(h/e^{2})$ as a function
of $B^{*}$, using discrete points $B^{*}=0,0.1,0.2,...,2$, for
$\overline{k_{\text{F}}}\ell=5,6$, and $7$. Linear interpolation allows us to
estimate the value of $\rho_{xx}$ for generic values of
$\overline{k_{\text{F}}}\ell$ between $5$ and $7$ and $B^{*}$ between $0$ and
$2$. The resulting curves can be translated into real units by inserting the
corresponding values of $\overline{k_{\text{F}}}$ and $\ell$, which allows us
to fit the data shown in Fig. 7(b) of the main text. The result of this
fitting is shown in Fig. 8(b) of the main text, and corresponds to
$\overline{k_{\text{F}}}\ell=6.5$, and $\overline{k_{\text{F}}}=0.88$ nm-1.
This value of $\overline{k_{\text{F}}}$ implies a typical electron/hole
density $2\pi\overline{k_{\text{F}}}^{2}=1.2\times 10^{13}$ cm-2, which is
consistent in order of magnitude with both our theoretical estimates and our
measurements of the Hall constant.

We extracted the dependence of the linear magnetoresistance slope
$d[\rho_{xx}/(h/e^{2})]/dB^{*}$ on the chemical potential by numerically
evaluating the curve $\rho_{xx}(B^{*})$ for different values of $\mu/\Gamma$,
and then fitting each curve to a line over the range $0.2<B^{*}<1$.

## References

  * Hsieh  _et al._ (2009) David Hsieh, Y Xia, Dong Qian, L Wray, J H Dil, F Meier, J Osterwalder, L Patthey, J G Checkelsky, N P Ong, _et al._ , “A tunable topological insulator in the spin helical dirac transport regime,” Nature 460, 1101–1105 (2009).
  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, “Colloquium,” Rev. Mod. Phys. 82, 3045–3067 (2010).
  * Qi  _et al._ (2009) Xiao Liang Qi, Rundong Li, Jiadong Zang, and Shou Cheng Zhang, “Inducing a magnetic monopole with topological surface states,” Science 323, 1184–1187 (2009).
  * Essin  _et al._ (2009) Andrew M. Essin, Joel E. Moore, and David Vanderbilt, “Magnetoelectric polarizability and axion electrodynamics in crystalline insulators,” Phys. Rev. Lett. 102, 146805 (2009).
  * König  _et al._ (2007) Markus König, Steffen Wiedmann, Christoph Brüne, Andreas Roth, Hartmut Buhmann, Laurens W Molenkamp, Xiao Liang Qi, and Shou Cheng Zhang, “Quantum spin hall insulator state in hgte quantum wells,” Science 318, 766–770 (2007).
  * Roth  _et al._ (2009) Andreas Roth, Christoph Brüne, Hartmut Buhmann, Laurens W Molenkamp, Joseph Maciejko, Xiao Liang Qi, and Shou Cheng Zhang, “Nonlocal transport in the quantum spin hall state,” Science 325, 294–297 (2009).
  * Liu  _et al._ (2010) Chao Xing Liu, HaiJun Zhang, Binghai Yan, Xiao Liang Qi, Thomas Frauenheim, Xi Dai, Zhong Fang, and Shou Cheng Zhang, “Oscillatory crossover from two-dimensional to three-dimensional topological insulators,” Phys. Rev. B 81, 041307 (2010).
  * Nandi  _et al._ (2018) D. Nandi, Inti Sodemann, K. Shain, G. H. Lee, K.-F. Huang, Cui-Zu Chang, Yunbo Ou, S. P. Lee, J. Ward, J. S. Moodera, P. Kim, and A. Yacoby, “Logarithmic singularities and quantum oscillations in magnetically doped topological insulators,” Phys. Rev. B 97, 085151 (2018).
  * Yazyev  _et al._ (2010) Oleg V. Yazyev, Joel E. Moore, and Steven G. Louie, “Spin polarization and transport of surface states in the topological insulators ${\mathrm{bi}}_{2}{\mathrm{se}}_{3}$ and ${\mathrm{bi}}_{2}{\mathrm{te}}_{3}$ from first principles,” Phys. Rev. Lett. 105, 266806 (2010).
  * Jiang  _et al._ (2012) Yeping Jiang, Yilin Wang, Mu Chen, Zhi Li, Canli Song, Ke He, Lili Wang, Xi Chen, Xucun Ma, and Qi-Kun Xue, “Landau quantization and the thickness limit of topological insulator thin films of ${\mathrm{sb}}_{2}{\mathrm{te}}_{3}$,” Phys. Rev. Lett. 108, 016401 (2012).
  * Kim  _et al._ (2011) Yong Seung Kim, Matthew Brahlek, Namrata Bansal, Eliav Edrey, Gary A Kapilevich, Keiko Iida, Makoto Tanimura, Yoichi Horibe, Sang Wook Cheong, and Seongshik Oh, “Thickness-dependent bulk properties and weak antilocalization effect in topological insulator bi 2 se 3,” Physical Review B 84, 073109 (2011).
  * Zhang  _et al._ (2010) Yi Zhang, Ke He, Cui Zu Chang, Can Li Song, Li Li Wang, Xi Chen, Jin Feng Jia, Zhong Fang, Xi Dai, Wen Yu Shan, _et al._ , “Crossover of the three-dimensional topological insulator bi2se3 to the two-dimensional limit,” Nature Physics 6, 584–588 (2010).
  * Skinner  _et al._ (2013) B. Skinner, T. Chen, and B. I. Shklovskii, “Effects of bulk charged impurities on the bulk and surface transport in three-dimensional topological insulators,” Journal of Experimental and Theoretical Physics 117, 579–592 (2013).
  * Stampfer  _et al._ (2009) C. Stampfer, J. Güttinger, S. Hellmüller, F. Molitor, K. Ensslin, and T. Ihn, “Energy gaps in etched graphene nanoribbons,” Phys. Rev. Lett. 102, 056403 (2009).
  * Lee  _et al._ (2015a) Gil-Ho Lee, Dongchan Jeong, Kee-Su Park, Yigal Meir, Min-Chul Cha, and Hu-Jong Lee, “Continuous and reversible tuning of the disorder-driven superconductor–insulator transition in bilayer graphene,” Scientific reports 5, 13466 (2015a).
  * Skinner  _et al._ (2014) Brian Skinner, B. I. Shklovskii, and M. B. Voloshin, “Bound state energy of a coulomb impurity in gapped bilayer graphene,” Phys. Rev. B 89, 041405 (2014).
  * Note (1) The theoretical estimates of Ref. skinner2013theory are slightly modified for our case, since the two opposite TI surfaces screen the disorder in tandem and therefore the density of states should be doubled relative to the case of a single surface.
  * Steinberg  _et al._ (2011) H. Steinberg, J.-B. Laloë, V. Fatemi, J. S. Moodera, and P. Jarillo-Herrero, “Electrically tunable surface-to-bulk coherent coupling in topological insulator thin films,” Phys. Rev. B 84, 233101 (2011).
  * Liu  _et al._ (2012) Minhao Liu, Jinsong Zhang, Cui-Zu Chang, Zuocheng Zhang, Xiao Feng, Kang Li, Ke He, Li-li Wang, Xi Chen, Xi Dai, Zhong Fang, Qi-Kun Xue, Xucun Ma, and Yayu Wang, “Crossover between weak antilocalization and weak localization in a magnetically doped topological insulator,” Phys. Rev. Lett. 108, 036805 (2012).
  * Tian  _et al._ (2014) Jifa Tian, Cuizu Chang, Helin Cao, Ke He, Xucun Ma, Qikun Xue, and Yong P. Chen, “Quantum and classical magnetoresistance in ambipolar topological insulator transistors with gate-tunable bulk and surface conduction,” Scientific Reports 4, 4859 EP – (2014).
  * Zhang  _et al._ (2013) Li Zhang, Merav Dolev, Qi I. Yang, Robert H. Hammond, Bo Zhou, Alexander Palevski, Yulin Chen, and Aharon Kapitulnik, “Weak localization effects as evidence for bulk quantization in bi2se3 thin films,” Phys. Rev. B 88, 121103 (2013).
  * Hikami  _et al._ (1980) Shinobu Hikami, Anatoly I Larkin, and Yosuke Nagaoka, “Spin-orbit interaction and magnetoresistance in the two dimensional random system,” Progress of Theoretical Physics 63, 707–710 (1980).
  * Wang and Lei (2012) C. M. Wang and X. L. Lei, “Linear magnetoresistance on the topological surface,” Phys. Rev. B 86, 035442 (2012).
  * Chen  _et al._ (2011) J. Chen, X. Y. He, K. H. Wu, Z. Q. Ji, L. Lu, J. R. Shi, J. H. Smet, and Y. Q. Li, “Tunable surface conductivity in bi2se3 revealed in diffusive electron transport,” Phys. Rev. B 83, 241304 (2011).
  * Ilan  _et al._ (2006) R. Ilan, N. R. Cooper, and Ady Stern, “Longitudinal resistance of a quantum hall system with a density gradient,” Phys. Rev. B 73, 235333 (2006).
  * Parish and Littlewood (2005) Meera M. Parish and Peter B. Littlewood, “Classical magnetotransport of inhomogeneous conductors,” Phys. Rev. B 72, 094417 (2005).
  * Kozlova  _et al._ (2012) N. V. Kozlova, N. Mori, O. Makarovsky, L. Eaves, Q. D. Zhuang, A. Krier, and A. Patanè, “Linear magnetoresistance due to multiple-electron scattering by low-mobility islands in an inhomogeneous conductor,” Nature Communications 3, 1097 EP – (2012).
  * Knap  _et al._ (2014) Michael Knap, Jay D. Sau, Bertrand I. Halperin, and Eugene Demler, “Transport in two-dimensional disordered semimetals,” Phys. Rev. Lett. 113, 186801 (2014).
  * Ping  _et al._ (2014) Jinglei Ping, Indra Yudhistira, Navneeth Ramakrishnan, Sungjae Cho, Shaffique Adam, and Michael S. Fuhrer, “Disorder-induced magnetoresistance in a two-dimensional electron system,” Phys. Rev. Lett. 113, 047206 (2014).
  * Ramakrishnan  _et al._ (2017) Navneeth Ramakrishnan, Ying Tong Lai, Silvia Lara, Meera M. Parish, and Shaffique Adam, “Equivalence of effective medium and random resistor network models for disorder-induced unsaturating linear magnetoresistance,” Phys. Rev. B 96, 224203 (2017).
  * Ruzin (1993) I. M. Ruzin, “Hall transport in nonuniform two-dimensional conductors,” Phys. Rev. B 47, 15727–15734 (1993).
  * Hans and Friedrich (1958) Joachim Hans and Kuhrt Friedrich, “zna,” (1958) Chap. Der Geometrieeinfluß auf den transversalen magnetischen Widerstandseffekt bei rechteckförmigen Halbleiterplatten, p. 462, 6.
  * Jensen and Smith (1972) H Hojgaard Jensen and H Smith, “Geometrical effects in measurements of magnetoresistance,” Journal of Physics C: Solid State Physics 5, 2867 (1972).
  * Abanin and Levitov (2008) Dmitry A. Abanin and Leonid S. Levitov, “Conformal invariance and shape-dependent conductance of graphene samples,” Phys. Rev. B 78, 035416 (2008).
  * Lee  _et al._ (2015b) Joon Sue Lee, Anthony Richardella, Danielle Reifsnyder Hickey, K. Andre Mkhoyan, and Nitin Samarth, “Mapping the chemical potential dependence of current-induced spin polarization in a topological insulator,” Phys. Rev. B 92, 155312 (2015b).
  * Fu and Kane (2008) Liang Fu and C. L. Kane, “Superconducting proximity effect and majorana fermions at the surface of a topological insulator,” Phys. Rev. Lett. 100, 096407 (2008).
  * Tanaka  _et al._ (2009) Yukio Tanaka, Takehito Yokoyama, and Naoto Nagaosa, “Manipulation of the majorana fermion, andreev reflection, and josephson current on topological insulators,” Phys. Rev. Lett. 103, 107002 (2009).
  * Potter and Fu (2013) Andrew C. Potter and Liang Fu, “Anomalous supercurrent from majorana states in topological insulator josephson junctions,” Phys. Rev. B 88, 121109 (2013).
  * Moore (2010) Joel E Moore, “The birth of topological insulators,” Nature 464, 194 (2010).
  * Williams  _et al._ (2012) J. R. Williams, A. J. Bestwick, P. Gallagher, Seung Sae Hong, Y. Cui, Andrew S. Bleich, J. G. Analytis, I. R. Fisher, and D. Goldhaber Gordon, “Unconventional josephson effect in hybrid superconductor-topological insulator devices,” Phys. Rev. Lett. 109, 056803 (2012).
  * (41) Leonid P Rokhinson, Xinyu Liu, and Jacek K Furdyna, “The fractional ac josephson effect in a semiconductor-superconductor nanowire as a signature of majorana particles,” Nature Physics 8, 795.
  * Kurter  _et al._ (2015) C Kurter, A D Finck, Y S Hor, and D J Van Harlingen, “Evidence for an anomalous current-phase relation in topological insulator josephson junctions,” Nature communications 6, 7130–7130 (2015).
  * Hart  _et al._ (2014) Sean Hart, Hechen Ren, Timo Wagner, Philipp Leubner, Mathias Mühlbauer, Christoph Brüne, Hartmut Buhmann, Laurens W Molenkamp, and Amir Yacoby, “Induced superconductivity in the quantum spin hall edge,” Nature Physics 10, 638–643 (2014).
  * Lee  _et al._ (2014) Jae Hyeong Lee, Gil Ho Lee, Joonbum Park, Janghee Lee, Seung Geol Nam, Yun Sok Shin, Jun Sung Kim, and Hu Jong Lee, “Local and nonlocal fraunhofer-like pattern from an edge-stepped topological surface josephson current distribution,” Nano letters 14, 5029–5034 (2014).