Semiconductor structures are the main source of actual technologies in the areas of sensing, computing, electronics, optical electronics and communications. Bandgap engineering is the main tool to obtain the necessary properties of semiconductor devices. Particularly, layered semiconductor structures containing quasi-two-dimensional electron gas become one of the most popular configurations for existing and perspective applications1. Therefore, the interest of scientific community to this kind of materials grows exponentially2.
Atomic impurities in semiconductor devices generally considered as an interference in an ideal crystallographic structure of semiconductor and are therefore one of the main reasons for the resistance of the whole structure. That is why modulation doping is mainly used outside the active areas as a source of free carriers. However, there are technologies that use impurities in semiconductor structures as a key feature3. In our previous works we proposed some ideas about possible applications of impurity-connected effects in quantum well (QW) type structures4,5,6,7.
The first work to calculate impurity states energies in a quantum well8used variational technique and gave birth to a series of subsequent publications9,10,11among others). This kind of method remains actual up to very recently12,13, etc.). Nevertheless, the information variational method can provide us is somewhat limited. It does not allow to reliably calculate the impurity wave function (WF), which implies that we cannot, for example, obtain matrix elements and, accordingly, the rates for all kinds of electronic transitions involving impurity states. Alternatively, expansion techniques, in which the impurity states wave function is expanded over QW subband wave functions, are used as well14,15,16. Those methods provide both wave functions and energy positions. However, they require significant processing powers and have their limitations. For instance, the well should be sufficiently deep to contain enough subbands to converge the results.
It is known that shallow impurity atom in/near the QW produces several states resembling Rydberg series. Some impurity states go out of the bandgap and appear in the background of the continuous energy spectrum of lower subbands. Such states are called resonant ones (with a reference to Fano resonance17) and they appear below every energy subband situated above the gap edge. Sometimes they described as broadened states.
In this work we describe in detail the numerical method based on the expansion mathematical model proposed by14. It can be used to numerically obtain energy positions and WFs of resonant and non-resonant impurity state in QW-like layered semiconductor heteroctructures. We use it to calculate some of the impurity states associated to a hydrogenic donor doped to the center of a rectangular GaAs/AlGaAs QW. Our results agree with an older work using the same mathematical model15. However, it contains much more data, which were missed by that work presumably due to the lack of access to necessary computational powers. Such data allowed us to propose the new classification of impurity states based on energy positions and the shapes of wave functions in QW-like semiconductor structures.
Mathematical model
To apply the method below to the 2D semiconductor nanostructure first we need to find the subband electronic description of the structure without the impurity atom. To find the subband energies and wave functions ξj(z) we must solve the one-dimensional time-independent Schrodinger equation18
$$\:{H}_{0}\xi\:\left(z\right)=E\xi\:\left(z\right)$$
(1)
For the simple case of one particle problem in a rectangular quantum well the Hamiltonian can be written as
$$\:{H}_{0}=-\frac{{\hbar}^{2}}{2{m}^{*}}\frac{{d}^{2}}{d{z}^{2}}+V\left(z\right).$$
(2)
Here the potential part \(\:V\left(z\right)\) is a well profile \(\:V\left(z\right)=\left\{\begin{array}{c}{V}_{b},\:\:\:\:\:\:\:\:\:\left|z\right|>L/2\\\:0,\:\:\:\:\:\:\:\:\:\:\:\left|z\right|\le\:L/2\end{array}\right\}\), L is a well width with Vb being a barrier height. In the more complex cases of external fields the Hamiltonian (2) is modified by adding respective terms. Such Hamiltonian is used when there is a small concentration of impurities or at low, usually helium, temperatures. In the case of elevated temperatures and high concentrations one can add a Hartree potential term, created by electrons and ionized impurities to the Hamiltonian and solve the Eq. (1) self-consistently along with Poisson and -maybe- electroneutrality ones. Anyway, that is a trivial task, which can be performed by available commercial19and free20 software available. For the solution we need to specify the structure geometry and some properties of the materials, particularly effective mass \(\:{m}^{*}\)(z). and dielectric permittivity ε(z), as function of structure growth direction coordinate z. As a result of such a solution we get al.l the necessary data for our problem, namely the wave functions of localized subbands \(\:{\xi\:}_{i}\left(z\right)\), and energy positions of subband bottoms \(\:{E}_{i}\).
As a base of the method to find impurity levels we used the one developed in14. The Hamiltonian to find the impurity states is written in cylindric coordinates r, Ѳ, z to take advantage of azimuthal symmetry of the system (assuming that in the directions perpendicular to z all the parameters are constant):
$$\:H=-\frac{{\hbar}^{2}}{2{m}^{*}}\frac{{d}^{2}}{d{z}^{2}}-\frac{{\hbar}^{2}}{2{m}^{*}}\left(\frac{{d}^{2}}{d{R}^{2}}+\frac{1}{R}\frac{d}{dR}+\frac{1}{{R}^{2}}\frac{{d}^{2}}{d{\theta\:}^{2}}\right)+V\left(z\right)-\frac{{e}^{2}}{4\pi\:\epsilon\:{\epsilon\:}_{0}\sqrt{{R}^{2}+{\left(z-{z}_{D}\right)}^{2}}}$$
(3)
where \(\:{m}^{*}\) is the electron effective mass, zD is z position of impurity center.
Accordingly, the WF \(\Psi\left(R,\theta,z\right)\) is sought as
$$\:{\Psi\:}\left(R,\theta\:,z\right)\equiv\:{e}^{im\theta\:}{\psi\:}_{m}\left(R,z\right),$$
(4)
were m is the azimuthal quantum number and i is the imaginary unit. \(\:{\psi\:}_{m}\left(R,z\right)\) is presented as an expansion over subband WFs:
$$\:{\psi\:}_{m}\left(R,z\right)=\sum\:_{j}{f}_{j}\left(R\right){\xi\:}_{j}\left(z\right).$$
(5)
Now using the solution of the Schrodinger equation with Hamiltonian \(\:{H}_{0}\) along with the equation with \(\:H\) after some math one can get the Eq.
$$\:\left(\frac{{d}^{2}}{d{R}^{2}}+\frac{1}{R}\frac{d}{dR}+\left({k}_{N}^{2}-\frac{{m}^{2}}{{R}^{2}}\right)\right){f}_{N}\left(R\right)={U}_{Nn}\left(R\right){f}_{n}\left(R\right),$$
(6)
were \(\:{k}_{N}^{2}=\frac{2{m}^{*}}{{\hbar}^{2}}(E-{E}_{N})\), \(\:{U}_{Nn}\left(R\right)=\frac{2{m}^{*}}{{\hbar}^{2}}\frac{1}{4\pi\:\epsilon\:{\epsilon\:}_{0}}\int\:{\xi\:}_{N}^{\text{*}}\left(z\right)\frac{{e}^{2}}{\sqrt{{R}^{2}+{\left(z-{z}_{D}\right)}^{2}}}{\xi\:}_{n}\left(z\right)dz\), and \(\:{E}_{N}\) (N = 1,2,3…) are energy solutions of (1).
The solution of (6) is found as14:
$$\:{f}_{N}\left(R\right)=\underset{0}{\overset{{\infty\:}}{\int\:}}{G}_{N}\left(R,{R}^{{\prime\:}},E\right)\sum\:_{n}{U}_{Nn}\left({R}^{{\prime\:}}\right){f}_{n}\left({R}^{{\prime\:}}\right)R{\prime\:}dR{\prime\:},$$
(7)
where \(\:{G}_{N}\left(R,{R}^{{\prime\:}},E\right)\) are Green functions given by:
$$\:{G}_{N}\left(R,{R}^{{\prime\:}},E\right)=\left\{\begin{array}{c}-{K}_{m}\left({k}_{N}{R}^{{\prime\:}}\right){I}_{m}\left({k}_{N}R\right),\:\:R<{R}^{{\prime\:}},\:E<{E}_{N}\\\:-{I}_{m}\left({k}_{N}{R}^{{\prime\:}}\right){K}_{m}\left({k}_{N}R\right),\:\:\:R\ge\:{R}^{{\prime\:}},\:E<{E}_{N}\\\:\frac{\pi\:}{2}{N}_{m}\left({k}_{N}{R}^{{\prime\:}}\right){J}_{m}\left({k}_{N}R\right),\:\:\:\:R<{R}^{{\prime\:}},E\ge\:{E}_{N}\:\\\:\frac{\pi\:}{2}{J}_{m}\left({k}_{N}{R}^{{\prime\:}}\right){N}_{m}\left({k}_{N}R\right),\:\:\:\:R\ge\:{R}^{{\prime\:}},E\ge\:{E}_{N}\end{array}\right.,$$
(8)
where J, N (I, K) are the Bessel functions of first (second) kind.
To solve the Eqs. (6–8) an overcomplicated version of shooting method was used in14. It was based on border conditions and continuity of WF expansion terms. Instead we employ a simpler way described below. Apart from this, we use the same system of equations and expect to get the same solutions within the numerical error.
Replacing in (7) integral with a sum of finite differences we get:
$$\:{f}_{N}\left({R}_{p}\right)=\sum\:_{q}\left({G}_{N}\left({R}_{p},{R}_{q},E\right){R}_{q}{\varDelta\:R}_{q}\sum\:_{n}{U}_{Nn}\left({R}_{q}\right){f}_{n}\left({R}_{q}\right)\right),$$
(9)
Where \(\:{R}_{p}\)and \(\:{R}_{q}\) get the same series of NR consequent magnitudes between 0 and some \(\:{R}_{max}\), \(\:{\varDelta\:R}_{j}={R}_{j+1}-{R}_{j}\:(j=p,q\)) \(\:p,q=1..{N}_{R}\). We did not find any numerical advantage of choosing variable \(\:\varDelta\:R\:\)and used the simple case of uniform step \(\:\varDelta\:R=const=\frac{{R}_{max}}{{N}_{R}-1}\), where \(\:\:{R}_{j}=(j-1)\times\:\varDelta\:R\). For given E, \(\:{R}_{p}\) and \(\:{R}_{q}\), \(\:{G}_{N}\left({R}_{p},{R}_{q},E\right)\) and \(\:{U}_{Nn}\left({R}_{q}\right)\) can be calculated, and (9) in fact is a system of linear algebraic equations with respect to variables \(\:f\left(R\right)\). The system consists of Ns×NR equations (Ns is the number of subbands used for expansion) and can be presented as:
$$\:M\varvec{f}=0,$$
(10)
Where f is a vector of Ns×NR elements and contains all tabulated functions fN, starting from 1st to \(\:{N}_{s}\)-th (\(\:{f}_{p+\left(N-1\right)*{N}_{R}}={f}_{N}\left({R}_{p}\right),\:N=1..{N}_{s},\:p=1..{N}_{R}\)), whilst M is a square matrix:
$$\:{M}_{p+\left(N-1\right)\text{*}{N}_{R},\:\:\:q+\left(n-1\right)\text{*}{N}_{R}}={\delta\:}_{pq}{\delta\:}_{nN}-{G}_{N}\left({R}_{p},{R}_{q},E\right){R}_{q}{\varDelta\:R}_{q}{U}_{Nn}\left({R}_{q}\right)$$
(11)
where \(\:N,n=1..{N}_{s},\:p,q=1..{N}_{R}\).
In21 it is shown that energy positions of both resonant and non-resonant impurity levels correspond to energy values for which the condition Det(M) = 0 fulfills. So, the solution of this equation gives impurity level energies and the solution of the system (9) provides with the set of Ns×NR numerical elements of column vector \(\:{f}_{j}\left(R\right)\), which allows to build the complete impurity WFs using Eqs. (4) and (5). Mathematically the equation system (9) has infinite number of solutions but owing to it is linear, we can just equate any of \(\:{f}_{j}\) scalar variables to any real constant, find other variables and then normalize the WF (4) which can be written in cylindric coordinates as:
$$\:{{\Psi\:}}_{norm}\left(R,\theta\:,z\right)=\frac{{\Psi\:}\left(R,\theta\:,z\right)}{\sqrt{2\pi\:\underset{z=-{\infty\:}}{\overset{{\infty\:}}{\int\:}}\underset{R=0}{\overset{{\infty\:}}{\int\:}}{{\Psi\:}(\text{R},\text{z})}^{2}RdRdz}}$$
(12)
We guess that this -or similar- simple method was not proposed in14 because at the time the author did not have access to computers with big enough memory to properly deal with very large matrices. For instance, a typical matrix appearing in this work contained about 576,000,000 real elements.
To characterize the impurity WFs we used transverse \(\:{a}_{R}\) and longitudinal \(\:{a}_{z}\) localization characteristics calculated as follows.
$$\:{a}_{R}=\underset{V}{\overset{}{\int\:}}{{\Psi\:}}^{2}RdV=2\pi\:\underset{z=-{\infty\:}}{\overset{{\infty\:}}{\int\:}}dz\underset{R=0}{\overset{{\infty\:}}{\int\:}}{{\Psi\:}(\text{R},\text{z})}^{2}{R}^{2}dR$$
(13)
$$\:{a}_{z}=\underset{V}{\overset{}{\int\:}}{{\Psi\:}}^{2}RdV=2\pi\:\underset{z=-{\infty\:}}{\overset{{\infty\:}}{\int\:}}\left|{z}_{cm}-z\right|dz\underset{R=0}{\overset{{\infty\:}}{\int\:}}{{\Psi\:}(\text{R},\text{z})}^{2}RdR$$
(14)
One may interpret \(\:{a}_{R}\) and \(\:{a}_{z}\) as some kind of quantities equivalent to Bohr radius, defined as the most probable transverse and longitudinal distances to the center of mass of the WF \(\:{z}_{cm}\). In this work we considered only the case of a center doped well \(\:{z}_{D}=0\), therefore the center of the WF coincides with the well center \(\:{z}_{cm}=0\) because of the well symmetry along z.
Numerical proceeding
First, we numerically find the subband WFs \(\:{\xi\:}_{j}\left(z\right)\) and respective energies \(\:{E}_{j}\) as eigenfunctions and eigenvalues of Hamiltonian \(\:{H}_{0}\) (Eq. (2)), employing effective mass \(\:{m}^{*}\)=0.067m0 (m0 is a free electron mass), dielectric constant \(\:\epsilon\:\)=12.565 and rectangular well profile: well width L=30 nm, and barrier potential Vb=50Ry (with effective Rydberg value 1Ry=\(\:\frac{{\hbar}^{2}}{2{m}^{*}{a}_{b}^{2}}=\)5.80meV, \(\:{a}_{b}=\frac{4\pi\:\epsilon\:{\epsilon\:}_{0}{\hbar}^{2}}{{m}^{*}{e}^{2}}\)). These parameters of our structure correspond to one structure used in15, in order to compare the outcoming results. The task above is trivial and we shall not discuss it.
In our reference system the center of the well is at z = 0 and we take z in units of L. Here we consider impurity position in the center of the well zD=0.
Now we form the matrix M as (10) and evaluate its discriminant D(E) for a given energy E. To find the energy positions of impurity we must solve the equation D(E) = 0 (that provides a series of solutions described below) and then to construct the impurity wave functions for the obtained energies, we solve the Eq. 9 obtaining the corresponding functions fj and form the full WFs according to Eqs. (4) and (5). Note that to correctly reproduce the WF the energy should be found with high numerical precision (we used 7 significative digits) for exactly the same set of numerical parameters.
The primary possible numerical error sources for the proceeding above are the following: (1) \(\:\varDelta\:R\) in finite differences, (2) number of subbands used for expansion Ns, (3) the maximum magnitude of R used to replace the infinity in the integral of Eq. 7 for the numerical purposes \(\:{{R}_{max}=N}_{R}\varDelta\:R\). Our calculations showed that for all impurity states we report in this work varying parameters \(\:\varDelta\:R\) and Ns provide good convergence of the results. Thus, for our calculations we stuck with Ns=7 and \(\:\varDelta\:R\)=0.07 nm, values such that negligible changes appear when refined further. It is worth noting here that among the solutions found, the most vulnerable to an increase in ∆R turned out to be the lowest energy one, corresponding to the wave function most localized in R, which is quite expectable.
A more complex situation occurs when varying Rmax. Obviously, Rmax should cover absolutely most of the WF, so the bigger Rmax the better. On the other hand, we find that when Rmax is too big in comparison to the WF extension, the calculated WF can exhibit an artificial peak at the border R = Rmax. The impurity level energy is not affected much by Rmax, which causes that peak, if we continue increasing Rmax, it remains stable but at some Rmax the energy solution can disappear. We associate this behavior with an equivalent of numerical integration error. On the other hand, decreasing \(\:\varDelta\:R\) does not have much effect on energy solution. So, for better reproduction of impurity WFs and to be able to find the energy of some specific impurity levels, the calculations should be performed with Rmax within the specific reasonable range that depends on the localization of respective state.
The Fig. 1 is an illustration of the series of solutions provided by the method. Here the modified determinant \(\:{D}_{M}\) of the matrix M is depicted as a function of energy for different Rmax. \(\:{D}_{M}\) is used for visualization purposes in the graph and calculated as follows
Modified determinant DM of the matrix M for GaAs/AlGaAs 30 nm-wide rectangular quantum well with barrier height Vb=50Ry as a function of energy in Ry for different values of Rmax. s-levels (m = 0), Rmax=120 (orange squares), 180 (grey diamonds), 240 (light blue triangles), 480 (blue minuses) and 720 nm (green circles).
$$\:{D}_{M}=\frac{\left|D\right|}{D}\times\:\left(\sqrt{{\text{l}\text{n}}^{2}\left({D}^{2}\right)+\beta\:}+\text{l}\text{n}\left({D}^{2}\right)\right),$$
(15)
where D is a matrix determinant in some arbitrary units, and β is a parameter used to adjust the visualization (here β = 100). This visualization concept has some advantages of semilogarithmic scale (particularly it conserves the smoothness and monotonousness of original determinant while allowing comparison of curve shapes with different orders of magnitude) and preserve the sign of original determinant which is important to see its zeros.
In Fig. 1 the solutions appear at the points where the curve intersects the horizontal axis. We can see that some of the solutions coincide for different Rmax or at least converge to some energy when Rmax grows. We shall refer to those solutions as to numerically stable ones and assume that a stable solution corresponds to a real physical impurity level. Other solutions are unstable; that is they change significantly with Rmax and do not converge within reasonable magnitudes of Rmax. Such numerically unstable solutions may not have any physical meaning and we will ignore them. As one can see, the stable solutions tend to form a ladder asymptotically approaching the bottom of each subband whereas unstable ones are either above the subband bottoms and below the stable solutions of the next subband; or just in the immediate proximity of subband bottom. The latter can be attributed to the upper steps of respective ladders and can be calculated more exactly with bigger Rmax. All the plots in Fig. 1 are for m = 0. We performed similar calculations for m = 1 and 2 and identified stable levels.
In the Fig. 2 we present DM graphs with curves for m = 0,1 and 2 for Rmax=240 nm (which we used for calculations). The graphs for m different than 0 have the same general behavior (stable levels just below respective subbands with somewhat bigger energies for zero determinant, unstable levels above the subbands) with one notable exception: an unstable level with m = 2 just below the first stable m = 0 level (marked in the figure).
Modified determinant DM GaAs/AlGaAs 30 nm-wide rectangular quantum well with barrier height Vb=50Ry as a function of energy in Ry for different m: m=0 (s, orange squares), m=1 (p, blue triangles), m=2 (d, grey circles), Rmax=240 nm.
Analyzing Figs. 1 and 2 one can conclude that the set of solutions provided by zero determinants resembles a kind of Rydberg series with three independent quantum numbers. In particular, we explicitly have the azimuthal one, m, which appears in exponent of the wave function (Eq. (4)), and then determines the order of Bessel functions in Eq. (8). Two other quantum numbers can be deduced from existing solutions for each m. One of them corresponds to the number of specific subband number and another one to the step number within the respective ladder. Naturally we expect that increasing each of those quantum numbers leads to increasing energy and less localized corresponding wave function.
Here and below, for a proper classification of impurity levels we will employ the following notation. Following15, we use letters “s”, “p”, “d”, “f” … to refer to m = 0, 1, 2, 3… accordingly; in loose analogy to the classification of electronic states in hydrogenic atom. However, it is necessary to have in mind that the analogy is far from perfect and can be confusing in some situations. In regard, a particular impurity state is denoted as “letter, n1, n2”, where n1 is a subband number where 1 corresponds to the ground subband, and n2 is a ladder number starting from 1 for the lowest level in the ladder. For example, s11 is a ground impurity state that corresponds to −0.606Ry and s21 is a first step of the ladder below the second subband (2.622Ry), see Fig. 1.
