Computation of the Spatial Distribution of Charge-Carrier Density in Disordered Media

The space- and temperature-dependent electron distribution n(r,T) determines optoelectronic properties of disordered semiconductors. It is a challenging task to get access to n(r,T) in random potentials, while avoiding the time-consuming numerical solution of the Schrödinger equation. We present several numerical techniques targeted to fulfill this task. For a degenerate system with Fermi statistics, a numerical approach based on a matrix inversion and one based on a system of linear equations are developed. For a non-degenerate system with Boltzmann statistics, a numerical technique based on a universal low-pass filter and one based on random wave functions are introduced. The high accuracy of the approximate calculations are checked by comparison with the exact quantum-mechanical solutions.


Introduction
2][3] The spatial and energy disorder creates a random potential, which decisively affects electron states.Among other effects, a disorder potential causes the spatial localization of electrons in the low-energy range.The knowledge of the space-and temperature-dependent electron distribution n(r, T ), particularly in localized states created by random potentials, is required to understand charge transport and light absorption/emission in disordered semiconductors.The distribution n(r, T ) is most straightforwardly obtained by solving the Schrödinger equation in the presence of a disorder potential.However, this procedure is extremely demanding with respect to computation facilities.It is hardly affordable for realistically large chemically complex systems.Therefore, it is highly desirable to develop theoretical tools to get access to n(r, T ) without solving the Schrödinger equation.
[6][7] In the LLT, the random potential is converted into some effective potential, which drastically simplifies all calculations.However, recent studies 8,9 revealed substantial problems of the LLT.For instance, the effective potential in the LLT lacks a temperature dependence that is necessary to describe n(r, T ) appropriately.Moreover, the LLT is seen to be equivalent to the Lorentzian filter applied to a random potential. 8Such a choice of the filter function is rather unfortunate.The Lorentzian filter yields a significantly larger number of localized states in a random potential than the number of such states obtained via the exact solution of the Schrödinger equation. 8Therefore, more developed computational techniques are desirable.
Here we develop two numerical techniques to reveal n(r, T ) in disordered systems under degenerate conditions controlled by Fermi statistics, avoiding the time-consuming numerical solution of the Schrödinger equation.One of the techniques is based on converting the Hamiltonian into a matrix, which, being subjected to several multiplications with itself succeeded by inverting the outcome, yields the distribution n(r, T ).The other technique replaces the operation of matrix inversion by solving a system of linear equations controlled by the matrix generated from the Hamiltonian.
We also describe two recently developed computational techniques for calculations of n(r, T ) in non-degenerate systems controlled by Boltzmann statistics. 8,9One algorithm is based on applying a temperature-dependent universal low-pass filter (ULF) to the random potential V (r).This yields a temperature-dependent effective potential, W (r, T ), that enables a quasiclassical calculation of particle density n(r, T ).The ULF algorithm employs Fast-Fourier Transformation for calculating the effective potential W , enabling the analysis of very large systems.
The other algorithm is based on the recursive application of the Hamiltonian to multiple sets of random wave functions (RWF) for a specific realization of the random potential V (r).
Following the repeated application of the thermal operator, the temperature-dependent electron density is determined by averaging the outcomes over different RWF sets.This procedure offers several advantages over the widely-used LLT.Unlike the LLT, which relies on an adjustable parameter that can only be determined through comparison with the exact solution, 8,9 the RWF scheme lacks any adjustable parameters and works across all temperatures.
Additionally, the accuracy of the RWF approach in computing n(r, T ) can systematically be improved, whereas the accuracy of the LLT is inherently limited.

Calculation of n(r, T ) in a degenerate system controlled by Fermi statistics
To be definite, we consider a disorder potential acting on electrons, characterized by Gaussian statistics ('white noise'), i.e., the potential obeys ⟨V (r )⟩ R = 0 with the auto-correlation where ⟨. ..⟩R indicates the average over many realizations R of the random potential and S is the strength of the interaction.This quantity S yields natural definitions for the characteristic length scale and for the characteristic energy scale in the form where d is the space dimensionality and k B is the Boltzmann constant.
We consider here a collection of non-interacting electrons in some external potential in the thermodynamic equilibrium that is characterized by the temperature T and the Fermi energy ε f .The goal is to develop an effective numerical method for the calculation of the electron density (concentration) n(r, T ) as a function of coordinates r in a degenerate system.
The electron density can be defined as In this equation, summation index a labels the electron eigenstates, i. e., solutions of the Schrödonger equation Ĥψ a = ε a ψ a , where ψ a (r) and ε a are the wavefunction and the energy of the eigenstate; f (ε) is the Fermi function, and factor 2 before the sum in Eq. ( 4) accounts for the two possible spin orientations.
Below we assume that the system under study is discretized with a finite-difference (or tight-binding) method.A wavefunction ψ(r) is therefore represented as a collection of probability amplitudes ψ 1 , ..., ψ L at L points (grid nodes) evenly distributed in space.The Hamiltonian Ĥ is a matrix of size L × L. Equation (4) for n(r, T ) in this discrete setting obtains the form where n i is the electron density at grid node i ∈ {1, ..., L}; ψ a,i is the value of eigenfunction ψ a at node i; ε a is the electron energy that corresponds this eigenfunction; and ∆V is the spatial volume per one grid node (in the one-dimensional case, ∆V is simply the distance between nodes).
Equation ( 6) represents a standard way for a numerical calculation of the electron density in degenerate systems.However this way requires the solution of the eigenvalue problem, which takes a large amount of computational resources for the large size L of the Hamiltonian matrix.Below we suggest two methods that allow one to speed up the calculation of n(r, T ).
In Method 1 (see Section 2.1), the numerical solution of the eigenvalue problem is replaced by a matrix inversion.The latter numerical task is much faster than the solution of the eigenvalue problem in the case of a sparse matrix, in which almost all entries are equal to zero.In Method 2 (see Section 2.2), we employ a numerical solution of a system of linear equation, which is even faster than the matrix inversion.Performance of Methods 1 and 2, as compared to the standard method based on the eigenvalue problem, is tested in Section 2.3 on an example of a one-dimensional disordered system with one occupied band.
The idea of Methods 1 and 2 is based on a simple observation that the shape of the function where a number N is large, resembles the shape of the Fermi function in the vicinity of point x = 1.2][13][14] Replacing x with an appropriate linear function of the Hamiltonian is the essence of Methods 1 and 2. The detailed justification is given in Appendix A.

Method 1: matrix inversion
The method is based on the Hamiltonian Ĥ (a matrix L × L), the Fermi energy ε f , and two additional parameters: a "reference energy" ε 0 and the number of iterations N .These parameters are related to the temperature T , Details for the choice of these parameters are given in Appendix B.
The method starts from composing the matrix Â0 from the Hamiltonian where Î is the L × L unit matrix.Then, the matrix Â0 is squared N times To the resulting matrix ÂN we add the unit matrix and invert the sum to get a new matrix Finally, the electron density is obtained from the diagonal elements B ii of the matrix B: 2.2 Method 2: solving a system of linear equation Similarly to Method 1, the input parameters are ε f , ε 0 and N , which are related to the temperature T by Eq. ( 8).In addition, this method requires one more parameter N C .The choice of N C is discussed in Appendix B. For a given N C , we compose a matrix Û of size L × N C that obeys the following conditions (see Fig. 1 for the shape of this matrix in the one-dimensional case): • each entry of matrix Û is equal to either 0 or 1; • in each row of matrix Û , exactly one entry is equal to 1; • in each column of matrix Û , the nodes with nonzero entries are placed spatially as far from each other as possible.For example, in the one-dimensional case, the unities in each column are separated by (N C − 1) zeros, as illustrated in Fig. 1.
The matrix ÂN is calculated, as done in Method 1, see Eqs. ( 9) and (10).However, in contrast to Method 1, no matrix inversion is necessary.Instead, a system of linear equations is to be solved with respect to an unknown matrix X of size L×N C .This is a computationally easier task in comparison to the matrix inversion.Finally, the electron density is calculated as

Numerical example: a one-dimensional disordered system with one occupied band
Let us compare the performance of different methods to calculate n(r, T ) in a simple onedimensional tight-binding model with energy bands and disorder.We consider a linear chain of L lattice nodes at a distance a = 0.1 from each other.We measure distances and energies in the units determined by Eqs. ( 2) and (3).In these units, the hopping integrals between neighboring nodes H i,i+1 and H i+1,i are equal to −1/(2a 2 ) = −50.The on-site energies H jj are chosen to be where ξ j are random numbers uniformly distributed in the range −1 < ξ j < 1.All other matrix elements of the Hamiltonian Ĥ are equal to zero.Periodic boundary conditions apply.
The constant term V 0 = 100 makes the lower boundary of the energy spectrum ε min to be near zero, and the higher boundary ε max to be ≈ 200.The amplitude of periodic variations The amplitude of the random-noise potential V 2 = √ 30 is chosen such that the "disorder strength" S is equal to unity.As an example, we show one realization of the on-site energies in Fig. 2a.As a consequence of the periodic potential (the second term in the r.h.s. of Eq. ( 17)), the electron energy spectrum consists of the four energy bands separated from each other by band gaps.We consider the situation when the lowest energy band is completely occupied by electrons, and three other bands are empty (quarter filling).The electron density in such a case is expressed as where summation is performed over the eigenstates of the occupied band.
The wave functions ψ a , that enter Eq. ( 18), can be obtained by diagonalizing the Hamil- tonian.An example of the calculated electron density distribution is shown in Fig. 2b by blue lines.
Also we show in Fig. 2b the approximated electron density obtained by Method 1 (red dots) and by Method 2 (orange circles).The parameters for these methods are: ε f = 28.5 (the middle of the lowest band gap), ε 0 = 10, N = 3, and N C = 30.One can see that both Methods provide quite accurate results for the electron density.In Fig. 3 we compare the time required by three ways of calculating the electron densitythe exact method based on the Hamiltonian diagonalization, Method 1, and Method 2-on a desktop PC with MATLAB used for the matrix manipulations.One can see that Method 1, which employs matrix inversion, is approximately one order of magnitude faster than the usual method of Hamiltonian diagonalization.Method 2 provides additional speedup by more than an order of magnitude.
6][17] For the sake of simplicity, we use here the standard MATLAB functions in Method 1.Even in such a non-optimized setting, this Method demonstrates a substantial speedup in comparison with matrix diagonalization.
3 Calculation of n(r, T ) in a non-degenerate system controlled by Boltzmann statistics 3.1 Low-pass filter (LF) approach

Motivation
Already in the 1960s, Halperin and Lax recognized that electrons in a random potential cannot follow very short-range potential fluctuations. 10This effect is illustrated in Fig. 4, where the random white-noise potential V (x) in one dimension is depicted by the green solid line.The detailed shape of V (x) in the region 275 ≤ x ≤ 325 (in the dimensionless units given by Eq. ( 2)) is compared with the shape of the wave function ψ for the lowest energy state in this spatial region.Apparently, the characteristic width ℓ wf of the wave function, even for low-energy localized states, is substantially larger than the spatial scale of the fluctuations of the disorder potential V (r).The latter scale in semiconductor alloys is of the order of the lattice constant a ≈ 0.5 nm.The strong inequality ℓ wf ≫ a suggests that electrons in localized states are affected only by the mean disorder potential, averaged over the space scale ℓ wf .It is, therefore, not necessary to solve exactly the Schrödinger equation with the real disorder potential in order to get access to the individual features of electron states.Instead, one can apply to the disorder potential V (r) a low-pass filter 10 (LF) that smoothes the spatial fluctuations of V (r).
Halperin and Lax suggested the square of the wave function as the filter function. 10e width of the filter function, ℓ ≈ ℓ wf was adjusted dependent on the energy, ε, of the localized state, ℓ ∝ |ε| −1/2 , where ε is counted from the band edge in the absence of disorder.
A variational approach was used to determine the shape of the low-energy density of states in a random potential. 10Baranovskii and Efros 18 addressed the same problem by a slightly different variational approach and confirmed the result of Halperin and Lax.
Neither of the two groups considered, however, the individual features of localized states, being focused solely on the structure of the density of states in the low-energy region. 10,18r aim here is, in contrast, to calculate the spatial distribution of electron density, n(r, T ), in a given disorder potential V (r).We start below with the definition of the low-pass filter and then formulate the algorithm for the calculation of n(r, T ).
Figure 4: Realization for the white-noise disorder potential on a one-dimensional strip.Insert: the wave function ψ(x) of the state with the lowest energy in the region 275 ≤ x ≤ 325.The coordinate x is dimensionless in the units given by Eq. (2).

Definiton of a low-pass filter (LF)
In one dimension, the low-pass filter (LF) is determined by the operation where the filter function Γ should contain the appropriate length scale ℓ.This operation converts the real disorder potential V (x) into the smooth effective potential W (x).For instance, one can try a Lorentzian function Γ L , [6][7] Halperin and Lax 10 suggested instead to use for LF the square of the wave function The shape of the wave functions ψ(x) for the low-energy states was determined by the optimal-fluctuation-approach that yields the filtering function where the characteristic length ℓ HL should depend on the state energy. 10,18Remarkably, it appears that a universal, energy-independent value for ℓ HL , can be introduced, 8 ℓ HL = 0.76 ℓ 0 as evidenced in Fig. 5a, where the effective potential yielded by the filtering function given by Eq. ( 21), is compared with the positions and energies of the eigenstates.The eigenstates for the given realization of disorder potential V (x) were obtained via a straightforward solution of the Schrödinger equation.In Fig. 5, the 30 eigenstates with the lowest energies are depicted by red points.The excellent agreement between the local minima of the effective potential (shown by the solid green line) and the positions and energies of the exactly calculated eigenstates justifies the filter function given by Eq. (21).In Fig. 5b, such a comparison is illustrated for the case of a Lorentzian filter, determined by Eq. ( 20) with ℓ L = 0.27 ℓ 0 chosen to mimic the LLT result. 8Evidently, the choice of a Lorentzian filter function is not satisfying.The number of the local minima in the effective potential W (x) (shown by the solid blue line) is significantly larger than the true number of the exactly calculated eigenstates.This happens because the Lorentzian function given by Eq. ( 20) is not smooth at x = 0.The cusp at x = 0 filters too many extrema from the real disorder potential V (x), preventing the identification of true localized electron states by searching the minima of the effective potential W (x). The filter function suggested by Halperin and Lax 10 does not possess such a deficiency.
Not only the energies and spatial positions of localized states in disorder potential V (r), discussed above, are of interest for the theory.In fact, the key quantity for the optoelectronic properties of disordered semiconductors is the space-and temperature-dependent electron distribution n(r, T ).Below we extend the LF approach to calculate n(r, T ).For that purpose, we introduce the temperature T into the filter function and, concomitantly, into the definition of the effective potential W (r, T ) that yields n(r, T ).

Universal filter function to determine n(r, T )
The T -dependent spatial distribution of electron density n(r, T ) is related to the quasiclassical effective potential W (r, T ) as where N c is the effective density of states in the conduction band and µ is the chemical potential.This equation serves as the definition of the quasi-classical effective potential W (r, T ).This effective potential is smooth in comparison to the initial disorder potential V (r) because W (r, T ) is derived from the electron density n(r, T ), which has the spatial scale of the electron wave functions, i.e., is broader than the scale of the short-range fluctuations of V (r).Let us, therefore, obtain W (r, T ) by subjecting V (r) to the action of a universal low-pass filter (ULF).
The key question is how to find out the appropriate T -dependent filter function Γ(r, T ) that can be used to extract the shape of the effective potential W (r, T ) for a given realization of the white-noise potential V (r).One can show that the Fourier image of this function has the shape 9 where erfi is the imaginary error function, and λ = h/ √ 2mk B T .
In order to reveal the electron density distribution n(r, T ) for a given realization of the white-noise potential V (r), one should first calculate the Fourier image V (k) of V (r) using a fast-Fourier-transform (FFT).This function V (k) should be then multiplied by Γ(|k|), given by Eq. ( 23).The inverse Fourier transform of the product V (k) Γ(|k|) by the FFT yields the effective potential W (r, T ) because the inverse Fourier transform converts a product into a convolution. 9Inserting W (r, T ) into Eq.( 22) gives the electron density n(r, T ).In Fig. 6 and Fig. 7, we compare the results for W (r, T ) of the above procedure with the effective potentials obtained via Eq.( 22) from the electron density n(r, T ) calculated using the exact solution of Schrödinger equation in one and two dimensions, respectively.Coordinate x and temperature T in the figures are measured in units given by Eqs. ( 2) and (3).
The data in Figs. 6 and 7 demonstrate the high accuracy of the approach based on the T -dependent low-pass filter.Notably, the FFT operation used in this approach does not need a considerable computer time, in contrast to the exact calculations on the basis of Schrödinger equation.3.2 Random-wave-functions (RWF) approach to calculate n(r, T )

Background
The idea of the RWF approach resembles the one suggested recently by Lu and Steinerberger 19 to search for the low-lying eigenfunctions of various linear operators.An iterative application of the operator leads to the increasing contributions of the low-energy regions to the state vector. 19A similar approach has been suggested by Krajewski and Parrinello 20 for the calculation of the thermodynamic potential.
Let us consider the action of the operator ĥ = exp[− Ĥ/(2k B T )] on an arbitrarily chosen wave function.The goal is to model the equilibrium distribution of electrons, which is described by the Boltzmann statistics in the nondegenerate case considered here.In Boltzmann statistics, states with energy ε contribute to the distribution with the probability The wave function is the probability amplitude, which explains the factor 1/2 in the operator ĥ.The wave function is always a linear combination of eigenfunctions that correspond to different energies.The action of the operator ĥ suppresses the contributions of high-energy eigenfunctions in favor of the contributions of low-energy eigenfunctions.By the application of the operator ĥ to a collection of the random wave functions, the average contributions of eigenfunctions corresponding to different energies approach their distribution in thermal equilibrium.The averaging here is performed over the set of the random wave functions.Physically, this procedure corresponds to the averaging of the electron density n(r, T ).The question arises on how to numerically subject a wave function to the action of the operator ĥ = exp[− Ĥ/(2k B T )].It can be done by recursively applying the Hamiltonian Ĥ to the wave function: with a natural number 9 M ≈ 1/(2αk B T ) and a small parameter α.A simple analysis justifies the choice 9 α = 1.5/ϵ max , where ϵ max is the estimate for the upper boundary of the energy distribution.In the case of a regular grid with the lattice constant a, ϵ max = h2 /(ma 2 ).
Below we describe how to realize this idea technically.

The RWF algorithm
Let us consider the RWF algorithm on a spatial lattice with the volume ∆V per lattice site.The value of the random wave function ψ on each lattice site is chosen independently as a random number extracted from a normal distribution with the average value zero and variance 1/∆V .The following transformation of the wave function ψ, is applied M times.Then an estimate of the reduced electron density ñR (r, T ) is The calculation of ñR (r, T ) is carried out for a large number N R of realizations R of the random wave function ψ(r).Then, the electron density n(r, T ) is the arithmetic mean of the functions ñR (r) obtained for different realizations R, multiplied by a chemical-potentialrelated factor e µ/(k B T ) , The larger the number of realizations N R , the more accurate is the calculated electron density n(r, T ).
In Fig. 8 and Fig. 9 the reduced electron densities, obtained via the RWF algorithm are compared with the exact ones calculated using the solution of the Schrödinger equation for one and three dimensions, respectively.Evidently, the RWF algorithm with 1000 iterations accurately yields the electron density in a wide range of temperatures.

Discussion
In this work, we introduce four theoretical tools to get access to the space-and temperaturedependent electron density n(r, T ) in disordered media with a random potential V (r), avoiding the time-consuming numerical solution of the Schrödinger equation.
For the case of degenerate conditions controlled by Fermi statistics, the Hamiltonian is converted into a matrix, which, being subjected to several multiplications with itself succeeded by inverting the outcome, yields the distribution n(r, T ).The other possible technique for the case of Fermi statistics replaces the operation of matrix inversion by solving a system of linear equations controlled by the matrix generated from the Hamiltonian.
For non-degenerate conditions with Boltzmann statistics, the universal low-pass filter (ULF) approach and the random-wave-function (RWF) algorithm are suggested for approximate calculations of n(r, T ).Both methods require far less computational resources than the complete solution of the Schrödinger equation.
The ULF approach employs the temperature-dependent effective potential W (r, T ).This technique is based on the Fast Fourier Transformation, which does not impose any demands on computational resources, such as processor time and memory.Therefore, it can be applied to mesoscopically large three-dimensional disordered systems.
Being superior to the widely used approximate methods, the RWF is computationally more costly than the ULF approach, if mesoscopically large three-dimensional systems at low temperatures are addressed.However, the accuracy of calculations based on the RWF algorithm can be unlimitedly improved by increasing the number of the RWF realizations.
Figure 10: Illustration of set S i for a given row index i of matrix Û .White and black squares represent zeros and ones, respectively.Row i and column a of the matrix are highlighted in green.
Note that i ∈ S i since U ia = 1.Therefore the sum in the right-hand side of Eq. ( 41) contains the term B ii .
Let us argue now that, at large enough N C , all other terms B ij in this sum with j ̸ = i are negligible.One can see from Eq. (36) that the matrix B is close to the projector P to the set of eigenstates with energies below the Fermi energy.Typically, especially in the situation when the filled states are separated by a band gap from the empty ones, the matrix elements P ij of this projector decrease with increasing the distance between nodes i and j, and become negligible at some distance.Hence, the matrix elements B ij in the right-hand side of Eq. (41) becomes negligible if nodes i and j are far away from each other, i. e., if i ̸ = j.The only non-negligible term is B ii , and therefore, Substituting there B ii from Eq. (38), one arrives at Eq. ( 15) of Method 2.
Similarly, the combination of Eqs.(39) and (42) justifies Eq. ( 16) of Method 2 in the case B Choice of parameters ε 0 , N and N C There are two restrictions on the choice of parameters ε 0 and N and, hence, on the temperature value T that can be achieved with Methods 1 and 2. First, the approximate equalities (34) and (35) must hold in the whole range of electron energies.And second, the matrix ( ÂN + Î) must not be ill-conditioned.
To consider the first restriction, let us define the function f (ε) As an example, Fig. 11 shows f (ε) along with the Fermi function f (ε) for N = 4 and two choices of the "reference energy" ε 0 : ε 0 < ε f (upper panel) and ε 0 > ε f (lower panel).
In both cases, there is a discrepancy between the functions f (ε) and f (ε) below the energy 2ε 0 − ε f in the first case, and above the energy 2ε 0 − ε f in the second case.Electron energy levels must not fall into these areas of discrepancy, otherwise the contributions of these levels into the electron density would not be accounted correctly in Methods 1 and 2. Therefore, in the case of ε 0 < ε f , the lowest electron energy ε min must be larger than 2ε 0 − ε f .Similarly, in the case of ε 0 > ε f , the highest electron energy ε max must be smaller than 2ε 0 − ε f .These conditions can be rewritten as restrictions to the "reference energy" ε 0 : Now let us consider the second restriction.The inversion of the matrix ( ÂN + Î) in Method 1, or solving a system of linear equations expressed by this matrix in Method 2, is possible when this matrix is not ill-conditioned.This means that the condition number, a ratio of the largest and the smallest eigenvalues of the matrix, is less than the maximal value of the order of 10 15 (a number that corresponds to the accuracy of representation of real numbers in the computer memory).The minimal eigenvalue of matrix ( ÂN + Î) is close to unity, and corresponds to energy levels near to ε 0 .The largest eigenvalue corresponds to the energy level farthest from ε 0 , i. e., either to ε min or to ε max , and is approximately the largest of two values [(ε min − ε 0 )/(ε f − ε 0 )] 2 N and [(ε max − ε 0 )/(ε f − ε 0 )] 2 N .Hence, the parameters ε 0 and N have to obey the following restrictions: The choice of parameters ε 0 and N is based on equations (45) and (46).We consider two different options: (i) when the temperature T is given, and (ii) when the goal is to obtain the sharpest possible boundary between filled and empty states.
In the first option, one should try the natural numbers in ascending order as values of N .
For each such number N , one finds ε 0 according to Eq. ( 8) and check whether conditions (45) and (46) are fulfilled.The smallest suitable value of N is the best choice.Indeed, the smaller N is, the more sparse matrix ÂN + Î is and, consequently, the faster is matrix inversion in Method 1 or solution of linear equations in Method 2.
In the second option, the largest possible N is desirable, in order to minimize the temperature T according to Eq. ( 8).To achieve the largest T , one has to choose the value ε 0 that maximizes the denominator |ε f − ε 0 | in Eq. ( 46): for the first and the second choice of ε 0 in Eq. (47), correspondingly.
Finally, let us consider the choice of the parameter N C in Method 2. The accuracy of Method 2 improves with N C .However, larger N C give rise to a longer computation time due to the increase of the size of the matrix Û in Eq. ( 14).Hence, there is a trade-off between the accuracy and the efficiency.The best value of parameter N C can be estimated as where d is dimensionality of the space, a is the distance between neighboring sites in the lattice, and ℓ is a characteristic length of decay of the matrix element B ij with increasing the distance between sites i and j.

Figure 1 :
Figure 1: Sketch of matrix Û for the one-dimensional case.White and black squares represent zeros and unities, respectively.

Figure 2 :
Figure 2: Numerical example of one-dimensional tight-binding model: (a) on-site energies H jj at different nodes; (b) electron density in the lowest energy band calculated by three methods: exact diagonalization of the Hamiltonian (Eq.(18), blue lines), Method 1 (red dots), and Method 2 (orange circles).

Figure 3 :
Figure 3: The average calculation time for the electron density in the one-dimensional tightbinding model shown in Fig. 2 as a function of the number of nodes L. Different curves correspond to different numerical methods.

Figure 6 :
Figure 6: Comparison between the exact effective potential (solid blue lines) and the filtered potential, (dashed red lines) for a one-dimensional sample with while-noise potential.Reprinted with permission from Nenashev et al., Phys.Rev. B 107, 064207 (2023).Copyright (2023) by the American Physical Society.

Figure 7 :
Figure 7: Comparison between the exact electron density n(x, y, T ) (upper part) and that obtained by using the universal filter function (lower part) in a two-dimensional white-noise potential.Reprinted with permission from Nenashev et al., Phys.Rev. B 107, 064207 (2023).Copyright (2023) by the American Physical Society.

Figure 8 :
Figure 8: Comparison between exact (solid lines) and approximate (symbols) reduced electron density ñ(x, T ) = e −µ/(k B T ) n(x, T ) in a one-dimensional white-noise potential at different temperatures T and N R = 1000 iterations.For clarity, the lowest three curves are scaled down by a multiplication by 10 −3 , 10 −2 and 10 −1 , as indicated in the legend.Reprinted with permission from Nenashev et al., Phys.Rev. B 107, 064207 (2023).Copyright (2023) by the American Physical Society.

Figure 9 :
Figure 9: Reduced electron density ñ(x, y, z, T ) in a three-dimensional sample with whitenoise potential.Compared are the exact density (solid orange line) with that obtained by the RWF algorithm (dashed green line for N R = 20 and dotted blue line for N R = 1000).Profiles along the x-axis with different values of coordinates y and z are shown (see inset).Sample size is 2 × 2 × 2 dimensionless units, the discretization grid parameter is a = 0.1, the temperature is T = T 0 .Periodic boundary conditions apply.For clarity, profiles are multiplied by different coefficients, as indicated in the plot.Reprinted with permission from Nenashev et al., Phys.Rev. B 107, 064207 (2023).Copyright (2023) by the American Physical Society.

) Methods 1
and 2 are based on the fact that this function is close to the Fermi function f (ε) in a vicinity of the Fermi energy ε f , see Eqs. (31), (34) and (35).Let us now consider the behavior of the function f (ε) in a broad range of energies ε.
and then one needs to choose the maximal number N that obeys the restrictions (46),