This article describes the use of the SACI package—a package for calculating the energy levels and wavefunctions of a multi-electron quantum dot modelled as a 2D harmonic well with electrons interacting through a Coulomb potential and under the influence of a perpendicular magnetic field.
|This article has not been updated for Mathematica 8.
Quantum dots are artificially fabricated atoms, in which charge carriers are confined in all three dimensions just like electrons in real atoms. Consequently, they exhibit properties normally associated with real atoms such as quantised energy levels and shell structures. These properties are described by the electron wavefunctions whose evolution is governed by the Schrödinger equation and the Pauli exclusion principle.
There are many methods available to solve the Schrödinger equation for multiple electrons. They roughly fall into the categories of the diagonalisation method, mean-field density-functional theory, and the self-consistent field approach. One of the first theoretical studies of quantum dots was by Pfannkuche et al. , who compared the results of Hartree-Fock self-consistent calculations and exact diagonalisation of the Hamiltonian for two electrons in a circularly symmetric parabolic potential. They found good agreement between the two methods for the triplet state but marked differences for the singlet state, indicating important spin correlations were not included properly in their Hartree-Fock model. This suggests that the proper treatment of electron spins is crucial for correctly obtaining the electronic structures in quantum dots.
Examples of self-consistent field approaches in the literature include Yannouleas and Landman [2, 3], who studied circularly symmetric quantum dots using an unrestricted spin-space Hartree-Fock approach, and McCarthy et al. , who developed a Hartree-Fock Mathematica package. Macucci et al.  studied quantum dots with up to 24 electrons using a mean-field local-density-functional approach, in which the spin exchange-correlation potential was approximated by an empirical polynomial expression. Lee et al.  also studied an -electron quantum dot using density-functional theory, where the generalised-gradient approximation was used for exchange-correlation potentials. Exchange interaction comes directly from the antisymmetrisation of wavefunctions as required by Pauli’s exclusion principle. In density-functional theory this is a major problem since the mathematical object is the electron density rather than the electron wavefunction, making evaluation of the exchange interaction intrinsically difficult.
Diagonalisation approaches in the literature include Ezaki et al. [7, 8], Eto [9, 10], Reimann et al. , and Reimann and Manninen , who each applied a brute force approach by numerically diagonalising the -electron Hamiltonian using Slater determinants composed of single-electron eigenstates as the basis functions. This approach, namely the configuration interaction (CI) method, takes into account the full interaction and correlation of the electrons in the system as long as the numerical results converge with an increasing number of basis functions. However, such an approach involves the calculation of a very large number of interaction integrals and the inversion of large matrices, which can be prohibitively expensive in terms of computer resources. Reimann et al. [11 employed matrices of dimensions up to 108,375 with 67,521,121 nonzero elements for a six-electron quantum dot. Calculations for any higher number of electrons were not considered numerically viable using the conventional CI formalism, even with state-of-the-art computing facilities.
We have recently developed a spin-adapted configuration interaction (SACI) method to study the electronic structure of multi-electron quantum dots . This method is based on earlier work by quantum chemist R. Pauncz , which expands the multi-electron wavefunctions as linear combinations of antisymmetrised products of spatial wavefunctions and spin eigenfunctions. The SACI method has an advantage over using Slater determinants in that a smaller basis is used. This reduces the computational resources required, allowing calculations for dots with more than six electrons on a desktop computer.
After some theoretical background, we present results from a Mathematica package which employs the SACI method to calculate the energy levels and wavefunctions of multiple electrons confined by a 2D harmonic well and under the influence of a perpendicular magnetic field. Mathematica enables the exact calculation of interaction integrals which greatly improves the speed and accuracy of calculations.
The SACI package encapsulates all the functionality needed to calculate the energy levels and wavefunctions for electrons in a 2D harmonic well potential with/without the influence of a perpendicular magnetic field.
Ensure that SACI.m is on your $Path and load the package.
You can also choose Input Get File Path to locate SACI.m on your computer. This package will be used in the Theory section to demonstrate results and produce examples, as well as in the Calculations section to calculate energy levels and wavefunctions.
When measuring the component of an electron’s spin angular momentum along an axis, the result is either or . Conventionally, we denote this axis as the -axis and work in units in which is one. The two spin eigenstates in which we are certain of the electron’s spin projection are denoted and , which are the eigenfunctions of the projection operator. To find projected spin eigenfunctions for multiple electrons we can simply multiply single electron spin eigenstates (e.g., or for five electrons). We call these simple multiplications of ’s and ’s elementary spin eigenfunctions. Quantum mechanics also allows linear combinations of projected spin eigenfunctions (e.g., ). As long as each term in the spin function has the same number of ’s and ’s, we will still have an eigenfunction of , that is , where is the spin projection quantum number given by half the number of spin up electrons minus half the number of spin down electrons ( in this case).
The square of the total spin is also quantised:
where is a multi-electron spin eigenfunction (spin eigenfunction) and is the total spin quantum number. Elementary spin eigenfunctions with the same number of ’s and ’s will be eigenfunctions of , each with the same projected spin quantum number , but they are not in general eigenfunctions of . Thus will in general be a linear combination of elementary spin eigenfunctions.
One way to calculate the spin eigenfunction is to use the Dirac identity to write the operator in terms of permutations
where is the number of electrons confined in the quantum dot and is the permutation operator. To form eigenfunctions of , we can calculate the action of on each of the elementary spin eigenfunctions. We convert into a matrix with the vector of coefficients of the elementary spin eigenfunctions in the linear combination. We can then solve the matrix eigenvalue problem to calculate all the spin eigenfunctions. There may be more than one eigenfunction for a given total spin quantum number, and so we carry out an orthonormalisation procedure to calculate an orthonormal spin eigenspace.
In the SACI package, we define and implicitly using ElectronList, in which 1 represents a spin up electron and 0 represents a spin down electron.
is the length of the list and is defined by the number of ones and zeros in the list. In this case and . We then define the total spin quantum number .
must be compatible with the number of electrons and to get meaningful results. With , , and in place we can calculate the spin basis and then display it.
We can see that we have a 2D spin eigenspace as there are two basis vectors produced. Calculating the spin eigenfunctions is a prerequisite to the rest of the calculations. Henceforth, we label spin eigenvectors , where is the number of electrons, and are the spin quantum numbers, and is an index to identify vectors in a multi-dimensional basis.
we get the energy eigenvalues
and the wavefunctions
where is the effective mass, the harmonic well constant, the cyclotron frequency, the effective harmonic well constant, the magnetic length, , the magnetic field, and the generalised Laguerre polynomials.
To coincide with experiments done on a circularly symmetric quantum, we will choose solid-state parameters appropriate for GaAs (i.e., and , where is the electron mass and is the vacuum permittivity).
When we are working in effective atomic units, the unit of length is the effective Bohr radius. This is given by the function lengthscale with the result in nanometers.
Here is a plot of the probability density for in GaAs in a harmonic well in which the gap between the first and second energy levels is 10 meV and there is no magnetic field. Fixing the first energy gap is a convenient way to define the harmonic well constant as it can be determined experimentally. We also set the magnetic field to 1 Tesla and plot six effective Bohr radii from the centre of the well.
The base scale is in nanometers and indicates the scale of the quantum dot.
Multiplying one-electron functions together (e.g., ) allows us to form multi-electron spatial wavefunctions that can be combined with the spin eigenfunctions. However, the combinations of spin and spatial wavefunctions must obey Pauli’s exclusion principle as described in the next section.
The Antisymmetry and Representation Matrices
The Pauli exclusion principle states that because electrons are indistinguishable fermions (half-integer spin particles), any wavefunction that describes the motion of electrons must change sign if both the spin and spatial coordinates of any pair of electrons are interchanged. Consequently in the -electron case, if we apply any permutation to the electrons, the following holds:
where is the multi-electron wavefunction and is the parity of the permutation. To antisymmetrise a wavefunction we apply the antisymmetrisation operator
where is the symmetric group (the set of all permutations of electrons under composition). A wavefunction must be an eigenfunction of the antisymmetrisation operator if it is to fully describe the physical properties of the electrons.
Since the permutation and spin operators commute, if we apply a permutation to a spin eigenfunction, the result is another spin eigenfunction. This new spin eigenfunction is not, in general, one of the functions already calculated but a linear combination of them. In fact, the coefficients of the linear combination form matrices which together form a representation of the symmetric group. Explicitly
where is the spin eigenfunction for electrons with total spin quantum numbers and in an eigenspace of dimension , and is the element of the representation matrix . It can be shown that for any two permutations and
where is the composite permutation derived by acting with first then . Equation (9) satisfies the condition for a representation of the symmetric group. These representation matrices and their properties are useful in simplifying the calculation of the Hamiltonian matrix.
We can see that our spin eigenfunctions are eigenfunctions of by the representation matrix of The permutation is given by a reordering of the list. For example, in a three-electron system, represents the transposition of electrons 1 and 2 and the transposition of electrons 2 and 3.
This shows spin eigenfunction two is unchanged by the permutation and spin eigenfunction one changes sign under the permutation. In general the action of a permutation can lead to a linear combination of the spin eigenfunctions.
This shows that under the permutation the spin eigenfunctions transform via a linear combination. The following demonstrates the relation given by equation (9).
To get an -electron spin-spatial basis function we can multiply single-electron wavefunctions, given by equation (5), combine them with a spin eigenfunction of and , and then antisymmetrise the result. The antisymmetrisation process restricts which spatial and spin combinations give nonzero results as discussed in detail in reference .
We use a five-electron quantum dot as an example, where the system wavefunction and there are two doubly occupied spatial orbitals. Two possible spin eigenfunctions for are
Both and are antisymmetric under permutation , but is antisymmetric under permutation , while is symmetric under permutation . This means vanishes but is nonzero.
The CalculateSpinBasis function applies a procedure to ensure the spin basis functions are also eigenfunctions of the permutations , … . This property is required for the Hamiltonian matrix element formulae, described in the Matrix Elements Rules section, to be valid. A basis whose spin eigenfunctions are simultaneously eigenfunctions of , … are guaranteed to exist as the permutation operators together with and pairwise commute.
The multi-electron Hamiltonian is given as the following:
This is simply the sum of the single-electron energies defined by equation (3) plus the pairwise Coulomb interactions. The fact that we are working in a solid state media is modelled by using the relative effective mass and the effective permittivity .
We wish to solve the following equation:
for the energy levels and wavefunctions . The standard diagonalisation method is to approximate as a sum of a finite sum of basis functions with the convergence improving as we increase the size of the basis. This converts the problem into a matrix eigenvalue problem. However, to get a good approximation to the solution we need to choose the basis wisely as there is a limit to the size of the matrix that is computationally feasible. We also need to be able to efficiently calculate the Hamiltonian matrix elements.
The SACI method uses normalised, antisymmetrised products of spin eigenfunctions and spatial wavefunctions given by products of Fock-Darwin solutions as basis elements. The elements of our Hamiltonian matrix are then for various spatial wavefunctions and , spin eigenfunctions and , and normalisation constants and .
The SACI method reduces the number of basis elements required by ensuring the wavefunction satisfies the Pauli exclusion principle from the outset and restricting the calculation to specific spin quantum numbers and .
Although each of the basis elements is a sum of terms, the properties of the antisymmetriser and the spin eigenfunctions allow enormous simplification of the calculation similar to the Slater-Condon rules for matrix elements between Slater determinates. In effect, only the product spatial wavefunction and spin eigenfunction need to be specified for each basis element with the normalisation constant and antisymmetrisation built into the rules for calculating the Hamiltonian elements.
We need to be careful in the selection of basis elements, however, as we need an orthonormal basis for the diagonalisation expansion to be valid. Remember, some spin and spatial combinations vanish and thus cannot be included in the basis. We must also make sure the set of basis elements is mutually orthogonal and normalised (i.e., if and , and 0 otherwise). To do this, we select only spatial wavefunctions that have doubly occupied orbitals in the first and second position, the third and fourth position, … to match the spin eigenfunctions symmetry or antisymmetry under the operators , , … . The properties of the antisymmetriser then guarantee the production of a complete orthonormal basis.
The SACI package automates the process of the orthonormal basis generation. As we can select only a finite number of basis elements out of an infinite set, we need to have some method of selecting appropriate basis elements. As we are usually interested in the ground and first few excited states of a quantum dot, we have chosen to use basis elements in increasing order of the sum of their component one-electron energies.
Elaborating, with no magnetic field the energy levels of the Fock-Darwin wavefunctions are given by . If there were no interaction between the electrons, the energy of a multi-electron state would simply be the sum of the one-electron energies, for example, the energy of would be . We select a basis consisting of all those elements with a noninteracting energy less than a certain cutoff.
The variable cutoffenergy is set to a cutoff energy of . We can then create a basis.
BasisCreate produces a list of basis elements with noninteracting energy less the specified cutoff. The spatial component of each basis element is given by a list of pairs of integers denoting the quantum numbers and for each orbital. The number at the end of the list denotes the index of the spin eigenfunction (its position when ShowSpinStatus runs) with which it is combined.
We notice that spin eigenfunction “1″ is antisymmetric with respect to the permutation , so it can be combined with the doubly occupied orbitals, but spin eigenfunction “2″ is symmetric with respect to , so it can only be combined with nondoubly occupied orbitals.
Matrix Elements Rules
As mentioned previously, there are general rules for calculating Hamiltonian matrix elements. These rules are coded in the function InteractionHamiltonianMatrixElement, and their derivation can be found in . In brief, when the spatial wavefunctions of the two basis elements differ by more than two orbitals (i.e., after being sorted so the orbitals are in maximum correspondence), then the matrix element is zero; when the spatial orbitals differ by two or less orbitals, a different rule applies for each of the three cases (i.e., zero, one, and two orbitals differ). These formulas are quite general and can be used for any interaction Hamiltonian that acts pairwise.
In the Hamiltonian element formulas we must calculate the two-electron interaction integral. For the case of electrons in a 2D harmonic well interaction through a Coulomb potential under the influence of a perpendicular magnetic field, the integral is in the following form:
We need to be able to do this integral for many combinations of the quantum numbers , and . We notice that when and , we have a singularity as the Coulomb interaction goes to infinity as the electrons move closer together.
However by using Mathematica and a little ingenuity, we can do this integral exactly in all cases needed in computation. The detailed derivation (breaking down and transforming the integral into a sum of components) can be found in . For example, we set the quantum numbers , and in order.
We also note that if (i.e., the sum or the orbital angular momentum quantum numbers are not equal), the interaction integral is zero.
These integrals are done with the effective harmonic well constant set to one. Values of the interaction integral for the effective harmonic well constant can easily be transformed to any value of . Thus when we do an interaction integral, it is done for and transformed to the correct value of afterwards. This allows us to store interaction integral values for various combinations of , and as function definitions and save them to use later. This process is done dynamically so that values are saved as definitions only if they are used.
Here are some examples of interaction Hamiltonian matrix elements indicating that they can be done exactly.
It can also be proved that the interaction element is zero if the sum of the orbital angular momentum quantum numbers are not the same, hence, the last calculation leads to zero.
The preceding calculations show that the spin eigenfunction with which a spatial wavefunction is combined makes a nontrivial contribution to the Hamiltonian matrix element calculation.
In this section we will demonstrate how the package can be used to do various calculations such as energy levels, wavefunctions, convergence plots, and energy level evolutions with magnetic fields.
To calculate energy levels, we need to choose a basis. For larger basis sizes the calculation will take longer and require more memory (the size of the matrix is the size of the basis squared) but the convergence will be better.
We choose a basis with a single total orbital angular momentum quantum number . The properties of the circularly symmetric quantum dot allow us to separate calculations for different values of , thus reducing the size of the basis needed for each calculation allowing for greater accuracy. We then calculate the interaction Hamiltonian using the function InteractionHamiltonian.
The calculation of this matrix can be time-consuming, so it is worthwhile to save the result, as the interaction Hamiltonian calculated for a particular number of electrons and spin quantum numbers and can then be used for any value of the harmonic well constant, magnetic field, relative effective mass, and relative permittivity. We have also loaded previously calculated integral definitions using the function LoadIntegralandHamiltonianData to speed up our calculation and saved the data afterward using SaveIntegralandHamiltonianData to speed up future calculations run in other sessions.
Next we set the magnetic field strength to 0 Tesla remembering that firstenergygap was set to 10 meV previously and relativeeffectivemass and relativepermittivity are set to the parameters for GaAs.
We then can calculate the Hamiltonian for GaAs with these particular values of the harmonic well constant and magnetic field and plot the energy levels.
The plot shows the first 10 energy levels for a GaAs quantum dot system of three electrons with spin quantum numbers , and total angular momentum quantum number in a zero magnetic field. The harmonic well constant is implicitly defined by a first energy gap of 10 meV.
As mentioned, the accuracy of our result as well as the computational time required increases with the size of the basis. We can test the convergence with the following series of functions.
As the interaction Hamiltonian for lower cutoff energies is a submatrix of interaction Hamiltonians with higher cutoff energies, we can do convergence tests using a single interaction Hamiltonian. We will use the interaction Hamiltonian from the previous section.
The Convergence function takes submatrices of the Hamiltonian using different basis sizes as specified by the input sizelist and calculates the energy of the levels specified by wantedlevels. The ground state is specified by 1, the first excited state by 2, and so on. This allows the convergence properties of the solution to be ascertained.
The list basis size can be specified manually, but a convenient method is to choose all the basis sizes that correspond to the different cutoff energies.
In the preceding code, we have specified wantedlevels to examine the ground state and first excited state. We set cutoffmin to be the first cutoff energy that will produce the highest energy level requested, as we can only approximate up to levels when we have basis elements.
As there may be large gaps in the basis sizes, we can also add in points at intermediate intervals. Here, we add a basis size of 85.
We can then calculate the energies of the ground state and first excited state as we increase the basis size. The energies are in meV.
The following plot shows convergence of energies of the two lowest states.
We can see that the convergence is quite good for relatively small basis sizes in this case.
Electron Density Plots
We can extract electron density plots for various energy levels by calculating the normalised energy eigenvectors of the full Hamiltonian and using the function ElectronDensity. We can then plot the electron density and use lengthscale to give the length dimensions in nanometers.
First we calculate the energy eigenvectors for the first six energy levels.
Next we use ElectronDensityFunction to convert the normalised energy eigenvector into a probability density function by combining it with the basis basisLzero. This is then simplified to cancel out the complex components. As the density function is a measurable quantity, we are guaranteed that the complex components will cancel. Chop is used to remove any complex components that remain due to rounding.
We can see that we end up with a numerical approximation to the wavefunction. Here is an example.
The following plot shows the electron probability density functions for the first six energy levels with , , and for the three-electron quantum dot system that has been the subject of our examples.
Energy Levels versus Magnetic Fields
The SACI package also allows the calculation of the evolution of the ground state with a magnetic field using a single interaction Hamiltonian.
We can use LevelVsField to calculate a list of energy levels for a range of magnetic fields defined by magmin, magmax, and step given in Tesla, and for the levels defined by wantedlevels. Here we choose the first six energy levels and plot for the magnetic field range of 0 to 4 Tesla at intervals of every 0.5 Tesla.
Calculating the energy levels for one value at a time has the benefit of producing more accurate results for a given matrix size and allows for easier separation of different energy states from the table of eigenvalues. The following plot shows the results.
Note that although we are plotting six levels, two pairs of these remain degenerate as the magnetic field changes, so we only see four lines in the plot. Running the program several times, it is possible to plot curves with different spin quantum numbers and total orbital angular momentum quantum numbers on the same graph. From this, it is possible to see the graphs for different quantum numbers intersecting as the magnetic field increases. This corresponds to magnetic transitions in the energy levels. Magnetic transitions in the ground state of quantum dots have been measured experimentally.
This article presents a Mathematica implementation of the SACI method for calculating the energies and wavefunctions of multi-electron quantum dot systems. The results obtained from this package are highly accurate with established confidence for up to eight electrons when using a PC with a Pentium IV 2.4GHz processor, which can be readily extended by using more powerful computers. Such numerically exact calculations provide important benchmarks against which one can test other approximate schemes developed to study more complex systems.
|||D. Pfannkuche, V. V. Gudmundsson, and P. A. Maksym, “Comparison of a Hartree, a Hartree-Fock, and an Exact Treatment of Quantum-Dot Helium,” Physical Review B (Condensed Matter), 47(4), 1993 pp. 2244-2250.|
|||C. Yannouleas and U. Landman, “Spontaneous Symmetry Breaking in Single and Molecular Quantum Dots,” Physical Review Letters, 82(26), 1999 pp. 5325-5328.|
|||C. Yannouleas and U. Landman, “Collective and Independent-Particle Motion in Two-Electron Artificial Atoms,” Physical Review Letters, 85(8), 2000 pp. 1726-1729.|
|||S. A. McCarthy, J. B. Wang, and P. C. Abbott, “Electronic Structure of an -Electron Quantum Dot,” Computer Physics Communications, 141, 2001 pp. 175-204.|
|||M. Macucci, K. Hess, and G. J. Iafrate, “Simulation of Electronic Properties and Capacitance of Quantum Dots,” Journal of Applied Physics, 77(7), 1995 pp. 3267-3276.|
|||I.-H. Lee, V. Rao, R. M. Martin, and J.-P, Leburton, “Shell Filling of Artificial Atoms within Density-functional Theory,” Physical Review B (Condensed Matter), 57(15), 1998 pp. 9035-9041.|
|||T. Ezaki, N. Mori, and C. Hamaguchi, “Electronic Structures in Circular, Elliptic, and Triangular Quantum Dots,” Physical Review B (Condensed Matter), 56(11), 1997 pp. 6428-6431.|
|||T. Ezaki, Y. Sugimoto, N. Mori, and C. Hamaguchi, “Electronic Properties in Quantum Dots with Asymmetric Confining Potential,” Semiconductor Science and Technology, 13(8a), 1998 pp. A1-A3.|
|||M. Eto, “Electronic Structures of Few Electrons in a Quantum Dot under Magnetic Fields,” Japanese Journal of Applied Physics, 36, 1997 pp. 3924-3927.|
|||M. Eto, “Many-Body States in a Quantum Dot under High Magnetic Fields,” Japanese Journal of Applied Physics, 38, 1999 pp. 376-379.|
|||S. M. Reimann, M. Koskinen, and M. Manninen, “Formation of Wigner Molecules in Small Quantum Dots,” Physical Review B (Condensed Matter), 62(12), 2000 pp. 8108-8113.|
|||S. M. Reimann and M. Manninen, “Electronic Structure of Quantum Dots,” Reviews of Modern Physics, 74(4), 2002 pp. 1283-1343.|
|||J. B. Wang, C. Hines, and R. D. Muhandiramge, “Electronic Structure of Quantum Dots,” The Handbook of Theoretical and Computational Nanoscience (M. Rieth and W. Schommers, eds.), Stevenson Ranch, CA: American Scientific Publishers, 2006.|
|||R. Pauncz, The Construction of Spin Eigenfunctions: An Exercise Book, New York: Kluwer Academic/Plenum Publishers, 2000.|
|||V. A. Fock, “Bemerkung zur Quantelung des Harmonischen Oszillators im Magnetfeld,” Zeitschrift fur Physics, 47(5-8), 1928 pp. 446-448.|
|||C. G. Darwin, “The Diamagnetism of the Free Electron,” Mathematical Proceedings of the Cambridge Philosophical Society, 27, 1930 pp. 86-90.|
|R. D Muhandiramge, and J. Wang, “Electronic Structure of Multi-Electron Quantum Dots,” The Mathematica Journal, 2012. dx.doi.org/10.3888/tmj.10.2-7.|
Available at www.mathematica-journal.com/data/uploads/2012/05/SACI.m.
About the Authors
Ranga D. Muhandiramge is a mathematics Ph.D. student at the University of Western Australia. In 2003, Muhandiramge did his honours project jointly in mathematics and physics modelling quantum dots, for which he received first-class honours and a prize for the most outstanding honours graduand. He currently works in the area of operations research on mine field path planning with the assistance of scholarships from the Hackett Foundation and the Defense Science and Technology Organisation.
Jingbo Wang is an associate professor of physics at The University of Western Australia. Her research areas and interests range widely from atomic physics, molecular and chemical physics, spectroscopy, acoustics, chaos, nanostructured electronic devices, and mesoscopic physics to quantum information and computation. Wang has published over 100 research papers in peer-reviewed journals and international conference proceedings. Wang uses Mathematica extensively in both her teaching and research, and she is also an editorial board member of the Journal of Computational and Theoretical Nanoscience.
Ranga D. Muhandiramge
School of Mathematics and Statistics, M019
University of Western Australia
35 Stirling Highway
Crawley, WA 6009, Australia
School of Physics
University of Western Australia
35 Stirling Highway
Crawley, WA 6009, Australia