Length scale of puddle formation in compensation-doped semiconductors and topological insulators

TL;DR

This paper investigates the length scale of electron-hole puddle formation in compensation-doped semiconductors and topological insulators, revealing a smaller-than-expected growth with the band gap and implications for TI thin film quality.

Contribution

It provides a numerical and scaling analysis of puddle length scales, challenging previous predictions and offering a method to estimate film thickness for improved topological insulator quality.

Findings

Puddle length scale grows approximately linearly with  for accessible parameters.

Numerical results suggest a smaller length scale than the ^2 scaling predicted by simple models.

Extrapolation indicates a ^2 /   scaling for large .

Abstract

In most semiconductors and insulators the presence of a small density of charged impurities cannot be avoided, but their effect can be reduced by compensation doping, i.e. by introducing defects of opposite charge. Screening in such a system leads to the formation of electron-hole puddles, which dominate bulk transport, as first recognized by Efros and Shklovskii. Metallic surface states of topological insulators (TI) contribute extra screening channels, suppressing puddles. We investigate the typical length $l_p$, which determines the distance between puddles and the suppression of puddle formation close to metallic surfaces in the limit where the gap $\Delta$ is much larger than the typical Coulomb energy $E_c$ of neighboring dopants, $\Delta \gg E_c$. In particular, this is relevant for three dimensional Bi-based topological insulators, where $\Delta/E_c \sim 100$. Scaling arguments predict $l_p \sim (\Delta/E_c)^2$. In contrast, we find numerically that $l_p$ is much smaller and grows in an extended crossover regime approximately linearly with $\Delta/E_c$ for numerically accessible values, $\Delta/E_c \lesssim 35$. We show how a quantitative scaling argument can be used to extrapolate to larger $\Delta/E_c$, where $l_p \sim (\Delta/E_c)^2/\ln(\Delta/E_c)$. Our results can be used to predict a characteristic thickness of TI thin films, below which the sample quality is strongly enhanced.

Equations

$$H=H_n+H_C=\frac{\Delta }{2}\sum _i{f}_i{n}_i+\frac{1}{2}\sum _{i\ne j}{V}_{ij}{q}_i{q}_j$$

$$q_i=\frac{f_i+1}{2}-n_i.$$

$$V_{ij}=\frac{e^2}{4\pi \epsilon {\epsilon }_0\sqrt{|{\mathbf{r}}_i-{\mathbf{r}}_j{|}^2+a_B^2}}=\frac{E_c}{\sqrt{|{\mathbf{x}}_i-{\mathbf{x}}_j{|}^2+1^2}}.$$

$$E_c=\frac{e^2{N}^{1/3}}{4\pi \epsilon {\epsilon }_0}$$

$${ϵ}_j=\frac{\Delta }{2}{f}_j-{\varphi }_j=\frac{\Delta }{2}{f}_j-\sum _{i\ne j}{V}_{ij}{q}_i.$$

$$\Delta E_{(\alpha ,\beta )}={ϵ}_{\beta }-{ϵ}_{\alpha }-V_{\alpha \beta }>0$$

$$\text{centering}{R}_g\sim {\Delta }^2\text{ and }{\rho }_g\sim {\Delta }^{-3}.\text{@add@centering}$$

$$C_{\rho \rho }(\mathbf{r},{\mathbf{r}}^{\prime })=⟨\rho (\mathbf{r})\rho ({\mathbf{r}}^{\prime })⟩=Q_0\delta (\mathbf{r}-{\mathbf{r}}^{\prime })+C_{\rho \rho }^{\text{nl}}(\mathbf{r}-{\mathbf{r}}^{\prime }),$$

$$Q_s(\Delta ,r)=4\pi {\int }_0^{r}ds{s}^2{C}_{\rho \rho }^{\text{nl}}(\Delta ,s).$$

$${\ell }_P\sim {\Delta }^{\gamma },\gamma \approx 1.1\pm 0.1$$

$$p_{L_z/2}^0(\Phi )=\frac{1}{\sigma (L_z)\sqrt{2\pi }}\exp \left[-\frac{1}{2}{\left(\frac{\Phi }{\sigma (L_z)}\right)}^2\right]$$

$$\sigma (L_z)$$

$$\approx 2.41\sqrt{L_z}-\frac{7.75}{\sqrt{L_z}}\text{for}{L}_z\gg 1$$

$$p_0(|\Phi |>\Delta /2)$$

$$\sigma (L_c^0{)}^2$$

$$L_c^0\sim \frac{{\Delta }^2}{46.5\ln [1/n_0]}\approx \frac{{\Delta }^2}{139\ln [\Delta ]}$$

$$n_0(d)=⟨\sum _i\delta (d-z_i){\delta }_{q_i,0}⟩,$$

$$Q_2(\Delta ,{\ell }_P)$$

$$n_0(L_z/2)$$

$$n_0({\ell }_s)$$

$${\ell }_s>c{\left(\frac{1}{\alpha }\right)}^{4/3}\text{for }⟨{\mu }^S⟩=0$$

$$\frac{⟨H_{\mathrm{C}}⟩}{V}$$

$$\frac{⟨H_{\mathrm{n}}⟩}{V}=n_{0}N\Delta$$

$$Q_0=\frac{⟨Q⟩}{V}$$

$$⟨{\varphi }^2⟩$$

$$C_{\rho \rho }^{\text{nl}}(\Delta ,s)={\Delta }^{-\beta }\overline{C_{\rho \rho }^{\text{nl}}}(s/{\Delta }^{\gamma }).$$

$${\Delta }^0$$

$${\Delta }^2$$

$$V^m({\mathbf{x}}_i,\mathbf{x})=V[\mathbf{x}-(x_i^1,x_i^2,x_i^3)]-V[\mathbf{x}-(x_i^1,x_i^2,-x_i^3)]$$

$$V^s({\mathbf{x}}_i,\mathbf{x})$$

$${\ell }_s\sim {\Delta }^{\gamma },\gamma \approx 1.1\pm 0.2$$

Full text

I Introduction

The formation of locally conducting puddles is a phenomenon caused by charged Coulomb disorder in insulators, semiconductors and Dirac-matter like graphene, topological surface states or Weyl semimetals. Efros and Shklovskii Efros and Shklovskii (1984) predicted that puddle formation is, in three dimensions, an unavoidable consequence of the long-range nature of the Coulomb interaction. Puddles are formed to screen large potential fluctuations exceeding the size of the gap $\Delta$.

In graphene it has been shown both theoretically and experimentally that puddles are necessary to understand most transport experiments close to charge neutrality Hwang et al. (2007); Adam et al. (2007); J. Martin et al. (2008); Zhang Yuanbo et al. (2009). Recently, they have been observed for the first time in the bulk of a three-dimensional topological insulator Borgwardt et al. (2016); Rischau et al. (2016). These materials, from the class of the Bi2-xSbxTe3-ySey compounds Ren et al. (2011), are almost perfectly compensated semiconductors with a band gap of order $250-300$ meV almost 2 orders of magnitude larger than the typical Coulomb energy $E_{c}$ of neighboring dopants Borgwardt et al. (2016). The relatively high density ($>10^{19}\rm{cm}^{-3}$) of dopants implies a strongly fluctuating Coulomb potential in the bulk. This leads to band bending and eventually to the formation of electron and hole puddles Skinner et al. (2012, 2013). The additional surface states in the topological materials induce an additional screening channel close to the surface. Here surface puddles form Skinner et al. (2013); Skinner and Shklovskii (2013) which are akin to puddles that form in graphene on a substrate which has charged impurities.

As the puddles are separated by insulating regions, they do not directly contribute to the DC conductivity. However, they do contribute to the optical conductivity at finite frequencies, which has been used to detect their presence and to measure the effective charge density in conducting regions Borgwardt et al. (2016). We have shown that screening from thermal excitations can efficiently suppress puddle formation leading to a characteristic temperature dependence of the optical response Borgwardt et al. (2016). Furthermore, in similar compounds a giant negative magnetoresistance was found experimentally and explained by merging of puddles driven by the Zeeman effect Breunig et al. (2017).

Surface puddles and puddles in graphene can be observed directly in real space by scanning tunneling microscopy (STM) J. Martin et al. (2008); Zhang Yuanbo et al. (2009); H. Beidenkopf et al. (2011). From the two-dimensional STM maps, the size of the potential fluctuations and the corresponding length scale can be directly read off. These agree well with theoretical results, where these quantities are calculated self-consistently Hwang et al. (2007); Adam et al. (2007); Skinner and Shklovskii (2013). However, nothing is known experimentally about the length scales of puddles in the bulk and the effect of surface screening on the bulk puddle formation.

In the following we demonstrate numerically that the length scales governing the distance of puddles, the suppression of (bulk) puddles close to surfaces of TIs, and the suppression of puddles in thin films grow much slower with $\Delta/E_{c}$ than expected from scaling arguments. First we introduce the model and consider the scaling behavior of the charge-charge correlation function. We show numerical results for the bulk, and demonstrate that the simple scaling theory fails. Then we additionally take into account the gapless surface states which provide an extra screening channel. The length scales governing the size of puddles on the surface is different, and independent of the bulk band gap Skinner and Shklovskii (2013). The bulk length describes, however, the size of a region where surface screening suppresses the formation of bulk puddles and is therefore important to understand the properties of thin topological insulator samples. We use scaling arguments to extrapolate our numerical results for $\Delta/E_{c}\lesssim 35$ to the experimentally relevant regime of $\Delta/E_{c}\sim 100$.

II Model and Simulations

Bi-based topological insulators typically have a very large dielectric constant $\varepsilon\approx 200$. Electron binding energies are therefore small. Thus, the bare energies of the dopants are located very close to the band edges and can be approximated by $+\Delta/2$ for the donors and $-\Delta/2$ for the acceptors. To model the non-linear screening of randomly placed charged impurities in such a system we use a simple classical model Basylko et al. (2000); Skinner et al. (2012, 2013); Borgwardt et al. (2016):

[TABLE]

where $f_{i}=\pm 1$ are random numbers with $f_{i}=+1$ for a donor states and $f_{i}=-1$ for the acceptor state at position $\bm{r}_{i}$. $V_{ij}$ denotes the Coulomb interaction between the dopants at positions $\bm{r}_{i}$ and $\bm{r}_{j}$. $n_{i}\in\{0,1\}$ denotes the electronic occupation of the $i-$th dopant and is determined by minimizing the Hamiltonian. It is related to its charge $q_{i}$ (in units of $|e|$ where $e$ is the electron charge) by

[TABLE]

A donor (acceptor) in its ground state is characterized by $f_{i}=1$, $n_{i}=0$ and $q_{i}=1$ ($f_{i}=-1$, $n_{i}=1$, $q_{i}=-1$). Somewhat counter-intuitively, screening occurs when the Coulomb interaction drives donors or acceptors into a neutral states with $q_{i}=0$. Several neutral donor states close by form an electron puddle, while neighboring neutral acceptor states form hole puddles. The Coulomb energy is modeled by

[TABLE]

Here the short-distance cutoff $a_{B}=\tfrac{4\pi\varepsilon_{0}\varepsilon}{m^{*}e^{2}}$ was introduced by Skinner et. al. Skinner et al. (2012, 2013) to take into account that the wave function of the bound state is smeared over a length scale set by the effective Bohr radius of the impurity state. Skinner et. al. Skinner et al. (2012, 2013) argued that due to the large dielectric constant in Bi based topological insulators, $a_{B}$ is large and of similar size as the typical distance of dopants. We use $a_{B}=N^{-1/3}$ where $N=N_{A}=N_{D}$ is the density of dopants where we assume a perfectly compensated system where the density of donors equals the density of acceptors $N_{A}=N_{D}$. For the last equality in Eq. (3) we expressed all distance in units of the average dopant distance $N^{-1/3}$. Here

[TABLE]

is the typical energy scale describing the Coulomb interaction of neighboring dopants. A large $\varepsilon\sim 200$ leads to a small energy scale $E_{c}\sim 3.3\,\rm{meV}\sim 40\,K$, (assuming a typical density $N=10^{20}\,\rm{cm}^{-3}$) about 2 orders of magnitude smaller than typical band gaps $\Delta$. Indeed, in Ref. Borgwardt et al. (2016) we used the temperature dependence of the optical response to determine $E_{c}$ and found $\Delta/E_{c}\approx 150$, similar parameters have also been found in Ref. Breunig et al. (2017). In the following we assume $T\ll E_{c}\ll\Delta$ and consider therefore only properties at $T=0$.

The model (1) describes how donor and acceptor states interact with each other. It does not include the states in the electronic bands. This turns out Borgwardt et al. (2016) to be well justified in the limit $\Delta/E_{c}\gg 1$ as the density of the relevant electronic states is much smaller than the density of dopants in this limit.

To find the true ground state of the model in Eq. (1) is an exponentially hard problem, but there is an algorithm to find an approximate ground state, called a pseudo ground state, in polynomial time Efros and Shklovskii (1984); Skinner et al. (2012). The physical properties of a pseudo ground state are expected to be indistinguishable from that of the true ground state. The single particle energies are defined as

[TABLE]

In a pseudo ground state

[TABLE]

has to be fulfilled for all pairs with $n_{\beta}=0,n_{\alpha}=1$. This state can be reached by exchanging electrons between states where this condition is not met. The algorithm is described in detail in Refs. Skinner et al. (2012, 2013). Simulations are performed in a cubic volume $V=L^{3}$ with periodic boundary conditions with $2L^{3}$ dopants, typically we use $L=50$ or $L=60$ corresponding to $250000$ or $432000$ dopants. Numerical results shown below are averaged over $200-800$ disorder realizations, i. e. random configurations of the dopant positions. We have checked Lux (2016) that our code reproduces published results (e.g. on the Coulomb gap in the density of states) from other groups Skinner et al. (2012) on the same model in all quantitative details. In the following we use dimensionless units where all length are measured in units of $N^{-1/3}$, and all energies are measured in units of $E_{c}$ defined in Eq. (4). In these units the only free parameter of our model is $\Delta$ besides the (dimensionless) system sizes considered in Sec. V.

III Length scales and scaling

One of the main questions that we will address is the following: What is the typical distance between electron and hole puddles or, equivalently, on what length scale does the potential typically change by an amount of $\Delta$? It turns out that this length scale also characterizes screening properties on average, as discussed in more detail below.

A simple scaling argument by Efros and Shklovskii Shklovskii and Efros (1972) suggests that the corresponding length scales as $R_{g}\sim\Delta^{2}$. The argument is as follows: in a volume of size $V\sim R^{3}$ there are on average $NR^{3}$ positive and negative charges where $N$ is the density of dopants. But these two numbers are not exactly equal, instead the typical charge of the region is (in the uncorrelated state) $Q_{R}\sim\pm\sqrt{NR^{3}}$. This implies a typical potential of order $\phi_{R}\sim Q_{R}/R\sim\sqrt{R}$ within that region. The fact that this potential diverges for $R\to\infty$ shows that this situation is unstable and the huge potential fluctuations have to be screened. The potential can be screened when the Coulomb potential is sufficiently strong to change the charging state of the dopant. This is possible for $\phi\sim\pm\Delta/2$. Using that $\phi\sim\sqrt{R}$, this strongly suggests that the typical length scale $R_{g}$, describing both the screening length and the length scale where the potential changes by $\pm\Delta$, is proportional to $\Delta^{2}$. Accordingly, the typical charge density in a volume $V=R_{g}^{3}$ is $\rho_{g}\sim Q_{R_{g}}/V\sim\sqrt{R_{g}^{3}}/R_{g}^{3}=R_{g}^{-3/2}\sim 1/\Delta^{3}$. To summarize, this scaling argument suggests

[TABLE]

Restoring dimensionfull units, these equations read $R_{g}\sim N^{-1/3}(\Delta/E_{c})^{2}$ and $\rho_{g}\sim\pm eN(E_{c}/\Delta)^{3}$. We will show below that our numerical results for $\Delta/E_{c}\lesssim 35$ show a much slower growth of length scales with $\Delta/E_{c}$. We will attribute this to a huge crossover regime and the presence of logarithmic corrections obtained from a refined version of the scaling argument in section V.

In Fig. 1 we compare several one-dimensional cuts of the potential in the uncorrelated state (dashed lines) and in the correlated ground state (solid lines) obtained from numerical simulations. The potential fluctuations of the uncorrelated state are much larger than $\pm\Delta/2$ (shown here for $L=50$) triggering screening. For the correlated ground state, in contrast, the potential fluctuations are strongly reduced and lie within the band gap. Puddles are formed in the tiny regions (shaded in gray) where $|\phi|$ slightly exceeds $\Delta/2$. One finds that in these regions $|\phi|-\Delta/2\sim E_{c}$.

The $3d$ plots in panel e) and f) show directly the puddles. Neutral donors constitute electron puddles and are colored in orange, while neutral acceptor are part of hole puddles, colored in blue. The snapshots and the cuts suggest that not only a single, but several length scales govern puddle formation Baranovskii S. D. et al. (1984); Lee et al. (1999); Chen and Skinner (2016). The short one, ${\ell_{P}}$ governs the closest distance of puddles and rapid fluctuations of the potential as shown in panel b). Much longer length scales govern the formation of lengthy, anisotropic cluster structures and also regions without puddles, see panel d).

To obtain more quantitative results, one can study the statistical properties of either the potential $\phi(\bm{r})$ or directly of the charge distribution $\rho(\bm{r})$, since both are related by the Poisson equation $\nabla^{2}\phi=-\rho$ (up to the short-distance cutoff $a_{B}$ introduced above). In the following, we will mainly discuss the charge-charge correlation function $C_{\rho\rho}$. We split this into a local part $\sim\delta(\bm{r}-\bm{r}^{\prime})$ and a non-local part $C^{\textrm{nl}}_{\rho\rho}$:

[TABLE]

where we used charge neutrality $\langle\rho\rangle=0$. Here and in the following the expectation value $\langle\cdot\rangle$ denotes a disorder average. After disorder averaging all correlation functions only depend on the distance $r=|\bm{r}-\bm{r}^{\prime}|$. Thanks to charge neutrality we know that $\int d^{3}\bm{r}\,C^{\textrm{nl}}_{\rho\rho}(\bm{r})=-Q_{0}$. The weight of the $\delta$-peak $Q_{0}$ corresponds to $2N(1-n_{0})$ where $n_{0}$ is the fraction of neutral dopants.

IV Screening in the bulk

Screening in insulating charged Coulomb systems is a highly non-local and non-linear mechanism. Early work by Baranovskii, Shklovskii and Efros Baranovskii S. D. et al. (1984) (see also a lucid discussion in Ref. Lee et al. (1999)) pointed out that adding a single charge can trigger an avalanche of discrete changes of the charge of dopants not only in the neighborhood of the charge but also at large distances. This appears to be a highly anisotropic, non-local (and perhaps fractal Chen and Skinner (2016)) process. The change of the potential at larger distances is random in sign but does not decay rapidly. In contrast to a metal, there is therefore no true screening (as is also obvious from the fact that the system is characterized by a Coulomb gap). In the following we will not track these changes but focus on the shorter length scale ${\ell_{P}}$ which governs the impurity-averaged charge correlations $C^{\textrm{nl}}_{\rho\rho}(\Delta,s)$, but also controls the typical ‘nearest’ distance of oppositely charged puddles. Later we will argue that the same length scale also governs the impact of metallic surface states on puddle formation.

Instead of studying directly the charge-charge correlation function $C^{\textrm{nl}}_{\rho\rho}(\Delta,s)$, we find it more convenient to investigate the distance dependence of the ‘screening charge’ defined by

[TABLE]

The advantage of this quantity is that it has a direct physical interpretation: it describes the charge accumulated – on average – around a dopant within the radius $r$ multiplied with the charge of that dopant and the density of dopants. As negative charges accumulate around a positive charge and vice versa, the screening charge is always negative. Total charge neutrality requires that around a positive (negative) charge exactly the charge $-1$ ($+1$) accumulates for $r\to\infty$. As neutral dopants do not contribute, one therefore obtains $Q_{s}(\Delta,r\to\infty)=-2N(1-n_{0})=-Q_{0}$. This also follows directly by integrating Eq. (8) over $\bm{r}$ in a charge-neutral system. In our simulations we use boxes of size $L$ with periodic boundary conditions. For $r>L/2$ we therefore have to replace in the integral in Eq. (9) the factor $4\pi s^{2}$ by $W(s)=\int\delta(s-|r|)d^{3}r$. This does not affect the scaling plots discussed below but is useful to check for overall charge neutrality.

We show numerical results for $Q_{s}$ in Fig. (2). On a rather short length scale (see inset) the screening charge reaches the value $-Q_{0}=-2N(1-n_{0})\approx-2$ (the plot uses units where $N=1$ and the fraction of neutral dopants, $n_{0}$, is less than $2\%$ for all shown values of $\Delta$). The scaling plot (main figure) suggests that the length scale ${\ell_{P}}$, on which the screening occurs, grows almost linear in $\Delta$ in the numerically accessible regime

[TABLE]

Below we will argue that the exponent $\gamma$ is not a ‘true’ asymptotic exponent but only an effective parameter describing an extended crossover regime.

For values of $r\gtrsim 0.4\Delta^{\gamma}$, the screening charge exceeds $-Q_{0}\approx-2$. This implies that there is a substantial amount of overscreening in the system: on average too much charge of opposite sign accumulates around each charged dopant.

We have checked that other observables, for example the potential correlation function or the typical distance of neutral dopants of different type, show similar scaling behaviors, see Appendix B for an example. Most importantly, they all consistently show the importance of the length scale ${\ell_{P}}$ which governs not only screening, but also the length scale on which the dominant short-distance fluctuations of the potential occur. ${\ell_{P}}$ therefore also determines the distance of puddles of opposite charge.

V Screening by metallic surface states

Topological insulators differ from ordinary insulators or semiconductors because topology enforces the existence of conducting surface states. These states are of interest in the context of our discussion, because they provide an extra channel for screening. STM measurements of surface states can also be used to obtain quantitative information on the strength and length scale of potential fluctuation at the surface H. Beidenkopf et al. (2011). Most importantly, the suppression of puddles in thin slabs of topological insulators is expected to lead to a substantial reduction of the bulk conductivity and should therefore enhance the quality of devices based on topological insulators substantially. A major goal of this section is therefore to estimate how thin a topological insulator has to be so that puddle formation is effectively suppressed. Note that such a suppression can occur even in the absence of metallic surface states, as has been discussed heuristically by Mitin Mitin (2010) for semiconductor heterostructures.

The surface states of a $3d$ topological insulator can, generically, be described by a Dirac equation, and thus have asymptotically a density of states proportional to the doping level. Their electronic properties can be characterized by the surface doping $\mu^{S}$ and the effective fine structure constant $\alpha=\tfrac{e^{2}}{4\pi\varepsilon_{\rm{surf}}\varepsilon_{0}\hbar v_{F}}$, where, in vacuum, $\varepsilon_{\rm{surf}}=\tfrac{\varepsilon_{\rm{bulk}}+1}{2}\sim 100$. Typical values for $\alpha$ in Bi-based topological insulators are in the range of $\alpha\approx 0.1\dots 0.2$ (using, e.g., $v_{F}$ taken from ARPES data Xia Y. et al. (2009)). In Ref. Skinner and Shklovskii (2013), Skinner, Chen and Shklovskii develop a detailed analytic theory on how bulk impurity states affect the surface. We will instead investigate the question how the screening from surface states feeds back on bulk properties using some of their results.

If the surface possesses a finite doping, described by a finite chemical surface potential $\mu^{S}$, it can screen charges on a length scale described by the surface screening length $\ell^{S}_{s}\sim v_{F}/(\alpha|\mu^{S}|)$. We first consider the limit that $\ell^{S}_{s}$ is smaller than the distance of bulk impurities, $\ell^{S}_{s}\lesssim N^{-1/3}$ or, equivalently, $|\mu^{S}|\gtrsim E_{c}/\alpha^{2}$. In this case the surface state of the topological insulator acts effectively like a perfect metal. Then, screening of a dopant with charge $q_{i}$ at distance $z$ from the surface is described by positioning a mirror charge with charge $-q_{i}$ at the same distance on the opposite side of the surface. This simple screening mechanism can be implemented in a straightforward way into the model described in Sec. II. To model a thin slab of a topological insulator with two metallic surface states, one formally needs an infinite sequence of mirror charges. As described in the appendix D, an accurate and numerically efficient description is obtained by using just a single mirror charge and a linear correction term setting the potential to zero at both surfaces.

Besides its importance for applications of TI materials, the problem of an infinitely large TI slab of finite thickness $L_{z}$ has also a technical advantage which we will use in the following: the ’bare’ potential $\Phi_{0}(\bm{r})$ arising from randomly placed impurities remains finite for finite $L_{z}$ even in the thermodynamic limit (while it would diverge in the absence of surface screening). This is the potential one obtains in the absence of puddle formation when all donors (acceptors) have charge $+1$ ($-1$). One can easily calculate the distribution $p(\Phi_{0}(z))$ of this potential by averaging over the position of dopants, here $z$ is the coordinate perpendicular to the surfaces. In the following we will focus our discussion for simplicity on the distribution in the middle of the sample, $z=L_{z}/2$. Due to the central limit theorem, this initial distribution (before puddle formation) is Gaussian

[TABLE]

where the width of the distribution $\sigma(L_{z})$ can simply be computed from $\sqrt{\langle\Phi_{0}(z=L_{z}/2)^{2}\rangle}$ and, in units of $E_{c}$, is given by

[TABLE]

where the potential $V^{s}$ is defined in the appendix D and the second line is a fit to the numerical integral valid for large $L_{z}$. The prefactor of the leading term depends only on the geometry of the setup but is otherwise universal, the subleading term is linear in the chosen cutoff $a_{B}$ defined in Eq. (3). For small $L_{z}$, the width $\sigma(L_{z})$ is much smaller than the gap, implying that puddle formation cannot take place. As discussed in the introduction, Efros and Shklovskii Shklovskii and Efros (1972) estimated the length scale triggering puddle formation (for a different geometry) from the condition $\sigma(L_{z})\sim\Delta$, leading to a length scale proportional to $\Delta^{2}$ as discussed in Eq. (7). We will need in the following a more quantitative version of this argument. We will take into account, that for the successful screening of the potential in the limit $L_{z}\to\infty$, it is not necessary to redistribute $O(1)$ charges. Instead we can use the total density of neutral dopants in the thermodynamic limit, $n_{0}$, as an estimate of the volume fraction where the bare potential triggers puddle formation by becoming larger than $\Delta/2$. The characteristic width of the slab $L_{c}^{0}$ below which puddle formation is suppressed, is therefore estimated from the condition

[TABLE]

with

[TABLE]

solved by

[TABLE]

for small $n_{0}$. Using Eq. (12), we find for $\Delta\to\infty$ as an estimate for the characteristic width $L_{c}^{0}$

[TABLE]

with (sizable) relative corrections of the order of $\ln[\ln\Delta]/\log[\Delta]$. In the last line we uses that in the asymptotic regime, $n_{0}\sim 1/\Delta^{3}$,  see Eq. (7). This formula misses logarithmic corrections to $n_{0}$ which give, however only subleading terms beyond the precision of Eq.(16). While our equation can only be an order-of-magnitude estimate, we have kept multiplicative numerical prefactors to indicate their rather large numerical value.

Note that Eq. (13) and therefore also (16) was obtained only by considering properties of the bare potential, not including any self-consistent screening effects. The formulas can therefore only be viewed as a crude estimate of the relevant length scale of the problem obtained by extrapolating from the bare potential. The result clearly suggests the presence of logarithmic corrections to scaling but we cannot exclude that a resummation of logarithms leads to a modification of the logarithm, e.g., $L_{c}\sim\Delta^{2}/\ln^{\alpha}[\Delta]$. A scaling $L_{c}\sim\Delta^{\gamma}$ for $\Delta\to\infty$ with an exponent $\gamma$ smaller than $1$, can, however, be excluded as in this case the regions where the bare potential can trigger puddle formation, are exponentially suppressed. In the following we will compare the estimate from condition Eq. (13) to the full numerical solution obtained for moderately large values of $\Delta\lesssim 35$ and find that the formula nevertheless reproduces the approximately linear $\Delta$ dependence ($\gamma\approx 1.1\pm 0.1$) in the regime, $\Delta\lesssim 35$.

The surface screening will suppress potential fluctuations and the formation of puddles close to the surface. Therefore, all donors will have charge $+1$, all acceptors have charge $-1$ and the density of neutral dopants vanishes close to the surface. This physics can be captured by computing the density of neutral dopants, having a charge [math], as a function of distance from the surface

[TABLE]

where $z_{i}$ is the (dimensionless) distance of dopant $i$ from the surface and $\delta_{i,j}$ denotes the Kronecker delta. In contrast to charge and potential, the density of neutral dopants is a quantity not fluctuating in sign and therefore its average is both easier to compute (statistical fluctuations are much weaker) and to interpret.

In Fig. 3 (upper panel) we show $n_{0}(d)$ for different values of $L_{z}$. As expected, the puddle formation and therefore $n_{0}(d)$ is suppressed close to the two metallic surfaces. For sufficiently thin $L_{z}$, $n_{0}(d)$ becomes small even in the center of the sample. The lower panel of Fig. 3 therefore shows the density of neutral dopants in the center of the slab, $n_{0}(L_{z}/2)$, a function of $L_{z}$.

In Fig. 4 we show four different length scales, $\ell_{P},\ell_{c},\ell_{s}$ and $L_{c}^{0}$, as function of $\Delta$. The first three length scales have been extracted from our numerics and the equations

[TABLE]

They describe the characteristic length scale $\ell_{P}$ on which - on average - a charge is screened in the bulk (see Fig. 2), the characteristic width $\ell_{c}$ of a slab of a topological insulator below which the density of puddles drops to half the bulk value, and the length scale on which puddle formation is suppressed close to the metallic surface of a thick slab of a topological insulator (see appendix E).

All three curves show an approximately linear behavior with $\Delta$ quite different from the $\ell\sim\Delta^{2}$ (with logarithmic corrections) expected from scaling arguments. All curves are well described by fits of the form $\ell_{i}=a_{i}+c_{i}\Delta^{1.1}$. Remarkably, the same behavior is also obtained from the estimate $L_{c}^{0}$, which was obtained from the condition in Eq. (13), i.e., from properties of the bare potential (before puddle formation). Therefore $\ell_{i}$ is well described by a linear fit to $L_{c}^{0}$ as shown by the solid lines in Fig. 4. The dashed red line in Fig. 4 shows a power-law fit, $L_{c}^{0}\sim\Delta^{1.1}+const.$, which works remarkably well. We therefore conclude that (i) the average bulk screening and the surface screening are governed by the same length scale, and that (ii) one can use the ‘naive’ scaling argument, (13), to obtain this length scales (up to multiplicative factors of $O(1)$ and a small offset of $O(1)$).

As we have shown above, the asymptotic behavior for $L_{c}^{0}$ is according to Eq. (16) given by $\Delta^{2}/\ln[\Delta]$ and definitely not by an exponent close to $1$. The apparent power-law behavior with an exponent close to one therefore reflects only an extended crossover regime: there is no ‘true’ power-law with an exponent smaller than $2$ (see, e.g., Ref. Löhneysen et al. (2007) for a discussion on the definition and determination of exponents). The approximately linear behavior for $10\lesssim\Delta\lesssim 50$ arises from the interplay of logarithmic corrections at large $\Delta$ and subleading corrections for small $\Delta$, see Eq. (12).

As also the numerically determined length scales $\ell_{i}$ show the same behavior, we conclude that also in this case the numerics probes the same extended crossover regime for numerically accessible vales of $\Delta\lesssim 35$. Below, we will argue that one can use the results for $L_{c}^{0}$ to estimate the value of $\ell_{i}$ for larger values of $\Delta\sim 100$, relevant for Bi-based topological insulators.

Above we assumed a perfectly metallic surface state, $|\mu^{S}|\gtrsim E_{c}/\alpha^{2}$, which has a screening length that is short compared to the mean distance of dopants $N^{-1/3}$. Using the result given above, that the screening by bulk states sets in only at a parametrically larger scale, $\ell_{s}$, we can relax this requirement. Our results should be valid as long as the surface screening length $\ell_{S}^{S}\sim v_{F}/(\alpha|\mu^{S}|)$ is small compared to $\ell_{s}$, or $|\mu^{S}|\gg v_{F}/(\alpha\ell_{s})$.

Here $\mu^{S}$ denotes an effective chemical surface potential. Even if the chemical potential of the surface state is exactly at the Dirac point, $\langle\mu^{S}\rangle=0$, disorder will induce a finite density of states allowing for screening. Due to the charged dopants metallic puddles will form on the surface which can, in turn, screen bulk charges. To estimate the effect of these surface puddles (not to be confused with puddles in the bulk) we use the results of Ref. Skinner and Shklovskii (2013) (similar results in the context of graphene have, e.g., been obtained in Refs. Hwang et al. (2007); Adam et al. (2007)). For the computation of the resulting surface screening length $\ell_{S}^{S}$ for $\langle\mu^{S}\rangle=0$ the authors of Ref. Skinner and Shklovskii (2013) did not take into account any bulk-screening effects, which is justified as long as $\ell_{S}^{S}\ll\ell_{s}$. Under these conditions, Skinner, Chen and Shklovskii Skinner and Shklovskii (2013) found that $|\mu^{S}|\sim E_{c}/\alpha^{2/3}$ or $\ell_{S}^{S}\sim N^{-1/3}/\alpha^{4/3}$. From the condition $\ell_{S}^{S}\ll\ell_{s}$, we obtain (using our dimensionless units)

[TABLE]

where $c\approx 0.6$ according to Ref. Skinner and Shklovskii (2013), where the authors estimate $\alpha\approx 0.24$ for a Bi-based topological insulator, which results in the condition $\ell_{s}>4$ for this class of systems.

If the condition (21) is fulfilled, the surface state of a topological insulator provides sufficient screening to suppress efficiently the formation of puddles within the distance $\ell_{s}$.

VI Discussion and Outlook

We have investigated the influence of charge dopants in (topological) insulators, focusing on the case of perfect compensation with equal densities of donors and acceptors. Motivated by the physics of Bi-based topological insulators, we studied the limit where the gap $\Delta$ is large compared to the Coulomb energy $E_{c}$ of neighboring dopants with $\Delta/E_{c}\sim 100$ as a typical value Borgwardt et al. (2016). In our numerical simulations we are not able to such large values of $\Delta$. Therefore analytical estimates are needed to extrapolate to larger values of $\Delta$.

Our main focus has been the investigation of the length scales governing the formation (and destruction) of puddles. As has been pointed out before in the literature Baranovskii S. D. et al. (1984); Lee et al. (1999), due to the long-ranged nature of the Coulomb interaction and the highly non-linear screening effects, there is more than one such length scale. We have, however, found that the size of a (average) screening cloud around an impurity, the typical distance of electron- and hole puddles, and - most importantly - the length scales governing the suppression of puddle formation in the bulk due to metallic surface states, are all similar and show an approximately linear increase with $\Delta$ for $\Delta\lesssim 35$ (a fit gives $\ell_{i}\sim\Delta^{1.1}$). We have found a simple analytic estimate of such length scales based on properties of the bare potential, which reproduces this behavior in an extended crossover regime but predicts $\ell_{i}\sim\Delta^{2}/\ln[\Delta]$ for $\Delta\to\infty$.

One can use this analytic estimate to obtain a quantitative extrapolation of the numerically determined results to larger values of $\Delta$. This is shown in the lower panel of Fig. 4. The fit $\ell_{i}=a^{\prime}_{i}+c^{\prime}_{i}L_{c}^{0}$ gives an excellent fit to the numerically determined length scales $\ell_{i}$ (see figure caption for details and fit parameters). Using this extrapolation, we can estimate the corresponding length scale for large values of $\Delta$.

Assuming, for example, $\Delta/E_{c}\sim 100$ and $N\approx 10^{19}\,{\rm cm}^{-3}$, our best estimates for the dimensionless length scales are $\ell_{c}\approx 72.9\pm 15.0$, $\ell_{s}\approx 16.3\pm 2.8$ where errors have been estimated based on the use of different extrapolations of $n_{0}$, see Fig. 4. In physical units this implies that the width of the region close to the metallic surface where puddle formation is inhibited is about $62.6-88.6\,$nm. Puddle formation in the center of a thin slab of a topological insulator is predicted to be suppressed by metallic surface states by at least a factor 2 if the slab is thinner than $268.7-407.9\,$nm. As puddles largely control the bulk conduction at low temperatures by reducing the energy gap for transport processes Efros and Shklovskii (1984); Skinner et al. (2012, 2013), the suppression of puddle formation in the bulk is expected to be accompanied by a strong suppression of bulk conduction. More precisely, at least three effects will contribute to the increase of bulk resistivity the screening from metallic surfaces, the suppression of Coulomb fluctuations due to the dimensional crossover (even without metallic surfaces), and the crossover from a 3d to a 2d percolation problem of electrons moving in a correlated potential Efros and Shklovskii (1984); Mitin (2010). As we have shown that surface and bulk effects are governed by similar length scales proportional to each other, we expect that all effects are of similar importance. In compensation doped Bi-based compounds the bulk conductance is expected to be suppressed considerably (i.e. much faster than to be expected from geometric factors) when the slab becomes thinner than, e.g., $270$ nm. It will be interesting to develop a quantitative theory for transport in the future which combines numerical calculations for smaller values of $\Delta/E_{c}$ with analytic extrapolation schemes for large $\Delta$ similiar to the ones used in this paper.

Acknowledgments

We thank Y. Ando, O. Breuning, N. Borgwardt, M. Grüninger, B. Skinner, and, especially, B. Shklovskii for insightful discussions and useful comments. The numerical simulations have been performed on the CHEOPS cluster at RRZK Cologne. We thank the DFG for financial support within project C02 of CRC 1238.

Appendix A Sum rules and scaling arguments

From the definition of $C^{\textrm{nl}}_{\rho\rho}$, Eq. (8), and the Poisson equation one can derive a set of exact sum-rules

[TABLE]

Here $\langle H_{\rm C}\rangle$ is the disorder average of the Coulomb energy, $\langle H_{\rm n}\rangle$ are the single particle energies of the dopants, see Eq. (1), $\langle Q\rangle$ is the number of ionized dopants (not counting the neutral ones), and $\langle\phi^{2}\rangle$ is the expectation value of the square of the potential (all expressed in our dimensionless units).

We can use these sum-rules to obtain a more rigorous version of the scaling argument given above. We start from the assumption that the physics of the system is governed by a single length scale large compared to the average distance of impurities. In this case $C^{\textrm{nl}}_{\rho\rho}(\Delta,s)$ can be written as

[TABLE]

We will show that from this assumption alone Eq. (7) can be derived. Later, we will conclude that the scaling ansatz is not fully valid: there are substantial subleading corrections even for values of $\Delta\lesssim 35$ and logarithmic corrections in the $\Delta\to\infty$ limit. In this appendix, we will, however, explore only the consequences of the scaling ansatz.

The bulk fluctuations of the potential are of the order of the bandgap $\Delta$, see Fig. 1, which implies $\langle\phi^{2}\rangle\sim\Delta^{2}$. Furthermore, as the fraction of neutral atoms vanishes for large $\Delta$, $n_{0}\to 0$ for $\Delta/E_{c}\to\infty$, the density of charged dopants, $Q_{0}$ is of order $\Delta^{0}$ with only subleading corrections. To leading order we therefore obtain from Eqns. (24) and (25)

[TABLE]

Therefore, the scaling ansatz predicts $\beta=3\gamma$ and $2=-\beta+4\gamma$ or, equivalently, $\gamma=2$ and $\beta=6$, implying a length scale $\sim\Delta^{2}$ and a typical charge density $\sim 1/\Delta^{3}$. This is just a refined version of the argument presented above in Eq. (7).

From Eq. (22) and Eq. (27) we further deduce that the Coulomb energy density $\langle H_{\rm C}\rangle/V\sim-\Delta^{-\beta+2\gamma}=-\Delta^{-\gamma}$ ( $C^{\textrm{nl}}_{\rho\rho}$ is negative as it describes screening, for example the accumulation of negative charge around a positive one). This implies that the Coulomb energy is minimized by choosing the screening length $R_{s}\sim\Delta^{\gamma}$ as small as possible. However, this minimization competes with the increasing $H_{\rm n}\sim n_{0}N\Delta$, see Eq. (23), and of course has to respect the constraints, in particular Eq.(28).

Our numerical result are in strong disagreement with the scaling result. Several factors play a role: First, even for $\Delta\sim 35$ the asymptotic scaling regime is not yet reached. It can be seen from Fig. 4, indicating that both $\ell_{s}$ and $\ell_{P}$ are well below $10$ in this regime. Second, a more quantitative estimate based on properties of the bare potential strongly suggested the presence of logarithmic corrections, see Eq. 16. This may indicate that in Eq. (28), $\int ds\,s^{3}\,C^{\textrm{nl}}_{\rho\rho}(\Delta,s)$ obtains logarithmic corrections from a slow decay $\sim 1/s^{4}$ of $C^{\textrm{nl}}_{\rho\rho}(\Delta,s)$.

Appendix B Correlations of the potential

As show in Fig. 1 the potential in the bulk of a compensation-doped insulator fluctuates in space. It is approximately restricted to the range $[-\Delta/2,\Delta/2]$ and exceeds $\pm\Delta/2$ by an amount of order $E_{c}$ only in the region where puddles form. The correlation function $\langle\phi(r)\phi(0)\rangle$ shows on which length scale the characteristic potential fluctuations occur.

In Fig. 5 we show $\langle\phi(r)\phi(0)\rangle$ normalized to $\langle\phi(0)\phi(0)\rangle\sim\Delta^{2}$. At short distances (of the order of the distance of impurities), this is governed by the autocorrelation of the potential of a single charge and decays on a length scale set by $a_{B}$. As can be seen in the upper panel of Fig. 5, the normalized correlation function is independent of $\Delta$ in this regime. The next-largest length scale, the screening length, on which the correlations decay is more interesting. As expected, we find that this length scale is governed by the bulk screening length ${\ell_{P}}\sim\Delta^{\gamma}$, see Eq. 10. This is shown by the scaling plot in the lower panel of Fig. 5. Clearly the same length scale $\approx 0.2N^{-1/3}(\Delta/E_{c})^{1.1}$ (including prefactors) determines the screening radius and the dominant length scale of potential fluctuations. Note that scaling does not hold at the short distances ($\lesssim a_{B}$ and/or impurity distance $N^{-1/3}$) and that we had to subtract a short distance cutoff to obtain a reasonable scaling collapse.

At larger length scales, the correlation function becomes negative. This physics is, however, not governed by ${\ell_{P}}$ as follows from the absence of a scaling collapse in this regime. As discussed in the main text, the physics in the second regime is related to overscreening and occurs on a length scale which we cannot resolve with our numerical simulations.

Appendix C Density of neutral dopants

To determine the density of neutral dopants in the thermodynamic limit, we have simulated boxes of size $L^{3}$ with periodic boundary conditions. The resulting density for $\Delta=25$ is shown in the upper panel of Fig. 6 as function of $1/L$. As we do not know the analytic $1/L$ dependence of that quantity and as the numerical result is consistent with different interpolating functions, we use the following procedure. To estimate the density in the thermodynamics limit, we use a quadratic extrapolation scheme (green line), assuming that finite size effects are of order $1/L^{2}$. Within the statistical error bars this is approximately equivalent to the value obtained for the largest system size used in our numerics. A linear extrapolation in $1/L$ (black line) gives a higher value for $n_{0}(L\to\infty)$. The error bar is obtained from the largest and smallest one-standard-deviation values obtained from both interpolation schemes. It therefore reflects not only the statistical uncertainty of our data but also the much larger systematic error related to the unknown $L$ dependence of finite size effects.

The lower panel of Fig. 6 shows how the density of neutral dopants drops for increasing $\Delta$. In the crossover regime accessible to our numerics, $n_{0}(\Delta)$ decays much slower than the $1/\Delta^{3}$ law expected up to logarithmic corrections from the scaling arguments, see Eq. (7). This slow decay is, however, consistent with the slow rise of the length scales characterizing screening, which increase approximately linear instead of quadratically with $\Delta$ in the same parameter regime.

Appendix D Screened Coulomb potential in the presence of two metallic layers

The effective potential of a single impurity with (dimensionless) coordinate $\bm{x}_{i}$ with distance $x^{3}_{i}$ from a metallic layer at $x^{3}=0$ is described by (using the same conventions and cutoffs as in Eq. (3))

[TABLE]

A ”mirror charge” guarantess that the potential vanishes on the metallic surface. In the presence of two metallic surfaces located at $z=0$ and $z=L_{z}$, an infinite sequence of mirror charges is required. For efficient simulations it is mandatory that the screened potential can be computed rapidly. We have found that a potential containing only the mirror charge to the closest metallic surface in combination with a linear correction term which sets the potential to zero at both surfaces is sufficiently accurate and numerically efficient, see Fig. 7. We therefore use in our simulations for $z_{i}\leq L_{z}/2$

[TABLE]

For $z_{i}>L_{z}/2$ we use a mirror image of the potential given above with mirror plane $z=L/2$. Due to this construction the derivative of the potential with respect to $z_{i}$ has a small jump at $z_{i}=L_{z}/2$. We have found that this leads to a tiny, hardly visible bump at $z=L_{z}/2$ in the density of neutral dopants, $n_{0}(z)$. When we determining the density of neutral dopants in the center of the slab, we fit a parabola to $n_{0}(z)$ for $0.4L_{z}<z<0.6L_{z}$ omitting a tiny region of width $0.04L_{z}$ around $z=L_{z}/2$. In practice, the tiny bump (and the small correctiont to the fit described above) does, however, have no qualitative or quantitative influence on our results.

Appendix E Screening from a single surface

In this appendix we briefly discuss the suppression of puddle formation close to an metallic surface (for a system much thicker than $\ell_{c}$). The inset of Fig. 8 shows $n_{0}(d)$ for values of $\Delta$ ranging from 10 to 26. On a relatively short length scale, the bulk value of $n_{0}(d)$ is reached. In the lower panel of Fig. 8 we plot the width $\ell_{s}$ of the zone, where surface screening suppresses puddle formation, defined by $n_{0}(\ell_{s})=(n_{0})_{\rm bulk}/2$. After subtracting the offset $a_{B}=1$ we obtain numerically an approximate power law relation in the numerically accessible regime

[TABLE]

We also performed simulation with several other values of $a_{B}$ ($a_{B}=0,0.5,1.5,2$) and have checked that subtracting $a_{B}$ results in the same curve (for a fixed value of $\Delta$). Also the scaling plot in the upper panel of Fig. 8 confirms that $\ell_{d}$ governs the size of ‘dead zone’, where puddle formation is suppressed.