Nonlinear screening and stopping power in two-dimensional electron gases

We have used density functional theory to study the nonlinear screening properties of a two-dimensional (2D) electron gas. In particular, we consider the screening of an external static point charge of magnitude Z as a function of the distance of the charge from the plane of the gas. The self-consistent screening potentials are then used to determine the 2D stopping power in the low velocity limit based on the momentum transfer cross-section. Calculations as a function of Z establish the limits of validity of linear and quadratic response theory calculations, and show that nonlinear screening theory already provides significant corrections in the case of protons. In contrast to the 3D situation, we find that the nonlinearly screened potential supports a bound state even in the high density limit. This behaviour is elucidated with the derivation of a high density screening theorem which proves that the screening charge can be calculated perturbatively in the high density limit for arbitrary dimensions. However, the theorem has particularly interesting implications in 2D where, contrary to expectations, we find that perturbation theory remains valid even when the perturbing potential supports bound states.


I. INTRODUCTION
The screening response is a fundamental property of an electron gas in arbitrary dimensions. The situation in two dimensions is of particular interest because of the possible realization of quasi-two-dimensional systems in semiconductor heterostructures 1 , image and band-gap surface states at metal surfaces 2 , electrons at the surface of liquid helium 3 , metallic overlayers on insulating substrates and layered materials 4 . Examples of problems in which an understanding of the screening of an external charge is important include impurity-limited electron transport in two-dimensional electron gases 1 , the electronic structure of intrinsic or photo-induced defects in semiconductor heterostructures 5 and dynamic interactions with moving charges 6,7 . In each of these cases, the strong interaction of mobile electrons with the charged impurity leads to significant modifications of the local electronic structure.
One of the first attempts to deal with 2D screening was the work of Stern and Howard 8 which was motivated by measurements of electron mobilities in Si inversion layers. To perform explicit calculations they considered a model in which the two-dimensional electron gas (2DEG) is represented as a sheet of zero thickness, the ideal 2D limit, with a charged impurity situated a distance d from the sheet. The screened impurity potential was then obtained within a Thomas-Fermi approximation and was used to investigate the possibility of impurity bound states. However no attempt was made to perform fully self-consistent calculations of the screening charge density and the screened potential. This was taken up later in a series of papers by Vinter 9 for a Si inversion layer. This work is especially notable for its attempt to realistically account for the underlying bandstructure of the semiconductor and the finite thickness of the 2DEG. A detailed discussion of screening and bound state formation was presented for this particular heterostructure and the results were shown to be consistent with observed impurity-limited mobilities.
Also of interest is the dependence of screening on the dimensionality of the system. That important differences might arise is indicated by the well known fact that any purely attractive potential in 2D always has at least one bound state 10 , in sharp contrast with the situation in 3D. Extending previous nonlinear screening calculations from 3D 11 to 2D and exploring these possible differences is part of the motivation for the work described in this paper. To do this we perform self-consistent nonlinear screening calculations for an ideal 2DEG within the context of density functional theory. Our studies as a function of density reveal some interesting properties of the screening in 2D that had not been appreciated before. The interpretation of these numerical results is facilitated by the proof of a general theorem regarding the nature of screening in the high density limit. A preliminary account of some of this work was presented earlier 12 .
As a specific application of the results we consider the low velocity stopping power for a heavy projectile moving parallel to the plane of the 2DEG. This problem has been treated previously in both linear and nonlinear response approximations. Of the linear theories, we mention the early work of Fetter 13 using a hydrodynamic model which was followed by a calculation based on the random phase approximation (RPA) by Horing et al. 14 . In this latter work a finite distance, d, between the projectile and the plane of the gas was considered. The case of a projectile moving in the plane of the gas has also been studied, not only within the RPA 15 , but also with the inclusion of local field corrections 16 and for a gas at finite temperatures 17 . Of course the interaction cannot always be considered as weak and a fully nonlinear theory of the screening is in general necessary. One step in this direction is the RPA quadratic response formulation 18 which was used by Bergara et al. 19 in 2D to provide the lowest order nonlinear correction to the stopping power. Although this theory is perturbative and therefore limited in its range of validity, it has the merit that the nonlinear correction is determined for arbitrary projectile velocities.
The stopping power at low projectile velocities can also be formulated in terms of the scattering of electrons from the screened potential. For qualitative purposes it is sometimes useful to consider model potentials 20 whose arbitrariness can be limited by imposing physical constraints such as the Friedel sum rule. This was the approach taken by Wang and Ma 21 to determine proton and antiproton stopping powers in a 2DEG. A somewhat more fundamental approach based on the nonlinear screening theory of Sjölander and Stott 22 was used by Krakovsky and Percus 23 for negatively charged projectiles. Their results for this case are consistent with the model potential results of Wang and Ma 21 . None of these calculations, however, are truly self-consistent. Our results based on self-consistent density functional theory calculations eliminate this source of uncertainty in the determination of 2D stopping powers.

II. SELF-CONSISTENT SCREENING
The electrons of our ideal 2DEG are assumed to have an isotropic effective mass m * and move in the presence of a uniform neutralizing positive background. We also imagine the entire 2D system to be immersed in an extended dielectric with permittivity ε as might occur in a semiconductor heterostructure. Other dielectric arrangements can also be envisaged, for example the surface of liquid helium, but these will not be considered here. All lengths will be expressed in units of the effective Bohr radius a 0 = εh 2 /m * e 2 , and all energies in units of the effective Hartree, H = e 2 /εa 0 . The mean density of the gas, n 0 , is characterized by the density parameter r s = 1/ √ πn 0 . Our objective is to determine the nonlinear screening of a stationary point charge, Z, located a distance d from the plane of the 2D gas. The methodology follows that used in earlier calculations for a three-dimensional gas 11,24 and is suitable for both positive and negative charges. This screening response is determined by solving self-consistently the two-dimensional Kohn-Sham equations where the effective potential is given by The external potential v ext (r) in the plane of the gas is −Z/ √ d 2 + r 2 , and ∆v H (r) is the Hartree potential due to the electronic screening density, ∆n(r) = n(r) − n 0 . The change in the exchange-correlation potential, ∆v xc (r) = v xc [n(r)] − v xc [n 0 ], is defined in the local density approximation (LDA) using the parameterization of the 2D exchange-correlation energy given by Tanatar and Ceperley 25 . The total screening density is given by where the first sum extends over all bound states of the effective potential, and the second extends over all occupied continuum states up to the Fermi level E F . We assume that each spatial orbital is doubly-occupied for spin. The free-particle solutions ψ 0 i (r) are obtained in the absense of the external potential. Due to the cylindrical symmetry of the problem, the solutions of Eq. (1) have the form ψ i (r) ≡ e imφ R nm (r) where the angular momentum quantum number m takes on integral values and n distinguishes the different radial solutions 8,26 . The continuum states behave asymptotically as R km (r) ∼ r −1/2 cos(kr − 1 2 |m|π − 1 4 π + η m ) where k = √ 2E is the wavevector and η m (k) is the 2D scattering phase shift. For the free-particle solutions, η m = 0. These scattering phase shifts are related to the total screening density according to the 2D Friedel sum rule (FSR) 8 At self-consistency, this quantity must equal the charge Z of the external impurity. The solutions of Eq. (1) are obtained using standard numerical techniques. The calculation of ∆v H (r), the only nontrivial step in determining the effective potential, is performed by making use of a fast fourier transform technique to go successively between real and wavevector space. Details of the method are given in an Appendix.

III. NUMERICAL RESULTS
Self-consistent calculations were performed for a range of densities and impurity charges Z, including negative values. In addition, calculations were done with and without exchange-correlation effects included in order to compare with earlier linear calculations which were performed within the random phase approximation (RPA).
To begin we show results for a charge Z = 1 in the plane of the 2DEG (d = 0). Fig. 1 presents the self-consistent effective potential ∆v eff for r s values ranging between 1 and 10 in steps of 1. This range spans low to high densities and in particular includes densities of experimental interest. Also shown for comparison is the Thomas-Fermi (TF) screened potential which, as will be explained shortly, is the r s → 0 limit of the DFT calculations. As such, we see that the potential for r s = 1 is already quite close to the TF limit.
The variation with r s is systematic, evolving with increasing r s into a potential with a strongly repulsive region beyond r ≃ 0.5 a.u. This behaviour can be explained in terms of the underlying electronic structure of the screened impurity. For all the densities shown, the screened potential supports one bound state which is occupied by two electrons. Since the total screening density must integrate to unity according to the FSR in Eq. (5), the continuum states must themselves contribute a total charge of +1 in order to compensate for the overscreening provided by the two bound electrons. This is illustrated in Fig. 2 which shows the bound and continuum densities for two representative r s values, r s = 2 and 10. It can be seen that the bound state distribution is very similar despite the very different background densities. This is consistent with the similarity of the inner region of the potentials in Fig. 1 which is primarily responsible for the shape of the bound state. At large distances from the impurity the densities in Fig. 2 exhibit characteristic Friedel oscillations with wavelength π/k F . The negative portion of the screening densities in Fig. 2 corresponds to a local charge density which is positive and provides the necessary compensation of the negative bound state electronic charge. However, the distribution in space of this compensating charge is quite different in the two cases and leads to the different behaviour of the r s = 2 and r s = 10 potentials for r > 0.5 a.u. The form of the r s = 10 continuum density is particularly interesting. In the range 0 ≤ r ≃ 10 a.u., ∆n cont (r) has approximately a constant value of −n 0 (hence r∆n cont (r) behaves linearly as seen in Fig. 2). In other words, the Z = 1 impurity with 2 bound electrons can be viewed as an H − ion which sits at the center of a positively charged disc of radius R ≃ r s . Returning to Fig. 1, we can now understand the repulsive part of the potentials seen for large r s . The LDA potential of a 2D H − ion itself has a +1/r Coulomb tail arising from the Hartree potential. In the 2DEG environment, this H − potential is screened by the positive background and as a result, the potential goes to zero at approximately R ≃ r s .
Another interesting quantity is the bound state energy shown in Fig. 3 as a function of r s . The solid curve is the result of the full DFT calculation including xc, while the dashed curve is the corresponding result obtained when ∆v xc is set to zero. We refer to the latter as the Hartree approximation. The Hartree result is monotonically decreasing with decreasing r s and goes to a limiting value of E 0 = −0.2862 H at r s = 0 which, as we shall see, is the binding energy for the TF screened potential. The xc result likewise decreases with decreasing r s , but reaches a mininum near r s ≃ 0.5 and then increases to the same limiting value at r s = 0 as in the Hartree approximation. The more negative xc eigenvalue reflects the stronger binding due to the xc potential which is attractive in the vicinity of the positive impurity as a result of the pile-up of the electron screening density. To illustrate this we compare the self-consistent potentials with and without xc for r s = 2 in Fig. 4. It can be seen from this figure that ∆v eff with xc is more attractive in the core region which is important for the determination of the bound state. Associated with this is a relative phase shift of the long range Friedel oscillations between the two cases.
The decreasing trend in the xc eigenvalue is consistent with the behaviour of the potentials in Fig. 1 but is opposite to what is found in the analogous 3D calculations 11 . There the bound state energy increases with decreasing r s , that is, the bound state becomes shallower with increasing density. The explanation for the 3D behaviour is that the screening of the impurity potential is more effective with increasing density, and as a result, the bound state eventually ceases to exist beyond a certain critical density. This behaviour is consistent with the expected applicability of perturbation theory in the high density limit when the Fermi energy E F is much larger than the magnitude of the screened potential |∆v eff (r)| 27 . If one were to use the same reasoning in 2D, the screening density in wavevector space would be given by the linear response result (Sec. IV provides further details) where χ 0 (q) is the Fourier transform of the static noninteracting density response function of the 2D gas. Within linear response theory, the screening density arises from the perturbation of the plane wave states of the uniform gas and there cannot be a bound state contribution. However, as we have just seen, bound states do persist in 2D even in the high density limit. This is a peculiarity of 2D and is associated with the fact that a purely attractive potential always has at least one bound state regardless of the strength of the interaction 10 . This effect is obviously missed by linear response theory and would seem to invalidate a perturbative approach. Nevertheless, it turns out that the results of linear response theory are indeed correct even though the theory neglects bound states entirely. We defer discussion of this point to the following Section. The opposite limit r s → ∞ is also of interest. In the Hartree approximation, the bound state eigenvalue approaches zero near r s = 8. Whether a very shallow bound state persists beyond this value could not be confirmed due to difficulties in obtaining converged self-consistent solutions. This was less of a problem with the inclusion of exchange and correlation since these additional interactions stabilize the bound state relative to the Hartree calculation. In the range 8 ≤ r s ≤ 12 the bound state eigenvalue could be fit approximately to the expression E 0 ≃ 0.023 − 0.7/r s H. Extrapolation of the numerical data beyond r s = 12 suggests that the bound state eigenvalue goes to zero at some large, but finite, r s value. However, it was not possible to confirm this since it was difficult to obtain converged solutions for large values of r s . Thus we cannot state unequivocally what the large r s behaviour is, and in particular, whether or not a free H − ion is stable in 2D within the LDA. In 3D it is known 28 that the H − ion is not stable within the LDA and our data suggest that this might also be the case in 2D. This remains an interesting question to explore in the future.
As a final comment about the eigenvalue obtained with xc, we note that the H − ion is stabilized in the range studied by the fact that it sits at the center of a positively charged disc of radius R ≃ r s . The electrostatic potential due to this charge has the value −2/r s H at the centre of the disc relative to infinity. The depletion of charge also gives rise to a shift in the xc potential of −v xc (n 0 ) which is positive and therefore a destabilizing effect. For large r s , −v xc (n 0 ) ≃ +1.6/r s H which, when combined with the electrostatic shift, gives a potential shift of about −0.4/r s H. It is clear that this represents a significant contribution to the xc eigenvalue at the larger r s values we have considered. However the existance of a bound state cannot be determined without a detailed knowledge of the self-consistent potential near the impurity. In the case of the Hartree approximation, the stabilizing effect of the electrostatic potential of the positively charged disc (−2/r s H) is clearly insufficient to compensate for the destabilizing effect of the Hartree interaction between the two bound electrons in the large r s limit.
We consider next the situation of a negatively charged impurity such as an antiproton or acceptor state in a semiconductor 29 . Fig. 5 illustrates the variation of the screening density with r s . The results found here in 2D are similar to those found previously in 3D 24 . The negative impurity repels electrons from its vicinity and roughly speaking exposes a positively charged disc of radius R ≃ r s which neutralizes the impurity. In other words, the impurity sits at the center of a hole in the electron gas. The similarity to the situation for a positive impurity is emphasized in Fig. 6 where we superpose the screening densities for Z = 1 and Z = −1 at r s = 10. The Z = −1 screening density is in fact very similar to the continuum part of the Z = 1 screening cloud. Thus, at sufficiently large r s , the H − ion for Z = 1 effectively acts as an external Z = −1 impurity as far as the rest of the electron gas is concerned. This emphasizes that the H − configuration exists as a well-defined entity in this range of r s values.
In all cases we find that the screened Z = −1 potential does not support any bound states. We mention this since it has previously been claimed 30 that the introduction of a negative test charge into a 2D gas could give rise to potentials which would bind an electron (or rather, a second negative test charge). To make contact with this earlier work we compare in Fig. 7 the nonlinearly screened potentials with those obtained on the basis of linear response theory. One curve shows the screened potential as obtained in the RPA. This potential is defined in Fourier space as where the RPA dielectric function is given by 31 At this level of approximation, only Hartree interactions are relevant and the screened potential defined above is the potential experienced by a Z = −1 test charge moving in the vicinity of the impurity. A second curve shows the linearly screened potential, v LF C sc = v(q)/ǫ LF C , obtained when local field corrections (LFC) are included in the definition of the dielectric function. Within the LDA, this is given by where v ′ xc (n 0 ) = dv xc (n 0 )/dn 0 is the LDA local field correction, i.e., in standard notation, v(q)G(q) = −v ′ xc (n 0 ). The potential v LF C sc is simply the Hartree potential obtained with a screening charge density that includes the effects of LFC. As can be seen in Fig. 7 for r s = 4, this potential has a large negative region in real space and supports a bound state for a unit negative test charge of one electron mass. The authors of Ref. 30 then argue that the potential also binds an electron in the 2DEG and that this binding may be a relevant pairing mechanism. However this conclusion is invalid for several reasons. First, the Hartree potential for the nonlinearly screened impurity (the short-dashed curve in Fig. 7) has a much shallower attractive portion as compared to the corresponding result of the linear calculation (long-dashed curve). A similar result was found analytically for a point charged screened by a uniformly charged disk 32 . But more importantly, an electron, as opposed to a negative test charge, also feels the effect of the induced xc potential. With this contribution included, the full self-consistent nonlinear potential ∆v eff defined in Eq. (2) has a very shallow attractive part, and there is no tendency for bound state formation, as confirmed numerically. At a more fundamental level, the question of pairing involves the effective interaction between pairs of electrons at the Fermi level which typically have a high relative velocity. From this point of view it is unclear to what extent a statically screened potential for a test charge can be used as an estimate of the pairing interaction.
We next consider the effect of moving the impurity charge out of the plane of the 2DEG. This is the situation corresponding to an external charge incident on the 2DEG from the outside, or a remote ionized impurity in a semiconductor heterostructure. 1,5 Fig. 8 shows the screening density for a few values of d. It can be seen that the amplitude of the induced density decreases rapidly with increasing d and the density becomes more spread out along the plane. The reason for this behaviour is that the potential loses its Coulomb singularity as soon as the impurity is out of the plane, and therefore becomes smoother and weaker. In Fig. 9 we show the nonlinear screening density for the case r s = 2 and d = 2 a.u. together with the bound state and continuum contributions. The bound state in this case is very shallow and extended, in fact more extended than the total density itself. This is possible since the continuum contribution to the density has a similar extent but is of opposite sign, so that there is a cancellation of the densities in the asymptotic region.
Also shown in Fig. 9 is the screening density as determined by linear response theory. Despite the fact that the nonlinear screening response has a bound state contribution, the linear response result is seen to agree very well. The reason for this agreement is not obvious and will be explained in detail in Sec. IV. In addition, the density is already fairly close to the classical image theory result (d/2π)(r 2 + d 2 ) −3/2 . This expression for d = 2a 0 gives a density at the origin of 0.04 a.u. as compared to the nonlinear result of 0.035 a.u. With further increases in d the nonlinear density rapidly approaches the classical screening charge density.

D. Stopping Power
As discussed in the Introduction, there have been a number of calculations of 2D stopping powers. Our purpose here is to present results which are based on the full nonlinear screening charge densities and potentials. Within the so-called kinetic theory framework, the stopping power at low velocities for a projectile moving with velocity v in the plane of the 2DEG is given by the expression 33 where σ tr (E F ) is the momentum-transfer cross section defined in terms of the scattering phase shifts by 8 To leading order in the velocity, it is sufficient to determine the scattering phase shifts using the static nonlinearly screened potentials calculated in the present paper.
In Fig. 10 we show the stopping power as a function of the projectile charge Z for d = 0 and r s = 2. For small Z, S has the expansion S = S 1 Z 2 + S 2 Z 3 + . . . where the first two terms are the linear and quadratic response results, respectively. To emphasize the deviations from linear response, we present the results in the form S/(vZ 2 ). In this representation, the stopping power including the quadratic response correction appears as a straight line with slope S 2 /v. This correction was previously calculated within quadratic RPA 19 and is shown in Fig. 10 as the straight line. The fact that the line is tangent to the nonlinear Hartree curve at Z = 0 shows that the dynamic quadratic response formulation is consistent with the kinetic theory approach for small values of Z and for low velocities. However, the validity of the quadratic response result is mainly limited to negative charges down to Z ≃ −1 and breaks down completely for Z ≃ 1. It is thus apparent that nonlinearities are generally very important. Furthermore, it can be seen that the inclusion of xc leads to a large enhancement of the stopping power in the range −1 ≤ Z ≤ 1, even to a greater extent than found in 3D. 34 This emphasizes that properties associated with the scattering of continuum states will be strongly affected by the xc part of the self-consistent potential.
The behaviour of the stopping power as a function of density is also of interest and in Fig. 11 we show r s S/v vs. r s for Z = 1. The quantity S/v has the physical interpretation of a friction coefficient 33 and is actually a monotonically decreasing function of r s with a finite limiting value of πZ 2 for r s → 0. This result follows 20 from the 2D momentumtransfer cross section σ tr (E F ) = Z(2π/v 2 F ) tanh(πZ/v F ) which is obtained when scattering from a bare Coulomb potential, −Z/r, is considered. This is the exact result in the high density limit and simply reflects the fact that screening is not important for high energy scattering. As a result of this finite limiting value, the quantity r s S plotted in Fig. 11 shows a maximum at r s ≃ 1 before going to zero for r s → 0.
The results obtained previously by Wang and Ma 21 for Z = 1 are qualitatively similar to our results although there are some quantitative differences. We recall that they approximated the screened potential by the linearized TF potential and adjusted a screening parameter in order to satisfy the FSR. This potential is everywhere attractive and misses the repulsive region shown in Fig. 1 which is associated with H − formation at low densities. These differences become more apparent when the stopping powers for Z = 1 and Z = −1 are compared. This is shown in the inset to Fig. 11 where the ratio of these stopping powers is given. The ratio tends to unity for r s → 0 since linear response is valid in this limit. However, the ratio also tends to unity for large r s since as we have seen, the H − configuration is effectively the same as an external negative charge insofar as the scattering of continuum states is concerned (see Fig. 6). In contrast, for an assumed model potential, Wang and Ma find the ratio of stopping powers to be close to 2 for all r s . This points to the danger of using model potentials which do not properly account for the true nature of the electronic screening. Another difficulty of the model potential approach is that differences associated with different screening approximations are no longer apparent. This is emphasized in Fig. 11 by the plot of the stopping power in the Hartree approximation. The results are qualitatively similar to those with xc, but there are important differences. Clearly additional input is needed in a model potential approach to capture these differences.
Finally, we compare with the results of Krakovsky and Percus 23 for Z = −1. As mentioned in the Introduction, they obtain the screening charge density using the Sjölander-Stott integral equation 22 which is a different way of introducing exhange and correlation. This induced density is shown in a second paper 35 and is qualitatively similar to our results in Fig. 5. They then use this density to construct an effective scattering potential for which the scattering phase shifts are determined. However they give no details about how their potential is constructed nor what interaction effects, such as exchange and correlation, are included. Nevertheless, the stopping powers they obtain for Z = −1 in 2D are similar to our results. We find in fact that our Hartree and xc stopping powers bracket their results with deviations which are typically less than 20%.

IV. HIGH DENSITY LIMIT
We briefly commented earlier on the behaviour of the screening cloud in the high density limit (E F → ∞, r s → 0). In this Section we investigate this question in more detail and in particular, establish the connection between linear and nonlinear screening for a 2DEG. This problem in 3D was touched on in an early paper by Butler 36 which dealt with the possibility of positronium formation in metals. Using a square well model potential, it was shown that the total electron density at the origin is an analytic function of the strength of the potential. This in fact is a special case of a more general result obtained by Kohn and Majumdar 37 concerning the analyticity of physical properties of the entire electron gas. Butler then used his result regarding the analyticity of the density to argue that the total density at the position of a positron in a metal could be obtained by perturbation theory (to all orders) even if the screened positron potential had a bound state. This assertion was confirmed for the specific example of the square well model where lowest order perturbation theory already provided a good estimate of the exact total density in the limit of high electron gas densities.
We shall now prove a theorem which demonstrates that some of Butler's conclusions have a much broader range of validity. It is unnecessary to specify either the form of the externally imposed potential nor the dimensionality of the system. The physical situation we consider is a noninteracting Fermi gas in the presence of an external potential λV (r). The effect of electron-electron interactions will be dealt with later. The potential is assumed to be smooth and to approach zero sufficiently rapidly as |r| → ∞. The parameter λ is an arbitrary coupling constant whose physical value is unity.
The problem is to determine the induced density ∆n(r) which is due to the introduction of λV (r). This quantity is conveniently defined in terms of the Green's operator where H 0 is the unperturbed hamiltonian. In a spatial representation, the Green's function is given by where ψ i (r) is an energy eigenstate of H. From this we see that Thus the total density is given by where the upper limit of the energy integration is the Fermi energy E F below which all states are occupied. The change in density from the unperturbed situation is clearly ∆n(r) = n(r) − n 0 (r) where G 0 (r, r, z) is the Green's function for λ = 0. It is the screening density ∆n(r) that we wish to investigate in the high density limit. Since the potential is multiplied by the coupling constant λ, the density as given by Eq. (16) depends implicitly on this parameter. The variation of the induced density with λ, holding E F fixed, is simply By making use of the Dyson equation we find ∂Ĝ(z, λ) ∂λ =Ĝ(z, λ)VĜ(z, λ) .
Thus, the variation of the induced density can be expressed in the form We note at this point that this quantity still depends implicitly on λ through the full Green's function G(r, r ′ , z). We now make use of the spectral representation of the Green's function in Eq. (13). The product ψ i (r)ψ * i (r ′ ) can be treated as real since a complex wave function at the energy E i must have a time-reversed pair which is its complex conjugate. Thus, we have This expression is formally exact and is the crucial step in the argument needed to arrive at the final conclusion. It is clear that the only terms which give a nonvanishing contribution are those with E i < E F and E j > E F , or vice-versa. It is these terms which give a finite result in the E F → ∞ limit. If this limit is taken first for a given (i, j), the result is zero. In other words, it is important to sum over all states before the E F → ∞ limit is taken. We also note that the i = j terms in Eq. (21) do not present any difficulty despite the appearance of the energy denominator. From Eq. (20), these terms contribute which is well-defined. The energy integral in Eq. (21) can be written alternatively as Using this in Eq. (21), we find that Eq. (20) is equivalent to This is a more useful form for taking the E F → ∞ limit since now the energy argument of the Green's function is large. If E F ≫ |V (r)| everywhere (we assume here a bounded potential 27 ), the effect of the potential on the propagation of the electron can be neglected and the Green's function G can be replaced by the free-particle propagator G 0 . We then have where we have simply reversed the earlier argument to change the limits of integration in the last line. Several comments are in order: 1. The replacement of G by G 0 is only valid if the energy argument is large; these Green's functions are certainly different at low energies. In fact, the perturbed and unperturbed systems have completely different state spectra and H may contain bound states which do not appear in H 0 . Nevertheless, since Eq. (20) can be written in the equivalent form (24), the replacement of G by G 0 in Eq. (25) is justified in the E F → ∞ limit.
2. The final line in Eq. (25) is identical to the result obtained by linear response theory. Since it is independent of λ, it can be integrated from λ = 0 to λ = 1 to give 3. The final result in Eq. (26) does not depend on the assumed dimensionality and is valid even when the potential V (r) supports bound states. These bound states of course do not appear in the spectrum of G 0 but are nevertheless included in ∆n(r). The result given in Eq. (26) will be referred to as the high density screening theorem.
As an application of Eq. (26) we consider an ideal uniform gas in D dimensions. The use of the spectral representation for the Green's function gives in this case where the factor of 2 is for spin, L is the size of the sample in a given direction, f (E) is the zero-temperature Fermi function and χ 0 (q) is the static free particle density response function Substitution of Eq. (27) into Eq. (26) gives This is a rigorous statement that linear response gives the correct induced density in the high density limit in any dimension.
The static response functions in the various dimensions are given by the expressions In each case, the q → 0 limit of χ 0 (q) is the density of states g(E F ) at the Fermi level If the potential is smooth we can suppose that V (q) is effectively zero for q > q max with q max ≪ k F . In this case, Eq. (29) reduces to which is recognized as the linearized TF approximation for the density. It should be noted that 2D is special for a number of reasons. First, since E F is proportional to the density n, the full TF approximation is itself linear. Second, the free-particle response function has the constant value χ 0 (q) = 1/π for all wavevectors up to 2k F . Thus, the replacement of χ 0 (q) by χ 0 (0) is in some sense a less restrictive approximation than it is in other dimensions. In particular, if the only finite Fourier components of V (r) occur for q < 2k F (which is always the case when k F → ∞), we have The analogous result in 3D is ∆n(r) = 0 while the 1D result does not have a finite limit. Thus the induced density only has a nontrivial limiting value in 2D. The result for the 2D screening charge in Eq. (33) is all the more remarkable if we note that the potential V (r) can in fact support bound states, as it must in 2D if the potential is purely attractive. It is worth verifying this explicitly by means of a numerical example. We consider the following model potential which has no physical basis but is chosen simply to illustrate the point. It is a repulsive double barrier potential for V 0 > 0 and an attractive double well potential for V 0 < 0. With this potential we calculate the induced density by solving the Schrödinger equation for all states below E F . In Fig. 12 we show results for a repulsive potential with V 0 = 0.125 H and r 0 = 5 a.u. as a function of gas density. The screening charge is due to continuum states in this case and it can be seen that by r s = 0.5, the induced density is very close to the prediction of Eq. (33). Similarly, if an attractive potential is placed in a gas with r s = 0.5 (E F = 4.0 H), we obtain the bound and continuum charge densities shown in Fig. 13. Here, the magnitude of |V 0 | was adjusted to vary the number of bound states. Two cases were considered: (a) V 0 = −0.125 H for which a single m = 0 bound state is found and (b) V 0 = −0.25 H, for which two m = 0 bound states exist. The bound and continuum densities are quite different for these two cases, but the total density agrees very well with the result in Eq. (33), even though the high density limit has only marginally been reached. These examples confirm very nicely the behaviour stipulated by the high density screening theorem.
We can now apply the theorem to explain some of the results found in the Z = 1 calculations of Sec. III. Of course interactions were included there, but the Kohn-Sham system of equations correspond to a system of noninteracting electrons moving in the total effective potential in Eq. (2). In Fourier space, we have ∆v eff (q) = v ext (q) + v(q)∆n(q) + v ′ xc (n 0 )∆n(q) .
We have here linearized the xc potential since, in the high density limit of interest, the induced density is small compared to the background density. Identifying ∆v eff (q) as given by Eq. (35) with V (q) in Eq. (29), we have ∆n(q) = −χ 0 (q)∆v eff (q) (36) which implies This is simply the induced density obtained within linear response theory, treating Hartree and xc interactions at the mean-field level. An example of the linear response result was already given in Fig. 9 for the case of a Z = 1 charge outside the plane of the 2DEG. The close agreement of this density with the nonlinear screening density is confirmation of the high density screening theorem. For high densities, v xc is in fact dominated by the exchange potential v x = −C x /r s and we have v ′ xc (n 0 ) ≃ −(πC x /2)r s . Thus in the limit r s → 0, the xc term can be dropped and the induced density is simply given by the result The corresponding effective potential is The real space quantities are obtained from these expressions by an inverse Fourier transform. We now consider the particular example of a point charge Z in the plane of the gas. With v(q) = 2π/q in 2D, Eq.(39) becomes and in real space we obtain the TF potential It is a monotonic potential which behaves as −Z/r for r → 0 and falls off as −Z/4r 3 for r → ∞. Being negative definite, it will support at least one bound state for any value of Z. For Z = 1 the bound state energy is found to be E 0 = −0.2862 H which agrees with the value found by Stern and Howard 8 . This is the limiting value of the bound state eigenvalue shown in Fig. 3.
We thus come to the conclusion that the r s → 0 limit of the Kohn-Sham effective potential is the TF potential in Eq. (41). To illustrate this we show in Fig. 14 the induced density for Z = 1 and a very high density corresponding to r s = 0.2. Also shown are the bound and continuum components which add to the total. The latter is also displayed as r∆n(r) in Fig. 15 together with the TF density as obtained from Eq. (38). It can be seen that the DFT density simply oscillates about this limiting value. With decreasing r s , the amplitude of these oscillations decrease and the DFT density uniformly approaches the TF density.
For completeness, we note that the result corresponding to Eq. (40) in 3D is with q T F = 4k F /π the Thomas-Fermi wavevector. The real-space potential in this case is With increasing k F , the range of the potential decreases. Thus at some point the potential will no longer support a bound state, and all states in the presence of the external charge will be continuum scattering states. As alluded to earlier, the screening of the external potential becomes more effective with increasing density. This is in marked contrast to the situation in 2D where the screened potential has a density independent limit which continues to support bound states in the high density limit. In summary, if a potential V (r) is introduced into a 2D gas and the induced density is constructed by summing up over all eigenstates of the hamiltonian H (including all possible bound states), the induced density will be given exactly by Eq. (33) in the high density limit. Remarkably, this density is the same as what one would obtain by treating the potential as a weak perturbation and applying perturbation theory.

V. CONCLUSIONS
In this paper we have provided a detailed discussion of the nonlinear screening of an external point charge by a two-dimensional electron gas. Nonlinear effects are most evident when the charge lies in the plane of the 2DEG but rapidly decrease in importance as the charge is moved out of the plane. As an application of the theory we have considered the problem of 2D stopping power and have compared our results with earlier work. Quantitative differences point to the importance of performing self-consistent nonlinear screening calculations.
We also derive what we refer to as the high density screening theorem. In 2D, the theorem has the interesting consequence that the screened potential has a density independent limiting form in the high density limit and that the screening density is directly proportional to this potential. Explicit calculations for model potentials confirm these conclusions and demonstrate the surprising fact that linear response theory is valid even when the potential acting on the gas supports bound states. These nonintuitive results highlight some of the peculiarities of electronic screening in two dimensions. The solid curve is the nonlinearly screened potential including exchange-correlation, vext + ∆vH + ∆vxc; the short-dashed curve is the Hartree potential, vext + ∆vH , calculated with the same screening density as used for the solid curve ; the chain curve is the linearly screened potential within the RPA (i.e., the Hartree approximation); and the long-dashed curve is the Hartree potential, vext + ∆vH , obtained with a linear screening density that includes local-field effects (i.e., exchange-correlation).