{'fraction_non_alphanumeric': 0.04686448912474774, 'fraction_numerical': 0.015247776365946633, 'mean_word_length': 4.24294670846395, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Largest Lyapunov exponents for lattices of interacting classical spins

We investigate how generic the onset of chaos in interacting many-body classical systems is in the context of lattices of classical spins with nearest neighbor anisotropic couplings. Seven large lattices in different spatial dimensions were considered. For each lattice, more than 2000 largest Lyapunov exponents for randomly sampled Hamiltonians were numerically computed. Our results strongly suggest the absence of integrable nearest-neighbor Hamiltonians for the infinite lattices except for the trivial Ising case. In the vicinity of the Ising case, the largest Lyapunov exponents exhibit a power-law growth, while further away they become rather weakly sensitive to the Hamiltonian anisotropy. We also provide an analytical derivation of these results.

The concept of microscopic chaos in many-particle systems plays an essential role in the foundations of statistical physics [1][2][3]. In classical systems, chaos is defined by the appearance of exponential instabilities with respect to infinitesimal perturbations of the initial conditions. The spectrum of these instabilities is characterized by a set of eigenvalues -Lyapunov exponentsand the corresponding eigenvectors. In general, interacting, many-particle, classical systems are expected to be chaotic. The largest Lyapunov exponents and the whole Lyapunov spectra have been calculated numerically for classical many-body systems, such as gases of hard-core particles [4,5], fluids with soft interactions [6], and lattice two-dimensional rotators [6,7], and analytically in a few cases [8][9][10][11][12][13][14]. Nevertheless, little is known about the universality of the features of the Lyapunov spectra, in particular when it comes to the systems with smooth dynamics. It is also not clear how generic the onset of chaos actually is, and what happens in many-particle systems in the vicinity of integrable, i.e. nonchaotic, limits of the microscopic Hamiltonians. This Letter deals with the above issues on the basis of a systematic numerical and analytical investigation of the systems of interacting classical spins. These systems have been extensively studied, e.g., in the context of the spin diffusion problem [15][16][17], but their Lyapunov spectra have not yet been computed. Finite classical spin systems are also known to exhibit nontrivial integrable limits [18]. Another aspect of our motivation is the connection to the quantum case. There is a growing appreciation that generic interacting many-particle quantum systems exhibit relaxation behavior similar to their classical chaotic counterparts, in particular as far as nuclear spin decays in solids are concerned [19][20][21][22][23][24][25][26][27]. In a broader context, the important conceptual development in the present study is the massive character of the numerical investigation covering an entire class of Hamiltonians.
In this Letter, we consider several large lattices of interacting classical spins. For each lattice, we compute the largest Lyapunov exponent λ max for several thousand randomly selected microscopic Hamiltonians. Thereby, we obtain the dependence of the largest Lyapunov exponent on the anisotropy of the spin-spin interaction. Our results give a strong indication of the absence of integrable Hamiltonians for infinite spin lattices with nearestneighbor interaction besides the trivial Ising case (explained below). We also find that, to the extent afforded by our numerical accuracy, the system becomes chaotic in the immediate vicinity of the integrable Ising case, with λ max exhibiting a universal power-law scaling. Further away from the Ising limit, λ max becomes only weakly dependent on the Hamiltonian anisotropy, especially for bipartite lattices. The letter is concluded with a simple analytical derivation that describes the above-mentioned numerical results.
We investigate the seven lattices shown in Fig. 1: (L1) a chain, (L2) a rectangular ladder, (L3) a square lattice, (L4) a bilayer of square lattices, (L5) a cubic lattice, (L6) a triangular ladder and (L7) a triangular lattice. The interaction Hamiltonian for each lattice is of the nearestneighbor (NN) type with periodic boundary conditions: where (S ix , S iy , S iz ) ≡ S i are the three projections of the ith classical spin normalized by the condition S 2 i = 1. We numerically integrate the equations of motion associated with the Hamiltonian (1):Ṡ i = S i × h i [28], where h i is the local field given by the expression Here e x , e y and e z are the unit vectors along the respective directions, and j(i) implies the summation over the nearest neighbors of the i-th lattice site. We use the fourth-order Runge-Kutta algorithm and choose a time step of 0.005, sufficiently small so that, on the time scales of our simulations, energy is conserved. The initial conditions of each trajectory are chosen randomly on the energy shell with zero total energy, which corresponds to infinite temperature [29].
In order to obtain λ max , we numerically calculate a phase space trajectory and, at every time step, track the evolution of a tangent space vector that defines an infinitesimal perturbation of it [22,30]. The asymptotic average growth rate of this vector is equal to the largest Lyapunov exponent. In all cases, λ max becomes sizeindependent for sufficiently large lattices.
For each lattice, we computed λ max for many combinations of the coupling constants randomly selected on the "interaction sphere" J 2 x + J 2 y + J 2 z = 1. Specifically, 8000 combinations selected from the isotropic distribution were selected for the lattices (L1), (L3), (L5), and 2000 combinations for each of the remaining lattices. In addition, in order to investigate the scaling behavior near the Ising limit, we have processed more combinations in the vicinity of J z = 1: 1000 for the lattices (L1), (L3), and (L5) and 250 for the others.
Our main numerical findings, namely λ max as a function of the parameter J max ≡ max(|J x |, |J y |, |J z |), are presented in Fig. 2. The maximum value J max = 1 is realizable only in the Ising case, and thus represents the integrable limit with λ max = 0. The minimum value of J max is 1/ √ 3. It corresponds to either Heisenberg case J x = J y = J z or "anti-Heisenberg" case J x = J y = −J z (or equivalent cases). In Fig. 3, we present the rescaled plots λ max (J max )/λ max (1/   (i) For all lattices considered, no integrable cases besides the Ising case J max = 1 were found. This indicates that, for infinite lattices of classical spins with the nearest-neighbor interaction, the existence of another integrable case is highly unlikely.
(ii) The value of λ max is mainly controlled by J max , especially for the bipartite lattices.
(iii) The dependence λ max (J max ) for bipartite lattices (L1-L5) has nearly universal form, as illustrated by the rescaling shown in Fig. 3.
(iv) The above dependence is nearly flat below J max ≈ 0.85; i.e., away from the integrable limit λ max is very weakly sensitive to the details of microscopic interaction.
(vi) As can be seen from Fig. 2(a), the nonbipartite lattices (L6) and (L7) show a fork-shaped spread of λ max as J max approaches 1/ √ 3. The upper and the lower tips of the fork correspond to the anti-Heisenberg and Heisenberg cases, respectively.
Gross features of the above results can be reproduced by rather simple analytical estimates. In particular, away from the Ising case, the plateau values of λ max seen in Fig. 2(a) at J max < 0.85 can be estimated with a factor-of-two accuracy as the typical frequency of onespin motion given by the root-mean-squared value of the local field h i at the infinite temperature: λ max ∼ [n(J 2 x + J 2 y + J 2 z )/3] 1/2 , where n is the number of the nearest neighbors. This expression also predicts that the heights of these plateaus scale as n 1/2 , while the inset of Fig. 3 indicates that the actual scaling is close but somewhat steeper. Now we turn to the approximation for the dependence λ max (J max ) as the system approaches the Ising limit. We assume that J z ≫ J x , J y , which implies that J max = J z . We also introduce variable J ⊥ ≡ [(J 2 x + J 2 y )/2] 1/2 = [(1 − J 2 max )/2] 1/2 to denote the typical value of the transverse coupling. We consider two phase space trajectories In order to estimate λ max , we linearize the equations of motion (2) with respect to small δS i (t) and keep the leading order in terms of J ⊥ /J z : where δS iz and δS iϕ are the projection of vectors δS i on the directions of vectors e z and e z × S i (t), respectively. The third projection of δS i does not appear in Eqs. (3) and (4), because it can be expressed in terms of δS iz -a consequence of the constraint S 2 i = 1. The parameters C ij z (t) and C ij ϕ (t) are determined by S i (t) and S j (t) and have characteristic fluctuation times 1/J ⊥ and 1/J z , respectively. Now, we make an assumption justified by the final result [Eq. (12)] that J ⊥ ≪ λ max ≪ J z . We estimate the growth of the typical values of δS iz and δS iϕ over time τ such that λ max ≪ 1/τ ≪ J z . On the timescale τ , the parameters C ij z (t) stay nearly constant, while C ij ϕ (t) strongly fluctuate, so that C ij ϕ (t) τ ≈ 0. We first write In this problem, a relatively slow growth of S iϕ is coupled to a random-walk-like growth of S iz . In order to extract λ max , we, therefore, look at the leading terms in the growth of δS 2 iϕ and δS 2 iz : In Eq.(8), we neglected the term linear in J ⊥ , because C ij ϕ (t) τ ≈ 0.
Since we are only concerned with the scaling of λ max with J max and n, we convert Eqs. (7) and (8) into an order-of-magnitude estimate for the typical growth of δS 2 iϕ and δS 2 iz . The estimate includes (i) dropping lattice index in Eqs. (7) and (8), (ii) estimating the instantaneous values of parameters C ij z (t) and C ij ϕ (t) by 1, (iii) replacing the sum in Eq. (7) by √ n (a consequence of the random sign of the n terms in that sum), and (iv) approximating the integral term in Eq.
n/J z . The factor √ n/J z in the latter estimate is due to n terms in the sum, each producing a contribution to the integral of the order 1/h iz ∼ 1/(J z √ n). Finally, after dividing thus simplified Eqs. (7) and (8) by (δS ϕ ) 2 and (δS z ) 2 , respectively, we obtain Since the parameters δS 2 ϕ and δS 2 z should grow at the same rate, this implies that Substituting Eq. (11) into Eq. (10) and comparing the right-hand side with the expression 1+2λ max τ , we finally obtain the estimate: which is consistent with the numerically observed 1/3power law shown in Fig. 2(b). In Ref. [14] the same power law was obtained for weakly interacting dilute gases via a perturbation expansion around the integrable ideal gas. The prefactor scaling as √ n is also compatible with the numerical results, as illustrated in the inset of Fig. 2(b).
Equation (11) predicts further that the components of the Lyapunov eigenvector corresponding to λ max depend systematically on whether they are parallel to the z-direction or not. As illustrated in Fig. 4, this prediction also agrees with our numerical results.
Our derivation of Eq. (12) is based on the assumption of the random fluctuations of the sum in Eq. (6). However, in the Ising limit, the Fourier transform of this sum contains only a finite number of frequencies. Therefore, in the vicinity of the Ising limit, recurrences may occur if correlations in S iz do not decay sufficiently fast. Such recurrences, in turn, would contradict our assumption of the fast decay of the time correlations of the above sum. This would be most problematic for lattices with smaller numbers of nearest neighbors, such as spin chains. Indeed, we observe in Fig. 2(b), that spin chains exhibit larger deviations from the 1/3-power law than other lattices. The above recurrences may, in fact, be the route to breaking down the chaotic nature of the spin dynamics around the Ising limit. However, our numerical results indicate that, if such a breakdown occurs, it occurs at values of |J ⊥ /J z | < 10 −4 .
Finally, we mention that the above estimate may be repeated for the case of J x and J y smaller but not much smaller than J z . In this a case, τ must be chosen much shorter than 1/J ⊥ and thus both C ij z (t) and C ij ϕ (t) can be assumed constant on the time scale of τ . This estimate would then give λ max ∼ = (nJ max ) 1/2 (1 − J 2 max ) 1/4 , which, as illustrated in Fig. 3, exhibits a good overall agreement with the numerical results for bipartite lattices over the entire range of J max .
In conclusion, we have presented a systematic study of largest Lyapunov exponents λ max for a very large variety of classical spin systems. Our findings strongly suggest the absence of integrable nearest-neighbor Hamiltonians for the type of lattices considered except for the Ising case. As far as the behavior of λ max is concerned, a number of universal features enumerated above as (ii)-(vi) are observed. We have also analytically derived the scaling of λ max with the anisotropy of the Hamiltonian.
The numerical part of this work was performed at the bwGRiD computing cluster at the University of Heidelberg. The authors are grateful to T. A. Elsayed and G. P. Morriss for discussions. A.S.W.'s work is financially supported by an Unga Forskare grant from the Swedish Research Council.