
\section{\label{problems}Three problems}

     The choice of problems here reflects an early motivation for this
work, which was to implement some QCD codes in a functional language.
{\em Quantum Chromodynamics} is the theory of the strong nuclear
force.  It has a reputation as a ``grand challenge'' numerical
problem.  In the present work (as it turned out) the emphasis is more
on elegant specification of the underlying algorithms than high
performance.

     Through a series of mathematical tricks
\cite{feynman,wilson1,wilson2,wilson3} QCD can be expressed as a {\em
lattice gauge theory}, in which form it is amenable to Monte Carlo
simulation.  In this article, lattice gauge theory will be introduced
from a more algorithmic point of view, assuming no particular knowledge
of the physical systems involved.  First two simpler problems which
share some of its features are introduced.

\subsection{\label{poisson}The 2-dimensional Poisson equation}

     A fundamental problem in electrostatics is calculating the
electric field produced by a given distribution of charges.  The
electrostatic potential is obtained by solving the associated Poisson
equation.  For simplicity we describe the two-dimensional case---a
third dimension introduces nothing essentially new.  The given
charge distribution is a function $c (x, y)$ of the coordinates.  The
electrostatic potential $u (x, y)$ is the solution of the equation
\begin{equation}
\label{continuum}
(\mbox{\boldmath $\nabla$}^2 u) (x, y) = - c (x, y).
\end{equation}
Having solved for $u$, the electric field is the gradient
$\mbox{\boldmath $F$} (x, y) = - (\mbox{\boldmath $\nabla$} u) (x, y)$.

     The most simple-minded approach to solving (\ref{continuum}),
is to replace it by the finite difference approximation
\begin{equation}
\label{discrete}
u (x + 1, y) + u (x - 1, y) +
u (x, y + 1) + u (x, y - 1) - 4 u (x, y) = - c (x, y).
\end{equation}
The coordinates $x$, $y$ are now integers---assume that they and $c$
have been rescaled to eliminate appearance of the grid spacing parameter.
Rearranging (\ref{discrete}) as
\begin{equation}
\label{relaxation}
u (x, y) = \frac{1}{4} (\mbox{sum neighbouring $u$ values} + c (x, y)),
\end{equation}
immediately suggests the simplest iterative algorithm for solving this
system.  In the {\em relaxation} algorithm, starting from some initial guess
for $u$, each iteration replaces the value of $u$ at a point by the average
value of its neighbours (plus an offset determined by the local
charge).  Variants differ in whether the $u$ values inserted in the RHS
of (\ref{relaxation}) are those from the previous sweep of the grid (Jacobi
relaxation) or whether they are the most recent values, which may
include some values already updated in the current sweep
(Gauss-Seidel relaxation).  The last version has further variations
depending on the order in which the grid is traversed---a
popular scheme is ``red-black relaxation'' in which all even sites are
visited first, then all odd sites.  Straight relaxation is
a slow algorithm, but it does form the core update of more powerful
methods---see for example \cite{multigrid}.

\subsection{\label{xy}The XY Model}

     Statistical physics is largely concerned with the study of
phase transitions in complex systems.  Many simplified mathematical
models have been introduced which exhibit transitions similar to those
of physical systems.  The {\em XY model} is an example \cite{xyModel}.
This kind of model can be motivated by referring
back to the Poisson equation.  Equations (\ref{continuum}) or
(\ref{discrete}) can be obtained from variational principles.
Solving the system (\ref{discrete}) is equivalent to minimising the
functional
%\begin{equation}
%\label{electroEnergy}
%  E (u) = \sum_{x, y} \frac{1}{2} \left((u (x + 1, y) - u (x, y))^2 + 
%				        (u (x, y + 1) - u (x, y))^2\right) -
%	  c (x, y) u (x, y).
%\end{equation}
\begin{eqnarray}
  E (u) & = & \sum_{x, y} \frac{1}{2} \left((u (x + 1, y) - u (x, y))^2 + 
		      (u (x, y + 1) - u (x, y))^2\right) - \nonumber \\
\label{electroEnergy}
	&   & \hspace{1in} c (x, y) u (x, y).
\end{eqnarray}
$E (u)$ is the electrostatic energy of the configuration $u$; the
lowest energy state is the one ``preferred by nature''.

     In statistical physics, states with higher energy are also taken into
account.  The aim of a simulation is to generate a sequence of
configurations of the system, in which the frequency of appearance of a
particular configuration $u$ with energy $E (u)$ is weighted by the
Boltzmann probability distribution:
\begin{displaymath}
  P (u) \propto \exp (- E (u) / \theta).
\end{displaymath}
Here, $\theta$ is the temperature of the system.

     In the XY model $u (x, y)$ is a unit-length two-dimensional vector
rather than a simple real number.  This vector can be represented as a
complex number with unit magnitude.  The electrostatic energy function,
(\ref{electroEnergy}), is replaced by
\begin{equation}
\label{XYEnergy}
  E (u) = - \Re \sum_{x, y} u^* (x, y) u (x + 1, y) +
                            u^* (x, y) u (x, y + 1),
\end{equation}
where $u^*$ means complex conjugate of $u$, and $\Re$ means real part.

     Algorithms for generating sequences of configurations with Boltzmann
distribution in this type of system resemble relaxation.  They are
iterative algorithms, involving repeated sweeps across the grid, and the
update of a single site normally only needs knowledge of the current
values at neighbouring sites.

     In the {\em Metropolis} algorithm \cite{metropolis} 
the update of a site is as
follows:  choose a candidate new value for the local $u$ by some
pseudo-random procedure.  The procedure used for this choice should give equal
probabilities for rotating the vector by $+\delta$ and $-\delta$, for
all angles $\delta$.  Also it should be ``ergodic''---by combining
some sequence of allowed updates it must be possible to get
from any state to any other.  Calculate the change in energy, $\Delta E$,
which would be produced by replacing the old value of $u$ by the
candidate new one (this depends on the state of the neighbouring
sites).  If $\Delta E$ is negative, accept the candidate update.  If
$\Delta E$ is positive, accept the candidate with probability $P = exp
(-\Delta E / \theta)$.  An operational procedure for doing this is to
generate a pseudo-random value $x$ uniformly in the interval $[0, 1)$,
accepting the candidate if $x < P$.  If the candidate is rejected,
retain the original value of $u$.

The proof of this algorithm is omitted.  The crux
is to show that, when the probability of a system state is given by the
Boltzmann formula, the expected frequency of transitions
(produced by the algorithm) from any state $A$ to any other state $B$
is the same as that from $B$ to $A$ (``detailed balance'') so
that the Boltzmann distribution represents equilibrium.

\subsection{\label{lgt}Lattice gauge theory}

     This section describes the so-called $U(1)$ lattice gauge
theory.  $U(1)$ is the group of unit magnitude complex numbers,
with complex multiplication as the product.  This is the ``gauge
group'' of {\em Quantum Electrodynamics} rather than QCD, but
the two theories are closely related.
\begin{figure}[hbt]
\centerline{
\begin{picture}(80, 90)(-20, -30)
\multiput(-10,  0)( 0, 20){3}{\line(1, 0){60}}
\multiput( 0, -10)(20,  0){3}{\line(0, 1){60}}
\put(0, 0){\circle*{6}}
\put(10, -30){Site}
\end{picture}
\begin{picture}(80, 90)(-20, -30)
\multiput(-10,  0)( 0, 20){3}{\line(1, 0){60}}
\multiput( 0, -10)(20,  0){3}{\line(0, 1){60}}
\put(2,  0){\rule[-1pt]{16pt}{3pt}}
\put(10, -30){Link}
\end{picture}
\begin{picture}(80, 90)(-20, -30)
\multiput(-10,  0)( 0, 20){3}{\line(1, 0){60}}
\multiput( 0, -10)(20,  0){3}{\line(0, 1){60}}
\put(2,  2){\rule{16pt}{16pt}}
\put(0, -30){Plaquette}
\end{picture}
}
\caption{\label{cells}Some elementary cells of a lattice.}
\end{figure}
     The main difference between $U(1)$ gauge theory and the XY
model is geometric.  In the XY model the variables $u$ are
associated with the {\em sites} of the lattice.  The energy
function (\ref{XYEnergy}) is a sum over the {\em links} of the
lattice---it involves multiplying together the values of the
variables at the ends of each link.

     In gauge theory the variables are associated with the links
of the lattice, and the energy function is computed on the
next-higher dimension cells of the lattice---the elementary
squares, or {\em plaquettes} (figure \ref{cells}).

     We can associate a group element (a unit-magnitude complex number)
$U_P$ with any path $P$ in the lattice by multiplying 
together the link variables along the path.  In particular, we
can associate a path with each elementary plaquette---the four-link
path bounding that plaquette.  The energy function is now
\begin{displaymath}
E (u) = - \sum_{P \in \mbox{\em plaquettes}} \Re (U_P)
\end{displaymath}
where {\em plaquettes} is the set of elementary plaquette-paths in the
lattice.  (In the XY model energy, one of the factors in the link product
has to by conjugated.  In gauge theory, variables associated
with links that have negative incidence on the plaquette are conjugated.
This will be illustrated more fully in the next two sections.)

     The Metropolis algorithm can be carried over from the XY model.
The local update is identical.  The difference lies in how the
neighbouring values that affect the local contribution to the
energy are collected.  In the XY model, the local contribution
to the energy was determined by the values of site variables
at neighbouring sites.  In gauge theory it is determined
by link variables on the ``staples'' illustrated (in the 3-dimensional
case) in figure \ref{staples}.  For real physics, the lattice should
be four-dimensional.  It is usually taken as a 4d cubic grid, but the
theory has sometimes been formulated on other lattice geometries;
the implementations given here will be largely independent of
the details of this geometry.
\begin{figure}[hbt]
\centerline{
\begin{picture}(120, 100)(0, -20)
\thicklines
\put(60, 0){\line(0, 1){60}}
\thinlines
\put(65, 5){\line(1, 1){20}}
\put(65, 65){\line(1, 1){20}}
\put(85, 25){\line(0, 1){60}}
\put(55, -5){\line(-1, -1){20}}
\put(55, 55){\line(-1, -1){20}}
\put(35, -25){\line(0, 1){60}}
\put(0, 0){\line(0, 1){60}}
\put(0, 0){\line(1, 0){50}}
\put(0, 60){\line(1, 0){50}}
\put(120, 0){\line(0, 1){60}}
\put(70, 0){\line(1, 0){50}}
\put(70, 60){\line(1, 0){50}}
\end{picture}
}
\caption{\label{staples}Surrounding links, contributing to the energy of the central one.}
\end{figure}

