Simulations of the Electron Spectrum of Quantum Wires in n-Si of Arbitrarily Doping Profile by Thomas-Fermi Method

Electron spectrum in doped n-Si quantum wires is calculated by the Thomas-Fermi (TF) method under finite temperatures. The many-body exchange corrections are taken into account. The doping profile is arbitrary. At the first stage, the electron potential energy is calculated from a simple two-dimensional equation. The effective iteration scheme is proposed there that is valid for multidimensional problems. Then the energy levels and wave functions of this quantum well are simulated from the Schrödinger equations. The expansion by the full set of eigenfunctions of the linear harmonic oscillator is used. The quantum mechanical perturbation theory can be utilized to compute the energy levels. Generally, the perturbation theory for degenerate energy levels should be used.


Introduction
Investigations of the electron spectra of low-dimensional and highly doped structures are important to many nanotechnology applications including nanoscale transistors, sensors and modern plasmonics [1]- [6].Quantum structures based on silicon are used in nanoscale transistors, quantum logic devices, and as building blocks in nanostructures [1] [2] [3] [5].The optical transitions between electron levels within δ-doped quantum structures are utilized as infrared and terahertz modulators, detectors and lasers [4] [5].
The combined Schrödinger-Poisson method is widely used to compute spectra of low-dimensional structures [4].Namely, the Schrödinger equation for the electron spectrum is utilized jointly with the Poisson equation for the electron potential energy; the density of the electric charge, or the electron concentration, is expressed through the computed wave functions.The many-body corrections are expressed through the electron concentration and are added to the electron potential energy.But there are some lacks when using this method directly.This method needs the calculations of the whole spectrum i.e. energy levels and wave functions, and the electron potential energy simultaneously, so the initial values of the electron spectrum should be determined quite accurately.Otherwise, the convergence may be poor or even absent.Moreover, in two-dimensional (2D) and 3D quantum structures, or quantum wires and dots, the degeneration of energy levels may occur; this fact also can affect the convergence.
Investigations of δ-doped quantum structures are possible with a simpler approach based on the statistical Thomas-Fermi (TF) method [7] [8] with the many-body corrections.The preference of TF method is a sequential calculation of the electron potential energy and the electron concentration; only then the electron spectrum is under simulations.While this method provides quite exact calculations itself, the final results can be specified then with the Schrödinger-Poisson method.
In this paper, it is used the TF method to calculate the electron spectrum in highly doped n-Si quantum wires under finite temperatures T, where the many-body effects like exchange-correlations are taken into account [8].The 2D shape of the doping profile is arbitrary.It differs from [8].An effective iteration algorithm is proposed for solving the 2D equation for the electron potential energy.This algorithm can be used also for 3D problems like quantum dots.
Then the electron eigenvalues and eigenfunctions of the quantum wire are calculated with an expansion by the harmonic oscillator functions.In a general case of 2D and 3D problems, the perturbation theory should be developed for degenerate energy levels.

Thomas-Fermi Method to Compute Electron Potential Energy
Consider a single δ-doped electron quantum wire of n-Si.The arbitrary high doping profile is considered in the plane x, y perpendicular to the axis of the wire OZ.The atomic units are used for distances , ν = 6 is the number of the lowest electron valleys in Si [9].In n-Si the lowest valleys are lateral, and the effective mass is anisotropic m || , m ⊥ .
At the first stage the electrostatic electron energy V is calculated by the TF method.The single equation for δ-doped electron quantum wire with dimensionless variables is [8]: where Here n is the total electron concentration; N 1d and N d0 are 1D and bulk donor concentrations respectively.Φ 1/2 (v) is the Fermi integral that has been approximated from [5].The donor levels are assumed shallow and single charged.The concentration of 1D donors is high N 1d0 ≥ 10 20 cm −3 ; they are fully ionized.The 1D doping is localized at the distances x 0 , y 0 ≈ 1 -5 nm and depends on two coordinates x, y.The value of V xc is the many-body correction to the electron energy due to exchange-correlations [4] [7].Results of simulations do not depend on the value of the critical electron concentration in the case of n c ≤ 10 18 cm −3 .Because in n-Si the achievable doping concentrations are N 1d ≤ 3 × 10 21 cm −3 , the exchange term is dominating in the many-body corrections in n-Si.When neglecting the correlation term with the logarithm in (2), the error is < 10%.Therefore, below a simpler correction term due to the exchange only is used, V x = −2n 1/3 × f(n).
The position of the Fermi level μ has been obtained from the condition of the total neutrality [4] [5]: [ ] Here E d is the donor energy with respect to the bottom of the conduction band.
In [8] the axially symmetric case was considered only.The method of solving Equation ( 1) used there cannot be used for a multidimensional case.Here Equation (1) has been solved by the general iteration method [10] [11] [12].The stationary problem has been replaced by a nonstationary one: [ ] ( ) The solving of Equation ( 4) should be realized till establishing stationary solution.Here the iterations have been applied based on the factorization, or Douglas-Rachford, method [11] [12].It is analogous to the classical Gummel's one known in electronics: [ ] where s is the number of the iteration.The iteration algorithm with two fractional steps is used to compute χ and then V (s+1) : .
The iteration scheme Equation (6) approximates Equation ( 4) at all the fractional steps.The rapid convergence of the method has been demonstrated even when the exchange energy has been taken into account.To obtain the accuracy < 0.1%, the number of iterations should be ~50, when the iteration parameter u is u = 0.1 -1.Note that a similar approach based on the Douglas-Rachford, or factorization, method can be applied for general multidimensional structures, whereas in 2D case also the iteration scheme based on the Peaceman-Rachford, or alternate directions, algorithm can be used.While the Peaceman-Rachford method sometimes results in a more rapid convergence, it is not applicable to 3D problems [11] [12].The summarized approximation method is not practically applicable for the establishing method [11] [12] [13].

Calculations of Electron Spectrum
After calculations of the electron potential energy V(x,y) and the exchange energy V x (x,y), the energy levels E j , the wave functions Ψ j (x,y) of the discrete spectrum of the quantum well, and the electron concentration n in each electron level have been computed from the Schrödinger equations [14] [15]: , , , ; There are two different orientations of electron valleys, as seen from Equations (7).
The main attention is paid here to the solving of more complicated Equation (7b).The preference of the proposed method is in the sequential realization of the finding of the electron potential energy and only then the solving of the Schrödinger equations.The second problem can be solved separately for instance by the standard finite element approach using COMSOL Multiphysics [16].But attention should be paid there to an approximation of boundary conditions and to selection of true wave functions.
Below a simple approximate solution method is realized, which is based on the parabolic approximation of the total potential energy with the exchange correction and the expansion of the wave function by the full set of the eigenfunctions of the linear harmonic oscillator.Then the standard quantum mechanical perturbation theory is applied.The discrete Fourier transform is unsatisfactory, because as the result a non-sparse matrix is formed, as our investigations have been demonstrated.Moreover, there is a problem how to select the functions of the discrete spectrum there; those functions should tend to zero at the infinity.
To solve Equation (7b), the following transformations of variables x, y have been applied: As a result, Equation (7b) takes the form: ( ) A solution of Equation ( 9) is searched as a series [14]: An arbitrary complete set of the functions can be used.Here an expansion with the orthonormal Hermite functions φ m (x) is applied [14] [15]: H m (x) are the Hermite polynomials [14].
Below the following notation is used: For lowest energy levels it is rather better to use the values of x c , y c that result in P ≈ 0 near the minimum point x = y = 0, see Figure 1, curve 2.
The approximate expression for the total electron potential energy near the minimum x = y = 0 is, see Figure 1, curve 2, where the cross-section y = 0 is shown: Therefore, to compute the minimum energy levels, the values of x c , y c are suitable: To get higher energy levels more accurately, it is possible to increase the values of x c , y c , to get a better approximation of the electron potential energy specifically near these energy values, see Figure 1, curves 3, 4, where the values of x c are bigger than those in Equation ( 14).Equation ( 9) is equivalent to the following matrix equation: Here 0 Ε > is the electron energy measured from the bottom of the electron potential energy with the exchange correction.The corresponding matrix elements are: The numerical integration has been done by means of the Gauss quadrature formula [10]: When it has been used the renumbering (m,n)  p, Equation ( 15) can be rewritten: ( ) 0 0 0, .
The full spectrum of the quantum well can be obtained from Equation ( 18).
The eigenvalues can be obtained directly or with using the iterations for the inverse matrix [8].Below the quantum mechanical perturbation theory [14] [15] is applied for this purpose.
Because 2D problem is investigated, a possible degeneration of energy levels may occur.Consider a situation when there are several close levels [14] 0 j j j p p p P const Ε + ≈ , the numbers of the energy levels are 1 , , Nj p p  .As usually, a few 2 -5 close energy levels, or resonant ones, take place.Let us rewrite Equation (18) for these resonant levels: For another, or nonresonant, levels Ε 0q Equation ( 18) is simplified: As a result, the following secular equation has been derived: Here the matrix elements are ( ) Thus, to find the eigenvalues and eigenfunctions, the problem of seeking the eigenvalues of the relatively small matrix M pjpk should be solved.Really for higher energy levels it is no more than 6 × 6.
In the nondegenerate case for lower energy levels, the standard perturbation method for the eigenvalue of the energy with the number p 0 is: Also the variation methods can be used to estimate the eigenvalues [8] [12].
They are effective for several lowest energy values.
Thus, at lower temperatures not all the donor levels are ionized.The results of simulations of the total electron potential energy W ≡ V(x,y) + V x (x,y), the exchange energy V x (x,y), and the total electron concentration n(x,y) are presented in Figures 2-4.The energy unit is Ry * ≈ 0.12 eV, the unit for distances x, y is   concentration n does not depend practically on T, see Figures 2-4, left columns, as well as the exchange energy correction, central columns.
After calculating the electron potential energy it is possible to simulate the electron levels of the well and the wave functions from Equation (15).The profiles of three lowest wave functions symmetrical with respect to x and y are presented in Figure 5 and Figure 6.The wave function of the lowest, or fundamental,  the errors is the using of optimal values of approximation parameters x c , y c , see Figure 1.For higher energy levels these values should be 1.3 -2 times bigger than those computed from Equation (14).
A preference of the Thomas-Fermi method is a possible choice of the functions of a preferred symmetry only.In the Schrödinger-Poisson method all the wave functions should be simulated simultaneously.
The dependencies of the corresponding values of the energy on temperature are given in Figure 5.It is seen that there is no dependence on temperature of the difference of the energies of several lowest levels.It is important for practical needs.In the deep well, when the doping level is high, the difference of lowest energy levels does not depend on the temperature within the interval from 20 to 300 K.When considering the higher energy levels, this is valid only approximately.
Also the simulations are presented in Figure 7 for the doping profile different from those in Figures 2-6, Figure 8.The results of the simulations are tolerant to changes of parameters of semiconductor within a wide range of parameters.

Conclusions
The Thomas-Fermi method with many-body corrections has been applied to calculate the electron spectrum in highly doped n-Si quantum wires of an arbitrary doping profile under finite temperatures T.An effective iteration algorithm is proposed for solving the two-dimensional equation for the electron potential energy.This iteration method can be used for arbitrary orientations of the axes of the quantum wire and for other multidimensional structures like quantum dots.The preference of the proposed method is in the sequential realization of

Figure 1 .
Figure 1.Different parabolic approximations of the electron potential energy V(x).Curve 1 is V(x), 2, 3, 4 are parabolic approximations at different values of x c .The curve 2 corresponds to approximation near the minimum of V(x), Equation (14), curves 3, 4 are at bigger values of x c .

Figure 2 .
Figure 2. The dependencies of the total electron concentration n, cm −3 , left columns, the exchange electron energy V x , Ry * units, central columns, and total electron potential energy W, Ry * , right columns, on the coordinates x, y when the exchange is taken into account.Part (a) is at T = 20 K, (b) is at T = 50 K, (c) is at T = 100 K, (d) is at T = 200 K, (e) is at T = 300 K.The maximum doping is N 1d0 = 5 × 10 20 cm −3 , 0 0 5 2.6 nm x a * = ≈ ,

Figure 3 .
Figure 3.The same as in Figure 2, but the maximum doping is N 1d0 = 10 21 cm −3 .

Figure 5 .
Figure 5. lowest wave functions related to Figure 2. Parts (a) are at T = 20 K, parts (b) are at T = 300 K.

Figure 6 .
Figure 6.Three lowest wave functions related to Figure 4. Parts (a) are at T = 20 K, parts (b) are at T = 300 K.

Figure 7 .
Figure 7. Part (a) is the dependencies of the total electron concentration n, cm −3 , left columns, the exchange electron energy V x , Ry * , central columns, and total electron potential energy W, Ry * , right columns, on the coordinates x, y at the temperature T = 300 K. Part (b) is the dependencies of the energy levels on the temperature.The maximum doping is N 1d0 = 3 × 10 21 cm −3 , but another doping profiles: 0