To put QCD on a computer, we proceed as follows [Wilson:74a],
[Creutz:83a]. The four-dimensional
space-time continuum is replaced by a four-dimensional hypercubic periodic
lattice of size 
, with the
quarks living on the sites and the gluons living on the links of the lattice.
is the spatial and 
is the temporal extent of the lattice. The
lattice has a finite spacing a. The gluons are represented by
complex SU(3) matrices associated with each link in the
lattice. The 3 in SU(3) reflects the fact that there are three colors of
quarks, and SU means that the matrices are unitary with unit determinant
(i.e., ``special unitary''). This link matrix describes how the color
of a quark changes as it moves from one site to the next. For example, as a
quark is transported along a link of the lattice it can change its color
from, say, red to green; hence, a red quark at one end of the link can
exchange colors with a green quark at the other end. The action functional
for the purely gluonic part of QCD is
where 
is a coupling constant and
is the product of link matrices around an elementary square or plaquette on
the lattice-see Figure 4.8. Essentially all of the time in
QCD simulations of gluons is spent multiplying these SU(3) matrices
together. The main component of this is the 
kernel, which
most supercomputers can do very efficiently. As the action involves
interactions around plaquettes, in order to satisfy detailed
balance we can update only half the links in any one
dimension simultaneously, as shown in Figure 4.9 (in two dimensions
for simplicity). The partition function for full-lattice QCD including quarks
is
where 
is a large sparse matrix the
size of the lattice squared. Unfortunately, since the quark or fermion
variables 
are anticommuting Grassmann numbers, there is no
simple representation for them on the computer. Instead, they must be
integrated out, leaving a highly non-local fermion determinant:
This is the basic integral one wants to evaluate numerically.
Figure 4.8: A Lattice Plaquette
Figure 4.9: Updating the Lattice
The biggest stumbling block preventing large QCD simulations with
quarks is the presence of the determinant 
in the
partition function. There have been many proposals for dealing with the
determinant. The first algorithms tried to compute the change in the
determinant when a single link variable was updated [Weingarten:81a].
This turned out to be prohibitively expensive. So instead, the approximate
method of pseudo-fermions [Fucito:81a] was used.
Today, however, the preferred approach is the so-called Hybrid Monte
Carlo algorithm [Duane:87a], which is exact.
The basic idea is to invent some dynamics for the variables in the system in
order to evolve the whole system forward in (simulation) time, and then do a
Metropolis accept/reject for the entire evolution on the basis of the total
energy change. The great advantage is that the whole system is updated in
one fell swoop. The disadvantage is that if the dynamics are not correct,
the acceptance will be very small. Fortunately (and this is one of very few
fortuitous happenings where fermions are concerned), good dynamics can be
found: the Hybrid algorithm [Duane:85a]. This is a neat combination
of the deterministic microcanonical method [Callaway:83a],
[Polonyi:83a] and the stochastic Langevin method [Parisi:81a],
[Batrouni:85a], which yields a quickly evolving, ergodic algorithm for
both gauge fields and fermions. The computational kernel of this algorithm
is the repeated solution of systems of equations of the form
where 
and 
are vectors that live on the sites of the
lattice. To solve these equations, one typically uses a conjugate
gradient algorithm or one of its cousins,
since the fermion matrix 
is sparse. For more
details, see [Gupta:88a]. Such iterative matrix algorithms have
as their basic component the 
kernel, so again
computers which do this efficiently will run QCD well.