\documentstyle[a4wide,12pt]{article}

\begin{document}
\section*{Introduction to numerical projects}


Here follows a brief recipe and recommendation on how to write a report for each
project.
\begin{itemize}
\item Give a short description of the nature of the problem and the eventual 
numerical methods you have used.
\item Describe the algorithm you have used and/or developed. Here you may find it convenient
to use pseudocoding. In many cases you can describe the algorithm
in the program itself.

\item Include the source code of your program. Comment your program properly.
\item If possible, try to find analytic solutions, or known limits
in order to test your program when developing the code.
\item Include your results either in figure form or in a table. Remember to
       label your results. All tables and figures should have relevant captions
       and labels on the axes.
\item Try to evaluate the reliabilty and numerical stability/precision
of your results. If possible, include a qualitative and/or quantitative
discussion of the numerical stability, eventual loss of precision etc. 

\item Try to give an interpretation of you results in your answers to 
the problems.
\item Critique: if possible include your comments and reflections about the 
exercise, whether you felt you learnt something, ideas for improvements and 
other thoughts you've made when solving the exercise.
We wish to keep this course at the interactive level and your comments can help
us improve it.
\item Try to establish a practice where you log your work at the 
computerlab. You may find such a logbook very handy at later stages
in your work, especially when you don't properly remember 
what a previous test version 
of your program did. Here you could also record 
the time spent on solving the exercise, various algorithms you may have tested
or other topics which you feel worthy of mentioning.
\end{itemize}



\section*{Format for electronic delivery of report and programs}
%
The preferred format for the report is a PDF file. You can also
use DOC or postscript formats. 
As programming language we prefer that you choose between C/C++, Fortran or Python.
The following prescription should be followed when preparing the report:
\begin{itemize}
\item Use Fronter to hand in your projects, log in  at 
blyant.uio.no and choose 'fellesrom fys3150'.
Thereafter you will see an icon to the left with 'hand in' or 'innlevering'.
Click on that icon and go to the given project. 
There you can load up the files within the deadline.
\item Upload {\bf only} the report file and the source code file(s) you have developed.
The report file should include all of your discussions and a list of the codes you have developed. 
Do not include library files which are available at the course homepage, unless you have
made specific changes to them.
\item Comments  from us on your projects, approval or not, corrections to be made 
etc can be found under
your Fronter domain and are only visible to you and the teachers of the course.

\end{itemize}
Finally, 
we do prefer that you work two and two together. Optimal working groups consist of 
2-3 students. You can then hand in a common report. 





\section*{Project 5, The time-dependent Schr\"odinger equation, deadline   November 30  12am (midnight)}


This project deals with the solution of the time-dependent Schr\"odinger method.
We need here the results from project 2 and you should not do this project if you did not solve the optional part of project 2. You will namely need the eigenvectors you generated in order to solve the time-dependent Schr\"odinger equation.

This project can lead, with extensions, to an eventual publication. However, in its present form it is meant to be finished by the deadline. The project is however demanding. 

We give first a general background to the problem. 

Given a complete set of of orthonormal basis vectors $\left| \phi_k \right\rangle\,$, we expand the state vector $\left|\Psi \right\rangle$ in this basis, and truncate after the first $N$ basis vectors

\begin{equation}
	\left| \Psi \right\rangle = \sum^N_{k=1} c_k \left| \phi_k \right\rangle, \quad c_k = \left\langle \phi_k | \Psi \right\rangle.
\end{equation}

The above coefficients $c_k$ are those which where computed in project 2.
Conversely, any coefficient vector $\mathbf{c}$ defines a unique wave function on this truncated basis. In the case where the basis vectors are labelled by several quantum numbers, we will have to define some mapping of these onto the single index $k$.  The action of operators can now be viewed as matrix multiplication on the coefficient vector. That is, the action of an operator $A\,$ on the state vector corresponds to $\mathbf{A} \mathbf{c}\,$, where the matrix elements of $A$ in the basis of vectors $\left| \phi_k \right\rangle\,$ are given by

\begin{equation}
	\left(\mathbf{A}\right)_{jk} = \left\langle \phi_j  | A | \phi_k \right\rangle.
\end{equation}

In this basis the time-dependent Schr\"odinger equation may be written

\begin{equation}
	i \frac{\partial}{\partial t} \mathbf{c} \left(t\right)= \mathbf{H} \mathbf{c} \left(t\right).
	\label{eq:OnceAgainSchroed}
\end{equation}

This is the general case of the semi-discrete Schr\"odinger equation. 
If we use the eigenvectors $\left|\psi_k\right\rangle$ of the Hamiltonian as a basis, the problem becomes trivial when $H$ is time-independent. In this basis $H$ is diagonal, with the diagonal elements being the eigenvalues of $H$. Equation~(\ref{eq:OnceAgainSchroed}) decouples

\begin{equation}
	i \frac{\partial}{\partial t} \mathbf{c} \left(t\right)= \mathbf{D} \mathbf{c} \left(t\right), \qquad \mathbf{D} 
	= \mathrm{diag} \left( E_1, E_2, \ldots, E_N \right),
\end{equation}
which is what was done in project 2.

The time development is given simply as $c_k\left(t\right)=\exp \left[ -i(t-t_0)E_k \right] c_k \left( t_0 \right)$. We now look at the case where $H$ has a time-dependent perturbation $H_1\left(t\right)$

\begin{equation}
H\left(t\right) = H_0 + H_1 \left( t \right).
\end{equation}

Here, $H_0$ is the stationary part of the Hamiltonian of an, as of yet, unspecified number of particles.  Here we will consider $H_0$ to be the full Hamiltonian of the two-electron problem, including the two-particle Coulomb terms, as it was solved in project 2. 
We will use the eigenvectors $ \left| \psi_k \right\rangle $ of $H_0$ 
as a basis to approximate the total Hamiltonian $H$. 
To find the eigenvectors of the many-body Hamiltonian, an approximation method that gives the ground state and the excited states up to some specified cut-off value must be used. The full diagonalization from project 2 is one possible approach. Using the $N$ vectors, obtained in project 2, 
$\left| \psi_k  \right\rangle$ of the many-body basis, the matrix version of the Hamiltonian becomes

\begin{equation}
 \mathbf{H}=\mathbf{H}_0+\mathbf{H}_1\left(t\right),
\end{equation}
where
\begin{equation}
 \mathbf{H}_0 = \mathrm{diag} \left(E_1,E_2,\ldots,E_N \right), \qquad 
 \left[ \mathbf{H}_1\left( t \right) \right] _{jk} = \left\langle \psi_j | H_1 | \psi_k \right\rangle.
\end{equation}

There is one potential difficulty however: To use the eigenfunctions $\left| \psi_k \right\rangle$ of the stationary part of the Hamiltonian, we have to be able to calculate the matrix elements $\left\langle \psi_j | H_1\left(t\right) | \psi_k \right\rangle$. In the case of the many-particle Hamiltonian this can be prohibitively expensive. In most cases, this certainly rules out calculating these matrix elements for each time step. In many cases, however, the time-dependent perturbation $H_1\left(t\right)$ can be written as a product of a function depending only on time, and a function of spatial coordinates. 
The form of perturbation of interest 
in this project is that of a system which is influenced by an electric field.
Here we will simply limit this perturbation to be time-dependent only for the  two particles
\begin{equation}
H_1 \left(\vec{x},t\right) = E_{0} \sin \left(\omega t\right).
\label{sdjfhkajlsdfhuithrt}
\end{equation}
To solve the time-dependent Schr\"odinger equation we will use a so-called multi-step method developed by Blanes and Moan (Physics Letters A, Vol.~265, pages 35 - 42, (2005)).

We start with our solution to the time-independent Schr\"odinger equation from project 2, rewritten as
\begin{equation}
H_0 \psi _j \left({\bf r}_1,{\bf r}_2\right) = E_j \psi _j \left({\bf r}_1,{\bf r}_2\right).
\end{equation} 
We have limited ourselves to the spatial part of the wave function. Our problem in project 2 does also not have an angular dependence.
Here, $H_0$ is the time-independent part of the Hamiltonian.

The fourth-order time stepping scheme of Blanes and Moan is given for a time
step $n+1$ as
\begin{equation}
\mathbf{\Psi}^{n+1} \approx\left[ \exp{-(\mathbf{M_2})}\exp{\mathbf{M_1}}\exp{\mathbf{M_2}}\right] \mathbf{\Psi}^n,
\label{eq:finals}
\end{equation}
where $\mathbf{M_1}$ and $\mathbf{M_2}$ are defined below
and consist of terms proportional 
to integrals of the form $\int \mathbf{H}\left(s\right)ds$ and  $\int s\mathbf{H}\left(s\right)ds$. 
The error goes like $O(\Delta t^5)$. The variable $n$ represents the number of time steps that have been performed.

When applying a time-dependent perturbation, which in our case is given by Eq.~\ref{sdjfhkajlsdfhuithrt}, 
the Hamiltonian matrix $\mathbf{H}$ is in the eigenvector basis given by
\begin{equation}
\mathbf{H} = \mathbf{H}_0 + \mathbf{H}_1 \left(t\right),
\end{equation}
where $\mathbf{H}_0$ is a diagonal matrix of the energy-eigenvalues, $\mathbf{H}_0 = \mathrm{diag} \left(E_1,E_2,\ldots,E_N \right)$, and $\mathbf{H}_1 \left(t\right)\,$ is given by, in our case, 
\begin{equation}
\left[\mathbf{H}_1 \left(t\right)\right]_{jk} = E \left( t\right) \int^{\infty}_{-\infty}\int^{\infty}_{-\infty} \psi^*_j \left({\bf r}_1,{\bf r}_2\right) \psi_k \left({\bf r}_1,{\bf r}_2\right) d{\bf r}_1 d{\bf r}_2.
\end{equation}

The two terms $M_1$ and $M_2$ are given as
\begin{equation}
	M_1 = -i \int ^{t_n+\Delta t}_{t_n} H \left( s \right) ds , 
\end{equation}
and 
\begin{equation}
	M_2 = \frac{i}{\Delta t} \int^{t_n+\Delta t}_{t_n} \left( s-\left(\frac{\Delta t}{2} -t_n \right)H \left( s \right) ds= \frac{i}{\Delta t} \left[ \int^{t_n+\Delta t}_{t_n} sH\left(s\right) ds \right] +\left(  \frac{\Delta t}{2}-t_n \right)M_1.
\end{equation}

Since our perturbation depends on time only, it is easy to calculate the above integrals. In case we need to compute a more involved expectation value (including for example a spatial dependence), 
then the exponential functions would be matrices and one would need to 
evaluate the exponential of a matrix. 
\begin{enumerate}
\item[a)]  Set up a program which solves Eq.~(\ref{eq:finals}) for the ground state and the first excited state from project 2. Allow for different values 
of the strength in front of the Coulomb interaction from project 2.
Test your program for the ground state of the non-interacting two-particle case with a frequency of the electric field given by the energy 
difference between the ground state and the first excited state. Set
$E_0=1$ and  choose a time interval with
\begin{center}
  \begin{tabular}{ l l}
  	$t_{\mathrm{final}}=600$ & $\Delta t = 0.05$ \\
  \end{tabular}    
\end{center}
A good test of your code is to test the one-particle case with 
no time-dependent perturbation. In this case it is possible to find a simple
closed form expression. Study the ground state and compare with your numerical
results.
\item[b)] It is often difficult to see the physics of a two-variable complex function (three, counting time). It is helpful to introduce the \emph{single-particle density} $\rho_i ({\bf r}_i,t)$ defined by integrating the two-particle density, $\left|\psi \left({\bf r}_1,{\bf r}_2,t\right)\right|\,$, over one of the particles, say particle one
\begin{equation}
\rho_1 \left( {\bf r}_1,t \right) = \int ^{\infty}_{-\infty} \left|\psi \left({\bf r}_1,{\bf r}_2,t\right)\right|^2 d{\bf r}_2.
\end{equation}
In our case we have no angular dependence of the wave function. 
Furthermore, we will only focus on the  density of the relative motion. The center of mass
motion is not of interest since the both the time-dependent part and the interacting part
do not depend affect the center of mass motion.
The density  of interest is therefore 
\begin{equation}
\rho\left( {\bf r},t \right) = \left|\psi \left({\bf r},t\right)\right|^2.
\end{equation}
Since there is no angular dependence (we have chosen $l=0$)
the radius ${\bf r}$ reduces simply to $r$. The angular dependence has been integrated out. 

Compute the above density for both the non-interacting
and the interacting systems.
Study the time-development of the density function for the ground state and
the first excited state using the same strength parameters in front of the
Coulomb interaction as in project 2. Use now the spacing between the ground state and the first excited state of the interacting system to define
the frequency $\omega$. Comment your results.  


\end{enumerate}
There are several tests which have been to be done in addition to those listed above. 
For example one needs to check
that unitarity is observed and that the norm is conserved. However, these topics
are beyond the aim of this project. Furthermore, a realistic perturbation has
a spatial dependence as well. In this case one needs to develop a function
which computes the exponential of a matrix. With these ingredients in place, this method can be extended to consider realistic quantum mechanical cases and could form the background for a possible scientific article.
One very hot topic is to study the temporal evolution of electrons in so-called quantum dots. 
\end{document}