{'fraction_non_alphanumeric': 0.0370469541933633, 'fraction_numerical': 0.01579803002773262, 'mean_word_length': 4.102068696330991, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Effect of bulk charged impurities on the bulk and surface transport in three-dimensional topological insulators

In the three-dimensional topological insulator (TI), the physics of doped semiconductors exists literally side-by-side with the physics of ultra-relativistic Dirac fermions. This unusual pairing creates a novel playground for studying the interplay between disorder and electronic transport. In this mini-review we focus on the disorder caused by the three-dimensionally distributed charged impurities that are ubiquitous in TIs, and we outline the effects it has on both the bulk and surface transport in TIs. We present self-consistent theories for Coulomb screening both in the bulk and at the surface, discuss the magnitude of the disorder potential in each case, and present results for the conductivity. In the bulk, where the band gap leads to thermally activated transport, we show how disorder leads to a smaller-than-expected activation energy that gives way to VRH at low temperatures. We confirm this enhanced conductivity with numerical simulations that also allow us to explore different degrees of impurity compensation. For the surface, where the TI has gapless Dirac modes, we present a theory of disorder and screening of deep impurities, and we calculate the corresponding zero-temperature conductivity. We also comment on the growth of the disorder potential as one moves from the surface of the TI into the bulk. Finally, we discuss how the presence of a gap at the Dirac point, introduced by some source of time-reversal symmetry breaking, affects the disorder potential at the surface and the mid-gap density of states.


I. INTRODUCTION
The three-dimensional (3D) topological insulator (TI) 1-5 has generated a great deal of excitement in the physics community because of its gapless surface states, which host a spectrum of quantum transport phenomena 6,7 . Unfortunately, while a number of crystals have been identified to be 3D TIs, most of them are not actually insulators, but instead have a relatively large bulk conductivity that shunts the surface conductivity for TI crystals of substantial thickness ( 10 µm). How to achieve a bulk-insulating state is a problem that is widely discussed in the current literature [8][9][10][11][12][13][14][15][16] .
Typically, as-grown TI crystals are heavily doped ntype semiconductors, so that the Fermi level resides in the bulk conduction band. In order to arrive at a bulk insulating state, such TIs are compensated by acceptors. With increasing compensation K = N A /N D , where N D and N A are the concentrations of monovalent donors and acceptors, respectively, the Fermi level shifts from the conduction band to inside the gap and then into the valence band. When compensation of donors is complete, K = 1, the Fermi level is in the middle of the gap and the most insulating state of the TI is reached. The hope is that for a TI with bulk band gap E g ∼ 0.3 eV (as, for example, in Bi 2 Se 3 ) the bulk resistivity should obey the activation law with activation energy ∆ = E g /2 ∼ 0.15 eV, so that at room temperatures and below the TI is well insulating. The typical experimental situation near K = 1, however, is frustrating 15 . In the range of temperatures between 100 K and 300 K the resistivity is activated, but with an activation energy that is roughly three times smaller than expected, ∆ ∼ 50 meV. At T ∼ 100 K the activated transport is replaced by variable range hopping (VRH) and the resistivity grows even more slowly with decreasing T . Finally, at even smaller temperature, T < 50 K, the resistivity saturates 17 at a value < 10 Ωcm.
In a recent paper 18 we showed that the unexpectedly large bulk conductivity of TIs at K = 1 can be explained as a consequence of the enormously-fluctuating random Coulomb potential created by randomly-positioned donor and acceptor impurities. In later papers we extended this analysis to the case of near complete compensation 19 , K < 1 and 1 − K 1, and we examined the effect of random Coulomb impurities on the surface disorder and transport properties 20 . In this mini-review our goal is to outline in a general way the effects of random, 3Ddistributed Coulomb impurities in TIs on both the bulk and surface properties. We describe the screening mechanisms for the random Coulomb potential both within the bulk of the TI and at the surface, and we present predictions for the magnitude of the disorder potential and the conductivity.
Our theoretical treatment is also motivated by the recent experiments of Ref. 21, where the random potential at the surface of typical TIs (Bi 2 Se 3 and Bi 2 Te 3 ) was studied directly by spectroscopic mapping with a scanning tunneling microscope. It was shown that near the Dirac energy random fluctuations of the potential have a Gaussian-like distribution with a width ∼ 20 -40 meV that can be attributed to deep impurity charges. We show below that such fluctuations are consistent with disorder produced by three-dimensionally distributed bulk Coulomb impurities that are screened by the gapless TI surface.
Crucial to our theoretical description throughout this paper is the assumption of a random spatial distribution of impurities. This assumption is readily justified for TI samples made by cooling from a melt, where the distribution of impurities in space is a snapshot of the distribution that impurities have at higher temperature, when their diffusion practically freezes 22 . In 3D TIs, as in conventional narrow band gap semiconductors, at this temperature there is a concentration of intrinsic carriers larger than the concentration of impurities. Intrinsic carriers thus screen the Coulomb interaction between impurities, so that impurities remain randomly distributed in space. When the temperature is lowered to the point where intrinsic carriers recombine, the impurities are left in random positions 23,24 . If the diffusion of impurities freezes at T ∼ 1000 K it is reasonable to assume that impurities are randomly positioned for semiconductors with bulk band gap E g ≤ 0.3 eV. Throughout this paper we deal with such narrow band gap TIs, such as Bi 2 Se 3 , for which our description of randomly-positioned impurities is accurate. We also assume everywhere that donor and acceptor energy levels are shallow, meaning that their binding energy is much smaller than E g . The remainder of this paper can be divided into two parts. In the first part, comprising Secs. II -IV, we focus on bulk properties, essentially treating the TI as a strongly-or completely-compensated semiconductor and ignoring the surface states. In Sec. II we give a conceptual explanation of the bulk disorder potential and the origin of the anomalously small bulk resistivity. Sec. III formulates a numerical model of the TI bulk and uses it to calculate the corresponding electron density of states (DOS). In Sec. IV we present our algorithm for computing the thermally activated conductivity, analyze our results, and arrive at an expression for the unusually small bulk activation energy. We also evaluate the localization length of states close to the Fermi energy and estimate the characteristic temperature associated with variablerange hopping.
The second part of this paper, comprising Secs. V -VIII, deals with the effects of Coulomb impurities on the properties of the TI surface. In Sec. V we describe a selfconsistent theory of the screened disorder potential at the TI surface and compare it with experiment. Sec. VII uses this theory to calculate the conductivity of surface electrons. Sec. VI briefly discusses how the amplitude of the disorder potential transitions from its large bulk value to its smaller value at the surface. Finally, Sec. VIII discusses an extension of our analysis to the case where the TI surface has a gap introduced by some source of time-reversal symmetry breaking. Where applicable, the major results of each section are summarized at the beginning of the section.

II. ORIGIN OF THE ENHANCED BULK CONDUCTIVITY
As mentioned in the Introduction, randomlypositioned impurities create a disordered Coulomb landscape in the bulk of the TI, which has the effect of reducing the activation energy ∆ relative to what one would naively expect by thinking about flat valence and conduction bands. In this section we explain this idea more fully, focusing first on the case of complete compensation, where the bulk transport can be described using the theory of a completely compensated semiconductor (CCS) 24,25 .
This theory is based on the idea that at K = 1, when almost all donors and acceptors are charged, random spatial fluctuations of the local concentration of impurities result in large fluctuations of charge. Their potential is poorly screened, because of the vanishing average concentration n = N D −N A of electrons, and therefore has huge fluctuations. These fluctuations bend the conduction and valence band edges and in some places bring them to the Fermi level, creating electron and hole puddles that in turn non-linearly screen the random potential. As a result, the amplitude of potential fluctuations is limited by E g /2, so that the ground state, illustrated schematically in Fig. 1, resembles a network of p-n junctions 24,25 . The characteristic size of these p-n junctions is 18 which can be thought of as the correlation length of the random potential. For the typical parameters E g ∼ 0.3 eV, N D = 10 19 cm −3 , and dielectric constant κ = 30, this length scale R ≈ 150 nm N −1/3 D = 4.6 nm. That is, we deal with a very long range potential.
As a result of these long range fluctuations, the resistivity can be dramatically different from the naive expectation based on thinking about flat valence and conduction bands. First, at relatively high temperatures activated conductivity is due to electrons and holes activated from the Fermi level to their corresponding classical percolation levels (classical mobility edges), E e and E h , in the conduction and the valence bands. According to numerical modeling 18 at K = 1, the activation energy ∆ 0.15E g , meaning that E e and E h are substantially closer to the Fermi level µ than to the unperturbed bottom of the conduction band, E c , or ceiling of the valence band, E v (Fig. 1a). (E c and E v are the energies of the conduction and valence bands, respectively, as they would be in the absence of a random potential.) Thus, one can think of the universal small factor ∆/E g ≈ 0.15 as corresponding to a percolation threshold associated with percolation through the potential created by random Coulomb impurities in 3D.
Second, at sufficiently low temperatures electrons and holes can hop (tunnel) directly between puddles, so that activated transport is replaced by VRH. In Ref. 18 we showed that with decreasing temperature the activated resistivity crosses over directly to the Efros-Shklovskii (ES) law 26 where T ES = Ce 2 /k B κξ, e is the electron charge, ξ is the localization length of electron states with energy close to the Fermi level, and C ≈ 4.4 is a numerical coefficient. Together our results for the activated and VRH resistivity established the universal upper limit of the bulk resistivity ρ(T ) for a 3D TI compensated by shallow impurities. In a Ref. 19, we expanded our focus to consider not just the maximum possible bulk resistivity that appears at K = 1, but to address the more practical question of the dependence of the bulk resistivity on the degree of compensation K at 1 − K 1. Indeed, with existing methods of growth of TI samples one cannot get K = 1 exactly, and it is important to know how the results for a CCS, where K = 1, are extended to the case of a strongly compensated semiconductor (SCS), for which 0 < 1 − K 1. For example, one can ask at which value of 1 − K does the activation energy ∆ become twice smaller than at K = 1. For definiteness we consider an n-type SCS, where the concentration of electrons n = N D − N A N D . We model numerically the ground state of such a SCS and its resistivity using algorithms similar to those of Ref. 18. We find that, in agreement with analytic theory 24 , when 1 − K grows the screening of the random potential improves and the correlation length R of the random potential decreases. The amplitude of the random potential decreases as well; hole puddles shrink and eventually vanish; and the chemical potential µ moves up, so that E c − µ decreases. One can say that with increasing (1 − K) screening happens by bending of the conduction band only, while all acceptors remain occupied by electrons and negatively charged. All these changes are illustrated by the transition from a) to b) in Fig. 1.
As a result of these changes with growing 1 − K, the activation energy ∆ decreases. We find that the relation ∆ = 0.3(E c − µ) obtained in Ref. 18 for K = 1 remains valid for 1 − K 1 as well (see Fig. 6 below). [In p-type semiconductors, where K = N D /N A , a similar relationship holds: ∆ = 0.3(µ − E v ).] By K = 0.97 the activation energy ∆ is already several times smaller than at K = 1. This result shows that achieving the maximum bulk resistivity, with ∆ = 0.15E g , is not easy. It also helps to explain the origin of the large scatter in the magnitude of ∆ among different TI samples 15 .
Our prediction that ∆ = 0.3(E c − µ) can in principle be directly compared with experiments in TIs. Indeed, for each K the position of the Fermi level, (E c − µ), can be found via measurements of the concentration of electrons in the surface states using Shubnikov-de-Haas oscillations.
At lower temperatures the activated bulk conduction crosses over to ES VRH. In Sec. IV we study this crossover numerically and also show how T ES , which is correlated with ∆, decreases with 1 − K.
It should be mentioned that these results for the bulk conductivity are also applicable to other narrow gap semiconductors, for example, to InSb. Historically, a large effort was made to make InSb insulating via strong compensation, with the goal of improving the performance of InSb-based photodetectors. Results were again frustrating: the dark resistivity was too small. Our results are in reasonable agreement with transport experimental data for InSb 27,28 .

III. MODEL OF BULK IMPURITIES AND THE DENSITY OF STATES
In order to study numerically the bulk properties of a heavily doped SCS, we introduce a model of the bulk donors and acceptors. In this section we first describe our numerical model and then use it to calculate the position of the Fermi level relative to the band edges as a function of compensation, K, and to evaluate the density of states of impurity states. Our major results are shown below in Figs. 2 and 3. Specifically, we model the bulk as a cube containing a large number of randomly-positioned donors and acceptors. We numerate all donors and acceptors by the index i and use n i = 0 or 1 to denote the number of electrons residing on a donor or acceptor. We also introduce the binary variable f i to discriminate between donors (for which f i = 1) and acceptors (f i = −1). The Hamiltonian of our system is then is the interaction energy between two like-charged impurities at a distance r, and all energies are defined relative to the Fermi level. The first term of Eq. (4) contains the difference between the energies of donors and acceptors, which for the case of shallow impurities is very close to the semiconductor gap E g . The second term of H represents the total interaction energy of charged impurities. Note that Eq. (4) does not include the kinetic energy of electrons and holes in the conduction and valence bands and, therefore, aims only at a description of the low temperature physics of SCS (k B T E c − µ). The form of the interaction law V (r) requires some consideration. For two impurities at a distance r a B , where a B is the effective Bohr radius of impurity states, one can use for V (r) the normal Coulomb interaction V (r) = e 2 /κr. For example, one can consider a pair of empty and distant donors. In such a donor pair one donor shifts the energy of the electron level on the other by V (r) = −e 2 /κr. This classical form for V (r) is good for a lightly doped SCS. In a heavily doped SCS, on the other hand, where a B > N −1/3 D , most impurities have at least one neighbor at distance r < a B , and quantum mechanical averaging over the electron wave function becomes important. (This is why an uncompensated heavily doped semiconductor is a good metal.) For example, a pair of donors cannot create an electron energy state deeper than that of the helium-like ion, which has binding energy 2e 2 /κa B is the binding energy of the shallow donor state. The interaction law V (r) should therefore be "softened" at short distances r < a B to reflect quantum mechanical effects. We model this behavior by continuing to use the classical Hamiltonian of Eq. (4) with a truncated Coulomb potential V (r) = e 2 /κ(r 2 + a 2 B ) 1/2 . Below it is convenient to express energies in units of e 2 N 1/3 D /κ. In these units, a typical TI with band gap 0.3 eV has E g ≈ 30. We unfortunately could not model E g = 30 directly, since in this case the very large correlation length of the random potential, R, leads to large size effects. Instead, we present results for the more modest value E g = 15, for which the size effect requires extrapolation 18 only for K = 1. Results for the smaller E g = 10 are largely identical 18 .
In our numerical simulations, we first randomly place donors and acceptors within the simulation volume; results presented below correspond to M = 20000 donors and 20000K acceptors. We then search for the arrangement of electrons (or equivalently, the set of electron occupation numbers {n i }) that minimizes H, and we use this set to calculate the DOS and the conductivity. We begin our search from the state where all M K acceptors are populated by electrons and negative (n i = 1, q i = −1), and where an equal number of randomly chosen donors are empty and positive (n i = 0, q i = 1), while the remaining M (1 − K) donors are filled and neutral (n i = 1, q i = 0). The charged donors and acceptors in this initial state create a random potential whose magnitude exceeds E g , and as a result the system's energy is well above that of the ground state. In order to bring the system closer to its ground state, we attempt sequentially to transfer electrons from an occupied impurity (either a neutral donor or a negatively charged acceptor) to an unoccupied one (a positively charged donor or a neutral acceptor). If the proposed move lowers the total system energy H, then it is accepted, otherwise it is rejected. To check whether H goes down with each proposed move, for a given set of electron occupation numbers {n i } it is convenient to introduce the single-electron energy state, ε i , at a given impurity i: In the ground state, single electron energies must satisfy the ES criterion for all i, j with n i = 1 and n j = 0. We use our numerical simulation to loop through all pairs of impurity sites i, j and enforce this criterion; if a given pair does not satisfy Eq. (6), then we move the electron from impurity i to j and recalculate all ε i . This process is continued until no single-electron transfers are possible that lower H. The final arrangement of electrons can be called a pseudo-ground state, since higher order stability criteria of the true ground state (corresponding to simultaneously changing three or more electron numbers) are not checked. Such pseudo-ground states are known to describe the properties of real ground states with a high degree of accuracy 24,29 . Results below are obtained at E g = 15 and a B = N bottom with growing 1 − K. At 1 − K > 0.01 this dependence is in reasonable agreement with the prediction of single band theory (which ignores the valence band and acceptors) 24 1. In the latter case the short range Coulomb interaction at distances r N −1/3 D leads to an additional contribution to µ of the same order of magnitude.
The resulting DOS of impurities is shown in Fig. 3 for K = 1 and K = 0.95. g * (ε) is the DOS in the units of (1 + K)N D /(e 2 N 1/3 D /κ) and is normalized to unity. At K = 1, the nearly constant and symmetric DOS between ε = −E g and ε = E g reflects the practically uniform distribution of the random potential from −E g /2 to E g /2 and, correspondingly, of the band edges E c and E v between 0 to E g and between 0 to −E g , respectively (see On the other hand, at K < 1 the DOS of impurity states loses the donor-acceptor symmetry it has at K = 1. As described in Sec. II (see Fig. 1), with growing 1 − K hole puddles are eliminated so that acceptors become disengaged from screening. The acceptor DOS (leftmost peak of Fig. 3) therefore splits from the donor one, which in turn develops two peaks separated by the Fermi level at ε = 0. The large right peak belongs to empty donors, while the small and narrow left peak belongs to occupied donors (electron puddles). These two donor peaks are separated by the ES Coulomb gap. , as a function of electron energy ε calculated from the Fermi level. Results are plotted for K = 0.95 (blue) and K = 1 (red) using Eg = 15. Impurity states with ε < 0 are occupied and those with ε > 0 are empty. At K = 1 the total DOS of impurities has donoracceptor symmetry, which is lost with growing 1 − K.

IV. NUMERICAL MODELING OF THERMALLY ACTIVATED CONDUCTIVITY
In the previous section we described our procedure for finding the energy levels of donor and acceptor impurities in the pseudo-ground state. We now discuss how these results can be used to calculate the bulk conductivity of a SCS, and we present results for the conductivity both in the high-temperature, activated regime and in the low-temperature, VRH regime. Our major results are twofold. First, we find that in the activated regime the activation energy decreases as the chemical potential approaches the conduction band according to ∆ ≈ 0.3(E c − µ) [see Fig. 6]. Second, we study how the characteristic temperature T ES in the VRH regime depends on compensation, and we find that Our process for numerically calculating the resistivity is as follows. Once the energies {ε i } are known (as calculated using the procedure described in Sec. III), we evaluate the resistivity using the approach of the Miller-Abrahams resistor network 24,30 . In this description each pair of impurities i, j is said to be connected by a link where the activation energy ε ij is defined 24 as follows: The resistivity of the system as a whole is found using a percolation approach 24 . Specifically, we find the minimum resistance R c such that if all links with resistance R ij > R c are cut, then there still exists a percolation pathway connecting opposite faces of the simulation volume. This approach captures the exponential dependence of the resistivity on the temperature, and we ignore details of the prefactor. Below we plot the temperature in the dimensionless units T * = 2k B T κ/e 2 N 2/3 D ξ and the resistivity ρ using the dimensionless quantity (ln ρ) * = (ξN 1/3 D /2) ln R c /R 0 . These dimensionless units eliminate any explicit dependence on the localization length ξ.
In Fig. 4 the resulting resistivity is plotted as a function of (T * ) −1/2 over the huge range of temperature 200 > T * > 0.03 for four different values of the compensation K. The resulting linear dependence at 0.3 > T * > 0.03 indicates that at low temperatures the resistivity is well described by the ES law [see Eq. (3)]. The higher temperature range 200 > T * > 1 is plotted separately as a function of 1/T * in Fig. 5. Here the linear slope suggests a well-defined activation energy that depends on the compensation K. At extremely high T * 50, which generally corresponds to unrealistically large temperatures, the conduction is dominated by activation of carriers across the band gap, which is not captured by our model.
Extracting the slope of the curves in Fig. 5 (dashed lines) gives an estimate of the activation energy ∆ as a function of compensation K. Combining this result with the values for the chemical potential µ(K) calculated in Sec. III yields the data shown in Fig. 6, where ∆ is plotted as a function of (E c − µ) for all the studied values of compensation K = 1, 0.99, 0.98, 0.97, 0.96, 0.95. One can see that the equation ∆ 0.3(E c − µ) holds reasonably well for all K in this interval.
So far we have emphasized results that do not explicitly depend on the localization length ξ. In fact, ξ determines the magnitude of T ES , and therefore determines the value of temperature at which the conduction transitions from activated to VRH behavior. We argue now that in a TI ξ is quite large, leading to a prominent role for VRH. To see this, one can imagine an electron with energy close to the Fermi level tunneling from one electron puddle to another, distant one. If such an electron were to tunnel along the straight line connecting the two puddles it would tunnel through high barriers and its wave function would decay sharply, with a decay length ξ a B . However, this straight line does not constitute the path of least action for the tunneling electron. Instead, a tunneling electron can use the same geometrical path as a classical percolating electron, which has energy ∆ above the Fermi level, and thereby avoid large barriers. One can roughly estimate the tunneling decay length by assuming that along such a "percolating" tunneling path the potential energy barriers V are uniformly distributed in the range 0 ≤ V ≤ ∆ and neglecting the additional contribution to the action associated with curvature of this path. Integration over V then gives a localization length ξ ∼ / √ m∆ and k B T ES = 4.4(m∆) 1/2 (e 2 /κ ). For a TI with a B = N The dependence T ES ∝ √ ∆ implies that when ∆ increases ∼ 2.5 times, as in Fig. 6 corresponding to the difference between K = 0.95 and K = 1, the ES temperature T ES increases by ∼ 60%. For a TI with κ = 30 and N D = 10 19 cm −3 , this corresponds to a variation in T ES from 500 to 800 K. The regime of ES VRH in TIs can be studied experimentally, but such a study requires sufficiently thick samples that the bulk conduction provides a larger contribution to the total conductance than the TI surfaces.

V. SELF-CONSISTENT THEORY OF THE SURFACE DISORDER POTENTIAL
In the first part of this paper we showed how the bulk conduction is strongly influenced by the presence of random Coulomb impurities, which produce large bending of the bulk conduction and valence bands. We now turn our attention to the problem of how these same impurities affect the surface transport provided by the Dirac-like surface states. For this problem we adopt the same model of monovalent Coulomb impurities that are randomly distributed throughout the bulk of the TI, and we focus our attention on the case of complete (or nearly-complete) compensation N D = N A ≡ N , where the Fermi level lies within the bulk band gap. As we show below, for determining the properties of the surface one can safely ignore the weak nonlinear screening by electron and hole puddles formed in the bulk (illustrated in Fig. 1).
In this section we present a self-consistent theory for the magnitude of the disorder potential at the TI surface, following Ref. 20. Our primary result is an expression for the amplitude of fluctuations of the electric potential energy, Γ, at the TI surface as a function of the chemical potential, µ, measured relative to the Dirac point. In particular, for µ = 0 we show below that Here α = e 2 /κ s v is the effective fine structure constant, where κ s is the effective dielectric constant at the surface and v is the Dirac velocity. This expression describes screening of the disorder potential via the formation of electron and hole puddles at the TI surface. The characteristic size of these puddles is given by and the corresponding total number of electrons (or holes) per unit area in surface puddles is given by Eqs. (8) -(10) are derived below, along with results corresponding to large µ. Below we also derive a simple relation for the autocorrelation function of the potential at the TI surface, which has an unusually slow decay and can be used to verify the bulk origin of disorder. These results were confirmed by numerical simulation in Ref.

20.
Our primary tool for describing screening of the electric potential is the Thomas-Fermi (TF) approximation, which applies in the limit where the potential φ(r) varies slowly compared to the characteristic Fermi wavelength of electrons at the surface. Specifically, the TF approximation gives where E f (n) = v 4π|n| sgn(n) = (e 2 /ακ s ) 4π|n| sgn(n) is the local Fermi energy and n(r) is the 2D electron concentration at the point r on the surface. The TF approximation is justified whenever α 1, as we show below. In TIs such small α can be seen as the result of the large bulk dielectric constant κ 30. We note here that for describing the properties of the surface state, which exists at a dielectric discontinuity, one should use for the effective dielectric constant κ s the arithmetic mean of the internal and external dielectric constants. If the TI is in vacuum, then κ s = (κ + 1)/2 κ/2.
When the chemical potential is large enough in magnitude that µ 2 e 2 φ 2 , where ... denotes averaging over the TI surface, the relation E f (n) can be linearized to read E f [n(r)] µ+δn(r)/ν(µ). Here δn(r) = n(r)−n 0 is the difference in the electron concentration relative to the state with zero electric potential, n 0 = α 2 κ 2 s µ 2 /(4πe 4 ), and ν(µ) = α 2 κ 2 s |µ|/(2πe 4 ) is the density of states at E f = µ. From this density of states one can define a screening radius r s = κ s /2πe 2 ν = e 2 /α 2 κ s µ that characterizes the distance over which fluctuations in the Coulomb potential are screened by the surface. The TF approximation is valid when the Fermi wavelength λ f ∼ n −1/2 0 ∼ e 2 /ακ s µ is much smaller than r s , which gives the condition α 1. One can understand qualitatively the magnitude of the potential fluctuations, Γ, using the following simple argument. For a given point on the TI surface, one can say that only impurities within a distance R r s contribute to the potential; those impurities at a distance R r s are effectively screened out (one can say that they are screened by their image charges in the "metallic" TI surface). Impurities with R < r s , on the other hand, are essentially unscreened. There are ∼ N r 3 s such impurities, and their net charge is of order Q ∼ e N r 3 s , with a random sign. The absolute value of the potential at the surface is then ∼ Q/κ s r s , so that Γ ∼ eQ/κ s r s ∼ (e 2 N 1/3 /κ s )(N r 3 s ) 1/6 ∼ e 2 N/κ s ν ∼ e 4 N/α 2 κ 3 s |µ|. In order to more accurately derive the value of Γ, one can start by considering the potential created by a single impurity charge +e. When such an impurity charge is placed a distance z from the TI surface (say, above the origin), it creates a potential φ 1 (r; z) that within the TF approximation is given by 31 where J 0 (x) is the zeroth order Bessel function of the first kind. At large z/r s , Eq. (12) can be expanded to give A simple physical derivation of Eq. (13) is based on the notion 32 that for a distant impurity, such that z r s , a surface with screening radius r s effectively plays the role of a metallic surface positioned below the real surface at a distance z = −r s /2. Equation (13) can then be viewed as the sum of the potentials created by the original charge at a distance z above the plane and its opposite image charge at a distance z + r s below the plane, expanded to lowest order in r s /z.
The total potential at the origin is where the index i labels all impurity charges, q i is the sign of impurity i, and r i and z i are the radial and azimuthal coordinates of its position. Under the assumption that all impurity positions are uncorrelated and randomly-distributed throughout the bulk of the TI, the average of φ 2 is given by Here, the quantity 2N d 2 r dz describes the probability that the volume element d 2 r dz contains an impurity charge, and the integration is taken over the semi-infinite volume of the bulk of the TI. The width of the disorder potential at the TI surface, Γ, is defined by Γ 2 = e 2 φ 2 . Inserting Eq. (12) into Eq. (14) and taking the integral then gives Eq. (15) is correct so long as the fluctuations in the Coulomb potential energy are small compared to the chemical potential, or Γ |µ|; this gives the condition written in parentheses.
On the other hand, when |µ| is very small, the fluctuations in the Coulomb potential become large compared to the chemical potential, and one cannot talk about a spatially uniform local density of states ν or screening radius r s . Instead, the Fermi energy has strong spatial variations, and the random potential is screened by the formation of electron and hole puddles at the surface. Nonetheless, one can define an average density of states ν at the surface, which determines, self-consistently, the typical screening radius r s and the magnitude of the potential fluctuations at the TI surface. This value ν can be equated with the thermodynamic density of states of the system, dµ/d n , where n is the overall electron concentration of the surface.
Consider, for example, the case µ = 0, where by symmetry the average value of the potential φ = 0. At any given point r on the surface, the potential φ(r) is the sum of contributions from many individual impurity charges, provided that the characteristic screening radius r s = κ s /2πe 2 ν N −1/3 . This implies that, by the central limit theorem, the value of the potential across the surface is Gaussian-distributed with some variance φ 2 = Γ 2 /e 2 that remains to be calculated. Within the TF approximation the local density of states at the point r is ν[−eφ(r)] = eα 2 κ 2 s |φ(r)|/(2πe 4 ), so that one can calculate the average density of states as This result for ν can be inserted into the first equality of Eq. (15), Γ 2 = e 2 N/κ s ν , to give a self-consistent relation for the amplitude of potential fluctuations 33 . This procedure gives the result first announced at the beginning of this section, Eq. (8). Substituting Eqs. (8) and (16) into the expression for the screening radius, r s = κ s /2πe 2 ν , gives Eq. (9). One can also calculate the total concentration of electrons/holes in surface puddles, n p , implied by this result for Γ 2 . This is done by first inverting the TF relation, Eq. (11), at µ = 0 to give n(φ) = (α 2 κ 2 s /4πe 2 )φ 2 sgn(φ). Integrating this expression for n(φ) weighted by the Gaussian probability distribution for φ gives Substituting the result of Eq. (8) for Γ 2 then gives Eq. (10). One can also combine this result for the residual electron/hole concentration, n p , with the expression for the screening radius, r s , to arrive at an estimate for the number of electrons/holes per puddle: M p ∼ πn p r 2 s ∼ π/16α 2 . Apparently at small α puddles typically contain many electrons/holes, M p 1. Our primarily results, outlined in Eqs. (8) - (10), are valid within the TF approximation so long as the typical Fermi wavelength, λ f ∼ e 2 /ακ s Γ, is much smaller than the typical screening radius, r s ∼ e 2 /α 2 κ s Γ, which again gives the condition α 1.
As we mentioned above, at µ = 0 the screening radius r s describes the characteristic size of electron or hole puddles at the TI surface. More generally, r s plays the role of a length scale over which potential fluctuations at the surface are correlated. Such correlations can be discussed in a quantitative way by defining the potential auto-correlation function: where ... r denotes averaging over the spatial coordinate r , and where by symmetry the correlation function depends on |r| = r only. In the remainder of this section we derive approximate analytical results for C(r), and show that spatial correlations in the potential have an unusually slow decay. At r = 0, Eq. (17) reproduces the definition of φ 2 , so that C(0) = Γ 2 /e 2 . At small enough distances that r r s , one can expect that the value of C(r) is determined primarily by unscreened impurities that are within a distance r s from the surface, as explained above during the derivation of Γ 2 . On the other hand, at r r s correlations are produced primarily by impurities that are relatively far from the surface, as can be seen from the following scaling argument. Consider two surface points separated by a distance r r s . One can imagine drawing a cube of size r that extends into the bulk of the TI and which contains the two surface points on opposite edges of one of its faces. Such a cube contains ∼ N r 3 impurities, and has a net impurity charge with magnitude q ∼ e √ N r 3 and random sign. These impurity charges are located at a mean distance ∼ r r s above the surface and, therefore, by Eq. (13), contribute a net potential ∼ qr s /κ s r 2 ∼ (e/κ s ) N r 2 s /r to both surface points. The square of this potential roughly gives the autocorrelation of the potential, C(r) ∼ e 2 N r 2 s /κ 2 s r. A more careful expression for C(r) can be derived by writing similar to Eq. (14). Inserting the asymptotic expression of Eq. (13) for φ 1 and evaluating the integral gives This result was also confirmed by numerical simulation in Ref. 20. Eq. (19) implies an unusually slow decay of potential correlations at the surface, which, as explained above, arises from long-range fluctuations of the potential created by deep bulk impurities. This behavior can be contrasted with the much faster decay of C(r) that would result from a two-dimensional (2D) distribution of Coulomb impurities at a distance d from the surface: 34 C(r) ∼ e 2 n i dr 2 s /κ 2 s r 3 , where n i is the 2D impurity concentration. Thus, by studying C(r) experimentally by scanning tunneling microscopy, one can discriminate between disorder by bulk impurities and disorder by impurities located in a layer close to the surface.
We now discuss the magnitude of Γ and r s implied by these expressions for typical TIs, which generally have an impurity concentration N ∼ 10 19 cm −3 . Typical values of the Dirac velocity and fine structure constant for TIs can be taken from Ref. 21, which reports v = 1.3 eVÅ and estimates α = 0.24. Using these parameters gives Γ ∼ 30 meV and r s ∼ 20 nm at the Dirac point, µ = 0. At large |µ| 30 meV, both Γ 2 and r s decay as 1/|µ|.
As mentioned in the Introduction, the theory presented in this section provides a good description of the recent experimental results of Ref. 21, where the random potential at the surface of the 3D TIs Bi 2 Se 3 and Bi 2 Te 3 was studied using a scanning tunneling microscope 20 . Indeed, in these experiments it was found that the electric potential at the surface was well-characterized by a Gaussian distribution with a standard deviation Γ ∼ 10 -20 meV, and the characteristic length scale of potential fluctuations was estimated as r s ∼ 20 -30 nm. One can compare these measurements to our theoretical predictions by using the parameters listed above and inserting the measured chemical potential µ ∼ 100 meV into Eq. (15). This procedure gives Γ ∼ 18 meV, and the corresponding screening radius r s ∼ 5 nm, so that our theory is indeed in reasonably good agreement with experiment. Further, Ref. 21 found that the disorder potential at the surface was not correlated with the position of surface impurities, indicating that the surface disorder potential originates primarily from impurities deep below the TI surface, as we have described.
Throughout this section, we have worked within the assumption that bulk impurities are completely ionized, or in other words that there is no screening by conduction band electrons or valence band holes in the bulk. Such an assumption is valid when the chemical potential resides in the middle of a large bulk band gap. In this case donors or acceptors can only be neutralized by very large band bending discussed in Sec. II (see Fig. 1). Such fluctuations take place over a long length scale R that scales as the square of the distance between the Fermi level and the nearest band edge [see Eq. (2)] and is typically on the order of hundreds of nanometers for typical TIs 18 . On the other hand, near the surface of the TI the potential fluctuations are screened much more effectively and over a much shorter distance, r s , by the (ungapped) surface states. As shown above, r s is typically 20 nm, and the amplitude of surface potential fluctuations Γ ∼ 30 meV E g ∼ 300 meV. One can therefore safely assume that near the surface there is no large band bending and one can indeed treat bulk impurities as completely ionized.
The effect of bulk screening should appear only in the long-range behavior of the correlation function, r R, where the 1/r decay of C(r) is truncated and, as one can show, is replaced with C(r) ∼ e 2 N Rr 2 s /κ 2 s r 2 . Finally, we note that our theory ignores the possibility of screening by material outside the TI. For example, if the TI is placed next to a metal electrode or an ionic liquid 35 , then this external material can screen the large potential fluctuations created by the bulk, thereby decreasing Γ and r s .

VI. FROM SURFACE TO BULK
In Sec. II we showed that deep within the bulk of the TI the disorder potential has large fluctuations of order Γ ∼ E g . On the other hand, in Sec. V we showed that at the TI surface the disorder potential has a much smaller amplitude, Γ ∼ (e 2 N 1/3 /κα 2/3 ). In this section we elaborate briefly on the crossover between these two results, or in other words we describe how the amplitude of potential fluctuations grow as one moves from the surface of the TI into the bulk.
Generally speaking, as one moves a distance z > 0 into the bulk of the TI, the amplitude of the disorder potential increases in magnitude. In order to see quantitatively how Γ grows as a function of z, one can assume, for the moment, that the TI surface is equivalent to a perfect metallic plane. In this case, each impurity at position (r , z ) has a corresponding image charge at (r , −z ), and the total potential at (0, z) is equal to the sum of the potentials created by the original impurity and its image. One can calculate Γ 2 (z) by averaging the square of this potential over all possible positions of the impurity charge [as in Eq. (14)]. This calculation gives Γ 2 (z) = 8πN e 4 z/κ 2 . That is, Γ 2 (z) grows linearly with the distance z from the TI surface. This growth continues until Γ becomes large enough that Γ 2 (z) = (E g /2) 2 , at which point electron and hole puddles begin to form in the bulk and one arrives at the bulk screening picture described in Ref. 18. This distance corresponds to z = R/4; at smaller z the potential fluctuations are small enough that practically all donors and acceptors are charged.
One can now recall that the TI surface is not perfectly metallic, and that its screening length r s is finite, so that Γ 2 (z) should be somewhat larger. In fact, at z r s one can still use the formula above for Γ 2 (z) by introducing a small modification allowing for the fact that the metallic surface is effectively shifted to the position z = −r s /2 (as discussed in Sec. V). Making this adjustment gives Γ 2 (z) = 8πN e 4 (z + r s /2)/κ 2 at z r s , which does not significantly alter our conclusions.

VII. SURFACE CONDUCTIVITY
We now turn our attention to the problem of how the 3D-distributed Coulomb impurities within the TI bulk affect the surface conductivity. As discussed at the beginning of Sec. V, we limit our consideration to the case where the Fermi level resides within the bulk band gap, where one can safely assume that all relevant bulk impurities are ionized.
Our primary result is an expression for the electron conductivity σ of the surface as a function of the average 2D surface electron concentration n. In particular, when n n p , where n p is the typical puddle concentration at µ = 0 [see Eq. (10)], we find that the conductivity is given by where e 2 /h is the conductance quantum. At much smaller electron concentrations, n n p , the conductivity saturates at a value σ min , which we estimate as To derive these results, we first note that in the limit of large chemical potential µ, where the electron density is only weakly modulated by the disorder potential, one can show using the Boltzmann kinetic equation that for electrons with a massless Dirac spectrum the conductivity is given by [36][37][38][39] Here τ is the momentum relaxation time. In the limit of zero temperature, the scattering rate 1/τ can be found by integrating the squared scattering potential produced by a given impurity over all impurities and over all scattering angles. More simply, one can arrive at an expression for 1/τ by taking the result for the scattering rate of a 2D layer of impurities with concentration n i at distance z [for example, Eq. (38) of Ref. 37], replacing n i with 2N dz, and then integrating over all planes z containing impurities. This procedure gives (1−cos 2 θ).
Evaluating the integral of Eq. (23) at small α gives Inserting this result for τ into Eq. (22) and substituting µ = e 2 k f /ακ s and k f = √ 4πn yields the result for conductivity announced at the beginning of the section, Eq. (20).
Equation (20) can be contrasted with the widely-used result for the 2D model of charge impurities 36,37,39,40 , for which the conductivity is linearly proportional to the electron density: σ/(e 2 /h) ∼ (1/α 2 )(n/n i ). This difference can be understood conceptually by noting that, for large angle scattering, only those impurities at a distance smaller than the Fermi wavelength, λ f ∼ n −1/2 , contribute significantly to scattering. One can therefore define, roughly speaking, an effective 2D concentration of scattering impurities as N λ f ∼ N/n 1/2 . Inserting N/n 1/2 for n i gives σ ∝ (1/α 2 )(n 3/2 /N ), similar to Eq. (20). The remaining factor 1/ ln(1/α) in Eq. (20) is related to low-angle scattering by distant impurities with z λ f . So far we are unaware of any transport data for TIs that shows σ ∝ n 3/2 . Recent conductivity measurements on ultra-thin TIs (with thickness ∼ 10 nm λ f ) suggest 41 σ ∝ n, consistent with the 2D model of impurities.
Our 3D model also produces a distinct result for the minimum conductivity σ min that appears in the limit of small average electron concentration. At small enough chemical potential that µ e 2 N 1/3 /α 2/3 κ s , the surface breaks into electron and hole puddles, and one can think that the effective carrier concentration saturates at ∼ n p [see Eq. (10)]. An estimate of σ min can therefore be obtained by setting n ∼ n p in Eq. (20), which gives 20 the result of Eq. (21). 2D models of disorder impurities also produce a minimum conductivity that is independent of the impurity concentration, but which has a different dependence on α. Specifically, at small α such models give 36,42 σ min ∼ (e 2 /h) ln(1/α). Our model suggests a minimum conductivity that is larger by a factor ∼ [α ln 2 (1/α)] −1 .

VIII. TI SURFACE WITH A GAP
In Secs. V and VI, we discussed the disorder potential created by Coulomb impurities at a gapless TI surface, whose massless spectrum is protected by time-reversal symmetry. On the other hand, one can open a gap at the TI surface if one introduces some source of timereversal symmetry breaking, such as an external magnetic field 4,43,44 , proximity to a magnetic material or magnetic impurities 45,46 , the proximity effect from an adjacent superconductor 47 , or electron tunneling between two nearby TI surfaces 48,49 (see also the review of Ref. 6). The resulting gapped spectrum is illustrated schematically in Fig. 7.
In this final section we briefly discuss how the presence of a gap with magnitude U affects the disorder potential at the surface and the mid-gap density of states. We focus our discussion around the case where the chemical potential µ = 0, which roughly corresponds to the largest disorder potential and the minimum in the thermodynamic density of states. Again, we limit our consideration to the case where the TI is sufficiently thick that one can describe impurities as three-dimensionally distributed.
In the absence of a gap, U = 0, the disorder potential is well-described by the results of the previous section.
In particular, the disorder potential width Γ = Γ 0 ≡ (2 1/6 √ π/α 2/3 )(e 2 N 1/3 /κ s ) [see Eq. (8)] and the average density of states ν = ν 0 ≡ (α 4/3 /2 1/3 π)(κ s N 1/3 /e 2 ) [see Eq. (16)]. If the gap U is small enough that U Γ 0 , then the disorder potential at the surface is essentially unaffected by the gap, since local fluctuations in the Fermi level are much larger than U . For example, if the gapless surface spectrum is replaced with a "massive" dispersion relation, as plotted in Fig. 7, then one can estimate the first-order effect of the gap by carrying out the same self-consistent procedure outlined in Sec. V. In particular, the gapped dispersion relation of Eq. (25) has a corresponding density of states where Θ(x) is the Heaviside step function. At U/Γ 0 1, one can assume a Gaussian distribution of the Coulomb potential φ with some unknown variance Γ 2 , integrate this distribution over φ multiplied by ν(−eφ; U ) to produce the thermodynamic density of states ν , and then use the self-consistency relation Γ 2 = e 2 N/κ s ν to arrive at a value for Γ [see Eq. (15)]. Expanding the result of this procedure for small U/Γ 0 gives for the disorder potential width a slightly enhanced value Similarly, the thermodynamic density of states is slightly depleted: On the other hand, if U is much larger than Γ 0 , then the surface screens poorly and the disorder potential grows. In this case screening of the disorder potential by the surface happens only nonlinearly, through the formation of electron and hole puddles at locations where the magnitude of the Coulomb potential energy reaches the gap energy U/2. This is similar to the bulk nonlinear screening discussed in Sec. II, and naturally produces Γ(U ) ∼ U . The typical correlation length of the disorder potential at the surface (the nonlinear screening length) is given by as in Eq. (2), with E g → U .
One can estimate the corresponding concentration of electrons/holes in surface puddles, n p , by noting that a square area of size R 2 U at the surface should contain enough electrons/holes to neutralize the net charge of Coulomb impurities in the adjacent cubic volume R 3 U of the TI bulk. This gives n p R 2 U ∼ N R 3 U , or in other words The corresponding thermodynamic density of states can be estimated by noting that when the chemical potential µ is raised by an amount ∼ U/2, surface hole puddles should dry up and be replaced by a correspondingly increased number of electron puddles. This suggests ν = dµ/dn ∼ n p /U , which gives ν ∼ e 2 N κ s U 2 .
Notice that if the gap U is reduced to the point where U ∼ Γ 0 , then ν → ν 0 , as can be seen by comparing Eq. (31) with Eqs. (8) and (16). Of course, these estimates assume that the surface gap U is smaller than the bulk band gap E g , and consequently that R U R, so that impurities near the surface are not screened by bending of the bulk bands. If the surface gap U is larger than E g , then the disorder potential variance is truncated at Γ(U ) ∼ E g due to bulk screening.
We note that Eqs. (30) and (31) were first derived in Ref. 50 in the context of semiconductor heterostructures in a transverse magnetic field, where a 2D electron gas experiences disorder from adjacent 3D impurities and the gap U in the kinetic energy spectrum is provided by the Landau level spacing ω c . These authors also showed how the disorder potential is reduced and the density of states increased as the chemical potential µ is increased from zero 50 . Specifically, Γ ∼ U − 2µ and ν ∼ e 2 N/κ s (U − 2µ) 2 , provided that U − 2µ e 2 N 1/3 /κ s . Of course, for TIs the effect of a transverse magnetic field goes beyond simply opening a single gap at the Dirac point 43,44 . We do not consider here the full problem of screening of Coulomb impurities in the presence of a magnetic field, but in principle this problem can be dealt with along the lines of Ref. 50.