# Geometric approach to chaos in the classical dynamics of abelian lattice gauge theory

## I Introduction

Classical dynamical aspects of lattice gauge theories have recently attracted some interest . The classical limit of a lattice gauge theory is interesting both from the point of view of classical dynamical system theory and in the perspective of its quantum counterpart. As far as the latter aspect is concerned, the interest mainly resides in the fact that very few non-perturbative tools are available to study quantum gauge field theories, while in the classical limit it is in principle possible to exactly simulate the real time evolution of the system at any energy. From the point of view of the theory of classical dynamical systems, lattice gauge theories are highly non-trivial many-degree-of-freedom Hamiltonian systems which exhibit a rich and interesting phenomenology. The classical Hamiltonian dynamics of such systems is known to be chaotic ; however, a precise characterization of the different chaotic regimes which may occur in these systems is still lacking.

Many particular aspects of this general problem have been considered in the literature, concerning the properties of pure gauge theories and of theories coupled with matter (mainly Higgs fields), e.g., the systematic study of the Lyapunov spectra , the study of thermalization processes , and the relation between Lyapunov exponents and observable quantities like the plasmon damping rate .

A particular problem which is still open is the dependence of the largest Lyapunov exponent $`\lambda `$ of the pure Yang-Mills lattice theory on the energy density $`\epsilon `$ — energy per plaquette, or energy per degree of freedom — particularly at low $`\epsilon `$ . First, the $`\epsilon `$ scaling of $`\lambda `$ seems different according to the fact that the theory is abelian — $`U(1)`$ — or non-abelian — $`SU(2)`$ and $`SU(3)`$; in the latter case two different scalings have been measured, namely $`\lambda \epsilon ^{1/4}`$  and $`\lambda \epsilon `$ , while in the former case a rapid decrease of $`\lambda `$ at low $`\epsilon `$ was observed . As we shall see in the following, our results suggest that in the $`U(1)`$ case the power law $`\lambda \epsilon ^2`$ holds. As pointed out by Müller and Trayanov  and subsequently by Nielsen et al. , such a problem is interesting because it is tightly related with the problem of the relevance of the chaotic lattice dynamics for the continuum limit of the gauge theory. In fact, intrinsic dimensional arguments can be used to show that the lattice spacing $`a`$, which can be set equal to one in all the numerical simulations after a convenient choice of units, enters the relation between $`\lambda `$ and $`\epsilon `$ as follows:

$$a\lambda (a)=f(a\epsilon (a)),$$
(1)   
hence if one observes numerically a power law $`\lambda \epsilon ^k`$, the latter can be read as

$$\lambda (a)a^{k1}\epsilon (a).$$
(2)   
This means that considering a continuum limit $`a0`$ in which $`\epsilon (a=0)`$ is finite, corresponding to a finite temperature of the resulting field, then the Lyapunov exponent is finite in the limit $`a0`$ only for the particular exponent $`k=1`$. Larger exponents $`k>1`$ would imply that $`lim_{a0}\lambda (a)=0`$, thus leading to a regular continuum theory, while exponents $`k<1`$ would mean that in the continuum theory the Lyapunov exponent diverges. The linear behavior plays then a very special role. As regards non-abelian theories, in Ref.  some evidence is reported supporting the fact that the correct scaling is the linear one, $`\lambda \epsilon `$, the other one ($`\lambda \epsilon ^{1/4}`$) being a spurious result due to finite-time effects. According to these results the continuum limit, for small but finite energy densities, of the Lyapunov exponent of non-abelian lattice gauge theories is finite. In fact, extracting reliable informations from numerical simulations in the low-energy regime is very difficult, mainly because of finite-time effects, which become very important at small energy densities as the characteristic instability time scales grow, and of finite-size effects which are known to be rather large in typical lattice gauge theory simulations.

In the present work we apply a recently proposed formalism , which is based on a Riemannian geometrization of Hamiltonian dynamics, to the classical dynamics of lattice gauge theories. Such a formalism allows one to relate chaotic dynamics and curvature properties of suitable manifolds, and to obtain an analytic formula for the Lyapunov exponent in terms of average curvature fluctuations . The quantities entering this formula are statistical averages which can be computed regardless of the knowledge of the dynamics, either by MonteCarlo or molecular dynamics simulations, or analytically in some cases . As a first step, we apply this formalism to the abelian — $`U(1)`$ — lattice gauge theory, which is the simplest one, leaving the non-abelian case to future work. In the case of a $`U(1)`$ gauge theory we perform a precise numerical measurement of the Lyapunov exponent by simultaneous integration of the Hamilton’s equations and of the tangent dynamics equations using a precise symplectic algorithm  for several lattice sizes. In these simulations we measure also the relevant geometric observables that allow for a characterization of the dynamics and that enter the above-mentioned theoretical expression for the Lyapunov exponent. We find that the analytic estimate compares very well with the outcomes of the numerical simulations for the Lyapunov exponents. Moreover, we find that the theoretical estimate is almost free from finite-size effects; for small energies we find the law $`\lambda \epsilon ^2`$ already when inserting in the formula the values of the geometric observables obtained with small ($`4\times 4\times 4`$) lattices, while the numerical values of the Lyapunov exponents are affected by large finite-size effects.

The paper is organized as follows: in Section II we review very briefly the geometric theory of Hamiltonian chaotic dynamics; in Section III we describe the model and the observables studied. Section IV is devoted to the presentation and the discussion of the results, and in Section V we draw some conclusions and we outline some future developments.

## II Geometry and chaotic dynamics

Let us now recall very briefly the main points about the geometric theory of Hamiltonian chaos. Details can be found in Refs. . Despite the fact that this theory is still at its beginning, it has already proved useful not only in its original context, but also in connection with the problem of the relationship between dynamics and statistical mechanics (in particular phase transitions) .

Hamiltonian dynamics can be rephrased in geometrical terms owing to the fact that the trajectories of a dynamical system with quadratic kinetic energy can be seen as geodesics of a suitable Riemannian manifold. There are several choices for the ambient manifold as well as for the metric tensor. As already discussed in Ref.  a particularly useful ambient space is the enlarged configuration space-time $`M\times ^2`$, i.e. the configuration space $`\{q^1,\mathrm{},q^i,\mathrm{},q^N\}`$ with two additional real coordinates $`q^0`$ and $`q^{N+1}`$. In the following $`q^0`$ will be identified with the time $`t`$. For standard Hamiltonians $`=T+V(𝐪)`$ where $`T=\frac{1}{2}a_{ij}\dot{q}^i\dot{q}^j`$, this manifold, equipped with Eisenhart’s metric $`g_E`$, has a semi-Riemannian (Lorentzian) structure ($`detg_E=1`$). The arc-length is given by

$$ds^2=a_{ij}dq^idq^j2V(𝐪)(dq^0)^2+2dq^0dq^{N+1},$$
(3)   
where both $`i`$ and $`j`$ run between $`1`$ and $`N`$. Let us restrict to geodesics whose arc-length parametrization is affine, i.e. $`ds^2=c_1^2dt^2`$; simple algebra shows that the geodesic equations

$$\frac{d^2q^\mu }{ds^2}+\mathrm{\Gamma }_{\nu \lambda }^\mu \frac{dq^\nu }{ds}\frac{dq^\lambda }{ds}=0\mu ,\nu ,\lambda =0\mathrm{}N+1,$$
(4)   
become Newton equations (without loss of generality $`a_{ij}=\delta _{ij}`$ is considered)

$$\frac{d^2q^i}{dt^2}=\frac{V}{q_i}$$
(5)   
for $`i=1\mathrm{}N`$, together with two extra equations for $`q^0`$ and $`q^{N+1}`$ which can be integrated to yield

$`q^0`$ $`=`$ $`t`$ (7)   
$`q^{N+1}`$ $`=`$ $`{\displaystyle \frac{c_1}{2}}t+c_2{\displaystyle _0^t}L(𝐪,\dot{𝐪})𝑑t`$ (8)   
where $`L(𝐪,\dot{𝐪})`$ is the Lagrangian, and $`c_1`$, $`c_2`$ are real constants. In the following we set $`c_1=1`$ in order that $`ds^2=dt^2`$ on the physical geodesics. As stated by Eisenhart theorem , the dynamical trajectories in configuration space are projections on $`M`$ of the geodesics of $`(M\times ^2,g_E)`$.

In the geometrical framework, the stability of the trajectories is mapped on the stability of the geodesics, hence it can be studied by the Jacobi equation for geodesic deviation

$$\frac{D^2J}{ds^2}+R(\dot{\gamma },J)\dot{\gamma }=0,$$
(9)   
where $`R`$ is the Riemann curvature tensor, $`\dot{\gamma }`$ is the velocity vector along the reference geodesic $`\gamma (s)`$, $`D/ds`$ is the covariant derivative and $`J`$, which measures the deviation between nearby geodesics, is referred to as the Jacobi field. The stability — or instability — of the dynamics, and thus deterministic chaos, originates from the curvature properties of the ambient manifold. In local coordinates, Eq. 9 is written as

$$\frac{D^2J^\mu }{ds^2}+R_{\nu \rho \sigma }^\mu \frac{dq^\nu }{ds}J^\rho \frac{dq^\sigma }{ds}=0,$$
(10)   
and in the case of Eisenhart metric it simplifies to

$$\frac{d^2J^i}{dt^2}+\frac{^2V}{q_iq^j}J^j=0,$$
(11)   
which is nothing but the usual tangent dynamics equation for standard Hamiltonians. The Lyapunov exponents are usually computed evaluating the rate of exponential growth of $`J`$ by means of a numerical integration of Eq. (11) .

In the particular case of constant curvature manifolds, Eq. (9) becomes very simple 

$$\frac{D^2J^\mu }{ds^2}+KJ^\mu =0,$$
(12)   
and has bounded oscillating solutions $`J\mathrm{cos}(\sqrt{K}s)`$ or exponentially unstable solutions $`J\mathrm{exp}(\sqrt{K}s)`$ according to the sign of the constant sectional curvature $`K`$, which is given by

$$K=\frac{K_R}{N1}=\frac{}{N(N1)},$$
(13)   
where $`K_R=R_{\mu \nu }\dot{q}^\mu \dot{q}^\nu `$ is the Ricci curvature and $`=R_\mu ^\mu `$ is the scalar curvature; $`R_{\mu \nu }`$ is the Ricci tensor. Manifolds with $`K<0`$ are considered in abstract ergodic theory (see e.g. Ref. ). Krylov  originally proposed that the presence of some negative curvature could be the mechanism actually at work to make chaos in physical systems, but in realistic cases the curvatures are neither found constant nor everywhere negative, and the straightforward approach based on Eq. (12) does not apply. This is the main reason why Krylov’s ideas remained confined to abstract ergodic theory with few exceptions.

In spite of these major problems, some approximations on Eq. (9) are possible even in the general case. The key point is that negative curvatures are not strictly necessary to make chaos, and that a subtler mechanism related to the bumpiness of the ambient manifold is actually at work. Upon an assumption of quasi-isotropy of the ambient manifold, i.e., that the manifold can be obtained as a small deformation of a constant-curvature space (see Ref.  for details), Eq. (9) can be approximated by an effective scalar equation which reads

$$\frac{d^2\psi }{dt^2}+K(t)\psi =0,$$
(14)   
where $`\psi `$ is a generic component of the vector $`J`$ (in this approximation all the components are considered equivalent), and $`K(t)`$ is a stochastic process which models the curvature along the geodesic curve. Such a stochastic model is defined by

$$K(t)=k_R+\delta ^2k_R^{1/2}\eta (t),$$
(15)   
where $`k_R=K_R/N`$, $``$ stands for an average taken along a geodesic, which, for systems in thermal equilibrium, can be substituted with a statistical average taken with respect to a suitable probability measure (e.g. the micro-canonical or the canonical measure); $`\eta (t)`$ is a stationary $`\delta `$-correlated Gaussian stochastic process with zero mean and variance equal to one. Using Eisenhart metric, and for standard Hamiltonians, the non-vanishing components of the Riemann tensor are $`R_{0i0j}=_{q_i}_{q_j}V`$, hence the Ricci curvature has the remarkably simple form

$$k_R=\frac{1}{N}^2V,$$
(16)   
where $`^2`$ is the Euclidean Laplacian operator. Equation (14) becomes a stochastic differential equation, i.e. the evolution equation of a random oscillator . It is worth noticing that Eq. (14) is no longer dependent on the dynamics, since the random process depends only on statistical averages. The estimate of the Lyapunov exponent $`\lambda `$ is then obtained through the evolution of the second moments of the solution of (14) as

$$\lambda =\underset{t\mathrm{}}{lim}\frac{1}{2t}\mathrm{log}\frac{\psi ^2(t)+\dot{\psi }^2(t)}{\psi ^2(0)+\dot{\psi }^2(0)}.$$
(17)   
As shown in Ref. , this yields the following expression for $`\lambda `$:

$$\lambda (k,\sigma _k,\tau )=\frac{1}{2}\left(\mathrm{\Lambda }\frac{4k}{3\mathrm{\Lambda }}\right),$$
(18)   
where

$$\mathrm{\Lambda }=\left(\sigma _k^2\tau +\sqrt{\frac{64k^3}{27}+\sigma _k^4\tau ^2}\right)^{1/3},$$
(20)   
$$\tau =\frac{\pi \sqrt{k}}{2\sqrt{k(k+\sigma _k)}+\pi \sigma _k};$$
(21)   
in the above expressions $`k`$ is the average Ricci curvature $`k=k_R`$ and $`\sigma _k`$ stands for the mean-square fluctuation of the Ricci curvature, $`\sigma _k=\delta ^2k_R^{1/2}`$.

The advantages in using the geometric approach to Hamiltonian chaos are thus evident. In fact, it is possible to give reliable estimates of the Lyapunov exponent without actually computing the time evolution of the system: the estimate (18) of $`\lambda `$ depends only on statistical averages which can be either computed analytically in some cases (for instance in the case of the FPU model ) or, in general, extracted from a Monte Carlo or a dynamical simulation, as it is the case of the model studied in the present work.

The behavior of the average geometric observables as the control parameter (e.g. the energy density or the temperature) is varied conveys an information which goes beyond the possibility of computing the Lyapunov exponent. In fact, one can look at the random oscillator equation (14) as an effective Jacobi equation for a geodesic flow on a surface $`M`$ whose Gaussian curvature is given by the random process $`K(t)`$. As long as nonlinear coupled oscillators are considered, the average Ricci curvature is positive, hence $`M`$ can be regarded as a sphere with a fluctuating radius. In the limit of vanishing fluctuations, one recovers the bounded evolution of the Jacobi field associated with integrable dynamics. Chaos suddenly appears as curvature fluctuations are turned on, nevertheless it it will be “weak” as long as $`\sigma _kk`$, i.e. as long as $`M`$ can be considered as a weakly perturbed sphere. On the contrary as the size of curvature fluctuations becomes of the same order of the average curvature, $`\sigma _kk`$, $`M`$ can no longer resemble a sphere, and the dynamics will no longer “feel” the integrable limit. Hence we expect the dynamics to be strongly chaotic. This is by no means a deep explanation of the existence of weakly and strongly chaotic regimes in Hamiltonian dynamics. Nevertheless it shows how the simple geometric concepts which enter the Riemannian description of Hamiltonian chaos, besides providing effective computational tools, are also useful in helping one’s physical intuition with images and analogies which would be difficult to find elsewhere.

## III Model and dynamical observables

The dynamical system that we are now considering is the classical lattice gauge theory based on the abelian gauge group $`U(1)`$. The Hamiltonian of such a system can be derived in two ways, either as the dual (Wegner model) of a lattice planar spin system , or starting from the Kogut-Susskind Hamiltonian for a generic gauge group $`G`$  and then specializing to the group $`U(1)`$. In the latter case we obtain the same Hamiltonian as in the former case by choosing the $`SO(2)`$ real matrix representation of the group $`U(1)`$.

The lattice Lagrangian, obtained from the Wilson action by fixing the temporal gauge (all the temporal components of the gauge fields are set to zero) and then by taking the continuum limit in the temporal direction, is

$$=\frac{ag^2}{2}\underset{\mathrm{links}}{}\dot{U}_{x,\mu },\dot{U}_{x,\mu }\frac{1}{ag^2}\underset{\mathrm{plaquettes}}{}\left(1\frac{1}{2}\mathrm{Tr}U_{\mu \nu }\right),$$
(22)   
where $`U_{x,\mu }G`$ is a group element defined on a link of a $`d`$-dimensional cubic lattice, labeled by the site index $`x`$ and the oriented lattice direction $`\mu `$, $`g^2`$ is a coupling constant, $`a`$ is the lattice spacing, $`,`$ stands for the scalar product between group elements, defined as

$$A,B=\frac{1}{2}\mathrm{Tr}(AB^{}),$$
(23)   
and $`U_{\mu \nu }`$ is a shorthand notation for the plaquette operator,

$$U_{\mu \nu }=U_{x,\mu }U_{x+\mu ,\nu }U_{x+\mu +\nu ,\mu }U_{x+\nu ,\nu }.$$
(24)   
We can pass to a standard Hamiltonian formulation by putting

$$P=\frac{}{\dot{U}}=ag^2\dot{U},$$
(25)   
thus obtaining

$$ag^2=\frac{1}{2}\underset{\mathrm{links}}{}P_{x,\mu },P_{x,\mu }+\underset{\mathrm{plaquettes}}{}\left(1\frac{1}{2}\mathrm{Tr}U_{\mu \nu }\right).$$
(26)   
The parameters $`a`$ and $`g^2`$ can be scaled out, so we set $`g=a=1`$.

The Hamiltonian (26) is the classical Hamiltonian for a lattice gauge theory with a generic gauge group $`G`$. Let us now specialize to the abelian group $`G=U(1)`$. Choosing the representation

$$U=\left(\begin{array}{cc}\mathrm{cos}\phi & \mathrm{sin}\phi \\ \mathrm{sin}\phi & cos\phi \end{array}\right)$$
(27)   
we have

$$P=\dot{U}=\dot{\phi }\left(\begin{array}{cc}\mathrm{sin}\phi & \mathrm{cos}\phi \\ \mathrm{cos}\phi & \mathrm{sin}\phi \end{array}\right)$$
(28)   
and we find

$$\frac{1}{2}P,P=\frac{1}{2}\dot{\phi }^2.$$
(29)   
To write the plaquette operator, we use the fact that the group $`U(1)SO(2)`$ is abelian. Then the product of the four rotations is a rotation of the sum of the oriented angles along the plaquette (the discrete curl of the field $`\phi `$)

$$\phi _{x,\mu \nu }=\phi _{x,\mu }+\phi _{x+\mu ,\nu }\phi _{x+\nu ,\mu }\phi _{x,\nu },$$
(30)   
and the magnetic energy, i.e., the potential energy of the dynamical system can be written as

$$V=\underset{\mathrm{plaquettes}}{}\left(1\frac{1}{2}\mathrm{Tr}U_{\mu \nu }\right)=\underset{\mathrm{plaquettes}}{}\left(1\mathrm{cos}\phi _{x,\mu \nu }\right).$$
(31)   
Our canonical variables are then the angles $`\phi _{x,\mu }`$ and the angular momenta $`\pi _{x,\mu }=\dot{\phi }_{x,\mu }`$, and the Hamiltonian (26) becomes

$$=\frac{1}{2}\underset{\mathrm{links}}{}\pi _{x,\mu }^2+\underset{\mathrm{plaquettes}}{}\left(1\mathrm{cos}\phi _{x,\mu \nu }\right).$$
(32)   
Constant-energy (microcanonical) simulations have been performed on three-dimensional lattices — with lattice sizes ranging from $`4^3`$ to $`15^3`$ — by integrating the canonical equations of motion

$`\dot{\phi }_{x,\mu }`$ $`=`$ $`\pi _{x,\mu };`$ (34)   
$`\dot{\pi }_{x,\mu }`$ $`=`$ $`{\displaystyle \frac{V}{\phi _{x,\mu }}},`$ (35)   
where $`x`$ runs over the $`L^3`$ lattice sites and $`\mu =1,2,3`$, by using a precise third-order bilateral symplectic algorithm . We remind that symplectic algorithms are integration schemes which exactly preserve the canonical form of the equations of motion. The exact Hamilton equations are replaced, in the time discretization procedure, by a map that is symplectic, hence the discrete-time flow that approximates the true Hamiltonian flow generated by Hamilton’ s equations is still Hamiltonian. All the geometric constraints on the phase space trajectories which are enforced by the canonical form of the equations of motion are thus exactly preserved during the numerical integration procedure and the Hamiltonian flow is faithfully represented. As a by-product of this features, symplectic algorithms conserve very well the total energy of the system: in our simulations relative energy fluctuations were of the order of $`10^7÷10^8`$. In Eqs. (32) $`V`$ is given by Eq. (31), whose explicit expression on a three-dimensional lattice is

$$V=\underset{x}{}\underset{(\mu \nu )}{}\left(1\mathrm{cos}\phi _{x,\mu \nu }\right),$$
(36)   
where $`(\mu \nu )=12,13,23`$. The forces (rhs of Eq. 35) are given by

$$\frac{V}{\phi _{x,\mu }}=\underset{\delta =1,2}{}\mathrm{sin}\phi _{x\mu \delta ,\mu \mu +\delta }\mathrm{sin}\phi _{x,\mu \mu +\delta }.$$
(37)   
In order to compute the largest Lyapunov exponent $`\lambda `$ by the standard method , the tangent dynamics equations (11), which now reads as

$$\ddot{J}_{x,\mu }+\underset{y}{}\underset{\nu }{}\frac{^2V}{\phi _{x,\mu }\phi _{y,\nu }}J_{y,\nu }=0,$$
(38)   
have been integrated simultaneously with the Hamilton’s equations (32) by means of the same algorithm. The largest Lyapuonov exponent has then been computed according to the definition

$$\lambda =\underset{t\mathrm{}}{lim}\frac{1}{t}\mathrm{log}\frac{\left[|\dot{J}|^2(t)+|J|^2(t)\right]^{1/2}}{\left[|\dot{J}|^2(0)+|J|^2(0)\right]^{1/2}}.$$
(39)   
Here $`|J|^2=_x_\mu J_{x,\mu }^2`$ is the squared Euclidean norm of the tangent vector $`J`$.

According to the discussion of Sec. II, the relevant geometric observable which is able to characterize the chaotic dynamics of the model is the Ricci curvature $`k_R`$, computed with the Eisenhart metric, defined by Eq. (16), which now can be rewritten as

$$k_R=\frac{1}{L^3}\underset{x}{}\underset{\mu }{}\underset{\nu \mu }{}[\mathrm{cos}\phi _{x,\mu \nu }+cos\phi _{x\nu ,\mu \nu }].$$
(40)   
The average and the rms fluctuations of the Ricci curvature $`k_R`$ have been computed in all the simulations. The results are presented and discussed in the next Section.

## IV Results and discussion

The average and the fluctuations of the Ricci curvature (40) are plotted in Fig. 1 against the energy density $`\epsilon `$ for three different lattice sizes. The resulting patterns reveal that in the low-energy regime the average $`k`$ is much larger than the fluctuation $`\sigma _k`$, whereas in the high-energy regime the converse is true. These results are qualitatively the same as those that have been recently obtained for a planar (XY) spin model on a $`2d`$ cubic lattice . Such a fact is not surprising, since the present model, possessing a local gauge invariance under proper planar rotations, is expected to behave on a $`d`$-dimensional lattice — from a statistical-mechanical point of view — as its dual model (invariant under a global symmetry) on a $`d1`$-dimensional lattice . Moreover, this means that our results are consistent with the fact that the system undergoes a Berežinskij-Kosterlitz-Thouless (BKT) transition at a finite energy density (of order $`\epsilon 1`$ in our units). However, the analysis of the geometric quantities is not expected to give sharp indications of the presence of a transition in the case of a BKT transition. On the contrary in the case of second-order transitions the fluctuations of the Ricci curvature exhibit peculiar cusp-like behaviours  which are completely absent here.

From the arguments reported at the end of Sec. II we expect a crossover from weak to strong chaos where $`\sigma _kk`$. The results reported in Fig. 1 thus suggest that such a crossover should occur around $`\epsilon 1`$, perhaps at a somewhat higher value of $`\epsilon `$ in the case of a $`15^3`$ lattice. Actually, as it is shown in Fig. 2, in the energy region around $`\epsilon 1`$ the numerically computed values of the largest Lyapunov exponent rapidly grow. At larger energies the values of $`\lambda `$ decrease because as $`\epsilon \mathrm{}`$ the model becomes integrable. At variance with what is usually found in other models (where the $`\lambda (\epsilon )`$ pattern is rather stable when the number of degrees of freedom of the system is varied) sizeable variations in the $`\lambda (\epsilon )`$ pattern are observed here at different lattice sizes. Again, the phenomenology is very close to the one observed in the XY model on a $`2d`$ lattice. The fact that there is a change in the chaotic properties of the dynamics is well evident looking at Fig. 3, where the Lyapunov exponents are plotted against $`\epsilon `$ in a log-log scale. In fact, in an energy range around $`\epsilon 1`$, the behaviour of $`\lambda `$ as a function of the energy density deviates from a steep power law to a smoother one. In Fig. 3 the numerical values of $`\lambda `$ are compared with the theoretical estimates obtained according to the geometric theory outlined in Sec. II, i.e., obtained by substituting the computed values of $`k`$ and $`\sigma _k`$, shown in Fig. 1, into Eq. (18).

Two facts are immediately evident from Fig. 3. First, in the low-energy region (weak chaos) the theoretical estimates show a power-law behaviour $`\lambda \epsilon ^2`$ already using the geometric values computed using a very small ($`4^3`$) lattice, and the exponent of the power law remains the same for all the lattice sizes, whereas the numerical values of $`\lambda `$ show a less steep $`\epsilon `$-dependence for small lattices. Second, both the low-energy power law behaviour and the actual values of the numerical Lyapunov exponents are the closer to the theoretical estimates the larger the lattice size is. Let us mention that the actual value of $`\lambda `$ is theoretically underestimated in the transition region, as it is well evident from Fig. 3. Without entering the details of this mismatch, we mention that the reason can be understood in the light of a very similar situation found in $`1d`$, $`2d`$ and $`3d`$ XY models .

Moreover, the actual values of the theoretical estimates are almost free from finite-size effects, which are instead very large in the numerical values of the Lyapunov exponents. This fact is particulalrly evident at very low energy densities, as shown in Fig. 4, where the numerical and theoretical values of $`\lambda `$ are plotted against the lattice size.

## V Concluding remarks

The results presented in the previous section clearly show that the geometric approach to Hamiltonian chaotic dynamics outlined in Sec.II is particularly efficient to provide very good estimates of the largest Lyapunov exponent for a classical lattice gauge theory with abelian gauge symmetry. In particular, the theoretical estimates of the Lyapunov exponents are almost free from finite-size effects, especially at low energies. The theoretical values are very close to the numerical values that one extrapolates for very large lattice sizes even when the geometric quantities are computed using very small lattices. Also the dependence of $`\lambda `$ on the energy density in the weakly chaotic (low-energy) region is $`\lambda \epsilon ^2`$ already for a $`4^3`$ lattice, where the numerical values of the $`\lambda `$’s exhibit a much less steep scaling.

This is a confirmation of the fact that the geometric estimates of the Lyapunov exponents obtained by Eq. (18) are very good for large systems, as already observed in other cases . It is not surprising since the whole theory uses the large size of the systems as a hypothesis . Moreover, this feature is very promising in the perspective of applying the geometric theory to non-abelian lattice gauge theories. In fact, as already mentioned in the Introduction, in the non-abelian case there is still some doubt about the correct value of the exponent $`\alpha `$ of the power law $`\lambda \epsilon ^\alpha `$ in the low-energy region, and this has remarkable importance for the relevance of chaos in the continuum limit of the theory. The application of the geometric theory of chaos might perhaps help in obtaining reliable estimates of the exponent $`\alpha `$ already using small lattice sizes.

###### Acknowledgements.

We thank G. Pettini for fruitful discussions. LC acknowledges also useful discussions with A. Giansanti. The present work was mainly carried out during one of the authors’ (LC) stay at the Département de Physique Théorique de l’Université de Genève which is gratefully acknowledged for the kind hospitality. This work is part of the EEC project CHRXCT94/0579 (OFES 95.0200).