Solution of the 1D Schrödinger equation in semiconductor heterostructures using the immersed interface method

Journal of Computational Physics 231 (2012) 6173–6180

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

Solution of the 1D Schrödinger equation in semiconductor heterostructures using the immersed interface method E. Momox ⇑, N. Zakhleniuk, N. Balkan School of Computer Science and Electronic Engineering, University of Essex, Wivenhoe Park, Colchester CO4 3SQ, United Kingdom

a r t i c l e

i n f o

Article history: Received 12 January 2011 Received in revised form 9 September 2011 Accepted 13 May 2012 Available online 26 May 2012 Keywords: Schrödinger equation Immersed interface method Finite difference method Semiconductors Potential profiles

a b s t r a c t Due to the enormous progress in computer technology and numerical methods that has been achieved over the past several decades, the use of numerical simulation methods in all scientific disciplines gain more and more importance. In the physics field, these methods have provided remarkable numerical solutions to problems considered analytically intractable. The solution of the Schrödinger equation in semiconductor heterostructures is a good example. However, many of these numerical schemes are cumbersome to implement for nonexperts in numerical computing. With this reason as motivation, a novel method simple enough to implement yet powerful enough to solve Schrödinger equation in semiconductor devices with high accuracy is presented herein. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction In a large number of problems encountered in physics, and specifically in the modelling of semiconductor heterostructures, a fast and accurate solution of the one-dimensional time-independent Schrödinger equation has to be found in order to better understand the physical properties of such devices and to simulate their expected performance. The abovementioned is an effect of the characteristic dimensions of up-to-date devices that have progressed further than the submicron regime where quantum effects play a paramount role. However, analytic solutions to these types of problems are rarely encountered or even nonexistent. For this reason, several methods have been devoted to the numerical solution of such issues, including the transmission line method [1], spectral methods [2–4], the so called shooting method [5–8], matrix transfer method [9–11], finite element method [12–14], and of more immediate interest, the finite difference method [11,15,16]. Many of these techniques are cumbersome to implement thus knowledge in numerical computing might be required which constitutes an obstacle for nonexperts. As demonstrated in the forthcoming sections, the presented method has been readily and successfully implemented for the solution of eigenvalue problems in arbitrary potential profiles originated in semiconductor devices. The foundation of the herein presented algorithm is the immersed interface method (IIM) that has been developed in recent years predominantly designed for interface problems and which modifies the numerical schemes in the neighbourhood of or on the interfaces [17]. Its application in solving the Schrödinger equation in semiconductor heterostructures is described in Section 2 whereas some numerical examples appear in Section 3 to exemplify the generality and ease for computing the eigenstates of such profiles when using this novel numerical approach.

⇑ Corresponding author. E-mail address: [email protected] (E. Momox). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.05.017

6174

E. Momox et al. / Journal of Computational Physics 231 (2012) 6173–6180

2. Derivations For the analysis, we consider a quantum well (QW) of arbitrary potential profile. Under the effective mass approximation, the envelope function of an electron or a hole in the QW, is given by the one-dimensional time independent Schrödinger equation



! 2  d 1 d h þ UðzÞ WðzÞ ¼ EWðzÞ; 2 dz me ðzÞ dz

ð1Þ

where E is the energy of the particle, U(z) is the potential, me(z) is the effective mass of the particle and 2p⁄ is Planck’s constant. We will use normalized parameters to facilitate the derivation of the following equations. Hence, the energy E, the potential U(z), and the coordinate z, are normalized as e = E/E0, u(f) = U(z)/E0, f = z/d, respectively, where E0 = ⁄2 p2/ 2me(0)d2, me(0) is the effective mass of the particle at z = 0, and d the width of the potential U(z). With this notation, the Schrödinger equation, Eq. (1), can be written in the following form

  1 d 1 d þ uðfÞ  2 WðfÞ ¼ eWðfÞ; p df m ðfÞ df

ð2Þ

where m⁄(f) = me(z)/me(0). If we consider the effective mass of the particle to be a piecewise constant, then it can be taken outside the derivative. Therefore, Eq. (2) has the form

! 2 1 d  2  þ uðfÞ WðfÞ ¼ eWðfÞ p m ðfÞ df2

ð2aÞ

to further simplify Eq. (2a), we can define M ¼ MðfÞ ¼ p2 m1 ðfÞ so that Eq. (2a) can be rewritten as

M Wff þ uðfÞW ¼ eW;

ð3Þ

where the subscripts denote differentiation. In the finite difference method (FDM), the range of integration (a, b) is divided into N  1 equal subintervals, grid points, of length h represented as fi = a + (i  1)h, i = 1, 2, 3, . . . , N and h = (b  a)/(N  1). On the other hand, the piecewise constant function M is allowed to have discontinuities at some grid points. For the analysis, we can consider one of these discontinuities, i.e. f = f⁄ where the potential u(f) may also have an abrupt change at such a point. The point f⁄ might be an element of the set of grid points or will have two of them as its neighbours, i.e., fj 6 f⁄ < fj+1, as shown in Fig. 1. In order to solve Eq. (3) under the FDM, it is convenient to replace the second order derivative by a finite difference approximation such as the central divided difference scheme. This scheme will result in an eigenvalue problem with a large tridiagonal matrix which can then be solved with the use of conventional numerical libraries such as ARPACK [18]. A similar procedure is found in [16]. However, the direct use of the final equations reported in [16] leads to far-reaching numerical difficulties. The correct form of such equations is of paramount importance for its numerical implementation. By the same token, we derive here their correct final form. If the definition of the second order central finite difference approximation is substituted in Eq. (3), we have

 Mðfi Þ

1

h

2

 ½wðfiþ1 Þ  2wðfi Þ þ wðfi1 Þ þ uðfi Þwðfi Þ ¼ ewðfi Þ;

rearranging terms in Eq. (4) we get

Fig. 1. Arbitrary potential profile, continuous line, originated in a single heterojunction and grid points close to the interface or discontinuity.

ð4Þ

E. Momox et al. / Journal of Computational Physics 231 (2012) 6173–6180



Mðfi Þ h

2

 wðfiþ1 Þ þ

2Mðfi Þ h

2

 Mðf Þ þ uðfi Þ wðfi Þ  2 i wðfi1 Þ ¼ ewðfi Þ h

6175

ð4aÞ

or

b1 ðiÞwðfi1 Þ þ b2 ðiÞwðfi Þ þ b3 ðiÞwðfiþ1 Þ ¼ ewðfi Þ;

ð4bÞ

where we have defined

b1 ðiÞ ¼ b3 ðiÞ ¼ 

Mðfi Þ 2

h

and b2 ðiÞ ¼ 2b1 ðiÞ þ uðfi Þ:

ð5Þ

Eq. (4b) represents the problem we must solve. It has the following form, A w = ew, where A is a sparse tridiagonal matrix with dimensions (N  2)  (N  2), w the eigenvectors, and e the corresponding eigenvalues, its representation in matrix notation is,

10 1 1 0 b2 ð2Þ b3 ð2Þ 0 w2 w2 CB w C B b ð3Þ b ð3Þ b ð3Þ B w C 2 3 CB 3 C B 1 B 3 C CB C C B B CB w4 C B B w4 C b1 ð4Þ b2 ð4Þ b3 ð4Þ CB C C B B ¼ e CB . C B B . C: .. .. .. C C C B B B . . . . . CB . C B B . C CB C C B B @ @ wN2 A b1 ðN  2Þ b2 ðN  2Þ b3 ðN  2Þ A@ wN2 A wN1 wN1 0 b1 ðN  1Þ b2 ðN  1Þ 0

ð6Þ

Note that we have set w(a) and w(b) equal to zero as the wavefunction must be zero at ±1, a condition that computationally is not feasible. Clearly, equating the wavefunction to zero at finite distance points is comparable to inserting an infinite barrier (infinite potential) at those points. Moreover, the assumption that these infinite barriers will not perturb significantly the energy structure can be made. However, the discontinuities in the mass and potential must be taken into account for a valid solution of the eigenvalue problem, failing to do this will likely lead to loss of accuracy in the eigenstate calculation. By the same token, we introduce here the IIM to locally modify the scheme near or on the interface, f = f⁄. The grid points which are close to such a point are fi=j, and fi=j+1, and are called irregulars whereas the majority are called regulars. Consequently, it is necessary to determine the new values for the bs at all irregular points. Therefore, we use a Taylor series expansion about the point f⁄ for w(fi) = wi, where i = j  1, j, j + 1, and j + 2, i.e.,

1 wj1 ¼ w þ ðfj1  f Þwf þ ðfj1  f Þ2 wff ; 2 1 wj ¼ w þ ðfj  f Þwf þ ðfj  f Þ2 wff ; 2 1 þ  þ wjþ1 ¼ w þ ðfjþ1  f Þwf þ ðfjþ1  f Þ2 wþff ; 2 1 þ  þ wjþ2 ¼ w þ ðfjþ2  f Þwf þ ðfjþ2  f Þ2 wþff ; 2

ð7Þ

ð7aÞ ð7bÞ ð7cÞ

Notice the () and (+) superscripts on the functions defined at the left and right of the discontinuity, f = f⁄, respectively. On þ the other hand, Eqs. (7b) and (7c) involve the values of wþ ; wþ f , and wff . For the foregoing reason, relationships to connect the +  ( )-region to the ( )-region are required. We call these relationships, internal boundary conditions. At the interface f = f⁄, particle conservation requires the continuity of the wavefunction, i.e.,

wþ ¼ w ¼ w :

ð8Þ

Bastard [19] introduced another internal boundary condition which can be obtained assuming Eq. (8) and integrating Eq. (3), þ þ that is to say M w f ¼ M wf , or

wþf ¼

M  w : Mþ f

ð8aÞ

Eq. (8a) is the only internal boundary condition which ensures stationary eigenstates even if there is a jump in the effective mass [19]. The remaining internal boundary condition can be derived from Eqs. (3) and (8), namely,

wþff ¼ wff

M uþ  u þw : Mþ Mþ

ð8bÞ

The left hand side of Eq. (4b) can be written at fj , using Eq. (7), to get the next expression,

    1 1 b1 w þ ðfj1  f Þwf þ ðfj1  f Þ2 wff þ b2 w þ ðfj  f Þwf þ ðfj  f Þ2 wff 2 2   1 þ b3 wþ þ ðfjþ1  f Þwþf þ ðfjþ1  f Þ2 wþff 2

ð9Þ

6176

E. Momox et al. / Journal of Computational Physics 231 (2012) 6173–6180

or rearranging terms,

   b b w ½b1 þ b2  þ wþ b3 þ wf b1 ðfj1  f Þ þ b2 ðfj  f Þ þ wþf b3 ðfjþ1  f Þ þ wff 1 ðfj1  f Þ2 þ 2 ðfj  f Þ2 2 2 b þ wþff 3 ðfjþ1  f Þ2 : 2

ð9aÞ

Taking into account Eq. (8), the internal boundary conditions, the following expression has been derived,

   þ    1 u  u M þ wf b1 ðfj1  f Þ þ b2 ðfj  f Þ þ þ b3 ðfjþ1  f Þ w b1 þ b2 þ b3 1 þ ðfjþ1  f Þ2 þ 2 M M    b b M b þ wff 1 ðfj1  f Þ2 þ 2 ðfj  f Þ2 þ þ 3 ðfjþ1  f Þ2 : 2 2 M 2

ð9bÞ

Matching coefficients with Eq. (3), we get a 3  3 linear system of equations at fj,

  þ  1 u  u ¼ u ; b1 þ b2 þ b3 1 þ ðfjþ1  f Þ2 þ 2 M M b1 ðfj1  f Þ þ b2 ðfj  f Þ þ b3 ðfjþ1  f Þ þ ¼ 0; M b1 b2 b3 M  2  2 ðfj1  f Þ þ ðfj  f Þ þ ðfjþ1  f Þ2 þ ¼ M  : 2 2 2 M

ð10Þ

Proceeding in the same manner but rewriting the internal boundary conditions, the system of equations at fj+1 has the following form,

  1 u  uþ þ b2 þ b3 ¼ uþ ; b1 1 þ ðfj  f Þ2  2 M Mþ þ b2 ðfjþ1  f Þ þ b3 ðfjþ2  f Þ ¼ 0; M 1 Mþ 1 1 b1 ðfj  f Þ2  þ b2 ðfjþ1  f Þ2 þ b3 ðfjþ2  f Þ2 ¼ M þ : 2 2 2 M

b1 ðfj  f Þ

ð11Þ

The linear systems of equations, Eqs. (10) and (11), represent the values that the b coefficients must have before and after the interface or discontinuity point f⁄. In other words, for a problem with N interfaces, such as the potential profile originated in multiple quantum well devices, the two linear systems must be solved N-times and their solutions conveniently substituted in Eq. (6) for the correct solution of the eigenvalue problem. Next, we consider various potential profiles which appear in semiconductor heterostructures and employ the herein described method for the appraisal of the energy levels and wavefunctions. Note that in Figs. 2 and 5 the wavefunctions were shifted with respect to their energy in order to be properly plotted onto the potential profile, which is the most common layout found in the relevant literature. 3. Applications 3.1. Infinite parabolic quantum well As a first application of the presented algorithm, we consider a parabolic quantum well made of Gallium arsenide (GaAs). An infinite parabolic potential represents a harmonic oscillator commonly encountered in textbooks on quantum mechanics. Its eigenvalues are given by En = (n  1/2)⁄x0 where n = 1, 2, 3, . . .. On the other hand, the energy levels are equally spaced by ⁄x0 and the eigenfunctions exhibit even–odd parity behaviour. Computationally, it makes no sense to have a potential profile with infinite height. However, it is possible to emulate the infinite parabolic GaAs QW by setting the amplitude to a large value, e.g. 80 eV. In computing ⁄x0 = 5.39572 eV and Fig. 2, the following parameters were considered; DVCB = 0.67, DEg = 1.247 eV and m = 0.067 m0. Fig. 2 shows the 15 bounded energy levels and their squared wavefunctions. In Table 1, various solutions of the infinite parabolic well were computed against the number of grid points employed for the calculation of the eigenvalues. Likewise, Fig. 3 shows the absolute error for the lowest confined states as a function of the number of grid points. The lowest computed energy levels are in excellent agreement with the theoretical ones whereas the remaining deviate considerably due to the influence of the finiteness of the parabolic profile which indicates that intermediate levels ‘‘see’’ an infinite potential so they are evenly separated. In this first application, no discontinuities are present but the solution is in excellent agreement with the theory. Put another way, the presented method becomes the standard FDM if no discontinuities are present which demonstrates the generality of the herein described algorithm. It can be seen from Fig. 3 that about 104 grid points are enough to compute a solution with satisfactory accuracy.

6177

E. Momox et al. / Journal of Computational Physics 231 (2012) 6173–6180

Fig. 2. Infinite Parabolic Quantum Well. Energy levels and squared wavefunctions. The considered semiconductor material is GaAs. The wavefunctions were shifted with respect to their energy level.

Table 1 Comparison of the computed confined energy levels as a function of the number of grid points in the parabolic quantum well against the analytical solution. n

Exact solution ⁄x(n  1/2) [eV]

Computed solution grid points = 1  103

Computed solution grid points = 8  103

Computed solution grid points = 64  103

Computed solution grid points = 128  103

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

2.6978598 8.0935793 13.4892988 18.8850183 24.2807378 29.6764573 35.0721768 40.4678963 45.8636158 51.2593353 56.6550548 62.0507743 67.4464938 72.8422133 78.2379329

2.7058625 8.1172244 13.5278604 18.9377701 24.3469536 29.7554104 35.1631386 40.5701265 45.9763072 51.3813540 56.7838822 62.1787050 67.5480892 72.8328005 77.7981384

2.6992647 8.0977883 13.4963007 18.8948019 24.2932917 29.6917701 35.0902355 40.4886771 45.8870323 51.2849914 56.6812275 62.0707283 67.4362544 72.7201024 77.6955594

2.6980441 8.0941322 13.4902201 18.8863079 24.2823954 29.6784826 35.0745679 40.4706406 45.8666391 51.2622570 56.6561810 62.0434448 67.4069553 72.6894969 77.6667844

2.6980020 8.0940059 13.4900098 18.8860136 24.2820173 29.6780209 35.0740226 40.4700120 45.8659272 51.2614621 56.6553038 62.0424873 67.4059246 72.6884169 77.6657656

Fig. 3. Absolute error for each of the lowest states in the infinite parabolic quantum well as a function of the number of grid points.

6178

E. Momox et al. / Journal of Computational Physics 231 (2012) 6173–6180

3.2. Double asymmetric quantum well Furthermore, the implementation of the presented algorithm can be done so as to have a function for its further use as a subfunction. With this line of reasoning, we consider now a double asymmetric QW where the length of the larger well is varied from 7–100 nm, see Fig. 4; whereas the width of the first well is fixed to 7 nm. The barriers are made of Al0.3Ga0.7As and the wells of GaAs. On the other hand, both GaAs quantum wells are separated by a 5 nm slab of Al0.3Ga0.7As. This structure was proposed by Ferreira et al. [20] and Fig. 4a illustrates the energy level diagram of the asymmetrical quantum well structure. Clearly, whenever an eigenstate of the larger well matches the energy of the ground state ( 47 meV), mainly localized in the 7 nm thick well, an anticrossing due to bonding and antibonding states occurs. Fig. 4a is in excellent agreement with the one obtained in the aforesaid reference.

3.3. Quantum Cascade Laser (QCL) As a final application, we consider a QCL structure reported by Page et al. [21]. They proposed such a structure to improve the temperature dependence of the devices emitting at k  9 lm. The structure, GaAs/Al0.45Ga0.55As, has got the following layer sequence starting from the injection barrier 46, 19, 11, 54, 11, 48, 28, 34, 17, 30, 18, 28, 20, 30, 26, 30 Å. The normal scripts are GaAs wells and bold scripts Al0.45Ga0.55As barriers. The QCL is under an electric field of 48 kV/cm. Our simulation, Fig. 5, is in good agreement with Fig. 1 depicted in the abovementioned reference. It also displays, in bold without marker, the squared wavefunctions of the ground, lower, excited, and injector states.

Fig. 4. Illustration of the potential profile of a double asymmetric AlGaAs/GaAs quantum well heterostructure. L, length of the larger quantum well.

Fig. 4a. Energy level scheme of an asymmetrical AlGaAs/GaAs double quantum well heterostructure vs. the thickness of the larger quantum well.

E. Momox et al. / Journal of Computational Physics 231 (2012) 6173–6180

6179

Fig. 5. Conduction band diagram and squared wavefunctions for the QCL under an electric field of 48 kV/cm. The principal states are in bold. The wavefunctions were shifted with respect to their energy level.

4. Conclusion More often than not, the complexity of several numerical methods hampers its implementation by novices in numerical computing. As demonstrated in this paper, the described algorithm can be readily put into practice, as it ultimately copes with the presence of interfaces by solving a set of linear system of equations. It is used for the appraisal of the eigenstates in arbitrary potential profiles originated in semiconductor heterostructures by solving the Schrödinger equation. This novel numerical approach makes use of the Immersed Interface Method to locally modify the second order finite difference scheme to treat the presence of such interfaces or discontinuities. As for the computing time, it is a function of the number of grid points (N) employed and of the requested eigenstates, but in most cases about 104 grid points will suffice to obtain an entirely acceptable solution in a few seconds. Clearly, the use of sparse matrices accounts for the short computing time.

Acknowledgments The authors thank Mr. A. Boland-Thoms for proofreading an early version of the manuscript. E. Momox acknowledges the CONACyT for financial support.

References [1] C.D. Papageorgiou, A.D. Raptis, A method for the solution of the Schrodinger equation, Computer Physics Communications 43 (3) (1987) 325–328. [2] M.D. Feit, J.A. Fleck Jr., A. Steiger, Solution of the Schrodinger equation by a spectral method, Journal of Computational Physics 47 (3) (1982) 412–433. [3] D. Yevick, B. Hermansson, The accuracy of band-structure calculations based on the split-step fast Fourier transform method, Journal of Physics C: Solid State Physics 18 (22) (1985) 4303–4314. [4] Inho Kim, T.K. Gustafson, Lars Thylen, Analysis of quantum-confined structures using the beam propagation method, Applied Physics Letters 57 (3) (1990) 285–287. [5] Ali El-Hajj, Hafez Kobeissi, Mahmoud Korek, SSM: a set of subprograms for calculating eigenvalues for a diatomic molecule using a simplified shooting method, Computer Physics Communications 74 (2) (1993) 297–302. [6] H. Kobeissi, A. El-Hajj, M. Kobeissi, On the numerical integration of the Schrodinger equation for symmetrical potentials: a highly stable shooting method with the Numerov integrator, Journal of Physics A: Mathematical and General 23 (24) (1990) 5725–5732. [7] J. Killingbeck, Shooting methods for the Schrodinger equation, Journal of Physics A: Mathematical and General 20 (6) (1987) 1411–1418. [8] D. Indjin, Z. Ikonic, V. Milanovic, On shooting method variations for the 1-D Schrodinger equation and their accuracy, Computer Physics Communications 72 (2-3) (1992) 149–153. [9] B. Jonsson, S.T. Eng, Solving the Schrodinger equation in arbitrary quantum-well potential profiles using the transfer matrix method, Journal of Quantum Electronics 26 (11) (1990) 2025–2035. [10] Wayne W. Lui, Masao Fukuma, Exact solution of the Schrodinger equation across an arbitrary one-dimensional piecewise-linear potential barrier, Journal of Applied Physics 60 (5) (1986) 1555–1559. [11] Theodoros P. Horikis, Eigenstate calculation of arbitrary quantum structures, Physics Letters A 359 (5) (2006) 345–348. [12] K. Hayata, M. Koshiba, K. Nakamura, A. Shimizu, Eigenstate calculation of quantum well structures using finite elements, Electronics Letters 24 (10) (1988) 614–616. [13] K. Nakamura, A. Shimizu, M. Koshiba, K. Hayata, Finite-element analysis of quantum wells of arbitrary semiconductors with arbitrary potential profiles, Journal of Quantum Electronics 25 (5) (1989) 889–895. [14] Khai Q. Le, Finite element analysis of quantum states in layered quantum semiconductor structures with band nonparabolicity effect, Microwave and Optical Technology Letters 51 (1) (2009) 1–5. [15] I-H. Tan, G.L. Snider, L.D. Chang, E.L. Hu, A self-consistent solution of Schödinger-Poisson equations using a nonuniform mesh, Journal of Applied Physics 68 (8) (1990) 4071–4076.

6180

E. Momox et al. / Journal of Computational Physics 231 (2012) 6173–6180

[16] D.J. Costinett, T.P. Horikis, High-order eigenstates calculation of arbitrary quantum structures, Journal of Physics A: Mathematical and Theoretical 42 (23) (2009) 235201. [17] Randall J. Leveque, Zhilin Li, The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM Journal of Numerical Analysis 31 (4) (1994) 1019–1044. [18] R.B. Lehoucq, D.C. Sorensen, C. Yang, ARPACK User’s Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods, SIAM, Philadelphia, 1998. [19] G. Bastard, Superlattice band structure in the envelope-function approximation, Physical Review B 24 (10) (1981) 5693–5697. [20] R. Ferreira, G. Bastard, Evaluation of some scattering times for electrons in unbiased and biased single- and multiple-quantum-well structures, Physical Review B 40 (2) (1989) 1074–1086. [21] H. Page, C. Becker, A. Robertson, G. Glastre, V. Ortiz, C. Sirtori, 300 K operation of a GaAs-based quantum-cascade laser at k  9 lm, Applied Physics Letters 78 (22) (2001) 3529–3531.