d-Dimensional Classical Heisenberg Model with Arbitrarily-Ranged Interactions: Lyapunov Exponents and Distributions of Momenta and Energies

Abstract

We numerically study the first-principle dynamics and thermostatistics of a d-dimensional classical inertial Heisenberg ferromagnetic model ($d=1,2,3$) with interactions decaying with the distance $r_{ij}$ as $1/r_{ij}^{\alpha }$ ($\alpha \ge 0$), where the limit $\alpha =0$ ($\alpha \to \infty$) corresponds to infinite-range (nearest-neighbour) interactions, and the ratio $\alpha /d>1$ ($0\le \alpha /d\le 1$) characterizes the short-ranged (long-ranged) regime. By means of first-principle molecular dynamics we study: (i) The scaling with the system size $\mathcal{N}$ of the maximum Lyapunov exponent $\lambda$ in the form $\lambda \sim {\mathcal{N}}^{-\kappa }$, where $\kappa (\alpha /d)$ depends only on the ratio $\alpha /d$; (ii) The time-averaged single-particle angular momenta probability distributions for a typical case in the long-range regime $0\le \alpha /d\le 1$ (which turns out to be well fitted by q-Gaussians), and (iii) The time-averaged single-particle energies probability distributions for a typical case in the long-range regime $0\le \alpha /d\le 1$ (which turns out to be well fitted by q-exponentials). Through the Lyapunov exponents we observe an intriguing, and possibly size-dependent, persistence of the non-Boltzmannian behavior even in the $\alpha /d>1$ regime. The universality that we observe for the probability distributions with regard to the ratio $\alpha /d$ makes this model similar to the $\alpha$-XY and $\alpha$-Fermi-Pasta-Ulam Hamiltonian models as well as to asymptotically scale-invariant growing networks.

1. Introduction

Although nature is strictly quantum at the microscopic level, classical systems of particles have been studied for decades, usually considered as good approximations in some limiting cases. In what concerns magnetism, classical spin models have succeeded in providing good qualitative agreement with experiments on several materials [1,2,3]. Many of these investigations focused on equilibrium properties of models characterized by nearest-neighbour interactions (see, e.g., Refs. [4,5,6,7,8,9,10,11]). In contrast, a paradigmatic example dealing with infinite-range interactions, the so-called Hamiltonian Mean Field (HMF) model, consisting of N fully coupled classical XY rotators, was introduced in [12]. This was done through the addition of a kinetic term in the Hamiltonian of a classical XY model, so that equations of motion could be derived, leading to the possibility of investigating its dynamic properties by a direct integration of these equations. Since then, a vast amount of papers devoted to the study of the HMF model have put forward its anomalous properties [13,14,15,16,17,18,19,20,21,22,23]: existence of negative specific heat, quasistationary states (QSSs), nonmaxwellian velocity probability distributions, anomalous diffusion and ergodicity breaking, among others.

The HMF model was further generalized into the so-called $\alpha$-XY model [24], via the introduction of an interaction term decaying with the distance r between rotators as $r^{-\alpha }$. The parameter $\alpha \ge 0$ controls the range of the interaction, with the limits $\alpha =0$ and $\alpha \to \infty$ corresponding respectively to the HMF and nearest-neighbour models. For the $\alpha$-XY model defined on d-dimensional lattices, the mean potential energy per particle diverges in the thermodynamic limit for $0\le \alpha /d\le 1$. Consequently, serious mathematical difficulties emerge in the thermodynamic limit for the Boltzmann-Gibbs (BG) partition function. Nowadays, it has become evident that BG statistical mechanics fails in describing appropriately systems generically characterized by relevant long-range correlations (e.g., due to long-range interactions). In such cases, the classification of thermodynamical quantities as either extensive or intensive is not sufficient, since one frequently finds a third category, the one of nonextensive quantities.

These kinds of difficulties have led to generalizations of BG statistical mechanics, and the most widely used theory for dealing with such systems so far is nonextensive statistical mechanics [25]. This theory emerged from the proposal of a generalized entropy, known as $S_q$, characterized by a real index q [26],

$$S_q=k\sum _{i=1}^W{p}_i\left({ln}_q\frac{1}{p_i}\right),$$

(1)

where we have introduced the q-logarithm definition,

$${ln}_{q}u=\frac{u^{1-q}-1}{1-q};\left({ln}_{1}u=lnu\right).$$

(2)

Therefore, one recovers BG entropy as ${lim}_{q\to 1}{S}_q=S_{\mathrm{BG}}$, whereas in the microcanonical ensemble, where all microstates present equal probability, $p_i=1/W$, Equation (1) becomes,

$$S_q=k{ln}_{q}W.$$

(3)

Since then, a considerable progress was achieved, leading to generalized functions, distributions, and important equations of physics. In particular, the q-Gaussian distribution, which generalizes the standard Gaussian, appears naturally by extremizing the entropy $S_q$ [25],

$$P_q\left(u\right)=P_0{exp}_q(-\beta u^2),$$

(4)

with

$${exp}_q\left(u\right)={[1+(1-q)u]}_{+}^{1/(1-q)};\left({exp}_1\left(u\right)=exp\left(u\right)\right),$$

(5)

where $P_0\equiv P_q\left(0\right)$ and ${\left[y\right]}_{+}=y$, for $y>0$ (zero otherwise). Above, the q-exponential function appears precisely as the inverse function of the q-logarithm of Equation (2), i.e., ${exp}_q\left({ln}_{q}u\right)={ln}_q\left({exp}_q\left(u\right)\right)=u$.

The $\alpha$-XY model has also been largely studied [27,28,29,30,31,32,33,34] and, in the d-dimensional case, it was shown [29] that in its high-energy regime, the maximum Lyapunov exponent scales with the number $\mathcal{N}\equiv N^d$ of rotators as $\lambda \sim {\mathcal{N}}^{-\kappa }$, where the exponent $\kappa (\alpha ,d)$ depends on $\alpha$ and d through the ratio $\alpha /d$, as in the case of other quantities, with $\kappa >0$ in the long-range regime ($\alpha /d\le 1$) and $\kappa =0$ (thus yielding positive Lyapunov exponents) in the short-range regime ($\alpha /d>1$). In addition, q-Gaussians have been obtained for the time-averaged momenta probability distributions not only before, but even after the crossover to the very-large-time QSS, with a value of the entropic index $q\equiv q_p(\alpha /d)$ depending again on the ratio $\alpha /d$ [34]. The usual Maxwellian distributions are recovered as $\alpha /d$ increases above unity. A curious behavior occurs also for the ensemble-averaged momenta probability distributions. The same departure from the BG predictions is observed for the time-averaged energy distributions; in this case, instead of the usual BG exponential, q-exponential probability distributions are obtained, with $q\equiv q_E(\alpha /d)$ [34]. Similar universality features have been recently obtained in quite different contexts, such as a generalized Fermi-Pasta-Ulam model with long-range interactions [35,36,37,38,39], complex networks with preferential attachment growth [40,41], and a system of particles under overdamped motion interacting repulsively with power-law interactions [42].

By considering three-dimensional rotators (instead of planar ones), the corresponding classical inertial ferromagnetic Heisenberg model is obtained; this has not been very much explored in the literature [43,44,45,46,47]. In the fully coupled (infinite-range interactions) version, QSSs have also been observed, implying that, just as in the HMF model, the thermodynamic and time limits do not commute: if the $N\to \infty$ limit is taken first, the system will remain in the QSS, whereas if the $t\to \infty$ limit is considered first, the system may attain the Boltzmann temperature after a transient time, which was shown to increase with the system size [44,45,46,47]. Inspired by the $\alpha$-XY model, the classical inertial ferromagnetic Heisenberg model was also generalized by introducing an interaction term decaying with the distance r between rotators as $r^{-\alpha }$ ($\alpha \ge 0$) [47]. Such an $\alpha$-Heisenberg model, studied up to the moment only in dimension $d=1$, exhibited a clear QSS in the long-range-interaction regime, with its duration decreasing for increasing values of $\alpha$. In this paper, we extend the $\alpha$-Heisenberg model to dimensions $d=2$ and 3, and show that the scaling with the system size of the maximum Lyapunov exponent follows the same trend as in the above indicated $\alpha$-XY model. In addition, we calculate numerically the probability distributions for the time-averaged momenta, as well as the energy distributions, in order to get $\alpha /d$-dependent entropic indexes. The anomalous (non-Maxwellian for the momenta and nonexponential for the energies) distributions persist well beyond the $\alpha =d$ threshold and intriguingly coexist with positive values of the maximal Lyapunov exponent, at least for the sizes that have been numerically attained.

2. The Model

The $\alpha$-Heisenberg inertial model consists of a collection of $\mathcal{N}$ interacting rotators, represented by three-component rotating vectors ${\overrightarrow{S}}_i$ with unit norm, located at fixed positions on a d-dimensional hypercubic lattice of linear size N ($\mathcal{N}\equiv N^d$).The Hamiltonian is given by

$$\mathcal{H}=K+V_{\alpha }=\frac{1}{2}\sum _{i=1}^{\mathcal{N}}{\overrightarrow{L}}_i^2+\frac{1}{2\tilde{\mathcal{N}}}\sum _{i=1}^{\mathcal{N}}\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^{\mathcal{N}}\frac{1-{\overrightarrow{S}}_i·{\overrightarrow{S}}_j}{r_{ij}^{\alpha }},$$

(6)

where $\alpha \ge 0$ controls the interaction range, $r_{ij}$ denotes the distance between rotators i and j, the angular momenta ${\overrightarrow{L}}_i$ coincide with the angular velocities since we are considering unit moments of inertia. Moreover, the prefactor with $\tilde{\mathcal{N}}(\alpha ,d)\equiv {\sum }_{j\ne i}\frac{1}{r_{ij}^{\alpha }}$ in the potential energy is convenient in order to get a total potential energy $V_{\alpha }$ which scales linearly with $\mathcal{N}$ [25], like the total kinetic energy.

2.1. The $\alpha =0$ and $\alpha \to \infty$ Limiting Cases

It is easily seen that, for $\alpha =0$, $\tilde{\mathcal{N}}(0,d)=\mathcal{N}-1\sim \mathcal{N}$, so in the infinite-range case the total potential energy reduces to

$$V_0=\frac{1}{2\mathcal{N}}\sum _{i=1}^{\mathcal{N}}\sum _{j=1}^{\mathcal{N}}(1-{\overrightarrow{S}}_i·{\overrightarrow{S}}_j)=\frac{\mathcal{N}}{2}-\frac{1}{2\mathcal{N}}\sum _{i=1}^{\mathcal{N}}{\overrightarrow{S}}_i\sum _{j=1}^{\mathcal{N}}{\overrightarrow{S}}_j=\frac{\mathcal{N}}{2}(1-M^2),$$

(7)

where the restriction $j\ne i$ has been trivially dismissed since ${\overrightarrow{S}}_i^2=1$, and the order parameter M is the norm of the total magnetization vector

$$\overrightarrow{M}=\frac{1}{\mathcal{N}}\sum _{i=1}^{\mathcal{N}}{\overrightarrow{S}}_i.$$

(8)

In addition, following the usual prescription for the kinetic temperature, $T=\frac{2}{n\mathcal{N}}\langle K\rangle =\frac{1}{\mathcal{N}}\langle K\rangle$, where $n=2$ represents the number of degrees of freedom per spin (three components of the spin vector minus the constant norm restriction). Thus, considering the energy per particle $U\equiv \frac{\langle \mathcal{H}\rangle }{\mathcal{N}}$ in Hamiltonian (6) with $\alpha =0$ together with Equation (7), one gets

$$U=T+\frac{1}{2}(1-M^2),$$

(9)

which, together with the consistency equation for the magnetization [1,2,3],

$$M=\mathcal{L}(M/T),$$

(10)

with the Langevin function $\mathcal{L}\left(x\right)\equiv cothx-\frac{1}{x}$, allows us to determine a relation between the thermodynamical variables of the model, which undergoes a second-order phase transition for the critical values $(T_c,U_c)=(\frac{1}{3},\frac{5}{6})$. The caloric curve, as well as the $M-T$ curve and some other details of the model can be found in Ref. [47]. It can be shown that [27], within BG statistical mechanics, Equations (9) and (10) are not only valid for $\alpha =0$ but also throughout the long-range regime $0\le \alpha /d<1$.

No analytical results are available in the short-range regime except for the limiting case $\alpha \to \infty$, which corresponds to nearest-neighbors interactions, where $\tilde{\mathcal{N}}(\infty ,d)=2d$ (that is, the number of nearest-neighbour spins on a d-dimensional hypercubic lattice) and the potential energy reduces to

$$V_{\infty }=\frac{1}{2d}\sum _{\langle i,j\rangle }(1-{\overrightarrow{S}}_i·{\overrightarrow{S}}_j),$$

(11)

where $\langle i,j\rangle$ indicates that sites labelled i and j are nearest neighbours. The equation for the internal energy in the one-dimensional case reads

$$U=T+\frac{1}{2}(1-\mathcal{L}(1/2T)),$$

(12)

while $M=0,\forall T\ne 0$ (there is no phase transition in the $d=1$ short-range regime).

2.2. Equations of Motion

From the definition ${\overrightarrow{L}}_i={\overrightarrow{S}}_i\times {\dot{\overrightarrow{S}}}_i$ one trivially gets ${\dot{\overrightarrow{S}}}_i={\overrightarrow{L}}_i\times {\overrightarrow{S}}_i$. In addition, ${\dot{\overrightarrow{L}}}_i={\overrightarrow{S}}_i\times {\ddot{\overrightarrow{S}}}_i=-{\overrightarrow{S}}_i\times \frac{\partial V_{\alpha }}{\partial {\overrightarrow{S}}_i}$, so the equations of motion are given by

$$\begin{array}{cc}\hfill {\dot{\overrightarrow{S}}}_i& ={\overrightarrow{L}}_i\times {\overrightarrow{S}}_i\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\dot{\overrightarrow{L}}}_i& ={\overrightarrow{S}}_i\times {\overrightarrow{M}}_i\hfill \end{array}$$

(14)

where

$${\overrightarrow{M}}_i\equiv -\frac{\partial V_{\alpha }}{\partial {\overrightarrow{S}}_i}=\frac{1}{\tilde{\mathcal{N}}}\sum _{j\ne i}^{\mathcal{N}}\frac{{\overrightarrow{S}}_j}{r_{ij}^{\alpha }},$$

(15)

carries the contribution of the interactions of spin i with all $(\mathcal{N}-1)$ other spins. This quantity becomes the total magnetization vector of (8) when $\alpha =0$, and the sum $\frac{1}{2d}{\sum }_{\langle i,j\rangle }{\overrightarrow{S}}_j$, extended to the nearest neighbours of spin i, when $\alpha =\infty$. From Equations (13) and (14) not only follows the aforementioned conservation of the spin norms, since ${\dot{\overrightarrow{S}}}_i^2=2{\overrightarrow{S}}_i·{\overrightarrow{S}}_i=2{\overrightarrow{S}}_i·({\overrightarrow{L}}_i\times {\overrightarrow{S}}_i)=0$, but also the conservation of the total angular momentum $\overrightarrow{L}={\sum }_{i=1}^{\mathcal{N}}{\overrightarrow{L}}_i$, since $\dot{\overrightarrow{L}}={\sum }_{i=1}^{\mathcal{N}}{\dot{\overrightarrow{L}}}_i=\frac{1}{\tilde{\mathcal{N}}}{\sum }_{i=1}^{\mathcal{N}}{\sum }_{j\ne i}\frac{{\overrightarrow{S}}_i\times {\overrightarrow{S}}_j}{r_{ij}^{\alpha }}=\overrightarrow{0}$.

3. Results

We have numerically integrated the equations of motion (13) and (14) by means of the standard fourth-order Runge-Kutta algorithm with an integration step $h=0.02$. It is well known that numerical methods, like the one used herein, present an energy drift for long-time simulations [48,49,50]; using the integration step above, we have checked that the total energy was preserved within a relative precision ${10}^{-5}$.We have considered two different kinds of initial conditions. In the first one we take random sets of angular momenta $\left\{{\overrightarrow{L}}_i\right\}$, with components drawn from a symmetric uniform distribution and then redefined in order to yield zero total momentum. As for the spins $\left\{{\overrightarrow{S}}_i\right\}$, two of their components are randomly chosen, while the third one is obtained by the orthogonality restriction ${\overrightarrow{S}}_i·{\overrightarrow{L}}_i=0$$(\forall i)$. Then, we get spins of unit length dividing by the norm. Finally, the angular momenta are rescaled so as to obtain the desired total energy per particle U which remains virtually constant throughout our microcanonical ensemble simulations. This process yields a poorly magnetized initial state with $\overrightarrow{M}\left(0\right)\sim \overrightarrow{0}$ and $V_{\alpha }\left(0\right)\sim 0.5$. In the second kind of initial conditions, we take all the spins equal and normalized, so as to get a fully magnetized state with $M\left(0\right)=1$ and zero potential energy. Then we randomly choose the angular momenta in the perpendicular plane to finally translate and rescale them to get zero total angular momentum and the desired total energy. Large time-limit qualitative and quantitative results coincide with both kinds of initial conditions.

3.1. Time-Averaged Momenta and Energy One-Particle Distributions

As already shown in Ref. [47], for values of the total energy below the critical one but close enough to the transition, after a short initial transient, the kinetic temperature of the system, $T\left(t\right)=\frac{1}{\mathcal{N}}\langle K\left(t\right)\rangle$, decreases to a first plateau value given by $T=U-1/2$ (consider $M=0$ in Equation (9)), corresponding to a nonmagnetized state. Only after waiting long enough, the system abandons the QSS when the temperature as well as the magnetization rise to reach the values predicted by Equation (9) and stays forever in this second plateau, which is non Boltzmannian though its temperature coincides with that associated with BG statistics. It is in this (nearly) final QSS where we choose a time window $[t_0,t_f]$ of width $\Delta t=W\times h$ and calculate the averaged-in-time momenta as $\langle {\overrightarrow{L}}_i\rangle =\frac{1}{W+1}{\sum }_{k=0}^W{\overrightarrow{L}}_i(t_0+kh)$. Then we put together all the $3\mathcal{N}$ averaged-in-time angular momenta components to build the corresponding histograms. In Figure 1 we present the probability distribution functions of time-averaged angular momenta $P\left(\langle L_i\rangle \right)$ for one-, two-, and three-dimensional systems. These distributions were computed for systems with $\overrightarrow{M}\left(0\right)\sim 0$ initial conditions, total energy per particle $U=0.76<U_c=5/6$, size $\mathcal{N}=262144={\left(512\right)}^2={\left(64\right)}^3$ (corresponding to $d=1,2,$ and 3), and $\alpha /d=0.9$, that is, in the long-range-interaction regime. A clear collapse of all three distributions into a q-Gaussian (cf. Equation (4)) is exhibited, with $q_L(\alpha /d)=1.42$ and ${\beta }_L=4.44$; the maximum values $P_0\equiv P(\langle L_i\rangle =0)$ fluctuated slightly with the lattice dimensionality, being found as $P_0=0.27(d=1)$, $P_0=0.19(d=2)$, and $P_0=0.22(d=3)$.Further results exhibiting values of $q_L$ as a function of the ratio $\alpha /d$ will be published elsewhere.

In Figure 2 we exhibit the probability distribution functions of time-averaged angular momenta $P\left(\langle L_i\rangle \right)$ for $\alpha /d=1.6$, that is, in the short-range regime, with the same initial conditions and parameters used in Figure 1. Now, the collapse into a single curve occurs only around the origin; next we point out possible reasons for the vital differences between the probability distributions shown in Figure 1 and Figure 2. (i) The distinction in the time evolution of the kinetic temperature, shown in the left insets in both figures. Particularly, one should emphasize that the equilibrium temperature (along which time averages $\langle L_i\rangle$ were computed) in Figure 2 depends on $\alpha$ and d separately. (ii) Finite-size effects, which become more significant in the limit of short-range interactions, might be playing an important role. (iii) At the BG equilibrium, the corresponding short-range-interaction models present further relevant distinctions, e.g., in $d=1$ there is no phase transition, whereas an ordered state exists in $d=3$ [1,2,3]. (iv) The maximum Lyapunov exponents show significant differences in the cases investigated in Figure 1 and Figure 2, as we shall see below.

We have also calculated the probability distributions of time-averaged individual energies, i.e., $\langle E_i\rangle =\frac{1}{W+1}{\sum }_{k=0}^W{E}_i(t_0+kh)$, with $E_i=\frac{1}{2}{\overrightarrow{L}}_i^2+\frac{1}{2\tilde{\mathcal{N}}}{\sum }_{j\ne i}^{\mathcal{N}}\frac{1-{\overrightarrow{S}}_i·{\overrightarrow{S}}_j}{r_{ij}^{\alpha }}$. Results for the long-range-interaction regime, with $\alpha /d=0.9$ and the same initial conditions and parameter values considered in Figure 1, are shown in Figure 3. Again, all data collapse, now in the energy range $E_i\ge \mu$ (see details in Figure 3); since different dimensions present distinct densities of states, the data does not collapse near the origin. One should call attention to the fact that the value $q_E=1.3$ of Figure 3 coincides with the recent result for the $\alpha$-XY model, obtained also with $\alpha /d=0.9$ [34]. Moreover, as found in the study of the $\alpha$-XY model [34], a curious discrepancy in the values of the index q occurred in the estimates of Figure 1 and Figure 3, i.e., $q_L>q_E$. Whether this may be a real effect, or a consequence of finite-size effects, remains a point for further investigation. Moreover, similarly to Figure 2, the corresponding probability distributions $P\left(\langle E_i\rangle \right)$ in the short-range regime (for $\alpha /d=1.6$), do not reveal a collapse of the curves into a q-exponential, so that the significant distinctions in the short-range regime persist.

3.2. Nearest-Neighbour Limit

In order to shed more light on the short-range regime, we shall focus now on the limiting $\alpha \to \infty$ case. We have calculated the equilibrium temperature, as well as the magnetization as a function of the total energy per particle and the dimension. Analytical results are available only for the caloric curve in the one-dimensional case, as given by Equation (12) and shown in Figure 4, where we obtain an excellent agreement between the theoretical results and the numerical ones obtained in our microcanonical ensemble simulations with $M\left(0\right)=1$ initial conditions (which are mandatory in this case since, for small values of the constant total energy per particle we need to have vanishing initial potential energy and the potential energy for the $\overrightarrow{M}\left(0\right)\sim 0$ initial conditions is $V\sim 0.5$) and $\mathcal{N}=$ 46,656. The caloric curve corresponding to the long-range regime $0\le \alpha /d<1$ is also shown in Figure 4 where it can be seen that the temperature coincides in both limits $\alpha \to 0$ and $\alpha \to \infty$ in all dimensions within the energy range $0\le U\le 0.25$. For larger values of U, the equilibrium temperature decreases with the dimension in the nearest neighbours limit, which is to be expected since the number of connections and thus the potential energy increases with the dimension while the total energy is kept constant so the kinetic energy and thus the temperature must decrease. Finally, it is worth mentioning that in the energy range $U>U_c$, the temperature versus energy curves for $d>1$ are surprisingly well approximated by Equation (12) with the Langenvin function evaluated at $1/2dT$.

Figure 5 shows the M vs. U curve for the long-range regime compared with the nearest neighbours regime curve for 1, 2 and 3 dimensions with the same parameter values and initial conditions used in Figure 4. As expected, in one dimension one gets $M=0$ for any value of U since there is no phase transition for $d=1$. For $d=3$ the order parameter M vanishes at the critical value $U\equiv U_{3d}\sim 0.6$ for which the ferromagnetic-paramagnetic transition takes place. In dimension $d=2$, it is observed a less sharp decay of the magnetization, which vanishes for $U\equiv U_{2d}\sim 3.5$.

3.3. Size-Scaling of the Largest Lyapunov Exponent

The maximum Lyapunov exponent characterizes the rate of separation of two initially close orbits. It is defined as $\lambda ={lim}_{t\to \infty }\frac{1}{t}ln\frac{d\left(t\right)}{d\left(0\right)}$, where $d\left(t\right)$ is the distance between orbits, measured in phase space, at time t. In order to numerically calculate it, we consider the usual approach given in Ref. [51]. Figure 6 shows the obtained values of the maximum Lyapunov exponents for one-dimensional systems of sizes $\mathcal{N}=1000,2000,\cdots ,$ 10,000, energy per particle $U=5.0$ and interaction range $\alpha =0.0,0.2,\cdots ,1.6$. An average over five different realizations of disorder has been taken in order to improve statistics. In a log-log plot the relation $\lambda \sim {\mathcal{N}}^{-\kappa }$ takes the form of a straight line with slope $-\kappa$, which has then been obtained via a least squares fitting. Figure 7 and Figure 8 show the corresponding results for two- and three-dimensional systems for the same values of $\alpha /d$ as in Figure 6 and the system sizes specified in the captions. Figure 9 displays the values of $\kappa$ obtained in Figure 6, Figure 7 and Figure 8 versus $\alpha /d$, where a clear collapse of the data is observed, with $\kappa \left(0\right)\sim 0.3$ and $\kappa \sim 0$ for $\alpha /d\ge {\alpha }_c\equiv 1.6$. The hypothesis $\kappa =\kappa (\alpha /d)$ is thus confirmed. It is worth mentioning that we have an intermediate regime $1\le \alpha /d\le {\alpha }_c$ with zero Lyapunov exponent and short-range interactions. This curious anomaly demands further study to check whether it is a finite-size effect or something else.

As mentioned in the beginning of this section, we have used the standard fourth-order Runge-Kutta algorithm for the integration of the equations of motion, by monitoring that the total energy was preserved within a relative precision ${10}^{-5}$ (at least). However, further simulations using specially designed symplectic methods for Hamiltonian problems [48,49,50] are welcome.

4. Conclusions

We have here focused on the classical inertial Heisenberg model in $d=1,2,3$ dimensions with arbitrarily-ranged interactions, characterized by a power-law decay $r^{-\alpha }$ ($\alpha \ge 0$), with the distance r between rotators. Such an $\alpha$-Heisenberg model has been studied previously only in dimension $d=1$ [47], with a special emphasis on the time evolution of its kinetic temperature, where it was shown to exhibit a clear QSS in the long-range-interaction regime, with its duration decreasing for increasing values of $\alpha$. Herein, the properties of the model have been preliminary analyzed in $d=1,2,3$ dimensions, within first-principle numerical calculations, particularly the one-particle moment and energy distributions, as well as the largest Lyapunov exponent.

The Lyapunov exponent scales with the size through the exponent $\kappa$, which was shown to depend only on the ratio$\alpha /d$. However, intriguingly enough, the results differ from those already known for the quite analogous $\alpha$-XY and $\alpha$-Fermi-Pasta-Ulam models. Indeed, for the XY and Fermi-Pasta-Ulam models, $\kappa$ monotonically decreases for increasing $\alpha /d$ and vanishes at $\alpha /d=1$. This fact is in variance with the result obtained for the Heisenberg model, for which $\kappa$ decreases for increasing $\alpha /d$ but only appears to vanish above$\alpha /d=1$, nearly at $\alpha /d\simeq 1.6$.

We have illustrated the one-particle momenta and energy distributions in the long-range region with one example, namely for $\alpha /d=0.9$, and the results have shown a clear collapse of data and were well fitted with a q-Gaussian and a q-exponential respectively. As found in the recent study of the $\alpha$-XY model [34], an interesting discrepancy in the values of the index q occurred, leading to $q_L>q_E$. Whether this may be a real effect, or a consequence of finite-size effects, remains a point for further investigation.In addition, we have analyzed similar distributions in the short-range region, namely $\alpha /d=1.6$. In this case, the results are quite more complex, in the sense that no clear collapse, and consequently, no adequate fitting occurs.

The present numerical analysis is consistent with the theoretical prevision of a change of behavior at $\alpha /d=1$, in agreement with recent results obtained in other long-range-interaction models, such as a generalized Fermi-Pasta-Ulam model [35,36,37,38,39], complex networks with preferential attachment growth [40,41], and a system of particles under overdamped motion interacting repulsively [42]. In future works we intend to cover the entire $\alpha /d\ge 0$ region in order to clarify further this point, and particularly, to investigate why some of the Heisenberg results found herein sensibly differ from those available in the literature for the XY and Fermi-Pasta-Ulam models.