{'fraction_non_alphanumeric': 0.03577596131706194, 'fraction_numerical': 0.013124568270780566, 'mean_word_length': 4.54058363897305, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Electric Currents at Semiconductor Surfaces from the Perspective of Drift-Diffusion Equations

Surface sensitive electric current measurements are important experimental tools poorly corroborated by theoretical models. We show that the drift-diffusion equations offer a framework for a consistent description of such experiments. The current flow is calculated as a perturbation of an equilibrium solution depicting the space charge layer. We investigate the accumulation and inversion layers in great detail. Relying on numerical findings, we identify the proper length parameter, the relationship of which with the length of the space charge layer is not simple. If the length parameter is large enough, long-ranged modes dominate the Green's function of the current equation, leading to two-dimensional currents. In addition, we demonstrate that the surface behavior of the currents is ruled by only a few parameters. This explains the fact that simplistic conductivity models have proven effective but makes reconstructions of conductance profiles from surface currents rather questionable.

Surface sensitive electric current measurements are important experimental tools poorly corroborated by theoretical models. We show that the drift-diffusion equations offer a framework for a consistent description of such experiments. The current flow is calculated as a perturbation of an equilibrium solution depicting the space charge layer. We investigate the accumulation and inversion layers in great detail. Relying on numerical findings, we identify the proper length parameter, the relationship of which with the length of the space charge layer is not simple. If the length parameter is large enough, long-ranged modes dominate the Green's function of the current equation, leading to two-dimensional currents. In addition, we demonstrate that the surface behavior of the currents is ruled by only a few parameters. This explains the fact that simplistic conductivity models have proven effective but makes reconstructions of conductance profiles from surface currents rather questionable.

I. INTRODUCTION
Electric current measurements at semiconductor surfaces at (sub-) micron distances have become a topic of interest [1][2][3] . In such experiments, the electrostatic potential due to current injections from compact electrodes is locally probed at the surfaces. The progress in the construction of new multi-probe scanning tunneling microscopes 4,5 paves the way towards transport experiments at distances of hundreds (and perhaps tens) of nanometers. However, the understanding of the surface-sensitive current measurements is far from being satisfactory. In particular, there are several reports on the observation of two-dimensional current flows at semiconductor surfaces [6][7][8][9][10] . It is suggested that this is due to the formation of an inversion layer near the semiconductor surface. Alternatively, when two-dimensional currents are observed only at lower temperatures, decoupling from the bulk conductivity is argued, e.g., Ref. 10. A phenomenological approach widely explores the distinction between two-and threedimensional currents 11 to describe the experimental results. However, the mechanism behind the generation of the two-dimensional currents and their relationship to the three-dimensional current flows are poorly understood. In more analytical approaches, a two-step procedure is applied [1][2][3]12,13 . First, the variation of the carrier density related to the surface Fermi level pinning is calculated from the Poisson equation. Second, the current density is calculated using the (classical) drift current equation. Although this approach offers insight into the underlying physics, it imposes several difficulties. The most serious objection to the scheme results from the fact that the drift equation simply generalizes the Ohm's law. Consequently, it is not suitable for describing the p-n junctions present in inversion layers. Furthermore, the relation between the electrostatic potential behind band bending and the one behind the current flow remains elusive. An interpretation of their sum or difference has not been clarified. The last shortcoming of the scheme is that it does not take into account the recombination processes and it is unclear how they could be considered. This motivated us to model the surface sensitive transport within a more general theory, which is the system of the drift-diffusion equations 14 . These equations offer a unified description of the space charge layer and the transport experiments. In this paper, we restrict our model to two bands, i.e., electrons and holes. Extensions to other bands, e.g., surface bands, can be naturally worked out. In addition, we pin down a mathematical mechanism responsible for the appearance of the two-dimensional currents within three-dimensional models. Our analysis reveals that only a few parameters are necessary to capture the surface manifestation of the non-trivial depth-dependent conductance; two parameters are sufficient in the case of inversion or accumulation layers. This explains why simple three-layer models of the conductivity profile can be successfully used to describe the experimental results 15 . These findings are important as they shed light on the limitations of the reconstruction of the depth-dependent carrier profile from the surface-sensitive measurements 13 . The paper is organized as follows. In the next section, we show how the drift-diffusion equations describe the space-charge layer and the current transport experiments. The model involves the Poisson equation for the electrostatic potential and current equations for every considered carrier type. Here, we formulate the current equations using a quasi-Fermi level Ansatz. In sect. II B, we comment on the linearized Poisson equation and the reduction of the drift-diffusion equations to the drift equation. Next, in sect. III, we outline the presentation scheme for the surface profiles of the quasi-Fermi level and the voltage drop. The scheme facilitates the classification of solutions in terms of two-and three-dimensional currents as it takes advantage of the relevant scaling properties. Furthermore, the current equation is analyzed. After a brief summary of the known results (sect. IV A), we identify the parameters crucial for the surface currents and pinpoint the mechanism behind the appearance of two-dimensional currents (sect. IV B). These findings also enable comments on the reconstruction problem. In sect. V, we study the inversion layer at the Si(111)-Ag surface as reported in Ref. 9. Finally, we sum up our developments and point out a few open issues. Several technicalities are discussed in the Appendices The equation for quasi-Fermi levels was examined with mathematical rigor in the context of the so-called Calderon problem 16 . This theory puts the electrical impedance tomography 17 on firm ground. Despite some overlap between our work and these developments, these two problems are both physically and mathematically autonomous.

A. Model
Based on the Boltzmann transport equation for semiconductors, many calculation schemes of electric transport phenomena are derived 14,18 , among them the driftdiffusion equations. An intuitive derivation from scattering analysis can be found in Ref. 19. Notably, this theory is the first-choice scheme for device simulations 20 . The (low current) drift-diffusion model is a coupled system of equations for the electrostatic potential (multiplied by the elementary charge e) V (x) and electron n(x) and hole p(x) densities. The equations for stationary current flow read as where N (x) stands for the net density of charged ions, R for the recombination rate, f n (x) and f p (x) for external current sources of electrons and holes, respectively, ε 0 denotes the vacuum permittivity, ε r -the material dielectric constant, T -the temperature, and k B -the Boltzmann constant. The diffusion coefficient is absent in the equations due to the Einstein-Smoluchowski relation between the mobility µ and diffusion coefficient D, D = k B T µ/e. If the subscripts for the electrons and holes are not explicitly given, the relation holds for holes and electrons separately. On the right hand sides of eq. (2) and (3), other current sources, like currents generated by the electromagnetic waves, can be added. The electron j n and hole j p current densities are In equilibrium, all currents vanish: j n = j p = 0. For electrons, the condition can be rewritten in the form and hence The constant electron density n b is introduced to have a dimensionless expression under the logarithmic function. An analogous calculation can be performed for holes. The above calculation makes it evident that the equations (1)-(3) build on the Boltzmann statistics. The space-charge layer is modeled by the appropriate boundary condition imposed on V 21,22 . We denote the equilibrium potential and carrier densities with a subscript zero. These equilibrium quantities depend on the distance from the surface only and we assume that the carrier densities saturate in the bulk. Below, we will consider the convention that the surface is located at z = 0 and the semi-axis z > 0 corresponds to the crystal bulk. The surface transport experiments are believed to slightly perturb the structure of the space charge layer without any substantial damage to it. As such, we are interested in the linearization of the above equations around the V 0 , n 0 , and p 0 solutions. To this end, we write the electrostatic potential in the following form and parametrize the carrier densities with the quasi-Fermi levels ϕ, A quasi-Fermi level traces the deviations from the equilibrium occupancy of the related band. The linearized equations are obtained by plugging the above Ansatz into eq. (1)-(3) and keeping the linear terms in v, ϕ n , and ϕ p only. The final relations read as whereR stands for the linearized recombination rate, and A similar approach can be found in Ref. 22. In the following, we neglect band mixingR and this results in the decoupling of the three equations: the current equations can be calculated independently and then the Poisson equation can be solved. The boundary conditions imposed on any of the quasi-Fermi levels ϕ ensure that there are no currents flowing through the surface: We assume that the surface charge is not substantially changed by the transport experiments, and hence Upon the identification σ n = µ n n (and analogously for holes), the current equations (12)-(13) get the form of the (classical) drift current equation. Hence, previous results 12,13 can be straightforwardly adapted to the framework of the drift-diffusion equations. At this stage, one could incorporate surface bands into the theory by adding an additional band to eq. (1)- (3), but we leave this issue to a future study.
The system of drift-diffusion equations has been widely used in electronic device simulations, but one should keep in mind that it is an approximation. The model does not report on thermal effects, and so it is not suitable for systems with substantial heating of the sample. Also, the equations are based on the quasi-classical description of the electronic bands and the Boltzmann statistics for charge carriers. Hence, they may not be able to capture quantum effects and poorly perform for systems with Fermi levels close to the band edges. Models consistent with the Fermi-Dirac statistics 14 have rarely been considered in the literature. We believe that the qualitative analysis delivered below remains relevant even for the Fermi-Dirac statistical distribution of the charge carriers.

B. Linearized Poisson equation
The linearized Poisson equation has the structure of the Schrödinger operator with a source. The homogeneous part of the equation reads The potential is positive and reflects the screening efficiency of the free charge carriers. The inverse of λ(n 0 + p 0 ) is called the Debye-Hückel length. It is instructive to consider perturbations around a constant carrier density ̺. Then, the equation turns into the wellknown screened Poisson equation. A Fourier transformation of the equationP v = −λ̺ϕ yields the relatioñ Si Ge ̺ [cm −3 ] 100 K 300 K 100 K 300 K 10 13 720 nm 1244 nm 847 nm 1467 nm 10 16 23 nm 39 nm 27 nm 46 nm 10 18 2.3 nm 3.9 nm 2.7 nm 4.6 nm Table I. Debye-Hückel length (λ̺) −1/2 for silicon and germanium at various temperatures and carrier densities ̺.
where the transformed functions are denoted with tildes. If the source termφ(k) has a considerable amplitude for small k (k 2 ≪ λ̺) only, then the approximate solution reads asṽ = −φ, and hence v = −ϕ. Plugging this into the linearized current equation, we obtain where the sign of f depends on the band under consideration. It is the drift equation widely used in modeling the surface current measurements. Here, it emerges as a long-wavelength approximation to the drift-diffusion equation. It is valid far from the current source, where the quasi-Fermi levels are determined by small modes. This can also be seen from the Green's function in real space, which reads This equation shows that the potential v can be regarded as the smoothed-out quasi-Fermi level. For the slowly varying charge densities, the abovementioned Green's function can be used approximately upon the substitution ̺ → ̺(x) (adiabatic approximation). The characteristic length-scale dividing the short-and long-ranged variations is given by 1/ √ λ̺. For silicon and germanium, we give the values of the Debye-Hückle length in Tab. I. For large carrier densities, the parameter is in the range of several nanometers, and the drift equation appears as a legitimate approximation. At low carrier densities or quick variation of σ, this is not true and one needs to deal with the equation system. The intuition that the electrostatic potential smears the quasi-Fermi levels over a certain volume is well known. As such, one can expect that the notions of two-and three-dimensional currents and electrostatic potentials remain valid within the drift-diffusion equations. Indeed, if a current equation results in the quasi-Fermi level behaving at the surface as ln(r) or r −1 , then the electrostatic potential follows the same function with a different multiplicative factor (r stands for the distance from the current source). Notably, while the quasi-Fermi level associated with the logarithmic function is z-independent, the electrostatic potential will vary with the distance from the surface. More details are given in Appendix A.

III. PRESENTATION SCHEME
The distinction between two-and three-dimensional currents marks a qualitative difference. In the first case, ϕ and v at the surface vary as ln(r), in the latter as r −1 . For the sake of a qualitative discussion focused on this distinction, related scaling properties are helpful. The presentation scheme applied in recent experimental reports 6,7 makes these relations manifest. In those experiments, the current source and drain were located at the surface at (0, 0) and (D, 0), and two additional electrodes measured the voltage drop between points D 2 (1 + x, 0) and D 2 (1 − x, 0) for some 0 < x < 1. To illustrate this, we show the resistance R (voltage drop divided by the current flowing through the system) as a function of x for several D. Then 6 , the resistance R depends only on x in the case of two-dimensional currents. Both x and D are necessary to determine the resistance in the threedimensional case where the quantity R · D is independent of D. In both cases, the formulae for the resistance are analytic: for two-dimensional currents, where σ 2 is a twodimensional conductivity parameter, and for three dimensional-currents; σ 3 stands for the threedimensional conductivity. The numerical factors, σ 2 and σ 3 result from the current equation, while v 2 (0) and v 3 (0) from the Poisson equation, see Appendix A.
Obviously, if v 2 = v 3 = 1, one gets the formulae obtained from the drift equation.

A. General framework
Now, we consider the current equations (12) and (13). As already mentioned, we neglect the recombination and, in consequence, the current equations for different bands get decoupled. We concentrate on the relation where ϕ is the quasi-Fermi level, σ > 0 stands for the conductivity, and f is a source function. The boundary condition (15) prevents any current flow through the surface. A comparison of eq. (23) and eq. (12)-(13) makes it clear that the conductivity characterizes every single band. We aim for a description of the quasi-Fermi level close to the surface and the extraction of physical content. We begin with a brief résumé of the results published in Ref. 12. In that paper, it was pointed out that the theory may result in strange long-ranged behavior if d dz σ(0) = 0. We clarify this issue in Appendix B, concluding that there is no physical difference between d dz σ(0) = 0 and d dz σ(0) = 0. To arrive at the general solution of eq. (23), we look for solutions of the function ξ: which recasts eq. (23) into a Schrödinger-like form: The boundary condition (15) is The useful formula for the Green's function G is given in terms of solutions of the following one-dimensional equation where the potential U (z) is given by the formula {ψ(k; z)} k are normalized (generalized) eigenfunctions of the linear operatorL. The general solution of eq. (23) is where K 0 stands for the modified Bessel function of the second kind of the zeroth order. As demonstrated in Appendix C, operatorL has a positive spectrum only; in the case of semi-space, it corresponds to the continuous spectrum. We consider below the quasi-Fermi level due to a point source located at (x ′′ , y ′′ , z ′′ = 0) where I stands for the current supplied to the sample. This gives rise to the formula for the surface profile φ(r) = ϕ(r, z = 0): where r is the two-dimensional radius (x − x ′′ ) 2 + (y − y ′′ ) 2 . Below, we consider a onesource model, although it leads to sample electrical charging (no stationary solution). A complete model has to include both current sources and drains. Due to the linearity of the equations, this is a simple generalization of the one-source case. This paper deals with the elementary features of the model. In particular, we are interested in what can be seen on surfaces at some distance from the source. We do not examine the structure of the very contact and the details of the injection. This approach is justified by the local character of the screened Poisson equation and the fact that two-and three-dimensional currents depend on the long-ranged modes.

B. Mechanism behind the dimensional reduction
The following analysis builds on eq. (31). It clarifies that the conductivity profile impacts the surface current measurements through the function ψ 2 (k; 0). As such, we numerically investigate the possible shapes of that function for different profiles σ(z). In Fig. 1A, three close but clearly distinct conduction profiles 23 are shown, in Fig. 1B-resulting potentials U (z) and in Figs. 1C and D-functions ψ 2 (k; 0). The most important fact about ψ 2 (k; 0) is evident at the first glance, i.e., the functions have several pronounced peaks, in particular around k = 0. In general, the peaks seem to originate from resonances; shape resonances naturally appear in the potentials of the type shown in Fig. 1C with regions, from which a classical particle cannot escape unless it is given sufficient energy. It can be inferred that the resonances mark effective transmission channels between the surface and the bulk through the potential U . On a semi-axis, where the Laplace operator is not always a positive operator 24 , the zero-resonance is allowed even for a positive and monotonically decreasing potential U (z). We observe the zero-resonance if σ approaches its bulk value from above, i.e., there is some enhancement of the conductivity close to the surface. In the opposite case, there is a zero-antiresonance and ψ 2 (k; 0) has a local minimum for k = 0. We aim to describe experiments at large distances, so we concentrate on the peak for k = 0; however, analogous discussions can be conducted for any other peak. As shown in Fig. 1D, the k = 0 peaks calculated for three different functions σ are nearly indistinguishable and well approximated by a Lorentzian. Based on our experience, this is generic behavior: any well-defined peak of ψ 2 (k; 0) is satisfactorily reproduced by the Lorentz function where A stands for the peak amplitude, Γ for its width, and p for its position. Discrepancies are sometimes ob-served at the tails. Following formula (31), φ(r) emerges upon integration. Function K 0 (x) diverges at the origin and quickly vanishes starting from x ≈ 3. A thin peak in ψ 2 (k; 0) can dominate the integral so that we can approximate where q is an arbitrary parameter. This allows writing the quasi-Fermi level at the surface in the form and hence The constant term can be absorbed into the logarithm modifying q. This parameter has no physical meaning as it corresponds to shifting the potential by a constant value, and when considering both current source and drain, one can get rid of this parameter. Note that there is no such parameter in eq. (21), which corresponds to the physical measurements. So far we have taken into account the peaked structure of ψ 2 (k; 0) with no assumption as to its functional form. If ψ 2 (k; 0) is given by a Lorentzian, the integral can be performed, and the twodimensional (sheet) conductivity can be expressed as The above formula takes advantage of the fact that ψ 2 (0; 0) = σ(0)/σ(∞) 12 . If r is too large (0.1 Γ 2 r 2 for a Lorentzian), the peak becomes broad and the approximation (33) is no longer instrumental. Then, threedimensional conductivity prevails 12 . An Ansatz with a single peak for ψ 2 (k, 0) allows the modeling of the twoand three-dimensional character of the surface currents, as demonstrated in Fig. 2. Notably, parameter Γ −1 is the only characteristic length and its value is not simply related to the width of the space charge layer or any other length parameter to be identified in the system. The logarithmic solution of eq. (23) appears for a plate as a zero-mode with no z-dependence. It dominates the solution for r comparable to the plate thickness. We argue in Appendix D that the quasi-Fermi levels are approximately z-independent, even in the semi-space. It follows from the fact that the long-ranged modes closely resemble the zero-mode beneath the surface. The above considerations show how the confined currents appear in the system and allow calculations of the twodimensional conductivity from the function σ(z). The physical significance of the structure of the theory goes beyond this statement. Assuming the Lorentz shape of the peak around k = 0, we see that the surface conductivity far from the source is governed by two parameters (Γ, σ(∞)) only. At smaller distances, additional three parameters can appear: the position, width, and amplitude of the second peak. Typically, the second peak becomes important at distances well below 1 µm and the third one at distances where the theory is no longer valid. As a consequence, the reconstruction of the conductivity profile σ(z) can be doubted since many functions σ will lead to the same peak structure.
To demonstrate how the theory works, we comment on the current flow measurements on Si(111) √ 3 × √ 3-Ag as reported in Ref. 9. At room temperature, both ptype and n-type doped samples were investigated. For ptype doped samples, three-dimensional currents were observed, while for n-type doped samples, two-dimensional currents were seen. This was explained as a result of the surface Fermi level pinning 0.16 eV above the top of the valence band and the formation of inversion and accumulation layers for n-type and p-type doped samples, respectively. As such, C. Liu et al. assumed that the two-dimensional current was flowing through the surface channel. The system has been experimentally investigated in many ways in the context of the surface conductivity 2,25 . There is no doubt that the current through the surface states depends also on the bulk electronic structure, as evidenced by the fact that the two-dimensional conductance coincides with the appearance of the inversion layer. As a consequence, if the surface band is an active transport channel, it is not known what fraction of the current goes through. We have no good model of the phenomenon, so we assume a limiting case-the inversion layer as the dominant channel. Reports on Ge(001) and Ge(001):H 7 demonstrate that it is always a viable option. The bulk Fermi level in the sample is said 9 to be about 0.25-0.30 eV below the conduction band at room temperature. This is consistent with a low donor density, in calculations we assume 6 · 10 15 cm −3 . Four point probes were arranged in a line with an equal spacing of 20 µm.
In that geometry, the change to the resistance in the temperature range 120-300 K was measured. The resistance decreased along with the temperature from 300 K to about 140 K; however, at lower temperatures it remained roughly constant. The measurements were performed on samples of 500 µm in thickness. We begin by calculating the equilibrium quantities: V 0 , n 0 , and p 0 ; see sect. II A. The equilibrium densities of the carriers at 120 and 300 K are shown in Fig. 3. The hole density at the surface highly exceeds the electron one, so we neglect the electron current. Notably, the carrier density at the surface is not governed by the Debye-Hückel length, which shows how the electrostatic potential changes. Next, we calculate the function ψ 2 (k; 0). The peak for k ∼ 0 appears to be extremely thin. Due to the accumulation of numerical errors, we cannot accurately analyze its structure. The data for k > 10 −6 µm −1 suggest a single peak with Γ ∼ 10 −10 µm −1 ; more sophisticated numerical work could probably alter this value but cannot make Γ −1 commensurate with the sample thickness d. This makes it necessary to reanalyze the model in terms of discrete eigenvalues and eigenfunctions; see Appendix E. The result can be easily stated. What counts is the zero mode only. The contribution of the mode associated with the lowest non-zero eigenvalue k 2 1 is proportional to ψ 2 (k 1 , 0). Since k 2 1 ∼ d −2 , ψ 2 (k 1 ; 0) is orders of magnitude smaller than the zero mode. The electrostatic potential reads as v(r, z = 0) = Iv 2 (0) where q is a dummy parameter, as discussed in sect. IV B. The zero mode usually prevails at a distance from the source comparable to the sample thickness. Here, it is the conductivity profile that promotes that mode to the only measurable mode. Note that it explains why twodimensional currents are detected at distances from several micrometers to half a millimeter 9 with no substantial change to the resistivity. Following eq. (13) the current density j p reads µ p p 0 ∇φ. Under the assumption that the mobility is a single num- ber, i.e. it does not vary with z, the current density in the radial direction is proportional to p 0 (z)/r and vanishes in the other directions. So, the largest current density is close to the surface and it rapidly decreases with the distance from the surface; see Fig. 3. In our approach, the conductivity of the inversion layer is decoupled from the bulk conductivity. It is in stark contrast to the models using the drift equation that assume a continuous transformation of the conductivity between the hole conductivity at the surface and the electron one in the bulk 13 . The use of the drift-diffusion equations also has quantitative consequences. This is evident for the mobility values that we can calculate from the resistance R measured in the experiment 9 and our model. The relevant formula reads as where v 2 is defined by eq. (A3). The results are shown in Tab. II. The obtained mobility at room temperature is extremely small, but surprisingly, the values on the order of 10 cm 2 /(V·s) were reported in Ref. 26 for the surface channel on the Si(111)-√ 3 × √ 3-Ag. The variation of the mobility with temperature is reasonable within the drift-diffusion equations; the change calculated from the drift equation seems to be nonphysical. The mobility in inversion layers has been investigated for a long time 20 , mostly in the context of microelectronics 27 . A further study is necessary to determine if the mobility calculated here can be simply related to those results.

VI. CONCLUSIONS
In this paper, we described the surface-sensitive transport measurements by the drift-diffusion equations. The outlined calculation scheme offers a unified and logical framework. The space charge layer corresponds to the equilibrium solution, while the transport equations are regarded as a small perturbation of those solutions. We qualitatively characterize the screened Poisson equation needed to calculate the electrostatic potential. This allows us to draw a connection between our model and the models using the drift equation. We also analyze the mechanism leading to two-dimensional currents. It reveals that the surface currents are governed by very few parameters for the accumulation and inversion layers. It is unclear if the Poisson equation can help in reconstructing the conductivity profiles. Finally, we analyze the experimental data for the Si(111)-√ 3 × √ 3-Ag within a simple model of the inversion layer. In this case, the observed two-dimensional currents are due to the zero mode of the minority charge carriers. The work can be continued in many ways. To develop a complete model comparable to experiments, the incorporation of the surface bands is imperative. Another interesting issue is the universality of the Lorentzian shape of the function ψ 2 (k; 0). Both the mathematical mechanism resulting in this shape and an explicit formula for Γ are highly anticipated. The limits for the reconstruction of the conductivity profile, which are approximately outlined in our work, need further elaboration. Finally, the implication of the close-to-surface mobility (as calculated above) and its transferability between various systems is an open issue that can be highlighted by future experiments.  p b +n b at 300 K and 120 K. The bulk carrier density reads 6 · 10 15 and 5 · 10 15 cm −3 at 300 and 120 K, respectively. Lower panel: The number of holes located between the surface and the depth z as a fraction of all the positive current carriers in the sample. The exact value of the sample thickness d in the range 100-1000 µm has no impact on the graphs. These graphs demonstrate that half of the holes are located in the 5 nm zone beneath the surface, while nearly all positive charge carriers are located in the 20 nm subsurface zone at room temperature. At 120 K the holes are concentrated even closer to the surface.

ACKNOWLEDGMENTS
The author thanks Professor J. Konior for reading the manuscript, and anonymous referees for insightful discussions and useful references. Funding for this research has been provided by the EC under the Large-scale Integrating Project in FET Proactive of the 7th FP entitled "Planar Atomic and Molecular Scale devices" (PAMS). allows writing the Poisson equation in the form of eq. (A3). However, in this case it is an approximate equation, neglecting terms of the order z/(r 2 + z 2 ) 3/2 , and hence valid close to the surface only. The local nature of the screened Poisson equation suggests that the solution at the surfaces is close to the all-inclusive solution.

Appendix C: Absence of bound eigenstates ofL
In this section, we argue that the current operatorL, see eq. (25), admits no negative eigenvalues. Solutions of eq. (23) with no source correspond to extrema of the functional which is interpreted as the energy dissipated by the current in the source-free region of space per unit time 28 when ϕ stands for the voltage drop (within the drift equation). The integration is done over the whole region. For any given function ϕ, the functional results in a non-negative number. Upon substitution (24) the energy gains an equivalent form excluding any negative eigenvalue as well as the zero eigenvalue with a square-integrable eigenfunction. The equivalence between functional (C1) and (C2) can be shown if boundary condition (26) is taken into account. In the cases of interest here, there are two options. For the semi-space approximation R 2 × [0, ∞), the potential U (z) defined by eq. (28) asymptotically approaches zero, U (z) → 0 for z → ∞. Hence, no bound states are allowed and only a continuous (scattering) spectrum is present. Thus, the conclusions of Ref. 12 are valid in any case; the current will have the three-dimensional character at sufficiently small and large distances from the source. Therefore, the confinement close to the surface may only appear at limited distances. The case of a real sample of a finite thickness (R 2 × [0, d]) is described in Appendix E.
Numerical investigations show that for a class of physical σ(z) this approximation is valid for z 0 in the range of the surface layer. The difference between u(k; z) and the (generalized) eigenfunctions ψ(k; z) lies in normalization. The multiplicative factor can be easily found u(k; z) = σ(0) σ(∞) ψ(k; z) ψ(k; 0) .

Appendix E: Finite-thickness effects
For the current equation (23) considered for samples of finite thickness d, the expansion of the Green's function in terms of the eigenfunctions and eigenvalues reads as K 0 (k n r)ψ(k n ; z)ψ(k n ; z ′ ), (E1) where k 2 n stands for the n-th eigenvalue of the operator L; see eq. (25). In Sect. V, we consider the system that strongly enhances the modes close to k = 0. The width Γ of the peak in ψ 2 (k; 0) is orders of magnitude smaller than k 1 ∼ d −1 . Hence, the above modes cannot describe the surface currents. However, to any solution of eq. (23) zero modes can be added. Far from the source, the zero mode has the form where C is an arbitrary constant and q is the dummy parameter. Note that q and C are not independent. Heuristically, the logarithm in the above formula is due to sum-ming up all modes parallel to the surface with the zero eigenvalue in the normal direction, i.e., As such, we expect that only the zero mode ∼ ln qr + C contributes to the measured current. This approach supports the conclusions of Appendix D that quasi-Fermi levels associated with the logartithmic currents are zindependent. Furthermore, it allows for a non-rigorous estimation of parameter Γ by comparing the above coefficient of ln(·) function and eq. (36). It results in the relation which may be used if the resulting Γ does not depend on d over several orders of magnitude.