Numerical studies of Hamiltonian SU(N) lattice gauge theory in (3+1)
dimensions have shown that the gauge fields exhibit chaotic
behavior in the classical limit [1]. The numerical value of the
largest positive Lyapunov exponent 
has been obtained for SU(2)
and SU(3) [1,2].
For the SU(2) gauge theory the
complete spectrum of Lyapunov exponents was obtained on small lattices
[3]. These calculations, which follow the evolution of a classical
gauge field configuration in real time, also showed that the
energy density distribution on the lattice rapidly approaches a thermal
distribution [4]. This confirms the expectation of a finite
growth rate of the coarse-grained entropy density of the gauge field,
which follows from the observation that the sum over all positive
Lyapunov exponents at fixed energy density grows like the volume
[3]. Hence, at any given level of coarse-graining, the classical gauge
field ``self-thermalizes'' on a time scale of the order of the inverse
Lyapunov exponent.
It is of course very important to be sure that the analyzed chaotic dynamical behavior, especially the exponentially growing divergence of initially adjacent field configurations is of physical and not numerical origin, i.e. this phenomenon survives the continous time limit. With other words - since exponentially growing deviations in a computer simulation well can have their origin in certain, ultimately eventual properties of the applied numerical algorithms, - it is better to ensure that the numerical solution of the time dependent differential equations describing the evolution of the lattice gauge field theory themselves are stable against small perturbations which are always present due to the finite resolution of real numbers on a computer.
Furthermore --- since it is often more comfortable to use somewhat more
variables in the numerical simulation than the number of independent
physical degrees of freedom, like using the cartesian pair 
instead of an angle 
alone because this facilitates
the calculation
--- time independent constraints (like 
) have
to be fulfilled and conserved during the update procedure.
These requirements are not trivial, especially in the case of Hamiltonian
lattice gauge theories, where the satisfaction of the Gauss law and the
equations of motion in the continous time limit should commute,
but they do not if using
an arbitrary, naive recursion formula for the time update.
The purpose of the present article is to show a class of numerical algorithms which conserve Gauss' law and the unitarity of SU(2) group elements --- represented as real quaternions --- on each link and site of the lattice, respectively. Furthermore, considering a conservative Hamiltonian system, the conservation of the total energy is required.
First we discuss how to conserve the Noether charge by a numerical update procedure which coincides with the Hamiltonian equations of motion in the continous time (i.e. zero time step) limit. Then, we study the stability properties of some representants of such algorithms, tentatively called "odd even time" (OET) algorithm, "implicit midpoint" (IMP) algorithm and "half implicit half explicit endpoint" (HIHEP) algorithm, in the case of harmonic oscillator and exponential growth. Especially the last case is relevant for numerical studies of classically chaotic, constrained, conservative Hamiltonian systems. Here is namely essentially important, that we remain on the allowed surface in the phase space, even approximating the solution of the equations of the motion, or for that purpose, that small deviances orthogonal to this surface tend to decrease. The exact geometry of such a surface is, however, very complicated for lattice gauge field theories.
As a simple representant of constrained systems the motion of a point particle on a circle of unit radius is discussed. The method of Lagrange multipliers is presented in order to create algorithms which describe the above motion in cartesian coordinates and conserve the norm, the unit radius, exactly in each recursion step.
After these studies we turn to lattice gauge theory, but try to keep the discussion general in respect of the local symmetry group and the representation of its elements. Finally the representative algorithm OET is explicitely constructed for the classical SU(2) lattice gauge theory.