Eigenstate modelling in arbitrary shaped nanostructres with gradual heterointerfaces

We consider the problem of evaluating the energy levels and envelope wavefunctions of electrons in semiconductor nanostructures with a smooth confining potential, such as inter- diffused quantum dots or quantum wells. We show how Gaussian functions can be used to model various confinement potentials in 1D, 2D and 3D. We develop a fast and accurate numerical algorithm for solving the one-particle Schrodinger equation with a position-dependent effective mass for this type of potential. We examine the convergence and accuracy of our algorithm. This method can be used to analyze experimental data, such as the PL spectra of quantum wells, wires and dots and the charge carrier transport through various nanostructures. GaAs/Al0.3Ga0.7As system is used in calculations.


Introduction
The problem of accurate evaluation of bound states of electrons and holes in semiconductor nanostructures is not trivial. For practical purposes such as modelling the photoluminescence spectra, tunnelling probability or charging and ionization of quantum dots the computational algorithm should be accurate and easy to use without requiring excessive computing power.
The Schrodinger equation can be solved analytically only for a handful of potentials. In 2D and 3D they are either symmetric or separable in certain coordinate systems [1]. For semiconductor heterostructures with varying composition the effective mass is positiondependent further complicating the problem. Numerical algorithms for one-dimensional problems are fairly straightforward, such as the transfer matrix method [2]. In two and three dimensions however, it is not as easy to correctly define the boundary conditions and solve the eigenvalue problem for an arbitrary non-separable confining potential. Various methods continue to be developed [3,4] but most of them are limited either by required computing power or by the amount of information which can be obtained.
In this work we show that Gaussian functions can be used to model arbitrary shaped nanostructures in 1D, 2D or 3D and develop a fast and accurate numerical algorithm for solving the Schrodinger equation for this type of potential. We use the effective mass, one particle, and envelope function approximation. Nonparabolicity and particle interaction are not considered here, but the method can be modified to include these effects. We use GaAs/Al 0.3 Ga 0. 7 As system in all calculations in this paper. Both the conduction band offset and effective mass depend linearly on the AlAs mole fraction t for t < 0.45.

Gaussian potential
We propose the following general potential for modelling the confinement in real nanostructures: This potential is very convenient. It is smooth, local (vanishes at infinity) and each term can be represented as a product U (x)U (y)U (z), which makes the resulting Schrodinger equation easy to solve. This potential is also very flexible, since it does not rely on any type of symmetry and thus can be used to model real, asymmetric nanostructures. Due to its locality, it can also be used to model nanostructure complexes (superlattices, nanowire and quantum dot arrays). Figure 1 shows how the rectangular quantum well can be modelled very accurately by a sum of 9 Gaussian functions. The energy levels and wavefunctions in the latter case are calculated by our method described below. Their differences in this case are ∆E 1 ≈ 0.3 meV and ∆E 2 ≈ 1.5 meV .  For semiconductor heterostructures we also need to consider the position dependence of the effective mass. If its dependence on AlAs mole fraction is linear, the effective mass for a single Gaussian well is defined as follows (m * b -the effective mass of barrier material at infinity, m * wthe effective mass in the center of the well):

The method in general
We consider the Schrodinger equation with the position-dependent effective mass [5]: We then use the spectral method first proposed in [6] for nanowires and further developed in [7][8][9][10] for quantum dots with abrupt interfaces. We expand the wavefunction into a complete orthonormal basis of eigenfunctions of a simple potential, such as a quantum box or a harmonic oscillator, both of which are separable in Cartesian coordinates in 3D.
Then we can substitute (4) into (3), multiply (3) by Ψ * pqs on the left and integrate over the whole space, obtaining a matrix eigenvalue problem: If we evaluate matrix elements on the left side we can solve this problem numerically, obtaining the complete set of eigenvectors {C jkl } n and eigenvalues E n . As we can only diagonalize a finite matrix, we have to truncate the expansion (4) and the result now depends on the number of basis functions N x , N y , N z and their exact form (for example quantum box dimensions L x , L y , L z or harmonic oscillator frequencies ω x , ω y , ω z ). To obtain the best possible agreement with the exact result for certain N x , N y , N z , we choose the parameters variationally (so as to minimize the ground state energy).
Unlike other numerical methods, such as a widely used finite element method, this algorithm does not require discretization in space or energy. It outputs normed continuous wavefunctions and all energy levels, including degenerate ones. This method is general and can be used with any potential with addition of electric and magnetic fields. The goal of this work is to create a fast, simple and accurate algorithm for the chosen family of potentials. By evaluating all matrix elements analytically, we will ensure that there is no numerical integration, which can be costly and introduce additional error.

Matrix elements calculation
We will first obtain the explicit expressions for matrix elements in the case of a single Gaussian well in one dimension. It will be easy to extend them to 2 and 3 dimensions. In this section we also consider the effective mass to be constant in space. We will use two basis sets. First is the set of solutions of harmonic oscillator (HO) potential. They have the following form (H j -Hermite polynomials).
Here j = 0, 1, 2, ... Then matrix elements will have the form: Here we used the relation for H j (gx) from [11]. The second basis is the set of solutions of the quantum box (QB) problem: Here j = 1, 2, 3, ... To obtain matrix elements for the potential analytically we need to integrate: Now let us assume that L a. Then we can move the integration boundaries to infinity and (11) is reduced to the Euler-Poisson integral. While this assumption may lead to additional error, when we compare the result to the HO basis calculation (figure 3) they show very good agreement.
The algorithm was implemented in Wolfram Mathematica 10. In our experience, the most time-consuming operation is calculation of the matrix elements. Numerical matrix diagonalization is performed much faster. Both operations scale with number of basis functions as N 2d with d being the number of dimensions. In one dimension the calculations are fast and take several seconds with a personal computer even for N = 100. Thus, we were able to analyze the continuous dependence of the results on the parameter L for HO and QB bases (figure 3). We can see definite intervals of convergence where the energy does not depend on L. These intervals become longer as N increases. For two dimensions the region of convergence was found. For simple potentials it is almost square, meaning we can choose L x , L y and L z independently in their own intervals of convergence.
The accuracy of the algorithm can also be checked by comparing the results with analytically solvable potentials (figure 1). A good example is modified Poschl-Teller potential which has and is very close in shape to the Gaussian potential. Indeed, for a = 5 nm their ground state energies are almost equal (∆E 1 ≈ 1 meV ). For wider wells ∆E 1 decreases. We also used exactly solvable symmetrical potentials in 2D and 3D to verify the results of our method in several dimensions.

Position dependent effective mass
Now we allow the effective mass to be position dependent. For the QB basis the matrix elements can still be calculated semi-analytically (using the Taylor expansion). The sum in (16) converges very fast, giving at least 5 correct digits for 20 terms. We compare calculation results for the constant and position dependent effective mass (figure 4). We can see that in the latter case the energy levels decrease (∆E 1 ≈ −2 meV , ∆E 2 ≈ −3 meV , ∆E 3 ≈ −5 meV ) and electron becomes more localized. It can be also understood qualitatively, since m * b > m * w , thus the mean effective mass increases when we account for m * b . In conclusion, we have shown that Gaussian functions can be used to model arbitrary shaped confinement potentials in any number of dimensions. We have developed a fast and accurate algorithm for evaluating energy levels and wave functions for this family of potentials and examined its convergence. The electric and magnetic fields can be included in straightforward way. Note that this method can be modified to include the effects of nonparabolicity and particle interaction (for example, by an iterative algorithm).