Relaxation to equilibrium in models of classical spins with long-range interactions

TL;DR

This paper investigates how fluctuations influence the relaxation process of classical Heisenberg spins with long-range interactions, revealing that environmental noise can eliminate quasistationary states and lead to rapid equilibrium attainment.

Contribution

It demonstrates that environmental fluctuations cause the disappearance of quasistationary states in long-range spin systems, contrasting deterministic and stochastic dynamics.

Findings

Deterministic dynamics exhibit quasistationary states with diverging lifetimes.

Strong environmental noise induces rapid, size-independent relaxation to equilibrium.

Glauber Monte Carlo dynamics also show fast relaxation, confirming the effect of fluctuations.

Abstract

For a model long-range interacting system of classical Heisenberg spins, we study how fluctuations, such as those arising from having a finite system size or through interaction with the environment, affect the dynamical process of relaxation to Boltzmann-Gibbs equilibrium. Under deterministic spin precessional dynamics, we unveil the full range of quasistationary behavior observed during relaxation to equilibrium, whereby the system is trapped in nonequilibrium states for times that diverge with the system size. The corresponding stochastic dynamics, modeling interaction with the environment and constructed in the spirit of the stochastic Landau-Lifshitz-Gilbert equation, however shows a fast relaxation to equilibrium on a size-independent timescale and no signature of quasistationarity, provided the noise is strong enough. Similar fast relaxation is also seen in Glauber Monte Carlo dynamics of the model, thus establishing the ubiquity of what has been reported earlier in particle dynamics (hence distinct from the spin dynamics considered here) of long-range interacting systems, that quasistationarity observed in deterministic dynamics is washed away by fluctuations induced through contact with the environment.

Equations

$$H=-\frac{J}{2N}\sum _{i,j=1}^N{\mathbf{S}}_i\cdot {\mathbf{S}}_j+D\sum _{i=1}^N{S}_{iz}^{2n},$$

$$\{S_{i\alpha },S_{j\gamma }\}={\delta }_{ij}{ϵ}_{\alpha \gamma \delta }{S}_{j\delta },$$

$$\{A,B\}=\sum _{i=1}^N{\mathbf{S}}_i\cdot \frac{\partial A}{\partial {\mathbf{S}}_i}\times \frac{\partial B}{\partial {\mathbf{S}}_i}.$$

$$\mathbf{m}\equiv \frac{{\sum }_{i=1}^N{\mathbf{S}}_i}{N}=(m_x,m_y,m_z)$$

$${\dot{\mathbf{S}}}_i=\{{\mathbf{S}}_i,H\};i=1,2,\dots ,N,$$

$${\dot{S}}_{ix}$$

$${\dot{S}}_{iy}$$

$${\dot{S}}_{iz}$$

$${\dot{\theta }}_i$$

$${\dot{\varphi }}_i$$

$${\dot{\mathbf{S}}}_i={\mathbf{S}}_i\times {\mathbf{h}}_i^{\mathrm{eff}},$$

$${\mathbf{h}}_i^{\mathrm{eff}}=-\frac{\partial H}{\partial {\mathbf{S}}_i}=\mathbf{m}+{\mathbf{h}}_i^{\mathrm{aniso}};{\mathbf{h}}_i^{\mathrm{aniso}}\equiv (0,0,-2nD{S}_{iz}^{2n-1}).$$

$${\dot{\mathbf{S}}}_i={\mathbf{S}}_i\times \left({\mathbf{h}}_i^{\mathrm{eff}}+{\mathbf{\eta }}_i(t)\right)-\gamma {\mathbf{S}}_i\times \left({\mathbf{S}}_i\times \left({\mathbf{h}}_i^{\mathrm{eff}}+{\mathbf{\eta }}_i(t)\right)\right),$$

$$⟨{\eta }_{i\mu }(t)⟩=0,⟨{\eta }_{i\mu }(t){\eta }_{j\nu }(t^{\prime })⟩=2{\delta }_{ij}{\delta }_{\mu \nu }\mathcal{D}\delta (t-t^{\prime }).$$

$$\mathcal{D}=\gamma T/(1+{\gamma }^2),$$

$$h(\theta ,\varphi ,m_x,m_y,m_z)\equiv -m_x\sin \theta \cos \varphi -m_y\sin \theta \sin \varphi -m_z\cos \theta +D{\cos }^{2n}\theta .$$

$$f_{\mathrm{eq}}(\theta ,\varphi )\propto \exp (-\beta h(\theta ,\varphi ,m_x^{\mathrm{eq}},m_y^{\mathrm{eq}},m_z^{\mathrm{eq}})),$$

$$(m_x^{\mathrm{eq}},m_y^{\mathrm{eq}},m_z^{\mathrm{eq}})=\frac{\int \int \sin \theta \mathrm{d}\theta \mathrm{d}\varphi (m_x^{\mathrm{eq}},m_y^{\mathrm{eq}},m_z^{\mathrm{eq}})f_{\mathrm{eq}}(\theta ,\varphi )}{\int \sin \theta \mathrm{d}\theta \mathrm{d}\varphi f_{\mathrm{eq}}(\theta ,\varphi )}.$$

$$m_{\mathrm{eq}}=\frac{\int \mathrm{d}\theta \mathrm{d}\varphi {\sin }^2\theta \cos \varphi e^{\beta m_{\mathrm{eq}}\sin \theta \cos \varphi -\beta D{\cos }^{2n}\theta }}{\int \mathrm{d}\theta \mathrm{d}\varphi \sin \theta e^{\beta m_{\mathrm{eq}}\sin \theta \cos \varphi -\beta D{\cos }^{2n}\theta }}.$$

$$ϵ=-\frac{1}{2}{m}_{\mathrm{eq}}^2+D⟨{\cos }^{2n}\theta {⟩}_{\mathrm{eq}},$$

$$⟨{\cos }^{2n}\theta {⟩}_{\mathrm{eq}}=\frac{\int \mathrm{d}\theta \mathrm{d}\varphi \sin \theta {\cos }^{2n}\theta e^{\beta m_{\mathrm{eq}}\sin \theta \cos \varphi -\beta D{\cos }^{2n}\theta }}{\int \mathrm{d}\theta \mathrm{d}\varphi \sin \theta e^{\beta m_{\mathrm{eq}}\sin \theta \cos \varphi -\beta D{\cos }^{2n}\theta }}.$$

$$\displaystyle m_{\rm eq}\Big{(}\int{\rm d}\theta{\rm d}\phi\sin\theta~{}e^{-\beta D\cos^{2n}\theta}-~{}\beta\int{\rm d}\theta{\rm d}\phi\sin^{3}\theta\cos^{2}\phi~{}e^{-\beta D\cos^{2n}\theta}\Big{)}=0.$$

$$1-\frac{2}{\beta_{c}}=\frac{\Gamma(3/(2n))-\Gamma(3/(2n),\beta_{c}D)}{(\beta_{c}D)^{2/(2n)}\Big{[}\Gamma(1/(2n))-\Gamma(1/(2n),\beta_{c}D)\Big{]}}.$$

$${ϵ}_c=D\frac{\int \mathrm{d}\theta \mathrm{d}\varphi \sin \theta {\cos }^{2n}\theta e^{-\beta D{\cos }^{2n}\theta }}{\int \mathrm{d}\theta \mathrm{d}\varphi \sin \theta e^{-\beta D{\cos }^{2n}\theta }};$$

$$\epsilon_{c}=\frac{\Gamma(1/(2n))-2n\Gamma(1+1/(2n),\beta_{c}D)}{2n\beta_{c}\Big{[}\Gamma(1/(2n))-\Gamma(1/(2n),\beta_{c}D)\Big{]}}.$$

$$\frac{\partial P_0(\mathbf{S},t)}{\partial t}+\frac{\partial }{\partial \mathbf{S}}\cdot (\mathbf{S}\times {\mathbf{h}}^{\mathrm{eff},0})P_0=0,$$

$${\mathbf{h}}^{\mathrm{eff},0}\equiv {\mathbf{h}}^{\mathrm{eff},0}[P_0]=\mathbf{m}[P_0]+(0,0,-2nD{S}_z^{2n-1});\mathbf{m}[P_0]\equiv \int \mathrm{d}\mathbf{S}\mathbf{S}{P}_0(\mathbf{S},t).$$

$$\displaystyle\Big{[}(S_{y}h^{\rm eff}_{z}-S_{z}h^{\rm eff}_{y})(\cos\theta\cos\phi~{}\hat{\theta}-\sin\phi~{}\hat{\phi})+(S_{z}h^{\rm eff}_{x}-S_{x}h^{\rm eff}_{z})(\cos\theta\sin\phi~{}\hat{\theta}+\cos\phi~{}\hat{\phi})$$

$$\displaystyle-(S_{x}h^{\rm eff}_{y}-S_{y}h^{\rm eff}_{x})\sin\theta~{}\hat{\theta}\Big{]}\cdot\Big{[}\hat{\theta}\frac{\partial}{\partial\theta}+\frac{1}{\sin\theta}\hat{\phi}\frac{\partial}{\partial\phi}\Big{]}f(\theta,\phi,t)$$

$$\displaystyle=\Big{[}\cos\theta\cos\phi(S_{y}h^{\rm eff}_{z}-S_{z}h^{\rm eff}_{y})+\cos\theta\sin\phi(S_{z}h^{\rm eff}_{x}-S_{x}h^{\rm eff}_{z})-\sin\theta(S_{x}h^{\rm eff}_{y}-S_{y}h^{\rm eff}_{x})\Big{]}\frac{\partial f}{\partial\theta}$$

$$\displaystyle+\Big{[}-\sin\phi(S_{y}h^{\rm eff}_{z}-S_{z}h^{\rm eff}_{y})+\cos\phi(S_{z}h^{\rm eff}_{x}-S_{x}h^{\rm eff}_{z})\Big{]}\frac{1}{\sin\theta}\frac{\partial f}{\partial\phi}$$

$$\displaystyle=-\Big{[}m_{y}\cos\phi-m_{x}\sin\phi\Big{]}\frac{\partial f}{\partial\theta}+\Big{[}m_{x}\cot\theta\cos\phi+m_{y}\cot\theta\sin\phi-m_{z}+(2n)D\cos^{2n-1}\theta\Big{]}\frac{\partial f}{\partial\phi}.$$

$$\displaystyle\frac{\partial f}{\partial t}=\Big{(}m_{y}\cos\phi-m_{x}\sin\phi\Big{)}\frac{\partial f}{\partial\theta}-\Big{(}m_{x}\cot\theta\cos\phi+m_{y}\cot\theta\sin\phi-m_{z}+(2n)D\cos^{2n-1}\theta\Big{)}\frac{\partial f}{\partial\phi},$$

$$(m_x,m_y,m_z)\equiv (m_x,m_y,m_z)[f]=\int \sin {\theta }^{\prime }\mathrm{d}{\theta }^{\prime }\mathrm{d}{\varphi }^{\prime }(\sin {\theta }^{\prime }\cos {\varphi }^{\prime },\sin {\theta }^{\prime }\sin {\varphi }^{\prime },\cos {\theta }^{\prime })f({\theta }^{\prime },{\varphi }^{\prime },t).$$

Full text

Keywords: Stationary states, Metastable states, Stochastic particle dynamics

1. 1 Introduction and model of study

2. 2 Previous studies

3. 3 Our queries and summary of results obtained

4. 4 Equilibrium properties

5. 5 Analysis of the deterministic dynamics (11)

6. 5.1 Behavior in the limit $N\to\infty$

7. 5.2 Behavior for finite $N$

8. 5.3 Numerical results

9. 6 Analysis of the stochastic dynamics (13)

10. 6.1 Behavior in the limit $N\to\infty$

11. 6.2 Behavior for finite $N$

12. 6.3 An alternative to dynamics (13): A Monte Carlo dynamical scheme

13. 6.4 Numerical results

14. 7 Conclusions

15. 8 Acknowledgements

16. 9 Appendix A: Derivation of Eqs. (43), (44), and (45) of the main text

17. 10 Appendix B: Derivation of Eqs. (26) and (65) of the main text

18. 11 Appendix C: Numerical scheme for integrating Eq. (13)

1 Introduction and model of study

Stochasticity is present in any statistical system constituted by a finite number of interacting degrees of freedom, which is known to induce fluctuations in both static and time-dependent observables of the system, thereby affecting their statistical properties. Stochasticity may arise due to sampling of initial conditions and due to interaction with the external environment. It is evidently of interest to investigate how these two sources of stochasticity interplay in dictating the long-time state of the system, and in particular, in predicting the values of the macroscopic observables the system attains in the stationary state. In this work, we explore the aforementioned issues within the ambit of a model many-body interacting system comprising classical spins that are interacting with one another through an inter-particle potential that is long-ranged in nature. Namely, the interparticle potential decays rather slowly as a function of the separation $r$ between the particles, specifically, as $1/r^{\alpha};~{}0\leq\alpha\leq d$, with $d$ being the embedding spatial dimension of the system [1, 2, 3, 4, 5].

Long-range interacting (LRI) systems may be found at all length scales, from atomic to astrophysical, and are known to exhibit a range of physical phenomena that appear counterintuitive when viewed vis-à-vis short-range systems for which the interaction has a finite range. A basic property that distinguishes LRI systems from short-range ones is the violation of additivity, whereby a macroscopic LRI system cannot be divided into independent macroscopic subparts so that thermodynamic quantities referring to the subparts add up to yield the corresponding values for the composite system. While non-additivity results in such unusual features as inequivalent equilibrium ensembles and a non-concave entropy function, more striking is its consequence on dynamical properties, namely, that an isolated LRI system relaxes to equilibrium over a time that diverges with the number of degrees of freedom [3]. Consequently, the larger the system is, the longer is the time that it takes to attain equilibrium, resulting in slowly-evolving nonequilibrium states being directly accessible to experimental observations [3].

Our model of study consists of $N$ globally-coupled classical Heisenberg spins of unit length denoted by ${\mathbf{S}}_{i}=(S_{ix},S_{iy},S_{iz});~{}i=1,2,\ldots,N$. Expressed in terms of spherical polar angles $\theta_{i}\in[0,\pi]$ and $\phi_{i}\in[0,2\pi)$, one has $S_{ix}=\sin\theta_{i}\cos\phi_{i},~{}S_{iy}=\sin\theta_{i}\sin\phi_{i},~{}S_{iz}=\cos\theta_{i}$. The Hamiltonian of the system is given by

[TABLE]

where $n>0$ is an integer. The spin components satisfy

[TABLE]

where $\epsilon_{\alpha\gamma\delta}$ is the fully antisymmetric Levi-Civita symbol. Here and in the following, we use Roman indices to label the spins and Greek indices to denote the spin components. Noting that the canonical variables for the $i$-th spin are $\cos\theta_{i}$ and $\phi_{i}$, the Poisson bracket $\{,\}$ appearing in Eq. (2) is defined for two functions $A,B$ of the spins by [6] $\{A,B\}\equiv\sum_{i=1}^{N}(\partial A/\partial\phi_{i}\,\partial B/\partial(\cos\theta_{i})-\partial A/\partial(\cos\theta_{i})\,\partial B/\partial\phi_{i})$, which may be re-expressed as

[TABLE]

We now explain the various terms appearing in Eq. (1). Here, the first term with $J>0$ on the right hand side models a ferromagnetic mean-field interaction between the spins. On the other hand, the second term with $D>0$ on the right hand side accounts for local anisotropy: for example, $n=1$ (respectively, $n=2$) models quadratic (respectively, quartic) anisotropy, and will be referred to below as the quadratic (respectively, the quartic) model. Single-spin Hamiltonian of Heisenberg spins and involving quadratic and quartic terms has been considered previously in the literature, see, e.g., Ref. [7]. The anisotropy term in Eq. (1) lowers energy by having the magnetization vector

[TABLE]

pointing in the $xy$ plane. The length of the magnetization vector is given by $m\equiv\sqrt{m_{x}^{2}+m_{y}^{2}+m_{z}^{2}}$. The coupling constant $J$ in Eq. (1) has been scaled down by the system size $N$ to order to make the energy extensive, in accordance with the Kac prescription [8]. The system (1) is however intrinsically non-additive, since extensivity does not guarantee additivity, although the converse is true. In the following, we set $J$ to unity without loss of generality, and also take unity for the Boltzmann constant.

In dimensionless time, the dynamics of the system (1) is governed by the set of coupled first-order differential equations [6]

[TABLE]

where the dot denotes derivative with respect to time. Using Eq. (5), we obtain the dynamical evolution of the spin components as

[TABLE]

Taking the vector dot product of both sides of Eq. (5) with ${\bf S}_{i}$, it is easily seen that the dynamics conserves the length of each spin. Summing Eq. (8) over $i$, we find that $m_{z}$ is a constant of motion. Note that for the special case of no anisotropy ($D=0$), the total magnetization $m$ is a constant of motion. The total energy of the system is a constant of motion under the dynamical evolution (5), and as such, the latter models microcanonical dynamics of the system (1). We remark that the dynamical setting of Eq. (5) is very different from that involving particles characterized by generalized coordinates and momenta and time evolution governed by a Hamiltonian given by a sum of a kinetic and a potential energy contribution, e.g., that of the celebrated Hamiltonian mean-field (HMF) model [3]. As a result, none of the results, static or dynamic, derived for the latter may be a priori expected to apply to the model (1). From Eqs. (6)-(8), we obtain the time evolution of the variables $\theta_{i}$ and $\phi_{i}$ as

[TABLE]

Equations (6)-(8) may be interpreted as the precessional dynamics of the spins in an effective magnetic field ${\bf h}^{\rm eff}_{i}\equiv{\bf h}^{\rm eff}_{i}(\{{\bf S}_{i}\})$:

[TABLE]

where ${\bf h}^{\rm eff}_{i}$, the effective field for the $i$-th spin, is obtained from the Hamiltonian (1) as

[TABLE]

Thus, the effective magnetic field has a global and a local contribution, with the former being due to the magnetization set up in the system by the effect of all the spins, and the latter being due to the field ${\bf h}^{\rm aniso}_{i}$ set up for individual spins by the anisotropy term in the Hamiltonian (1).

The paper is organized as follows. In Section 2, we summarize previous studies of model (1) for $n=1$, which is followed in Section 3 by listing of our queries in this work, namely, the relaxation properties of the deterministic dynamics (5) and the corresponding stochastic dynamics constructed in the spirit of the stochastic Landau-Lifshitz-Gilbert (LLG) equation [9]. Here, we also give a summary of results obtained in this work. The following sections are then devoted to a derivation of these results. We start with a derivation of the equilibrium properties of the model (1) in Section 4. This is followed in Section  5 by a study of the deterministic dynamics (5) in the limit $N\to\infty$ in terms of the so-called Vlasov equation in Subsection 5.1. Here we also study linear stability of a representative stationary state of the Vlasov equation, and demarcate for two representative values of $n$ (namely, $n=1,2$) regions in the parameter space where the state is stable under the Vlasov evolution. Subsection 5.2 is devoted to a discussion on the behavior of the dynamics when $N$ is large but finite. The behavior of the stochastic dynamics in the limit $N\to\infty$ and in the case when $N$ is large but finite are taken up in Section 6. In this section, we also discuss a Monte Carlo scheme that serves as an alternative to the stochastic LLG scheme to study effects of noise on the deterministic dynamics (5). All throughout, we provide numerical checks of our theoretical predictions, considering $n=1,2$ in the Hamiltonian (1). We draw our conclusions in Section 7. Some of the technical details of our analytic and numerical analysis are collected in the three appendices.

2 Previous studies

The quadratic model was first considered in Ref. [6] that addressed the equilibrium and relaxational properties of the model. The system was shown to exhibit in Boltzmann-Gibbs microcanonical equilibrium a magnetized (equilibrium magnetization $m_{\rm eq}\neq 0$) phase at low values of the energy $\epsilon$ per spin and a nonmagnetized ($m_{\rm eq}=0$) phase at high values, with a continuous transition between the two occurring at a critical value $\epsilon_{c}$. It was established that within microcanonical dynamics and for a class of nonmagnetized initial states, there exists a threshold energy $\epsilon^{\star}<\epsilon_{c}$, such that in the energy range $\epsilon^{\star}<\epsilon<\epsilon_{c}$, relaxation to equilibrium magnetized state occurs over a time that scales superlinearly with $N$ [6, 10]. On the other hand, for energies $\epsilon<\epsilon^{\star}$, the dynamics shows a fast relaxation out of the initial nonmagnetized state over a time that scales as logarithm of $N$. The particular initial state that was considered was the so-called waterbag (WB) state, in which the spins have $\phi$’s chosen independently and uniformly over the interval $[0,2\pi)$ and the $\theta$’s chosen independently and uniformly over an interval symmetric about $\pi/2$, that is, over the interval $[\pi/2-a,\pi/2+a]$, with $a$ being a real positive quantity. These results, obtained on the basis of numerical integration of the equations of motion, were complimented by an analytical study in the limit $N\to\infty$ of the time evolution, à la a Vlasov-type equation, of the single-spin phase space distribution. The distribution counts the fraction of the total number of spins that have given $\theta$ and $\phi$ values. It was demonstrated that the distribution associated to the WB state is stationary under the Vlasov evolution, but is unstable for energies below $\epsilon^{\star}$ and stable for energies above. For finite $N$, the eventual relaxation to equilibrium observed for energies $\epsilon^{\star}<\epsilon<\epsilon_{c}$ was accounted as due to statistical fluctuations adding non-zero finite-$N$ corrections to the Vlasov equation that are at least of order greater than $1/N$ [11]. The WB state that for energies $\epsilon^{\star}<\epsilon<\epsilon_{c}$ is stationary and stable in an infinite system but which shows a slow evolution for finite $N$ exemplifies the so-called quasistationary states (QSSs) [3].

3 Our queries and summary of results obtained

Starting from the premises discussed in the previous section, we pursue in this work a detailed characterization of the relaxational dynamics and its ubiquity in the context of long-range spin models, by considering the model (1) for general values of $n$. We study for general $n$ the Boltzmann-Gibbs microcanonical equilibrium properties of the model, deriving in particular an expression for the continuous phase transition point $\epsilon_{c}(n)$, such that the system is in a magnetized phase for lower energies and in a nonmagnetized phase for higher energies. Though not guaranteed for LRI systems [3], by virtue of the model (1) exhibiting a continuous transition in equilibrium, we conclude by invoking established results [12] that microcanonical and canonical ensembles are equivalent in equilibrium. Consequently, one may associate to every value of the conserved microcanonical energy density $\epsilon$ a temperature $T$ of the system in canonical equilibrium that guarantees that the average energy in canonical equilibrium equals $\epsilon$ in the limit $N\to\infty$. This allows to also derive the phase diagram of the model (1) in canonical equilibrium.

The WB single-spin distribution is non-analytic at $\theta=\pi/2\pm a$, and one may wonder as to whether such a peculiar feature led to quasitationarity in the $n=1$ model reported in Refs. [6, 11] and summarized in the preceding section. As a counterpoint, and to demonstrate that quasistationarity is rather generic to the model (1), we consider as initial states suitably smoothened versions of the WB state, the so-called Fermi-Dirac (FD) state, for which the single-spin distribution is a perfectly analytic function, and study its evolution in time. A linear stability analysis of the FD state under the infinite-$N$ Vlasov dynamics establishes the existence of a threshold energy value $\epsilon^{\star}(n)$, such that the state is stationary but linearly unstable under the dynamics for energies $\epsilon<\epsilon^{\star}(n)$. For finite $N$, we establish that the relaxation to equilibrium occurs as a two-step process: an initial relaxation from the FD state to a magnetized QSS, followed by a relaxation of the latter over a timescale $\sim N$ to Boltzmann-Gibbs microcanonical equilibrium. The magnetized QSS has thus a lifetime $\sim N$. For energies $\epsilon^{\star}(n)<\epsilon<\epsilon_{c}(n)$, however, the FD state is dynamically stable under the Vlasov evolution, exhibiting for finite $N$ a relaxation towards equilibrium over a scale $\sim N^{\alpha}$, where the exponent $\alpha$ has an essential dependence on $n$. In this case, one concludes observing a nonmagnetized QSS with a lifetime $\sim N^{\alpha}$. As for $\alpha$, while one obtains for $n=1$ the value $\alpha=3/2$ (as opposed to the value $\alpha=2$ for the WB state reported in Ref. [10]), one observes a linear dependence ($\alpha=1$) for the quartic model. Note that the two-step relaxation process for energies $\epsilon<\epsilon^{\star}(n)$ was not reported in previous studies of the model (1), see Refs. [6, 11], and is being reported here for the first time. While magnetization $m$ turns out to be a useful macroscopic observable to monitor in order to establish the aforementioned relaxation scenario, it does not serve the purpose when considering energies $\epsilon>\epsilon_{c}(n)$ where both the initial FD and the final Boltzmann-Gibbs microcanonical equilibrium state are nonmagnetized. Here, by identifying a suitable observable (e.g., $\sum_{i=1}^{N}\cos^{4}\theta_{i}/N$ and $\sum_{i=1}^{N}\cos^{2}\theta_{i}/N$ for respectively the quadratic and the quartic model), we show that the relaxation to equilibrium occurs over a timescale that has an $N$ dependence distinct from what was observed for magnetization relaxation for energies $\epsilon^{\star}(n)<\epsilon<\epsilon_{c}(n)$. Namely, the relaxation time scales as $N^{2}$ for the quadratic model and as $N^{3/2}$ for the quartic model. We may thus conclude for energies $\epsilon>\epsilon_{c}(n)$ the existence of a nonmagnetized QSS with a lifetime that diverges with the system size. We stress that the existence of QSSs with lifetimes $\sim N^{2}$ was not discussed in previous studies [6, 11], and it is here that we report on such states for the first time.

Our next issue of investigation is the robustness of QSSs with respect to fluctuations induced through contact with the external environment. Modelling the environment as a heat bath, previous studies of Hamiltonian particle dynamics (e.g., that of the HMF model) have invoked a scheme of coupling to the environment that allows for energy exchange and consequent stochastic Langevin evolution of the system. These studies have suggested a fast relaxation to equilibrium over a size-independent timescale provided the noise is strong enough [13, 14, 15]. In the context of the model (1), in order to assess the effects of noise induced by the external environment, we study a stochastic version of the dynamics (11) that considers the effective field ${\bf h}_{i}^{\rm eff}$, see Eq. (11), to have an additional stochastic component due to interaction with the environment. The resulting dynamics, built in the spirit of the stochastic Landau-Lifshitz-Gilbert equation well known in studies of dynamical properties of magnetic systems (see Ref. [9] for a review), reads

[TABLE]

where the second term on the right represents dissipation with the real parameter $\gamma>0$ being the dissipation constant, and $\boldsymbol{\eta}_{i}(t)$ is a Gaussian white noise with independent components that satisfy

[TABLE]

Here, ${\cal{D}}>0$ is a real constant that characterizes the strength of the noise. Note that the stochastic dynamics (13) conserves the length of each spin, as does the deterministic dynamics (11). The former models dynamics within a canonical ensemble for which the energy is not conserved during the dynamical evolution, while, as already mentioned earlier, the latter models energy-conserving microcanonical dynamics.

The presence of noise in Eq. (13) makes the state of the system at a given time, characterized by the set of values $\{{\bf S}_{i}(t)\}$, to vary from one realization of the dynamics to another, even when starting every time from the same initial condition. Although Eq. (13) has the flavor of Langevin dynamics, it is different in that the noise and dissipation terms are incorporated in a way that it has the desirable feature of keeping the length of each spin to be unity at all times during the dynamical evolution. Since the noise terms in Eq. (13) depend on the state of the system, itself stochastic in nature, the noise is said to be multiplicative in common parlance. As we argue later in the paper, requiring the dynamics (13) to relax at long times to Boltzmann-Gibbs canonical equilibrium at a given temperature $T$ fixes the constant ${\cal D}$ to be related to $\gamma$ in the manner

[TABLE]

a choice we also consider in this work. Our numerical simulation of the dynamics (13) follows the scheme detailed in Appendix C. The results show that in presence of strong-enough noise, the system shows a fast relaxation to Boltzmann-Gibbs equilibrium on a size-independent time scale, with no existence of intermediate quasistationary states. We also implement a Monte Carlo scheme as an alternative to the dynamics (13) to study the effects of environment-induced noise on the dynamics (5). On implementing the scheme, we find similar to the study of the dynamics (13) a fast relaxation to Boltzmann-Gibbs equilibrium on a size-independent timescale. Our studies thus serve to reaffirm what has been observed earlier in particle dynamics of LRI systems, namely, that quasistationarity, observed in conservative dynamics, is completely washed away in presence of stochasticity in the dynamics.

4 Equilibrium properties

In this section, we investigate the properties of the system (1) in the thermodynamic limit $N\to\infty$ and in canonical equilibrium at temperature $T=1/\beta$. Note that model (1) is a mean-field system that describes the motion of a spin moving in a self-consistent mean-field generated by its interaction with all the spins, with the single-spin Hamiltonian given by

[TABLE]

Consequently, it is rather straightforward to write down exact expressions for the average magnetization and the average energy in equilibrium and in the thermodynamic limit. For example, the single-spin equilibrium distribution is given by

[TABLE]

with the equilibrium magnetization components $(m_{x}^{\rm eq},m_{y}^{\rm eq},m_{z}^{\rm eq})$ obeying the self-consistent equation

[TABLE]

With $D>0$, so that the system orders in the $xy$-plane, we may choose the ordering direction to be along $x$ without loss of generality, yielding $m_{x}^{\rm eq}\neq 0,m_{y}^{\rm eq}=m_{z}^{\rm eq}=0$. From Eq. (18), we thus obtain for $m_{\rm eq}\equiv m_{x}^{\rm eq}$ the equation [6]

[TABLE]

The average energy per spin equals [6]

[TABLE]

with

[TABLE]

From the fact that the model (1) with $n=1$ shows a continuous phase transition in magnetization across critical inverse temperature $\beta_{c}=1/T_{c}$ [6], we may anticipate that so is the case for general $n$. Consequently, we may consider Eq. (19) close to the critical point, i.e., for \beta\lower 2.0pt\hbox{,\buildrel\scriptstyle>\over{\scriptstyle\sim},}\beta_{c}, when $m_{\rm eq}$ is small so that the equation may be expanded to leading order in $m_{\rm eq}$, as

[TABLE]

With $m_{\rm eq}\neq 0$, one obtains $\beta_{c}$ as the value of $\beta$ that sets the bracketed quantity to zero; on performing the integrals, one obtains $\beta_{c}$ to be satisfying

[TABLE]

Here, $\Gamma(s)$ is the Gamma function, while $\Gamma(s,x)$ is the upper incomplete Gamma function. At the critical point, when we have $m_{\rm eq}=0$, one obtains the critical energy density as $\epsilon_{c}=D\langle\cos^{2n}\theta\rangle_{\rm eq}$, that is,

[TABLE]

on performing the integrals, we get

[TABLE]

Note that for $n=1$, one may check using the above expressions that $1-2/\beta_{c}=1/(2\beta_{c}D)-e^{-\beta_{c}D}/(\sqrt{\pi\beta_{c}D}\mathrm{Erf}[\sqrt{\beta_{c}D}])$ and that \epsilon_{c}=D\Big{(}1-2/\beta_{c}\Big{)}, where $\mathrm{Erf}[x]\equiv(2/\sqrt{\pi})\int_{0}^{x}{\rm d}t~{}e^{-t^{2}}$ is the error function, as was reported in Ref. [6].

Since the phase transition exhibited by the model (1) is a continuous one, the canonical and microcanonical ensemble properties in equilibrium would be equivalent [12], and hence, Eq. (25) also gives the conserved microcanonical energy density at the transition point. Figure 1 shows for $n=1,2$ the energy density $\epsilon_{c}$ as a function of $D$, obtained by first solving numerically for a given $D$ the transcendental equation (23) for $\beta_{c}$ and then using the obtained value of $\beta_{c}$ in Eq. (25). Moreover, one may construct a one-to-one mapping between a value of microcanonical equilibrium energy density $\epsilon$ and canonical equilibrium temperature $T$ by first solving Eq. (19) at a given $T$ to obtain the equilibrium magnetization $m_{\rm eq}$, then substituting in Eq. (20) to obtain the corresponding average energy in canonical equilibrium, and finally demanding that the latter is the conserved energy density in microcanonical equilibrium. On carrying out this program for $n=1$ and $D=5.0$, one obtains the results shown in Fig. 2, where we also show $m_{\rm eq}$ as a function of microcanonical energy density $\epsilon$.

5 Analysis of the deterministic

dynamics (11)

5.1 Behavior in the limit $N\to\infty$

We now discuss how the model (1) in the limit $N\to\infty$ and under the dynamical evolution (11) relaxes to equilibrium while starting far from it, by invoking the corresponding Vlasov equation. The relaxation may be characterized by monitoring the time evolution of the single-spin distribution function $P_{0}({\bf S},t)$ normalized as $\int{\rm d}{\bf S}~{}P_{0}({\bf S},t)=1~{}\forall~{}t$. As detailed in Appendix B, the time evolution of $P_{0}({\bf S},t)$ follows the Vlasov equation

[TABLE]

where we have

[TABLE]

For later purpose, it is convenient to consider the single-spin distribution $f(\theta,\phi,t)$, defined such that $f(\theta,\phi,t)\sin\theta{\rm d}\theta{\rm d}\phi$ is the probability to have a spin at time $t$ with its angles between $\theta$ and $\theta+{\rm d}\theta$ and between $\phi$ and $\phi+{\rm d}\phi$, and which is related to $P_{0}({\bf S},t)$ as $f(\theta,\phi,t)=P_{0}({\bf S},t)$, with the normalization $\int_{0}^{2\pi}{\rm d}\phi\int_{0}^{\pi}{\rm d}\theta~{}\sin\theta~{}f(\theta,\phi,t)=1~{}\forall~{}t$. To obtain the time evolution of $f$, let us express the second term on the left hand side of Eq. (26), which equals $({\bf S}\times{\bf h}^{{\rm eff},0})\cdot\partial P_{0}/\partial{\bf S}$, in spherical polar coordinates with unit vectors $(\hat{r},\hat{\theta},\hat{\Phi})$ and $r=\sqrt{S^{2}_{x}+S^{2}_{y}+S^{2}_{z}}=1$. We get

[TABLE]

Substituting in Eq. (26), we get for the time evolution of $f(\theta,\phi,t)$ the equation

[TABLE]

with

[TABLE]

Let us consider as a far-from-equilibrium initial condition a state $f_{0}(\theta,\phi)$ that is uniform in $\phi$ over $[0,2\pi)$ and uniform in $\theta$ over a symmetric interval about $\theta=\pi/2$:

[TABLE]

with $p(\pi/2-\theta)=p(\pi/2+\theta)$. Such a state is evidently nonmagnetized, i.e., with $m_{x}=m_{y}=m_{z}=0$. It then follows that such a state is a stationary solution of the Vlasov equation (29). We now study dynamical stability of this stationary state with respect of fluctuations. The method of analysis follows closely the one pursued in Ref. [6]. To this end, we linearize the Vlasov equation (29) with respect to small fluctuations $f_{1}(\theta,\phi,t)$ by expanding $f(\theta,\phi,t)$ as

[TABLE]

with $|\lambda|\ll 1$. The linearized Vlasov equation reads

[TABLE]

where $\widetilde{m}_{x}$ and $\widetilde{m}_{y}$ are linear in $f_{1}$: $(\widetilde{m}_{x},\widetilde{m}_{y})[f_{1}]\equiv\int{\rm d}\theta{\rm d}\phi\sin\theta~{}(\sin\theta\cos\phi,\sin\theta\sin\phi)~{}f_{1}(\theta,\phi,t)$.

Now, since $f_{1}(\theta,\phi,t)$ is $2\pi$-periodic in $\phi$, we may implement the following Fourier expansion:

[TABLE]

In the long-time limit, we may expect the linearized Vlasov dynamics to be dominated by the Fourier mode of frequency $\omega$ with the smallest imaginary part, so that one effectively has $f_{1}(\theta,\phi,t)=\sum_{k}~{}g_{k}(\theta,\omega)e^{i(k\phi+\omega t)}$, yielding

[TABLE]

with $I_{\pm}=\int{\rm d}\theta\sin^{2}\theta~{}g_{\pm 1}(\theta,\omega)$. It then follows that the relevant eigenmodes of Eq. (33) are those with $k=\pm 1$. Indeed, as follows from Eq. (33), modes $k\neq\pm 1$ only oscillate in time. This fact that only the long-wavelength (i.e., small-$k$) mode perturbations are the ones that determine the stability of stationary states holds in general for long-range systems. In this regard, the reader may refer to the phenomenon of Jeans instability in a prototypical long-range system, the gravitational systems [16].

Using Eq. (35) and the aforementioned expansion of $f_{1}$ in Eq. (33), we find that the coefficients $g_{\pm 1}(\theta,\omega)$ satisfy

[TABLE]

Multiplying both sides by $\sin^{2}\theta$ and then integrating over $\theta$, we find, by using the definition of the quantity $I_{\pm}$ and the fact that $I_{\pm}\neq 0$, that

[TABLE]

Let us consider as a representative example for $f_{0}(\theta,\phi)$ the form

[TABLE]

where $\mu\equiv\sin^{2}a$ with $a>0$ and $\beta_{\rm FD}>0$ being real parameters, and $A$ is the normalization constant. In the limit $\beta_{\rm FD}\to\infty$, it is easy to see that $p(\theta)$ is a uniform distribution over the range $\theta\in[\pi/2-a,\pi/2+a]$; correspondingly, the distribution (38) becomes

[TABLE]

and is thus identical to the WB state [6]. For finite but large $\beta_{\rm FD}$, the distribution is smoothened around the boundaries at $\theta=\pi/2\pm a$. Figure 3 shows $p(\theta)$ for different values of $\beta_{\rm FD}$ and for $\mu=0.5$, which makes it evident the similarity in the form of $p(\theta)$ to the Fermi-Dirac (FD) distribution. Henceforth, we will refer to the distribution (38) as the FD state. While it is not possible to derive analytical results for the FD distribution for general $\beta_{\rm FD}$, simplifications occur for large $\beta_{\rm FD}$ when exact expressions may be derived, as detailed below.

In Appendix A, we show that for large $\beta_{\rm FD}$, we have

[TABLE]

correct to order $1/\beta_{\rm FD}^{2}$, while to same order, the energy corresponding to the state (38) is given by

[TABLE]

Next, using Eq. (38) in Eq. (37), we get to order $1/\beta_{\rm FD}^{2}$ the equation

[TABLE]

where we have

[TABLE]

while $\mu$ is to be obtained by solving Eq. (44). The latter equation gives for $n=1$ two possible values of $\mu$ given by

[TABLE]

and for $n=2$ a single value given by

[TABLE]

Using Eqs. (47) and (48) in Eq. (45) and retaining terms up to order $1/\beta_{\rm FD}^{2}$, we finally obtain a relation connecting $\omega,D,\beta_{\rm FD}$ and $\epsilon$, which is correct to same order. For $n=1$, we obtain two equations:

[TABLE]

and

[TABLE]

where we have $\widetilde{\gamma}\equiv 3\epsilon/D$ and $\widetilde{\alpha}\equiv 4D^{2}\widetilde{\gamma}$. For $n=2$, Eqs. (48) and (45) give to order $1/\beta_{\rm FD}^{2}$ a single equation:

[TABLE]

where we have $\gamma^{2}\equiv 5\epsilon/D$ and $\alpha\equiv 16D^{2}\gamma^{3}$. On physical grounds, we would want Eqs. (49), (50), and (51) to be valid for all $\omega$, including $\omega=0$. Equation (49) however gives for $\omega=0$ an inconsistent relation $8\pi^{2}D^{5}/(3\beta_{\rm FD}^{2})=0$, and hence, may be discarded.

In the limit $\beta_{\rm FD}\to\infty$, when the FD state (38) reduces to the WB state (42), Eqs. (50) and (51) reduce respectively to

[TABLE]

and

[TABLE]

Setting $\omega=0$ in Eqs. (52) and (53) yields the value $\epsilon^{\star}(D,\beta_{\rm FD}\to\infty)$ of the energy density such that these equations give only real roots for the frequency $\omega$ for $\epsilon>\epsilon^{*}$. We find that $\epsilon^{\star}$ satisfies the following equations:

[TABLE]

and

[TABLE]

It then follows that the WB state is linearly unstable under the Vlasov dynamics for energy density $\epsilon$ smaller than $\epsilon^{\star}$, and is linearly stable for $\epsilon>\epsilon^{\star}$. For $\epsilon<\epsilon^{\star}$, the perturbation $f_{1}(\theta,\phi,t)$ grows exponentially in time. Here, on setting $\omega^{2}=-\Omega^{2}$ with real $\Omega$, one gets

[TABLE]

For finite but large $\beta_{\rm FD}$, when Eqs. (50) and (51) are valid, we may expect on the basis of the above that there exists an energy threshold $\epsilon^{\star}(D,\beta_{\rm FD}~{}{\rm large})$ such that the FD state is linearly unstable and that the scaling (56) holds for energy $\epsilon<\epsilon^{\star}$, while the state is stable for energies $\epsilon>\epsilon^{\star}$. Such an $\epsilon^{\star}$ may be obtained by setting $\omega=0$ in Eqs. (50) and (51); there are more than one value of $\epsilon^{\star}$ for given $D$ and $\beta_{\rm FD}$ that one obtains in doing so, and we take for the physically meaningful $\epsilon^{\star}$ only the value that reduces to Eqs. (54) and (55) as one takes the limit $\beta_{\rm FD}\to\infty$. The result for $\epsilon^{\star}$ as a function of $D$ is shown in Fig. 1 for $n=1,2$.

5.2 Behavior for finite $N$

Equation (26) describes the time evolution in an infinite system, and here we ask: what happens when the system size $N$ is large but finite? Such a situation arises while studying the dynamics (11) numerically when obviously one has a finite $N$. In this case, as shown in Appendix B, the state of the system is described by a discrete single-spin density function $P_{d}({\bf S},t)$, which to leading order in $N$ may be expanded as

[TABLE]

For times $t\ll N$, the time evolution of $P_{0}$ is given by Eq. (26), with that for $\delta P$ given by

[TABLE]

with $\delta{\bf h}^{\rm eff}\equiv\delta{\bf h}^{\rm eff}[\delta P]$. Equivalent to Eq. (57), one may write

[TABLE]

where for times $t\ll N$, one has the time evolution of $f$ given by Eq. (29), while as was done in obtaining Eq. (29), one may show that the time evolution of $\delta f(\theta,\phi,t)\equiv\delta P({\bf S},t)$ is obtained from Eq (58) as

[TABLE]

Suppose we choose $f(\theta,\phi,0)$ to be $f_{0}(\theta,\phi)$ given by Eq. (31). It then follows from Eq. (29) that $f(\theta,\phi,t)=f(\theta,\phi,0)$, while Eq. (60) takes the same form as the linearized Vlasov equation (33):

[TABLE]

Based on our analysis in the preceding section, we may then conclude that for energies $\epsilon<\epsilon^{\star}$, when $f_{0}(\theta,\phi)$ is an unstable stationary solution of the Vlasov equation, Eq. (59) would give

[TABLE]

Using Eq. (56) that is a solution of an equation of the same form, Eq. (33), as Eq. (61), we thus obtain

[TABLE]

Thus, for $\epsilon<\epsilon^{\star}$, the relaxation time over which the magnetization acquires a value of $O(1)$ scales as $\log N$. On the other hand, for energies $\epsilon>\epsilon^{\star}$, when $f_{0}(\theta,\phi)$ is Vlasov-stationary and stable, the system would remain unmagnetized for times $t\ll N$. In this case, it is known that for longer times, the time evolution is described by (see Appendix B):

[TABLE]

Then, for $\epsilon>\epsilon^{\star}$, only for longer times of order $N$ when the dynamics (64) comes into play would there be an evolution of the initial unmagnetized state. Consequently, the state $f_{0}(\theta,\phi)$ manifests itself in a finite system as a long-lived QSS that evolves very slowly, that is, over a timescale that diverges with $N$.

5.3 Numerical results

Here, we discuss numerical results in support of our theoretical analysis of the preceding section. We present our results for two representative values of $n$, namely, $n=1,2$. In performing numerical integration of the dynamics (11), unless stated otherwise, we employ a fourth-order Runge-Kutta integration algorithm with timestep equal to $0.01$. In the numerical results that we present, data averaging has been typically over several hundreds to thousand runs of the dynamics starting from different realizations of the FD state (38).

We first discuss the results for $n=1$, for which we make the choice $D=5.0$ that yields the equilibrium critical energy $\epsilon_{c}\approx 0.2381$. Choosing as an initial condition the nonmagnetized FD state (38) with $\beta=1000$, for which one has the stability threshold $\epsilon^{\star}\approx 0.0795$ (see Fig. 1), Fig. 4(a) shows for energy $\epsilon<\epsilon^{\star}$ a fast relaxation out of the initial state on a timescale $\sim\log N$ (see Fig. 4(b)). This is consistent with the prediction based on Eq. (63), which is further validated by the collapse of the data for $\sqrt{N}m(t)$ vs. $t$ for different values of $N$ shown in Fig. 4(c); here, the growth rate $\Omega$ of $m(t)$ is obtained as the magnitude of imaginary part of the root of Eq. (50) for which the imaginary part is the largest in magnitude. Figure 5(a) shows that the relaxation observed in Fig. 4 out of the initial FD state is not to Boltzmann-Gibbs equilibrium but is to a magnetized QSS that has a lifetime that scales linearly with $N$, Fig. 5(b). Summarizing, for energy $\epsilon<\epsilon^{\star}$, relaxation of nonmagnetized FD state to Boltzmann-Gibbs equilibrium is a two-step process: in the first step, the system relaxes over a timescale $\sim\log N$ to a magnetized QSS, while in the second step, this QSS relaxes over a timescale $\sim N$ to Boltzmann-Gibbs equilibrium.

For energies $\epsilon^{\star}<\epsilon<\epsilon_{c}$, Fig. 6(a) shows that consistent with our analysis, the initial FD state appears as a nonmagnetized QSS that relaxes to Boltzmann-Gibbs equilibrium over a time which by virtue of the data presented in Fig. 6(b) may be concluded to be scaling with $N$ as $N^{3/2}$. For energies $\epsilon>\epsilon_{c}$ too is the initial FD state a stable stationary solution of the Vlasov equation, and is expected to show up as a QSS. However, here the magnetization is not the right quantity to monitor since both the FD state and Boltzmann-Gibbs equilibrium are nonmagnetized. Consequently, we choose $\langle\cos^{4}\theta\rangle=(1/N)\sum_{i=1}^{N}\cos^{4}\theta_{i}$ to monitor as a function of time (note that for $\epsilon>\epsilon_{c}$, the quantity $\langle\cos^{2}\theta\rangle$ is strictly a constant for infinite $N$, showing fluctuations about this constant value for finite $N$). Figure 7 shows that indeed the initial FD state does show up as a QSS that has a lifetime that scales quadratically with $N$. In all cases reported above and in the following when we observe an initial QSS with zero magnetization relaxing eventually to a magnetized state in equilibrium, it may be noted that due to the symmetry of the Hamiltonian (1) under spin rotation about the $z$-axis, the particular direction the magnetization chooses in equilibrium may depend on the particular realization of the QSS under study. The equilibrium magnetization vector may even have some rotation in time, and only the application of an external field may select a given orientation of the vector.

To demonstrate that the aforementioned relaxation scenario is quite generic to the model (1), we now present in Figs. 8 – 11 results for another value of $n$, namely, $n=2$. As may be observed from the figures, one has the same qualitative features of the relaxation process as that discussed above for $n=1$. Note that for energies $\epsilon>\epsilon_{c}$, one has in contrast to the $n=1$ case the quantity $\langle\cos^{4}\theta\rangle$ a constant in time and consequently one monitors $\langle\cos^{2}\theta\rangle$ as a function of time, see Fig. 11. Differences from the $n=1$ case appear in specific scalings of QSSs: the nonmagnetized QSS occurring for energies $\epsilon^{\star}<\epsilon<\epsilon_{c}$ has a lifetime scaling as $N$, while the one occurring for energies $\epsilon>\epsilon_{c}$ has a lifetime growing with $N$ as $N^{3/2}$.

6 Analysis of the stochastic dynamics (13)

6.1 Behavior in the limit $N\to\infty$

The stochastic dynamics (13) in the limit $N\to\infty$ may be studied by considering the time evolution of the single-spin distribution function $P_{0}({\bf S},t)$ derived in Appendix B as

[TABLE]

Note that the state (38), or more generally, the state (31), is not a stationary solution of Eq. (65), while, as already discussed, they both solve the energy-conserving Vlasov dynamics (26) in the stationary state. From the structure of the above equation, it follows that for times $t\ll 1/\gamma$, one may neglect the right hand side, and consequently, the time evolution is governed solely by the left hand side set to zero, which is nothing but the Vlasov equation (26). As a result, the energy is conserved for times $t\ll 1/\gamma$, and, based on the analysis presented in Section 5.1, the state (38) appears as an unstable stationary state for energies $\epsilon<\epsilon^{\star}$ and as a stable stationary state for energies $\epsilon>\epsilon^{\star}$. For times of order $1/\gamma$, we may however not neglect the right hand side of Eq. (65), and hence, we would observe the energy to be changing over times of $O(1/\gamma)$ and the state (38) to be evolving to relax to the stationary state of Eq. (65), which is nothing but the Boltzmann-Gibbs equilibrium, see Appendix B.

6.2 Behavior for finite $N$

In this case, finite-$N$ corrections need to be added to Eq. (65), and as shown in Appendix B, the time evolution is instead given by

[TABLE]

Then, based on our previous analysis, we may conclude that for a given $N$, when the noise is strong enough that $1/\gamma\ll N$, the dynamics (66) would be dominated by the first term on the right hand side. As a result, over times $t\sim 1/\gamma$, the state (38) would relax to the Boltzmann-Gibbs equilibrium state, and no size-dependent relaxation and hence QSSs should be expected. What happens in the opposite limit, that is, for $1/\gamma\gg N$? Then, over times of $O(N)$, the QSS observed for times $t\ll N$, would start evolving towards Boltzmann-Gibbs equilibrium. The relaxation would be further assisted by the effects of noise that come into effect over times of order $1/\gamma$. On the basis of the foregoing, we may expect that for a given $N$, as one tunes $\gamma$ from very small to very large values, one should see a cross-over behavior, from a size-dependent relaxation at small $\gamma$ to a size-independent one at large $\gamma$.

6.3 An alternative to dynamics (13): A Monte

Carlo dynamical scheme

An alternative way of modeling the effect of environment-induced noise on the dynamics (5) is to invoke a Monte Carlo update scheme of the spin values that guarantees that the long-time state of the system is Boltzmann-Gibbs equilibrium. In this scheme, randomly selected spins attempt to rotate by a stipulated amount (which itself could be random) with a probability that depends on the change in the energy of the system as a result of the attempted update of the state of the system [17, 18]. Specifically, to perform the Monte Carlo dynamics at temperature $T=1/\beta$, one implements the following steps [19]:

One starts with a spin configuration in the nonmagnetized FD state. 2. 2.

Next, one selects a spin at random and attempts to change its direction at random, that is, choose a value of $\theta$ uniformly in $[0,\pi]$ and a value of $\phi$ uniformly in $[0,2\pi)$ and assign these values to the spin. 3. 3.

One then computes $\Delta E$, the change in the energy of the system that this attempted change of spin direction results in. 4. 4.

If $\Delta E<0$, that is, the system energy is lowered by the change of spin direction, the change is accepted. 5. 5.

On the other hand, if the energy increases by changing the spin direction, that is, $\Delta E>0$, one computes the Boltzmann probability $p=\exp(-\beta\Delta E)$. Next, if a random number $r$ chosen uniformly in $[0,1]$ satisfies $r<p$, the change in spin direction is accepted; otherwise, the attempted change is rejected and the previous spin configuration is retained. 6. 6.

Time is measured in units of Monte Carlo steps (MCS), where one step corresponds to $N$ attempted changes in spin direction. 7. 7.

At the end of every MCS, one computes the desired physical quantities such as the magnetization. In practice, one repeats steps (ii) – (v) to obtain values as a function of time of these physical quantities averaged over a sufficient number of independent configurations.

Note that unlike the deterministic dynamics (11), the above Monte Carlo scheme does not conserve energy.

6.4 Numerical results

Here, we first discuss for $n=1$ results obtained from numerical integration of the stochastic dynamics (13) on implementing the algorithm discussed in Appendix C. For the results reported in this work, we take $\gamma=0.05$ and integration timestep equal to $10^{-3}$. Data averaging has been typically over several hundreds to thousand runs of the dynamics starting from different realizations of the FD state (38). Our aim here is to compare stochastic dynamics results with those from deterministic dynamics observed at a given energy $\epsilon$. By virtue of equivalence of microcanonical and canonical ensembles in equilibrium, we choose the temperature $T$ in the stochastic dynamics to have a value that ensures that one has in equilibrium the same value of magnetization as the one observed for the deterministic dynamics at energy $\epsilon$; this is done by using plots such as those in Fig. 2. Figure 12(a) shows that under stochastic dynamics with $1/\gamma\ll N$, the initial FD state shows a fast relaxation to Boltzmann-Gibbs equilibrium on a size-independent timescale and there is no sign of quasistationarity during the process of relaxation. Figure 12(b) shows that at the chosen value of $T$, the average energy of the system in equilibrium does coincide with the conserved energy of the deterministic dynamics, as it should due to our choice of $T$. Figure 12(c) shows for $N=10000$ the evolution of energy under the stochastic dynamics (13) for four values of the dissipation parameter $\gamma$. Scaling collapse of the data suggests relaxation of the initial state over the timescale $\sim 1/\gamma$, consistent with our analysis in Section 6.2. Relaxation on a size-independent timescale is also observed for energies $\epsilon^{\star}<\epsilon<\epsilon_{c}$ (see Fig. 6(a)) and for energies $\epsilon>\epsilon_{c}$ (see Fig. 7(a)); note that in these cases too we have $1/\gamma\ll N$. The expected cross-over in the relaxation behavior as one tunes for a fixed $N$ the value of $\gamma$ from low to high values is verified by the plot in Fig. 13. Similar results as for $n=1$ are also observed for $n=2$, see Figs. 15, 10, and 11.

Next, we show in Figs. 14 and 16 the results from Glauber Monte Carlo simulation of the system (1) for $n=1,2$. Here too one observes a fast relaxation to equilibrium over a size-independent time scale. This may be explained based on the fact that size-dependent relaxation is a feature of energy-conserving Vlasov dynamics, as is evident from the discussions in Section 5. While within the scheme (13), energy conservation is violated on the timescale $\sim 1/\gamma$, the same within the Monte Carlo scheme happens on the scale of one time unit. Hence, obviously, within the later scheme, size-dependent relaxation will not be observed anyhow, while in the former, fast relaxation requires choosing the noise to be strong enough that $1/\gamma\ll N$.

7 Conclusions

In this work, we wanted to assess the effects of stochasticity, such as those arising from the finiteness of system size or those due to interaction with the external environment, on the relaxation properties of a model long-range interacting system of classical Heisenberg spins. Under deterministic spin precessional dynamics, we showed for a wide range of energy values a slow relaxation to Boltzmann-Gibbs equilibrium over a timescale that diverges with the system size. The corresponding stochastic dynamics, modeling interaction with the environment and constructed in the spirit of (i) the stochastic Landau-Lifshitz-Gilbert equation, and (ii) the Glauber Monte Carlo dynamics, however shows a fast relaxation to equilibrium on a size-independent timescale, with no signature of quasistationarity. Our work establishes unequivocally how quasistationarity observed in deterministic dynamics of long-range systems is washed away by fluctuations induced through contact with the environment.

In the light of results on slow relaxation to equilibrium reported in this work, it would be interesting to address the issue of how the system (1) prepared either in Boltzmann-Gibbs equilibrium or in QSSs responds to an external field. One issue of particular relevance is when the field is small, and one has for short-range systems in equilibrium a linear response to the field that may be expressed in terms of fluctuation properties of the system in equilibrium. While investigation of similar fluctuation-response relations for LRI systems has been pursued in the context of particle dynamics (e.g., that of the HMF model [20, 21]) and strange scaling of fluctuations in finite-size systems has been reported [22], it would be interesting to pursue such a study for the spin model (1). Investigations in this direction have been reported in Ref. [23].

8 Acknowledgements

The work of DD is supported by UGC-NET Research Fellowship Sr. No. 2121450744, dated 29-05-2015, Ref. No. 21/12/2014(ii) EU-V. SG thanks Yoshiyuki Y. Yamaguchi for fruitful discussions, for suggesting to study the model (1) for $n=2$, and for pointing out Ref. [7]. The manuscript was finalized while S.G. was visiting the International Centre for Theoretical Physics (ICTP), Trieste during May 2019 as a Regular Associate of the Quantitative Life Sciences section, and would like to acknowledge the support and hospitality of the ICTP.

9 Appendix A: Derivation of Eqs. (43), (44),

and (45) of the main text

The normalization $A$ satisfies

[TABLE]

with $f_{\rm FD}(x)=1/(1+e^{\beta_{\rm FD}(x-\mu)})$, and where in obtaining the last equality, we have performed integration by parts. When $\beta_{\rm FD}$ is large, the first term on the right hand side of Eq. (67) drops out. In order to evaluate the second term, using the fact for large $\beta_{\rm FD}$, $\partial f_{\rm FD}(x)/\partial x=-\delta(x-\mu)$, we Taylor expand $\sqrt{x}$ about $\mu$, which on substituting in Eq. (67) gives for large $\beta_{\rm FD}$ that

[TABLE]

with

[TABLE]

Here, we have considered the limit of large $\beta_{\rm FD}$ in evaluating all the three integrals $\bar{I}_{0},\bar{I}_{1},\bar{I}_{2}$. Using Eq. (LABEL:eq:FD-integrals) in Eq. (68), we obtain for large $\beta_{\rm FD}$ the following result correct to order $1/\beta_{\rm FD}^{2}$:

[TABLE]

which is Eq. (43) of the main text.

The energy corresponding to the state (38) is given by

[TABLE]

where we have used integration by parts to arrive at the last equality. For large $\beta_{\rm FD}$, the first term on the right hand side of Eq. (71) drops out, while noting that we have $\partial f_{\rm FD}(x)/\partial x=-\delta(x-\mu)$, we evaluate the second term by Taylor expanding $x^{(2n+1)/2}$ about $x=\mu$. We finally get for large $\beta_{\rm FD}$ that

[TABLE]

with $\bar{I}_{0},\bar{I}_{1},\bar{I}_{2}$ given by Eq. (LABEL:eq:FD-integrals). Using the latter, we get for large $\beta_{\rm FD}$ that

[TABLE]

correct to order $1/\beta_{\rm FD}^{2}$. On using Eq (70), we finally get

[TABLE]

correct to order $1/\beta_{\rm FD}^{2}$. Equation (74) is Eq. (44) of the main text.

Our next job is to show how using Eq. (38) in Eq. (37) leads to Eq. (45) of the main text. It may be straightforwardly shown by using Eqs. (38) and Eq. (37) that

[TABLE]

Noting that for large $\beta_{\rm FD}$, one has $\partial f_{\rm FD}(x)/\partial x=-\delta(x-\mu)$, we may expand $g(x)$ in a Taylor series about $x=\mu$, and evaluate the right hand side. One gets to order $1/\beta_{\rm FD}^{2}$ the result

[TABLE]

Using Eq. (70), we get to order $1/\beta_{\rm FD}^{2}$ the equation

[TABLE]

where $\mu$ is to be obtained by solving Eq. (74). Equation (77) is Eq. (45) of the main text.

10 Appendix B: Derivation of Eqs. (26) and

(65) of the main text

Here, we derive Eqs. (26) and (65) of the main text. We start with the equation of motion (13):

[TABLE]

where $\boldsymbol{\eta}_{i}(t)$ is a Gaussian white noise with

[TABLE]

and

[TABLE]

In terms of components, Eq. (78) reads

[TABLE]

where we have used Einstein summation convention for repeated indices, and

[TABLE]

Let us define $F_{d}({\bf S},t)$, the discrete single-spin density function, as

[TABLE]

For a given noise realization $\{\boldsymbol{\eta}_{i}\}$, let us obtain the time evolution equation for $F_{d}$. To this end, differentiating both sides of the last equation with respect to time, using Eq. (81), and the property $a\delta(a-b)=b\delta(a-b)$, one gets

[TABLE]

with

[TABLE]

Averaging Eq. (85) over the noise statistics (79), one gets for the averaged distribution $P_{d}({\bf S},t)$ the equation [24, 25]

[TABLE]

Using Eq. (88), we have

[TABLE]

so that

[TABLE]

where we have used the fact that $\epsilon_{\alpha\lambda\beta}$ is completely antisymmetric with respect to the indices. Consequently, Eq. (89) gives

[TABLE]

Next, Eq. (88) gives $g^{\alpha\beta}g^{\lambda\beta}=(1+\gamma^{2})\epsilon_{\alpha\sigma\beta}\epsilon_{\lambda\rho\beta}S^{\sigma}S^{\rho}$, so that the right hand side of Eq. (92) now reads

[TABLE]

Consequently, Eq. (92) now reads

[TABLE]

where note that $\partial/\partial{\bf S}\cdot({\bf S}\times{\bf h}^{{\rm eff}})P_{d}=({\bf S}\times{\bf h}^{{\rm eff}})\cdot\partial P_{d}/\partial{\bf S}$.

Let us define an averaged one-spin density function $P({\bf S},t)$ as the average of $P_{d}({\bf S},t)$ with respect to a large number of initial conditions close to the same macroscopic state. To this end, we have to leading order the expansion

[TABLE]

where the deviation $\delta P$ between $P_{d}$ and $P_{0}$, which is of order $N^{0}$, is of zero average: $\langle\delta P({\bf S},t)\rangle=0$. Substituting Eq. (95) in Eq. (94), and using ${\bf h}^{\rm eff}={\bf h}^{{\rm eff},0}[P_{0}]+1/\sqrt{N}~{}\delta{\bf h}^{\rm eff}[\delta P]$, we get

[TABLE]

Using $\langle\delta P\rangle=0$, implying $\langle\delta{\bf h}^{\rm eff}\rangle=0$, then yields

[TABLE]

In the limit $N\to\infty$, when the last term on the right hand side drops out, requiring that the Boltzmann-Gibbs equilibrium state $P_{0}({\bf S})={\cal N}\exp(-\beta(-{\bf S}\cdot{\bf m}[P_{0}]+DS_{z}^{2n}))$, with ${\cal N}$ being the normalization, solves Eq. (97) in the stationary state, we must have

[TABLE]

Using

[TABLE]

we see that Eq. (98) is satisfied provided ${\cal D}(1+\gamma^{2})\beta=\gamma$. Consequently, Eq. (97) may be rewritten as

[TABLE]

In the limit $N\to\infty$, one gets Eq. (65) of the main text.

Let us consider the case of deterministic dynamics ($\gamma=0$). Then, in the limit $N\to\infty$, one gets from Eq. (102) the Vlasov equation, Eq (26), of the main text:

[TABLE]

Alternatively, Eq. (103) describes for finite $N$ the time evolution for times $t\ll N$, with that for $\delta P$ obtained from Eq. (96) as

[TABLE]

The time evolution for times of order $N$ is obtained from Eq. (102) as

[TABLE]

11 Appendix C: Numerical scheme for integrating

Eq. (13)

Here we summarize a method [25] to numerically integrate the dynamics (13) for given values of $\gamma$, $T$ and $N$. To integrate the dynamics over a time interval $[0:{\cal T}]$, we first choose a time step size $\Delta t\ll 1$, and set $t_{n}=n\Delta t$ as the $n$-th time step of the dynamics, with $n=0,1,2,...,{\cal N}_{t}$, and ${\cal N}_{t}={\cal T}/\Delta t$. One step of the update scheme from $t_{n}$ to $t_{n+1}=t_{n}+\Delta t$ involves the following updates of the dynamical variables for $i=1,2,\ldots,N$ and $\mu,\nu,\ldots=x,y,z$:

[TABLE]

where Einstein summation convention is implied. Here, $\Delta W_{i}^{\nu}$ is a Gaussian distributed random number with zero mean and variance equal to $2\gamma T$.