The aim of this project is to develop your own code for solving eigenvalue problems. The matrix to diagonalize is the same as the one we encountered in project one, a so-called tridiagonal Toeplitz matrix. This matrix has analytical eigenpairs (eigenvalues and eigenvectors) and gives us an excellent testing ground for our algorithms. In this project we will develop an eigenvalue solver based on Jacobi's method. The project will also introduce you to units tests and we will compare our results against other eigenvalue solvers (from LAPACK and/or numpy).
This project aims also at introducing to you the concept of scaling of equations. This means often either making various variables dimensionless or introducing units which are more convenient.
We will start with the two-point boundary value problem of a buckling beam or a spring fastened at both ends. This is one of the problems which has analytical solutions. Thereafter, by simply adding a new variable along the diagonal elements, we can study quantum mechanical problems. In particular, we will study the harmonic oscillator problem in three dimensions, with one or two electrons. For the latter case we can study the role of the repulsive Coulomb interaction and extract interesting physics results. For selected frequencies, this interacting two-electron problem exhibits analytical solutions, one of the few cases of an interacting system where wecan find analytical solutions. See M. Taut, Phys. Rev. A 48, 3561 (1993) for the derivation of analytical expressions for the eigenpairs..
Electrons confined in small areas in semiconductors, so-called quantum dots, form a hot research area in modern solid-state physics, with applications spanning from such diverse fields as quantum nano-medicine to the contruction of quantum gates. The article on quantum computing with quantum dots by Loss and DiVincenzo is an excellent read for those interested in this exciting topic.
We start with the following differential equation, namely $$ \gamma \frac{d^2 u(x)}{dx^2} = -F u(x), $$ where \( u(x) \) is the vertical displacement of the beam in the \( y \) direction. The beam has length \( L \), \( x\in [0,L] \) and \( F \) is a force applied at \( (L,0) \) in the direction towards the origin. The parameter \( \gamma \) is a constant defined by properties like the rigidity of the beam. We apply again so-called Dirichlet boundary conditions and set \( u(0)=u(L)=0 \).
In this specific case two of the parameters \( \gamma \), \( F \) and \( L \) are known. As an example, assume we know \( F \) and \( L \). Then the eigenvalue problem we set up below will allow us to find \( \gamma \).
We define a dimensional variable $$ \rho = \frac{x}{L}, $$ meaning that we have \( \rho \in [0,1] \). By reordering the equation as $$ \frac{d^2 u(\rho)}{d\rho^2} = -\frac{FL^2}{\gamma} u(\rho)=-\lambda u(\rho), $$ with \( \lambda= FL^2/\gamma \) we have an equation that when discretized, becomes an eigenvalue problem. We use the same expression for the second derivative of a function \( u \) as we did in project 1, namely $$ \begin{equation} u''=\frac{u(\rho+h) -2u(\rho) +u(\rho-h)}{h^2} +O(h^2), \tag{1} \end{equation} $$ where \( h \) is our step. Next we define minimum and maximum values for the variable \( \rho \), \( \rho_{\mathrm{min}}=0 \) and \( \rho_{\mathrm{max}}=1 \), respectively. With a given number of mesh points, \( N \), we define the step length \( h \) as, with \( \rho_{\mathrm{min}}=\rho_0 \) and \( \rho_{\mathrm{max}}=\rho_N \), $$ \begin{equation*} h=\frac{\rho_N-\rho_0 }{N}. \end{equation*} $$ The value of \( \rho \) at a point \( i \) is then $$ \rho_i= \rho_0 + ih \hspace{1cm} i=1,2,\dots , N. $$ We can rewrite the differential equation for a value \( \rho_i \) as $$ -\frac{u(\rho_i+h) -2u(\rho_i) +u(\rho_i-h)}{h^2} = \lambda u(\rho_i), $$ or in a more compact way as $$ -\frac{u_{i+1} -2u_i +u_{i-1} }{h^2} = \lambda u_i. $$ Following our approach from project 1, we can rewrite this equation in a more a general form, but now as an eigenvalue problem, that is $$ \begin{equation} \begin{bmatrix} d& a & 0 & 0 & \dots &0 & 0 \\ a & d & a & 0 & \dots &0 &0 \\ 0 & a & d & a &0 &\dots & 0\\ \dots & \dots & \dots & \dots &\dots &\dots & \dots\\ 0 & \dots & \dots & \dots &a &d & a\\ 0 & \dots & \dots & \dots &\dots &a & d\end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \\ \dots \\ u_{N-2} \\ u_{N-1}\end{bmatrix} = \lambda \begin{bmatrix} u_1 \\ u_2 \\ u_3 \\ \dots \\ u_{N-2} \\ u_{N-1}\end{bmatrix} . \tag{2} \end{equation} $$ As in project 1, we have not included the endpoints \( u_0 \) and \( u_N \). We have defined \( d=2/h^2 \) and the non-diagonal ones as \( a=-1/h^2 \). This eigenvalue problem has analytical eigenpairs, with eigenvalues given as (note that we define the matrix to run from \( i=1 \) to \( i=N-1 \)) $$ \lambda_j = d+2a\cos{(\frac{j\pi}{N})}, \hspace{0.2cm} j=1,2,\dots N-1. $$
The eigenvectors are $$ \mathbf{u}_j = [ \sin(\frac{j\pi}{N}), \sin(\frac{2j\pi}{N}), ..., \sin(\frac{(N-1)j\pi}{N}) ]^T, \hspace{0.2cm} j=1,2,\dots N-1. $$
The lecture notes of Tom Lyche provide an excellent description of this matrix. Chapters 1.2 and 1.3 of this text present a thourough discussion of tridiagonal matrices and two-point boundary value problems.
Your task now is to write a function which implements Jacobi's rotation algorithm (see Lecture notes chapter 7) in order to solve Eq. (2). However, the first simple check is to set up the matrix to diagonalize for a given \( N \) and use either armadillo's or numpy's functions for diagonalizing a matrix. You can then check that you obtain the analytical eigenvalues.
For Jacobi's method, we define the quantities \( \tan\theta = t= s/c \), with \( s=\sin\theta \) and \( c=\cos\theta \) and $$ \begin{equation*}\cot 2\theta=\tau = \frac{a_{ll}-a_{kk}}{2a_{kl}}. \end{equation*} $$ We can then define the angle \( \theta \) so that the non-diagonal matrix elements of the transformed matrix \( a_{kl} \) become non-zero and we obtain the quadratic equation (using \( \cot 2\theta=1/2(\cot \theta-\tan\theta) \) $$ \begin{equation*} t^2+2\tau t-1= 0, \end{equation*} $$ resulting in $$ \begin{equation*} t = -\tau \pm \sqrt{1+\tau^2}, \end{equation*} $$ and \( c \) and \( s \) are easily obtained via $$ \begin{equation*} c = \frac{1}{\sqrt{1+t^2}}, \end{equation*} $$ and \( s=tc \).
How many similarity transformations are needed before you reach a result where all non-diagonal matrix elements are essentially zero? Try to estimate the number of transformations and extract a behavior as function of the dimensionality of the matrix. Compare your results with the analytical ones.
Plot also the eigenvector for the lowest eigenvalue and compare this eigenvector with the analytical solution.
You can check your results against the Armadillo function for solving eigenvalue problems. The armadillo function eig$_$sym can be used to find eigenvalues and eigenvectors. A Python program which solves this part of the project is available under the project writing slides.
Comment your results (here you could for example compute the time needed for both algorithms for a given dimensionality of the matrix).
In this project (and later ones as well) we will implement so-called unit tests. Our unit tests are mainly meant to test mathematical properties of our algorithm. During the development phase of a program it is useful to devise tests that your program should pass. One of these is to make sure that for a given simple test matrix (say a \( 5\times 5 \) matrix) our algorithm for searching for the largest non-diagonal element always returns the correct answer. Furthermore, for a known simple matrix, irrespective of changes made, we should always get the same and correct eigenvalues. Another test is to make sure that the orthogonality shown in exercise (a) is preserved. Try to figure out other tests your code should pass, based either on the mathematical properties of the algorithms or more program specific tests. Implement at least two unit tests in this project.
Here we will assume that these electrons move in a three-dimensional harmonic oscillator potential (they are confined by for example quadrupole fields) and repel each other via the static Coulomb interaction. We assume spherical symmetry. You don't need to think of quantum mechanics when solving this project, look at it as another diagonalization problem.
We are first interested in the solution of the radial part of Schroedinger's equation for one electron. This equation reads $$ \begin{equation*} -\frac{\hbar^2}{2 m} \left ( \frac{1}{r^2} \frac{d}{dr} r^2 \frac{d}{dr} - \frac{l (l + 1)}{r^2} \right )R(r) + V(r) R(r) = E R(r). \end{equation*} $$ In our case \( V(r) \) is the harmonic oscillator potential \( (1/2)kr^2 \) with \( k=m\omega^2 \) and \( E \) is the energy of the harmonic oscillator in three dimensions. The oscillator frequency is \( \omega \) and the energies are $$ \begin{equation*} E_{nl}= \hbar \omega \left(2n+l+\frac{3}{2}\right), \end{equation*} $$ with \( n=0,1,2,\dots \) and \( l=0,1,2,\dots \).
Since we have made a transformation to spherical coordinates it means that \( r\in [0,\infty) \). The quantum number \( l \) is the orbital momentum of the electron. Then we substitute \( R(r) = (1/r) u(r) \) and obtain $$ \begin{equation*} -\frac{\hbar^2}{2 m} \frac{d^2}{dr^2} u(r) + \left ( V(r) + \frac{l (l + 1)}{r^2}\frac{\hbar^2}{2 m} \right ) u(r) = E u(r) . \end{equation*} $$ The boundary conditions are \( u(0)=0 \) and \( u(\infty)=0 \).
We introduce a dimensionless variable \( \rho = (1/\alpha) r \) where \( \alpha \) is a constant with dimension length and get $$ \begin{equation*} -\frac{\hbar^2}{2 m \alpha^2} \frac{d^2}{d\rho^2} u(\rho) + \left ( V(\rho) + \frac{l (l + 1)}{\rho^2} \frac{\hbar^2}{2 m\alpha^2} \right ) u(\rho) = E u(\rho) . \end{equation*} $$ We will set in this project \( l=0 \). Inserting \( V(\rho) = (1/2) k \alpha^2\rho^2 \) we end up with $$ \begin{equation*} -\frac{\hbar^2}{2 m \alpha^2} \frac{d^2}{d\rho^2} u(\rho) + \frac{k}{2} \alpha^2\rho^2u(\rho) = E u(\rho) . \end{equation*} $$ We multiply thereafter with \( 2m\alpha^2/\hbar^2 \) on both sides and obtain $$ \begin{equation*} -\frac{d^2}{d\rho^2} u(\rho) + \frac{mk}{\hbar^2} \alpha^4\rho^2u(\rho) = \frac{2m\alpha^2}{\hbar^2}E u(\rho) . \end{equation*} $$ The constant \( \alpha \) can now be fixed so that $$ \begin{equation*} \frac{mk}{\hbar^2} \alpha^4 = 1, \end{equation*} $$ or $$ \begin{equation*} \alpha = \left(\frac{\hbar^2}{mk}\right)^{1/4}. \end{equation*} $$ Defining $$ \begin{equation*} \lambda = \frac{2m\alpha^2}{\hbar^2}E, \end{equation*} $$ we can rewrite Schroedinger's equation as $$ \begin{equation*} -\frac{d^2}{d\rho^2} u(\rho) + \rho^2u(\rho) = \lambda u(\rho) . \end{equation*} $$ This is the first equation to solve numerically. In three dimensions the eigenvalues for \( l=0 \) are \( \lambda_0=3,\lambda_1=7,\lambda_2=11,\dots . \)
We define minimum and maximum values for the variable \( \rho \), \( \rho_{\mathrm{min}}=0 \) and \( \rho_{\mathrm{max}} \), respectively. You need to check your results for the energies against different values \( \rho_{\mathrm{max}} \), since we cannot set \( \rho_{\mathrm{max}}=\infty \).
With a given number of mesh points, \( N \), we define the step length \( h \) as, with \( \rho_{\mathrm{min}}=\rho_0 \) and \( \rho_{\mathrm{max}}=\rho_N \), $$ \begin{equation*} h=\frac{\rho_N-\rho_0 }{N}. \end{equation*} $$ The value of \( \rho \) at a point \( i \) is then $$ \rho_i= \rho_0 + ih \hspace{1cm} i=1,2,\dots , N. $$ We can rewrite the Schroedinger equation for a value \( \rho_i \) as $$ -\frac{u(\rho_i+h) -2u(\rho_i) +u(\rho_i-h)}{h^2}+\rho_i^2u(\rho_i) = \lambda u(\rho_i), $$ or in a more compact way $$ -\frac{u_{i+1} -2u_i +u_{i-1}}{h^2}+\rho_i^2u_i=-\frac{u_{i+1} -2u_i +u_{i-1} }{h^2}+V_iu_i = \lambda u_i, $$ where \( V_i=\rho_i^2 \) is the harmonic oscillator potential.
We define first the diagonal matrix element $$ \begin{equation*} d_i=\frac{2}{h^2}+V_i, \end{equation*} $$ and the non-diagonal matrix element $$ \begin{equation*} e_i=-\frac{1}{h^2}. \end{equation*} $$ In this case the non-diagonal matrix elements are given by a mere constant. All non-diagonal matrix elements are equal. With these definitions the Schroedinger equation takes the following form $$ \begin{equation*} d_iu_i+e_{i}u_{i-1}+e_{i}u_{i+1} = \lambda u_i, \end{equation*} $$ where \( u_i \) is unknown. We can write the latter equation as a matrix eigenvalue problem $$ \begin{equation} \begin{bmatrix}d_1 & e_1 & 0 & 0 & \dots &0 & 0 \\ e_1 & d_2 & e_2 & 0 & \dots &0 &0 \\ 0 & e_2 & d_3 & e_3 &0 &\dots & 0\\ \dots & \dots & \dots & \dots &\dots &\dots & \dots\\ 0 & \dots & \dots & \dots &\dots e_{N-3} &d_{N-2} & e_{N-2}\\ 0 & \dots & \dots & \dots &\dots &e_{N-2} & d_{N-1} \end{bmatrix} \begin{bmatrix} u_{1} \\ u_{2} \\ \dots\\ \dots\\ \dots\\ u_{N-1} \end{bmatrix}=\lambda \begin{bmatrix} u_{1} \\ u_{2} \\ \dots\\ \dots\\ \dots\\ u_{N-1} \end{bmatrix}. \tag{3} \end{equation} $$
Note here that we have not included the endpoints since we assume that they are knowm (as we did in project 1). The matrix is also symmetric in this case.
Add the harmonic oscillator potential to your tridiagonal Toeplitz matrix from 2a-2c and diagonalize the matrix. Study the results as functions of the number of integration points \( N \) and your approximation to \( \rho_{\mathrm{max}} \). The analytical results with our scaling for the one-electron energies are \( \lambda=3,7,11,15,\dots \). How many integration points do you need in order to reproduce the analytical results with say four leading digits after the decimal point?
You can reuse the code you wrote for part (b), but you need to add the potential \( \rho^2 \) to the diagonal elements.
Note: if you are not interested in the explicit quantum physics case here, you can replace this part with either 2g or 2h below. The latter ones are optional exercises if you do this.
We will now study two electrons in a harmonic oscillator well which also interact via a repulsive Coulomb interaction. Let us start with the single-electron equation written as $$ \begin{equation*} -\frac{\hbar^2}{2 m} \frac{d^2}{dr^2} u(r) + \frac{1}{2}k r^2u(r) = E^{(1)} u(r), \end{equation*} $$ where \( E^{(1)} \) stands for the energy with one electron only. For two electrons with no repulsive Coulomb interaction, we have the following Schroedinger equation $$ \begin{equation*} \left( -\frac{\hbar^2}{2 m} \frac{d^2}{dr_1^2} -\frac{\hbar^2}{2 m} \frac{d^2}{dr_2^2}+ \frac{1}{2}k r_1^2+ \frac{1}{2}k r_2^2\right)u(r_1,r_2) = E^{(2)} u(r_1,r_2) . \end{equation*} $$
Note that we deal with a two-electron wave function \( u(r_1,r_2) \) and two-electron energy \( E^{(2)} \).
With no interaction this can be written out as the product of two single-electron wave functions, that is we have a solution on closed form.
We introduce the relative coordinate \( \mathbf{r} = \mathbf{r}_1-\mathbf{r}_2 \) and the center-of-mass coordinate \( \mathbf{R} = 1/2(\mathbf{r}_1+\mathbf{r}_2) \). With these new coordinates, the radial Schroedinger equation reads $$ \begin{equation*} \left( -\frac{\hbar^2}{m} \frac{d^2}{dr^2} -\frac{\hbar^2}{4 m} \frac{d^2}{dR^2}+ \frac{1}{4} k r^2+ kR^2\right)u(r,R) = E^{(2)} u(r,R). \end{equation*} $$
The equations for \( r \) and \( R \) can be separated via the ansatz for the wave function \( u(r,R) = \psi(r)\phi(R) \) and the energy is given by the sum of the relative energy \( E_r \) and the center-of-mass energy \( E_R \), that is $$ \begin{equation*} E^{(2)}=E_r+E_R. \end{equation*} $$
We add then the repulsive Coulomb interaction between two electrons, namely a term $$ \begin{equation*} V(r_1,r_2) = \frac{\beta e^2}{|\mathbf{r}_1-\mathbf{r}_2|}=\frac{\beta e^2}{r}, \end{equation*} $$ with \( \beta e^2=1.44 \) eVnm.
Adding this term, the \( r \)-dependent Schroedinger equation becomes $$ \begin{equation*} \left( -\frac{\hbar^2}{m} \frac{d^2}{dr^2}+ \frac{1}{4}k r^2+\frac{\beta e^2}{r}\right)\psi(r) = E_r \psi(r). \end{equation*} $$ This equation is similar to the one we had previously in (b) and we introduce again a dimensionless variable \( \rho = r/\alpha \). Repeating the same steps as above, we arrive at $$ \begin{equation*} -\frac{d^2}{d\rho^2} \psi(\rho) + \frac{1}{4}\frac{mk}{\hbar^2} \alpha^4\rho^2\psi(\rho)+\frac{m\alpha \beta e^2}{\rho\hbar^2}\psi(\rho) = \frac{m\alpha^2}{\hbar^2}E_r \psi(\rho) . \end{equation*} $$ We want to manipulate this equation further to make it as similar to that in (a) as possible. We define a new 'frequency' $$ \begin{equation*} \omega_r^2=\frac{1}{4}\frac{mk}{\hbar^2} \alpha^4, \end{equation*} $$ and fix the constant \( \alpha \) by requiring $$ \begin{equation*} \frac{m\alpha \beta e^2}{\hbar^2}=1 \end{equation*} $$ or $$ \begin{equation*} \alpha = \frac{\hbar^2}{m\beta e^2}. \end{equation*} $$ Defining $$ \begin{equation*} \lambda = \frac{m\alpha^2}{\hbar^2}E, \end{equation*} $$ we can rewrite Schroedinger's equation as $$ \begin{equation*} -\frac{d^2}{d\rho^2} \psi(\rho) + \omega_r^2\rho^2\psi(\rho) +\frac{1}{\rho} = \lambda \psi(\rho). \end{equation*} $$ We treat \( \omega_r \) as a parameter which reflects the strength of the oscillator potential.
Here we will study the cases \( \omega_r = 0.01 \), \( \omega_r = 0.5 \), \( \omega_r =1 \), and \( \omega_r = 5 \) for the ground state only, that is the lowest-lying state.
With no repulsive Coulomb interaction you should get a result which corresponds to the relative energy of a non-interacting system. Make sure your results are stable as functions of \( \rho_{\mathrm{max}} \) and the number of steps.
We are only interested in the ground state with \( l=0 \). We omit the center-of-mass energy.
You can reuse the code you wrote for part (b), but you need to add the potential \( \omega_r^2\rho^2+1/\rho \) to the diagonal elements.
Comment the results for the lowest state (ground state) as function of varying strengths of \( \omega_r \).
For specific oscillator frequencies, the above equation has answers in an analytical form, see the article by M. Taut, Phys. Rev. A 48, 3561 (1993).
This exercise is a continuation of the previous and adds more quantum physics to the analysis. In this exercise we want to plot the wave function for two electrons as functions of the relative coordinate \( r \) and different values of \( \omega_r \). With no Coulomb interaction you should have a result which corresponds to the non-interacting case. Plot either the function itself or the probability distribution (the function squared) with and without the repulsion between the two electrons. Varying \( \omega_r \), the shape of the wave function will change.
We are only interested in the wave function for the ground state with \( l=0 \) and omit again the center-of-mass motion.
The eigenvectors are normalized. Plot then the normalized wave functions for different values of \( \omega_r \) and comment the results.
Here follows a brief recipe and recommendation on how to write a report for each project.
The preferred format for the report is a PDF file. You can also use DOC or postscript formats or as an ipython notebook file. As programming language we prefer that you choose between C/C++, Fortran2008 or Python. The following prescription should be followed when preparing the report: