
\section{\label{problems}Three problems}

     The choice of problems here reflects strongly an original
motivation for this work, namely to ``do QCD'' in a functional
language.  {\em Quantum Chromodynamics}, the theory of the strong nuclear
force, is a relativistic quantum field theory.  Through a series of
mathematical tricks\footnote{The Feynman path-integral formulation of
quantum mechanics \cite{feynman}; analytic continuation
to Euclidean space-time (Wick rotation); and Wilson's lattice
discretisation of gauge theories \cite{wilson1,wilson2,wilson3}}, which
reduce $d$-dimensional quantum theories to $(d + 1)$-dimensional
statistical mechanics problems, it can be expressed as a {\em lattice
gauge theory}.  In this 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 calculation the
electric field produced by a given distribution of charges.  For an
isolated point charge, the electric field is described (in 3-dimensions)
by the inverse-square law.  For a more complicated system of charges,
with non-trivial boundary conditions (for example, the charges may be
enclosed in a metal box) there is no simple solution.  The
electrostatic potential is obtained by solving the associated Poisson
equation.  For simplicity we describe the two-dimensional case---a
third dimension doesn't introduce anything 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 a $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-Siedel 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 \cite{multigrid}.

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

     Statistical physics is largely concerned with the study of {\em
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}.
It has an unusual
transition driven by ``condensation of vortex configurations''.
The kind of model of interest here 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\footnote{Taking the partial derivative of $E (u)$ with respect to
$u (x, y)$, for some fixed $(x, y)$, and setting the result to zero,
recovers equation (\ref{discrete}).}
\begin{displaymath}
  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{displaymath}
$E (u)$ is the electrostatic energy of the configuration
$u$; the lowest energy state is, so to speak, the one ``preferred by
nature''\footnote{The squared terms in $E$ are discrete difference
approximations to the squared electric field, $\mbox{\boldmath $F$}^2$.}.

     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\footnote{As $\theta
\rightarrow 0$, states with energy above the minimum are exponentially
suppressed, and the simulation yields the ground state with probability
approaching 1---an observation underlying the {\em simulated
annealing} method of minimisation.}.  A phase transition generally
manifests itself in discontinuous behaviour of some expectation-values
as $\theta$ passes through a critical value.

     In the XY model $u (x, y)$ isn't a single real number.
Instead it is a unit-length two-dimensional vector.  
This vector can be represented as a complex number with unit
magnitude.  The energy function can then be written as
\begin{displaymath}
\label{XY}
  E (u) = - \Re \sum_{x, y} u^* (x, y) u (x + 1, y) +
                            u^* (x, y) u (x, y + 1),
\end{displaymath}
where $u^*$ means complex conjugate of $u$, and $\Re$ means
real part\footnote{Alternatively, each vector can be parametrised
by its angle from some fixed direction.  This gives the
equivalent expression $E (u) = - \sum \cos (\theta (x + 1, y) - \theta (x, y))
+ \cos (\theta (x, y + 1) - \theta (x, y))$}.
It isn't difficult to see that $E$ is minimised when
every vector is pointing in the same direction as its neighbours.
In other words, when {\em all} vectors point in the same direction.

     Some more interesting configurations are also very stable in the XY
model.  A few of these ``vortex'' configurations, which are local
(rather than absolute) minima of the energy function, are illustracted
in figure \ref{vortices}.
\begin{figure}[hbt]
\centerline{
\begin{picture}(120, 100)(-20, -20)
\multiput(-10,   0)( 0, 20){4}{\line(1, 0){80}}
\multiput(  0, -10)(20,  0){4}{\line(0, 1){80}}
\thicklines
\put( 0,  0){\line(-1, 1){10}}
\put( 0, 20){\line(-1, 2){10}}
\put( 0, 40){\line( 1, 2){10}}
\put( 0, 60){\line( 1, 1){10}}
\put(20,  0){\line(-2, 1){10}}
\put(20, 20){\line(-1, 1){10}}
\put(20, 40){\line( 1, 1){10}}
\put(20, 60){\line( 2, 1){10}}
\put(40,  0){\line(-2,-1){10}}
\put(40, 20){\line(-1,-1){10}}
\put(40, 40){\line( 1,-1){10}}
\put(40, 60){\line( 2,-1){10}}
\put(60,  0){\line(-1,-1){10}}
\put(60, 20){\line(-1,-2){10}}
\put(60, 40){\line( 1,-2){10}}
\put(60, 60){\line( 1,-1){10}}
\end{picture}
\begin{picture}(120, 100)(-20, -20)
\multiput(-10,   0)( 0, 20){4}{\line(1, 0){80}}
\multiput(  0, -10)(20,  0){4}{\line(0, 1){80}}
\thicklines
\put( 0,  0){\line(-1,-1){10}}
\put( 0, 20){\line(-2,-1){10}}
\put( 0, 40){\line(-2, 1){10}}
\put( 0, 60){\line(-1, 1){10}}
\put(20,  0){\line(-1,-2){10}}
\put(20, 20){\line(-1,-1){10}}
\put(20, 40){\line(-1, 1){10}}
\put(20, 60){\line(-1, 2){10}}
\put(40,  0){\line( 1,-2){10}}
\put(40, 20){\line( 1,-1){10}}
\put(40, 40){\line( 1, 1){10}}
\put(40, 60){\line( 1, 2){10}}
\put(60,  0){\line( 1,-1){10}}
\put(60, 20){\line( 2,-1){10}}
\put(60, 40){\line( 2, 1){10}}
\put(60, 60){\line( 1, 1){10}}
\end{picture}
}
\caption{\label{vortices}Vortex configurations of the XY model}
\end{figure}

     Algorithms for generating sequences of configurations with Boltzmann
distribution for this type of system resemble relaxation.  They are
iterative algithms, 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\footnote{The procedure 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 $exp
(-\Delta E)$.  If the candidate is rejected, retain the original value
of $u$.\footnote{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\footnote{$U(1)$ is the group of unit magnitude complex numbers,
with complex multiplication as the product}, associated
with QED ({\em Quantum Electrodynamics}).  QED is closely
related to QCD; the modifications needed to go to QCD
are outlined at the end of the section.
\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{XY}) 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.  The energy function is
\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.\footnote{In the XY model energy, one of the factors in the link product
has to by conjugated.  In gauge theory, variables associated
with links which 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 which effect 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}.
\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}

     For real physics, the lattice should be four dimensional.
The procedure described here simulates the ``gauge-particles''
of a theory (photons or gluons).  It doesn't take account of
matter particles.  The modification required to handle QCD rather
than QED is that the link variables should be elements in
the group, $SU (3)$, of unit-determinant, $3 \times 3$, unitary matrices,
rather than simply unit-magnitude complex numbers.

