The most successful numerical approach to quantum field theory
begins with a formulation of quantum mechanics developed by
Feynman in which a quantum amplitude is described as a weighted
integral over all possible paths (not necessarily obeying the
classical equations) which start at the system's initial state
and end at the final state. For single particle quantum mechanics
the quantum amplitude 
for a
transition from position 
at time 
to position 
at
time 
is written as:
where 
is the classical action for the path 
given
by
This is a sophisticated Wiener integration over function space and is typically an awkward formalism for analytic calculation. However, it is nicely suited for numerical work since it replaces the normal operator/Hilbert space formalism of quantum mechanics with an explicit integral.
The path integral appropriate for quantum field theory is similar to the equation above except that the integration must be performed over all possible time evolutions of field configurations rather than particle trajectories. In our physical problem, a field configuration specifies both the quark and gluon fields as particular functions of space. A particular time evolution then specifies these fields as functions of space and time. This problem is easily put in a numerically tractable from by replacing the space-time continuum by a grid or lattice of points, conventionally a uniform, four-dimension mesh.
The field theory analogue of the single-particle action given in the equation above is a similar polynomial in the field variables and their derivatives, integrated over space-time. Thus, the corresponding discrete field theory action will be a four-dimensional sum of a local density which depends on the lattice field variables at a specific lattice site and its nearest neighbors.
The actual integration appropriate for the lattice QCD evaluation
of an observable 
is typically performed as a Monte Carlo
average,
over an ensemble of configurations for the gluon fields 
,
.
Each configuration assigns a specific 
complex matrix
U to each link connecting neighboring sites in the lattice. The
ensemble used above is generated by a Metropolis or molecular
dynamics algorithm so as to be distributed according to the
positive definite statistic weight:
The sum is over all elementary squares or plaquettes, 
,
that can be constructed out of four lattice links and 
is the ordered product of the corresponding U matrices
associated with those links. The quark fields correspond to
anti-commuting classical variable and cannot be treated
numerically as an integral but instead are represented by the
determinant above. Here D is a nearest neighbor difference
operator and m, the quark mass matrix. Typically, the force
generated by the determinant is computed using a noisy estimator
which can be done using, for example, a conjugate gradient method
for computing the inverse of the sparse matrix D+m. In addition
to the quark mass matrix m, the coupling strength, related to
the parameter 
, is the only other free parameter in the
calculation.
This is an ideal formulation for massively parallel computing. A
typical large-scale lattice calculation might work with a 
hypercubic lattice. If each processor in a parallel machine is
assigned a 
subvolume, 4K processors would be required. The
most computationally demanding part of a conjugate-gradient
iteration requires about 500 floating point operations per
lattice site or 128K flops/processor. A 3-component complex
vector must be transferred both in and out of the processor for
each link that joins the 
subcube to its neighbors, or 
or 6K words total. This suggests a
reasonably favorable 20:1 computation to communications load for
each processor.