{'fraction_non_alphanumeric': 0.03221482507484027, 'fraction_numerical': 0.008846789687681515, 'mean_word_length': 4.61938237509415, 'pattern_counts': {'":': 0, '<': 4, '<?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': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Truly self-consistent solution of Kohn-Sham equations for extended systems with inhomogeneous electron gas

The density functional approach in the Kohn-Sham approximation is widely used to study properties of many-electron systems. Due to the nonlinearity of the Kohn-Sham equations, the general self-consistent solution method for infinite systems involves iterations with alternate solutions of the Poisson and Schrödinger equations. One of problems with such an approach is that the charge distribution, updated by solving the Schrodinger equation, may be incompatible with the boundary conditions of the Poisson equation for Coulomb potential. The resulting instability or divergence manifests itself most appreciably in the case of infinitely extended systems because the corresponding boundary-value problem becomes singular. In this work the stable iterative scheme for solving the Kohn-Sham equations for infinite systems with inhomogeneous electron gas is described based on eliminating the long-range character of the Coulomb interaction, which causes the tight coupling of the charge distribution with the boundary conditions. This algorithm has been previously successfully implemented in the calculation of work function and surface energy of simple metals in the jellium model. Here it is used to calculate the energy spectrum of quasi-two-dimensional electron gas in the accumulation layer at the semiconductor surface n-InAs. The electrons in such a structure occupy states that belong to both discrete and continuous parts of the energy spectrum. This causes the problems of convergence in the usually used approaches, which do not exist in our case. Because of the narrow bandgap of InAs, it is necessary to take the nonparabolicity of the conduction band into account; this is done by means of a new effective mass method. The calculated quasi-two-dimensional energy bands correspond well to experimental data measured by the angle resolved photoelectron spectroscopy technique.


Introduction
The density functional approach in the Kohn-Sham approximation is widely used to study properties of many-electron systems. The set of equations to describe the nonuniform electron gas of infinite system in the Kohn-Sham approximation should be satisfied by the self-consistent distribution of electron density N (r) and Coulomb potential U (r). Due to the nonlinearity of the Kohn-Sham equations, the general method of self-consistent solution involves iterations with the alternative solving of the Poisson and Schrödinger equations [1]. One of problems of such an approach is that the charge distribution renewed by means of the Schrödinger equation solution may be incompatible with boundary conditions on Poisson equation for Coulomb potential, as it has long been known since the articles of Bardeen [2] and Lang and Kohn [3]. The resulting instability or even divergence of iterations manifests itself most appreciably in the case of infinitely extended systems because the corresponding boundary-value problem becomes singular. It was shown in [1,4], the published attempts to deal with this difficulty are reduced in fact to cancellation of iterative method declared initially, and its replacing somehow approximate scheme of calculation. In general, these schemes are usually semi-empirical and do not permit to evaluate the extent of deviation from the exact solution.
In [1] the iterative algorithm of solving the Kohn-Sham equations for unbounded systems with inhomogeneous electron gas has been formulated, which is based on eliminating the longrange character of Coulomb interaction being a cause of tight coupling of the charge distribution with boundary conditions. This algorithm was applied to a semi-infinite degenerate electron gas, limited infinitely high potential wall or self-consistent potential barrier [4]. The work function and surface energy of simple metals in the jellium model have been calculated and successfully compared with the published experimental and computational results [5,6].
Here we present the results of applying our algorithm to the calculation of the energy spectrum of quasi-two-dimensional (2D) electron gas in accumulation layer on the semiconductor surface. In such a structure, the states electrons are occupied belong both discrete and continuous parts of the energy spectrum. Therefore, the high density of electrons that bound near the surface and degenerate electron gas in the bulk of semiconductor coexist simultaneously. The stable convergence of iterations was obtained even in this case, in which the self-consistency previously considered unattainable [7]. To assess the adequacy of the method at the real structures, we calculated the energy spectrum of the quasi-2D electron gas in the case of the accumulation layer at the surface of n-InAs. The calculated energy levels correspond well to experimental data measured by the angle resolved photoelectron spectroscopy (ARPES) method [8]. In addition, the difference has been estimated between self-consistent solutions in the Hartree approximation and with taking into account the exchange-correlation interaction.

Description of the computational scheme
In this section, we describe a theoretical framework of our analysis. The atomic system of units |e| = m = = 1 is used everywhere based on the charge |e| and mass m of the free electron, unless otherwise specified. Therein the unit of length is the Bohr radius a B , and the unit of energy is Hartree Ha = e 2 /a B .
In the Kohn-Sham approximation, the gas of interacting electrons is described by the set of equations for single-particle wave function Ψ ε Here ε is the energy eigenvalue of single-particle state. The effective potential U eff (r) is the sum of Coulomb U and exchange-correlation U xc potentials The Coulomb potential energy of electrons obeys Poisson equation and boundary conditions corresponding to the posed problem where N + (r) is the density of positive charge, N (r) is the electron density. The exchange-correlation potential U xc is taken in the local approximation of the density functional where ϵ xc = ϵ x + ϵ c is the sum of exchange and correlation energy per particle. The density of electrons is determined through the single-particle wave functions by the expression where factor 2 takes into account the spin degeneracy of single-particle states.
As it was pointed in Introduction, the boundary conditions to the Poisson equation impose certain requirements on the spatial distribution of electron density. The result of integration in Eq. (5) does not necessarily compatible with these requirements until the self-consistency is reached. However, the incompatibility problem of boundary conditions with the right-hand side of Poisson equation does not arise, if to present the total density of electrons as the sum Here the induced charge density N ind depends explicitly on the unknown potential. As a result, the long-range Coulomb interaction is replaced with the screened interaction of finite range. The induced density of electrons is related with the effective potential by an approximate local expression that is usual for the semiclassical degenerate Fermi gas The function N qu (r) in Eq. (6) determines a quantum correction to the semiclassical N ind (r). For greater transparency from here on we restrict ourselves to the approximation of the Hartree's self-consistent potential U . The full description of the algorithm with allowance for the exchange-correlation potential U xc can be found in [4]. However all numerical results presented here were obtained by using the complete system of the Kohn-Sham equations including U xc .
The separation of the induced density transforms the Poisson equation (3) to the form The quantum correction N qu (r) to the electron density introduced here by means of Eqs. (5) -(7) is defined as The algorithm of consecutive steps i = 0, 1, 2,... for iterative solving of equations (1) and Eq. (8) can be now formulated as follows 1. One neglects small-scale quantum variations of electron density at the initial step of iterations 2. Solve the nonlinear Poisson equation with the known right-hand side and the imposed boundary conditions 3. Find the total electron density self-consistently with U (r) 4. Find the eigenfunctions of the single-particle Kohn-Sham Hamiltonian 5.Find anew the total electron density, using the obtained eigenvalues and eigenfunctions 6. Compute the density N qu (r) of quantum correction for the next step of iteration, using the found total electron density 7. Return to solving the Poisson equation with the new function N i+1 qu (r) in the right-hand side The convergence criterion for this iteration process is the value of δ, which measures the maximal difference between the electron density N i s (r), Eq. (12), obtained after solving the Poisson equation, and the density N i (r), Eq. (14), found by means of the Schrödinger equation

Work function of simple metals
The efficiency of the proposed algorithm has been successfully tested with the self-consistent calculations of the electronic properties of a number of typical unbounded systems [4]. From the results published out there we will present here only the part related to the calculation of the work function of simple metals in the jellium model. In the approximation of a homogeneous positive background, the work function of simple metals was calculated in many papers. However the published data reveal sometimes not only quantitative, but also qualitative discrepancy. As a reason we believe the complete selfconsistency of electron density and effective potential have not been reached in these works. The purpose of this part of the present work is to demonstrate the truly self-consistent solution to one of the classical problems of the theory of inhomogeneous electron gas, and appreciate the importance of self-consistency by comparing the results with those obtained by other methods in the same formulation of the problem. We have had no aim to improve the jellium model or the density functional theory for calculating the surface properties of simple metals. Our goal was to find an exact solution of the known problem in its conventional statement.
Let the positive background occupy the half-space z > 0, where 0 is coordinate at the metalvacuum boundary. Following [4], we introduce dimensionless variables The boundary condition for the Poisson equation were taken as Hence, the neutral in a whole metal is considered. Wave functions Ψ ε of single-particle states exponentially damped at ζ → −∞ for states with ε ≤ ε F . In the bulk of metal the eigenfunctions Ψ k have the asymptotic form where dimensionless quantum number k = √ 2ε is normalized to k F . The eigenfunctions (18) are normalized to delta function δ(k − k ′ ).
The self-consistent electron density and effective potential have been obtained after convergence the offered iteration process at various values of Wigner-Seitz radius 0.3 ≤ R s < 5 [4]. To check the degree of self-consistency of the obtained solutions, we have used the exact strict relationship for the self-consistent electrostatic potential derived in [9]. It relates the equilibrium potential difference between the bulk of metal and the edge of the positive background with a derivative of the energy per electron with respect to the background density where ϵ = 3 5 ε 0 F + ε xc is the total mean energy per electron in the homogeneous electron gas. Contrary to the work [6] where the criterion (19) was explicitly used in the process of calculations, this relationship was not taken into account in our procedure. So, we can use Eq. (19) to verify if the calculated potential is actually self-consistent with the electron density. Good agreement between ∆ = u(0) − u(∞) and ∆ BV =∆ BV /ε 0 F should prove there was obtained the truly self-consistent solution of Kohn-Sham equations.   [4] with published experimental and theoretical data is shown in Figure 1. Here the dependence W (R s ) found in [4] (curve 3) is compared to similar ones published in the papers [5] and [6]. For R s > 2 there is semi-quantitative agreement with those obtained in [5] (curve 1) and [6] (curve 2. One can see from the figure that both curves 2 and 3 have maxima at R s ≈ 1.7, however, at small R s there is a qualitative difference. According to [6], the work function should have had a local minimum at R s ≈ 1, which is not confirmed by our calculation. This local minimum has possibly resulted from poor convergency of iterations during the solution of Kohn-Sham equations at small R s , which has been declared by authors of the article [6]

Energy spectrum of quasi-two-dimensional electron gas in accumulation layer at the semiconductor surface of n-InAs.
So far, we have had to deal with systems where only the states of continuous spectrum were occupied by electrons. Now let us consider the infinite systems, in which the localized states of the discrete spectrum coexist with delocalized states of the continuous spectrum. An example is the electron gas in the accumulation layer at the surface of semiconductor. As it was noted in the Sec. 3.1 of our paper [4], approximate calculation methods, that have been reported in the literature for such systems (see [10] page 463 and [7,11]), do not lead to the self-consistent solution. Thus, it was natural to see whether our algorithm can be applied in such a problem. First of all, it was necessary to make sure that the algorithm works for large values of N qu , which, according to (9), have to originate from the large discrepancy between the number N ind of semiclassical and N 2D quasi-two-dimensional electrons inside the quantum well. Second, the problems of this kind are typical of the physics of semiconductor nanostructures. Therefore, it was necessary to ensure that the self-consistent calculations in such a statement can be made.
One of such systems is the structure of n-InAs with an accumulation layer at the semiconductor surface. In these structures there is an increased concentration of charge carriers near the surface, the movement of which is limited in the direction normal to the surface. We denote this direction as z-axis. In the directions parallel to the surface (x-y plane) carriers are  Fig. 2). Hence the induced density, defined by (7), can reach sufficiently high values, which may vary widely from the true self-consistent density. What could break the convergence of the proposed iterative algorithm. The first direct observation of surface-bound electronic states in structures based on n-InAs with accumulation layer was made by Tsui [12], using the tunneling spectroscopy method. However, such data are not very useful for comparison with our calculations, since the bias voltage on the structure varies during the measurement process and, accordingly, there is changing the electron energy spectrum. The development of the ARPES has enabled direct observation of these bound states and the two-dimensional electron dispersion law [8]. To make it easier to compare the results of calculations with measurements we have taken the numerical parameters of structures in accordance with the data in this paper. Figure 2 illustrates the band diagram of structure with an accumulation layer at the semiconductor surface. It represents schematically the downward bending of the conduction band and positions of the size-quantized energy levels. The shaded area corresponds to the states of continuous spectrum, which are occupied by electrons of the degenerate gas.
Since the InAs is a narrow-bandgap semiconductor, it is important to take account for the nonparabolicity of the conduction band within calculations of energy spectrum of quasi-2D electron states. In the two-band approximation of the Kane's kp-theory [13,14], the energy spectrum ε(k) of narrow-bandgap semiconductors in the neighborhood of Γ-points of the Brillouin zone is determined by the equation where E g is the energy gap of the semiconductor, m * 0 is the effective mass at the bottom of the conduction band, and k = (k ∥ , k z ) is the quasi-wavevector of electrons. Two roots of this equation describe the energy of electrons in the non-parabolic conduction band and that of light holes in the valence band, respectively. At E g → ∞ these expressions become the usual one-band quadratic spectra of energy for the electrons and light holes. If the periodic potential of the semiconductor crystal is perturbed by smooth non-periodic potential U (r), the states of electrons and holes in this case are defined as eigenstates of the effective Hamiltonian and obey the effective mass equation. In the derivation of the effective Hamiltonian, either Wannier function expansion (see e.g. [15], Ch.6) or Luttinger-Kohn function expansion (see [16] and, e.g., [17], Ch.21) of the sought-for wave function is used.
In the one-band approximation, both methods lead to the same result. The effective mass equation can be obtained from the dispersion equation and considering left-hand-side as the characteristic polynomial of the operator ε 1 (−i∇ r ) ≡ k 2 /2m * 0 with the problem of eigenvalues and eigenfunctions in the form we can apply the same replacement rule (22) to the two-band dispersion equation (20). This leads to a scalar equation of the effective mass theory 1 instead of matrix equation for two envelope functions, which are usually considered [18]. In explicit form we obtain the following approximate equation for the eigenfunctions that can be seen directly from Eq. (20). In deriving this equation from Eq. (23), the commutator ) has been omitted. It is of the order of the square of the small parameter, (L g /L) 2 , where L g = 1/ √ m * 0 E g and L is the characteristic scale of U (r). Thus, the applicability of the effective mass theory in the two-band case requires additional condition L g /L < 1 for smoothness of the potential U , besides conventional a/L << 1, where a is the lattice constant of semiconductor. In the case of InAs, a ≃ 0.6 nm and the ratio L g /a ≤ 5, so the quantum well of the accumulation layer with L ≥ 30 nm should be good described in the effective mass approximation.
In order to find the boundary conditions for the solution of equation (24) one needs to know how that solution relates with the exact eigenstate Φ ε (r) of the conduction band in the presence of the potential U (r). The successive derivation of Eq. (24) is based on the expansion of Φ ε (r) in terms of the Bloch functions ϕ ck (r) of the conduction band and gives rise to the expression where w c (r) is the Wannier function of the conduction band. The expression (25) is simply the sliding smoothing of the function Φ ε (r) with a window in the form of Wannier function. The characteristic scale of the latter is of order of the lattice constant of semiconductor. Therefore, the boundary and the normalizing conditions for the smoothed solution Ψ ε of (24) should be the same as those for the exact solution Φ ε . The one-dimensional inhomogeneity of the accumulation layer allows one to write the threedimensional wave function Ψ(r) of electrons as the product of a plane wave that depends on the coordinate r ∥ parallel to the surface, and the function ψ ε (z) Further we use the dimensionless units n = N/N + , ζ = k F z,k = k/k F , ϵ = ε/ε 0 F , u eff = U eff /ε 0 F with the choice of parameters in the LDA of U xc as in [4]. The function ψ(z) satisfies to the one-dimensional Schrödinger equation of the form For local states of the discrete spectrum we have The effects of the nonparabolic conduction band should be included for solving Poisson equation as well. Eq. (8) can rewritten as where c n = (8/3π)(4/9π) 1/3 R s and induced density n ind is given by Here the Fermi energy ϵ F ̸ = 1 due to the U xc contribution and nonparabolicity of InAs conduction band. Potential u(ζ) must satisfy the boundary condition Writing of the Poisson equation in the form of (31) corresponds to the condition that the ion charge forms the homogeneous positive background of density N + (z) = N + for z > 0. Self-consistent electron density and effective potential has been obtained by an iterative solution of the equations (27) and (31). The proposed algorithm has provided the convergence of the process. The total density of electrons was given by the expression where k ∥max is a root of equation ε j (k 2 ∥ ) = ε F . The first term in (33) is the density of delocalized states n 3D and the second one is the density of quasi-two-dimensional electrons n 2D . Fig. 3 demonstrates the spatial dependence of n 3D (solid curve) and n 2D (dashed curve) for the case when there are two size-quantized bands in the quantum well. It can be seen, that in the region of the quantum well the density n 2D of localized electrons, which is poorly approximated by the semiclassical expression (32), is much higher than the density of electrons in the states of the continuous spectrum. Nevertheless, the algorithm under study has shown the stable convergence of iterations. Fig. 4 shows the energy dispersion of electronic state as a function of k || . We calculated the dispersion law E j (k) of quasi-two-dimensional subbands for the structure with the materials parameters presented in [8]. Circles in the Fig. 4 represent the measured energy of electrons for different values of k ∥ . They were extracted from the insert in Fig. 2 of the article [8]. Authors of this work have assumed that clean surface of the samples corresponds to the condition of flat zones. In this case the measured band bending at the surface covered with some layer of Cs has been estimated as U 0 = 400 meV. However, the best match of our results with experimental dispersion curve was achieved with U 0 = 435 meV (see Fig. 4).
We believe the band bending ∆U 0 = 35 meV for a clean surface could exist, but to be unobservable for two reasons: 1) the possible absence of a bound state in a shallow potential well, and 2) since this value is within the stated experimental error of order of 35 meV.  In the end we recall that all presented results of calculations were obtained with account the exchange-correlation interaction. It turned out, however, that for the semiconductor structures with parameters specified in [8], contribution to the energy of size-quantized subbands due to u xc is small and is approximately 3 meV.