% periodic del2
% eigenstates
% expectation

% * intuition for q and p:  weighted pos'n and change of prob'ty
% * F.T. is sensitive to phase ... it isn't ignorable.  need another example
%   on importance of phase.

% * P < 1:  early on, xref fwd to fock spaces.  Note we *don't* always
%   assume vectors are normalized.

% Notation:  \oplus on elements to replace ordered pairs.  Note that this
%   avoids confusion (if not ambiguity) with inner product.

% * sth about hermite fcns ... even if only as wikipedia ref

% * harmonic oscillator:  classical and QM Hamiltonians?

% xxx does unitary imply self-adjoint?

% xxx I do need uncountable prb spaces ... :/

%% !!! EACH SECTION SHOULD *START* WITH A SUMMARY.  NO MYSTERY NOVEL !!!

%% ================================================================
\section{Introduction}

\begin{center}
\emph{
Everything should be as simple as possible, and no simpler.
} \\
--- Albert Einstein\index{Einstein, Albert} (1879-1955).
\end{center}

The perennial question authors pose in their introductions is:  Why write
another treatment of an old subject?  My reasons for this paper are as follows.

I want to explain basic principles of quantum mechanics to myself in a way that
I can understand and remember --- for general interest as well as for specific
background reference for my research in quantum statistical mechanics.  In a
recent graduate mathematical physics course (just a few months ago!), given my
insufficient level of preparation, I found the level of generality too
overwhelming for myself as a newcomer.  I became lost in all the technicalities
which arise (necessarily!) in the infinite-dimensional case: unbounded
operators, domain restrictions, singular measures, abstract spectral theory,
etc.  Working in the greatest possible generality is the prerogative of
experience, but it is the bane of instruction.

Quantum mechanics is not a single theory, but a \emph{recipe} --- it asserts
some basic postulates and tells you how to write down a PDE describing the
behavior of a system.  The two most common systems to which that recipe is
applied at the introductory level are the infinite well and the harmonic
oscillator.  Yet, these are not the only choices.  In this paper, I choose some
finite-dimensional systems which are easy to explicitly compute with and
visualize.  Infinities notoriously introduce computational challenges;
finite-dimensional systems (especially two-level systems) provide the
opportunity to experiment very concretely and very computably.  These things
are worth working out and worth understanding for their own sakes --- but,
moreover, a solid, computational grounding in the finite-dimensional case makes
the more general, abstract treatments of the subject digestible.

This paper (and this talk) are aimed at mathematics graduate students --- I am
one at present.  However, I've attempted to keep the level of presentation
accessible to an undergraduate who has seen calculus, differential equations,
and linear algebra.  It is for this reason that I downplay measure theory in my
discussion of probability spaces.

I find \cite{bib:griffiths} accessible, and wholeheartedly recommend it as a
prequel to a graduate math-physics course.  Some of my added value in this
paper is that I put as much focus on numerics as on algebraic solutions.  The
presentation in \cite{bib:nc} is similar in spirit to my own, but with
different goals.  See \cite{bib:griffiths}, \cite{bib:mermin}, and
\cite{bib:nc} for more information on quantum mechanics, finite-dimensional and
otherwise.

I acknowledge a large debt of gratitude to Itai Seggev, whose presentation
\cite{bib:seggev} of the solar-neutrino problem (which can be viewed as a
two-level quantum system) inspired many of the computational results of this
paper.  As well, the theoretical statements here are due mainly to Janek Wehr's
mathematical physics course \cite{bib:wehr}.

%% ================================================================
\newpage
\section{The cartoon version}

A friend of mine who is a graduate student in another department --- someone
non-mathematical but well-educated --- recently asked me, ``What does the
\emph{quantum} mean in phrases such as \emph{quantum mechanics}, \emph{quantum
computation}, \emph{quantum information}, etc.?''  I believe that we, as
mathematical experts, should be able to give concise, accurate answers to good
layman questions such as this.  This section is my current answer to this
question; it comes in three parts.

%% ================================================================
\subsection{First part:  granularity}
\label{sec:cartoon_1}

When we look around us, we see things that look continuous, e.g.  the surface
of a pool of water, the top of a table, etc.  But when we look with a
sufficiently powerful magnifying glass, things become discrete.
Electromagnetic energy is radiated from atoms only with certain energy levels;
atoms have most of their mass concentrated in their nuclei with mostly empty
space all around; the periodic table is indexed by integers; there is no such
thing as half an electron.  These little packets of things are called
\emph{quanta}, which is the plural of \emph{quantum}.  This word comes from the
same Latin root as \emph{quantity}.

\begin{center}
\includegraphics{figures/magnifying_glass.eps}
\end{center}

%% ================================================================
\subsection{Second part:  indeterminacy}
\label{sec:cartoon_2}

At the same small distance scales where nature reveals its graininess, we also
see \emph{indeterminacy}.   You can read up (say at Wikipedia) on the
\plainidx{double-slit experiment}.  What this showed is that an electron (this
is true for other particles as well) is really not a localized particle ---
it's sort of a \emphidx{probability cloud} (or \emphidx{wave function} as we
will discuss below) where the electron, before being observed, has various
probabilities of being found in one place or another:

\begin{center}
\includegraphics{figures/cartoon_probty.eps}
\end{center}

When the particle is observed, it is found to be in one specific place.  (This
is called the \plainidx{collapse of the wave function}, which we will discuss
below.)

\begin{center}
\includegraphics{figures/cartoon_p1.eps}
\end{center}

But if we were to go back in a time machine and repeat our observation (or if
we repeat the experiment with many identical electrons), we might find it
another place:

\begin{center}
\includegraphics{figures/cartoon_p2.eps}
\end{center}

(Making a second observation of the \emph{same} particle right after the first
observation has a different effect.  We'll see why in the next section.)

If we keep repeating that experiment and tabulate all the observed positions
and make a \plainidx{histogram} of our results, we will start to see the
histogram taking the form of the probability distribution which was specified
by the wave function:

\begin{center}
\includegraphics{figures/cartoon_pn.eps}
\end{center}

Those $s,p$ orbital figures you saw in high-school chemistry were nothing more
than the 95\%-confidence-interval boundaries of three-dimensional probability
distributions.  (Above, I sketched one-dimensional probability distributions.)

This is the connection between \plainidx{probability and statistics}:
\begin{itemize}

\item Probability\index{probablity} predicts the future; it specifies the likelihoods of various
events occurring.

\item Statistics\index{statistics} describes the past; it tabulates the
proportions of various events which actually occurred.  We use statistical
inference to try to discover what probability model underlied the events
we observed.

\end{itemize}

Nature, observed at the macroscopic scale, is \plainidx{deterministic}: when we
drop a rock, it falls without doubt; planets orbit the sun in very predictable
ways.  In particular, with the advent of Kepler's laws and Newton's mechanics,
the motions of celestial objects could be completely explained without recourse
to the supernatural.  Many people came to believe that we live in a
\emphidx{clockwork universe}:  If one could know exactly the initial conditions
of the universe and the laws which govern time evolution of those conditions,
one could predict the future completely.  This sword is double-edged:  on one
hand, we would have full knowledge of the universe; on the other hand, there
would be no place for free will.  This latter is called the \emph{nightmare of
determinism}.

Quantum indeterminacy was greeted with some skepticism in the early twentieth
century; Albert Einstein famously said, ``God does not play dice with the
universe.''  However, indeterminacy was forced upon us by experiment, and
quantum mechanics is a theory which predicts experimental results quite well.
This sword is also double-edged:  on one hand, nature appears somehow random,
but on the other hand, we are relieved of the nightmare of determinism.

%% ----------------------------------------------------------------
\subsection{Third part:  time evolution}
\label{sec:cartoon_3}

The last piece of quantum theory is that the probability distributions of
as-yet-unobserved particles evolve with time:

\begin{center}
\includegraphics{figures/cartoon_probty_evol.eps}
\end{center}

Once a particle is observed, its position is fixed --- the probability
distribution peaks at the observed location when the wave function collapes ---
but this distribution continues to evolve with time, typically
\plainidx{spreading}:

\begin{center}
\includegraphics{figures/cartoon_collapse.eps}
\end{center}

The law governing the time evolution of wave functions is called
\emphidx{Schr\"odinger's equation}.

That's the qualitative cartoon version.  My goal in the remainder of this paper
is to show how these notions are made quantitative.

The last question one asks is, \emph{Who cares?}  Why study this stuff at all?
The best and shortest answer is simple:  Without quantum mechanics there would
be no Information Age, period.

%% ================================================================
\newpage
\section{State spaces}

In this section we will see how probability spaces (appendix \ref{appx:prb})
and complex vector spaces (appendix \ref{appx:la}) are combined, in the
language of quantum mechanics, into the notion of state spaces.

%% ----------------------------------------------------------------
\subsection{The first postulate:  wave functions}

Consider a single toss of a fair coin.  The coin may land heads-up or tails-up
on any given toss; we don't know until we carry out the experiment.  While the
coin is spinning in the air, we say it's neither heads-up nor tails-up, but the
fairness of the coin is flying along with it.  Only when the coin hits the
table do we say that the result of the flip has been determined.

The coin, once it lands, can be in one of two positions, namely,  heads or
tails, and each occurs half the time --- that is, with \plainidx{probability}
1/2.
We can write this graphically, as in the cartoon section.  However, instead of
$x$ with all its range of values, there are only two choices on the horizontal
axis:  $H$ and $T$.

We can also tabulate these probabilities using an ordered pair of numbers:
$$
	\colvectwo{0.5}{0.5}.
$$
(Using the terminology of appendix \ref{appx:prb}, this is a
\emphidx{probability mass function} for the probability space whose two
elements are heads and tails.)

Now suppose the coin isn't fair:  it lands heads-up with probability 0.64 and
tails-up with probability 0.36.  We can tabulate the event probabilities as
another ordered pair:
$$
	\colvectwo{0.64}{0.36}.
$$
In general, we have
$$
	\colcvectwo{p}{1-p}.
$$

Physical particles (for example, electrons) likewise can be observed to be in a
variety of physical states:
\begin{itemize}
\item maybe just one,
\item maybe two (e.g. electron spin up or down),
\item countably many (e.g. bound-electron energy levels), or
\item uncountably infinitely many (e.g. free-electron energy levels).
\end{itemize}
These situations, however many there may be, are measured with various
probabilities.  One of the bizarre things about quantum mechanics is that,
until the measurement, the particle is like the coin still flipping in mid-air:
all we can say about the particle's state is to tabulate the various
probabilities of what the state \emph{will} be when it is measured.

That much is perhaps not surprising, but one of the even more bizarre things
about quantum mechanics is that the probabilities are \emph{factored} into
complex-conjugate pairs.  (This fact is not obvious and it was not easily
discovered; it is justified only by history and laboratory experiment.  See
\cite{bib:griffiths} for some history.)

For the coin, we can write down a two-state \bfidx{wave function}
$$
	\psi =
	\colvectwo{a}{b}
$$
where $a$ and $b$ are complex numbers.  That is,
$$
	\psi \in \C^2.
$$
The probability of the coin landing heads-up (the first state) is
$$
	P(H) = |a|^2 = a^* a.
$$
Likewise, the probability of the coin landing tails-up (the second state) is
$$
	P(T) = |b|^2 = b^* b.
$$
Since it is certain that the coin lands either one way or the other
(i.e. with probability 1), we insist on the \bfidx{normalization}
$$
	P(H) + P(T) = |a|^2 + |b|^2 = 1.
$$
Then we have
$$
	|a|^2 = p, \quad |b|^2 = 1-p.
$$
For example,
$$
	\colvectwo{a}{b} =
	\colvectwo{\sqrt{0.5}}{\sqrt{0.5}}
$$
is a wave function for the fair coin.  Note however that
$$
	\colvectwo{a}{b} =
	\colvectwo{\sqrt{0.5}i}{\sqrt{0.5}i}
$$
also works since
$$
	(-\sqrt{0.5}i \cdot \sqrt{0.5}i) +
	(-\sqrt{0.5}i \cdot \sqrt{0.5}i) = 0.5+0.5 = 1.
$$
Likewise,
$$
	\colvectwo{a}{b} =
	\colvectwo{\sqrt{0.5}e^{i\alpha}}{\sqrt{0.5}e^{i\beta}}
$$
works for real $\alpha,\beta$.  These $e^{i\alpha}$ and $e^{i\beta}$ are called
\bfidx{phase factors}\index{complex phase}.

For the unfair coin, we could have wave functions
$$
	\colcvectwo{a}{b} =
	\colvectwo{0.8}{0.6}, \quad
	\colvectwo{0.8i}{0.6i}, \quad
	\colvectwo{0.8e^{i\alpha}}{0.6e^{i\beta}},
$$
etc.  In general,
$$
	\colcvectwo{a}{b} =
	\colcvectwo{\sqrt{p}e^{i\alpha}}{\sqrt{1-p}e^{i\beta}},
$$
is a wave function for the $p$-weighted coin.

As we will see in section \ref{sec:time_evolution}, the $a$ and $b$ may change
with time.  If the complex magnitudes change, then the probabilities change,
but if only the complex phases change with time, then the probabilities $|a|^2$
and $|b|^2$ do not change with time.

\begin{rem}
\label{rem:phase_not_obs}
That is to say: \emph{The complex phase factors disappear in the probability
computations; they have no observable significance.}
\end{rem}

One might ask, why do we need the complex part?  It turns out that it is the
complex wave function which participates in the differential equation (the
\plainidx{Schr\"odinger equation}, section \ref{sec:SE}) which
governs the time development of wave functions.

% xxx pure state of coin -- not flipping but parachuting?

Extensions of this concept:

\begin{itemize}

\item If we had not a coin but a six-sided die, then the wave function would be
not in $\C^2$ but $\C^6$.  The six numbers $|c_1|^2$ through $|c_6|^2$
represent the probability of the die landing in states $\phi_1$ through
$\phi_6$, namely, rolls of 1 through 6 respectively.

\item If we were considering countably infinitely many energy levels (e.g. an
elecrtron bound to an atom), then the wave function would be
infinite-dimensional, and the elements $|c_i|^2$ for $i=0,1,2,3,\ldots$ would
form a discretely probability density function.  Here
$$
	P(\psi=\phi_j) = |c_j|^2,
$$
i.e. the numbers $|c_i|^2$ tabulate the probability of $\psi$ being found in
the $j$th state, $\phi_j$, upon measurement.  The normalization condition
forces $\psi$ to be a \bfidx{square-summable
sequence}, i.e. $\psi \in \ell^2(\C)$:
$$
	\sum_{j=0}^\infty P(\psi = \phi_j) =
	\sum_{j=0}^\infty |c_j|^2 = 1.
$$

\item If we were considering uncountably infinitely many energy levels (e.g. a
free electron) then we would need to use the theory of continuous probability
distributions.  (See \cite{bib:kerl} for more information on this topic.)  The
normalization criterion forces $\psi$ to be \bfidx{square integrable}, i.e.
$\psi \in L^2(\R^d,\C)$:
$$
	\int_{\R^d} |\psi(\vecx)|^2 d\vecx = 1.
$$
In fact, this is the historical origin of the study of $L^2$ spaces.

\end{itemize}

These examples motivate the following definition.

\begin{defn}
For a finite probability space $\Omega$, with $\#\Omega = N$, the \bfidx{state
space} $V$ is $\C^N$.  Each element $\psi$ of $V$ with $\|\psi\|=1$ is called a
\bfidx{state vector}, and represents a probability distribution (definition
\ref{defn:pmf}) on $\Omega$.
\end{defn}

I summarize for emphasis: \textbf{Each state vector encodes a probability
distribution.}

The \bfidx{first postulate of quantum mechanics}\index{postulate, first} is
(following \cite{bib:nc}) simply the following: to any isolated physical system
is associated a state space, where the state of the system is described by a
unit vector in that space.

The \bfidx{fourth postulate of quantum mechanics}\index{postulate, fourth} ---
a topic which is not pursued here --- is that the state of a multi-component
system is represented by vectors in the tensor product of the state spaces for
each component system.

%% ----------------------------------------------------------------
\subsection{Expectations}

If we denote the heads-up state as $\phi_0$ and the tails-up state as $\phi_1$,
then we already saw that the wave function
$$
	\psi
	= \colcvectwo{a}{b}
	= \colcvectwo{\sqrt{p}e^{i\alpha}}{\sqrt{1-p}e^{i\beta}}
$$
encodes the probability of the coin landing heads-up or tails-up.

We can write down the \bfidx{mean} (or \bfidx{expectation} or \bfidx{expected
value}) of various \emph{functions} of the states, much as in section
\ref{sec:mean}.  Suppose you and I are flipping a $p$-weighted coin with the
agreement that if the coin lands heads-up, I win \$2 but if it lands tails-up I
lose \$3.  Writing 0 for heads and 1 for tails, my winnings are
$$
	w(0) = 2, \quad w(1) = -3.
$$
Then my expected earnings over the long run are
$$
	\expn{w(j)}
	= \sum_j w_j P(\psi=\phi_j)
	= w(0) p + w(1) (1-p)
	= 2p - 3(1-p)
	= 5p-3.
$$
(For a fair coin, I'm losing 50 cents per flip on average.  As discussed in the
appendix, I'd be winning money in the long run as long as this expectation is
positive --- i.e. if $p>0.6$.)

Now remember that wave functions factor the probabilities into
complex-conjugate pairs.  If the result of a coin-flipping experiment depends
on the heads-or-tails result, i.e. $f(\psi)$ where $\psi$ is the ordered pair
$$
	\psi = \colcvectwo{\sqrt{p}e^{i\alpha}}{\sqrt{1-p}e^{i\beta}},
	= \colcvectwo{a}{b},
$$
then equation \ref{eqn:expn_fj} becomes
\begin{align*}
	\expn{f(\psi)} = \sum_j f(\phi_j) P(X=j)
	&= f(\phi_0) |a|^2 + f(\phi_1) |b|^2 \\
	&= f(\phi_0) a^*a + f(\phi_1) b^*b \\
	&= a^* f(\phi_0) a + b^* f(\phi_1) b.
\end{align*}
In general, we put the complex-conjugate part of the probability factor
to the left of $f$ and the non-conjugated part to the right:
$$
	\expn{f(\psi)} = \sum_j c_j^* f(\phi_j) c_j.
$$

%% ----------------------------------------------------------------
\subsection{The second postulate:  observables and measurements}
\label{sec:second_postulate}

The \bfidx{second postulate of quantum mechanics}\index{postulate, second},
which again is justified by history and experiment, is as follows:

\begin{itemize}

\item An observation corresponds to a Hermitian operator $A:V\to V$ on the
state space $V$.

\item Recall from section \ref{sec:spectral_theorem} that, since $A$ is
Hermitian, $V$ has an orthonormal basis $\{\phi_1, \ldots, \phi_N\}$ of
eigenvectors of $A$, with respective eigenvalues $\{a_1, \ldots, a_N\}$.  If
$\psi$ is the wave function of our system, then it is a linear combination of
the eigenbasis:
$$
	\psi = \sum_{j=1}^N c_j \phi_j
$$
where
$$
	c_j = \braket{\phi_j}{\psi}.
$$

\item Measurement of the quantity to which $A$ corresponds has as its possible
outcomes the numbers $a_j$ with probabilities $|c_j|^2$.

\end{itemize}

If this is to work, we would want the probabilities to be non-negative (which
we have by positive-definiteness of the inner product), and we would want them
to sum to 1.  In fact they do, as shown in the following proposition.

\begin{prop}
We have
$$
	\sum_j |\braket{\phi_j}{\psi}|^2 = 1.
$$
\begin{proof}
First, we know we've normalized $\psi$ so that
$$
	\|\psi\|^2 = \braket{\psi}{\psi} = 1.
$$
Since the eigenfunctions $\phi_j$ are an orthonormal spanning set, we have
$$
	\psi=\sum_j c_j \phi_j
$$
with
$$
	c_j = \braket{\phi_j}{\psi}.
$$
Then
\begin{align*}
	1 &= \braket{\psi}{\psi}
		= \left\langle
			\sum_i c_i\,\phi_i \,\Big{|}\, \sum_j c_j \,\phi_j
		\right\rangle
		= \sum_i \sum_j \braket{c_i \phi_i}{ c_j \phi_j } \\
	&= \sum_i \sum_j c_i^* c_j \delta_{ij} \\
	&= \sum_j |c_j|^2
\end{align*}
by the orthonormality of the $\phi_j$'s.  But the $c_j$'s were just shorthand
for $\braket{\phi_j}{\psi}$ so we have
$$
	\sum_j |\braket{\phi_j}{\psi}|^2 = 1
$$
which is what we wanted to show.
\end{proof}
\end{prop}

%% ----------------------------------------------------------------
\subsection{Expectations, continued}
\label{sec:expn_braCket}

The expected value of $A$ is
$$
	\expn{A} = \sum a_j P(a_j).
$$
From the second postulate, those probabilities $P(a_j)$ are
\begin{align*}
	P(a_j) &= |\braket{\phi_j}{\psi}|^2 \\
	&= \braket{\psi}{\phi_j} \braket{\phi_j}{\psi}.
\end{align*}
Regrouping (since matrix multiplication is associative), we have
\begin{align*}
	P(a_j) &= \bra{\psi}\; \bigg(\ket{\phi_j}\bra{\phi_j}\bigg) \; \ket{\psi}.
\end{align*}
Then
$$
	\expn{A} = \sum \bra{\psi}\;
		\bigg(a_j \ket{\phi_j}\bra{\phi_j}\bigg)
		\; \ket{\psi}.
$$
Given the \plainidx{spectral decomposition} of $A$ (see section
\ref{sec:spectral_theorem}) as
$$
	A = \sum a_j \ket{\phi_j}\bra{\phi_j},
$$
we have
$$
	\expn{A} = \sum \braCket{\psi}{A}{\psi}.
$$

This\index{expected value} is the \plainidx{average value} or \plainidx{mean}
of all the observations, averaged over repeated measurements on identical
systems.

%% ================================================================
\newpage
\section{Time evolution}
\label{sec:time_evolution}

%% ----------------------------------------------------------------
\subsection{The Schr\"odinger equation}
\label{sec:SE}

As verified by experiment, wave functions evolve according to the
\bfidx{Schr\"odinger equation}:
$$
	i \hbar \frac{\D\psi(t,\vecx)}{\Dt}
	= \frac{-\hbar^2}{2m} \nabla^2 \psi(t,\vecx) + V(t,\vecx) \psi(t,\vecx)
	= H \psi(t,\vecx)
$$
i.e.
$$
	\frac{\D\psi}{\Dt}
	= \frac{-i H}{\hbar} \psi
$$
where $H$ is the Hermitian operator
$$
	H = \frac{-\hbar^2}{2m} \nabla^2 + V.
$$
The wave function varies in position and time, as in the cartoon version; the
\bfidx{potential function}, conventionally written with the letter $V$, may
also vary in position and time.  The potential represents external interactions
with the particle.  The quantity $\hbar$ is best thought of as being very
small.

We can write down a solution
\begin{align}
\label{eqn:SE}
	\frac{d\psi}{dt} &= \frac{-i}{\hbar} \, H \psi \\
	\psi(t) &= e^{-iHt/\hbar} \psi(0).
\end{align}
One checks this solution by differentiating:
\begin{align*}
	\frac{d\psi}{dt} &= \frac{-iH}{\hbar} \, e^{-iHt/\hbar} \psi(0)
	= \frac{-iH}{\hbar} \, \psi(t)
	= \frac{-i}{\hbar} \, H \psi(t).
\end{align*}

Recall that the Hamiltonian is hermitian.  As shown in proposition
\ref{prop:exp_hmn_unitary}, the imaginary exponential of a hermitian matrix is
\bfidx{unitary}.  The importance of this is that unitary matrices $U$ are
\bfidx{norm-preserving}.  That is, if $\|\psi\|=1$ then $\|U\psi\|$ is still 1;
probability functions remain probability functions as they evolve in time.

In fact, \cite{bib:wehr} makes the following grand summary:

\begin{center}
\emph{Quantum mechanics says that the energy operator generates the group of
time-evolution operators.}
\end{center}

The group structure is addressed in section \ref{sec:neutrino_numerics}.

%xxx compose $U(s)$ and $U(t)$.

%Mean/expectation.

%Real spectrum.  Existence of orthonormal basis of eigenvectors.
%Why they \emph{must} be orthonormal.

%% ----------------------------------------------------------------
\subsection{The third postulate:  wave-function collapse and subsequent
evolution}
\label{sec:third_postulate}

The \bfidx{third postulate of quantum mechanics}\index{postulate, third} is
that immediately after measurement of an observable $A$, the state of the
system is described \emph{only} by an eigenstate $\phi_k$ of $A$ (using the
notation of section \ref{sec:second_postulate}).  This is called the
\emphidx{collapse of the wave function}.  The evolution of the single-particle
system thereafter is described by the Schr\"odinger equation, with new initial
conditions.

\beginex
Let $A$ be an observable with orthonormal eigenbasis $\{\phi_1, \ldots,
\phi_N\}$.  We can write our wave function $\psi$ in terms of this basis as
$$
	\psi = \sum_{j=1}^N c_j \phi_j.
$$

If $N=2$ then maybe $c_1=0.6$ and $c_2=0.8$.  Then $\psi$ in this
basis is
$$
	\psi=\colvectwo{0.6}{0.8}.
$$
The probability of the particle being observed in the first or second state is
0.36 or 0.64, respectively.  Suppose the former.  Then $\psi$ is now
$$
	\psi=\colvectwo{1}{0}.
$$
After that, as the wave function evolves, probability may leak back into
the second slot --- or not, as shown in the next section.
\mendex

%% ----------------------------------------------------------------
\subsection{Eigenstates of the Hamiltonian}

Suppose a state $\psi_k$ is an eigenstate of the Hamiltonian (not just an
observable $A$, but the same Hamiltonian $H$ which participates in the
Schr\"odiner equation) with eigenvalue $E_k$.  Then
$$
	H \psi_k = E_k \psi_k.
$$
Then equation \ref{eqn:SE} becomes
\begin{align*}
	\psi(t) &= e^{-i E_k t/\hbar} \psi(0).
\end{align*}
Notice that the quantity $e^{-i E_k t/\hbar}$ is a complex scalar with
magnitude 1.  In particular, with a two-dimensional state space,
$$
	\psi(t) = \colvectwo{a(t)}{b(t)}
$$
and so
$$
	\colvectwo{a(t)}{b(t)}
	= e^{-i E_k t/\hbar} \colvectwo{a(0)}{b(0)}
	= \colvectwo
		{e^{-i E_k t/\hbar} a(0)}
		{e^{-i E_k t/\hbar} b(0)}.
$$
That is, the coefficients $a(t)$ and $b(t)$ evolve in time with only their
complex phases being modified; their complex magnitudes do not change in time.

%% ----------------------------------------------------------------
\subsection{Inhibition of time evolution}

There are two ways to make a system which doesn't evolve in time:

\begin{itemize}
\item Construct a system with a zero Hamiltonian, or at least one with a small
	enough Hamiltonian such that $e^{-iHt/\hbar}$ is sufficiently close to
	identity for the duration of the experiment, or
\item Construct a system in which the initial state is an eigenstate of the
	Hamiltonian.
\end{itemize}

%% ================================================================
\newpage
\section{Numerical results with a two-state system}

As described in \cite{bib:seggev}, a particular two-state quantum system may be
described wherein a neutrino oscillates between one of two flavors.  The
parameters (which I simplify from Seggev's more physically realistic
description) are \plainidx{masses} $a$ and $b$ and \plainidx{mixing angle}
$\theta$.  The interpretation of the mixing angle is that if it is zero, no
oscillation between neutrino flavors happens; if it is $\pi/4$, then full
oscillation occurs.

In this system, the wave function involves $t$ only; there is no notion of
position dependence.

The Hamiltonian is (given here without explanation):
$$
	H = \pcmatrix{
		a^2 \cos^2\theta + b^2 \sin^2\theta &
		(a^2-b^2)\cos\theta\sin\theta \\
		(a^2-b^2)\cos\theta\sin\theta &
		a^2 \sin^2\theta + b^2 \cos^2\theta \\
	}.
$$
Note that this is real and symmetric, hence Hermitian.

%% ----------------------------------------------------------------
\subsection{Algebra}

Compute the determinant of $H-\lambda I$.  After some algebra, we get
eigenvalues
$$
	\lambda = a^2, b^2
$$
with respective eigenvectors
$$
	\psi=
	\colrvectwo{\cos\theta}{\sin\theta},
	\colrvectwo{\sin\theta}{-\cos\theta}.
$$

%% ----------------------------------------------------------------
\subsection{Numerics}
\label{sec:neutrino_numerics}

Time evolution (with $\hbar=1$):
\begin{align*}
	U(t) &= e^{-i\,H\,t}.
\end{align*}
Composition:
\begin{align*}
	U(s) &= e^{-i\,H\,s} \\
	U(s+t) &= e^{-i\,H\,s} e^{-i\,H\,t} \\
		&= e^{-i\,H\,(s+t)}.
\end{align*}
Iterative approximation of $\psi(t)$:
\begin{align*}
	\psi(t) &= e^{-i\,H\,t}\, \psi(0) \\
	\psi(t+\Delta t)
		&= e^{-i\,H\,(t+\Delta t)}\, \psi(0) \\
		&= e^{-i\,H\,\Delta t}\, e^{-i\,H\,t}\, \psi(0) \\
		&= e^{-i\,H\,\Delta t}\, \psi(t).
\end{align*}

Matlab code:

\begin{verbatim}
a = 0.4; b = 0.5; theta = 0.6 * pi/4;
dt = 0.1; niter = 1000;
H = [ a^2*cos(theta)^2 + b^2*sin(theta)^2, (a^2-b^2)*cos(theta)*sin(theta);
    (a^2-b^2)*cos(theta)*sin(theta),      a^2*sin(theta)^2 + b^2*cos(theta)^2 ];
Udt = expm(-i * H * dt);
% Other print statements here ...
psit = [1.0; 0.0];
measurement_iter = round(4*niter/5);
for iter = 1:niter
    t     = iter*dt;
    psit  = Udt * psit;

    if (iter == measurement_iter) % Measurement
        u = rand(1)
        if (u < prob(1))
            psit = [1;0];
        else
            psit = [0;1];
        end
    end

    % Other computations / plots ... .
end
\end{verbatim}

In black is the probability $|\psi_0|^2$ for the neutrino to be found in state
0, and in blue is the probability $|\psi_1|^2$ for it to be found in state 1.
Note that the probabilities always sum to 1 as they vary.  After measurement,
the wave function collapses onto one of the eigenstates of the observable (not
an eigenstate of the Hamiltonian) --- i.e. state 0 or state 1 --- and then
continues to evolve in time.

\begin{center}
\includegraphics{code/neutrinos.eps}
\end{center}

Program output (see \texttt{code/neutrinos.m}):

\begin{verbatim}
H:
    0.1686   -0.0265
   -0.0265    0.2414

U(dt):
   0.9999 - 0.0169i   0.0001 + 0.0026i
   0.0001 + 0.0026i   0.9997 - 0.0241i

Eigenvalues of H:
    0.1600
    0.2500

Eigenvectors of H (columns):
   -0.9511   -0.3090
   -0.3090    0.9511

psi(0):
     1
     0

Is H Hermitian?  Here are H* and H:
    0.1686   -0.0265
   -0.0265    0.2414

    0.1686   -0.0265
   -0.0265    0.2414

Is U(dt) unitary?  Here are U(dt)*, U(dt), and U(dt)* U(dt):
   0.9999 + 0.0169i   0.0001 - 0.0026i
   0.0001 - 0.0026i   0.9997 + 0.0241i

   0.9999 - 0.0169i   0.0001 + 0.0026i
   0.0001 + 0.0026i   0.9997 - 0.0241i

   1.0000            -0.0000 + 0.0000i
  -0.0000 - 0.0000i   1.0000          

\end{verbatim}

%% ================================================================
\newpage
\section{Numerical results with a discrete 2D system}

Here I want to quantify the qualitative picture given in the cartoon version
(section \ref{sec:cartoon_1}).  Namely, let's write down a position-dependent
wave function and see how it evolves in time.  We can do that numerically (as
always in numerics) by sampling the $x$ axis (or $x,y$ axes, etc.) along a
\plainidx{discrete mesh} --- just like a graphing calculator does:

\begin{center}
\includegraphics{figures/wave_fcn_sampled.eps}
\end{center}

Above, we found a closed-form solution for the Schr\"odinger equation:
given
$$
	i \hbar \frac{\D\psi}{\Dt}
	= \frac{-\hbar^2}{2m} \nabla^2 \psi + V \psi
	= H \psi,
$$
we had
$$
	\psi(t) = e^{-iHt/\hbar} \psi(0).
$$
Here we instead solve the Schr\"odinger equation directly, using approximation
methods.

%% ----------------------------------------------------------------
\subsection{The time derivative}

Recall from freshman calculus that
$$
	\frac{\D\psi}{\Dt} \approx
		\frac{
			\psi(t+\Delta t, \vecx) -
			\psi(t,          \vecx)
		}{\Delta t}.
$$
Given a sequence
$$
	t_k = t_0 + k\Delta t,
$$
we have
$$
	\frac{\D\psi}{\Dt} \approx
		\frac{
			\psi(t_{k+1}, \vecx) - \psi(t_k, \vecx)
		}{\Delta t}.
$$

The \bfidx{Euler method} for solving a single-variable ODE is to use
this first-order approximation.  Namely, for
$$
	\frac{dy}{dt} = f(t,y),
$$
we use the approximation
$$
	\frac{y(t_{k+1}) - y(t_k)}{\Delta t} \approx f(t_k, y(t_k))
$$
and solve for $y(t_{k+1})$ to obtain
$$
	y(t_{k+1}) \approx y(t_k) + f(t_k, y(t_k)) \Delta t.
$$
Here, we start with
\begin{align*}
	i \hbar \frac{\D\psi}{\Dt}
		&= \frac{-\hbar^2}{2m} \nabla^2 \psi + V \psi \\
	\frac{\D\psi}{\Dt}
		&= i \left(
			\frac{\hbar}{2m} \nabla^2 \psi - \frac{V}{\hbar} \psi
		\right) \\
	\frac{\D\psi}{\Dt}
		&= i\left(c \nabla^2 \psi - V \psi\right)
\end{align*}
where I put
$$
	c = \frac{\hbar}{2m}
$$
and (with abusive notation) set $\tilde V = V/\hbar$, then replace
$\tilde V$ with $V$.  Then the Euler method gives
$$
	\psi(t_{k+1},\vecx) \approx \psi(t_k, \vecx) + i (c \nabla^2\psi - V \psi)
		\Delta t.
$$
The only question is how to compute the Laplacian on the right-hand side,
which I'll do below.

You can read up on the fourth-order \bfidx{Runge-Kutta method} in any
textbook on numerical methods.  In short, given
$$
	\frac{dy}{dt} = f(t,y),
$$
the approximation is
$$
	y(t_{n+1}) = y(t_n) + h/6(k_1 + 2 k_2 + 2 k_3 + k_4)
$$
where
\begin{align*}
	k_1 &= f(t_n,       y_n) \\
	k_2 &= f(t_n + h/2, y_n + k_1 h/2) \\
	k_3 &= f(t_n + h/2, y_n + k_2 h/2) \\
	k_4 &= f(t_n + h,   y_n + k_3 h).
\end{align*}
If this is confusing, just think of it as accomplishing what the Euler
method does, but with less error.

%% ----------------------------------------------------------------
\subsection{The one-dimensional Laplacian}

Also from freshman calculus, the one-dimensional Laplacian is approximately
$$
	\frac{\D^2\psi}{\Dx^2} \approx
		\frac{
			\psi(t, x + \Delta x) -
			2 \psi(t, x) +
			\psi(t, x - \Delta x)
		}{(\Delta x)^2}.
$$
With a mesh of $M$ points $x_0, \ldots, x_{M-1}$, we have
$$
	\frac{\D^2\psi}{\Dx^2}\pickpipe{k} \approx
		\frac{
			\psi(t, x_{k+1}) -
			2 \psi(t, x_k) +
			\psi(t, x_{k-1})
		}{(\Delta x)^2}.
$$

Computing the discrete Laplacian at $x_k$ involves sampling $\psi$
there and at two neighbor points:

	\setlength{\unitlength}{1.0in}
	\begin{picture}(6.5,0.5)
	%\put(0,0){\framebox(6.5,0.5)}

	\put(0.60,0.10){\line(1,0){5.0}}

	\put(1.60,0.10){\circle*{.05}}
	\put(2.10,0.10){\circle*{.05}}
	\put(2.60,0.10){\circle*{.05}}

	\put(1.00,0.30){\ldots}
	\put(1.50,0.30){$x_{k-1}$}
	\put(2.00,0.30){$x_k$}
	\put(2.50,0.30){$x_{k+1}$}
	\put(3.00,0.30){\ldots}
	\end{picture}

The question is, what to do for $\psi(t, x_{k-1})$ at the left edge $k=0$,
and for $\psi(t, x_{k+1})$ at the right edge $k=M-1$.

	\setlength{\unitlength}{1.0in}
	\begin{picture}(6.5,0.5)
	%\put(0,0){\framebox(6.5,0.5)}

	\put(0.60,0.10){\line(1,0){5.0}}

	\put(0.60,0.10){\circle{.05}}
	\put(1.10,0.10){\circle*{.05}}
	\put(1.60,0.10){\circle*{.05}}

	\put(0.50,0.30){$???$}
	\put(1.00,0.30){$x_0$}
	\put(1.50,0.30){$x_1$}

	\put(4.60,0.10){\circle*{.05}}
	\put(5.10,0.10){\circle*{.05}}
	\put(5.60,0.10){\circle{.05}}

	\put(2.50,0.30){\ldots}

	\put(4.50,0.30){$x_{M-2}$}
	\put(5.00,0.30){$x_{M-1}$}
	\put(5.50,0.30){$???$}

	\end{picture}

This dilemma is resolved by a choice of \plainidx{boundary conditions}:
\begin{itemize}
\item Dirichlet-zero\index{Dirichlet boundary conditions}, or
	\emphidx{clamped}, boundary conditions, taking
	$$
		\psi(t,x_{-1}) = \psi(t,x_{M})= 0.
	$$
\item Periodic\index{periodic boundary conditions} boundary conditions, taking
	$$
		\psi(t,x_{-1}) = \psi(t,x_{M-1})
		\qqand
		\psi(t,x_{M}) = \psi(t,x_{0}).
	$$

\end{itemize}

Note that in the 1D case, it is particularly easy to view the discrete
Laplacian (and hence the Hamiltonian) as a real symmetric matrix.
With Dirichlet-zero boundary conditions we have
$$
	\Delta \psi = \prmatrix{
		-2 &  1 &     {} &     {} &     {} & {} & {} \\
		1 & -2 &     -1 &     {} &     {} & {} & {} \\
		{} &  1 &     -2 &     -1 &     {} & {} & {} \\
		{} & {} & \ddots & \ddots & \ddots & {} & {} \\
		{} & {} &     {} &      1 &     -2 &  1 & {} \\
		{} & {} &     {} &     {} &      1 & -2 &  1 \\
		{} & {} &     {} &     {} &     {} &  1 & -2 \\
	}
	\pcmatrix{
		\psi(t, x_0) \\
		\psi(t, x_1) \\
		\psi(t, x_2) \\
		\vdots \\
		\psi(t, x_{M-3}) \\
		\psi(t, x_{M-2}) \\
		\psi(t, x_{M-1}) \\
	}
$$
With periodic boundary conditions, this is
$$
	\Delta \psi = \prmatrix{
		-2 &  1 &     {} &     {} &     {} & {} &  1 \\
		1 & -2 &     -1 &     {} &     {} & {} & {} \\
		{} &  1 &     -2 &     -1 &     {} & {} & {} \\
		{} & {} & \ddots & \ddots & \ddots & {} & {} \\
		{} & {} &     {} &      1 &     -2 &  1 &  {} \\
		{} & {} &     {} &     {} &      1 & -2 &  1 \\
		1  & {} &     {} &     {} &     {} &  1 & -2 \\
	}
	\pcmatrix{
		\psi(t, x_0) \\
		\psi(t, x_1) \\
		\psi(t, x_2) \\
		\vdots \\
		\psi(t, x_{M-3}) \\
		\psi(t, x_{M-2}) \\
		\psi(t, x_{M-1}) \\
	}
$$
The Hamiltonian is written by adding $V(t,x_i)$ along the main diagonal.

Using Matlab, you can pass this matrix to the \texttt{eig} routine to find out
the eigenvalues and eigenstates of the Hamiltonian.

%% ----------------------------------------------------------------
\subsection{The two-dimensional Laplacian}

In two dimensions, with an $M$-point mesh along $x$ and an $N$-point mesh along
$y$, we have
\begin{align*}
	\Delta \psi\pickpipe{k,\ell} &=
	\frac{\D^2\psi}{\Dx^2}\pickpipe{k,\ell} +
	\frac{\D^2\psi}{\Dy^2}\pickpipe{k,\ell} \\
	&\approx
%
		\frac{
			\psi(t, x_{k+1}, y_\ell) -
			2 \psi(t, x_k, y_\ell) +
			\psi(t, x_{k-1}, y_\ell)
		}{(\Delta x)^2} \\
%
	&+
		\frac{
			\psi(t, x_k, y_{\ell+1}) -
			2 \psi(t, x_k, y_\ell) +
			\psi(t, x_k, y_{\ell-1})
		}{(\Delta y)^2}.
%
\end{align*}
Computing the discrete Laplacian at $(x_k,y_\ell)$ involves sampling $\psi$
there and at four neighbor points:

	\setlength{\unitlength}{1.0in}
	\begin{picture}(6.5,1.7)
	%\put(0,0){\framebox(6.5,1.7)}

	\put(0.60,0.35){\line(1,0){3.0}}
	\put(0.60,0.60){\line(1,0){3.0}}
	\put(0.60,0.85){\line(1,0){3.0}}
	\put(0.60,1.10){\line(1,0){3.0}}
	\put(0.60,1.35){\line(1,0){3.0}}

	\put(1.10,0.10){\line(0,1){1.5}}
	\put(1.60,0.10){\line(0,1){1.5}}
	\put(2.10,0.10){\line(0,1){1.5}}
	\put(2.60,0.10){\line(0,1){1.5}}
	\put(3.10,0.10){\line(0,1){1.5}}

	\put(2.10,0.60){\circle*{.05}}

	\put(1.60,0.85){\circle*{.05}}
	\put(2.10,0.85){\circle*{.05}}
	\put(2.60,0.85){\circle*{.05}}

	\put(2.10,1.10){\circle*{.05}}

	\end{picture}

Again, the question of boundary conditions comes up, and some possible
solutions are as before:  Dirichlet-zero, or periodic (i.e. toric).

%% ----------------------------------------------------------------
\subsection{PDE solver}

The Matlab code is straightforward, using the \texttt{del2} function
for the discrete Laplacian with Dirichlet-zero boundary conditions:

\begin{verbatim}
xlo = -10; xhi =  10; dx  =   0.5;   %% x mesh
ylo = -10; yhi =  10; dy  =   0.5;   %% y mash
[x, y] = meshgrid(xlo:dx:xhi, ylo:dy:yhi);

h   =   0.01; c   =   1.0;
dt  =   0.01; niter = 1000;          %% t mesh

psi0 = exp(-x.^2 + -y.^2);           %% Initial condition
psi0 = psi0 / sqrt(sum(sum(conj(psi0).*psi0))); % Normalize

V = 0.0 * x * y;                     %% Potential
surf(x, y, abs(psi0))                %% Plot the initial wave function

psie = psi0; psir = psi0; t = 0;

for iter = 1:niter
    % Euler:
    psie = psie + i * dt * (c * del2(psie) - V.*psie);

    % RK4:
    k1 = i * (c * del2(psir)          - V .*  psir);
    k2 = i * (c * del2(psir + k1*h/2) - V .* (psir + k1*h/2));
    k3 = i * (c * del2(psir + k2*h/2) - V .* (psir + k2*h/2));
    k4 = i * (c * del2(psir + k3*h)   - V .* (psir + k3*h));
    psir = psir + dt/6*(k1 + 2*k2 + 2*k3 + k4);

    psie = psie / sqrt(sum(sum(conj(psie).*psie))); % Renormalize for
    psir = psir / sqrt(sum(sum(conj(psir).*psir))); %   round-off error

    t = t + dt;
end
surf(x, y, abs(psir))                %% Plot the final wave function
\end{verbatim}

Here's the initial wave function:
\begin{center}
\includegraphics{code/se2d0.eps}
\end{center}
Here's the final wave function (it looks about the same using Euler or
Runge-Kutta):
\begin{center}
\includegraphics{code/se2dr.eps}
\end{center}

Much more could be done with this example, e.g. varying the potential.
Nonetheless, my main point is that the second-derivative term tends to spread
out peaks over time, leading to delocalization of the wave function after
measurement.

%% ================================================================
\newpage
\section{The Heisenberg uncertainty principle}

The following is not central to the talk (and in fact, was not presented in the
talk).  However, I feel that the \plainidx{Heisenberg uncertainty principle} is
a result that many of us have heard of but which is often not precisely
understood.  There is a beautiful discussion in \cite{bib:nc}, which I
elaborate on and modify here.

There is some background which will not be proved here.  First, take it as
given that many observables come in \bfidx{conjugate pairs} --- technically,
this means that one is the Fourier transform of the other.  These include,
among others, position and momentum.

Second, Heisenberg showed that such pairs of observables (call them, say, $A$
and $B$) satisfy the \bfidx{canonical commutation relation}
$$
	[A,B] := AB-BA = i\hbar.
$$

Third, as \cite{bib:nc} point out, there is a common misconception that
measurement of position disturbs a particle, making measurement of its momentum
uncertain.  Disturbance can occur during measurement, but this is not the
content of the uncertainty principle.  Rather, it quantifies the products of
the \emphidx{standard deviations} of the measurements of $A$ and $B$ --- where
the standard deviation of a measurement describes its variability over
measurements made in a large number of experiments.

Fourth, recall from probability that the standard deviation of a quantity $X$
is
$$
	\Delta X := \sqrt{\expn{(X-\expn{X})^2}} = \sqrt{\expn{X^2} - \expn{X}^2}.
$$
If it so happens that $X$ has mean zero, then this reduces to
$$
	\Delta X = \sqrt{\expn{X^2}}.
$$

Fifth, recall from section \ref{sec:expn_braCket} that
the expected value of $X$ is
$$
	\expn{X} = \braCket \psi X \psi.
$$

\begin{prop}[Heisenberg uncertainty principle]
If $A$ and $B$ are a conjugate pair of observables, then
$$
	\Delta A \Delta B \ge \frac{\hbar}{2}.
$$
\begin{proof}
Suppose the system is in the state $\psi$.
To simplify the algebra a bit, let
$$
	C = A - \expn{A} \qqand D = B - \expn{B}.
$$
These have mean zero.  Since $A$ and $A$ are observables, they are Hermitian;
note that $C$ and $D$ are also Hermitian.  Also recall that shifting by
a constant doesn't affect the standard deviation, so
$$
	\Delta C = \Delta A
	\qqand
	\Delta D = \Delta B.
$$
The expected value of the commutator is
\begin{align*}
	\braCket{\psi}{[C,D]}{\psi}
		&= \braCket{\psi}{CD}{\psi} - \braCket{\psi}{DC}{\psi}.
\end{align*}
By the triangle inequality,
\begin{align*}
	|\braCket{\psi}{[C,D]}{\psi}|
		&\le |\braCket{\psi}{CD}{\psi}| + |\braCket{\psi}{DC}{\psi}|.
\end{align*}
The \plainidx{Cauchy-Schwartz inequality} gives us
\begin{align*}
	|\braCket{\psi}{CD}{\psi}| &\le
		%\left(\braCket{\psi}{C^2}{\psi} \braCket{\psi}{D^2}{\psi} \right)^{1/2}
		\sqrt{ \braCket{\psi}{C^2}{\psi} \braCket{\psi}{D^2}{\psi} }
	= \sqrt{ \expn{C^2} \expn{D^2} }
	= \Delta C \Delta D
\end{align*}
and likewise
\begin{align*}
	|\braCket{\psi}{DC}{\psi}| &\le \Delta C \Delta D.
\end{align*}
Then
$$
	\Delta C \Delta D \ge \frac{|\braCket{\psi}{[C,D]}{\psi}|}{2}
$$
but
$$
	\braCket{\psi}{[C,D]}{\psi} = i\hbar \braket{\psi}{\psi} = i\hbar
$$
so
$$
	\Delta C \Delta D \ge \frac{\hbar}{2}.
$$
Since $\Delta A = \Delta C$ and $\Delta B = \Delta D$, we have
$$
	\Delta A \Delta B \ge \frac{\hbar}{2}.
$$

\end{proof}
\end{prop}
