\documentclass[12pt]{thesis}

\usepackage{pslatex}  % Use Times and Symbol fonts
\usepackage{amsmath}  % The amsmath package includes lots of useful
                      % mathematical stuff
\usepackage{graphicx}  % graphics stuff
\usepackage{dcolumn}  % So that the math column stuff works
\usepackage{array}    % More flexibility in table formats
\newcolumntype{d}[1]{D{.}{.}{#1}}
\newcolumntype{L}{>{$}l<{$}}
\newcolumntype{R}{>{$}r<{$}}
\newcolumntype{C}{>{$}c<{$}}
\usepackage{booktabs} % Better-looking tables
\usepackage{cite}     % Improves citation handling
\usepackage{graphicx} % For convenient inclusion of figures
\usepackage{psfrag}

\bibliographystyle{ss}

%  Any special definitions (\newcommand) that you need

\newcommand\kjm{{\text{kJ mol$^{-1}$}}}
\newcommand\ind{{\text{ind}}}
\newcommand\disp{{\text{disp}}}
\newcommand\tot{{\text{tot}}}
\newcommand\elec{{\text{elec}}}
\newcommand\rep{{\text{rep}}}
\newcommand\etc{{\textit{etc.}}}
\newcommand\mbf{{\mathbf}}

\newlength{\fntxvi} \newlength{\fntxvii}
\newcommand{\chemical}[1]
{{\fontencoding{OMS}\fontfamily{cmsy}\selectfont
  \fntxvi\the\fontdimen16\font
  \fntxvii\the\fontdimen17\font
  \fontdimen16\font=3pt\fontdimen17\font=3pt
  $\mathrm{#1}$
  \fontencoding{OMS}\fontfamily{cmsy}\selectfont
  \fontdimen16\font=\fntxvi \fontdimen17\font=\fntxvii}}

% Define \lt and \gt for `less than' and `greater than'
\mathchardef\lt="313C \mathchardef\gt="313E
% Set up `<' and `>' as angle brackets
\mathcode`<="4268 \mathcode`>="5269

\setlength{\floatsep}{0.4cm}
\setlength{\textfloatsep}{0.2cm}
\setlength{\intextsep}{0.2cm}
\renewcommand{\textfraction}{0.0}
\frenchspacing

\begin{document}

\title{Almost Exactly Solvable Models of Confined Interacting Electrons}
\author{David Thompson\\Churchill College, Cambridge}
\date{June 2001}

\maketitle

\frontmatter

This thesis submitted in June 2001 in partial fulfilment of the requirements for the Certificate of Postgraduate Study, describes work undertaken by myself in the theoretical section of the department of chemistry under the supervision of Dr.~A.~Alavi.  The work described is my own original research, except as indicated by references or otherwise.

I wish to thank my supervisor for his endless support and guidance throughout this year, my colleagues in room 246 and the EPSRC for financial support.

\newpage

\tableofcontents

\mainmatter

\chapter{Introduction}

The study of exactly or almost exactly solvable systems has always been a useful starting point for much of theoretical physics and chemistry.  The subject of Quantum Mechanics is introduced to a student through the solution of exact systems and the solutions to these models provide a useful first approximation for the description of modes of molecular motion. 

One particularly interesting class of exactly solvable problems are the confined electron problems.  In these, in the most basic example, we confine a single electron to a one-dimensional region of space through an infinite external potential.  Upon applying the equations of Quantum Mechanics to this problem we find that the solutions for the energy are quantised - we observe a discrete spectrum of energy levels where each level is labelled by an integer quantum number.  From this simple model we observe one of the most important phenomenon of quantum systems, that of quantization, it is this and problems like this one that highlight the importance of these model systems.

The work detailed in this thesis will be concerned with the solution of a particular confined electron problem.  I will set up a system such that I have an infinte external potential which surrounds a spherical region of space; inside this excluded volume I will have two interacting electrons. I will study the interaction of these electrons as a function of ${r_s}$ where ${r_s}$ is defined as 

\begin{equation}
\label{rs}
r_s=\left(\frac{3}{4n\pi}\right)^{\frac{1}{3}},
\end{equation}

and is the characteristic volume of space that contains one electron.  The variable ${n}$ is the density of the system and is linked to the total number of electrons and the volume such that for this particular problem \ref{rs} becomes

\begin{equation}
r_s=\frac{R}{2^{\frac{1}{3}}},
\end{equation} 

with ${R}$, the radius of the confining sphere.  Before begining a detailed analysis of this problem I wish to discuss some of the work which has been done in this field and to draw some links to some relevant and highly interesting physically relevant systems that can be described using these electron confinement models. 

\section{Quantum dots}

The understanding of how electrons interact in physical systems is of key importnace to chemists and physicists alike, and it has recently become possible to probe this interaction experimentally through the study of Quantum Dots.

Recent progress in Semiconductor manufacture has led to structures that are so small that they can contain as few as one electron and because of the properties of these devices the number of electrons within them can be varied exactly.  These devices have been dubbed quantum dots.  These quantum dot systems show such a range of interesting physics \cite{Ashoori96} \cite{VandV95} \cite{PalaciosJJ95} that the alternative terminology 'artificial atoms' has been adopted by many theoreticians and experimentalists; this is a reference to the close analogy of their properties to atomic systems.  In a normal atomic system the electrons are attracted to a central location  - the positive nucleus of the atom, in the case of the artificial atoms there is also an attraction towards a central location and it has been likened to a bowl-like parabolic confining potential where the electrons have the tendency of falling towards the bottom of the bowl.  

These artifical atoms are formed inside a semi-conducting material through quite ingenious techniques: for example it is possible using modern techniques in surface manufacture to form Gallium Arsenide (\chemical{GaAs}) one atomic layer at a time \cite{ChoAY91}.  It is also possible to replace some of the Gallium with Aluminium, thus forming Aluminium Gallium Arsenide (\chemical{AlGaAs}).  This material acts as an insulator, but electrons are free to move in the \chemical{GaAs}; sandwiching a 10nm thickness of \chemical{GaAs} between \chemical{AlGaAs} layers one is able to confine electrons in a \chemical{GaAs} "quantum well".  The electron can move laterally within the well, but have no freedom to move perpendicularly to it.  By placing electrostatic gates at the surface of such a structure we are laterally confining a two-dimensional electron gas and creating a quantum dot.  Through variation of the potential across such a device we are able to control the number of electrons confined in the quantum well.

One of the key physical characteristics of these dots and what led them to being dubbed artificial atoms is that their size dominates much of their physics; these artificial atoms are typically much larger (${\sim}$100nm on a side) than natural atoms, and the electron orbits do not scale simply with the size.  Consider an atom containing many electrons which has a size that varies continuously.  At large sizes the Coulomb energy is low as the average separation between electrons is large, however there is another energy scaling that one must consider - that of the spacing of the different electron orbital energies.  In the artificial atom, as the atomic size increases the difference in orbital energy decreases, and it is found to decrease much faster than the Coulomb energy.  Thus in a large atom the effects of electron-electron interaction become much more important than in a small atom. 
 
How then does one go about modeling such systems?  As is the way with the theoretical sciences we seek the simplest system that best resembles our own, in this
 case the simplest example of our problem would be a single electron confined to move in the ${x}$,${y}$ direction in a spherically symmetric parabolic potential ${1/2m\omega_0r^2}$ with ${\omega_0}$ the frequency of the oscillation around the centre 
of the parabolic bowl, and ${r^2=x^2+y^2}$ the distance from the centre of the dot.  The eigenstates of this problem were studied early on by Fock \cite{FockV28} and Darwin \cite{DarwinCG30} and it was found that the corresponding energies of the system had the form

\begin{equation}
E_{nl}=\hbar\omega_0\left(2n+|l|+1\right).
\end{equation}

The two quantum numbers ${n}$ and ${l}$ are indicative of the motion of the electron in two dimenions;  ${n}$ is a positive integer and is linked to the radial solution to the problem, ${l}$ can be any integer, positive or negative.  It is observed that the ${l}$ quantum numbers have circular orbits associated with them; the larger the value of ${|l|}$ the larger the radius of the orbit.  It should be noted that in the absence of a magnetic field both orientations of the spin of the electron have the same energy.

This small piece on the basic experimental and theoretical background of quantum dots should serve as part introduction to the modeling of confined electron systems  and also as part justification of this work.  It is important to realise that now there do exist experimental systems such that electron-electron interaction can (to a good approximation) be studied exactly.

%
\chapter{Theory}
%
Here we shall begin to describe how the theory of two confined interacting electrons is formulated; it is easiest to begin with the theory of one electron in a sphere and to use these solutions as a basis for the two electron problem.

%
\section{The one electron problem}
%

Consider a single electron confined to a spherical region of space by a potential such that

\begin{equation}
\label{sphpot}
V(r)=
\begin{cases}
0& \text{$r < R$} \\
\infty& \text{$r > R$}
\end{cases}.
\end{equation}

With $R$ the radius of the sphere.  The corresponding Schr\"odinger equation for this problem has the form 

\begin{equation}
\label{1schr}
\left\{-\frac{1}{2}\nabla^2+V(r)\right\}\Psi(\mathbf{r})=E\Psi(\mathbf{r}),
\end{equation}

with $V(r)$ as \ref{sphpot}.  We look for solutions to \ref{1schr} which are separable in the natural co-ordinates of the problem - in this case spherical polar such that 

\begin{equation}
\Psi(\mathbf{r})=R(r)\Theta(\theta)\Phi(\phi).
\end{equation}

We can see that the potential $V(r)$ will only enter into the radial equation so for the azimuthal angle,$\phi$, we find an equation for $\Phi(\phi)$ which has the form

\begin{equation}
\label{phi}
\frac{\text{d}^2\Phi}{\text{d}\phi^2}+m^2\Phi=0.
\end{equation}

We denote $m$, a constant, the orbital angular momentum quantum number.  It labels the solutions to \ref{phi} such that  

\begin{equation}
\Phi_m(\phi)=\frac{1}{\sqrt{2\pi}}e^{im\phi}.
\end{equation}
%

For well behaved $\Phi_m(\phi)$ we see that $m=0,\pm1,\pm2,\cdots$.  The corresponding eqution for $\Theta(\theta)$ has the form


\begin{equation}
\label{theta}
\frac{1}{\sin(\theta)}\frac{\text{d}}{\text{d}\theta}\left(\sin(\theta)\frac{\text{d}\Theta}{\text{d}\theta}\right)+\left[l(l+1)-\frac{m^2}{{\sin}^2(\theta)}\right]\Theta=0.
\end{equation}

We find solutions to \ref{theta} for only certain values of $l$.  It is found that $l=0,1,2\cdots$ and that this places restrictions on the allowed values of $m$, thus $m$ can only take the values $m=0,\pm1,\pm2\cdots\pm\l$.  

The products of $\Theta_{lm}(\theta)$ and $\Phi_m(\phi)$ are known as the Spherical Harmonics and are labelled for any given $\theta,\phi,m$ and $l$ as $Y_{lm}(\theta,\phi)$; indeed these are the angular solutions to the Schr\"odinger equaton for any central potential $V(r)$.

Having solved the angular part of the problem all that is left to tackle is the radial equation, which must be solved for each $V(r)$.  This equation has the form

\begin{equation}
\label{radial}
-\frac{1}{2}\left(\frac{\text{d}^2R}{\text{d}r^2}+\frac{2}{r}\frac{\text{d}R}{\text{d}r}\right)+\left[V(r)+\frac{l(l+1)}{2}\right]R=ER.
\end{equation}

The $l(l+1)$ term is physically interpreted as an addition to the potential that keeps the electrons from the origin when $l$ is greater than zero.  In solving this equation we consider the form of the potential given by \ref{sphpot} and also the boundary condition

\begin{equation}
R(r)=0 \qquad r>R
\end{equation}

where we pay attention to note that $R$ is formally the radius of the confining sphere and that $R(r)$ is the radial function of $r$ that is being determined.  Given these conditions it is trivial to find the corresponding $R(r)$ and we relate them to the well known functions of mathematical physics - the spherical Bessel functions, denoted $j_l(\alpha_{nl}r)$.  The $\alpha_{nl}$ is such that

\begin{equation}
\label{root}
\alpha_{nl}=\frac{k_{nl}}{R},
\end{equation}

with $k_{nl}$ the $n^{\text{th}}$ root of the $l^{\text{th}}$ order spherical Bessel function. It is important to realise that $n$ is a new quantum number introduced upon solving \ref{radial} but that $l$ is the same quantum number as was found from solving \ref{theta}.  Thus the final solution for $\Psi(\mathbf{r})$ becomes

\begin{equation}
\label{unbasis}
\Psi(\mathbf{r})=j_l(\alpha_{nl}r)Y_{lm}(\theta,\phi).
\end{equation} 

There are some interesting features that should be noted about these solutions, the first is that we observe degenerate energy levels as $\Psi(\mathbf{r})$ with different $m$ quantum number all have the same energy.  In fact it is trivial to see that as $m$ is restricted to $m=0,\pm1,\pm2\cdots,\pm\l$ there are $2l+1$ possible $Y_{lm}(\theta,\phi)$ for any given $l$, which means that the corresponding degeneracy of each orbital is $2l+1$.  We observe a similar degeneracy to electron orbitals in atoms.

Having found the eigenstates of a single electron confined to a spherical region we can construct the eigenstates for the two interacting electrons as a linear combination of the one electron basis.  It is to the two electron problem we now turn.

\section{The two electron problem}

It should be noted that we follow, in methodology at least, the work of Dr. A. Alavi \cite{AlaviA00} who has recently solved the analogous problem of two interacting electrons confined to a three dimensional box.

Consider a two body wavefunction $\Psi(\mathbf{x_1,x_2})$, with $\mathbf{x}=(\mathbf{r}\sigma)$, which we can write in terms of a product of a spatial and spin wavefunction as

\begin{equation}
\label{utwobody}
\Psi(\mathbf{x_1,x_2})=\Psi^\pm(\mathbf{r_1,r_2})\chi^\mp(\sigma_1,\sigma_2).
\end{equation}

We associate the two spin states, singlet ($\chi^-$ : S=0) and triplet ($\chi^+$ : S = 1) with the symmetric ($\Psi^+$) and anti-symmetric ($\Psi^-$) spatial wavefunctions respectively. We can then expand the spatial wavefunctions $\Psi^\pm(\mathbf{r_1,r_2})$ in terms of \ref{unbasis} such that

\begin{equation}
\label{twobody}
\Phi^\pm(\mathbf{r_1,r_2})=N_{nlmpqr}\left\{\eta_{nlm}(\mathbf{r_1})\eta_{pqr}(\mathbf{r_2})\pm\eta_{pqr}(\mathbf{r_1})\eta_{nlm}(\mathbf{r_2})\right\}. 
\end{equation}

Where $\eta_{nlm}$ has the form 

\begin{equation}
\label{basis}
\eta_{nlm}=N_{nl}j_l(\alpha_{nl}r)Y_{lm}(\theta,\phi),
\end{equation}

with

\begin{equation}
\label{norm}
N_{nl}=\frac{2}{R^3}\left(j_{l+1}(k_{nl})\right)^{-1},
\end{equation}

and 

\begin{equation}
\label{Ng}
N_{nlmpqr}=\sqrt\frac{1}{\sqrt{2+2\delta_{np}\delta_{lq}\delta_{mr}}}.
\end{equation}

Equation \ref{norm} is simply a normalisation factor such that

\begin{equation}
\iint\eta_{nlm}(\mathbf{r_1})\eta_{pqr}(\mathbf{r_2})\,\text{d}\mathbf{r_1}\,\text{d}\mathbf{r_2}=\delta_{np}\delta_{lq}
\end{equation}

The notation we are using can be simplified somewhat if we think of the six quantum numbers $n,l,m,p,q$ and $r$ which specify any of our eigenstates completely, as a six dimensional vector $\mathbf{G}$, such that

\begin{equation}
\mathbf{G}=\left\{n l m p q r\right\}.
\end{equation}

Using this new vector notation for the quantum numbers we can thus recast equation \ref{utwobody} in a rather compact linear form such that the two body wavefunction can be written as a linear combination of the spatial wavefunctions 

\begin{equation}
\label{wavefunc}
\Psi^\pm(\mathbf{r_1,r_2})=\sum_{\mathbf{G}}c^\pm_{\mathbf{G}}\Phi^\pm_{\mathbf{G}}(\mathbf{r_1,r_2}).
\end{equation}

The $c^\pm_{\mathbf{G}}$ are real and must be determined, the $\Phi^\pm_{\mathbf{G}}(\mathbf{r_1,r_2})$ are a complete set and are orthonormal, and the $\mathbf{G}$ run over all possible positive integer combinations.

Havin established the most general basis we now turn to the evaluation of the eigenstates of the problem.  The two electron Hamiltonian for this problem is

\begin{equation}
\label{hamilt}
\hat{H}=\hat{T}+\hat{U}+\hat{V}.
\end{equation}

Where 

\begin{equation}
\label{kinetic}
\hat{T}=-\frac{1}{2}\sum_{i=1}^2\nabla_i^2,
\end{equation}

\begin{equation}
\label{indiv}
\hat{U}=\frac{1}{|\mathbf{r_1-r_2}|}, 
\end{equation}

and

\begin{equation}
\label{indivpot}
\hat{V}=\sum_V(\mathbf{r_i}). 
\end{equation}

We identify the $\hat{V}$ in \ref{indivpot} with the potential described by \ref{sphpot}.  We see that in solving our problem  we are going to have to calculate Hamiltonian matrix elements and looking at \ref{hamilt} we see that this will reduce to calculating individual contributions from the kinetic and Coulombic terms, given through \ref{kinetic} and \ref{indiv} respectively.  The kinetic term is reasonably straightforward and it can readily be shown that it is diagonal in the basis such that

\begin{equation}
\langle\Phi_\mathbf{G}|\hat{T}|\Phi_\mathbf{G'}\rangle=\frac{1}{2}\left(\alpha_{nl}^2+\alpha_{pq}^2\right)\delta_{\mathbf{GG'}}
=\frac{1}{2R^2}\left(k_{nl}^2+k_{pq}^2\right)\delta_{\mathbf{GG'}}, 
\end{equation}

the $k$'s in the above expression are the same as those given by \ref{root} and are the roots of the spherical Bessel functions.  The Coulomb term requires additional thought and it is to this that we shall devote a little time.

%
\section{Coulomb interaction}
%

From \ref{hamilt} we see that the form of the Coulomb matrix element to be evaluated is of the form

\begin{equation}
\label{coulint}
\hat{U}_\mathbf{GG'}=\iint\frac{\eta_{nlm}^*(\mathbf{r_1})\eta_{pqr}^*(\mathbf{r_2})\eta_{nlm}(\mathbf{r_1})\eta_{pqr}(\mathbf{r_2})}{|\mathbf{r_1-r_2}|}\text{d}\mathbf{r_1}\text{d}\mathbf{r_2}.
\end{equation}

This is a non trivial integral to evaluate.  As it is written above it is a six dimensional expression which will take an appreciable time to calculate.  In order to make it's evaluation slightly more feasible we have to find some way to simplify this expression.  If we think about the form of the basis used in the above expression we see that we can use the addition theorem for spherical harmonics   

\begin{equation}
\label{addthe}
\frac{1}{|\mathbf{r_1-r_2}|}=4\pi\sum_{l=0}^\infty\sum_{m=-l}^l\frac{1}{2l+1}\frac{r_1^l}{r_2^{l+1}}Y_{lm}^*(\theta_1,\phi_1)Y_{lm}(\theta_2,\phi_2); \qquad r1<r2.
\end{equation}

Inserting this into \ref{coulint} we find that we can separate the integral into a product of angular and radial terms

\begin{align}
\label{totint}
\begin{split}
\hat{U}_\mathbf{GG'}
&=4\pi\sum_{l''=0}^\infty\sum_{m''=-l''}^{l''}\frac{1}{2l''+1}\int Y_{lm}^*(\theta_1,\phi_1)Y_{l'm'}(\theta_1,\phi_1)Y_{l''m''}^*(\theta_1,\phi_1)\text{d}\Omega_1 \times \\
&\int Y_{qr}^*(\theta_2,\phi_2)Y_{q'r'}(\theta_2,\phi_2)Y_{l''m''}^*(\theta_2,\phi_2)\text{d}\Omega_2 \times \\
&N_{nl}N_{n'l'}N_{pq}N_{p'q'} \times \\
&\biggl[\int_{r_2=0}^R \int_{r_1=0}^{r_2} j_l^*(\alpha_{nl}r_1)j_q^*(\alpha_{pq}r_2)j_{l'}(\alpha_{n'l'}r_1)j_{q'}(\alpha_{p'q'}r_2)
\frac{r_1^{l''}}{r_2^{l''+1}}
r_1^2 r_2^2\text{d}r_1\text{d}r_2 + \\ 
&\int_{r_1=0}^R \int_{r_2=0}^{r_1} j_l^*(\alpha_{nl}r_1)j_q^*(\alpha_{pq}r_2)j_{l'}(\alpha_{n'l'}r_1)j_{q'}(\alpha_{p'q'}r_2) 
\frac{r_2^{l''}}{r_1^{l''+1}}
r_1^2 r_2^2\text{d}r_1\text{d}r_2\biggr]. 
\end{split}
\end{align}

Having separated out the co-ordinates we look for analytic solutions to the separate angular and radial parts.  Looking firstly at the angular product, we find that what we are looking at is an analytic expression for the evaluation of the integral of products of spherical harmonics.  We find \cite{BrinkDM68} \cite{Zare88} that a useful relationship exists between the so-called normalised spherical harmonics and the Wigner-3j symbols such that

\begin{equation}
\int C_{a\alpha}(\theta,\phi) C_{b\beta}(\theta,\phi) C_{c\gamma}(\theta,\phi)\text{d}\Omega = 4\pi
\begin{pmatrix} a & b & c \\ \alpha & \beta & \gamma\end{pmatrix}
\begin{pmatrix} a & b & c \\ 0 & 0 & 0\end{pmatrix},
\end{equation}

where there exists a tringular inequality between $a,b$ and $c$ such that if they do not obey

\begin{equation}
|a-c|\le b \le |a+c|
\end{equation}

and also 

\begin{equation}
\alpha+\beta+\gamma = 0
\end{equation}

the Wigner-3j symbol is identically zero.  The $C_{ij}(\theta,\phi)$ are linked to the ordinary spherical harmonics through the relationship

\begin{equation}
\label{simp1}
Y_{lm}(\theta,\phi)=\left(\frac{2l+1}{4\pi}\right)^{\frac{1}{2}} C_{lm}(\theta,\phi).
\end{equation}

Another useful relationship is one that exists for the ordinary spherical harmonics and relates the complex conjugate to the real form

\begin{equation}
\label{simp2}
Y_{lm}^*(\theta,\phi)=\left(-1\right)^m Y_{l-m}(\theta,\phi).
\end{equation}

Using \ref{simp1} and \ref{simp2} we can integrate out the angular dependence in \ref{totint}.  Doing this we find that \ref{totint} becomes

\begin{align}
\label{coulintfin}
\begin{split}
\hat{U}_\mathbf{GG'}
&=\sum_{l''=0}^\infty\sum_{m''=-l''}^{l''}(-1)^{m+r+2m''}
\left(2l''+1\right)\left[(2l+1)(2l'+1)(2q+1)(2q'+1)\right]^\frac{1}{2} \times \\
&\begin{pmatrix} l & l' & l'' \\ -m & m' & -m''\end{pmatrix}
\begin{pmatrix} q & q' & q'' \\ -r & r' & -r''\end{pmatrix}
\begin{pmatrix} l & l' & l'' \\ 0 & 0 & 0\end{pmatrix}
\begin{pmatrix} q & q' & q'' \\ 0 & 0 & 0\end{pmatrix} \times \\ 
&N_{nl}N_{n'l'}N_{pq}N_{p'q'} \times \\
&\left[2\int_{r_2=0}^R \int_{r_1=0}^{r_2} j_l^*(\alpha_{nl}r_1)j_q^*(\alpha_{pq}r_2)j_{l'}(\alpha_{n'l'}r_1)j_{q'}(\alpha_{p'q'}r_2)
\frac{r_1^{l^{''}}}{r_2^{l^{''}+1}}
r_1^2 r_2^2\text{d}r_1\text{d}r_2\right]
\end{split}
\end{align}

The factor of two in front of the radial term indicates that the sum of the terms coming from the expressions for the $r_1<r_2$ and $r_2<r_1$ are the same.  We have made some important progress, the original six dimensional given by \ref{totint} has been reduced to a two dimensional integral.  Also due to the high symmetry of the Wigner-3j symbols we expect a reasonable number of the Coulomb matrix elements to be zero.  In practice we find that although this simplifies the calculation and indeed makes it practicable for the system being studied, the majority of the computational time is spent on the evaluation of the Coulomb matrix elements and may be a limiting factor when extending the system to more electrons.  In practice, code was written such that the double radial integral is evaluated numerically using a Gauss-Chebyshev \cite{NexCMM79} quadrature over a radial mesh of 512 points.  The inner integral which is formally a function of the upper limit of integration $r_2$ is calculated once and for all at the beginning of the calculation and the results are stored in an array over the values of $r_2$.  The outer integral is then calculated using these tabulated values which it calls as and when it needs them.
The Gauss-Chebyshev method is employed as the pivots for the numerical integration are at the roots of trignometric functions, which is useful given that we know the spherical Bessel functions are closely linked to the trigonometric functions. 

In order to determine the effectiveness of the algorithm we have developed we need alternative methods of calculating some of the elements and comparing.  If we consider the integral $\hat{U}_{\mathbf{GG'}}$ with $\mathbf{G=G'}$ and $\mathbf{G}=\left\{100100\right\}$ we find that the integral is analytic. Also a direct "brute force" method was written.  We thus have three different numerical answers which should in theory  be the same.  

\begin{table}
\begin{center}
\caption{Comparison of the numerical and analytic values of the $\langle100100|\hat{U}|100100\rangle$ Coulomb matrix element}
\label{compint}
\begin{tabular}{CCd9} \\ \toprule
\text{Method of Evaluation} & \multicolumn{1}{c}{Value} \\ \toprule 
  Gauss-Chebyshev & 1.786073080 \\
Direct evaluation & 1.786176860 \\ 
Analytic Solution & 1.786073450 \\ \bottomrule
\end{tabular}  
\end{center}
\end{table} 

We see that the Gauss-Chebyshev method differs from the analytic solution by $2.07\times10^{-5}\%$, well within the numerical accuracy needed for the evaluation of these matrix elements.

\begin{figure}
\begin{center}
\includegraphics[scale=0.40,angle=0]{meshfig.pstex}
\caption{Convergence of direct evaluation as a function of mesh size}
\label{fig:mesh}
\end{center}
\end{figure}

Figure \ref{fig:mesh} below shows the convergence of the brute force method as a function of the mesh-size of the radial grid.  It should be noted that the time to do these calculations increased with increasing mesh size, such that the calculation of the 179  mesh point integral took several minutes.  It is for this reason, and that the direct evaluation was extremely slowly convergent, that the Gauss-Chebyshev approach was used, the Gauss-Chebyshev method giving excellent agreement and being extremely quick, with integrals being evaluated in hundredths of seconds.

Having now described how to succesfully calculate both the kinetic and Coulombic matrix elements we are in a position to solve the system completely.  We can write an expression for the energy of the system as an eigenvalue expression in the $c_{\mathbf{G}}$ such that

\begin{equation}
\sum_{\mathbf{G}}\hat{H}_{\mathbf{GG'}}c_{\mathbf{G'}}^\pm=E^\pm c_{\mathbf{G}}^\pm.
\end{equation}

This is in principle exact, with the limit on accuracy coming from the number of basis functions used in evaluating $\Psi$, the two-body wavefunction.  The energies and the convergence of energies as a function of $\mathbf{G}$ are discussed in the next chapter.

The next property of the system that is of extremely important interest is the ground state density.  It is worthwhile spending a reasonable amount of time on the calculation of the density for this system as there are some fascinating insights into the form of the correlation between the electrons as a function of the size of the system.

\section{Calculation of the ground state density} 

One of the most interesting properites of this two electron system is that the density functional theory (DFT), when applied to it, can be greatly simplified.  We know that for a two electron system we can write

\begin{equation}
\label{density}
n(\mathbf{r_1})=2\int |\Psi^\pm(\mathbf{r_1,r_2})|^2 \text{d}\mathbf{r_2}.
\end{equation}

With $n(\mathbf{r_1})$ being the electron density of particle one at position $\mathbf{r_1}$.  Now we can substitute \ref{wavefunc} into \ref{density} such that 

\begin{equation}
\label{odense}
n(\mathbf{r_1})=2\sum_{\mathbf{GG'}} c_{\mathbf{G}} c_{\mathbf{G'}} N_{\mathbf{G}} N_{\mathbf{G'}} S_{\mathbf{GG'}},
\end{equation}

the $N_{\mathbf{G}}$ are normalisation coefficients given by \ref{Ng}, the $c_{\mathbf{G}}$ are clearly the wavefunction coefficients and

\begin{equation}
S_{\mathbf{GG'}}=\int \Phi_{\mathbf{G}}(\mathbf{r_1,r_2})\Phi_{\mathbf{G'}}(\mathbf{r_1,r_2}) \text{d}\mathbf{r_2}. 
\end{equation}

Upon calculating the density through \ref{density} we find that it is spherically symmetric - although after a little thought it is easy to see this result;  We know that the density is a scalar quantity and is thus independent upon our choice of axes, it is also independent upon rotations about the chosen axes.  Given that this is the case we are able to integrate out the angular dependance of the density.  Formally we integrate \ref{odense} with respect to the solid angle $\Omega_1$ such that

\begin{equation}
\begin{split}
n(r_1)&=\int n(\mathbf{r_1})\text{d}\Omega_1 \\
&=2\sum_{\mathbf{GG'}} c_{\mathbf{G}} c_{\mathbf{G'}} N_{\mathbf{G}} N_{\mathbf{G'}} \int S_{\mathbf{GG'}} \text{d}\Omega_1 \\
&=\frac{1}{2\pi}\sum_{\mathbf{GG'}} c_{\mathbf{G}} c_{\mathbf{G'}} N_{\mathbf{G}} N_{\mathbf{G'}} \int S_{\mathbf{GG'}} \text{d}\Omega_1. 
\end{split}
\end{equation}

From a computational point of view we find that the density is simply a sum of products of spherical Bessel functions with integrals involving spherical Bessel functions and Kronecker deltas that depend upon the angular momentum (and projection) quantum numbers.

We calculate  the integrals over all possible combinations of the quantum numbers, over a mesh of 512 radial points which we store as an array.  We do this once and for all at the initialization of the density routine and call elements from the array as and when we need them.  The fact that the density is also just a function of one co-ordinate is extremely advantageous and the calculation of the density for this system is by far the quickest part of the whole calculation for any given basis set.  Having a density that is extremely easy to calculate is an attractive feature of the system we are studying, as all properties of the system can be derived from the density, if we know the density everything about the system is known.

\section{The exchange-correlation hole}

One of the interesting properties that is amenable to calculation once the density is known is the physical exchange-correlation hole.  Formally in DFT this is defined in terms of the density-density correlation function evaluated with the two body wavefunction \ref{utwobody} such that

\begin{equation}
\label{xchole}
n_{xc}(\mathbf{r_1,r_2})=\frac{1}{n(\mathbf{r_1})}\langle\Psi|\delta\hat{n}(\mathbf{r_1})\delta\hat{n}(\mathbf{r_2})|\Psi\rangle-\delta(\mathbf{r_1-r_2}).
\end{equation}

Where

\begin{equation}
\delta\hat{n}(\mathbf{r})=\hat{n}(\mathbf{r})-n(\mathbf{r})
\end{equation}

and $\hat{n}(\mathbf{r})$ is the number density operator such that,

\begin{equation}
\hat{n}(\mathbf{r})=\sum_i \delta(\mathbf{r-r_i}).
\end{equation}

As we have seen already for the density, these expressions in their most general form in DFT, can be greatly simplified for the two electron system being studied.  We find that the physical exchange-correlation hole given by \ref{xchole} can be written solely in terms of the density and the two-body wavefunction 

\begin{equation}
n_{xc}(\mathbf{r_1,r_2})=\frac{2|\Psi(\mathbf{r_1,r_2})|^2}{n(\mathbf{r_1})}-n(\mathbf{r_2}).
\end{equation}  

We find that the exchange correlation hole given by the expression above is normalised to -1.  Having this expression for the exchange-correlation hole means that we are able to investigate both the exchange and correlation components independently since we know that we can write

\begin{equation}
n_{xc}(\mathbf{r_1,r_2})=n_x(\mathbf{r_1,r_2})+n_c(\mathbf{r_1,r_2}).
\end{equation}

Furthermore we know that for a two electron system the $n_x(\mathbf{r_1,r_2})$ becomes

\begin{equation}
n_x(\mathbf{r_1,r_2})=-\frac{n(\mathbf{r_2})}{2},
\end{equation}

from which we can see that the correlation term has the form

\begin{equation}
\label{chole}
n_c(\mathbf{r_1,r_2})=\frac{2|\Psi(\mathbf{r_1,r_2})|^2}{n(\mathbf{r_1})}-\frac{n(\mathbf{r_2})}{2}.
\end{equation}

It is readily seen that the correlation hole is normalised to zero.  Now what does the expression \ref{chole} mean and how do we go about calculating it?  In practice we fix the position of one of the electrons, for example $r_1$, and then move the other electron over the whole $r_2$ space.  We will thus generate a function which tells us how likely it is that there will be an electron at $r_2$ given that there is one at $r_1$.  Now as was the case with the density we are able to use the high symmetry of the system to simplify our calculations tremendously.
As it is written \ref{chole} is a six-dimensional function, but upon fixing the position of electron one to be some radial distance $r_1$ from the origin (in a cartesian representation of the system) then we can allow the second electron to be at a distance $r_2$ from the origin and at an angle $\theta$ to $r_1$.  This describes the whole system and reduces it to a two variable problem as we vary only $\theta$ and $r_2$.  This simplification is due to the symmetry of the problem, reflected by the fact that the density is a function only of the radial distance.  Formally the form of \ref{chole} is 

\begin{equation}
n_{xc}(r_1,0,0;r_2,\theta,0)=\frac{2|\Psi(r_1,0,0;r_2,\theta,0)|^2}{n(r_1)}-n(r_2),
\end{equation}

$0 \le r_{1,2} \le R$ and $0 \le \theta \le \pi$.  Figure \ref{fig:xchole} gives a geometrical interpretation of what we are doing; it may also be a useful diagram for the reader to consult when looking at different plots of \ref{chole} when it comes to a discussion of the results.

\begin{figure}
\begin{center}
\includegraphics[scale=0.60,angle=0]{xc.eps}
\caption{Figure showing the position of the two electrons in calculating the physical exchange correlation hole}
\label{fig:xchole}
\end{center}
\end{figure}

\section{The Kohn-Sham potential}

We also realise that in DFT we can relate the density to the Kohn-Sham (KS) one electron orbitals such that for this two electron system

\begin{equation}
\label{ks}
n(r)=2\left(\varphi(r)\right)^2.
\end{equation}

The $\varphi(r)$ are the KS orbitals that satisfy the KS equation

\begin{equation}
\label{kseqn}
\left[-\frac{1}{2}\nabla^2 + v_s(r)\right]\varphi(r)=\epsilon\varphi(r).
\end{equation}

As we know that we can invert \ref{ks} and find the KS orbitals in terms of the density - a quantity we know - we can also invert \ref{kseqn} such that we can calculate the KS potential $v_s(r)$.  Physically this corresponds to the calculation of the local potential that the non-interacting electrons move in.  We find that upon inversion of \ref{kseqn} 

\begin{equation}
\label{vsr}
v_s(r)=\epsilon + \frac{1}{2}\frac{\nabla^2 n^{\frac{1}{2}}(r)}{n^{\frac{1}{2}}(r)}.
\end{equation} 

As we know the density we can calculate the KS potential to within an additive constant.

\chapter{Results}

In this section the results collated thus far will be shown and in the next chapter a detailed discussion of the results will be presented along with aims for research in the future - both short and long term.

\section{Energy convergence}

In beginning an energy calculation we must specify a basis set size.  This is done through the variation of input parameters in an input file that is read through into the main body of code.  These two input paramters we term $n_{max}$ and $l_{max}$, where the $max$ refers to the maximum values of the quantum numbers $n$ and $l$ as defined by \ref{basis}.  As a word on nomenclature I will refer to a given basis as $\left\{n_{max},l_{max}\right\}$.  As an example of what this notation means in terms of the basis, consider $\left\{{2,1}\right\}$.  In the expansion \ref{wavefunc} we find that $\mathbf{G}$ runs over all combinations of integers such that $n_{max}=2$ and $l_{max}=1$.

Below figure \ref{fig:unitconv} shows the convergence of the energy as we vary the orbital angular momentum quantum number, $l$.  Along any given line the value of $n$ is constant.  In figure \ref{fig:lunitconv}, a plot of energy convergence as a function of the radial quantum number $n$, we see that the final converged energy in the case of the unit sphere is 11.592524685 a.u.  Similar calculations were also done for spheres of radii ten and fifty, the graphs which illustrate their convergence as a function of $l$ are shown in figures \ref{fig:tenconv} and \ref{fig:fiftyconv} respectively.

In all three cases we see, what could be termed, a two dimensional convergence.  In all of the graphs we see that the system saturates very quickly with respect to the amount of angular information needed for any constant $n$.  To illustrate this point we see that in all cases for $l>4$ at all sphere sizes considered there is very little gain in adding additional angular information in the basis; the system is very well converged in the angular sense, as is shown quantitatively in figure \ref{fig:lunitconv}.  It appears that the radial information in the basis is the significant contributor to the convergence of the energy.

Looking at the form of the convergence we can consider a  converged basis to be of the size $\left\{5,8\right\}$.  We can qualify this statement if we look at some individual data points of the figures \ref{fig:unitconv},\ref{fig:tenconv} and \ref{fig:fiftyconv}.  These are collated in table \ref{convsph}.  We see that in all cases the relative errors on approaching $\left\{5,8\right\}$ decrease markedly.

It is realised from physical arguments (for example see section on Quantum dots) that the energy will be dominated at different $r_s$ by either the kinetic or the Coulombic components.  In fact it is easy to see that at small $r_s$ the energy goes as $\sim\frac{A}{r_s^2}$, whilst at large $r_s$ as $\sim\frac{B}{r_s}$, with $A$ and $B$ constants.  Using the converged basis, $\left\{5,8\right\}$, we have calculated the ground state energy for spheres of various sizes (1,5,15,25,35,45,50) and have plotted $E.r_s^2$ against $r_s$.  This is shown in figure \ref{fig:ers2}.  As the energy has the form

\begin{equation}
\label{Eapprox}
E\sim\frac{A}{r_s^2}+\frac{B}{r_s}+\cdots,
\end{equation} 

we see that in the limit of $r_s\rightarrow 0$, $E.r_s^2$ should tend to a constant - as is shown in the figure.

\begin{figure}
\begin{center}
\includegraphics[scale=0.40,angle=0]{Unitconv.pstex}
\caption{Convergence of energy as a function of orbital angular momentum quantum number l, for the unit sphere}
\label{fig:unitconv}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.40,angle=0]{lunitconv.pstex}
\caption{Convergence of energy as a function of radial quantum number n, for the unit sphere}
\label{fig:lunitconv}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.40,angle=0]{tenconv.pstex}
\caption{Convergence of energy as a function of radial quantum number n, for the sphere of radius ten}
\label{fig:tenconv}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.40,angle=0]{fiftyconv.pstex}
\caption{Convergence of energy as a function of radial quantum number n, for the sphere of radius fifty}
\label{fig:fiftyconv}
\end{center}
\end{figure}

\begin{table}
\begin{center}
\caption{Converged energies as a function of sphere size}
\label{convsph}
\begin{tabular}{cccCCCd9} \\ \toprule
Radius of sphere & Basis set & No. of basis functions & \multicolumn{1}{r}{Energy (a.u.)} & \multicolumn{2}{r}{Relative error} \\ \midrule
	1 & 2,1 & 12 & 11.59948053 & 0.06\% \\
          & 3,3 & 96 & 11.59347065 & 0.008\% \\
          & 4,5 & 360 & 11.592748512 & 0.002\%  \\
          & 5,8 & 1215 & 11.59252468 &       \\ \midrule
       10 & 2,1 & & 0.243021280 & 0.109\% \\
          & 3,3 & & 0.242804743 & 0.020\% \\
          & 4,5 & & 0.242768237 & 0.005\% \\ 
          & 5,8 & & 0.242756465 &       \\ \midrule
       50 & 2,1 & & 0.028355601 & 3.786\% \\
          & 3,3 & & 0.027331533 & 0.040\% \\
          & 4,5 & & 0.027325240 & 0.015\% \\  
          & 5,8 & & 0.027321140 &  
\\ \bottomrule
\end{tabular}
\end{center}
\end{table} 

\begin{figure}
\begin{center}
\includegraphics[scale=0.40,angle=0]{ers2.pstex}
\psfrag{Ers2}{$E.r_s^2$}
\psfrag{rs2}{$r_s^2}
\caption{Plot of $E.r_s^2$ against $r_s$}
\label{fig:ers2}
\end{center}
\end{figure}

\section{Ground state density}

Given the expressions detailed in the previous chapter the ground state density was investigated for different values of $r_s$.  It was found that whilst writing the code, for debugging purposes, it was best to work with a very small basis, so all the code was tested for $\left\{2,1\right\}$ and along the series $\left\{2,n\right\}$ with $n=1,8$.  This is a reasonable approach to follow as we know that densities will have the same overall shape, we show below in figure \ref{fig:unitdensity} the shape of the density for the basis sets $\left\{2,n\right\}$, with ${n=1,2,4,8}$ for the unit sphere.  We observe that for all different basis set sizes there is no change in the shape of the density.

\begin{figure}
\begin{center}
\includegraphics[scale=0.40,angle=0]{unitdensity.pstex}
\caption{Plot of $n(r).r_s^3$ against radial distance from origin of the unit sphere for differing basis set sizes} 
\label{fig:unitdensity}
\end{center}
\end{figure}

In figure \ref{fig:wigone} the scaled  density is plotted for three different sphere sizes R=9, 10 and 11 with the smallest basis $\left\{2,1\right\}$, correspondingly figure \ref{fig:wigoneone} shows the scaled density for the sphere of radius 10 for the basis $\left\{2,1\right\}$ and $\left\{2,8\right\}$ - we observe several interesting things from these results.  Firstly, that upon going along the radii 9, 10 and 11 we see the development of a minimum in the electron density at the origin and secondly, that upon using a larger basis set the minimum becomes less pronounced although still evident.  From the latter point one could conclude that even at the lowest basis the density is qualitatively correct in form and that the true converged density will be very similar to the densities illustrated here.  The former point is an interesting consequence of the physics of the system and the formation of the minimum has been dubbed the onset of "Wigner crystallisation," in analogy to the formation of the Wigner crystal in an infinite electron system.  This change in the system proposed by Wigner \cite{Wigner34} in 1934 forms a monatomic three-dimensional antiferromagnetic crystal.  The analogue of this transition has been reported in many confined electron simulations \cite{CeperleyDM80} \cite{BryantGW87} \cite{CreffieldCE99} \cite{Maksym00}.  This crystallisation we can think of as resulting from the localisation of electrons due to the dominance of the Coulomb interaction over the kinetic interaction of the two particles. 

\begin{figure}
\begin{center}
%\psfrag{nrrs3}{$$n(r).r_s^3$}
\includegraphics[scale=0.40,angle=0]{wig91011.pstex}
\caption{Plot of $n(r).r_s^3$ against radial distance from origin for spheres of radii 9, 10 and 11}
\label{fig:wigone}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.40,angle=0]{wigoneone.pstex}
\caption{Plot of $n(r).r_s^3$ against radial distance from origin for a sphere of radius 10, for the basis {21} and {28}}
\label{fig:wigoneone}
\end{center}
\end{figure}

One consequence of the onset of this crystallisation is that as the sphere increases in size, beyond the value for which the minimum is initially observed, the maximum should move out, away from the origin and towards the wall of the sphere.  In figure \ref{fig:wigtwo} we see the scaled density for spheres of radii 25, 35 and 45, and as predicted we observe this outward migration of the maximum. 

\begin{figure}
\begin{center}
%\psfrag{nrrs3}{$$n(r).r_s^3$}
\includegraphics[scale=0.40,angle=0]{wig253545.pstex}
\caption{Plot of $n(r).r_s^3$ against radial distance from origin for spheres of radii 25, 35 and 45}
\label{fig:wigtwo}
\end{center}
\end{figure}

The densities shown are for the smallest basis sets considered in this research, it hasn't been possible to show densities for larger sized basis sets as a problem was found in the code.  The inclusion of the densities shown already is legitimate given that we expect the system to behave as has been observed.  Also we have seen that the density doesn't change appreciably upon changing the basis set size - such that we would expect the densities of systems described with larger basis sets to have the same form as those already shown. 

\section{The correlation hole}

We saw in the preceding chapter that it is possible to calculate the correlation hole if the density is known, this was shown in \ref{chole}.  The correlation hole is  an interesting quantity as changes in the shape of the correlation hole enable us, at least quantitatively, to understand the spatial correlation between the electrons. We can see this simply by examining the form of equation 2.42.

The term $|\Psi(\mathbf{r_1,r_2})|^2$ is simply the probability of finding one electron at $\mathbf{r_1}$ and the other at $\mathbf{r_2}$.  Division through by the half the density at $\mathbf{r_1}$ gives us a conditional probability of finding an electron at $\mathbf{r_2}$ given that one is at $\mathbf{r_1}$.  If this term is larger than the probability of finding the electron at $\mathbf{r_2}$ (given by $\frac{n(\mathbf{r_2})}{2}$) then $n_c(\mathbf{r_1,r_2})$ will be positive and if not, negative.  Now the probabilities will be determined by the position of the reference electron, if $\mathbf{r_2}$ is chosen close to $\mathbf{r_1}$ then $n_c(\mathbf{r_1,r_2})$ will be negative and is indicative of an unfavourable interaction.  If it is chosen to be a large distance from $\mathbf{r_1}$ then $n_c(\mathbf{r_1,r_2})$ will be positive and we could interpret this qualitatively as an indication of an favourable interaction.

In figure \ref{fig:choleunit} we have a plot of the correlation hole with the reference electron at $r_1=0$, at the origin of the sphere.  We show results for $\theta=0,\pi$ and $\frac{\pi}{2}$.  We see that for all angles the shape of the hole is the same.  If we consider figure \ref{fig:xchole} we see that we would expect no angular dependence for the case $r_1=0$ as motion along a radial vector $r_2$ at an angle $\theta$ made with $r_1$ will be the same irrespective of the value of $\theta$.  In figure \ref{fig:choleunitnzero} the reference electron is moved to a somewhat arbitrary value of $0.3R$ along the $r_1$ mesh, for this system we would expect the shape of the correlation hole to have some angular dependence.  We see that for $\theta=0$ we effectively have $r_2$ superimposed on $r_1$ and we can think of electron two as moving along towards electron one, this is an unlikely occurrence  and hence the correlation hole is strongly negative, with a minimum around fifteen, the position of electron one.  For $\theta=\pi$ we see a large positive correlation hole, this corresponds to motion of electron two which is diametrically opposite to the position of electron one; this situation is extremely favourable.  In the middle ground of these two extremes we have the case of $\theta=\frac{\pi}{2}$, this is not nearly as favourable as for $\theta=\pi$, but preferable to $\theta=0$ 
   
\begin{figure}
\begin{center}
%\psfrag{nrrs3}{$$n(r).r_s^3$}
\includegraphics[scale=0.40,angle=0]{choleunit.pstex}
\caption{Plot of $n_c(\mathbf{r_1,r_2})$ against the radial distance from origin for $\theta=0,\pi,\frac{\pi}{2}$ in the unit sphere for $\mathbf{r_1}=0$}
\label{fig:choleunit}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
%\psfrag{nrrs3}{$$n(r).r_s^3$}
\includegraphics[scale=0.40,angle=0]{choleunitnzero.pstex}
\caption{Plot of $n_c(\mathbf{r_1,r_2})$ against the radial distance from origin for $\theta=0,\pi,\frac{\pi}{2}$ in the unit sphere for $\mathbf{r_1}\ne 0$}
\label{fig:choleunitnzero}
\end{center}
\end{figure}

The form of the correlation hole was also probed as a function of the size of the sphere, figure \ref{fig:cholers} shows $n_c(\mathbf{r_1,r_2})$ for spheres of radii 1, 10, 20, 30 and 40 with the reference electron at $r_1=0$.  We know that for $R\ge10$ we have moved into the Wigner crystal phase, where the electrons are localised in distinct regions of space.  We would expect the correlation hole to be strongly positive around some characteristic value of the size of the Wigner crystal and this is what we observe in our results. For the case of the unit sphere the density is a maximum at the origin and we expect the correlation hole to a have a negative value corresponding with that - as we move away from the origin the density decreases and we see an increase in $n_c(\mathbf{r_1,r_2})$, it becomes slightly positive and then reaches zero at the boundary of the sphere.  For $R>10$ the density at the origin is a minimum and the value of $n_c(\mathbf{r_1,r_2})$ it corresponds to has a local maximum in the direction away from the origin of the sphere.  As the density increases, upon formation of the Wigner crystal we observe a minimum in $n_c(\mathbf{r_1,r_2})$.  As $r_2$ moves to the boundary of the sphere $n_c(\mathbf{r_1,r_2})$ becomes strongly positive and this is a clear indication of finding the other electron diametrically opposite at the boundary of the sphere.  It is interesting to note the enhancement of the correlation hole around the reference electron in the systems shown for $R>10$, so in the Wigner crystal.  We observe that for $R>20$ the correlation hole is further away from the reference electron and we rationalise this as we know that the density is developing a minimum at the centre of the sphere. 

\begin{figure}
\begin{center}
%\psfrag{nrrs3}{$$n(r).r_s^3$}
\includegraphics[scale=0.80,angle=0]{choleradii.pstex}
\caption{Plot of $n_c(\mathbf{r_1,r_2})$ against the radial distance from origin for spheres of radii 1, 10, 20, 30, 40 with $r_1=0$}
\label{fig:cholers}
\end{center}
\end{figure}

\section{The Kohn-Sham potential}

As was mentioned in the theory section knowing the ground state density we can also calculate the Kohn-Sham (KS) potential through use of \ref{vsr}.  In performing these calculations we make use of the fact that the density is solely a function of $r$ such that 

\begin{equation}
\label{simprad}
\nabla^2 n(r)=\frac{2}{r}\frac{\text{d}n(r)}{\text{d}r}+\frac{\text{d}^2 n(r)}{\text{d}{r^2}}.
\end{equation}

We use the method of finite differences to evaluate the derivatives, and the plots of various $v_s(r)$ are shown below in figure \ref{fig:ksp}.  Before discussing these results it is instructive to look at what the KS potential is.  In DFT as proposed by Hohenberg, Kohn and Sham \cite{HohenbergP64} \cite{KohnW65} the KS potential is a one-body, local potential which acts on the non-interacting electrons such that the non-interacting system has the same ground state density as the interacting one.  It is introduced as a mathematical construct with no physical properties.  It has been shown though that this may not be the case and that the KS potential may behave as an actual physical potential \cite{Perdew82} \cite{Levy84}.  These considerations are mentioned for interest alone as this is a topic of some debate.  

\begin{figure}
\begin{center}
\includegraphics[scale=0.5,angle=0]{kspot.pstex}
\caption{Plot of $v_s(r).r_s^2$ against radial distance from origin for spheres of radii 8, 9, 10 and 11}
\label{fig:ksp}
\end{center}
\end{figure}

In figure \ref{fig:ksp} we see a plot of the KS potential, $v_s(r)$, for a range of sphere radii - 8, 9, 10 and 11, where it should be appreciated that we have concentrated our investigations to around the formation of the Wigner crystal.  We see that in the case of the smallest sphere $v_s(r)$ is always negative, at radius 9 around the origin $v_s(r)=0$ and for $R>10$ there is a small, but increasing with radius, positive portion to the potential.  We also see that the tail of the potential becomes more negative.  The development of a positive region of the potential is the development of a repulsive region, whilst the tail becoming more negative at large distances from the origin is indicative of a more attractive region being formed in the potential.  Given that we know that our system, upon formation of the Wigner crystal, has the electrons localised in space this is the form of the potential that we would expect.  For $R>10$ we see that the repulsive region around the origin pushes the electrons away and towards the boundary of the sphere, where the potential becomes more negative - more attractive - upon increasing sphere size. 

\chapter{Conclusion}

This work has been concerned with the solution of two interacting electrons confined to a spherical volume of space by an infinite external potential.  We have calculated the ground state energy for this system in the singlet spin state (antiparallel electron spins) and for a number of different sized spheres.  The convergence of the energy was seen to be two dimensional, in that it depended on two quantum numbers $n$ and $l$.  There appeared to be more dependance upon the radial quantum number $n$ as the system quickly converged upon varying $l$.  We found that the energy varied according to the empirical relationship $E\sim\frac{A}{r_s^2}+\frac{B}{r_s}$; an expression which captures the important physics of the system and which seems to predict different regimes of size dependence.  For instance at low $r_s$, in small spheres, we observe domination of the kinetic energy term, whilst for large $r_s$ we see the dominance of the Coulomb term and subsequent formation of the Wigner crystal.

Whilst a study of this empirical relationship implied the existence of a "Wigner" like transition it isn't until we calculate the ground state density that we observe the formation of such a system.  Here we found that for $1\le R\le 10$ we observed a maximum in the density at the origin of the sphere.  This we can think of as being the un-correlated system, where the kinetic energy dominates.  For $R\ge10$ we find the development of a minimum at the origin of the sphere and the outward progression of a maximum as $R$ increases (figure \ref{fig:wigtwo}).  The development of a minimum for $R=10$ corresponds to an $r_s$ of $\sim8$, this is in agreement with the value of $r_s$ Alavi \cite{AlaviA00} found for the onset of Wigner crystallisation.  Unfortunately there are no other studies of electrons confined in three dimensions as most workers who model these confined systems are trying to explicitly model quantum dots and so are concerned with two dimensional confinement.  

The correlation hole was calculated and gave us an understanding of the correlated motion of the electrons, we saw that a hole was formed around the reference electron such that $n_c(\mathbf{r_1,r_2})\lt0$.  We also understood that regions where $n_c(\mathbf{r_1,r_2})\gt0$ correspond to the probable position of the other electron.  Upon varying $R$ we observe the enhancement of $n_c(\mathbf{r_1,r_2})$ close to the reference electron (corresponding to the development of a local maximum in the origin (figure \ref{fig:cholers})).  In the Wigner crystal the correlation hole is further away from the reference electron due to the fact that $n(\mathbf{r_1})$ develops a minimum at the sphere centre.  Similar results for the density and correlation hole are seen in \cite{AlaviA00}.

Another quantity that was calculated which showed indications of the development of the Wigner crystal was the Kohn-Sham potential.  Here for $R\gt10$ we saw the development of a positive, repulsive, region around the origin of the sphere - the position of the reference electron.  Coupled with this we saw the enhancement of an attractive region near the boundary of the sphere.

This seemingly simple model of two electrons confined to a region of space has shown a wide range of interesting physics and has provided us with an insight into the correlated motion of electrons, something which is of fundamental importance in chemistry.  What now do we see as a continuation of this research?

We can divide our research goals into short term (weeks), mid-term (months) and long term (years) objectives.  In the long term it would be hoped that we could extend the system to the study of more electrons, different symmetries of confinement, the use of a correlated basis set and some more exotic investigations; an example being the application of some time dependent external electric field across the system in order to study time dependent correlation properties.

At present, with the calculational scheme described in this work, these long term goals are simply unpractical.  In order to realise these long term goals we see that work in the mid-term should concentrate on the development of a calculational scheme such that evaluation of the converged ground state energy takes minutes and not hours.  We know that at present most of the time is taken up in evaluating the Coulomb matrix elements so this would seem a logical place to begin enhancing and optimising code.  Included in the appendix of this thesis is a method that has been developed which looks like it may have the potential to drastically speed up the calculation of the Coulomb terms.  The mid term aim will be implementation of this theory and if it is found to be unsuccessful the development of one which does speed up the evaluation of the appropriate terms.

In the short term we aim to find the problem in the code written to evaluate the density, which has prevented us from using densities for any basis greater than $\left\{2,8\right\}$.  The densities shown, as was explained in the discussion of the results, will be similar in form but it will interesting to see how the density changes as the basis set increases.  It was observed in the results presented that whatever the changes will be they will be small.

\chapter{Appendix}

It has been mentioned that one of the key problems with the solution posed above is that the calculations for a basis set of size $\left\{5,8\right\}$ take several hours, this is going to make calculations of systems with more than two electrons unpractical.  A method has been developed such that we can overcome this difficulty, it is a method which has it's basis in the problem as originally solved in \cite{AlaviA00}.  We know that we can write the Coulomb term as a Fourier sum of wave vectors, such that 

\begin{equation}
\frac{1}{\mathbf{|r_1-r_2|}}=\sum_{\mathbf{g}} v_{\mathbf{g}}e^{i\mathbf{g}.\mathbf{\left(r_1-r_2\right)}},
\end{equation}

where

\begin{equation}
\label{intcoeff}
v_{\mathbf{g}}=\frac{1}{\Omega'}\int_{\Omega'} \frac{e^{-i\mathbf{g}.{\mathbf{r}}}}{|\mathbf{r}|} \text{d}\mathbf{r}.
\end{equation}

Here the region of integration $\Omega'$ runs over the maximum possible relative displacements of the two particles.  Thus the Coulomb integral given by \ref{coulint} can be re-written as 

\begin{equation}
\label{sepint}
\text{I}_c=\sum_{\mathbf{g}} \tilde{v_{\mathbf{g}}} \int \eta_{\mathbf{n}}^*(\mathbf{r_1})\eta_{\mathbf{n'}}(\mathbf{r_1})e^{i\mathbf{g}.\mathbf{r_1}}\text{d}\mathbf{r_1} \int \eta_{\mathbf{p}}^*(\mathbf{r_2})\eta_{\mathbf{p'}}(\mathbf{r_2})e^{-i\mathbf{g}.\mathbf{r_2}}\text{d}\mathbf{r_2}.
\end{equation}

where we have adopted a simplification to the notation

\begin{equation}
\mathbf{G}=\mathbf{np}=\left\{n l m p q r\right\}.
\end{equation}

Here lies the power of this method; the integral above for $\text{I}_{c}$ is seperable, and as written each of the one particle integrals can be tabulated in arrays once and for all at the begining of the computation.  The elements of the arrays can be called as and when needed and \ref{sepint} becomes a summation of pre-tabulated numbers, which will be extremely quick to evaluate.

We can write $\text{I}_c$ as

\begin{equation}
\text{I}_c=\sum_{\mathbf{g}} \tilde{v_{\mathbf{g}}}\Lambda_{\mathbf{nn'}}(\mathbf{r_1})\Lambda_{\mathbf{pp'}}(\mathbf{r_2}),
\end{equation}

with

\begin{equation}
\Lambda_{\mathbf{nn'}}(\mathbf{r_1})=\int \eta_{\mathbf{n}}^*(\mathbf{r_1})\eta_{\mathbf{n'}}(\mathbf{r_1})e^{i\mathbf{g}.\mathbf{r_1}}\text{d}\mathbf{r_1}
\end{equation}

and

\begin{equation}
\Lambda_{\mathbf{pp'}}(\mathbf{r_2})=\int \eta_{\mathbf{p}}^*(\mathbf{r_2})\eta_{\mathbf{p'}}(\mathbf{r_2})e^{-i\mathbf{g}.\mathbf{r_2}}\text{d}\mathbf{r_2}.
\end{equation}

Using \ref{basis} we see that we can write the above integrals of $\Lambda$ as

\begin{equation}
\begin{split}
\Lambda_{\mathbf{nn'}}(\mathbf{r_1})
&=N_{nl}N_{n'l'} \int_0^R j_l(\alpha_{nl}r_1)j_{l'}(\alpha_{n'l'}r_1) \times \\
&\left[\int_0^{2\pi} \int_0^\pi Y_{lm}^*(\theta_1,\phi_1)Y_{l'm'}(\theta_1,\phi_1)\sin(\theta_1)e^{igr_1\cos(\theta_1)}\text{d}\theta_1\text{d}\phi{1}\right]r_1^2\text{d}r_1,
\end{split}
\end{equation}

with a corresponding term for $\Lambda_{\mathbf{pp'}}(\mathbf{r_2})$.  We see that we can integrate out the dependence upon $\phi_1$ in the above expression and upon making a change of variable ($x_1=\cos(\theta_1)$) we can simplify the integral for $\theta_1$ such that

\begin{equation}
\begin{split}
\Lambda_{\mathbf{nn'}}(\mathbf{r_1})
&=2\pi N_{nl}N_{n'l'} \int_0^R j_l(\alpha_{nl}r_1)j_{l'}(\alpha_{n'l'}r_1) \times \\
&\left[\delta_{mm'}W_{lm}W_{l'm'}\left\{\int_{-1}^1 P_{lm}(x_1)P_{l'm'}(x_1)e^{igr_1x_1}\text{d}x_1 \right\}\right]r_1^2\text{d}r_1,
\end{split}
\end{equation}

where

\begin{equation}
W_{lm}=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}
\end{equation}

and the $P_{lm}(x_1)$ are the associated Legendre polynomials, a similar analysis can be performed upon $\Lambda_{\mathbf{pp'}}(\mathbf{r_2})$.  The only thing to be calculated now is $\tilde{v_{\mathbf{g}}}$, which is given by \ref{intcoeff}.  We find that

\begin{equation}
\begin{split}
\tilde{v_{\mathbf{g}}}
&=\frac{1}{\Omega'} \int_0^{2\pi} \int_0^\pi \int_0^{2R} \frac{e^{igr\cos{\theta}}}{r}r^2\sin{\theta}\text{d}r\text{d}\theta\text{d}\phi \\
&=\frac{4\pi}{\Omega'g^2}\left[1-\cos{2gR}\right].
\end{split}
\end{equation} 

Thus \ref{sepint} becomes

\begin{equation}
\begin{split}
\text{I}_c
&=\sum_\mathbf{g} \frac{16\pi^3}{\Omega'g^2}\left[1-\cos{2gR}\right] \times \\  &N_{nl}N_{n'l'}N_{pq}N_{p'q'}\int_0^R j_l(\alpha_{nl}r_1)j_{l'}(\alpha_{n'l'}r_1) \times \\
&\left[\delta_{mm'}W_{lm}W_{l'm'}\left\{\int_{-1}^1 P_{lm}(x_1)P_{l'm'}(x_1)e^{igr_1x_1}\text{d}x_1 \right\}\right]r_1^2\text{d}r_1 \times \\
&\int_0^R j_q(\alpha_{pq}r_2)j_{q'}(\alpha_{p'q'}r_2) \times \\
&\left[\delta_{rr'}W_{qr}W_{q'r'}\left\{\int_{-1}^1 P_{qr}(x_2)P_{q'r'}(x_2)e^{-igr_2x_2}\text{d}x_2 \right\}\right]r_2^2\text{d}r_2
\end{split}
\end{equation}

The expression above looks like it may be extremely quick to evaluate; we will be able to calculate the integrals of the Legendre polynomials using a fast Fourier transfrom method, storing the results in an array.  We already have routines to calculate the integral of products of spherical Bessel functions so the overall scheme should be reasonably straightforward to implement and result in a major time saving.  This method, if it is as fast as is hoped, will be extendable to systems of more than two electrons, which is one of the long term aims of the research. 

%
\bibliography{macros,all} 
%
\end{document}

