HTML conversions sometimes display errors due to content that did not convert correctly from the source. This paper uses the following packages that are not yet supported by the HTML conversion tool. Feedback on these issues are not necessary; they are known and are being worked on.

  • failed: chemformula

Authors: achieve the best HTML results from your LaTeX submissions by following these best practices.

License: CC BY 4.0
arXiv:2312.13022v1 [cond-mat.dis-nn] 20 Dec 2023

Computation of the spatial distribution of charge-carrier density in disordered media

Alexey V. Nenashev Faculty of Physics, Philipps-UniversitΓ€t Marburg, Marburg 35032, Germany nenashev_isp@mail.ru    Florian Gebhard Faculty of Physics, Philipps-UniversitΓ€t Marburg, Marburg 35032, Germany    Klaus Meerholz Department fΓΌr Chemie, UniversitΓ€t zu KΓΆln, Greinstraße 4-6, 50939 KΓΆln, Germany    Sergei D. Baranovskii Department fΓΌr Chemie, UniversitΓ€t zu KΓΆln, Greinstraße 4-6, 50939 KΓΆln, Germany sergei.baranovski@physik.uni-marburg.de
Abstract

The space- and temperature-dependent electron distribution n⁒(r,T)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_T ) determines optoelectronic properties of disordered semiconductors. It is a challenging task to get access to n⁒(r,T)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_T ) in random potentials, 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 that 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.

[Uncaptioned image]
keywords:
Disordered materials, electron states in random potential
\SectionNumbersOn\altaffiliation

On leave of absence from Rzhanov Institute of Semiconductor Physics and the Novosibirsk State University, Russia \alsoaffiliationFaculty of Physics, Philipps-UniversitΓ€t Marburg, Marburg 35032, Germany \abbreviationsIR,NMR,UV

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. 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)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_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)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_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)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_T ) 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) 4, 5, 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)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_T ) appropriately. Moreover, the LLT is seen to be equivalent to the Lorentzian filter applied to a random potential 8. 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 8. Therefore, more developed computational techniques are desirable.

Here we develop two numerical techniques to reveal n⁒(r,T)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_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)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_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)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_T ) in non-degenerate systems controlled by Boltzmann statistics 8, 9. One algorithm is based on applying a temperature-dependent universal low-pass filter (ULF) to the random potential V⁒(r)π‘‰π‘ŸV(r)italic_V ( italic_r ). This yields a temperature-dependent effective potential, W⁒(r,T)π‘Šπ‘Ÿπ‘‡W(r,T)italic_W ( italic_r , italic_T ), that enables a quasiclassical calculation of particle density n⁒(r,T)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_T ). The ULF algorithm employs Fast-Fourier Transformation for calculating the effective potential Wπ‘ŠWitalic_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)π‘‰π‘ŸV(r)italic_V ( italic_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)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_T ) can systematically be improved, whereas the accuracy of the LLT is inherently limited.

2 Calculation of n⁒(r,T)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_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=0subscriptdelimited-βŸ¨βŸ©π‘‰π«R0\langle V(\mathbf{r}\,)\rangle_{\rm R}=0⟨ italic_V ( bold_r ) ⟩ start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT = 0 with the auto-correlation function 10

⟨V(𝐫)V(𝐫)β€²βŸ©R=SΞ΄(π«βˆ’π«)β€²,\langle V(\mathbf{r}\,)V(\mathbf{r}\,{}^{\prime})\rangle_{\rm R}=S\delta\left(% \mathbf{r}-\mathbf{r}\,{}^{\prime}\right)\;,⟨ italic_V ( bold_r ) italic_V ( bold_r start_FLOATSUPERSCRIPT β€² end_FLOATSUPERSCRIPT ) ⟩ start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT = italic_S italic_Ξ΄ ( bold_r - bold_r start_FLOATSUPERSCRIPT β€² end_FLOATSUPERSCRIPT ) , (1)

where βŸ¨β€¦βŸ©Rsubscriptdelimited-βŸ¨βŸ©β€¦R\langle\ldots\rangle_{\rm R}⟨ … ⟩ start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT indicates the average over many realizations RR{\rm R}roman_R of the random potential and S𝑆Sitalic_S is the strength of the interaction. This quantity S𝑆Sitalic_S yields natural definitions for the characteristic length scale and for the characteristic energy scale in the form

β„“0=(ℏ4m2⁒S)1/(4βˆ’d),subscriptβ„“0superscriptsuperscriptPlanck-constant-over-2-pi4superscriptπ‘š2𝑆14𝑑\ell_{0}=\left(\frac{\hbar^{4}}{m^{2}S}\right)^{1/(4-d)}\;,roman_β„“ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = ( divide start_ARG roman_ℏ start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT end_ARG start_ARG italic_m start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_S end_ARG ) start_POSTSUPERSCRIPT 1 / ( 4 - italic_d ) end_POSTSUPERSCRIPT , (2)
T0=1kB⁒(md⁒S2ℏ2⁒d)1/(4βˆ’d),subscript𝑇01subscriptπ‘˜π΅superscriptsuperscriptπ‘šπ‘‘superscript𝑆2superscriptPlanck-constant-over-2-pi2𝑑14𝑑T_{0}=\frac{1}{k_{B}}\left(\frac{m^{d}S^{2}}{\hbar^{2d}}\right)^{1/(4-d)}\;,italic_T start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT end_ARG ( divide start_ARG italic_m start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT italic_S start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG roman_ℏ start_POSTSUPERSCRIPT 2 italic_d end_POSTSUPERSCRIPT end_ARG ) start_POSTSUPERSCRIPT 1 / ( 4 - italic_d ) end_POSTSUPERSCRIPT , (3)

where d𝑑ditalic_d is the space dimensionality and kBsubscriptπ‘˜π΅k_{B}italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT 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𝑇Titalic_T and the Fermi energy Ξ΅fsubscriptπœ€π‘“\varepsilon_{f}italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT. The goal is to develop an effective numerical method for the calculation of the electron density (concentration) n⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_T ) as a function of coordinates 𝐫𝐫\mathbf{r}bold_r in a degenerate system. The electron density can be defined as

n⁒(𝐫,T)=2β’βˆ‘a|ψa⁒(𝐫)|2⁒f⁒(Ξ΅a).𝑛𝐫𝑇2subscriptπ‘Žsuperscriptsubscriptπœ“π‘Žπ«2𝑓subscriptπœ€π‘Žn(\mathbf{r},T)=2\sum_{a}|\psi_{a}(\mathbf{r})|^{2}f(\varepsilon_{a}).italic_n ( bold_r , italic_T ) = 2 βˆ‘ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT | italic_ψ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT ( bold_r ) | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_f ( italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT ) . (4)

In this equation, summation index aπ‘Žaitalic_a labels the electron eigenstates, i. e., solutions of the SchrΓΆdonger equation H^⁒ψa=Ξ΅a⁒ψa^𝐻subscriptπœ“π‘Žsubscriptπœ€π‘Žsubscriptπœ“π‘Ž\hat{H}\psi_{a}=\varepsilon_{a}\psi_{a}over^ start_ARG italic_H end_ARG italic_ψ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT = italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT italic_ψ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT, where ψa⁒(𝐫)subscriptπœ“π‘Žπ«\psi_{a}(\mathbf{r})italic_ψ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT ( bold_r ) and Ξ΅asubscriptπœ€π‘Ž\varepsilon_{a}italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT are the wavefunction and the energy of the eigenstate; f⁒(Ξ΅)π‘“πœ€f(\varepsilon)italic_f ( italic_Ξ΅ ) is the Fermi function,

f⁒(Ξ΅)=1exp⁑[(Ξ΅βˆ’Ξ΅f)/kB⁒T]+1,π‘“πœ€1πœ€subscriptπœ€π‘“subscriptπ‘˜π΅π‘‡1f(\varepsilon)=\frac{1}{\exp[(\varepsilon-\varepsilon_{f})/k_{B}T]+1}\,,italic_f ( italic_Ξ΅ ) = divide start_ARG 1 end_ARG start_ARG roman_exp [ ( italic_Ξ΅ - italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT ) / italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_T ] + 1 end_ARG , (5)

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 ψ⁒(𝐫)πœ“π«\psi(\mathbf{r})italic_ψ ( bold_r ) is therefore represented as a collection of probability amplitudes ψ1,…,ψLsubscriptπœ“1…subscriptπœ“πΏ\psi_{1},...,\psi_{L}italic_ψ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , … , italic_ψ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT at L𝐿Litalic_L points (grid nodes) evenly distributed in space. The Hamiltonian H^^𝐻\hat{H}over^ start_ARG italic_H end_ARG is a matrix of size LΓ—L𝐿𝐿L\times Litalic_L Γ— italic_L. Equation (4) for n⁒(r,T)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_T ) in this discrete setting obtains the form

ni=2Δ⁒Vβ’βˆ‘a=1L|ψa,i|2⁒f⁒(Ξ΅a),subscript𝑛𝑖2Δ𝑉superscriptsubscriptπ‘Ž1𝐿superscriptsubscriptπœ“π‘Žπ‘–2𝑓subscriptπœ€π‘Žn_{i}=\frac{2}{\Delta V}\sum_{a=1}^{L}|\psi_{a,i}|^{2}f(\varepsilon_{a}),italic_n start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = divide start_ARG 2 end_ARG start_ARG roman_Ξ” italic_V end_ARG βˆ‘ start_POSTSUBSCRIPT italic_a = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_L end_POSTSUPERSCRIPT | italic_ψ start_POSTSUBSCRIPT italic_a , italic_i end_POSTSUBSCRIPT | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_f ( italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT ) , (6)

where nisubscript𝑛𝑖n_{i}italic_n start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT is the electron density at grid node i∈{1,…,L}𝑖1…𝐿i\in\{1,...,L\}italic_i ∈ { 1 , … , italic_L }; ψa,isubscriptπœ“π‘Žπ‘–\psi_{a,i}italic_ψ start_POSTSUBSCRIPT italic_a , italic_i end_POSTSUBSCRIPT is the value of eigenfunction ψasubscriptπœ“π‘Ž\psi_{a}italic_ψ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT at node i𝑖iitalic_i; Ξ΅asubscriptπœ€π‘Ž\varepsilon_{a}italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT is the electron energy that corresponds this eigenfunction; and Δ⁒VΔ𝑉\Delta Vroman_Ξ” italic_V is the spatial volume per one grid node (in the one-dimensional case, Δ⁒VΔ𝑉\Delta Vroman_Ξ” italic_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𝐿Litalic_L of the Hamiltonian matrix. Below we suggest two methods that allow one to speed up the calculation of n⁒(r,T)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_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

y⁒(x)=1x𝒩+1,𝑦π‘₯1superscriptπ‘₯𝒩1y(x)=\frac{1}{x^{\mathcal{N}}+1}\,,italic_y ( italic_x ) = divide start_ARG 1 end_ARG start_ARG italic_x start_POSTSUPERSCRIPT caligraphic_N end_POSTSUPERSCRIPT + 1 end_ARG , (7)

where a number 𝒩𝒩\mathcal{N}caligraphic_N is large, resembles the shape of the Fermi function in the vicinity of point x=1π‘₯1x=1italic_x = 1. Other approaches for approximating the Fermi function are also possible. 11, 12, 13, 14 Replacing xπ‘₯xitalic_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.

2.1 Method 1: matrix inversion

The method is based on the Hamiltonian H^^𝐻\hat{H}over^ start_ARG italic_H end_ARG (a matrix LΓ—L𝐿𝐿L\times Litalic_L Γ— italic_L), the Fermi energy Ξ΅fsubscriptπœ€π‘“\varepsilon_{f}italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT, and two additional parameters: a β€œreference energy” Ξ΅0subscriptπœ€0\varepsilon_{0}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT and the number of iterations N𝑁Nitalic_N. These parameters are related to the temperature T𝑇Titalic_T,

kB⁒T=|Ξ΅fβˆ’Ξ΅0|2N.subscriptπ‘˜π΅π‘‡subscriptπœ€π‘“subscriptπœ€0superscript2𝑁k_{B}T=\frac{|\varepsilon_{f}-\varepsilon_{0}|}{2^{N}}\,.italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_T = divide start_ARG | italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT | end_ARG start_ARG 2 start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT end_ARG . (8)

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

The method starts from composing the matrix A^0subscript^𝐴0\hat{A}_{0}over^ start_ARG italic_A end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT from the Hamiltonian

A^0=H^βˆ’Ξ΅0⁒I^Ξ΅fβˆ’Ξ΅0,subscript^𝐴0^𝐻subscriptπœ€0^𝐼subscriptπœ€π‘“subscriptπœ€0\hat{A}_{0}=\frac{\hat{H}-\varepsilon_{0}\hat{I}}{\varepsilon_{f}-\varepsilon_% {0}}\,,over^ start_ARG italic_A end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = divide start_ARG over^ start_ARG italic_H end_ARG - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT over^ start_ARG italic_I end_ARG end_ARG start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG , (9)

where I^^𝐼\hat{I}over^ start_ARG italic_I end_ARG is the LΓ—L𝐿𝐿L\times Litalic_L Γ— italic_L unit matrix. Then, the matrix A^0subscript^𝐴0\hat{A}_{0}over^ start_ARG italic_A end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT is squared N𝑁Nitalic_N times

A^p=(A^pβˆ’1)2,p=1,2,…,N.formulae-sequencesubscript^𝐴𝑝superscriptsubscript^𝐴𝑝12𝑝12…𝑁\hat{A}_{p}=\left(\hat{A}_{p-1}\right)^{2},\quad p=1,2,\ldots,N.over^ start_ARG italic_A end_ARG start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT = ( over^ start_ARG italic_A end_ARG start_POSTSUBSCRIPT italic_p - 1 end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT , italic_p = 1 , 2 , … , italic_N . (10)

To the resulting matrix A^Nsubscript^𝐴𝑁\hat{A}_{N}over^ start_ARG italic_A end_ARG start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT we add the unit matrix and invert the sum to get a new matrix

B^=(A^N+I^)βˆ’1.^𝐡superscriptsubscript^𝐴𝑁^𝐼1\hat{B}=\left(\hat{A}_{N}+\hat{I}\right)^{-1}.over^ start_ARG italic_B end_ARG = ( over^ start_ARG italic_A end_ARG start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT + over^ start_ARG italic_I end_ARG ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT . (11)

Finally, the electron density is obtained from the diagonal elements Bi⁒isubscript𝐡𝑖𝑖B_{ii}italic_B start_POSTSUBSCRIPT italic_i italic_i end_POSTSUBSCRIPT of the matrix B^^𝐡\hat{B}over^ start_ARG italic_B end_ARG:

ni=2Δ⁒V⁒Bi⁒iifΞ΅0<Ξ΅f,formulae-sequencesubscript𝑛𝑖2Δ𝑉subscript𝐡𝑖𝑖ifsubscriptπœ€0subscriptπœ€π‘“n_{i}=\frac{2}{\Delta V}\,B_{ii}\quad\mathrm{if}\quad\varepsilon_{0}<% \varepsilon_{f}\,,italic_n start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = divide start_ARG 2 end_ARG start_ARG roman_Ξ” italic_V end_ARG italic_B start_POSTSUBSCRIPT italic_i italic_i end_POSTSUBSCRIPT roman_if italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT < italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT , (12)
ni=2Δ⁒V⁒(1βˆ’Bi⁒i)ifΞ΅0>Ξ΅f.formulae-sequencesubscript𝑛𝑖2Δ𝑉1subscript𝐡𝑖𝑖ifsubscriptπœ€0subscriptπœ€π‘“n_{i}=\frac{2}{\Delta V}\,(1-B_{ii})\quad\mathrm{if}\quad\varepsilon_{0}>% \varepsilon_{f}\,.italic_n start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = divide start_ARG 2 end_ARG start_ARG roman_Ξ” italic_V end_ARG ( 1 - italic_B start_POSTSUBSCRIPT italic_i italic_i end_POSTSUBSCRIPT ) roman_if italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT > italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT . (13)

2.2 Method 2: solving a system of linear equation

Similarly to Method 1, the input parameters are Ξ΅fsubscriptπœ€π‘“\varepsilon_{f}italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT, Ξ΅0subscriptπœ€0\varepsilon_{0}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT and N𝑁Nitalic_N, which are related to the temperature T𝑇Titalic_T by Eq. (8). In addition, this method requires one more parameter NCsubscript𝑁𝐢N_{C}italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT. The choice of NCsubscript𝑁𝐢N_{C}italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT is discussed in Appendix B. For a given NCsubscript𝑁𝐢N_{C}italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT, we compose a matrix U^^π‘ˆ\hat{U}over^ start_ARG italic_U end_ARG of size LΓ—NC𝐿subscript𝑁𝐢L\times N_{C}italic_L Γ— italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT that obeys the following conditions (see Fig. 1 for the shape of this matrix in the one-dimensional case):

  • β€’

    each entry of matrix U^^π‘ˆ\hat{U}over^ start_ARG italic_U end_ARG is equal to either 0 or 1;

  • β€’

    in each row of matrix U^^π‘ˆ\hat{U}over^ start_ARG italic_U end_ARG, exactly one entry is equal to 1;

  • β€’

    in each column of matrix U^^π‘ˆ\hat{U}over^ start_ARG italic_U end_ARG, 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 (NCβˆ’1)subscript𝑁𝐢1(N_{C}-1)( italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT - 1 ) zeros, as illustrated in Fig. 1.

The matrix A^Nsubscript^𝐴𝑁\hat{A}_{N}over^ start_ARG italic_A end_ARG start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 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

(A^N+I^)⁒X^=U^subscript^𝐴𝑁^𝐼^𝑋^π‘ˆ\left(\hat{A}_{N}+\hat{I}\right)\hat{X}=\hat{U}( over^ start_ARG italic_A end_ARG start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT + over^ start_ARG italic_I end_ARG ) over^ start_ARG italic_X end_ARG = over^ start_ARG italic_U end_ARG (14)

is to be solved with respect to an unknown matrix X^^𝑋\hat{X}over^ start_ARG italic_X end_ARG of size LΓ—NC𝐿subscript𝑁𝐢L\times N_{C}italic_L Γ— italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT. This is a computationally easier task in comparison to the matrix inversion.

Refer to caption
Figure 1: Sketch of matrix U^^π‘ˆ\hat{U}over^ start_ARG italic_U end_ARG for the one-dimensional case. White and black squares represent zeros and unities, respectively.

Finally, the electron density is calculated as

ni=2Δ⁒Vβ’βˆ‘a=1NCXi⁒a⁒Ui⁒aifΞ΅0<Ξ΅f,formulae-sequencesubscript𝑛𝑖2Δ𝑉superscriptsubscriptπ‘Ž1subscript𝑁𝐢subscriptπ‘‹π‘–π‘Žsubscriptπ‘ˆπ‘–π‘Žifsubscriptπœ€0subscriptπœ€π‘“n_{i}=\frac{2}{\Delta V}\,\sum_{a=1}^{N_{C}}X_{ia}U_{ia}\quad\mathrm{if}\quad% \varepsilon_{0}<\varepsilon_{f}\,,italic_n start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = divide start_ARG 2 end_ARG start_ARG roman_Ξ” italic_V end_ARG βˆ‘ start_POSTSUBSCRIPT italic_a = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT end_POSTSUPERSCRIPT italic_X start_POSTSUBSCRIPT italic_i italic_a end_POSTSUBSCRIPT italic_U start_POSTSUBSCRIPT italic_i italic_a end_POSTSUBSCRIPT roman_if italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT < italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT , (15)
ni=2Δ⁒V⁒(1βˆ’βˆ‘a=1NCXi⁒a⁒Ui⁒a)ifΞ΅0>Ξ΅f.formulae-sequencesubscript𝑛𝑖2Δ𝑉1superscriptsubscriptπ‘Ž1subscript𝑁𝐢subscriptπ‘‹π‘–π‘Žsubscriptπ‘ˆπ‘–π‘Žifsubscriptπœ€0subscriptπœ€π‘“n_{i}=\frac{2}{\Delta V}\,\left(1-\sum_{a=1}^{N_{C}}X_{ia}U_{ia}\right)\quad% \mathrm{if}\quad\varepsilon_{0}>\varepsilon_{f}\,.italic_n start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = divide start_ARG 2 end_ARG start_ARG roman_Ξ” italic_V end_ARG ( 1 - βˆ‘ start_POSTSUBSCRIPT italic_a = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT end_POSTSUPERSCRIPT italic_X start_POSTSUBSCRIPT italic_i italic_a end_POSTSUBSCRIPT italic_U start_POSTSUBSCRIPT italic_i italic_a end_POSTSUBSCRIPT ) roman_if italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT > italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT . (16)

2.3 Numerical example: a one-dimensional disordered system with one occupied band

Let us compare the performance of different methods to calculate n⁒(r,T)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_T ) in a simple one-dimensional tight-binding model with energy bands and disorder. We consider a linear chain of L𝐿Litalic_L lattice nodes at a distance a=0.1π‘Ž0.1a=0.1italic_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 Hi,i+1subscript𝐻𝑖𝑖1H_{i,i+1}italic_H start_POSTSUBSCRIPT italic_i , italic_i + 1 end_POSTSUBSCRIPT and Hi+1,isubscript𝐻𝑖1𝑖H_{i+1,i}italic_H start_POSTSUBSCRIPT italic_i + 1 , italic_i end_POSTSUBSCRIPT are equal to βˆ’1/(2⁒a2)=βˆ’5012superscriptπ‘Ž250-1/(2a^{2})=-50- 1 / ( 2 italic_a start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) = - 50. The on-site energies Hj⁒jsubscript𝐻𝑗𝑗H_{jj}italic_H start_POSTSUBSCRIPT italic_j italic_j end_POSTSUBSCRIPT are chosen to be

Hj⁒j=V0+V1⁒cos⁑(π⁒j/2)+V2⁒ξj,subscript𝐻𝑗𝑗subscript𝑉0subscript𝑉1πœ‹π‘—2subscript𝑉2subscriptπœ‰π‘—H_{jj}=V_{0}+V_{1}\cos(\pi j/2)+V_{2}\xi_{j}\,,italic_H start_POSTSUBSCRIPT italic_j italic_j end_POSTSUBSCRIPT = italic_V start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + italic_V start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT roman_cos ( italic_Ο€ italic_j / 2 ) + italic_V start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT italic_ΞΎ start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT , (17)

where ΞΎjsubscriptπœ‰π‘—\xi_{j}italic_ΞΎ start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT are random numbers uniformly distributed in the range βˆ’1<ΞΎj<11subscriptπœ‰π‘—1-1<\xi_{j}<1- 1 < italic_ΞΎ start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT < 1. All other matrix elements of the Hamiltonian H^^𝐻\hat{H}over^ start_ARG italic_H end_ARG are equal to zero. Periodic boundary conditions apply. The constant term V0=100subscript𝑉0100V_{0}=100italic_V start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 100 makes the lower boundary of the energy spectrum Ξ΅minsubscriptπœ€min\varepsilon_{\mathrm{min}}italic_Ξ΅ start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT to be near zero, and the higher boundary Ξ΅maxsubscriptπœ€max\varepsilon_{\mathrm{max}}italic_Ξ΅ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT to be β‰ˆ200absent200\approx 200β‰ˆ 200. The amplitude of periodic variations of the potential is V1=20subscript𝑉120V_{1}=20italic_V start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = 20. The amplitude of the random-noise potential V2=30subscript𝑉230V_{2}=\sqrt{30}italic_V start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = square-root start_ARG 30 end_ARG is chosen such that the β€œdisorder strength” S𝑆Sitalic_S is equal to unity. As an example, we show one realization of the on-site energies in Fig. 2a.

Refer to caption
Figure 2: Numerical example of one-dimensional tight-binding model: (a) on-site energies Hj⁒jsubscript𝐻𝑗𝑗H_{jj}italic_H start_POSTSUBSCRIPT italic_j italic_j end_POSTSUBSCRIPT 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).

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

ni=2Δ⁒Vβ’βˆ‘a∈OB|ψa,i|2,subscript𝑛𝑖2Δ𝑉subscriptπ‘ŽOBsuperscriptsubscriptπœ“π‘Žπ‘–2n_{i}=\frac{2}{\Delta V}\sum_{a\in\mathrm{OB}}|\psi_{a,i}|^{2},italic_n start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = divide start_ARG 2 end_ARG start_ARG roman_Ξ” italic_V end_ARG βˆ‘ start_POSTSUBSCRIPT italic_a ∈ roman_OB end_POSTSUBSCRIPT | italic_ψ start_POSTSUBSCRIPT italic_a , italic_i end_POSTSUBSCRIPT | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT , (18)

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

The wave functions ψasubscriptπœ“π‘Ž\psi_{a}italic_ψ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT, that enter Eq. (18), can be obtained by diagonalizing the Hamiltonian. 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.5subscriptπœ€π‘“28.5\varepsilon_{f}=28.5italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT = 28.5 (the middle of the lowest band gap), Ξ΅0=10subscriptπœ€010\varepsilon_{0}=10italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 10, N=3𝑁3N=3italic_N = 3, and NC=30subscript𝑁𝐢30N_{C}=30italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT = 30. One can see that both Methods provide quite accurate results for the electron density.

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

In Fig. 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 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.

Note that Method 1 can be further improved by using advanced techniques to calculate the diagonal part of the inverted matrix 15, 16, 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)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_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. 10 This effect is illustrated in Fig. 4, where the random white-noise potential V⁒(x)𝑉π‘₯V(x)italic_V ( italic_x ) in one dimension is depicted by the green solid line. The detailed shape of V⁒(x)𝑉π‘₯V(x)italic_V ( italic_x ) in the region 275≀x≀325275π‘₯325275\leq x\leq 325275 ≀ italic_x ≀ 325 (in the dimensionless units given by Eq. (2)) is compared with the shape of the wave function Οˆπœ“\psiitalic_ψ for the lowest energy state in this spatial region. Apparently, the characteristic width β„“w⁒fsubscriptℓ𝑀𝑓\ell_{wf}roman_β„“ start_POSTSUBSCRIPT italic_w italic_f end_POSTSUBSCRIPT 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)π‘‰π‘ŸV(r)italic_V ( italic_r ). The latter scale in semiconductor alloys is of the order of the lattice constant aβ‰ˆ0.5π‘Ž0.5a\approx 0.5italic_a β‰ˆ 0.5 nm. The strong inequality β„“w⁒f≫amuch-greater-thansubscriptβ„“π‘€π‘“π‘Ž\ell_{wf}\gg aroman_β„“ start_POSTSUBSCRIPT italic_w italic_f end_POSTSUBSCRIPT ≫ italic_a suggests that electrons in localized states are affected only by the mean disorder potential, averaged over the space scale β„“w⁒fsubscriptℓ𝑀𝑓\ell_{wf}roman_β„“ start_POSTSUBSCRIPT italic_w italic_f end_POSTSUBSCRIPT. 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)π‘‰π‘ŸV(r)italic_V ( italic_r ) a low-pass filter 10 (LF) that smoothes the spatial fluctuations of V⁒(r)π‘‰π‘ŸV(r)italic_V ( italic_r ).

Halperin and Lax suggested the square of the wave function as the filter function 10. The width of the filter function, β„“β‰ˆβ„“w⁒fβ„“subscriptℓ𝑀𝑓\ell\approx\ell_{wf}roman_β„“ β‰ˆ roman_β„“ start_POSTSUBSCRIPT italic_w italic_f end_POSTSUBSCRIPT was adjusted dependent on the energy, Ξ΅πœ€\varepsilonitalic_Ξ΅, of the localized state, β„“βˆ|Ξ΅|βˆ’1/2proportional-toβ„“superscriptπœ€12\ell\propto|\varepsilon|^{-1/2}roman_β„“ ∝ | italic_Ξ΅ | start_POSTSUPERSCRIPT - 1 / 2 end_POSTSUPERSCRIPT, where Ξ΅πœ€\varepsilonitalic_Ξ΅ 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 10. Baranovskii 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, 18. Our aim here is, in contrast, to calculate the spatial distribution of electron density, n⁒(r,T)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_T ), in a given disorder potential V⁒(r)π‘‰π‘ŸV(r)italic_V ( italic_r ). We start below with the definition of the low-pass filter and then formulate the algorithm for the calculation of n⁒(r,T)π‘›π‘Ÿπ‘‡n(r,T)italic_n ( italic_r , italic_T ).

Refer to caption
Figure 4: Realization for the white-noise disorder potential on a one-dimensional strip. Insert: the wave function ψ⁒(x)πœ“π‘₯\psi(x)italic_ψ ( italic_x ) of the state with the lowest energy in the region 275≀x≀325275π‘₯325275\leq x\leq 325275 ≀ italic_x ≀ 325. The coordinate xπ‘₯xitalic_x is dimensionless in the units given by Eq. (2).

3.1.2 Definiton of a low-pass filter (LF)

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

W⁒(x)=∫dx′⁒Γ⁒(xβˆ’xβ€²)⁒V⁒(xβ€²),π‘Šπ‘₯differential-dsuperscriptπ‘₯β€²Ξ“π‘₯superscriptπ‘₯′𝑉superscriptπ‘₯β€²W(x)=\int{\rm d}x^{\prime}\Gamma(x-x^{\prime})V(x^{\prime})\,,italic_W ( italic_x ) = ∫ roman_d italic_x start_POSTSUPERSCRIPT β€² end_POSTSUPERSCRIPT roman_Ξ“ ( italic_x - italic_x start_POSTSUPERSCRIPT β€² end_POSTSUPERSCRIPT ) italic_V ( italic_x start_POSTSUPERSCRIPT β€² end_POSTSUPERSCRIPT ) , (19)

where the filter function ΓΓ\Gammaroman_Ξ“ should contain the appropriate length scale β„“β„“\ellroman_β„“. This operation converts the real disorder potential V⁒(x)𝑉π‘₯V(x)italic_V ( italic_x ) into the smooth effective potential W⁒(x)π‘Šπ‘₯W(x)italic_W ( italic_x ). For instance, one can try a Lorentzian function Ξ“LsuperscriptΞ“L\Gamma^{\rm L}roman_Ξ“ start_POSTSUPERSCRIPT roman_L end_POSTSUPERSCRIPT,

Ξ“L⁒(x)=eβˆ’|x|/β„“L2⁒ℓL,superscriptΞ“Lπ‘₯superscript𝑒π‘₯subscriptβ„“L2subscriptβ„“L\Gamma^{\rm L}(x)=\frac{e^{-|x|/\ell_{\rm L}}}{2\ell_{\rm L}}\;,roman_Ξ“ start_POSTSUPERSCRIPT roman_L end_POSTSUPERSCRIPT ( italic_x ) = divide start_ARG italic_e start_POSTSUPERSCRIPT - | italic_x | / roman_β„“ start_POSTSUBSCRIPT roman_L end_POSTSUBSCRIPT end_POSTSUPERSCRIPT end_ARG start_ARG 2 roman_β„“ start_POSTSUBSCRIPT roman_L end_POSTSUBSCRIPT end_ARG , (20)

because using as LF a Lorentzian function with β„“L=0.27⁒ℓ0subscriptβ„“L0.27subscriptβ„“0\ell_{\rm{L}}=0.27\,\ell_{0}roman_β„“ start_POSTSUBSCRIPT roman_L end_POSTSUBSCRIPT = 0.27 roman_β„“ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT has recently been proven 8 equivalent to the popular LLT approach 4, 5, 6, 7.

Halperin and Lax 10 suggested instead to use for LF the square of the wave function Γ⁒(x)=ψ2⁒(x)Ξ“π‘₯superscriptπœ“2π‘₯\Gamma(x)=\psi^{2}(x)roman_Ξ“ ( italic_x ) = italic_ψ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_x ). The shape of the wave functions ψ⁒(x)πœ“π‘₯\psi(x)italic_ψ ( italic_x ) for the low-energy states was determined by the optimal-fluctuation-approach that yields the filtering function

Ξ“HL⁒(x)=12⁒ℓHL⁒1cosh2⁑(x/β„“HL),superscriptΞ“HLπ‘₯12subscriptβ„“HL1superscript2π‘₯subscriptβ„“HL\Gamma^{\rm HL}(x)=\frac{1}{2\ell_{\rm HL}}\frac{1}{\cosh^{2}(x/\ell_{\rm HL})% }\;,roman_Ξ“ start_POSTSUPERSCRIPT roman_HL end_POSTSUPERSCRIPT ( italic_x ) = divide start_ARG 1 end_ARG start_ARG 2 roman_β„“ start_POSTSUBSCRIPT roman_HL end_POSTSUBSCRIPT end_ARG divide start_ARG 1 end_ARG start_ARG roman_cosh start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_x / roman_β„“ start_POSTSUBSCRIPT roman_HL end_POSTSUBSCRIPT ) end_ARG , (21)

where the characteristic length β„“HLsubscriptβ„“HL\ell_{\rm HL}roman_β„“ start_POSTSUBSCRIPT roman_HL end_POSTSUBSCRIPT should depend on the state energy 10, 18. Remarkably, it appears that a universal, energy-independent value for β„“HLsubscriptβ„“HL\ell_{\rm HL}roman_β„“ start_POSTSUBSCRIPT roman_HL end_POSTSUBSCRIPT, can be introduced 8, β„“HL=0.76⁒ℓ0subscriptβ„“HL0.76subscriptβ„“0\ell_{\rm HL}=0.76\,\ell_{0}roman_β„“ start_POSTSUBSCRIPT roman_HL end_POSTSUBSCRIPT = 0.76 roman_β„“ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT 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)𝑉π‘₯V(x)italic_V ( italic_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).

Refer to caption
Figure 5: (a) Effective potential for a Halperin-Lax low-pass filter. (b) Effective potential for a Lorentzian low-pass filter. Reprinted with permission from Gebhard et al., Phys. Rev. B 107, 064206 (2023). Copyright (2023) by the American Physical Society.

In Fig. 5b, such a comparison is illustrated for the case of a Lorentzian filter, determined by Eq. (20) with β„“L=0.27⁒ℓ0subscriptβ„“L0.27subscriptβ„“0\ell_{\rm{L}}=0.27\,\ell_{0}roman_β„“ start_POSTSUBSCRIPT roman_L end_POSTSUBSCRIPT = 0.27 roman_β„“ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT chosen to mimic the LLT result 8. Evidently, the choice of a Lorentzian filter function is not satisfying. The number of the local minima in the effective potential W⁒(x)π‘Šπ‘₯W(x)italic_W ( italic_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π‘₯0x=0italic_x = 0. The cusp at x=0π‘₯0x=0italic_x = 0 filters too many extrema from the real disorder potential V⁒(x)𝑉π‘₯V(x)italic_V ( italic_x ), preventing the identification of true localized electron states by searching the minima of the effective potential W⁒(x)π‘Šπ‘₯W(x)italic_W ( italic_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⁒(𝐫)𝑉𝐫V(\mathbf{r})italic_V ( bold_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⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_T ). Below we extend the LF approach to calculate n⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_T ). For that purpose, we introduce the temperature T𝑇Titalic_T into the filter function and, concomitantly, into the definition of the effective potential W⁒(𝐫,T)π‘Šπ«π‘‡W(\mathbf{r},T)italic_W ( bold_r , italic_T ) that yields n⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_T ).

3.1.3 Universal filter function to determine n⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_T )

The T𝑇Titalic_T-dependent spatial distribution of electron density n⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_T ) is related to the quasi-classical effective potential W⁒(𝐫,T)π‘Šπ«π‘‡W(\mathbf{r},T)italic_W ( bold_r , italic_T ) as

n⁒(𝐫,T)=Nc⁒exp⁑[ΞΌβˆ’W⁒(𝐫,T)kB⁒T],𝑛𝐫𝑇subscriptπ‘π‘πœ‡π‘Šπ«π‘‡subscriptπ‘˜B𝑇n(\mathbf{r},T)=N_{c}\exp\left[\frac{\mu-W(\mathbf{r},T)}{k_{\rm B}T}\right]\;,italic_n ( bold_r , italic_T ) = italic_N start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT roman_exp [ divide start_ARG italic_ΞΌ - italic_W ( bold_r , italic_T ) end_ARG start_ARG italic_k start_POSTSUBSCRIPT roman_B end_POSTSUBSCRIPT italic_T end_ARG ] , (22)

where Ncsubscript𝑁𝑐N_{c}italic_N start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT is the effective density of states in the conduction band and ΞΌπœ‡\muitalic_ΞΌ is the chemical potential. This equation serves as the definition of the quasi-classical effective potential W⁒(𝐫,T)π‘Šπ«π‘‡W(\mathbf{r},T)italic_W ( bold_r , italic_T ). This effective potential is smooth in comparison to the initial disorder potential V⁒(𝐫)𝑉𝐫V(\mathbf{r})italic_V ( bold_r ) because W⁒(𝐫,T)π‘Šπ«π‘‡W(\mathbf{r},T)italic_W ( bold_r , italic_T ) is derived from the electron density n⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_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⁒(𝐫)𝑉𝐫V(\mathbf{r})italic_V ( bold_r ). Let us, therefore, obtain W⁒(𝐫,T)π‘Šπ«π‘‡W(\mathbf{r},T)italic_W ( bold_r , italic_T ) by subjecting V⁒(𝐫)𝑉𝐫V(\mathbf{r})italic_V ( bold_r ) to the action of a universal low-pass filter (ULF).

The key question is how to find out the appropriate T𝑇Titalic_T-dependent filter function Γ⁒(𝐫,T)Γ𝐫𝑇\Gamma(\mathbf{r},T)roman_Ξ“ ( bold_r , italic_T ) that can be used to extract the shape of the effective potential W⁒(𝐫,T)π‘Šπ«π‘‡W(\mathbf{r},T)italic_W ( bold_r , italic_T ) for a given realization of the white-noise potential V⁒(r)π‘‰π‘ŸV(r)italic_V ( italic_r ). One can show that the Fourier image of this function has the shape 9

Ξ“^⁒(k)=πλ⁒k⁒eβˆ’Ξ»2⁒k2/4⁒erfi⁒(λ⁒k/2),^Ξ“π‘˜πœ‹πœ†π‘˜superscript𝑒superscriptπœ†2superscriptπ‘˜24erfiπœ†π‘˜2\hat{\Gamma}(k)=\frac{\sqrt{\pi}}{\lambda k}\,e^{-\lambda^{2}k^{2}/4}\,\text{% erfi}(\lambda k/2)\;,over^ start_ARG roman_Ξ“ end_ARG ( italic_k ) = divide start_ARG square-root start_ARG italic_Ο€ end_ARG end_ARG start_ARG italic_Ξ» italic_k end_ARG italic_e start_POSTSUPERSCRIPT - italic_Ξ» start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_k start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT / 4 end_POSTSUPERSCRIPT erfi ( italic_Ξ» italic_k / 2 ) , (23)

where erfi is the imaginary error function, and Ξ»=ℏ/2⁒m⁒kB⁒Tπœ†Planck-constant-over-2-pi2π‘šsubscriptπ‘˜B𝑇\lambda=\hbar/\sqrt{2mk_{\rm B}T}italic_Ξ» = roman_ℏ / square-root start_ARG 2 italic_m italic_k start_POSTSUBSCRIPT roman_B end_POSTSUBSCRIPT italic_T end_ARG.

In order to reveal the electron density distribution n⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_T ) for a given realization of the white-noise potential V⁒(𝐫)𝑉𝐫V(\mathbf{r})italic_V ( bold_r ), one should first calculate the Fourier image V^⁒(𝐀)^𝑉𝐀\hat{V}(\mathbf{k})over^ start_ARG italic_V end_ARG ( bold_k ) of V⁒(𝐫)𝑉𝐫V(\mathbf{r})italic_V ( bold_r ) using a fast-Fourier-transform (FFT). This function V^⁒(𝐀)^𝑉𝐀\hat{V}(\mathbf{k})over^ start_ARG italic_V end_ARG ( bold_k ) should be then multiplied by Ξ“^⁒(|𝐀|)^Γ𝐀\hat{\Gamma}(|\mathbf{k}|)over^ start_ARG roman_Ξ“ end_ARG ( | bold_k | ), given by Eq. (23). The inverse Fourier transform of the product V^⁒(𝐀)⁒Γ^⁒(|𝐀|)^𝑉𝐀^Γ𝐀\hat{V}(\mathbf{k})\hat{\Gamma}(|\mathbf{k}|)over^ start_ARG italic_V end_ARG ( bold_k ) over^ start_ARG roman_Ξ“ end_ARG ( | bold_k | ) by the FFT yields the effective potential W⁒(𝐫,T)π‘Šπ«π‘‡W(\mathbf{r},T)italic_W ( bold_r , italic_T ) because the inverse Fourier transform converts a product into a convolution 9. Inserting W⁒(𝐫,T)π‘Šπ«π‘‡W(\mathbf{r},T)italic_W ( bold_r , italic_T ) into Eq. (22) gives the electron density n⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_T ). In Fig. 6 and Fig. 7, we compare the results for W⁒(𝐫,T)π‘Šπ«π‘‡W(\mathbf{r},T)italic_W ( bold_r , italic_T ) of the above procedure with the effective potentials obtained via Eq. (22) from the electron density n⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_T ) calculated using the exact solution of SchrΓΆdinger equation in one and two dimensions, respectively. Coordinate xπ‘₯xitalic_x and temperature T𝑇Titalic_T in the figures are measured in units given by Eqs. (2) and (3).

Refer to caption
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.

The data in Figs. 6 and 7 demonstrate the high accuracy of the approach based on the T𝑇Titalic_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.

Refer to caption
Figure 7: Comparison between the exact electron density n⁒(x,y,T)𝑛π‘₯𝑦𝑇n(x,y,T)italic_n ( italic_x , italic_y , italic_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.

3.2 Random-wave-functions (RWF) approach to calculate n⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_T )

3.2.1 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 19. A 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 h^=exp⁑[βˆ’H^/(2⁒kB⁒T)]^β„Ž^𝐻2subscriptπ‘˜π΅π‘‡\hat{h}=\exp[-\hat{H}/(2k_{B}T)]over^ start_ARG italic_h end_ARG = roman_exp [ - over^ start_ARG italic_H end_ARG / ( 2 italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_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 Ξ΅πœ€\varepsilonitalic_Ξ΅ contribute to the distribution with the probability proportional to exp⁑[βˆ’Ξ΅/(kB⁒T)]πœ€subscriptπ‘˜π΅π‘‡\exp[-\varepsilon/(k_{B}T)]roman_exp [ - italic_Ξ΅ / ( italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_T ) ]. The wave function is the probability amplitude, which explains the factor 1/2121/21 / 2 in the operator h^^β„Ž\hat{h}over^ start_ARG italic_h end_ARG. The wave function is always a linear combination of eigenfunctions that correspond to different energies. The action of the operator h^^β„Ž\hat{h}over^ start_ARG italic_h end_ARG suppresses the contributions of high-energy eigenfunctions in favor of the contributions of low-energy eigenfunctions. By the application of the operator h^^β„Ž\hat{h}over^ start_ARG italic_h end_ARG 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⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_T ). The question arises on how to numerically subject a wave function to the action of the operator h^=exp⁑[βˆ’H^/(2⁒kB⁒T)]^β„Ž^𝐻2subscriptπ‘˜π΅π‘‡\hat{h}=\exp[-\hat{H}/(2k_{B}T)]over^ start_ARG italic_h end_ARG = roman_exp [ - over^ start_ARG italic_H end_ARG / ( 2 italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_T ) ]. It can be done by recursively applying the Hamiltonian H^^𝐻\hat{H}over^ start_ARG italic_H end_ARG to the wave function:

h^=eβˆ’H^/2⁒kB⁒Tβ‰ˆ(1βˆ’Ξ±β’H^)M,^β„Žsuperscript𝑒^𝐻2subscriptπ‘˜π΅π‘‡superscript1𝛼^𝐻𝑀\hat{h}=e^{-\hat{H}/2k_{B}T}\approx(1-\alpha\hat{H})^{M}\;,over^ start_ARG italic_h end_ARG = italic_e start_POSTSUPERSCRIPT - over^ start_ARG italic_H end_ARG / 2 italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_T end_POSTSUPERSCRIPT β‰ˆ ( 1 - italic_Ξ± over^ start_ARG italic_H end_ARG ) start_POSTSUPERSCRIPT italic_M end_POSTSUPERSCRIPT , (24)

with a natural number 9 Mβ‰ˆ1/(2⁒α⁒kB⁒T)𝑀12𝛼subscriptπ‘˜π΅π‘‡M\approx 1/(2\alpha k_{B}T)italic_M β‰ˆ 1 / ( 2 italic_Ξ± italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_T ) and a small parameter α𝛼\alphaitalic_Ξ±. A simple analysis justifies the choice 9 Ξ±=1.5/Ο΅max𝛼1.5subscriptitalic-Ο΅max\alpha=1.5/\epsilon_{\text{max}}italic_Ξ± = 1.5 / italic_Ο΅ start_POSTSUBSCRIPT max end_POSTSUBSCRIPT, where Ο΅maxsubscriptitalic-Ο΅max\epsilon_{\text{max}}italic_Ο΅ start_POSTSUBSCRIPT max end_POSTSUBSCRIPT is the estimate for the upper boundary of the energy distribution. In the case of a regular grid with the lattice constant aπ‘Žaitalic_a, Ο΅max=ℏ2/(m⁒a2)subscriptitalic-Ο΅maxsuperscriptPlanck-constant-over-2-pi2π‘šsuperscriptπ‘Ž2\epsilon_{\text{max}}=\hbar^{2}/(ma^{2})italic_Ο΅ start_POSTSUBSCRIPT max end_POSTSUBSCRIPT = roman_ℏ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT / ( italic_m italic_a start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ). 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 Δ⁒VΔ𝑉\Delta Vroman_Ξ” italic_V per lattice site. The value of the random wave function Οˆπœ“\psiitalic_ψ on each lattice site is chosen independently as a random number extracted from a normal distribution with the average value zero and variance 1/Δ⁒V1Δ𝑉1/\Delta V1 / roman_Ξ” italic_V. The following transformation of the wave function Οˆπœ“\psiitalic_ψ,

Οˆβ†’Οˆβˆ’Ξ±β’H^β’Οˆβ†’πœ“πœ“π›Ό^π»πœ“\psi\rightarrow\psi-\alpha\hat{H}\psi\,italic_ψ β†’ italic_ψ - italic_Ξ± over^ start_ARG italic_H end_ARG italic_ψ (25)

is applied M𝑀Mitalic_M times. Then an estimate of the reduced electron density n~R⁒(𝐫,T)subscript~𝑛R𝐫𝑇\tilde{n}_{\rm R}(\mathbf{r},T)over~ start_ARG italic_n end_ARG start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT ( bold_r , italic_T ) is

n~R⁒(𝐫,T)=2⁒|ψ⁒(𝐫)|2.subscript~𝑛R𝐫𝑇2superscriptπœ“π«2\tilde{n}_{\rm R}(\mathbf{r},T)=2\,|\psi(\mathbf{r})|^{2}\,.over~ start_ARG italic_n end_ARG start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT ( bold_r , italic_T ) = 2 | italic_ψ ( bold_r ) | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT . (26)

The calculation of n~R⁒(𝐫,T)subscript~𝑛R𝐫𝑇\tilde{n}_{\rm R}(\mathbf{r},T)over~ start_ARG italic_n end_ARG start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT ( bold_r , italic_T ) is carried out for a large number NRsubscript𝑁RN_{\rm R}italic_N start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT of realizations RR{\rm R}roman_R of the random wave function ψ⁒(𝐫)πœ“π«\psi(\mathbf{r})italic_ψ ( bold_r ). Then, the electron density n⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_T ) is the arithmetic mean of the functions n~R⁒(𝐫)subscript~𝑛R𝐫\tilde{n}_{\rm R}(\mathbf{r})over~ start_ARG italic_n end_ARG start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT ( bold_r ) obtained for different realizations R, multiplied by a chemical-potential-related factor eΞΌ/(kB⁒T)superscriptπ‘’πœ‡subscriptπ‘˜π΅π‘‡e^{\mu/(k_{B}T)}italic_e start_POSTSUPERSCRIPT italic_ΞΌ / ( italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_T ) end_POSTSUPERSCRIPT,

n⁒(𝐫,T)=eΞΌ/(kB⁒T)⁒⟨n~R⁒(𝐫,T)⟩R.𝑛𝐫𝑇superscriptπ‘’πœ‡subscriptπ‘˜π΅π‘‡subscriptdelimited-⟨⟩subscript~𝑛R𝐫𝑇Rn(\mathbf{r},T)=e^{\mu/(k_{B}T)}\langle\tilde{n}_{\rm R}(\mathbf{r},T)\rangle_% {\rm R}\;.italic_n ( bold_r , italic_T ) = italic_e start_POSTSUPERSCRIPT italic_ΞΌ / ( italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_T ) end_POSTSUPERSCRIPT ⟨ over~ start_ARG italic_n end_ARG start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT ( bold_r , italic_T ) ⟩ start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT . (27)
Refer to caption
Figure 8: Comparison between exact (solid lines) and approximate (symbols) reduced electron density n~⁒(x,T)=eβˆ’ΞΌ/(kB⁒T)⁒n⁒(x,T)~𝑛π‘₯𝑇superscriptπ‘’πœ‡subscriptπ‘˜π΅π‘‡π‘›π‘₯𝑇\tilde{n}(x,T)=e^{-\mu/(k_{B}T)}n(x,T)over~ start_ARG italic_n end_ARG ( italic_x , italic_T ) = italic_e start_POSTSUPERSCRIPT - italic_ΞΌ / ( italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_T ) end_POSTSUPERSCRIPT italic_n ( italic_x , italic_T ) in a one-dimensional white-noise potential at different temperatures T𝑇Titalic_T and NR=1000subscript𝑁R1000N_{\rm R}=1000italic_N start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT = 1000 iterations. For clarity, the lowest three curves are scaled down by a multiplication by 10βˆ’3superscript10310^{-3}10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT, 10βˆ’2superscript10210^{-2}10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT and 10βˆ’1superscript10110^{-1}10 start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT, 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.
Refer to caption
Figure 9: Reduced electron density n~⁒(x,y,z,T)~𝑛π‘₯𝑦𝑧𝑇\tilde{n}(x,y,z,T)over~ start_ARG italic_n end_ARG ( italic_x , italic_y , italic_z , italic_T ) in a three-dimensional sample with white-noise potential. Compared are the exact density (solid orange line) with that obtained by the RWF algorithm (dashed green line for NR=20subscript𝑁R20N_{\rm R}=20italic_N start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT = 20 and dotted blue line for NR=1000subscript𝑁R1000N_{\rm R}=1000italic_N start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT = 1000). Profiles along the xπ‘₯xitalic_x-axis with different values of coordinates y𝑦yitalic_y and z𝑧zitalic_z are shown (see inset). Sample size is 2Γ—2Γ—22222\times 2\times 22 Γ— 2 Γ— 2 dimensionless units, the discretization grid parameter is a=0.1π‘Ž0.1a=0.1italic_a = 0.1, the temperature is T=T0𝑇subscript𝑇0T=T_{0}italic_T = italic_T start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT. 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.

The larger the number of realizations NRsubscript𝑁RN_{\rm R}italic_N start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT, the more accurate is the calculated electron density n⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_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.

4 Discussion

In this work, we introduce four theoretical tools to get access to the space- and temperature-dependent electron density n⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_T ) in disordered media with a random potential V⁒(𝐫)𝑉𝐫V(\mathbf{r})italic_V ( bold_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⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_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⁒(𝐫,T)𝑛𝐫𝑇n(\mathbf{r},T)italic_n ( bold_r , italic_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⁒(𝐫,T)π‘Šπ«π‘‡W(\mathbf{r},T)italic_W ( bold_r , italic_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.

Appendix

Appendix A Why do Methods 1 and 2 work

To understand why Method 1 works, let us apply to the Hamiltonian H^^𝐻\hat{H}over^ start_ARG italic_H end_ARG a unitary transformation that diagonalizes it. This transformation also diagonalizes all the matrices that appear in Method 1. The diagonal elements of matrices H^^𝐻\hat{H}over^ start_ARG italic_H end_ARG, A^0subscript^𝐴0\hat{A}_{0}over^ start_ARG italic_A end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, A^Nsubscript^𝐴𝑁\hat{A}_{N}over^ start_ARG italic_A end_ARG start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT and B^^𝐡\hat{B}over^ start_ARG italic_B end_ARG are, in accord to Eqs. (9) – (11),

Ha⁒asubscriptπ»π‘Žπ‘Ž\displaystyle H_{aa}italic_H start_POSTSUBSCRIPT italic_a italic_a end_POSTSUBSCRIPT =Ξ΅a,absentsubscriptπœ€π‘Ž\displaystyle=\varepsilon_{a}\,,= italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT , (28)
A0,a⁒asubscript𝐴0π‘Žπ‘Ž\displaystyle A_{0,aa}italic_A start_POSTSUBSCRIPT 0 , italic_a italic_a end_POSTSUBSCRIPT =Ξ΅aβˆ’Ξ΅0Ξ΅fβˆ’Ξ΅0,absentsubscriptπœ€π‘Žsubscriptπœ€0subscriptπœ€π‘“subscriptπœ€0\displaystyle=\frac{\varepsilon_{a}-\varepsilon_{0}}{\varepsilon_{f}-% \varepsilon_{0}}\,,= divide start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG , (29)
AN,a⁒asubscriptπ΄π‘π‘Žπ‘Ž\displaystyle A_{N,aa}italic_A start_POSTSUBSCRIPT italic_N , italic_a italic_a end_POSTSUBSCRIPT =(Ξ΅aβˆ’Ξ΅0Ξ΅fβˆ’Ξ΅0)2N,absentsuperscriptsubscriptπœ€π‘Žsubscriptπœ€0subscriptπœ€π‘“subscriptπœ€0superscript2𝑁\displaystyle=\left(\frac{\varepsilon_{a}-\varepsilon_{0}}{\varepsilon_{f}-% \varepsilon_{0}}\right)^{2^{N}}\,,= ( divide start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ) start_POSTSUPERSCRIPT 2 start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT , (30)
Ba⁒asubscriptπ΅π‘Žπ‘Ž\displaystyle B_{aa}italic_B start_POSTSUBSCRIPT italic_a italic_a end_POSTSUBSCRIPT =1(Ξ΅aβˆ’Ξ΅0Ξ΅fβˆ’Ξ΅0)2N+1.absent1superscriptsubscriptπœ€π‘Žsubscriptπœ€0subscriptπœ€π‘“subscriptπœ€0superscript2𝑁1\displaystyle=\frac{1}{\left(\frac{\varepsilon_{a}-\varepsilon_{0}}{% \varepsilon_{f}-\varepsilon_{0}}\right)^{2^{N}}+1}\,.= divide start_ARG 1 end_ARG start_ARG ( divide start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ) start_POSTSUPERSCRIPT 2 start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT + 1 end_ARG . (31)

Let us consider those values of Ξ΅asubscriptπœ€π‘Ž\varepsilon_{a}italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT that are close to the Fermi energy Ξ΅fsubscriptπœ€π‘“\varepsilon_{f}italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT. For them,

(Ξ΅aβˆ’Ξ΅0Ξ΅fβˆ’Ξ΅0)2N=(Ξ΅aβˆ’Ξ΅fΞ΅fβˆ’Ξ΅0+1)2Nβ‰ˆexp⁑(2N⁒Ρaβˆ’Ξ΅fΞ΅fβˆ’Ξ΅0).superscriptsubscriptπœ€π‘Žsubscriptπœ€0subscriptπœ€π‘“subscriptπœ€0superscript2𝑁superscriptsubscriptπœ€π‘Žsubscriptπœ€π‘“subscriptπœ€π‘“subscriptπœ€01superscript2𝑁superscript2𝑁subscriptπœ€π‘Žsubscriptπœ€π‘“subscriptπœ€π‘“subscriptπœ€0\left(\frac{\varepsilon_{a}-\varepsilon_{0}}{\varepsilon_{f}-\varepsilon_{0}}% \right)^{2^{N}}=\left(\frac{\varepsilon_{a}-\varepsilon_{f}}{\varepsilon_{f}-% \varepsilon_{0}}+1\right)^{2^{N}}\approx\exp\left(2^{N}\frac{\varepsilon_{a}-% \varepsilon_{f}}{\varepsilon_{f}-\varepsilon_{0}}\right).( divide start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ) start_POSTSUPERSCRIPT 2 start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT = ( divide start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT end_ARG start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG + 1 ) start_POSTSUPERSCRIPT 2 start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT β‰ˆ roman_exp ( 2 start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT divide start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT end_ARG start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ) . (32)

Taking Eq. (8) into account, one can rewrite this estimate as

(Ξ΅aβˆ’Ξ΅0Ξ΅fβˆ’Ξ΅0)2Nβ‰ˆexp⁑(Β±Ξ΅aβˆ’Ξ΅fkB⁒T),superscriptsubscriptπœ€π‘Žsubscriptπœ€0subscriptπœ€π‘“subscriptπœ€0superscript2𝑁plus-or-minussubscriptπœ€π‘Žsubscriptπœ€π‘“subscriptπ‘˜π΅π‘‡\left(\frac{\varepsilon_{a}-\varepsilon_{0}}{\varepsilon_{f}-\varepsilon_{0}}% \right)^{2^{N}}\approx\exp\left(\pm\,\frac{\varepsilon_{a}-\varepsilon_{f}}{k_% {B}T}\right),( divide start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ) start_POSTSUPERSCRIPT 2 start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT β‰ˆ roman_exp ( Β± divide start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT end_ARG start_ARG italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_T end_ARG ) , (33)

where sign β€œ+++” is chosen if Ξ΅0<Ξ΅fsubscriptπœ€0subscriptπœ€π‘“\varepsilon_{0}<\varepsilon_{f}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT < italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT, and sign β€œβˆ’--” otherwise. Substitution of Eq. (33) into Eq. (31) and comparison to expression (5) for the Fermi function f⁒(Ξ΅)π‘“πœ€f(\varepsilon)italic_f ( italic_Ξ΅ ) provides the following result:

Ba⁒aβ‰ˆf⁒(Ξ΅a)ifΞ΅0<Ξ΅f,formulae-sequencesubscriptπ΅π‘Žπ‘Žπ‘“subscriptπœ€π‘Žifsubscriptπœ€0subscriptπœ€π‘“B_{aa}\approx f(\varepsilon_{a})\quad\mathrm{if}\quad\varepsilon_{0}<% \varepsilon_{f}\,,italic_B start_POSTSUBSCRIPT italic_a italic_a end_POSTSUBSCRIPT β‰ˆ italic_f ( italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT ) roman_if italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT < italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT , (34)
1βˆ’Ba⁒aβ‰ˆf⁒(Ξ΅a)ifΞ΅0>Ξ΅f,formulae-sequence1subscriptπ΅π‘Žπ‘Žπ‘“subscriptπœ€π‘Žifsubscriptπœ€0subscriptπœ€π‘“1-B_{aa}\approx f(\varepsilon_{a})\quad\mathrm{if}\quad\varepsilon_{0}>% \varepsilon_{f}\,,1 - italic_B start_POSTSUBSCRIPT italic_a italic_a end_POSTSUBSCRIPT β‰ˆ italic_f ( italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT ) roman_if italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT > italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT , (35)

which is valid after the unitary transformation that diagonalizes the Hamiltonian. If we undo this transformation, Eq. (34) acquires the form

Bi⁒jβ‰ˆβˆ‘a=1Lf⁒(Ξ΅a)⁒ψa,i⁒ψa,j*ifΞ΅0<Ξ΅f,formulae-sequencesubscript𝐡𝑖𝑗superscriptsubscriptπ‘Ž1𝐿𝑓subscriptπœ€π‘Žsubscriptπœ“π‘Žπ‘–superscriptsubscriptπœ“π‘Žπ‘—ifsubscriptπœ€0subscriptπœ€π‘“B_{ij}\approx\sum_{a=1}^{L}f(\varepsilon_{a})\psi_{a,i}\psi_{a,j}^{*}\quad% \mathrm{if}\quad\varepsilon_{0}<\varepsilon_{f}\,,italic_B start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT β‰ˆ βˆ‘ start_POSTSUBSCRIPT italic_a = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_L end_POSTSUPERSCRIPT italic_f ( italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT ) italic_ψ start_POSTSUBSCRIPT italic_a , italic_i end_POSTSUBSCRIPT italic_ψ start_POSTSUBSCRIPT italic_a , italic_j end_POSTSUBSCRIPT start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT roman_if italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT < italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT , (36)

that is, the diagonal matrix elements of the matrix B^^𝐡\hat{B}over^ start_ARG italic_B end_ARG are equal to

Bi⁒iβ‰ˆβˆ‘a=1Lf⁒(Ξ΅a)⁒|ψa,i|2ifΞ΅0<Ξ΅f,formulae-sequencesubscript𝐡𝑖𝑖superscriptsubscriptπ‘Ž1𝐿𝑓subscriptπœ€π‘Žsuperscriptsubscriptπœ“π‘Žπ‘–2ifsubscriptπœ€0subscriptπœ€οΏ½οΏ½B_{ii}\approx\sum_{a=1}^{L}f(\varepsilon_{a})|\psi_{a,i}|^{2}\quad\mathrm{if}% \quad\varepsilon_{0}<\varepsilon_{f}\,,italic_B start_POSTSUBSCRIPT italic_i italic_i end_POSTSUBSCRIPT β‰ˆ βˆ‘ start_POSTSUBSCRIPT italic_a = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_L end_POSTSUPERSCRIPT italic_f ( italic_Ξ΅ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT ) | italic_ψ start_POSTSUBSCRIPT italic_a , italic_i end_POSTSUBSCRIPT | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_if italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT < italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT , (37)

or, due to Eq. (6),

Bi⁒iβ‰ˆΞ”β’V2⁒niifΞ΅0<Ξ΅f,formulae-sequencesubscript𝐡𝑖𝑖Δ𝑉2subscript𝑛𝑖ifsubscriptπœ€0subscriptπœ€π‘“B_{ii}\approx\frac{\Delta V}{2}\,n_{i}\quad\mathrm{if}\quad\varepsilon_{0}<% \varepsilon_{f}\,,italic_B start_POSTSUBSCRIPT italic_i italic_i end_POSTSUBSCRIPT β‰ˆ divide start_ARG roman_Ξ” italic_V end_ARG start_ARG 2 end_ARG italic_n start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT roman_if italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT < italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT , (38)

that justifies Eq. (12) of Method 1.

Similarly, from Eq. (35) one can get after undoing the unitary transformation:

1βˆ’Bi⁒iβ‰ˆΞ”β’V2⁒niifΞ΅0>Ξ΅f,formulae-sequence1subscript𝐡𝑖𝑖Δ𝑉2subscript𝑛𝑖ifsubscriptπœ€0subscriptπœ€π‘“1-B_{ii}\approx\frac{\Delta V}{2}\,n_{i}\quad\mathrm{if}\quad\varepsilon_{0}>% \varepsilon_{f}\,,1 - italic_B start_POSTSUBSCRIPT italic_i italic_i end_POSTSUBSCRIPT β‰ˆ divide start_ARG roman_Ξ” italic_V end_ARG start_ARG 2 end_ARG italic_n start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT roman_if italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT > italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT , (39)

which justifies Eq. (13) of Method 1.

Now let us consider Method 2. For simplicity, we restrict ourselves to the case Ξ΅0<Ξ΅fsubscriptπœ€0subscriptπœ€π‘“\varepsilon_{0}<\varepsilon_{f}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT < italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT. The opposite case is completely analogous.

The matrix X^^𝑋\hat{X}over^ start_ARG italic_X end_ARG can be expressed as X^=B^⁒U^^𝑋^𝐡^π‘ˆ\hat{X}=\hat{B}\hat{U}over^ start_ARG italic_X end_ARG = over^ start_ARG italic_B end_ARG over^ start_ARG italic_U end_ARG according to Eqs. (11) and (14). Hence, the sum in the right-hand side of Eq. (15) can be rewritten as

βˆ‘a=1NCXi⁒a⁒Ui⁒a=βˆ‘a=1NCβˆ‘j=1LBi⁒j⁒Uj⁒a⁒Ui⁒a.superscriptsubscriptπ‘Ž1subscript𝑁𝐢subscriptπ‘‹π‘–π‘Žsubscriptπ‘ˆπ‘–π‘Žsuperscriptsubscriptπ‘Ž1subscript𝑁𝐢superscriptsubscript𝑗1𝐿subscript𝐡𝑖𝑗subscriptπ‘ˆπ‘—π‘Žsubscriptπ‘ˆπ‘–π‘Ž\sum_{a=1}^{N_{C}}X_{ia}U_{ia}=\sum_{a=1}^{N_{C}}\sum_{j=1}^{L}B_{ij}U_{ja}U_{% ia}\,.βˆ‘ start_POSTSUBSCRIPT italic_a = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT end_POSTSUPERSCRIPT italic_X start_POSTSUBSCRIPT italic_i italic_a end_POSTSUBSCRIPT italic_U start_POSTSUBSCRIPT italic_i italic_a end_POSTSUBSCRIPT = βˆ‘ start_POSTSUBSCRIPT italic_a = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT end_POSTSUPERSCRIPT βˆ‘ start_POSTSUBSCRIPT italic_j = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_L end_POSTSUPERSCRIPT italic_B start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT italic_U start_POSTSUBSCRIPT italic_j italic_a end_POSTSUBSCRIPT italic_U start_POSTSUBSCRIPT italic_i italic_a end_POSTSUBSCRIPT . (40)

For a given i𝑖iitalic_i, there is only one value of aπ‘Žaitalic_a such that Ui⁒a=1subscriptπ‘ˆπ‘–π‘Ž1U_{ia}=1italic_U start_POSTSUBSCRIPT italic_i italic_a end_POSTSUBSCRIPT = 1, according to construction of the matrix U^^π‘ˆ\hat{U}over^ start_ARG italic_U end_ARG (see Section 2.2). For this value of aπ‘Žaitalic_a, there are only few values of j𝑗jitalic_j such that Uj⁒a=1subscriptπ‘ˆπ‘—π‘Ž1U_{ja}=1italic_U start_POSTSUBSCRIPT italic_j italic_a end_POSTSUBSCRIPT = 1, as illustrated in Fig 10. Let us denote the set of these values of j𝑗jitalic_j as Sisubscript𝑆𝑖S_{i}italic_S start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT. One can therefore rewrite Eq. (40) as

Refer to caption
Figure 10: Illustration of set Sisubscript𝑆𝑖S_{i}italic_S start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT for a given row index i𝑖iitalic_i of matrix U^^π‘ˆ\hat{U}over^ start_ARG italic_U end_ARG. White and black squares represent zeros and ones, respectively. Row i𝑖iitalic_i and column aπ‘Žaitalic_a of the matrix are highlighted in green.
βˆ‘a=1NCXi⁒a⁒Ui⁒a=βˆ‘j∈SiBi⁒j.superscriptsubscriptπ‘Ž1subscript𝑁𝐢subscriptπ‘‹π‘–π‘Žsubscriptπ‘ˆπ‘–π‘Žsubscript𝑗subscript𝑆𝑖subscript𝐡𝑖𝑗\sum_{a=1}^{N_{C}}X_{ia}U_{ia}=\sum_{j\in S_{i}}B_{ij}\,.βˆ‘ start_POSTSUBSCRIPT italic_a = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT end_POSTSUPERSCRIPT italic_X start_POSTSUBSCRIPT italic_i italic_a end_POSTSUBSCRIPT italic_U start_POSTSUBSCRIPT italic_i italic_a end_POSTSUBSCRIPT = βˆ‘ start_POSTSUBSCRIPT italic_j ∈ italic_S start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUBSCRIPT italic_B start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT . (41)

Note that i∈Si𝑖subscript𝑆𝑖i\in S_{i}italic_i ∈ italic_S start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT since Ui⁒a=1subscriptπ‘ˆπ‘–π‘Ž1U_{ia}=1italic_U start_POSTSUBSCRIPT italic_i italic_a end_POSTSUBSCRIPT = 1. Therefore the sum in the right-hand side of Eq. (41) contains the term Bi⁒isubscript𝐡𝑖𝑖B_{ii}italic_B start_POSTSUBSCRIPT italic_i italic_i end_POSTSUBSCRIPT.

Let us argue now that, at large enough NCsubscript𝑁𝐢N_{C}italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT, all other terms Bi⁒jsubscript𝐡𝑖𝑗B_{ij}italic_B start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT in this sum with jβ‰ i𝑗𝑖j\neq iitalic_j β‰  italic_i are negligible. One can see from Eq. (36) that the matrix B^^𝐡\hat{B}over^ start_ARG italic_B end_ARG is close to the projector P^^𝑃\hat{P}over^ start_ARG italic_P end_ARG 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 Pi⁒jsubscript𝑃𝑖𝑗P_{ij}italic_P start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT of this projector decrease with increasing the distance between nodes i𝑖iitalic_i and j𝑗jitalic_j, and become negligible at some distance. Hence, the matrix elements Bi⁒jsubscript𝐡𝑖𝑗B_{ij}italic_B start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT in the right-hand side of Eq. (41) becomes negligible if nodes i𝑖iitalic_i and j𝑗jitalic_j are far away from each other, i. e., if iβ‰ j𝑖𝑗i\neq jitalic_i β‰  italic_j. The only non-negligible term is Bi⁒isubscript𝐡𝑖𝑖B_{ii}italic_B start_POSTSUBSCRIPT italic_i italic_i end_POSTSUBSCRIPT, and therefore,

βˆ‘a=1NCXi⁒a⁒Ui⁒aβ‰ˆBi⁒i.superscriptsubscriptπ‘Ž1subscript𝑁𝐢subscriptπ‘‹π‘–π‘Žsubscriptπ‘ˆπ‘–π‘Žsubscript𝐡𝑖𝑖\sum_{a=1}^{N_{C}}X_{ia}U_{ia}\approx B_{ii}\,.βˆ‘ start_POSTSUBSCRIPT italic_a = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT end_POSTSUPERSCRIPT italic_X start_POSTSUBSCRIPT italic_i italic_a end_POSTSUBSCRIPT italic_U start_POSTSUBSCRIPT italic_i italic_a end_POSTSUBSCRIPT β‰ˆ italic_B start_POSTSUBSCRIPT italic_i italic_i end_POSTSUBSCRIPT . (42)

Substituting there Bi⁒isubscript𝐡𝑖𝑖B_{ii}italic_B start_POSTSUBSCRIPT italic_i italic_i end_POSTSUBSCRIPT 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 Ξ΅0>Ξ΅fsubscriptπœ€0subscriptπœ€π‘“\varepsilon_{0}>\varepsilon_{f}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT > italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT.

Appendix B Choice of parameters Ξ΅0subscriptπœ€0\varepsilon_{0}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, N𝑁Nitalic_N and NCsubscript𝑁𝐢N_{C}italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT

There are two restrictions on the choice of parameters Ξ΅0subscriptπœ€0\varepsilon_{0}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT and N𝑁Nitalic_N and, hence, on the temperature value T𝑇Titalic_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 (A^N+I^)subscript^𝐴𝑁^𝐼(\hat{A}_{N}+\hat{I})( over^ start_ARG italic_A end_ARG start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT + over^ start_ARG italic_I end_ARG ) must not be ill-conditioned.

To consider the first restriction, let us define the function f~⁒(Ξ΅)~π‘“πœ€\tilde{f}(\varepsilon)over~ start_ARG italic_f end_ARG ( italic_Ξ΅ )

f~⁒(Ξ΅)=1(Ξ΅βˆ’Ξ΅0Ξ΅fβˆ’Ξ΅0)2N+1ifΞ΅0<Ξ΅f,formulae-sequence~π‘“πœ€1superscriptπœ€subscriptπœ€0subscriptπœ€π‘“subscriptπœ€0superscript2𝑁1ifsubscriptπœ€0subscriptπœ€π‘“\tilde{f}(\varepsilon)=\frac{1}{\left(\frac{\varepsilon-\varepsilon_{0}}{% \varepsilon_{f}-\varepsilon_{0}}\right)^{2^{N}}+1}\quad\mathrm{if}\quad% \varepsilon_{0}<\varepsilon_{f}\,,over~ start_ARG italic_f end_ARG ( italic_Ξ΅ ) = divide start_ARG 1 end_ARG start_ARG ( divide start_ARG italic_Ξ΅ - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ) start_POSTSUPERSCRIPT 2 start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT + 1 end_ARG roman_if italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT < italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT , (43)
f~⁒(Ξ΅)=1βˆ’1(Ξ΅βˆ’Ξ΅0Ξ΅fβˆ’Ξ΅0)2N+1ifΞ΅0>Ξ΅f.formulae-sequence~π‘“πœ€11superscriptπœ€subscriptπœ€0subscriptπœ€π‘“subscriptπœ€0superscript2𝑁1ifsubscriptπœ€0subscriptπœ€π‘“\tilde{f}(\varepsilon)=1-\frac{1}{\left(\frac{\varepsilon-\varepsilon_{0}}{% \varepsilon_{f}-\varepsilon_{0}}\right)^{2^{N}}+1}\quad\mathrm{if}\quad% \varepsilon_{0}>\varepsilon_{f}\,.over~ start_ARG italic_f end_ARG ( italic_Ξ΅ ) = 1 - divide start_ARG 1 end_ARG start_ARG ( divide start_ARG italic_Ξ΅ - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ) start_POSTSUPERSCRIPT 2 start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT + 1 end_ARG roman_if italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT > italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT . (44)

Methods 1 and 2 are based on the fact that this function is close to the Fermi function f⁒(Ξ΅)π‘“πœ€f(\varepsilon)italic_f ( italic_Ξ΅ ) in a vicinity of the Fermi energy Ξ΅fsubscriptπœ€π‘“\varepsilon_{f}italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT, see Eqs. (31), (34) and (35). Let us now consider the behavior of the function f~⁒(Ξ΅)~π‘“πœ€\tilde{f}(\varepsilon)over~ start_ARG italic_f end_ARG ( italic_Ξ΅ ) in a broad range of energies Ξ΅πœ€\varepsilonitalic_Ξ΅.

As an example, Fig. 11 shows f~⁒(Ξ΅)~π‘“πœ€\tilde{f}(\varepsilon)over~ start_ARG italic_f end_ARG ( italic_Ξ΅ ) along with the Fermi function f⁒(Ξ΅)π‘“πœ€f(\varepsilon)italic_f ( italic_Ξ΅ ) for N=4𝑁4N=4italic_N = 4 and two choices of the β€œreference energy” Ξ΅0subscriptπœ€0\varepsilon_{0}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT: Ξ΅0<Ξ΅fsubscriptπœ€0subscriptπœ€π‘“\varepsilon_{0}<\varepsilon_{f}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT < italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT (upper panel) and Ξ΅0>Ξ΅fsubscriptπœ€0subscriptπœ€π‘“\varepsilon_{0}>\varepsilon_{f}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT > italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT (lower panel). In both cases, there is a discrepancy between the functions f⁒(Ξ΅)π‘“πœ€f(\varepsilon)italic_f ( italic_Ξ΅ ) and f~⁒(Ξ΅)~π‘“πœ€\tilde{f}(\varepsilon)over~ start_ARG italic_f end_ARG ( italic_Ξ΅ ) below the energy 2⁒Ρ0βˆ’Ξ΅f2subscriptπœ€0subscriptπœ€π‘“2\varepsilon_{0}-\varepsilon_{f}2 italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT in the first case, and above the energy 2⁒Ρ0βˆ’Ξ΅f2subscriptπœ€0subscriptπœ€π‘“2\varepsilon_{0}-\varepsilon_{f}2 italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT 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<Ξ΅fsubscriptπœ€0subscriptπœ€π‘“\varepsilon_{0}<\varepsilon_{f}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT < italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT, the lowest electron energy Ξ΅minsubscriptπœ€min\varepsilon_{\mathrm{min}}italic_Ξ΅ start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT must be larger than 2⁒Ρ0βˆ’Ξ΅f2subscriptπœ€0subscriptπœ€π‘“2\varepsilon_{0}-\varepsilon_{f}2 italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT. Similarly, in the case of Ξ΅0>Ξ΅fsubscriptπœ€0subscriptπœ€π‘“\varepsilon_{0}>\varepsilon_{f}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT > italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT, the highest electron energy Ξ΅maxsubscriptπœ€max\varepsilon_{\mathrm{max}}italic_Ξ΅ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT must be smaller than 2⁒Ρ0βˆ’Ξ΅f2subscriptπœ€0subscriptπœ€π‘“2\varepsilon_{0}-\varepsilon_{f}2 italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT. These conditions can be rewritten as restrictions to the β€œreference energy” Ξ΅0subscriptπœ€0\varepsilon_{0}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT:

Ξ΅0<Ξ΅f+Ξ΅min2orΞ΅0>Ξ΅f+Ξ΅max2.formulae-sequencesubscriptπœ€0subscriptπœ€π‘“subscriptπœ€min2orsubscriptπœ€0subscriptπœ€π‘“subscriptπœ€max2\varepsilon_{0}<\frac{\varepsilon_{f}+\varepsilon_{\mathrm{min}}}{2}\quad% \mathrm{or}\quad\varepsilon_{0}>\frac{\varepsilon_{f}+\varepsilon_{\mathrm{max% }}}{2}\,.italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT < divide start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT + italic_Ξ΅ start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG roman_or italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT > divide start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT + italic_Ξ΅ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG . (45)
Refer to caption
Figure 11: Comparison of the Fermi function f⁒(Ξ΅)π‘“πœ€f(\varepsilon)italic_f ( italic_Ξ΅ ) and function f~⁒(Ξ΅)~π‘“πœ€\tilde{f}(\varepsilon)over~ start_ARG italic_f end_ARG ( italic_Ξ΅ ). Above: the case of Ξ΅0<Ξ΅fsubscriptπœ€0subscriptπœ€π‘“\varepsilon_{0}<\varepsilon_{f}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT < italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT, below: the case of Ξ΅0>Ξ΅fsubscriptπœ€0subscriptπœ€π‘“\varepsilon_{0}>\varepsilon_{f}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT > italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT.

Now let us consider the second restriction. The inversion of the matrix (A^N+I^)subscript^𝐴𝑁^𝐼(\hat{A}_{N}+\hat{I})( over^ start_ARG italic_A end_ARG start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT + over^ start_ARG italic_I end_ARG ) 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 1015superscript101510^{15}10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT (a number that corresponds to the accuracy of representation of real numbers in the computer memory). The minimal eigenvalue of matrix (A^N+I^)subscript^𝐴𝑁^𝐼(\hat{A}_{N}+\hat{I})( over^ start_ARG italic_A end_ARG start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT + over^ start_ARG italic_I end_ARG ) is close to unity, and corresponds to energy levels near to Ξ΅0subscriptπœ€0\varepsilon_{0}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT. The largest eigenvalue corresponds to the energy level farthest from Ξ΅0subscriptπœ€0\varepsilon_{0}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, i. e., either to Ξ΅minsubscriptπœ€min\varepsilon_{\mathrm{min}}italic_Ξ΅ start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT or to Ξ΅maxsubscriptπœ€max\varepsilon_{\mathrm{max}}italic_Ξ΅ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT, and is approximately the largest of two values [(Ξ΅minβˆ’Ξ΅0)/(Ξ΅fβˆ’Ξ΅0)]2Nsuperscriptdelimited-[]subscriptπœ€minsubscriptπœ€0subscriptπœ€π‘“subscriptπœ€0superscript2𝑁\left[(\varepsilon_{\mathrm{min}}-\varepsilon_{0})/(\varepsilon_{f}-% \varepsilon_{0})\right]^{2^{N}}[ ( italic_Ξ΅ start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) / ( italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) ] start_POSTSUPERSCRIPT 2 start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT and [(Ξ΅maxβˆ’Ξ΅0)/(Ξ΅fβˆ’Ξ΅0)]2Nsuperscriptdelimited-[]subscriptπœ€maxsubscriptπœ€0subscriptπœ€π‘“subscriptπœ€0superscript2𝑁\left[(\varepsilon_{\mathrm{max}}-\varepsilon_{0})/(\varepsilon_{f}-% \varepsilon_{0})\right]^{2^{N}}[ ( italic_Ξ΅ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) / ( italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) ] start_POSTSUPERSCRIPT 2 start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT. Hence, the parameters Ξ΅0subscriptπœ€0\varepsilon_{0}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT and N𝑁Nitalic_N have to obey the following restrictions:

(Ξ΅minβˆ’Ξ΅0Ξ΅fβˆ’Ξ΅0)2N<1015and(Ξ΅maxβˆ’Ξ΅0Ξ΅fβˆ’Ξ΅0)2N<1015.formulae-sequencesuperscriptsubscriptπœ€minsubscriptπœ€0subscriptπœ€π‘“subscriptπœ€0superscript2𝑁superscript1015andsuperscriptsubscriptπœ€maxsubscriptπœ€0subscriptπœ€π‘“subscriptπœ€0superscript2𝑁superscript1015\left(\frac{\varepsilon_{\mathrm{min}}-\varepsilon_{0}}{\varepsilon_{f}-% \varepsilon_{0}}\right)^{2^{N}}<10^{15}\quad\mathrm{and}\quad\left(\frac{% \varepsilon_{\mathrm{max}}-\varepsilon_{0}}{\varepsilon_{f}-\varepsilon_{0}}% \right)^{2^{N}}<10^{15}\,.( divide start_ARG italic_Ξ΅ start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ) start_POSTSUPERSCRIPT 2 start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT < 10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT roman_and ( divide start_ARG italic_Ξ΅ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ) start_POSTSUPERSCRIPT 2 start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT < 10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT . (46)

The choice of parameters Ξ΅0subscriptπœ€0\varepsilon_{0}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT and N𝑁Nitalic_N is based on equations (45) and (46). We consider two different options: (i) when the temperature T𝑇Titalic_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𝑁Nitalic_N. For each such number N𝑁Nitalic_N, one finds Ξ΅0subscriptπœ€0\varepsilon_{0}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT according to Eq. (8) and check whether conditions (45) and (46) are fulfilled. The smallest suitable value of N𝑁Nitalic_N is the best choice. Indeed, the smaller N𝑁Nitalic_N is, the more sparse matrix (A^N+I^)subscript^𝐴𝑁^𝐼\left(\hat{A}_{N}+\hat{I}\right)( over^ start_ARG italic_A end_ARG start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT + over^ start_ARG italic_I end_ARG ) 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𝑁Nitalic_N is desirable, in order to minimize the temperature T𝑇Titalic_T according to Eq. (8). To achieve the largest T𝑇Titalic_T, one has to choose the value Ξ΅0subscriptπœ€0\varepsilon_{0}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT that maximizes the denominator |Ξ΅fβˆ’Ξ΅0|subscriptπœ€π‘“subscriptπœ€0|\varepsilon_{f}-\varepsilon_{0}|| italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT | in Eq. (46):

Ξ΅0β‰ˆΞ΅f+Ξ΅min2orΞ΅0β‰ˆΞ΅f+Ξ΅max2,formulae-sequencesubscriptπœ€0subscriptπœ€π‘“subscriptπœ€min2orsubscriptπœ€0subscriptπœ€π‘“subscriptπœ€max2\varepsilon_{0}\approx\frac{\varepsilon_{f}+\varepsilon_{\mathrm{min}}}{2}% \quad\mathrm{or}\quad\varepsilon_{0}\approx\frac{\varepsilon_{f}+\varepsilon_{% \mathrm{max}}}{2}\,,italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT β‰ˆ divide start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT + italic_Ξ΅ start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG roman_or italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT β‰ˆ divide start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT + italic_Ξ΅ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG , (47)

and then one needs to choose the maximal number N𝑁Nitalic_N that obeys the restrictions (46),

Nβ‰ˆlog2⁑(15log10⁑|Ξ΅maxβˆ’Ξ΅0Ξ΅fβˆ’Ξ΅0|)𝑁subscript215subscript10subscriptπœ€maxsubscriptπœ€0subscriptπœ€π‘“subscriptπœ€0N\approx\log_{2}\left(\frac{15}{\log_{10}\left|\frac{\varepsilon_{\mathrm{max}% }-\varepsilon_{0}}{\varepsilon_{f}-\varepsilon_{0}}\right|}\right)italic_N β‰ˆ roman_log start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ( divide start_ARG 15 end_ARG start_ARG roman_log start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT | divide start_ARG italic_Ξ΅ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG | end_ARG ) (48)

or

Nβ‰ˆlog2⁑(15log10⁑|Ξ΅minβˆ’Ξ΅0Ξ΅fβˆ’Ξ΅0|)𝑁subscript215subscript10subscriptπœ€minsubscriptπœ€0subscriptπœ€π‘“subscriptπœ€0N\approx\log_{2}\left(\frac{15}{\log_{10}\left|\frac{\varepsilon_{\mathrm{min}% }-\varepsilon_{0}}{\varepsilon_{f}-\varepsilon_{0}}\right|}\right)italic_N β‰ˆ roman_log start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ( divide start_ARG 15 end_ARG start_ARG roman_log start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT | divide start_ARG italic_Ξ΅ start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_Ξ΅ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG | end_ARG ) (49)

for the first and the second choice of Ξ΅0subscriptπœ€0\varepsilon_{0}italic_Ξ΅ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT in Eq. (47), correspondingly.

Finally, let us consider the choice of the parameter NCsubscript𝑁𝐢N_{C}italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT in Method 2. The accuracy of Method 2 improves with NCsubscript𝑁𝐢N_{C}italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT. However, larger NCsubscript𝑁𝐢N_{C}italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT give rise to a longer computation time due to the increase of the size of the matrix U^^π‘ˆ\hat{U}over^ start_ARG italic_U end_ARG in Eq. (14). Hence, there is a trade-off between the accuracy and the efficiency. The best value of parameter NCsubscript𝑁𝐢N_{C}italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT can be estimated as

NC≃(β„“/a)d,similar-to-or-equalssubscript𝑁𝐢superscriptβ„“π‘Žπ‘‘N_{C}\simeq(\ell/a)^{d},italic_N start_POSTSUBSCRIPT italic_C end_POSTSUBSCRIPT ≃ ( roman_β„“ / italic_a ) start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT , (50)

where d𝑑ditalic_d is dimensionality of the space, aπ‘Žaitalic_a is the distance between neighboring sites in the lattice, and β„“β„“\ellroman_β„“ is a characteristic length of decay of the matrix element Bi⁒jsubscript𝐡𝑖𝑗B_{ij}italic_B start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT with increasing the distance between sites i𝑖iitalic_i and j𝑗jitalic_j.

{acknowledgement}

S.D.B. and K.M. acknowledge financial support by the Deutsche Forschungsgemeinschaft (Research Training Group β€œTIDE”, RTG2591) as well as by the key profile area β€œQuantum Matter and Materials (QM2)” at the University of Cologne. K.M. further acknowledges support by the DFG through the project ASTRAL (ME1246-42).

References

  • Baranovski 2006 Baranovski, S. D., Ed. Charge Transport in Disordered Solids with Applications in Electronics; John Wiley and Sons, Ltd, Chichester, 2006
  • Masenda et al. 2021 Masenda, H.; Schneider, L. M.; Adel Aly, M.; Machchhar, S. J.; Usman, A.; Meerholz, K.; Gebhard, F.; Baranovskii, S. D.; Koch, M. Energy Scaling of Compositional Disorder in Ternary Transition-Metal Dichalcogenide Monolayers. Adv. Electron. Mater. 2021, 7, 2100196
  • Baranovskii et al. 2022 Baranovskii, S. D.; Nenashev, A. V.; Hertel, D.; Gebhard, F.; Meerholz, K. Energy Scales of Compositional Disorder in Alloy Semiconductors. ACS Omega 2022, 7, 45741
  • Arnold et al. 2016 Arnold, D. N.; David, G.; Jerison, D.; Mayboroda, S.; Filoche, M. Effective Confining Potential of Quantum States in Disordered Media. Phys. Rev. Lett. 2016, 116, 056602
  • Filoche et al. 2017 Filoche, M.; Piccardo, M.; Wu, Y.-R.; Li, C.-K.; Weisbuch, C.; Mayboroda, S. Localization Landscape Theory of disorder in semiconductors. I. Theory and modeling. Phys. Rev. B 2017, 95, 144204
  • Piccardo et al. 2017 Piccardo, M.; Li, C.-K.; Wu, Y.-R.; Speck, J. S.; Bonef, B.; Farrell, R. M.; Filoche, M.; Martinelli, L.; Peretti, J.; Weisbuch, C. Localization Landscape Theory of disorder in semiconductors. II. Urbach tails of disordered quantum well layers. Phys. Rev. B 2017, 95, 144205
  • Li et al. 2017 Li, C.-K.; Piccardo, M.; Lu, L.-S.; Mayboroda, S.; Martinelli, L.; Peretti, J.; Speck, J. S.; Weisbuch, C.; Filoche, M.; Wu, Y.-R. Localization Landscape Theory of disorder in semiconductors. III. Application to carrier transport and recombination in light emitting diodes. Phys. Rev. B 2017, 95, 144206
  • Gebhard et al. 2023 Gebhard, F.; Nenashev, A. V.; Meerholz, K.; Baranovskii, S. D. Quantum states in disordered media. I. Low-pass filter approach. Phys. Rev. B 2023, 107, 064206
  • Nenashev et al. 2023 Nenashev, A. V.; Baranovskii, S. D.; Meerholz, K.; Gebhard, F. Quantum states in disordered media. II. Spatial charge carrier distribution. Phys. Rev. B 2023, 107, 064207
  • Halperin and Lax 1966 Halperin, B. I.; Lax, M. Impurity-Band Tails in the High-Density Limit. I. Minimum Counting Methods. Phys. Rev. 1966, 148, 722–740
  • Goedecker 1999 Goedecker, S. Linear scaling electronic structure methods. Rev. Mod. Phys. 1999, 71, 1085–1123
  • Ceriotti et al. 2008 Ceriotti, M.; KΓΌhne, T. D.; Parrinello, M. An efficient and accurate decomposition of the Fermi operator. J. Chem. Phys. 2008, 129, 024707
  • Lin et al. 2009 Lin, L.; Lu, J.; Car, R.; E, W. Multipole representation of the Fermi operator with application to the electronic structure analysis of metallic systems. Phys. Rev. B 2009, 79, 115133
  • Bowler and Miyazaki 2012 Bowler, D. R.; Miyazaki, T. O⁒(N)𝑂𝑁O(N)italic_O ( italic_N ) methods in electronic structure calculations. Rep. Progr. Phys. 2012, 75, 036503
  • Lin et al. 2009 Lin, L.; Lu, J.; Ying, L.; Car, R.; E, W. Fast algorithm for extracting the diagonal of the inverse matrix with application to the electronic structure analysis of metallic systems. Commun. Math. Sci. 2009, 7, 755–777
  • Tang and Saad 2012 Tang, J. M.; Saad, Y. A probing method for computing the diagonal of a matrix inverse. Numerical Linear Algebra with Applications 2012, 19, 485–501
  • Li et al. 2013 Li, S.; Wu, W.; Darve, E. A fast algorithm for sparse matrix computations related to inversion. J. Comput. Phys. 2013, 242, 915–945
  • Baranovskii and Efros 1978 Baranovskii, S. D.; Efros, A. L. Band edge smearing in solid solutions. Sov. Phys. Semicond. 1978, 12, 1328–1330
  • Lu and Steinerberger 2018 Lu, J.; Steinerberger, S. Detecting localized eigenstates of linear operators. Res. Math. Sci. 2018, 5, 33
  • Krajewski and Parrinello 2005 Krajewski, F. R.; Parrinello, M. Stochastic linear scaling for metals and nonmetals. Phys. Rev. B 2005, 71, 233105