Lyapunov instabilities in lattices of interacting classical spins at infinite temperature
We numerically investigate Lyapunov instabilities for one-, two- and three-dimensional lattices of interacting classical spins at infinite temperature. We obtain the largest Lyapunov exponents for a very large variety of nearest-neighbor spin-spin interactions and complete Lyapunov spectra in a few selected cases. We investigate the dependence of the largest Lyapunov exponents and whole Lyapunov spectra on the lattice size and find that both quickly become size-independent. Finally, we analyze the dependence of the largest Lyapunov exponents on the anisotropy of spin-spin interaction with the particular focus on the difference between bipartite and nonbipartite lattices.
Investigations of the Lyapunov instabilities in many-particle systems are often motivated by the role of chaos in the foundations of statistical physics [1, 2, 3], which is still not fully understood. In general, interacting many-particle classical systems are expected to be chaotic. The largest Lyapunov exponents and properties of the entire Lyapunov spectra have been calculated numerically for classical many-body systems, such as gases of hard-core particles [4, 5], fluids with soft interactions , and lattice two-dimensional rotators [6, 7, 8, 9, 10], and analytically in a few cases [11, 12, 13, 14, 15, 16, 17, 8, 9, 10]. In the present paper, we focus on lattices of interacting classical spins.
Classical spins often appear in the theoretical studies as the large-spin limit of quantum spins. Moreover, even when one deals with lattices of spins 1/2, parallels between classical and quantum spin dynamics still remain. These parallels have recently received much attention, in particular, in the context of the asymptotic exponential-oscillatory behavior of nuclear spin decays in solids [18, 19, 20, 21, 22, 23, 24]. These decays were identified in Refs. [18, 20] with chaotic eigenmodes in both classical and quantum many-spin systems.
Although it appears very likely a priori that lattices of interacting classical spins exhibit chaotic dynamics, no systematic investigation of the chaotic properties of these lattice was undertaken until our previous work , which presented a survey of the largest Lyapunov exponents for a very large variety of spin lattices and Hamiltonian anisotropies. The principal finding of Ref.  was that all Hamiltonians considered, with the exception of the Ising case, led to chaotic dynamics as evidenced by the nonzero value of the largest Lyapunov exponent. We also obtained both analytically and numerically the power-law scaling of the largest Lyapunov exponent in the vicinity of the integrable Ising limit.
In the present paper, we complement the findings of Ref.  in several respects. Namely, we compute complete Lyapunov spectra for a few selected spin lattices and show their dependence on the lattice size, and also present a more extensive investigation of the lattice size dependence of the largest Lyapunov exponent. Finally, we discuss the dependence of the largest Lyapunov exponent on the Hamiltonian anisotropy with particular emphasis on the difference between bipartite and nonbipartite lattices.
2 General formulation
2.1 Spin model
We consider periodically closed spin lattices with the nearest-neighbor (NN) interaction Hamiltonian of the following kind:
where are the three projections of the classical spin vector of unit length on the th lattice site (i.e. ), and , , are the coupling constants, which we also normalize by condition . Below, we often mention Ising, Heisenberg and “anti-Heisenberg” limits of the Hamiltonian (1). The Ising Hamiltonian corresponds to , , Heisenberg Hamiltonian to , and, finally, anti-Heisenberg Hamiltonian to . The total number of spins in the lattice is denoted as . The phase space of such a lattice has dimensionality .
We consider seven lattices shown in Fig. 1 and labeled as (L1-L7). Lattices (L1-L5) are bipartite, which means that they can be divided into two sublattices such that all the interacting neighbors for a spin on one sublattice belong to the other sublattice. Lattices (L6,L7) are nonbipartite.
The equations of motion associated with the Hamiltonian (1) can be obtained in the Poisson-bracket formalism [33, 34, 35]: , where the index admits values , or representing the projections , , or , respectively. The primary Poisson brackets are: , where is the Kronecker delta and the Levi-Civita symbol. The evaluation of involves the following Poisson brackets: . The resulting equations of motion are
where is the local field given by the expression
Here , and are the unit vectors along the respective directions, and implies the summation of the nearest neighbors of the -th lattice site.
In this work, we restrict ourselves to the Lyapunov instabilities on the zero energy shell, which corresponds to infinite temperature in the microcanonical sense.
2.2 Analytical considerations
Since the dynamics of this system are fully time-reversible, positive and negative Lyapunov exponents are expected to form conjugate pairs with equal absolute values. In the general case of different , , and , there should be only two zero Lyapunov exponents, corresponding to the energy and the time shift directions. In the case , of which the anti-Heisenberg case is an example, the -component of the total spin polarization becomes an integral of motion, and hence one more pair of Lyapunov exponents should also become equal to zero. In the Heisenberg case, , all three components of the total spin polarization become integrals of motion. However, these three integrals of motion are not dynamically independent, because one of them can be obtained as the Poisson bracket of the two others. The two independent integrals of motion thus imply two additional pairs of zero Lyapunov exponents. In the Ising case, the -component of each spin is an integral of motion. Hence, the system is integrable and all Lyapunov exponents are equal to zero.
Now we discuss the connections between the Lyapunov spectra of different anisotropic Hamiltonians. We first note that infinite-temperature Lyapunov spectra are expected to be identical for the Hamiltonians with the coupling constants and , because the zero-energy shells in the two cases are identical, while the change of the sign of the coupling constants amounts to the operation of time-reversal and flips the sign of the energy. We further note that, for bipartite lattices, the infinite temperature Lyapunov spectra for the coupling constants and should also be identical. In order to see this, one should make the transformation and for one of the two sublattices forming the bipartite lattice and then examine the resulting equations of motion.
The two observations made above also imply that, for a given bipartite lattice, the Lyapunov spectra for the Heisenberg and the anti-Heisenberg Hamiltonians are identical to each other, which means that, in the latter case, the Lyapunov spectrum has three pairs of zero Lyapunov exponents instead of two expected for the generic case of . The extra pair of the zero exponents originates from the following two (not independent) integrals of motion: or , where is the index taking value 0 for one sublattice and 1 for the other. The same considerations also imply that any Hamiltonian with , or equivalent, on a bipartite lattice can be converted to the axially symmetric Hamiltonian with . Therefore, the original Hamiltonain has an extra pair of zero Lyapunov exponents corresponding to the integral of motion .
Small spin clusters may have additional nontrivial integrals of motion. In particular, any 4-spin periodic chain with the general anisotropic Hamiltonian of form (1) is fully integrable. This is because the first and third spins in this chain rotate in the same local field , while the second and the fourth spins rotate in another local field . Therefore, there are two additional integrals of motion, namely: and . Finally, one can also check that is also an integral of motion. As a result, the number of integrals of motion (including energy) becomes 4, while the dimensionality of the phase space is 8, i.e. the problem is fully integrable. In the case of the isotropic Heisenberg Hamiltonian, a much larger variety of small spin clusters with nontrivial integrals of motion were cataloged in Ref. .
2.3 Numerical simulations
The equations of motion were integrated using a fourth-order Runge-Kutta algorithm with a time-step . This is sufficiently small so that on the time scale of our simulations, energy is conserved with the 6-digit accuracy. Simulations were run for 20000 time units, which was sufficient for accurate convergence of the smaller exponents.
The Lyapunov exponents were obtained by using a version of the standard reorthonormalization algorithm [4, 27, 28]. In order to obtain the first largest Lyapunov exponents, we numerically propagate a reference trajectory and initially orthogonal perturbation vectors . At every time step, we propagate the perturbations using the linear tangent space map that is obtained by numerically taking derivatives,
After each time interval , we hierarchically reorthogonalize the perturbation vectors using the Gram-Schmidt procedure, and then renormalize their lengths back to the initial values. The renormalization factors are recorded for the th time interval as . Finally, we compute the th Lyapunov exponent using the formula , where is the number of the renormalization time intervals.
As a test of the accuracy of our numerical routine, we have checked that the Lyapunov exponents form conjugate pairs and that the number of zero Lyapunov exponents is equal to the number expected from the symmetry of the Hamiltonian as discussed in Section 2.2. We have also checked the symmetries between Heisenberg and anti-Heisenberg Lyapunov spectra expected based on the arguments presented in Section 2.2.
In order to obtain initial conditions at zero total energy, we first choose random orientations for all spins. This produces a total energy close to zero, but with fluctuations of the order of . Then, in order to arrive at zero total energy, we evolve the dissipative equations of motion:
Depending on the sign in front of the right-hand side, this increases or decreases the total energy associated with the Hamiltonian (1). Once the zero value of the total energy is reached, we additionally assure the randomness of the initial conditions on the energy shell by performing sequential rotations of random spins by random angles around the directions of their respective local fields given by Eq.(3). These rotations preserve the total energy.
3 Results and discussion
3.1 Lyapunov spectra
We have computed the full Lyapunov spectra for the linear chain (L1), cubic lattice (L5) and triangular lattice (L7) with the Heisenberg Hamiltoninan, and also for the triangular lattice (L7) with the anti-Heisenberg Hamiltonian. The results are presented in Fig. 2. The spectra (i.e. the Lyapunov exponents as a function the exponent’s index) are typically weakly convex, regardless of the lattice. In the case of a cubic lattice with Heisenberg interaction, the spectrum is actually very close to linear. As explained in Section 2.2, the Lyapunov spectra for the bipartite lattices (L1) and (L5) with anti-Heisenberg Hamiltonian are identical to those already presented in Fig. 2 for the Heisenberg Hamiltonian. The effect of the lattice size on the Lyapunov spectrum for the linear chain with Heisenberg Hamiltonian is shown in Fig. 3. The size dependence of the spectrum appears to be very weak for and becomes virtually unobservable for chains containing more than 32 spins.
Comparing the Lyapunov spectra of classical spins with the Lyapunov spectra of other many-particle systems, we first note, that the spectra presented in Fig. 2 do not exhibit an offset from zero for the smallest positive exponents. Such an offset was observed in gases of particles  and high-dimensional billiards [29, 15] but not in systems with sufficiently soft interactions . The classical spin lattices obviously belong to the latter group.
We further remark on the existence of the delocalized Lyapunov-Goldstone modes, which were observed in dilute gases [4, 30, 31, 14] and in some other extended systems . If these modes exist in the spin systems, the projections on single spins of the Lyapunov vectors corresponding to the smallest nonzero Lyapunov exponents should exhibit a sinusoidal dependence on the positions of spins. In the present work, we did not investigate the properties of the Lyapunov vectors systematically. However, our several attempts to find the Lyapunov-Goldstone modes for the infinite temperature energy shells did not produce any positive evidence: the projections of the Lyapunov vectors for all exponents were strongly localized on the lattice and not sinusoidal. There is also no indication of a dependence of the smalles non-vanishing exponent on the system length, as is characteristic for the Lyapunov-Goldstone modes. The Lyapunov-Goldstone modes may still exist in classical spin systems at low temperatures but this is a subject beyond the scope of the present paper.
3.2 Largest Lyapunov exponents: finite size effects
Accurate numerical calculations of full Lyapunov spectra are very demanding, because the computational cost grows as . For this reason, the lattices investigated in the preceding subsection were relatively small. In this subsection, we focus only on the largest Lyapunov exponents, . The cost of computing grows only as , which allows us to investigate much larger lattices and a much greater variety of Hamiltonians. An extensive investigation of this kind was already reported by us in Ref. . In the present and the next subsections we present some results and analysis that were not included in Ref. .
In this subsection, we investigate the dependence of on the lattice size for all seven lattices shown in Fig. 1 with Heisenberg and the anti-Heisenberg Hamiltonians. The results are shown in Fig. 4. They indicate that, within our numerical accuracy, becomes size-independent for sufficiently large lattices. In particular, on the basis of these results even slow logarithmic growth of with the lattice size can be excluded.
The saturation of with the lattice size can be explained on the basis of the following consideration. The exponential growth of the perturbation vector in many-spin phase space with rate implies that the projection on the subspace of each individual spin should also, on average, grow exponentially with the same rate. The instantaneous growth rate of the perturbations of the coordinates of a given spin can be obtained from the linearized equations of motion, which, in turn can be obtained from equations (2) and (3),
Here, indicates a perturbation in the -th spin, while indicates a perturbation in its th component. From these equations, one can see that the instantaneous growth rate is limited from above by the value on the order of the maximum value of the local field , where is the number of nearest neighbors for each lattice site. Since this constraint does not depend on the lattice size, the growth of the perturbation belonging to the largest Lyapunov exponent must saturate as the lattice size increases. In principle, it might also be possible for the largest Lyapunov exponent to oscillate with the lattice size, but in a system with exponential decay of spatial correlation, this is extremely unlikely.
3.3 Largest Lyapunov exponents: dependence on the Hamiltonian anisotropy
In Ref. , we conducted a systematic survey of the dependence of on the Hamiltonian anisotropy on the “interaction sphere” constrained by the condition . We have found that the principal parameter controlling this dependence is . In particular, this parameter quantifies the approach to the integrable Ising case corresponding to . The main plot of Ref.  is reproduced in Fig. 5.
Here we focus on the difference between the anisotropy dependence of for bipartite lattices (L1-L5) and nonbipartite lattices (L6, L7). As can be seen in Fig. 5, the bipartite lattices (L1-L5) show a nearly universal dependence , which scales only with the number of interacting neighbors. The small spread of the sampled values of at a given value of indicates that depends only very little on the ratio of the two coupling constants that have smaller absolute values. Our investigation of the nonbipartite lattices (L6,L7) was, in fact, motivated by the impression that the number of interacting neighbors alone determines the entire anisotropy dependence of . We wanted to compare for the lattices (L6) and (L3), where, in the both cases, each site has four nearest neighbors, and for the lattices (L7) and (L5), each site has six nearest neighbors.
We found, however, that, as seen in Fig. 5, the nonbipartite lattices (L6) and (L7) exhibit a noticeable fork-like spread of as approaches . The upper and the lower tips of the fork correspond to the anti-Heisenberg and Heisenberg Hamiltonians, respectively. In the anti-Heisenberg case, the value of is close to that of a bipartite lattice with the same number of nearest neighbors. In the Heisenberg case, the value of is closer to the bipartite lattice with one fewer nearest neighbor. This spread indicates that the knowledge of alone is insufficient to determine . The two-dimensional anisotropy dependence behind this spread is shown in Fig. 6, where we present it for lattices (L5) and (L7) in the form of the color density plots as a function of and .
Less obvious from Fig. 5, is the fact that the significant majority of the sampled values of for the lattices (L6) and (L7) agree very well with the dependence for the bipartite lattices (L3) and (L5), respectively. The comparison between Figs. 6(a) and (b) clearly shows that the difference between lattices (L5) and (L7) is only pronounced, when all three interaction constants have the same sign and roughly the same value—or, in other words, when they approach the Heisenberg limit. This implies that for the anti-Heisenberg Hamiltonian has a more typical value than for the Heisenberg Hamiltonian, and that the expected universality of for the same number of nearest neighbors still roughly holds even for nonbipartite lattices. It thus appears that the conservation of the total spin, i.e. , in the Heisenberg case leads to the reduction of the effective number of the nearest neighbors by roughly one, as far as the value of is concerned.
We suspect that the situation here is similar to the origin of the frustrated low-temperature magnetism for the Heisenberg model on nonbipartite lattices. The interaction energy for a spin pair is minimal when the two spins are antiparallel to each other, but, on a nonbipartite lattice, such an antiparallel configuration cannot simultaneously exist for all pairs of interacting spins. Hence the ground state of a nonbipartite lattice is frustrated. In the case of the Lyapunov instabilities, the conservation of the total spin polarization implies that the perturbation vector corresponding to the largest Lyapunov exponent does not grow along the direction of the total spin polarization, i.e . At the same time, grows, on average, exponentially, which means that, for , . This, in turn, implies that, in the leading order, . As a result, when a given projection grows, this growth needs to be compensated by the growth of for other spins in the opposite direction. However, since the interaction is local, these “other spins” can only be the nearest neighbors. Achieving the maximum growth of the perturbation, and hence the largest value of the Lyapunov exponent presumably requires that the perturbation vector maximizes the anti-alignment of for the adjacent sites. For the bipartite lattices, this anti-alignment can, in principle, be made perfect, but for the nonbipartite lattices this is impossible. Such an explanation is consistent with the fact that for lattices (L6,L7) in the Heisenberg limit is smaller than in the anti-Heisenberg limit. It seems also to be connected to the fact that the above reduction approximately leads to the value of for a bipartite lattice with roughly one fewer nearest neighbor per site.
We finally remark that the overall small spread of values of for bipartite lattices at a given value of is related to the symmetries described in Section 2.2, which imply that the Lyapunov exponents are identical for eight combinations of the coupling constants characterized by the same value of , namely: , , , , , , , . The values of cannot change much between these eight points.
4 Summary and conclusions
We have investigated numerically the Lyapunov spectra of systems of many classical spins for a variety of lattices and coupling constants at infinite temperature. The possibility of varying the coupling constants, from the highly symmetric isotropic Heisenberg model, through partially symmetric couplings such as the anti-Heisenberg model, to the completely anisotropic integrable case, makes these systems particularly interesting. We have presented: (i) calculations of the Lyapunov spectra for selected lattices of interacting classical spins; (ii) investigations of the lattice-size dependence of the Lyapunov spectra; (iii) investigations the largest Lyapunov exponents for a broader group of coupling constants, lattices and large lattice sizes; and (iv) discussion of the difference between the largest Lyapunov exponents for bipartite and nonbipartite lattices. The computed Lyapunov spectra were found to be weakly convex. We have observed no finite offset of the smallest positive exponents, and, to the extent that we have searched, we have not encountered any evidence of Lyapunov-Goldstone modes. Both the largest Lyapunov exponents and the whole Lyapunov spectra were found to become independent of the lattice sizes for sufficiently large lattices. We have given an analytical argument explaining this finding. In addition, we have found that the largest Lyapunov exponents for bipartite and nonbipartite lattices depend differently on the anisotropy of the coupling. This is due to the special symmetry of the bipartite lattices with respect to the sign change of the two out of three coupling constants.
A.S.dW’s work is financially supported by an Unga Forskare grant from the Swedish Research Council. B.V.F. acknowledges the hospitality of the Kavli Institute for Theoretical Physics at the University of California, Santa-Barbara, where a part of this paper was written, and support by the National Science foundation under Grant No. NSF PHY11-25915. The numerical part of this work was performed at the bwGRiD computing cluster at the University of Heidelberg.
-  J. W. Gibbs. Elementary Principles in Statistical Mechanics. Princeton University Press, Princeton, 1979.
-  N. S. Krylov. Works on the Foundations of Statistical Physics. Yale University Press, New Haven, 1902.
-  P. Gaspard. Chaos, Scattering and Statistical Mechanics. Cambridge University Press, Cambridge, 1998.
-  H. A. Posch and R. Hirschl. Simulations of billiards and of hard body fluids. In D. Szasz, editor, Hard Ball Systems and the Lorentz Gas, volume 101 of Encyclopedia of Mathematical Sciences., page 280. Springer-Verlag, New York, 2000.
-  C. Dellago and H. A. Posch. Kolmogorov-Sinai entropy and Lyapunov spectra of a hard-sphere gas. Physica A, 240:68, 1997.
-  W. G. Hoover, H. A. Posch, C. Forster, C. Dellago, and M. Zhou. Lyapunov modes of two-dimensional many-body systems; soft disks, hard disks, and rotors. J. Stat. Phys., 109:765, 2002.
-  Vito Latora, Andrea Rapisarda, and Stefano Ruffo. Lyapunov instability and finite size effects in a system with long-range forces. Phys. Rev. Lett., 80:692–695, 1998.
-  R. O. Vallejos and C. Anteneodo. Theoretical estimates for the largest Lyapunov exponent of many-particle systems. Phys. Rev. E, 66:021110, 2002.
-  C. Anteneodo, R. N. P. Maia, and R. O. Vallejos Lyapunov exponent of many-particle systems: Testing the stochastic approach. Phys. Rev. E, 68:036120, 2003.
-  R. O. Vallejos and C. Anteneodo. Largest Lyapunov exponent of long-range XY systems. Physica A, 340:178, 2004.
-  R. van Zon, H. van Beijeren, and C. Dellago. Largest Lyapunov exponent for many particle systems at low densities. Phys. Rev. Lett., 80:2035, 1998.
-  H. van Beijeren, J. R. Dorfman, H. A. Posch, and C. Dellago. The Kolmogorov-Sinai entropy for dilute gases in equilibrium. Phys. Rev. E, 56:5272, 1997.
-  A. S. de Wijn and H. van Beijeren. Goldstone modes in Lyapunov spectra of hard sphere systems. Phys. Rev. E, 70:016207, 2004.
-  S. McNamara and M. Mareschal. On the origin of the hydrodynamic Lyapunov modes. Phys. Rev. E, 64:051103, 2001.
-  A. S. de Wijn. Lyapunov spectra of billiards with cylindrical scatterers: comparison with many-particle systems. Phys. Rev. E, 72:026216, 2005.
-  A. S. de Wijn and H. van Beijeren. A radius of curvature approach to the kolmogorovâsinai entropy of dilute hard particles in equilibrium. J. Stat. Mech.: Theory and Experiment, P08012, 2011.
-  D. M. Barnett, T. Tajima, K. Nishihara, Y. Ueshima, and H. Furukawa. Lyapunov exponent of a many body system and its transport coefficients. Phys. Rev. Lett., 76:1812, 1996.
-  B. V. Fine. Long-Time Relaxation on Spin Lattice as a Manifestation of Chaotic Dynamics. Int. J. Mod. Phys. B, 18:1119–1159, 2004.
-  Boris V. Fine. Universal long-time relaxation on lattices of classical spins: Markovian behavior on non-markovian timescales. J. Stat. Phys., 112:319–327, 2003.
-  B. V. Fine. Long-time behavior of spin echo. Phys. Rev. Lett., 94:247601, 2005.
-  S. W. Morgan, B. V. Fine, and B. Saam. Universal long-time behavior of nuclear spin decays in a solid. Phys. Rev. Lett., 101:067601, 2008.
-  E. G. Sorte, B. V. Fine, and B. Saam. Long-time behavior of nuclear spin decays in various lattices. Phys. Rev. B, 83:064302, 2011.
-  E. G. Sorte, B. V. Fine, and B. Saam. Phase relationship between the long-time beats of free induction decays and spin echoes in solids. Phys. Rev. B, 85:174425, 2012.
-  B. Meier, J. Kohlrautz, and J. Haase. Eigenmodes in the long-time behavior of a coupled spin system measured by nuclear magnetic resonance. Phys. Rev. Lett., 108:177602, 2012.
-  A. S. de Wijn, B. Hess, and B. V. Fine. Largest lyapunov exponents for lattices of interacting classical spins. Phys. Rev. Lett., 109:034101, 2012.
-  R. Steinigeweg and H.-J. Schmidt. Heisenberg-integrable spin systems. Math. Phys. Anal. Geom., 12:19–45, 2009.
-  G. Benettin, L. Galgani, A. Giorgilli, and J. M. Strelcyn. Lyapunov exponents for smooth dynamical systems and hamiltonian systems; a method for computing all of them, part i: Theory; part ii: Numerical application. Meccanica, 15:9–20, 1980.
-  T. A. Elsayed, B. Hess, and B. V. Fine. arXiv:1112.3626.
-  A. S. de Wijn and H. van Beijeren. Lyapunov spectrum of the many-dimensional dilute random Lorentz gas. Phys. Rev. E, 70:036209, 2004.
-  T. Taniguchi and G. P. Morriss. Localized behavior in the lyapunov vectors for quasi-one-dimensional many-hard-disk systems. Phys. Rev. E, 68:046203, 2003.
-  A. S. de Wijn and H. van Beijeren. Goldstone modes in Lyapunov spectra of hard sphere systems. Phys. Rev. E, 70:016207, 2004.
-  Hong-liu Yang and Günter Radons. When can one observe good hydrodynamic lyapunov modes? Phys. Rev. Lett., 100:024101, Jan 2008.
-  M. Steiner, J. Villain, and C. G. Windsor. Theoretical and experimental studies on one-dimensional magnets. Adv. Phys., 25:87, 1976.
-  K.-H. Yang and J. O. Hirschfelder. Generalizations of classical Poisson brackets to include spin. Phys. Rev. A, 22:1814â1816, 1980.
-  B. Hess. Investigations of Lyapunov spectra for classical spin lattices. Diploma thesis, University of Heidelberg, 2010. http://www.thphys.uni-heidelberg.de/ fine/Hess_Diploma_2010.pdf.