Instruction FACTA UNIVERSITATIS Series: Electronics and Energetics Vol. 28, N o 1, March 2015, pp. 103 - 111 DOI: 10.2298/FUEE1501103G THOMAS-FERMI METHOD FOR COMPUTING THE ELECTRON SPECTRUM AND WAVE FUNCTIONS OF HIGHLY DOPED QUANTUM WIRES IN n-Si  Volodymyr Grimalsky, Outmane Oubram, Svetlana Koshevaya, Christian Castrejon-Martinez CIICAp and the Faculty of Sciences on Chemistry and Engineering, Autonomous University of State Morelos (UAEM), Av. Universidad 1001, 62209 Cuernavaca, Mor., Mexico Abstract. The application of the Thomas-Fermi method to calculate the electron spectrum in quantum wells formed by highly doped n-Si quantum wires is presented under finite temperatures where the many-body effects, like exchange, are taken into account. The electron potential energy is calculated initially from a single equation. Then the electron energy sub-levels and the wave functions within the potential well are simulated from the Schrödinger equation. For axially symmetric wave functions the shooting method has been used. Two methods have been applied to solve the Schrödinger equation in the case of the anisotropic effective electron mass, the variation method and the iteration procedure for the eigenvectors of the Hamiltonian matrix. Key words: quantum wires, Thomas-Fermi method, exchange, wave functions, variation method, inverse matrix iterations 1. INTRODUCTION Investigations of the electron spectrum of low-dimensional and highly doped structures are central to many nanotechnology applications [1,2]. Quantum devices based on silicon have been the subject of a concentrated recent interest, both experimental and theoretical, including the recent proposals on quantum computing [1-4]. The infrared transitions between the electron sub-levels within -doped quantum wells are perspective for using in optoelectronics, especially for infrared modulators, detectors, and lasers [5,6]. The electron spectrum of -doped quantum wells can be calculated from solving Schrödinger equation jointly with the Poisson one (SP) [5,6]. There exist several difficulties for simulations of quantum structures in silicon, namely, anisotropy of the Received July 8, 2014; received in revised form November 22, 2014 Corresponding author: Volodymyr Grimalsky CIICAp and the Faculty of Sciences on Chemistry and Engineering, Autonomous University of State Morelos (UAEM), Av. Universidad 1001, 62209 Cuernavaca, Mor., Mexico (e-mail: v_grim@yahoo.com) 104 V. GRIMALSKY, O. OUBRAM, S. KOSHEVAYA, C. CASTREJON-MARTINEZ effective electron mass and slow convergence of SP method in the case of an arbitrary initial approximation. The investigation of -doped quantum structures is possible with a simpler approach based on the statistical Thomas-Fermi (TF) method [6-8]. The preference is the separation of the complex problem into the sequential ones, where the wave functions are computed after the found solution of the potential energy. In the case of axially symmetric high doping, the electron potential energy depends of the radius only. Also the combined method can be applied, where the final result of TF simulations is used as a starting one for SP [9]. Moreover, a comparison between SP method and TF one shows that the simple TF method gives a good approximation for the electron energy sub-levels and the total electron concentration within the -doped quantum wells [9]. In this paper the application of the TF method to calculate the electron spectrum in the quantum wells formed by highly doped n-Si quantum wires is presented under finite temperatures T, and many-body effects, like exchange, are taken into account [5,7]. The electron potential energy and the total electron concentration are calculated from a single equation solved by the Newton method. Then the wave functions and values of the electron energy sub-levels are computed from the Schrödinger equation where two possible orientations of electron valleys are considered. The peculiarities of solving Schrödinger equation in the case of the anisotropic electron mass are pointed out. 2. BASIC EQUATIONS Consider a single -doped electron quantum wire within n-Si. Below the atomic units are used for distances a0 * =2/(mce 2 ) 0.52 nm and for energy Ry * = e 2 /(2a0 * )  0.12 eV, where mc =  2/3 (m 2 m||) 1/3  1.06 me  10 -27 g,  = 6 is the number of the lowest electron valleys in Si. In the case of n-Si the lowest valleys are lateral and the effective mass is anisotropic: m||, m. With non-dimensional variables the basic equation of TF method for the -doped electron quantum wire is [8]: 1 0 1 1 ( ) 8 { [ ] (2 exp( ) 1) ( )}; ( 0) 0, ( ) 0; d d d E Vd dV r n V N N r r dr dr T dV r V r dr                 (1) where 3 / 2 1/ 2 1/ 3 1/ 2 1/ 2 1/ 2 0 2 17 3 1 1 0 0 1 [ ] ; 4 tanh(( / ) 1), ; 2 ( ); ( ) 0, ; 2 ( ) ; 1 exp( ) ( ) exp( ( / ) ); 10 . x c c x c d d c V VT n V T n n n n V n f n f n n n u du v u v N r N r r n cm                                  (2) Here V(r) is the electrostatic electron energy, n is the total electron concentration; N1d and Nd0 are 1D and 3D donor concentrations, respectively. Vx is the many-body correction Thomas-Fermi Method for Computing the Electron Spectrum snd Wave Functions... 105 to the electron energy due to the exchange [7]. Eq. (1) is the Poisson equation for the electrostatic electron energy, where the electron concentration n[V] is calculated from the equilibrium statistical Fermi distribution. Note that at the 1D donor concentrations N1d0  10 21 cm -3 the exchange energy Vx is comparable with the electrostatic one. The donor levels are assumed shallow and single charged. The concentration of 1D donors is high N1d0  10 20 cm -3 ; they are fully ionized. The 1D doping is localized at the distances r ~ r0  1 – 5 nm. The finite size of the highly doping region r0 is considered, because the distance unit a0 * in silicon is comparable with the size of the lattice cell. Moreover, 1D doping cannot be approximated by the -function directly, because this approximation leads to the logarithmic singularity of the electron potential energy at r ~ 0. The results of simulations do not depend on the value of the critical electron concentration when nc  10 18 cm -3 . The position of the Fermi level  has been obtained from the condition of the total neutrality of the semiconductor [6]: 1 0 [ 0] (2 exp( ) 1)d d E n V N T       (3) Here Ed is the donor energy with respect to the bottom of the conduction band. Eq. (1) has been solved by the Newton method [8]. .0|,0| )};()(][{8)( 1 )(8)( 1 ; 0 1 1            Rrr d s d s s dss dr d rNVNVn dr dV r dr d r Q V N V n dr d r dr d r VV       (4) Note that in the derivative (n/V) the exchange correction Vx does not vary. In the boundary conditions the parameter R is an enough big radius. In Eq. (4) the parameter Q  1 is chosen to provide better convergence [9]. The rapid convergence of the method has been demonstrated, even when the exchange energy has been taken into account. After calculation of the electron potential energy V(r), the energy sub-levels Ej, the wave functions (WF) j(x,y) of the discrete spectrum of the well, and then the electron concentration n in each electron sub-level within the quantum well have been computed from the following Schrödinger equations: ;])([)(ˆ )1()1()1( 2 )1(2 2 )1(2 )1( 1 jjjx jjc j EVrV yxm m H         (5a) ;])([ˆ )2()2()2( 2 )2(2 || 2 )2(2 )2( 2 jjjx jcjc j EVrV ym m xm m H         (5b) .)(;1)( 2/1222 yxrdxdy j       WF can be chosen as real, because the Hamiltonians are real. There are two different orientations of electron valleys in silicon, as seen from Eqs. (5). Namely, Eq. (5a) is for the isotropic case of the effective mass components, Eq. (5b) is for the anisotropic case. 106 V. GRIMALSKY, O. OUBRAM, S. KOSHEVAYA, C. CASTREJON-MARTINEZ 3. SIMULATIONS OF POTENTIAL ENERGY AND ELECTRON CONCENTRATION The results of simulations of the electron potential energy V(r) and the electron concentration n(r) are presented in Figs. 1,2. For all cases the volume doping is Nd0 = 1·10 16 cm -3 , r0 = 2a0 *  1.04 nm. The previous simulations [10] demonstrated that the exchange correction is important for the doping levels N1d  10 21 cm -3 . But the total electron potential energy W = V + Vx and the electron concentration n are practically the same as without this many-body correction. The potential energy depends on temperature T, as seen in Fig. 2. This is due to the partial ionization of volume donors Nd0 at low temperatures, as seen from Eq. (3). But the electron concentration n does not depend on T. Some difference is only at the periphery r >> r0. a) b) Fig. 1 Part a) is dependence of electron potential energy jointly with the exchange energy V+Vx on the radius r. Part b) is the dependence of the total electron concentration n(r). The values of the maximum doping are N1d0 = 310 21 cm -3 (curve 1), 10 21 cm -3 (curve 2), 310 20 cm -3 (curve 3), 510 19 cm -3 (curve 4), T = 300 K, r0 = 2a0 *  1.04 nm. The corresponding exchange energies Vx are also presented there in the upper part of the part a). a) b) c) Fig. 2 Part a) is dependence of electron potential energy jointly with the exchange energy V+Vx on the radius r for different temperatures T. Part b) is the dependence of the electron concentration for different temperatures; part c) is the same as b) in details. Curve 1 is for T = 300 K, 2 is for 200 K, 3 is for 150 K, 4 is for 100 K, 5 is for 50 K, 6 is for 20 K. The maximum doping is N1d0 = 310 21 cm -3 . Thomas-Fermi Method for Computing the Electron Spectrum snd Wave Functions... 107 4. WAVE FUNCTIONS AND ENERGY SUB-LEVELS After calculating the electron potential energy it is possible to simulate the electron energy sub-levels in the quantum well and the corresponding WF. To compute WF for the isotropic case (5a), where the effective masses are the same, the shooting method has been applied [11]. The axially symmetric WF (r) for the isotropic case, Eq. (5a), are presented in Fig. 3, a. The maximum doping level is N1d0 = 310 21 cm -3 , T = 300 K, as in Fig. 1, a, curve 1. The dependence of the electron energy sub-levels, two lowest ones E1,2, on the maximum doping is presented in Fig. 3, b, for T = 300 K. The dependence of the electron energy sub-levels on temperature is given in Fig. 3, c. One can see that the difference E2 – E1 depends on the temperature T there. a) b) c) Fig. 3 Part a) is the axially symmetric WF for the case of isotropic effective mass; part b) is the dependence of the energy sub-levels on the maximum doping, T = 300 K; c) is the dependence of the energy sub-levels on the temperature T for the maximum doping N1d0 = 310 21 cm -3 . The anisotropic case, Eq. (5b), with different effective masses has been solved by two simple methods, which are realized in the Cartesian coordinate frame XOY. The first one is the variation method [12]. Namely, the problem of the minimization of the functional of the electron energy is considered:                              dxdy dxdyH E 2 2 )( ˆ min (6) WF possesses different types of symmetry or antisymmetry in the plane XOY, due to the symmetry of the Hamiltonian, Eq. (5b). The probing functions for the symmetric case (±x, ±y)=(x,y) are chosen as: 2 2 1 01 01 2 2 2 2 2 02 02 02 02 exp( ( ) ( ) ); exp( ( ) ( ) )(1 ( ) ( ) ); x y x y x y x y a b x y x y           (7) 108 V. GRIMALSKY, O. OUBRAM, S. KOSHEVAYA, C. CASTREJON-MARTINEZ In the case of antisymmetry with respect to x (-x, ±y) = -(x,y) the probing functions are: );)()(1)()()(exp( );)()(exp( 2 02 2 02 2 02 2 02 2 2 01 2 01 1 y y b x x a y y x x x y y x x x   (8) Analogously, it is possible to write down the probing functions for other types of symmetry or antisymmetry, i.e. with the multipliers y or xy. Therefore, for the lowest WF there are two variation parameters x01, y01. For the second WF there are 3 independent variation parameters, because of the imposing orthogonality relation:        0 21 dxdy . (9) The second method is the search of the eigenvalues of the matrix of the Hamiltonian by means of the iteration procedure [13,14]. For this purpose, WF has been expanded by the truncated Fourier series. Zero boundary conditions have been used at the periphery x = ±Lx, y = ±Ly, where the boundaries Lx, Ly are chosen enough large. Namely, WF is represented by the vector, or the column of the coefficients of the expansion; the Hamiltonian has been represented by the matrix. Then the following iteration procedure has been applied [13,14]: ),( ),( ;)ˆ( 10 1 1 02 1        ss ss s ss EE EH (10) Here ),( 1  ss is the scalar product of vectors, s in the number of iterations, E0 < 0 is the parameter that has been chosen from the condition of maximum convergence. Usually, E0 is close to the lowest energy sub-level computed from the variation method. It is important that the matrix inversion can be realized in the simple manner, because of the diagonal domination of the shifted Hamiltonian matrix )ˆ( 02 EH  . When the second WF is searched, it should be orthogonal to the first WF 1: ),/(),( 111 1 1 11   sss . After each iteration it is better to normalize the vector: .1),(  ss The initial values of the vectors  s=0 can be chosen as ones found earlier from the variation method. The iterations with the direct Hamiltonian matrix )ˆ( 02 EH  diverge and cannot be applied. The profiles of WF for the two lowest sub-levels are presented in Fig. 4 for the temperature T = 300 K and the maximum doping N1d0 = 310 21 cm -3 . The dependencies of the two lowest energy sub-levels in the quantum well on the maximum doping concentration for T = 300 K and on the temperature T for the maximum doping N1d0 = 310 21 cm -3 are given in Fig. 5 for symmetric WF. Thomas-Fermi Method for Computing the Electron Spectrum snd Wave Functions... 109 a) b) c) d) e) f) g) h) Fig. 4 Wave functions computed for the case of anisotropic effective masses, T = 300 K, N1d0 = 310 21 cm -3 . Part a) is the first symmetric wave function; the left panel is computed from the matrix iteration method, the right panel is from the variation method. Part b) is the same for the second symmetric wave function. Parts c) and d), e) and f), g) and h) are the first and the second wave functions correspondingly computed from the variation method for different types of symmetry or antisymmetry. 110 V. GRIMALSKY, O. OUBRAM, S. KOSHEVAYA, C. CASTREJON-MARTINEZ a) b) Fig. 5 The dependence of two lowest energy sub-levels in the quantum well on the maximum doping concentration for T = 300 K (a) and on the temperature for N1d0 = 310 21 cm -3 (b). The case of anisotropic effective mass, Eq. (5b), symmetric WF, is considered. The solid lines are the data obtained from the matrix iteration method, the dot lines are ones from the variation method. The variation method yields accurate values of the energy sub-levels. For instance, at T = 300 K and N1d0 = 310 21 cm -3 the values of the energy for the symmetric case calculated from the iteration procedure are E1 = -4.09 Ry * , E2 = -1.94 Ry * . The same values computed from the variation method are E1 = -4.075 Ry * , E2 = -1.895 Ry * . The profiles of WF are calculated from the variation method approximately; there is some difference at the periphery from those computed from the matrix iteration procedure. In the report [10] the electron spectrum has been calculated from the shooting method applied in the polar coordinate system. There is coincidence of the energy sub-levels with the data presented here, but that numerical realization of the shooting method is more complicated. The difference of the lowest energetic sub-levels E2 – E1 does not depend on the temperature for the anisotropic case, symmetric WF, see Fig. 5, b. For the isotropic case, Eq. (5a), this is not valid, see Fig. 3, c. This can be explained by higher values of the electron sub-levels |E1,2| for the anisotropic case. It is possible to calculate WF more accurately also by means the standard simulators based on finite element methods, like COMSOL Multiphysics [15]. 5. CONCLUSIONS An application of TF method to the electron spectrum of quantum wires in n-Si can be subdivided into two problems. The first one is the simulation of the electron potential energy from the simple ordinary differential equation. The iteration procedure demonstrates rapid convergence even when the many-body effects, like exchange, are taken into account. Then it is possible to solve the Schrödinger equations for the wave functions and the energy sub-levels. Because of the anisotropy of the effective electron mass in silicon, this problem is generally two-dimensional. Two simple methods have been proposed. The variation method yields accurate values of the energy sub-levels, whereas the profiles of the electron wave functions are approximate at the periphery. The method based on the inverse matrix iterations is more accurate both for the eigenvalues and the eigenfunctions. Thomas-Fermi Method for Computing the Electron Spectrum snd Wave Functions... 111 Acknowledgement: The authors would like to thank to SEP-CONACyT (Mexico) for a partial support of our work. REFERENCES [1] D. W. Drumm, A. Budi, M. C. Per, S.P. Russo, and L.C. L. Hollenberg. “Ab initio calculation of valley splitting in monolayer δ-doped phosphorus in silicon”, Nanoscale Research Lett., vol. 8, no 1, pp. 111- 121, Jan. 2013. [2] B. Weber, S. Mahapatra, H. Ryu, S. Lee, A. Fuhrer, T. C. G. Reusch, D. L. Thompson, W. C. T. Lee, G. Klimeck, L. C. L. Hollenberg, and M. Y. Simmons, “Ohm’s law survives to the atomic scale”, Science, vol. 335, no 6064, pp. 64-67, Jan. 2012. [3] F. J. Ruess, W. Pok, T. C. G. Reusch, M. J. Butcher, Kuan Eng J. Goh, L. Oberbeck, G. Scappucci, A. R. Hamilton, and M. Y. Simmons, “Realization of atomically controlled dopant devices in silicon”, Small, vol. 3, pp. 563 – 567, Apr. 2007. [4] D.K. Ferry, S.M. Goodnick, and Jonathan Bird, Transport in Nanostructures, Cambridge: Cambridge Univ. Press, 2009. [5] L. Ramdas Ram-Mohan, Finite Element and Boundary Element Applications in Quantum Mechanics, Oxford: Oxford University Press, 2002. [6] Y. Fu and M. Willander, Physical Models of Semiconductor Quantum Devices, Dordrecht: Kluwer, 1999. [7] L. M. Gaggero-Sager, “Exchange and correlation via functional of Thomas-Fermi in delta-doped quantum wells”, Modelling Simul. Mater. Sci. Eng., vol. 9, pp. 1-5, Jan. 2001. [8] V. Grimalsky, L. M.Gaggero-S., S.Koshevaya, and A.Garcia-B., “Electron spectrum of -doped quantum wells by Thomas – Fermi method at finite temperatures”, in Proc. 27 th International Conference on Microelectronics (MIEL 2010), Nis, Serbia, 2010, pp. 119-122. [9] C. Castrejon-Martinez, V. Grimalsky, L. M. Gaggero-Sager, and S. Koshevaya, “Combined method for simulating electron spectrum of delta-doped quantum wells in n-Si with many-body corrections”, Progress In Electromagnetics Research M, vol. 31, pp. 215-229, Aug. 2013. [10] V. Grimalsky, O. Oubram, S. Koshevaya, and C. Castrejon-M., “Thomas-Fermi method for computing the electron spectrum of highly doped quantum wires in n-Si”, In Proc. 29 th International Conference on Microelectronics (MIEL 2014), Belgrade, Serbia, 2014, Paper 049. [11] V. A. Ilyina and P. K. Silaev, Numerical Methods for Theoretical Physicists, vol. 2, Moscow: Institute for Computing Research Publ., 2004 (in Russ.). [12] K. T. Hecht, Quantum Mechanics, Springer, N.Y., 2000. [13] W. H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, Numerical Recipes in Fortran, Cambridge University Press, Cambridge, 1997. [14] J. Stoer and R. Burlish, Introduction to Numerical Analysis. N.Y.: Springer, 2002. [15] S. M. Musa, ed., Computational Finite Element Methods in Nanotechnology. Boca Raton, CA: CRC Press, 2013 (www.comsol.com). https://www.sciencemag.org/search?author1=B.+Weber&sortspec=date&submit=Submit https://www.sciencemag.org/search?author1=S.+Mahapatra&sortspec=date&submit=Submit https://www.sciencemag.org/search?author1=H.+Ryu&sortspec=date&submit=Submit https://www.sciencemag.org/search?author1=S.+Lee&sortspec=date&submit=Submit https://archive.today/o/Zfl0/http:/www.sciencemag.org/search?author1=A.+Fuhrer&sortspec=date&submit=Submit https://archive.today/o/Zfl0/http:/www.sciencemag.org/search?author1=T.+C.+G.+Reusch&sortspec=date&submit=Submit https://archive.today/o/Zfl0/http:/www.sciencemag.org/search?author1=D.+L.+Thompson&sortspec=date&submit=Submit https://archive.today/o/Zfl0/http:/www.sciencemag.org/search?author1=W.+C.+T.+Lee&sortspec=date&submit=Submit https://archive.today/o/Zfl0/http:/www.sciencemag.org/search?author1=G.+Klimeck&sortspec=date&submit=Submit https://archive.today/o/Zfl0/http:/www.sciencemag.org/search?author1=G.+Klimeck&sortspec=date&submit=Submit https://archive.today/o/Zfl0/http:/www.sciencemag.org/search?author1=L.+C.+L.+Hollenberg&sortspec=date&submit=Submit https://archive.today/o/Zfl0/http:/www.sciencemag.org/search?author1=M.+Y.+Simmons&sortspec=date&submit=Submit http://www.comsol.com/