Numerical solutions for a two dimensional quantum dot model

TL;DR

This paper revisits a 2D quantum dot model, comparing polynomial solutions with new numerical calculations, and extends the spectrum of eigenfunctions and eigenvalues obtainable for the system.

Contribution

It introduces a more general numerical approach using the Numerov method, enabling calculation of a broader spectrum of eigenfunctions and eigenvalues for the quantum dot model.

Findings

Good qualitative agreement between polynomial and numerical solutions

Extended the range of external harmonic frequencies for eigenfunction calculation

Predicted the existence of bound states for the planar system with l=0

Abstract

In this paper, a quantum dot mathematical model based on a two-dimensional Schr\"odinger equation assuming the 1/r inter-electronic potential is revisited. Generally, it is argued that the solutions of this model obtained by solving a biconfluent Heun equation have some limitations. The known polynomial solutions are confronted with new numerical calculations based on the Numerov method. A good qualitative agreement between them emerges. The numerical method being more general gives rise to new solutions. In particular, we are now able to calculate the quantum dot eigenfunctions for a much larger spectrum of external harmonic frequencies as compared to previous results. Also the existence of bound state for such planar system, in the case l=0, is predicted and its respective eigenvalue is determined.

Equations

$$\left\{-\frac{1}{2}({\nabla }_1^2+{\nabla }_2^2)+\frac{{\Omega }^2}{2}(r_1^2+r_2^2)+\frac{1}{|{\vec{r}}_1-{\vec{r}}_2|}\right\}\Psi ({\vec{r}}_1,{\vec{r}}_2)=E_T\Psi ({\vec{r}}_1,{\vec{r}}_2)$$

$$\left[-{\nabla }_{\vec{r}}^2+{\omega }^2{r}^2+\frac{1}{r}\right]\psi (\vec{r})=\eta \psi (\vec{r})$$

$$\psi (\vec{r})=r^{-1/2}u(r)e^{\pm i\ell \theta }$$

$$\frac{\text{mbox}{d}^{2}u(r)}{\text{mbox}d{r}^2}+\left[\eta -\frac{1}{r}-{\omega }^2{r}^2-\frac{({\ell }^2-1/4)}{r^2}\right]u(r)=0$$

$$x{y}^{\prime \prime }(x)+[1+\alpha -2x^2]y^{\prime }+[-\delta /2+(\gamma -\alpha -2)x]y(x)=0$$

$$\alpha =2\ell ;\gamma =\frac{\eta }{\omega };\delta =\frac{2}{\sqrt{\omega }}$$

$$u(r)=r^{\ell +1/2}{e}^{-\omega r^2/2}y(\sqrt{\omega }r)$$

$$u(r\pm \delta )=u(r)\pm \delta u^{\prime }(r)+\frac{{\delta }^2}{2}{u}^{\prime \prime }(r)\pm \frac{{\delta }^3}{6}{u}^{\prime \prime \prime }(r)+\frac{{\delta }^4}{24}{u}^{iv}(r)$$

$$u(r+\delta )+u(r-\delta )=2u(r)+{\delta }^2{u}^{\prime \prime }(r)+\frac{{\delta }^4}{12}{u}^{iv}(r)$$

$$\frac{u(r+\delta )+u(r-\delta )-2u(r)}{{\delta }^2}=\left(1+\frac{{\delta }^2}{12}\frac{\text{mbox}{d}^2}{\text{mbox}d{r}^2}\right)u^{\prime \prime }(r)$$

$$\frac{\text{mbox}{d}^{2}u}{\text{mbox}d{r}^2}+k^2(r)u(r)=0$$

$$\left(1+\frac{{\delta }^2}{12}\frac{\text{mbox}{d}^2}{\text{mbox}d{r}^2}\right)u^{\prime \prime }(r)=-k^2(r)u(r)-\frac{{\delta }^2}{12}\frac{\text{mbox}{d}^2}{\text{mbox}d{r}^2}\left[k^2(r)u(r)\right]$$

$$\left[1+\frac{{\delta }^2}{12}{k}^2(r+\delta )\right]u(r+\delta )=2\left[1-\frac{5{\delta }^2}{12}{k}^2(r)\right]u(r)-\left[1+\frac{{\delta }^2}{12}{k}^2(r-\delta )\right]u(r-\delta )$$

$$\gamma -\alpha -2=2n\Rightarrow {\eta }_{n\ell }=2\left(n+\ell +1\right)\omega$$

Full text

1 Introduction

Modern technics in nanometer-scale semiconductor manufacturing enable the creation of quantum confinement of only a few electrons. These few-body systems are often called quantum dots Reimann . They can be described by a model where the electrons move in an external harmonic oscillator potential of frequency $\Omega$, exhibiting a two dimensional behavior Sikorski . Previous numerical calculations suggest that the harmonic oscillator potential can be successfully employed to describe two-electron quantum dot Merkt . Thus, a model for this kind of system can be described by the time independent Schrödinger equation, in atomic units ($\hbar=m=e=1$), with energies and frequencies in hartree (Ha) units,

[TABLE]

The subscripts 1 and 2 refer to each one of the electrons. The $\vec{r}_{i}$ are two-dimensional vectors with length $r_{i}=|\vec{r}_{i}|$. Introducing the usual relative and center of mass coordinates, $\vec{r}=\vec{r}_{1}-\vec{r}_{2}$ and $\vec{R}=(\vec{r}_{1}+\vec{r}_{2})/2$, Eq. (1), with the choice $\Psi(\vec{r}_{1},\vec{r}_{2})=\chi(\vec{R})\psi(\vec{r})$, gives rise to the following relative coordinate equation:

[TABLE]

defining the frequency $\omega=\Omega/2$. The total energy is given by $E_{T}={\epsilon}+\eta$, with $\epsilon$ being the center of mass amount of energy. The 2$D$ radial Schrödinger equation for the relative coordinate can be obtained by introducing the polar coordinates $(r,\theta)$ and putting its solution in the form

[TABLE]

being $\ell$ the integer angular momentum quantum number of the two-dimensional system. The radial function $u(r)$ should satisfy the following equation:

[TABLE]

The problem of two electrons in an external oscillator potential is studied in three dimensions Taut , and it is shown that the above radial equation is quasi-exactly solvable Turbiner_1 ; Turbiner_2 ; Usheridze , which means that it is possible to find exact simple solutions for some, but not all, eigenfunctions, corresponding to a certain infinite set of discrete oscillator frequencies. In a recent paper Caruso ; Silveira , it was shown that it is possible to determine exactly and in a closed form a finite portion of the energy spectrum and the associated eigenfunctions for the Schrödinger equation describing the relative motion of a two-electron system, by putting Eq. (4) into the form of a biconfluent Heun equation, like

[TABLE]

where

[TABLE]

and the relation between the functions $u(r)$ and $y(x)$, with $x=\sqrt{\omega}\,r$, is given by

[TABLE]

This method, indeed, did give rise to polynomial solutions only for certain frequencies, $\Omega$. Each solution is obtained for a specific frequency value. The same is true for other studies Taut_94 ; Taut_2000 ; Taut_2010 . This is a significative limitation. Therefore, it is natural to wonder if a numerical analysis of this problem could give rise to a broad class of solutions well defined for any chosen external frequency value $\Omega$ as it is expected in a real experimental situation. This is the scope of this paper.

The first part of this paper is aimed to confirm that the Numerov method, applied to Eq. (4), is able to predict the previous analytical results. Doing this, we are ready to find new solutions for arbitrary $\Omega$ values, that were not possible with the polynomial method.

2 Numerical method

In this section, a short summary for those who are not familiarized with the method is given.

This numerical method was initially developed to determine solutions of eigenvalue problems associated with second order ordinary differential equations that did not contain the first order derivative term. But using Eq. (3), we can eliminate the term of the first order derivative.

In this method, the solution is considered to be known at two subsequent points in the interval $[a,b]$, for example in $u(r-\delta)$ and $u(r)$, $\delta$ being a parameter arbitrarily small. The next step is to determine the solution at the next point $u(r+\delta)$. Then, we expand the term $u(r\pm\delta)$ in Taylor series to the fourth order obtaining

[TABLE]

and then one can add the terms $u(r+\delta)$ and $u(r-\delta)$, resulting in

[TABLE]

In this case, only even order derivatives remain. Then, Eq. (8) can be written as

[TABLE]

Eq. (4) can be written in a convenient way as

[TABLE]

where $k^{2}=\eta-1/r-\omega^{2}r^{2}-(\ell^{2}-1/4)/r^{2}$. Applying the operator $\displaystyle\left(1+\frac{\delta^{2}}{12}\frac{\mbox{d}^{2}}{\mbox{d}r^{2}}\right)$ on both sides of Eq. (10) we obtain

[TABLE]

Finally, Comparing Eq. (9) to Eq. (11), we were able to write Numerov’s difference formula:

[TABLE]

To solve this formula we need to start from an initial hint for the eigenvalue that determine the $a$ and $b$ points in which the wave function is practically zero ($u(a)=u(b)=0$). This can be done graphically by analyzing the effective potential. Then it is needed to match the left and right solutions at one of the classical turning points.

The error involved here is smaller, $\mathcal{O}(\delta^{6})$ than that in other methods based on lower-order expansion, such as Runge-Kutta, $\mathcal{O}(\delta^{4})$.

For a better understanding of the method we suggest reading the references Numerov ; Numerov2 ; Numerov3 ; Numerov4 ; Numerov5 ; Numerov6 , The algorithm used in this paper was implemented in a program developed by the authors using $C^{++}$ language. Both calculations and graphics shown here were done by using the CERN/ROOT package.

3 First numerical results

As was said in the introduction, the solutions obtained by solving the Biconfluent Heun equation, Eq. (5), present a strong limitation. Indeed, each polynomial solution is valid just for one specific frequency. Therefore, it is difficult to compare the theoretical prediction with a particular experiment for a given external harmonic oscillator potential of frequency $\Omega=2\omega$. This limitation motivates us to search for numerical solutions of the radial Schrödinger equation, Eq. (4) as done in reference Caruso-Helayel .

First of all, we aimed to reproduce the linear relation, given by Eq. (12), which was deducted in Caruso ; Silveira , between the energy $\eta_{n\ell}$ and the quantum numbers $n$ and $l$. Such linear relation, for the particular case when $\ell=0$, is shown in Fig. 2. The frequency for this numerical result was fixed as $\omega=0.01$ Ha. Remember that since the quantum dot is submitted to a microwave external excitation, we can use values for the frequency $\Omega$ in the range 0.01 Ha $<\Omega<$ 1 Ha.

Another test that confirms that it is possible to obtain numerically the available analytical predictions was done by running the numerical program for the frequency value associated with each polynomial solution found in Caruso ; Silveira . For simplicity, we give here just the comparison considering the states $l=0$. The theoretical value of the energy, computed by Eq. (12), and the respective output from the numerical method are shown in Table 1, for $1\leq n\leq 5$.

[TABLE]

A simple inspection of Table 1 show a good qualitative agreement between previous results and our numerical values.

In Fig. 2 some theoretical polynomial eigenfunctions $u(r)$ calculated with the results of the references Caruso ; Silveira are compared to the respective numerical eigenfunctions.

For the $n=3$ case, the numerical eigenfunctions is slightly different from the analytic one, nevertheless, the energy numerically found for this state is almost identical to the analytical one as we can see in the Table 1.

In general, we can conclude that we also got a good agreement for the wave-functions.

4 Numerical eigenfunctions for any external frequency

Now that we know that the numerical program is able to retrieve the analytic solutions, we can look for new solutions that are not possible to obtain analytically with the chosen polynomial method.

For example, considering an external frequency $\omega=0.01$ Ha, we can actually find numerically several solutions. These results are given in Table 2, where Eq. (12) is used to calculate the analytical values.

Again, Table 2 shows how the results are close.

The numerical eigenfunctions for these states are shown in Fig. 3.

To the best of our knowledge, it is the first time that a set (of quantum numbers $n$) of quantum dot wave-functions are obtained for the same value of the external frequency $\Omega$ applied to the quantum dot.

Note that in this section we chose to work with the frequency $\omega=0.01$ Ha. However, we could have chosen a variety of other frequencies and obtained similar results.

5 Bound state solution

First of all, we have to understand that, in a strictly planar system, bound states can exist only for $\ell=0$. This peculiar fact depends on the nature of the effective potential of Eq. (4) since, only in two spatial dimensions, the so called “centrifugal potential” becomes indeed attractive just for the value $\ell=0$. Otherwise, the sum of Coulombic and centrifugal potentials are always repulsive. Notice that this kind of solution cannot be predicted by the polynomial methods.

We were able to numerically find just one state with energy value of $\eta=-63.92$ Ha, also with $\omega=0.01$ Ha. This will not be different for other choices of $\omega$. The reason for this is that, for small values of $r$ ($0<r<0.5$), the well shape is not at all bias by the $\omega$ choice (at least in the range we are considering in this paper). If this value is compared to the ground state energy of the hydrogen atom in three dimensions, $\eta_{H}=-0.5$ Ha, we see that both differ from two orders of magnitude.

The eigenfunction is given in Fig. 4.

This kind of discontinuity in the first derivative of the wave function, seen in Fig. 4, is typical of a $\delta(r)$-potential. It should be stressed that in this case and even in that of a very deep well with a small characteristic width (which is indeed our case) only one energy eigenstate is expected in the framework of non-relativistic quantum mechanics, as we have found.

6 Conclusion

The Numerov numerical method, applied in this paper to a mathematical quantum dot model, was able to reproduce the energies and eigenfunctions previously found analytically, with good accuracy. In addition, we were able to find new solutions that were not predicted by the polynomial method applied to solve the Biconfluent Heun equation describing this model. We also managed to find a bound state solution that can not be obtained analytically, since equation Eq. (12) only predicts positive energies.

The numerical calculation carried out here has also the advantage to be able to compute the quantum dot wave-function for more values of the external frequency that keeps the quantum dot confined.

As a last remark, it is important to stress that all the calculations made in this work start from a mathematical model for a quantum dot where the inter-electronic potential is supposed to be given by the Coulombic potential $1/r$. Thus, someone (a theoretician) must say that this is a toy model, since it is well known that, mathematically, the electric charge conservation, in a strictly two dimensional space, requires a logarithmic type potential. However, one cannot neglect the fact that there is a huge number of papers in the literature that still use a $1/r$ potential to describe the quantum dot model. The justification of this choice may be found in what an experimentalist would say, namely, that there is not a real “2D” system. Indeed any $2D$ system is an idealization that actually should be immersed in a true $3D$ space, where the potential between two electric charges is proportional to $1/r$. Which choice should we made? This question and their implications were discussed in another paper by the authors Silveira2 . Our point of view is that anyway one should expect to have experimental data in order to compare them to the predictions of both models.