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

TL;DR

This paper introduces numerical methods to calculate electron distribution in disordered semiconductors, avoiding complex quantum solutions.

Contribution

The paper introduces novel numerical techniques for calculating charge-carrier density in disordered media under different statistical conditions.

Findings

Two numerical approaches for degenerate systems with Fermi statistics were developed and validated.

Two methods for non-degenerate systems with Boltzmann statistics were introduced and tested.

The approximate calculations showed high accuracy when compared to exact quantum-mechanical solutions.

Abstract

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.

Figures

Full text

1. Introduction

Disordered materials, such as amorphous organic and inorganic semiconductors and semiconductor alloys, play an important role in modern optoelectronics for computing, communications, photovoltaics, sensing, and light emission [1,2,3,4,5,6,7,8]. The spatial and energy disorder creates a random potential, that decisively affects electron states. Among other effects, the disorder potential causes spatial localization of electrons in the low-energy range. The knowledge of the space- and temperature-dependent electron distribution $[eqn]$ , particularly in localized states, created by random potentials, is required to understand charge transport and light absorption/emission in disordered semiconductors. The distribution $[eqn]$ is most straightforwardly obtained by solving the Schrödinger equation in the presence of 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 $[eqn]$ without solving the Schrödinger equation.

One of the currently mostly used theoretical tools to reveal the individual features of localized states in a random potential is the so-called localization-landscape theory (LLT) [9,10,11,12]. In the LLT, the random potential is converted into some effective potential, drastically simplifying the calculations. However, recent studies [13,14] revealed substantial problems of the LLT. For instance, the effective potential in the LLT lacks the temperature dependence, which is necessary to describe $[eqn]$ appropriately. Moreover, the LLT has been proven equivalent to the Lorentzian filter applied to a random potential [13]. Such 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 [13]. Therefore, more-developed computational techniques are desirable.

Here, we develop two numerical techniques to reveal $[eqn]$ 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 $[eqn]$ . 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 $[eqn]$ in non-degenerate systems controlled by Boltzmann statistics [13,14]. One algorithm is based on applying a temperature-dependent universal low-pass filter (ULF) to the random potential $[eqn]$ . This yields a temperature-dependent effective potential, $[eqn]$ , that enables a quasiclassical calculation of particle density $[eqn]$ . 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 a recursive application of the Hamiltonian to multiple sets of random wave functions (RWFs) for a specific realization of the random potential $[eqn]$ . 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 [13,14], the RWF scheme does not involve adjustable parameters, and it works across all temperatures. Additionally, the accuracy of the RWF approach in computing $[eqn]$ can systematically be improved, whereas the accuracy of the LLT is inherently limited.

2. Calculation of n(r,T) in a Degenerate System Controlled by Fermi Statistics

To be definite, we consider a disorder potential characterized by Gaussian statistics (‘white noise’), i.e., the potential obeys $[eqn]$ with the auto-correlation function [15]

$$[eqn]$$

where $[eqn]$ indicates the average over many realizations $[eqn]$ of the random potential and S is the strength of the interaction. The quantity S yields natural definitions for the characteristic length scale and for the characteristic energy scale in the form

$$[eqn]$$ $$[eqn]$$

where d is the space dimensionality and $[eqn]$ is the Boltzmann constant.

We consider a collection of non-interacting electrons in some external potential in the thermodynamic equilibrium characterized by the temperature T and the Fermi energy $[eqn]$ . The goal is to develop an effective numerical method to calculate the electron density (concentration) $[eqn]$ as a function of coordinates $[eqn]$ in a degenerate system. The electron density can be defined as

$$[eqn]$$

In this equation, summation index a labels the eigenstates, i.e., solutions of the Schrödonger equation $[eqn]$ , where $[eqn]$ and $[eqn]$ are the wavefunction and the energy of the eigenstate, $[eqn]$ is the Fermi function,

$$[eqn]$$

and factor 2 in front of the sum in Equation (4) accounts for the two possible spin orientations.

We assume that the system under study is discretized with a finite-difference (or tight-binding) method. The tight-binding method is widely used for modeling of various transport phenomena in the presence of disorder, ranging from Anderson localization to disordered topological insulators. As a few examples, let us mention studies of one-dimensional disordered chains [16,17,18], tight-binding models of topological insulators in amorphous systems [19], tight-binding models for lead halide perovskites [20], and nucleic acid sequences [21]. A wavefunction $[eqn]$ is, therefore, represented as a collection of probability amplitudes $[eqn]$ at L points (grid nodes) evenly distributed in space. The Hamiltonian $[eqn]$ is a matrix of size $[eqn]$ . Equation (4) for $[eqn]$ , in this discrete setting, attains the form

$$[eqn]$$

where $[eqn]$ is the electron density at the grid node $[eqn]$ , $[eqn]$ is the value of the eigenfunction $[eqn]$ at node i; $[eqn]$ is the electron energy that corresponds to this eigenfunction; and $[eqn]$ is the spatial volume per one grid node (in the one-dimensional case, $[eqn]$ is simply the distance between nodes).

Equation (6) represents a standard way to calculate numerically 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 $[eqn]$ . In Method 1 (see Section 2.1), the numerical solution of the eigenvalue problem is replaced by the 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 equations that is even faster than the matrix inversion. The performances of Methods 1 and 2, as compared to the standard method based on the eigenvalue problem, are tested in Section 2.3 on the 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

$$[eqn]$$

where a number $[eqn]$ is large, resembles the shape of the Fermi function in the vicinity of point $[eqn]$ . Other approaches for approximating the Fermi function are also possible [22,23,24,25]. The essence of Methods 1 and 2 is the replacement of x with an appropriate linear function of the Hamiltonian. The detailed justification is provided in Appendix A.

It is implicitly assumed in this section that the Fermi energy is inside the energy gap far from the band edges, as the model Hamiltonians exclude two-particle interaction. In a situation, when the two-particle interaction is unavoidable, it can be included into the model via the exchange-correlation potential of the density functional theory. The latter issue is out of scope in our study.

2.1. Method 1: Matrix Inversion

The inputs of Method 1 are the Hamiltonian $[eqn]$ (a matrix $[eqn]$ ), the Fermi energy $[eqn]$ , and two additional parameters, a “reference energy” $[eqn]$ and the number of iterations N. The three latter parameters are related to temperature T via the equality

$$[eqn]$$

The details for the choice of these parameters are given in Appendix B.

In Method 1, one first composes the matrix $[eqn]$ from the Hamiltonian,

$$[eqn]$$

where $[eqn]$ is the $[eqn]$ unit matrix. This matrix $[eqn]$ is then squared N times

$$[eqn]$$

Afterwards, the unit matrix is added to matrix $[eqn]$ and the sum is inverted to yield a new matrix

$$[eqn]$$

Finally, the electron density is obtained from the diagonal elements $[eqn]$ of the matrix $[eqn]$ :

$$[eqn]$$ $$[eqn]$$

2.2. Method 2: Solving a System of Linear Equation

Similarly to Method 1, the input parameters are $[eqn]$ , $[eqn]$ and N, which are related to temperature T by Equation (8). In addition, Method 2 requires one extra parameter $[eqn]$ . The choice of $[eqn]$ is discussed in Appendix B. For a given $[eqn]$ , a matrix $[eqn]$ of size $[eqn]$ is to be composed, that fulfills the following conditions (see Figure 1 for the shape of this matrix in the one-dimensional case):

• each entry of matrix $[eqn]$ is equal to either 0 or 1;
• in each row of matrix $[eqn]$ , exactly one entry is equal to 1;
• in each column of matrix $[eqn]$ , 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 $[eqn]$ zeros, as illustrated in Figure 1.

The matrix $[eqn]$ is calculated, as done in Method 1, see Equations (9) and (10). However, in contrast to Method 1, no matrix inversion is necessary. Instead, a system of linear equations

$$[eqn]$$

is to be solved with respect to the yet unknown matrix $[eqn]$ of size $[eqn]$ . This is a computationally easier task than the matrix inversion used in Method 1.

Finally, the electron density is calculated as

$$[eqn]$$ $$[eqn]$$

2.3. A Numerical Example: One-Dimensional Disordered System with a Single Occupied Band

Let us compare the performances of different methods in calculating $[eqn]$ on a simple one-dimensional tight-binding model with energy bands and disorder. We consider a linear chain of L lattice nodes at the distance $[eqn]$ from each other. The distances and energies are measured in the units determined by Equations (2) and (3). In these units, the hopping integrals between the neighboring nodes $[eqn]$ and $[eqn]$ are equal to $[eqn]$ . The on-site energies $[eqn]$ are chosen to be

$$[eqn]$$

where $[eqn]$ are the random numbers uniformly distributed in the range $[eqn]$ . All the other matrix elements of the Hamiltonian $[eqn]$ are equal to zero. Periodic boundary conditions apply. The constant term $[eqn]$ makes the lower boundary of the energy spectrum $[eqn]$ close to zero, and the higher boundary is $[eqn]$ . The amplitude of the periodic variations of the potential is $[eqn]$ . The amplitude of the random-noise potential $[eqn]$ is chosen to set the ’disorder strength’ S equal to unity. This relation can be understood as follows. In a discretized 1D model, the Dirac delta function $[eqn]$ turns into $[eqn]$ , where $[eqn]$ is a Kronecker delta, and a is the distance between sites. Hence, Equation (1), reformulated with the discrete variables, has the form $[eqn]$ . Inserting potentials $[eqn]$ , $[eqn]$ yields $[eqn]$ . Since $[eqn]$ , one obtains $[eqn]$ for the assumed values $[eqn]$ and $[eqn]$ . As an example, we show one realization of the on-site energies in Figure 2a.

As a consequence of the periodic potential (the second term in the right-hand side of Equation (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

$$[eqn]$$

where summation is performed over the eigenstates of the occupied band.

The wave functions $[eqn]$ , that enter Equation (18), can be obtained by diagonalizing the Hamiltonian. An example of the calculated electron density distribution is shown in Figure 2b by blue lines.

Also we show in Figure 2b the approximated electron density obtained by Method 1 (red dots) and by Method 2 (orange circles). The parameters for these methods are: $[eqn]$ (the middle of the lowest band gap), $[eqn]$ , $[eqn]$ , and $[eqn]$ . One can see that both Methods provide quite accurate results for the electron density.

In Figure 3, we compare the time required by three ways of calculating the electron density—the exact method based on the Hamiltonian diagonalization, Method 1, and Method 2—on a desktop PC with MATLAB (version R2023b, by MathWorks) 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 increase in speed by more than an order of magnitude.

Note that Method 1 can be further improved by using advanced techniques to calculate the diagonal part of the inverted matrix [26,27,28]. 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.

So far, we considered here degenerate systems with Fermi statistics. In Section 3, we address a complementary case of a non-degenerate system with Boltzmann statistics, following two recent studies [13,14]. While in those studies, hundreds of equations were used to treat the non-degenerate case, we demonstrate in this paper that the essence of the theory for Boltzmann statistics can be introduced using only several equations.

3. Calculation of n(r,T) in a Non-Degenerate System Controlled by Boltzmann Statistics

3.1. Low-Pass Filter (LF) Approach

3.1.1. Motivation

Already in the 1960s, Halperin and Lax recognized that electrons in a random potential cannot follow very short-range potential fluctuations [15]. This effect is illustrated in Figure 4, where the random white-noise potential $[eqn]$ in one dimension is depicted by the green solid line. The detailed shape of $[eqn]$ in the region $[eqn]$ (in the dimensionless units given by Equation (2)) is compared with the shape of the wave function $[eqn]$ for the lowest energy state in this spatial region. Apparently, the characteristic width $[eqn]$ of the wave function, even for low-energy localized states, is substantially larger than the spatial scale of the fluctuations of the disorder potential $[eqn]$ . The latter scale in semiconductor alloys is of the order of the lattice constant $[eqn]$ nm. The strong inequality $[eqn]$ suggests that electrons in localized states are affected only by the mean disorder potential, averaged over the space scale $[eqn]$ . 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 $[eqn]$ a low-pass filter [15] (LF) that smooths the spatial fluctuations of $[eqn]$ .

Halperin and Lax suggested the square of the wave function as the filter function [15]. The width of the filter function, $[eqn]$ was adjusted dependent on the energy, $[eqn]$ , of the localized state, $[eqn]$ , where $[eqn]$ 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 [15]. Baranovskii and Efros [29] 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 [15,29]. Our aim here is, in contrast, to calculate the spatial distribution of electron density, $[eqn]$ , in a given disorder potential $[eqn]$ . We start below with the definition of the low-pass filter and then formulate the algorithm to calculate $[eqn]$ . Afterwards, we will formulate an even more accurate technique for calculating $[eqn]$ based on the random-wave-function approach.

3.1.2. Definition of a Low-Pass Filter (LF)

In one dimension, the low-pass filter (LF) is determined by the operation

$$[eqn]$$

where the filter function $[eqn]$ should contain the appropriate length scale ℓ. This operation converts the real disorder potential $[eqn]$ into the smooth effective potential $[eqn]$ . For instance, one can try a Lorentzian function $[eqn]$ ,

$$[eqn]$$

because using as LF a Lorentzian function with $[eqn]$ has recently been proven [13] to be equivalent to the popular LLT approach [9,10,11,12].

Halperin and Lax [15] suggested instead to use for LF the square of the wave function $[eqn]$ . The shape of the wave functions $[eqn]$ for the low-energy states was determined by the optimal-fluctuation-approach that yields the filtering function

$$[eqn]$$

where the characteristic length $[eqn]$ should depend on the state energy [15,29]. Remarkably, it appears that a universal, energy-independent value for $[eqn]$ , can be introduced [13], $[eqn]$ as evidenced in Figure 5a, where the effective potential yielded by the filtering function given by Equation (21), is compared with the positions and energies of the eigenstates. The eigenstates for the given realization of disorder potential $[eqn]$ were obtained via a straightforward solution of the Schrödinger equation. In Figure 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 Equation (21).

In Figure 5b, such a comparison is illustrated for the case of a Lorentzian filter, determined by Equation (20) with $[eqn]$ chosen to mimic the LLT result [13]. Evidently, the choice of a Lorentzian filter function is not satisfying. The number of the local minima in the effective potential $[eqn]$ (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 Equation (20) is not smooth at $[eqn]$ . The cusp at $[eqn]$ filters too many extrema from the real disorder potential $[eqn]$ , preventing the identification of true localized electron states by searching the minima of the effective potential $[eqn]$ . The filter function suggested by Halperin and Lax [15] does not possess such a deficiency.

Not only the energies and spatial positions of localized states in disorder potential $[eqn]$ , 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 $[eqn]$ . Below we extend the LF approach to calculate $[eqn]$ . For that purpose, we introduce the temperature T into the filter function and, concomitantly, into the definition of the effective potential $[eqn]$ that yields $[eqn]$ .

3.1.3. Universal Filter Function to Determine n(r,T)

The T-dependent spatial distribution of electron density $[eqn]$ is related to the quasi-classical effective potential $[eqn]$ as

$$[eqn]$$

where $[eqn]$ is the effective density of states in the conduction band and $[eqn]$ is the chemical potential. This equation serves as the definition of the quasi-classical effective potential $[eqn]$ . The effective potential is smooth in comparison to the initial disorder potential $[eqn]$ because $[eqn]$ is derived from the electron density $[eqn]$ , which has the spatial scale of the electron wave functions, i.e., is broader than the scale of the short-range fluctuations of $[eqn]$ . Let us, therefore, obtain $[eqn]$ by subjecting $[eqn]$ to the action of a universal low-pass filter (ULF).

The key question is how to find out the appropriate T-dependent filter function $[eqn]$ that can be used to extract the shape of the effective potential $[eqn]$ for a given realization of the white-noise potential $[eqn]$ . One can represent this filter function as a functional derivative

$$[eqn]$$

This derivative can be calculated analytically via the first-order perturbation theory in the approximation of vanishing function $[eqn]$ , i.e., for the free electron gas. This perturbation approach yields the following shape of the filter function’s Fourier image [14]:

$$[eqn]$$

where erfi is the imaginary error function, and $[eqn]$ .

In order to reveal the electron density distribution $[eqn]$ for a given realization of the white-noise potential $[eqn]$ , one should first calculate the Fourier image $[eqn]$ of $[eqn]$ using a fast-Fourier-transform (FFT). This function $[eqn]$ should be then multiplied by $[eqn]$ , which is given by Equation (23). The inverse Fourier transform of the product $[eqn]$ by the FFT yields the effective potential $[eqn]$ because the inverse Fourier transform converts a product into a convolution [14]. Inserting $[eqn]$ into Equation (22) gives the electron density $[eqn]$ . In Figure 6 and Figure 7, we compare the results for $[eqn]$ of the above procedure with the effective potentials obtained via Equation (22) from the electron density $[eqn]$ 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 the units given by Equations (2) and (3).

The data in Figure 6 and Figure 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 computation time, in contrast to the exact calculations on the basis of Schrödinger equation.

Below, we introduce an even more accurate technique to calculate electron density $[eqn]$ while avoiding solving the Schrödinger equation. This technique is based on the random-wave-functions approach.

3.2. Random Wave Functions (RWF) Approach to Calculate n(r,T)

3.2.1. Background

The idea of the RWF approach resembles the one suggested recently by Lu and Steinerberger [30] 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 [30]. A similar approach has been suggested by Krajewski and Parrinello [31] for the calculation of the thermodynamic potential.

Let us consider the action of the operator $[eqn]$ 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 $[eqn]$ contribute to the distribution of electrons with the probability proportional to $[eqn]$ . The wave function is the probability amplitude, which explains the factor $[eqn]$ in the operator $[eqn]$ . The wave function is always a linear combination of eigenfunctions that correspond to different energies. The action of the operator $[eqn]$ suppresses the contributions of high-energy eigenfunctions in favor of the contributions of low-energy eigenfunctions. By the application of the operator $[eqn]$ 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 $[eqn]$ . The question arises on how to numerically subject a wave function to the action of the operator $[eqn]$ . This can be done by recursively applying the Hamiltonian $[eqn]$ to the wave function:

$$[eqn]$$

with a natural number [14] $[eqn]$ and a small parameter $[eqn]$ . A simple analysis justifies the choice [14] $[eqn]$ , where $[eqn]$ is the estimate of the upper boundary of the energy distribution. In the case of a regular grid with the lattice constant a, $[eqn]$ . Below, we describe how to realize this idea technically.

3.2.2. The RWF Algorithm

Let us consider the RWF algorithm on a spatial lattice with the volume $[eqn]$ per lattice site. The value of the random wave function $[eqn]$ on each lattice site is chosen independently as a random number extracted from a normal distribution with the average value zero and variance $[eqn]$ . The following transformation of the wave function $[eqn]$ ,

$$[eqn]$$

is applied M times. Then, an estimate of the reduced electron density $[eqn]$ is

$$[eqn]$$

The calculation of $[eqn]$ is carried out for a large number $[eqn]$ of realizations $[eqn]$ of the random wave function $[eqn]$ . Then, the electron density $[eqn]$ is the arithmetic mean of the functions $[eqn]$ obtained for different realizations R, multiplied by a chemical-potential-related factor $[eqn]$ ,

$$[eqn]$$

The larger the number of realizations $[eqn]$ is, the more accurate the calculated electron density $[eqn]$ will be.

In Figure 8, the reduced electron densities obtained via the RWF algorithm are compared with the exact one calculated using the solution of the Schrödinger equation for three dimensions. Evidently, the RWF algorithm with 1000 iterations accurately yields the electron density.

4. Conclusions

In this work, we introduce four theoretical tools to get access to the space- and temperature-dependent electron density $[eqn]$ in disordered media with a random potential $[eqn]$ , 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 $[eqn]$ . 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 $[eqn]$ . 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 $[eqn]$ . 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 when 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.