# Chaotic properties of spin lattices near second-order phase transitions

## Abstract

We perform a numerical investigation of the Lyapunov spectra of chaotic dynamics in lattices of classical spins in the vicinity of second-order ferromagnetic and antiferromagnetic phase transitions. On the basis of this investigation, we identify a characteristic of the shape of the Lyapunov spectra, the “G-index”, which exhibits a sharp peak as a function of temperature at the phase transition, provided the order parameter is capable of sufficiently strong dynamic fluctuations. As a part of this work, we also propose a general numerical algorithm for determining the temperature in many-particle systems, where kinetic energy is not defined.

###### pacs:

05.45.Jn, 75.10.Hk 75.30.Kz 05.20.-y 05.45.Pq 05.50.+q## I Introduction

The notion of chaos is often invoked in statistical physics to justify the ergodicity assumption. However, the relation between the primary characteristics of chaos, namely the Lyapunov exponents, and the equilibrium properties of many-particle systems still remains elusive. A particularly interesting issue in this regard is the temperature dependence of Lyapunov exponents near phase transitions. Lyapunov exponents characterize the sensitivity of phase-space trajectories with respect to small deviations of initial conditions, while phase transitions are often accompanied by strong fluctuations and the appearance of long-range correlations. Because of the striking nature of phase transitions, one may wonder if such dramatic changes involve chaos as well and if there are universal features to be found in the Lyapunov exponents close to phase transitions.

Most of the relevant investigations so far have been limited to the largest Lyapunov exponents and often reported dramatic signatures of phase transitions in their temperature dependences Butera and Caravati (1987); Mehra and Ramaswamy (1997); Kwon and Park (1997); Caiani et al. (1997, 1998a, 1998b); Firpo (1998); Latora et al. (1998); Casetti et al. (1999); Barre and Dauxois (2001); Cerruti-Sola et al. (2000); Casetti et al. (2000); Posch et al. (1990). It should be noted, however, that, for some systems, these signatures likely originate from the infinite range of particle-particle interactions Latora et al. (1998); Firpo (1998), while, for others, as we explain later, they are not intrinsic to Lyapunov exponents but rather reflect the nonanalytic behavior of the temperature with respect to the total energy near a phase transition and would disappear if Lyapunov exponents are plotted as functions of the total energy. At the same time, the investigations of Refs.Caiani et al. (1997, 1998a, 1998b) (reviewed in Ref.Casetti et al. (2000) ) indicated that a quantity closely related to Lyapunov exponents, namely, the curvature of the configuration space in the geometrical formulation of the dynamics, exhibits sharply increasing fluctuations at phase transitions.

In general, Hamiltonian dynamics in an -dimensional phase space generates not one but Lyapunov exponents organized in pairs of equal absolute values and opposite signs. The entire Lyapunov spectra have been investigated so far only across first-order phase transitions Dellago and Posch (1996); Dellago and Posch (1997); Bosetti and Posch (2014). In this paper, we present a detailed investigation of Lyapunov spectra as a function of temperature for lattices of classical spins with nearest-neighbor interaction in the vicinity of ferromagnetic (FM) and antiferromagnetic (AF) second-order phase transitions. We introduce a characteristic of the shape of the Lyapunov spectra, namely the “-index”, which exhibits a peak at the phase transition both as a function of temperature and energy, provided the order parameter is capable of sufficiently strong dynamic fluctuations. The present work builds on our earlier investigations of Lyapunov instabilities in classical spin lattices at infinite temperature de Wijn et al. (2012, 2013), where, in particular, we showed that the lattices are all chaotic with the exception of the Ising case.

## Ii General formulation and numerical aspects

We consider cubic lattices of classical spins with periodic boundary conditions and the nearest-neighbor interaction Hamiltonian:

(1) |

where are the three projections of the th classical spin normalized by the condition , and and are the coupling constants, which we choose such that . The notation indicates the nearest neighbors of the -th spin. Various alternatives to the Hamiltonian in Eq. (1) will be discussed in Sec. VI.

Our procedure for computing the spectrum of Lyapunov exponents is described in Ref. de Wijn et al. (2013). It follows the standard approach of Ref. Benettin et al. (1980). Index in the above notation orders the Lyapunov exponents in decreasing order, with being the largest positive Lyapunov exponent. Due to the demanding nature of the numerical calculations of the full Lyapunov spectrum (order per time step combined with long convergence times for small exponents), we have had to restrict ourselves to lattices of . For this system size, the finite-size effects on the Lyapunov exponents are already small de Wijn et al. (2013) (for more details, see Appendix B).

We numerically integrate the equations of motion associated with the Hamiltonian (1), , where is the local field. Here , and are orthogonal unit vectors. We use a fourth-order Runge-Kutta algorithm with time step 0.005. During the time of our simulations, typically equal to 20000, the total energy is conserved with absolute accuracy better than . The initial conditions corresponding to a given value of the total energy of the system are selected using the routine described in Ref. de Wijn et al. (2013), which draws them from a uniform distribution of the energy shell.

## Iii Determining temperature and identifying phase transitions

The total energy determines the temperature of the system. However, this temperature cannot be found using the average kinetic energy per particle, because the spin Hamiltonian cannot be decomposed into a quadratic-in-momentum kinetic energy and a momentum-independent potential energy Ben (). Below we describe a more general algorithm applicable to any system with smooth dynamics and short-range interactions. (A system-specific alternative to this algorithm would be to simulate a thermal bath or to extract temperature from the correlations between spin polarizations and local fields.)

Our algorithm is based on the definition , where is the entropy of the system, which, is, in turn, defined (after setting ) as . Here is the -dimensional volume of the energy shell in the -dimensional many-particle phase space. The above definitions lead to

(2) |

which implies that the volume of the energy shell changes nearly exponentially as a function of energy with the characteristic constant equal to the inverse temperature. For this reason, obtaining the above constant by random Monte-Carlo sampling of the entire many-particle phase space is not feasible. Instead, our algorithm consists of the following three steps: (i) It locates one point on any given energy shell using a dissipative dynamics routine introduced in Ref. de Wijn et al. (2013). (ii) It randomly samples that energy shell using sequential energy-conserving rotations of randomly chosen spins around the directions of their local fields by random angles. (iii) Finally, it explores the vicinity of each thus obtained point on the energy shell by tiny energy non-conserving rotations of each spin around a randomly chosen axis perpendicular to spin’s direction. The small angles for these rotations are drawn from a Gaussian distribution around zero. (For the lattices considered, we used a standard deviation of rad.) Since grows exponentially with energy, an energy increase as a result of step (iii) is more likely than an energy decrease. We recover the value of temperature by first obtaining the mean and the mean-squared changes of energy, and , respectively, and then substituting them into the formula

(3) |

which is derived in Appendix A.

After obtaining , we find the specific heat as . In the thermodynamic limit, exhibits lambda-point singularity at the FM and AF phase transitions. Since this singularity is washed out by the finite-size effects, we identify the phase-transition temperature with the maximum of — see figure 1. In particular, for an lattice with the Heisenberg Hamiltonian, we thereby obtain , which is the same as the appropriately rounded thermodynamic value that can be extracted from Refs. Deng et al. (2005); Brown and Ciftan (2006). The size dependence of the specific heat is illustrated in Appendix B.

In our systems, corresponds to infinite temperature, while and correspond to positive and negative temperatures respectively. Cubic spin lattices with nearest-neighbor interactions are bipartite, in the sense that they can be divided into two sublattices such that spins of one sublattice interact only with the spins of the other sublattice. The reversal of all spin coordinates for one sublattice changes the sign of while leaving the volume of the corresponding phase space elements the same. As a result, the volumes of energy shells are symmetric with respect to , i.e. . This symmetry implies that, if an AF transition occurs at temperature , then, in the same system, an FM transition occurs at temperature . We define the order parameters as , where the FM order implies all signs , while the AF order implies +1 and -1 alternating between adjacent lattice sites.

Despite the above symmetry of , the Lyapunov spectra are, in general, not symmetric with respect to , because the reversal of all three projections of a spin does not preserve their Poisson brackets (see Ref. de Wijn et al. (2013)) and hence changes the character of the dynamics. The only symmetric case is the -interaction characterized by and . In this case, one can reverse only - and -components for the spins of one of the two sublattices without reversing their -components, thereby protecting the Poisson brackets and, at the same time, reversing the energy. We further note, that, as illustrated in Appendix C, the Lyapunov spectra of bipartite lattices do not change under the simultaneous sign reversal of energy and the sign of one of the three coupling constants.

## Iv Numerical results: G-index

We now turn to the results of our simulations for the case of the Heisenberg interaction at positive temperatures. Fig. 2 shows the energy and temperature dependences of and the Kolmogorov-Sinai entropy (equal to the sum of all positive Lyapunov exponents) together with the specific heat and the order parameter. Several examples of complete Lyapunov spectra are presented in Fig. 3.

Comparing Figs. 2 (a) and (b), we observe that both and exhibit a steep change across the AF phase transition as functions of temperature but not as functions of energy. This behavior does not change with system size (see Appendix B). In general, such behavior is expected for any smooth function of energy , which is then converted to a function of temperature . For the latter function, . Since exhibits a singularity at the phase transition, so does . In other words, the steep changes of and around as such indicate only the change of the energy-temperature relation rather than an intrinsic sensitivity of Lyapunov instabilities to phase transitions.

The examples of spectra shown in Fig. 3, nevertheless indicate that the phase transition influences the shape of the Lyapunov spectra: the closer the temperature to , the more curved the spectrum. We quantify this shape change by a simple ratio, which we call the “-index”:

(4) |

It represents the ratio of the total area between the spectrum and the diagonal line extending in Fig. 3 from to , divided by the area under the spectrum. The -index is plotted in Fig. 2. It exhibits a sharp peak at the phase transition as a function of temperature and also a clear maximum at the corresponding energy. The size dependence of is illustrated in Appendix B.

Motivated by the above finding, we have systematically investigated the -index for other interaction parameters. Fig. 4 presents , , and at positive and negative temperatures for the Heisenberg interaction [(a) and (b)], generic anisotropic interaction [(c) and (d)], a less anisotropic interaction [(e) and (f)], and the -interaction [(g) and (h)]. The AF transitions in all these cases occurred at positive temperatures, while the FM transitions occurred at negative temperatures. Typically, as shown in Figs. 4 (a)–(f) exhibited a sharp peak at either FM or AF phase transitions but never at both. The -interaction was the only case in which we observed no peak in at either of the two phase transitions. A further example illustrating the symmetry of the -index with respect to the simulataneous sign change of the total energy and one of the coupling constants is given in Appendix C.

The behavior of away from the phase transition, in particular the appearance of humps of in Figs. 4 (b), (g), and (h), may also be of interest, but it extends beyond the scope of the present work. Here we only make two remarks: (i) This kind of humps should be distinguished from the “peaks” of the -index associated with the phase transition. In the thermodynamic limit, the “peaks” are expected to have discontinuous first derivatives and hence be cusp-like. This is a consequence of the earlier general argument about the conversion from energy to temperature dependencies near the second-order phase transitions. On the contrary, the “humps” away from the second-order phase transitions are expected to remain broad and smooth in the thermodynamics limit. While the above distinction is reasonably supported by our numerical results, the computational resources available to us were not sufficient to check the scaling of the -index peaks near the phase transition beyond the results presented in Fig. 10 of Appendix B. (ii) Our calculations far into the ordered phases for generic anisotropic couplings of the type presented in Figs. 4 (c), (d), (e), and (f) exhibited very slow convergence — probably because of the formation of magnetic domains. In these two cases, we were not able to check whether humps similar to those seen in Figs. 4 (b), (g), and (h) exist at sufficiently low temperatures.

## V Relation between the -index and the Lyapunov vectors

Now we turn to explaining the presence or the absence of the peaks of at . In general, all positive Lyapunov exponents tend to decrease with decreasing , because the phase space volume available to the system becomes smaller. The function given by Eq. (4) is sensitive to the difference between the temperature dependences of and the average positive Lyapunov exponent . Let us follow the behavior of starting from infinite temperature and then decreasing . As can be seen in Figs. 4 (a), (d), and (f), exhibits a peak at when initially decreases more slowly than and then drops faster around , thereby catching up with . We now propose an argument, which we later substantiate by examples, that the above behavior of is due to the fact that the order parameter is capable of strong dynamical fluctuations. In such a case, the Lyapunov vector corresponding to seeks the directions in the phase space corresponding to the faster-than-average dynamics, which are, in turn, correlated with the combinations of variables contributing to . In the opposite case, when is not capable of sufficiently strong dynamical fluctuations, the Lyapunov vector corresponding to ignores the respective directions in the phase space. In such a case, and exhibit very similar behavior over the entire range of temperatures seen in Figs. 4 (b), (c), (e), (g) and (h), and, as a result, does not have a peak at .

In order to exemplify the notion of strong dynamical fluctuations of the order parameter, let us assume that the magnetic order sets in along the -axis. [This is the only possibile direction for the interaction used for Figs. 4 (c), (d), (e), and (f), or one of a continuous set of possible directions for Figs. 4 (a), (b), (g) and (h).] Let us then decompose the Hamiltonian as

(5) |

where we use spin raising and lowering variables and , which are analogous to the raising and lowering quantum spin operators Fine (1997, 2004). We refer to the first two terms in the right-hand-side of Eq. (5) as “double-flip” and “flip-flop” terms respectively. The first of them changes the -projections of the two spins in the same direction, while the second one changes them in the opposite directions. The flip-flop term makes AF order fluctuate, while conserving the FM order. The double-flip term has the opposite effect.

For the Heisenberg Hamiltonian, . Therefore, the double-flip term is zero, while the flip-flop term dominates. As a result, the AF order strongly fluctuates in time near the phase transition, which, according to our argument, leads to the peak of seen in Fig. 4 (a). On the contrary, the FM order that sets in at negative temperatures does not fluctuate in time. Accordingly, does not exhibit a peak at in Fig. 4 (b).

For the Hamiltonian corresponding to Figs. 4 (c) – (f), . Therefore, the double-flip term dominates. This leads to the peak of at the FM transition and no peak at the AF transition.

For the XX-interaction corresponding to Figs. 4 (e, d), , i.e. the flip-flop and the double-flip terms have equal stength. This implies that the Lyapunov vector corresponding to does not particularly seek either FM or AF correlations. As a result, there are no peaks of at either FM or AF transition.

The above interpretation is supported by our Fourier analysis of the components of Lyapunov vectors de Wijn et al. (2013). Here are the cubic lattice indices. We compute the function , where

(6) |

and are the wave numbers of the discrete Fourier components.

for the Heisenberg case is presented in Fig. 5. The left column of this figure shows for the Lyapunov vector corresponding to at the temperature of the AF phase transition. For comparison, the middle column represents for the Lyapunov exponent taken from the middle of the positive side of the Lyapunov spectrum at the same temperature and the right column corresponds to but at the (negative) temperature of the FM transition. In the first case, the bright spots in Fig. 5 around indicate strong AF correlations. In the latter two cases, no correlations of AF or FM type are apparent.

## Vi Possible generalisations

The question arises as to how general are the results obtained in this article. One possible generalization of the system considered here is a lattice of interacting rigid rotators. In this case, the Hamiltonian would depend on the rotation angle and the angular momentum of each rotator. The dynamics of such a system would include several aspects different from the classical spin lattices considered in this work. Firstly, the dimensionality of the phase space and hence the number of the Lyapunov exponents would be two times larger. Secondly, since the kinetic energy of the rigid rotator system is not limited from above, it cannot have negative temperature. We do not expect the -index of the rigid rotator lattices to be the same as for the classical spin lattices considered in this article. At the same time, lattices of rigid rotators are also likely to exhibit a second-order phase transitions. We would then expect that the -index in this case, would also exhibit a peak at the phase transition for some but not for all possible interaction models. Such an investigation, however, extends beyond the scope of the present article. Another interesting potential investigation would be to compute the -index for the classical liquid-gas system near the critical point on the pressure-temperature phase diagram, where the line of the first-order phase transition ends.

## Vii Conclusions

To summarize, we have identified a characteristics of Lyapunov spectra of many-spin systems — the -index — which, as a function of temperature exhibits a clear peak at magnetic phase transitions, provided the variable associated with the order parameter is capable of strong dynamic fluctuations. We expect similar behavior near second-order phase transitions in other many-particle systems with short-range interactions. As a part of this work, we have also developed an algorithm for determining microcanonical temperatures of general Hamiltonian systems.

Note added: Recently, we discovered that a significant part of the justification of our temperature-determining algorithm [namely, roughly that up to Eq. (A6) in Appendix A] was done in Ref. Rugh (1997).

###### Acknowledgements.

The authors are grateful to T. A. Elsayed for helpful discussions during the initial stage of this work. A.S.dW’s work is financially supported by an Unga Forskare grant from the Swedish Research Council. The numerical part of this work was performed at the bwGRiD computing cluster at the University of Heidelberg.## Appendix A Derivation of Eq.(3) for temperature associated with a given energy shell

Let us start by mentioning that an intuitive insight in the forthcoming general derivation can be gained by considering an example of -dimensional Euclidean phase space, and assuming that the energy is given by the distance to the origin of a Cartesian coordinate system in this space. In this case, the family of energy shells becomes a continuous set of -dimensional hyperspherical surfaces with a common center.

Turning to the general case, let us denote the complete set of coordinates in -dimensional many-particle phase space as . An energy shell corresponding to energy is defined by condition

(7) |

The corresponding -dimensional phase-space volume is . Let us further consider a small element of volume on this energy shell. If the energy changes by value , the above element can be bijectively mapped onto an element of the new energy shell by moving in the direction orthogonal to the original energy shell. The change of the coordinates in this case is

(8) |

where

(9) |

Here is the vector orthogonal to the original energy shell at a given point. The individual components of this vector are . The volume of the above element of the energy surface after transformation (8) is

(10) |

where is the Kronecker delta. Strictly speaking, the determinant in Eq. (10) represents the growth of an -dimensional rather than -dimensional volume element, because it includes the energy direction itself. However, since this is only one of directions, the error in the final result associated with the above approximation is of the order of . Now, we write explicitly

(11) |

and then observe that for a system with short-range interactions, the first term in the above equation is of the order of , while the second term is of the order of and hence can be neglected. We further notice that the leading (first-order) contributions to the determinant in Eq. (10) in terms of come only from the diagonal elements of the matrix. Taking into account the above two considerations, we finally obtain that, in the limit ,

(12) |

where . The total change of the volume of the energy shell is then

(13) |

where the notation implies the average of the entire energy shell.

Eq. (13) together with Eq. (2) implies that . Since both and contain additive small contributions associated with the uncorrelated remote parts of a large system, the distributions for each of them are narrowly peaked around the respective average values (according to the central limit theorem). Therefore, in the limit ,

(14) |

In our simulations, both and are obtained from the perturbations of the phase space vector associated with the small random spin rotations introduced in the main part of the paper and characterized by mean-squared values . These perturbations are illustrated schematically in Fig. 6. The energy change for each perturbation is

(15) |

from which it follows that, in the limit , , while . Substituting the latter two formulas into Eq. (14), we obtain Eq. (3).

## Appendix B Finite-size effects

Here we include four figures illustrating the dependence of several quantities computed for the Heisenberg model on lattice sizes: Fig. 7 shows the size dependence of the specific heat, Figs. 8 and 9 show the largest Lyapunov exponent and the Kolmogorov-Sinai entropy as functions of energy and temperature, respectively, and, finally, Fig. 10 shows the -index.

## Appendix C Simultaneous sign reversal of the total energy and one of the coupling constants

Figure 11 shows plots similar to Figs. 4 (c) and (d), but for a different sign of the couplings. The comparison of Fig. 11 (a) with Fig. 4 (d), and Fig. 11 (b) with Fig. 4 (c) demonstrates the symmetry of Lyapunov spectra of bipartite spin lattices with respect to the simultaneous sign reversals of the total energy and of one of the coupling constants.

### References

- P. Butera and G. Caravati, Phys. Rev. A 36, 962 (1987).
- V. Mehra and R. Ramaswamy, Phys. Rev. E 56, 2508 (1997).
- K.-H. Kwon and B.-Y. Park, J. Chem. Phys. 107, 5171 (1997).
- L. Caiani, L. Casetti, C. Clementi, and M. Pettini, Phys. Rev. Lett. 79, 4361 (1997).
- L. Caiani, L. Casetti, C. Clementi, G. Pettini, M. Pettini, and R. Gatto, Phys. Rev. E 57, 3886 (1998a).
- L. Caiani, L. Casetti, and M. Pettini, J. Phys. A: Math. Gen. 31, 3357 (1998b).
- M.-C. Firpo, Phys. Rev. E 57, 6599 (1998).
- V. Latora, A. Rapisarda, and S. Ruffo, Phys. Rev. Lett. 80, 692 (1998).
- L. Casetti, M. Cerruti-Sola, M. Modugno, G. Pettini, M. Pettini, and R. Gatto, Rivista del Nuovo Cimento 22 1–74 (1999).
- J. Barre and T. Dauxois, Europhys. Lett. 55, 164 (2001).
- M. Cerruti-Sola, C. Clementi, and M. Pettini, Phys. Rev. E 61, 5171 (2000).
- L. Casetti, M. Pettini, and E. Cohen, Phys. Rep. 337, 237 (2000).
- H. A. Posch, W. G. Hoover, and B. L. Holian, Ber. d. Bunsenges. f. Phys. Chem. 94, 250 (1990).
- C. Dellago and H. A. Posch, Physica A 230, 364 (1996).
- C. Dellago and H. Posch, Physica A 237, 95 (1997).
- H. Bosetti and H. A. Posch, Communications in Theoretical Physics 62, 451 (2014).
- A. S. de Wijn, B. Hess, and B. V. Fine, Phys. Rev. Lett. 109, 034101 (2012).
- A. S. de Wijn, B. Hess, and B. V. Fine, Journal of Physics A: Mathematical and Theoretical 46, 254012 (2013).
- G. Benettin, L. Galgani, A. Giorgilli, and J. M. Strelcyn, Meccanica 15, 9 (1980).
- B. Hess, Investigations of Lyapunov spectra for classical spin lattices, Diploma Thesis, University of Heidelberg (2010), www.thphys.uni-heidelberg.de/fine/Hess_Diploma_2010.pdf.
- Y. Deng, H. W. J. Blöte, , and M. P. Nightingale, Phys. Rev. E 72, 016128 (2005).
- R. G. Brown and M. Ciftan, Phys. Rev. B 74, 224413 (2006).
- B. V. Fine, Phys. Rev. Lett. 79, 4673 (1997).
- B. V. Fine, Int. J. Mod. Phys. B 18, 1119 (2004).
- H. H. Rugh, Phys. Rev. lett. 78, 772 (1997).