Equilibration in the Nosé–Hoover Isokinetic Ensemble: Effect of Inter-Particle Interactions

Abstract

We investigate the stationary and dynamic properties of the celebrated Nosé–Hoover dynamics of many-body interacting Hamiltonian systems, with an emphasis on the effect of inter-particle interactions. To this end, we consider a model system with both short- and long-range interactions. The Nosé–Hoover dynamics aim to generate the canonical equilibrium distribution of a system at a desired temperature by employing a set of time-reversible, deterministic equations of motion. A signature of canonical equilibrium is a single-particle momentum distribution that is Gaussian. We find that the equilibrium properties of the system within the Nosé–Hoover dynamics coincides with that within the canonical ensemble. Moreover, starting from out-of-equilibrium initial conditions, the average kinetic energy of the system relaxes to its target value over a size-independent timescale. However, quite surprisingly, our results indicate that under the same conditions and with only long-range interactions present in the system, the momentum distribution relaxes to its Gaussian form in equilibrium over a scale that diverges with the system size. On adding short-range interactions, the relaxation is found to occur over a timescale that has a much weaker dependence on system size. This system-size dependence of the timescale vanishes when only short-range interactions are present in the system. An implication of such an ultra-slow relaxation when only long-range interactions are present in the system is that macroscopic observables other than the average kinetic energy when estimated in the Nosé–Hoover dynamics may take an unusually long time to relax to its canonical equilibrium value. Our work underlines the crucial role that interactions play in deciding the equivalence between Nosé–Hoover and canonical equilibrium.

1. Introduction

Often, one needs in studies in nonlinear dynamics and statistical physics to investigate the dynamical properties of a many-body interacting Hamiltonian system evolving under the condition of a constant temperature. For example, one might be interested in studying the dynamical properties of the system in canonical equilibrium at a certain temperature T, with the temperature being proportional to the average kinetic energy of the system by virtue of the Theorem of Equipartition (In this work, we measure temperatures in units of the Boltzmann constant). To this end, one may devise a dynamics having a temperature $T_{\mathrm{target}}$ as a dynamical parameter that is designed to relax an initial configuration of the system to canonical equilibrium at temperature $T_{\mathrm{target}}$, and then make the choice $T_{\mathrm{target}}=T$. A common practice is to employ a Langevin dynamics, i.e., a noisy, dissipative dynamics that mimics the interaction of the system with an external heat bath at temperature $T_{\mathrm{target}}$ in terms of a deterministic frictional force and an uncorrelated, Gaussian-distributed random force added to the equation of motion [1]. In this approach, one then tunes suitably the strength of the random force such that the Langevin dynamics relaxes at long times to canonical equilibrium at temperature $T_{\mathrm{target}}$. The presence of dissipation renders the dynamics to be irreversible in time. A complementary approach to such a noisy, dissipative dynamics was pioneered by Nosé and Hoover, in which the dynamics is fully deterministic and time-reversible, while achieving the same objective of relaxing the system to canonical equilibrium at the desired temperature $T_{\mathrm{target}}$ [2,3]; for a review, see [4,5]. The time evolution under the condition of relaxation at long times to canonical equilibrium at a given temperature is said to represent isokinetic ensemble dynamics when taking place according to the Nosé–Hoover equation of motion and to represent Langevin/canonical ensemble dynamics when taking place following the Langevin equation of motion.

To illustrate in detail the distinguishing feature of the Nosé–Hoover vis-à-vis Langevin dynamics, consider an interacting N-particle system characterized by the set $\{q_j,{\pi }_j\}$ of canonical coordinates and conjugated momenta. The particles, which we take for simplicity to have the same mass m, interact with one another via the two-body interaction potential $\Phi \left(\{q_j\}\right)$. In the following, we consider $q_j$’s and ${\pi }_j$’s to be one-dimensional variables for reasons of simplicity. Our analysis, however, extends straightforwardly to higher dimensions. The Hamiltonian of the system is given by

$${\mathcal{H}}_{\mathrm{system}}=\sum _{j=1}^N\frac{{\pi }_j^2}{2m}+\Phi \left(\{q_j\}\right),$$

(1)

where the first term on the right-hand side stands for the kinetic energy of the system.

In the approach due to Langevin, the dynamical equations of the system are given by

$$\begin{array}{c}\frac{\mathrm{d}{q}_j}{\mathrm{d}t}=\frac{{\pi }_j}{m},\frac{\mathrm{d}{\pi }_j}{\mathrm{d}t}=-\gamma \frac{{\pi }_j}{m}-\frac{\partial \Phi \left(\{q_j\}\right)}{\partial q_j}+{\eta }_j\left(t\right),\hfill \end{array}$$

(2)

where t denotes time, $\gamma >0$ is the dissipation constant, while ${\eta }_j\left(t\right)$ is a Gaussian, white noise satisfying

$$\overline{{\eta }_j\left(t\right)}=0,\overline{{\eta }_j\left(t\right){\eta }_k\left(t^{\prime }\right)}=2D{\delta }_{jk}\delta (t-t^{\prime }).$$

(3)

Here, the overbars denote averaging over noise realizations, while $D>0$ characterizes the strength of the noise. The dynamics (2) are evidently not time-reversal invariant. Choosing $D=\gamma T_{\mathrm{target}}$ ensures that the dynamics (2) relaxes at long times to the canonical distribution at $T_{\mathrm{target}}$ given by [1]

$$P\left(\{q_j,{\pi }_j\}\right)\propto exp(-{\mathcal{H}}_{\mathrm{system}}/T_{\mathrm{target}}),$$

(4)

in which the kinetic energy density of the system fluctuates around the average value $T_{\mathrm{target}}/2$.

In the approach due to Nosé and Hoover, a degree of freedom s augmenting the set $\{q_j,{\pi }_j\}$ is introduced, which is taken to characterize an external heat reservoir that interacts with the system through the momenta ${\pi }_j$’s. The Hamiltonian of the combined system is given by

$$\mathcal{H}=\sum _{j=1}^N\frac{{\pi }_j^2}{2m{s}^2}+\Phi \left(\{q_j\}\right)+\frac{p_s^2}{2Q}+(N+1)T_{\mathrm{target}}lns,$$

(5)

where Q is the mass and $p_s$ is the conjugated momentum of the additional degree of freedom. The dynamics of the system is given by the following Hamilton equations of motion:

$$\begin{array}{c}\frac{\mathrm{d}{q}_j}{\mathrm{d}t}=\frac{{\pi }_j}{m{s}^2},\frac{\mathrm{d}{\pi }_j}{\mathrm{d}t}=-\frac{\partial \Phi \left(\{q_j\}\right)}{\partial q_j},\hfill \\ \\ \frac{\mathrm{d}s}{\mathrm{d}t}=\frac{p_s}{Q},\frac{\mathrm{d}{p}_s}{\mathrm{d}t}={\displaystyle \sum _{j=1}^N}\frac{{\pi }_j^2}{m{s}^3}-(N+1)\frac{T_{\mathrm{target}}}{s}.\hfill \end{array}$$

(6)

It may be easily checked that unlike dynamics (2), dynamics (6) is invariant under time reversal. In terms of new variables

$$p_j\equiv \frac{{\pi }_j}{s},\zeta \equiv \frac{p_s}{Q},$$

(7)

and rescaled time

$$\tilde{t}\equiv \frac{t}{s},$$

(8)

one obtains from the Hamilton Equations (6) the following dynamics:

$$\begin{array}{c}\frac{\mathrm{d}{q}_j}{\mathrm{d}\tilde{t}}=\frac{p_j}{m},\hfill \end{array}$$

(9)

$$\begin{array}{c}\frac{\mathrm{d}{p}_j}{\mathrm{d}\tilde{t}}=-\frac{\partial \Phi \left(\{q_j\}\right)}{\partial q_j}-\zeta p_j,\hfill \end{array}$$

(10)

$$\begin{array}{c}\frac{\mathrm{d}s}{\mathrm{d}\tilde{t}}=\zeta s,\hfill \end{array}$$

(11)

$$\begin{array}{c}\frac{\mathrm{d}\zeta }{\mathrm{d}\tilde{t}}=\frac{1}{Q}\left({\displaystyle \sum _{j=1}^N}\frac{p_j^2}{m}-(N+1)T_{\mathrm{target}}\right)=\frac{1}{{\tau }^2}\left(\frac{K\left(P\right)}{K_0}-1\right),\hfill \end{array}$$

(12)

where $K\left(P\right)\equiv {\sum }_{j=1}^N p_j^2/\left(2m\right)$ is the kinetic energy, while we have defined

$$K_0\equiv (N+1)\frac{T_{\mathrm{target}}}{2},{\tau }^2\equiv \frac{Q}{2K_0}.$$

(13)

From Equations (9)–(12), we observe that a complete description of the time evolution of the system is given in terms of Equations (9), (10), and (12), without any reference to Equation (11) for s, so that, as far as the description of the system is concerned, the variable s is an irrelevant one that may be ignored. We note in passing that a different, but closely related, Hamiltonian giving directly the Nosé-Hoover equations of motion but without any time scaling, as in Equation (8), is discussed in [6]. We will from now on drop the tilde over time in order not to overload the notation. Let us note that, in terms of the variables $p_j$’s, the Hamiltonian (5) takes the form

$$\mathcal{H}={\displaystyle \sum _{j=1}^N}\frac{p_j^2}{2m}+\Phi \left(\{q_j\}\right)+\frac{Q{\zeta }^2}{2}+(N+1)Tlns.$$

(14)

From Equation (12), we find that, in the stationary state $(\mathrm{d}\zeta /\mathrm{d}t=0$), the kinetic energy of the system equals $(N+1)T_{\mathrm{target}}/2$ (the extra factor of unity takes care of the presence of the additional degree of freedom s). For large $N\gg 1$, we then have the desired result: an ensemble of initial conditions under the evolution given by Equations (9), (10), and (12) evolves at long times to a stationary state in which the average kinetic energy density has the value $T_{\mathrm{target}}/2$. The quantity $\tau$ in Equation (12) denotes a relaxation timescale over which the kinetic energy relaxes to its target value. Beyond the average kinetic energy, it has been demonstrated by invoking the phase space continuity equation that the distribution

$$f\propto exp\left[-\left(\sum _{j=1}^N\frac{p_j^2}{2m}+\Phi \left(\{q_j\}\right)+Q{\zeta }^2/2\right)/T_{\mathrm{target}}\right]$$

(15)

is a stationary state of the Nosé–Hoover dynamics [3]. It then follows that the corresponding stationary distribution for the system variables $\{q_j,p_j\}$ is the canonical equilibrium distribution:

$$P\left(\{q_j,p_j\}\right)\propto exp\left[-\left(\sum _{j=1}^N\frac{p_j^2}{2m}+\Phi \left(\{q_j\}\right)\right)/T_{\mathrm{target}}\right],$$

(16)

normalized as $\int \left({\prod }_{j=1}^N\mathrm{d}{q}_j\mathrm{d}{p}_j\right)P\left(\{q_j,p_j\}\right)=1$. Thus, the dynamics (9)–(12) that includes the additional dynamical variable s nevertheless preseves the canonical equilibrium distribution of the system. A general formalism for constructing modified Hamiltonian dynamical systems that preserve a canonical equilibrium distribution on adding a time evolution equation for a single additional thermostat variable is discussed in [7].

Equation (16) implies that the single-particle momentum distribution $P\left(p\right)$, defined such that $P\left(p\right)\mathrm{d}p$ gives the probability that a randomly chosen particle has its momentum between p and $p+\mathrm{d}p$, is a Gaussian distribution with mean zero and width equal to $T_{\mathrm{target}}$:

$$P\left(p\right)=\frac{1}{\sqrt{2\pi m{T}_{\mathrm{target}}}}exp\left(-\frac{p^2}{2m{T}_{\mathrm{target}}}\right).$$

(17)

Consequently, the moments $\langle p^n\rangle \equiv {\int }_{-\infty }^{\infty }\mathrm{d}p{p}^{n}P\left(p\right)$, with $n=1,2,3,\dots$, satisfy $\langle p^4\rangle /{\langle p^2\rangle }^2=3$.

In the above backdrop, the principal objective of this work is to answer the question: what is the effect of inter-particle interactions on the relaxation properties of the Nosé–Hoover dynamics? More specifically, considering a system embedded in a d-dimensional space, we ask: do systems with long-range interactions, in which the inter-particle interaction decays slower than $1/r^d$, behave in a similar way to short-range systems that have the inter-particle interaction decaying faster than $1/r^d$? How does the timescale over which the phase space distribution relaxes to its canonical equilibrium form behave in the two cases, and, in particular, is there a system-size dependence in the timescale for long-range systems with respect to short-range ones? Studying these issues is particularly relevant and timely in the wake of recent surge in interest across physics in long-range interacting (LRI) systems.

LRI systems may display a notably distinct thermodynamic behavior with respect to short-range ones [8,9,10,11,12]. These systems are characterized by a two-body interaction potential $V\left(r\right)$ that decays asymptotically with inter-particle separation r as $V\left(r\right)\sim r^{-\alpha }$, with $0\le \alpha \le d$ in d spatial dimensions. The limit $\alpha \to 0$ corresponds to the case of mean–field interaction. Examples of LRI systems are self-gravitating systems, plasmas, fluid dynamical systems, and some spin systems. One of the striking dynamical features resulting from long-range interactions is the occurrence of non-equilibrium quasi-stationary states (QSSs) during relaxation of LRI systems towards equilibrium. These states have lifetimes that diverge with the number of particles constituting the system, so that, in the thermodynamic limit, the system remains trapped in QSSs and does not attain equilibrium. Only for a finite number of particles do the QSSs eventually evolve towards equilibrium. Even in equilibrium, LRI systems may exhibit features such as ensemble inequivalence and a negative heat capacity in the microcanonical ensemble that are unusual for short-range systems.

In this work, we address our aforementioned queries within the ambit of a model system comprising classical $XY$-spins occupying the sites of a one-dimensional periodic lattice and interacting via a long-range (specifically, a mean–field interaction in which every spin interacts with every other and a short-range (specifically, a nearest-neighbor interaction in which every spin interacts with its left and right neighbors) interaction. With an aim to study the equilibrium properties as well as relaxation towards equilibrium, we simulate the Nosé–Hoover dynamics of the model by integrating the corresponding equations of motion in time. A signature of canonical equilibrium is a single-particle momentum distribution that is Gaussian (see Equation (17)). We find that the equilibrium properties of our model system evolving under the Nosé–Hoover dynamics coincide with those within the canonical ensemble. As regards relaxation towards canonical equilibrium, we observe that starting from out-of-equilibrium initial conditions, the average kinetic energy of the system relaxes to its target canonical-equilibrium value over a size-independent timescale. However, quite surprisingly, our results indicate that under the same conditions and with only long-range interactions present in the system, the momentum distribution relaxes to its Gaussian form in equilibrium over a scale that diverges with the system size. On adding short-range interactions, the relaxation is found to occur over a timescale that has a much weaker dependence on system size. This system-size dependence vanishes when only short-range interactions are present in the system. An implication of such an ultra-slow relaxation when only long-range interactions are present in the system is that macroscopic observables other than the average kinetic energy when estimated in the Nosé–Hoover dynamics may take an unusually long time to relax to its canonical equilibrium value. Our work underlines the crucial role that interactions play in deciding the equivalence between Nosé–Hoover and canonical equilibrium.

The paper is organized as follows. In Section 2, we describe the model of study. In Section 3, we obtain the so-called caloric curve of the model within the canonical ensemble, which we eventually invoke in later parts of the paper to decide on the equivalence of the equilibrium properties of the Nosé–Hoover dynamics and canonical equilibrium. In Section 4, we present results from simulations of the Nosé–Hoover dynamics of the model, and discuss the implications and relevance of the results. The paper ends with conclusions in Section 5.

2. Model of Study

Our system of study comprises a one-dimensional periodic lattice of N sites. Each site of the lattice is occupied by a unit-inertia rotor characterized by its angular coordinate ${\theta }_j\in [0,2\pi )$ and the corresponding conjugated momentum $p_j$, with $j=1,2,\dots ,N$. One may also think of the rotors as representing classical $XY$-spins. Note that both of the ${\theta }_j$’s and the $p_j$’s are one-dimensional variables. There exist both a long-range (specifically, a global or a mean–field) coupling and a short-range (specifically, nearest-neighbor) coupling between the rotors. Thus, a rotor on site j interacts with strength $J/\left(2N\right)$ with rotors on all the other sites and with strength K with the rotor occupying the $(j-1)$-th and the $(j+1)$-th site. The Hamiltonian of the system is given by [13,14]

$$\begin{array}{c}\hfill \mathcal{H}={\displaystyle \sum _{j=1}^N}\frac{p_j^2}{2}+\frac{J}{2N}{\displaystyle \sum _{j,k=1}^N}\left[1-cos({\theta }_j-{\theta }_k)\right]+K{\displaystyle \sum _{j=1}^N}\left[1-cos({\theta }_{j+1}-{\theta }_j)\right];{\theta }_{N+1}\equiv {\theta }_1,p_{N+1}\equiv p_1.\end{array}$$

(18)

Note that, for $K=0$, the Hamiltonian (18) reduces to that of the widely-studied Hamiltonian mean–field (HMF) model [15], which is regarded as a paradigmatic model to study statics and dynamics of LRI systems [10]. On the other hand, for $J=0$, the model (18) reduces to a short-ranged $XY$ model in one dimension.

In the following, we take both the mean–field coupling J and the short-range coupling K to be positive, thereby modeling ferromagnetic global and nearest-neighbor couplings. Consequently, both the long-range and the short-range coupling between the rotors favor an ordered state in which all the rotor angles are equal, thereby minimizing the potential energy contribution to the total energy. Such a tendency is, however, opposed by the kinetic energy contribution whose average in equilibrium may be characterized by a temperature by invoking the Theorem of Equipartition. Noting that, for a given N, the total potential energy is bounded from above while the total kinetic energy is not, one expects the system to show in equilibrium an ordered/magnetized phase at low energies/temperatures and a disordered/unmagnetized phase at high energies/temperatures. This scenario holds even with $K=0$.

The amount of order in the system is characterized by the $XY$ magnetization

$$\mathbf{m}\equiv \frac{1}{N}\left(\sum _{j=1}^{N}cos{\theta }_j,\sum _{j=1}^{N}sin{\theta }_j\right),$$

(19)

which is a vector whose length m has the thermodynamic value in equilibrium denoted by $m^{\mathrm{eq}}$ that is nonzero in the ordered phase and zero in the disordered phase. For $K=0$, the corresponding HMF model is known to display a second-order phase transition between a high-temperature unmagnetized phase and a low-temperature magnetized phase at the critical temperature $T_c=J/2$, with the corresponding critical energy density being $u_c=3J/4$ [10]. On the other hand, invoking the Landau’s argument for the absence of any phase transition at a finite temperature in a one-dimensional model with only short-range interactions, one may conclude for $J=0$ that the corresponding short-ranged $XY$ model does not display any phase transitions, though it has been shown to have interesting dynamical effects [16]. For general $J\ne 0,K\ne 0$, when both long-range and short-range interactions are present, the model displays a second-order phase transition between an ordered and a disordered phase [13,14]. Note that all the mentioned phase transitions are continuous. Although ensemble equivalence is not guaranteed for LRI systems, it has been argued that inequivalence arises when one has a first-order phase transition in the canonical ensemble, and not when one has a second-order transition [17]. Consequently, we may regard the phase diagram of model (18) to be equivalent within microcanonical and canonical ensembles. For an explicit demonstration of ensemble equivalence for the model (18), one may refer to [14].

In the following section, we will obtain the caloric curve of model (18) that relates the equilibrium internal energy with the equilibrium temperature of the system.

3. The Caloric Curve within the Canonical Ensemble

As mentioned in the preceding section, model (18) is known to have equivalent microcanonical and canonical ensemble descriptions in equilibrium. Consequently, in obtaining the caloric curve of the model, which will be invoked to decide the equivalence between the equilibrium properties of the Nosé–Hoover dynamics and canonical equilibrium, it will suffice to restrict our analysis to the canonical ensemble description of the model.

The Langevin/canonical ensemble dynamics (2) for the model (18) comprises the set of time-evolution equations

$$\begin{array}{c}\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}=p_j,\hfill \\ \\ \frac{\mathrm{d}{p}_j}{\mathrm{d}t}=-\gamma p_j+\frac{J}{N}{\displaystyle \sum _{k=1}^N}sin({\theta }_k-{\theta }_j)+K\left[sin({\theta }_{j+1}-{\theta }_j)+sin({\theta }_{j-1}-{\theta }_j)\right]+{\eta }_j\left(t\right),\hfill \end{array}$$

(20)

with the properties of the noise ${\eta }_j\left(t\right)$ given by Equation (3) with $D=\gamma T$. Within the microcanonical ensemble description of the system, the time evolution of the variables $\{{\theta }_j,p_j\}$ is given by Hamilton equations obtained from Equation (20) by setting $\gamma$ to zero. The Nosé–Hoover dynamics of the variables $\{{\theta }_j,p_j\}$ is obtained from Equations (9) and (10) as

$$\begin{array}{c}\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}=p_j,\hfill \\ \\ \frac{\mathrm{d}{p}_j}{\mathrm{d}t}=\frac{J}{N}{\displaystyle \sum _{k=1}^N}sin({\theta }_k-{\theta }_j)+K\left[sin({\theta }_{j+1}-{\theta }_j)+sin({\theta }_{j-1}-{\theta }_j)\right]-\zeta p_j,\hfill \end{array}$$

(21)

where the time evolution of the variable $\zeta$ is given by Equation (12).

In order to derive the desired caloric curve of model (18) within the canonical ensemble, we start with the canonical partition function of the system at temperature T given by $Z_N\equiv \int \left({\prod }_{j=1}^N\mathrm{d}{\theta }_j\mathrm{d}{p}_j\right)exp[-\beta \mathcal{H}\left(\{{\theta }_j,p_j\}\right)]$, with $\beta \equiv 1/T$. Using Equation (18), we get

$$Z_N={\left(\frac{2\pi }{\beta }\right)}^{N/2}{e}^{-\beta JN/2-\beta KN}\int \left({\prod }_{j=1}^N\mathrm{d}{\theta }_j\right)exp\left[\frac{\beta J}{2N}\left\{{\left({\sum }_{j=1}^{N}cos{\theta }_j\right)}^2+{\left({\sum }_{j=1}^{N}sin{\theta }_j\right)}^2\right\}+\beta K{\sum }_{j=1}^{N}cos({\theta }_{j+1}-{\theta }_j)\right].$$

(22)

Using the Hubbard–Stratonovich transformation $exp\left(a{x}^2\right)=1/\left(\sqrt{4\pi a}\right){\int }_{-\infty }^{\infty }\mathrm{d}zexp\left(-\frac{z^2}{4a}+zx\right),a>0$ in Equation (22), we obtain

$$\begin{array}{c}{Z}_N={\left(\frac{2\pi }{\beta }\right)}^{N/2}{e}^{-\beta JN/2-\beta KN}\frac{N\beta J}{2\pi }{\int }_{-\infty }^{\infty }\mathrm{d}{z}_1{\int }_{-\infty }^{\infty }\mathrm{d}{z}_2\int \left({\displaystyle \prod _{j=1}^N}\mathrm{d}{\theta }_j\right)exp[-\frac{N\beta J}{2}(z_1^2+z_2^2)\hfill \\ +\beta J{z}_1{\displaystyle \sum _{j=1}^N}cos{\theta }_j+\beta J{z}_2{\displaystyle \sum _{j=1}^N}sin{\theta }_j+\beta K{\displaystyle \sum _{j=1}^N}cos({\theta }_{j+1}-{\theta }_j)].\hfill \end{array}$$

(23)

Writing $z_1=zcos\varphi ,z_2=zsin\varphi$, with real $z={(z_1^2+z_2^2)}^{1/2}>0$ and $\varphi \in [0,2\pi )$ given by $\varphi ={tan}^{-1}(z_2/z_1)$, we get

$$\begin{array}{c}{Z}_N={\left(\frac{2\pi }{\beta }\right)}^{N/2}{e}^{-\beta JN/2-\beta KN}\frac{N\beta J}{2\pi }{\int }_0^{2\pi }\mathrm{d}\varphi {\int }_0^{\infty }\mathrm{d}zz\int \left({\displaystyle \prod _{j=1}^N}\mathrm{d}{\theta }_j\right)exp[-\frac{N\beta J}{2}{z}^2\hfill \\ +\beta Jz{\displaystyle \sum _{j=1}^N}cos({\theta }_j-\varphi )+\beta K{\displaystyle \sum _{j=1}^N}cos({\theta }_{j+1}-{\theta }_j)].\hfill \end{array}$$

(24)

In view of the invariance of the Hamiltonian (18) under rotation by an equal amount of all the ${\theta }_j$’s, we get [18]

$$Z_N={\left(\frac{2\pi }{\beta }\right)}^{N/2}{e}^{-\beta JN/2-\beta KN}N\beta J{\int }_0^{\infty }\mathrm{d}zz\int \left({\prod }_{j=1}^N\mathrm{d}{\theta }_j\right)exp\left[-\frac{N\beta J}{2}{z}^2+\beta Jz{\sum }_{j=1}^{N}cos{\theta }_j+\beta K{\sum }_{j=1}^{N}cos({\theta }_{j+1}-{\theta }_j)\right].$$

(25)

In order to proceed further, we consider separately the cases $K=0$ and $K\ne 0$ in the following.

3.1. $K=0$

For $K=0$, Equation (25) yields

$$Z_N={\left(\frac{2\pi }{\beta }\right)}^{N/2}N\beta J{\int }_0^{\infty }\mathrm{d}zzexp\left[-N\left\{\frac{\beta J}{2}(1+z^2)-ln\left({\int }_0^{2\pi }\mathrm{d}\theta exp(\beta Jzcos\theta )\right)\right\}\right].$$

(26)

In the thermodynamic limit, $Z_N$ may be approximated by invoking the saddle-point method to perform the integration in z on the right-hand side; one gets

$$Z_N={\left(\frac{2\pi }{\beta }\right)}^{N/2}N\beta J{z}_{s}exp\left[-N\left\{\frac{\beta J}{2}(1+z_s^2)-ln\left({\int }_0^{2\pi }\mathrm{d}\theta exp(\beta J{z}_{s}cos\theta )\right)\right\}\right],$$

(27)

where the saddle-point value $z_s$ solves the equation

$$z_s=\frac{I_1\left(\beta J{z}_s\right)}{I_0\left(\beta J{z}_s\right)},$$

(28)

with $I_n\left(x\right)=(1/\left(2\pi \right)){\int }_0^{2\pi }\mathrm{d}\theta exp(xcos\theta )cos\left(n\theta \right)$ being the modified Bessel function of first kind and of order n. It may be shown by following the arguments given in [18] that $z_s$ is nothing but the stationary magnetization $m^{\mathrm{eq}}$. Equation (28) has a trivial solution $m^{\mathrm{eq}}=0$ valid at all temperatures, while a non-zero solution exists for $\beta \ge {\beta }_c=2/J$ [10]. In fact, the system shows a continuous transition, from a magnetized phase ($m^{\mathrm{eq}}\ne 0$) at low temperatures to an unmagnetized phase ($m^{\mathrm{eq}}=0$) at high temperatures at the critical temperature $T_c=J/2$ [10].

In the thermodynamic limit, the internal energy density of the system $u=-{lim}_{N\to \infty }(1/N)\mathrm{d}ln{Z}_N/\mathrm{d}\beta$ is obtained by using Equations (27) and (28) as

$$u=\frac{1}{2\beta }+\frac{J}{2}\left(1-{\left(m^{\mathrm{eq}}\right)}^2\right);m^{\mathrm{eq}}=\frac{I_1\left(\beta J{m}^{\mathrm{eq}}\right)}{I_0\left(\beta J{m}^{\mathrm{eq}}\right)},$$

(29)

yielding the critical energy density

$$u_c=\frac{3J}{4}.$$

(30)

Equation (29) gives the caloric curve of the model (18) at canonical equilibrium for $J\ne 0,K=0$.

3.2. $K\ne 0$

For $K\ne 0$, Equation (25) gives

$$\begin{array}{c}{Z}_N={\left(\frac{2\pi }{\beta }\right)}^{N/2}N\beta J{\int }_0^{\infty }\mathrm{d}zzexp\left[-\frac{N\beta J}{2}(1+z^2)-\beta KN\right]{\mathcal{Z}}_N,\hfill \end{array}$$

(31)

$$\begin{array}{c}{\mathcal{Z}}_N\equiv \int \left({\displaystyle \prod _{j=1}^N}\mathrm{d}{\theta }_j\right)exp\left[\beta Jz{\displaystyle \sum _{j=1}^N}cos{\theta }_j+\beta K{\displaystyle \sum _{j=1}^N}cos({\theta }_{j+1}-{\theta }_j)\right],\hfill \end{array}$$

(32)

where we may identify the factor ${\mathcal{Z}}_N$ with the canonical partition function of a $1d$ periodic chain of N interacting angle-only rotors, where a rotor on each site interacts with strength K with the rotor on the left nearest-neighbor and the right nearest-neighbor site, and also with an external field of strength $Jz$ along the x direction.

One may evaluate ${\mathcal{Z}}_N$ by rewriting it in terms of a transfer operator $T(\theta ,{\theta }^{\prime })$ as

$$\begin{array}{c}{\mathcal{Z}}_N=\int \left({\displaystyle \prod _{j=1}^N}\mathrm{d}{\theta }_j\right)\mathcal{T}({\theta }_1,{\theta }_2)\mathcal{T}({\theta }_2,{\theta }_3)\dots \mathcal{T}({\theta }_N,{\theta }_1),\hfill \end{array}$$

(33)

$$\begin{array}{c}\mathcal{T}({\theta }_j,{\theta }_{j+1})\equiv exp\left[\beta Jz\left\{\frac{cos{\theta }_j+cos{\theta }_{j+1}}{2}\right\}+\beta Kcos({\theta }_{j+1}-{\theta }_j)\right].\hfill \end{array}$$

(34)

Let $\{{\lambda }_m\}$ denote the set of eigenvalues of the transfer operator $\mathcal{T}(\theta ,{\theta }^{\prime })$. In other words, denoting the eigenfunctions of $\mathcal{T}(\theta ,{\theta }^{\prime })$ as $f_m\left(\theta \right)$, we have $\int \mathrm{d}{\theta }^{\prime }\mathcal{T}(\theta ,{\theta }^{\prime })f_m\left({\theta }^{\prime }\right)={\lambda }_m{f}_m\left(\theta \right)$. In terms of $\{{\lambda }_m\}$, we obtain

$${\mathcal{Z}}_N=\sum _m{\left[{\lambda }_m\left(\beta Jz,\beta K\right)\right]}^N.$$

(35)

For large N, the sum in Equation (35) is dominated by the largest eigenvalue ${\lambda }_{\max }={\lambda }_{\max }\left(\beta Jz,\beta K\right)$, yielding

$${\mathcal{Z}}_N={\lambda }_{\max }^N.$$

(36)

Substituting Equation (36) in Equation (31), and approximating the integral on the right-hand side of the latter by the saddle-point method, one gets

$$Z_N={\left(\frac{2\pi }{\beta }\right)}^{N/2}N\beta J{z}_{s}exp\left[-N\left\{\frac{\beta J}{2}(1+z_s^2)+\beta K-ln{\lambda }_{\max }\left(\beta J{z}_s,\beta K\right)\right\}\right],$$

(37)

where $z_s$ solves the saddle-point equation $z_s\equiv {sup}_z\tilde{\varphi }(\beta ,z)$, with $\tilde{\varphi }(\beta ,z)$ being the free-energy function:

$$-\tilde{\varphi }(\beta ,z)\equiv -\frac{1}{2}ln\beta -\frac{\beta J}{2}(1+z^2)-\beta K+ln{\lambda }_{\max }\left(\beta Jz,\beta K\right).$$

(38)

The saddle-point equation may thus be written as

$$z_s={\frac{\partial ln{\lambda }_{\max }\left(\beta Jz,\beta K\right)}{\partial \left(\beta Jz\right)}|}_{z=z_s}.$$

(39)

Equation (37) gives the dimensionless free energy per rotor, $\varphi \left(\beta \right)\equiv -{lim}_{N\to \infty }(ln{Z}_N)/N$, as $-\varphi \left(\beta \right)={sup}_z\left[-\tilde{\varphi }(\beta ,z)\right]$, where we have suppressed the dependence of $\varphi \left(\beta \right)$ on K. We thus have

$$-\varphi \left(\beta \right)\equiv -\frac{1}{2}ln\beta -\frac{\beta J}{2}(1+z_s^2)-\beta K+ln{\lambda }_{\max }\left(\beta J{z}_s,\beta K\right).$$

(40)

Note that the free energy at a given temperature has a definite value given by Equation (40), and is obtained by substituting the saddle-point solution $z_s$ into the expression for the free-energy function $\tilde{\varphi }(\beta ,z)$.

In the thermodynamic limit, the internal energy density of the system $u=-{lim}_{N\to \infty }(1/N)\mathrm{d}ln{Z}_N/\mathrm{d}\beta$ is obtained as

$$u=\frac{1}{2\beta }+\frac{J}{2}(1+z_s^2)+\beta J{z}_s\frac{\mathrm{d}{z}_s}{\mathrm{d}\beta }+K-{\frac{\mathrm{d}ln{\lambda }_{\max }(\beta Jz,\beta K)}{\mathrm{d}\beta }|}_{z=z_s}.$$

(41)

Using Equation (39), and the fact that, as for $K=0$, the quantity $z_s$ is nothing but the stationary magnetization $m^{\mathrm{eq}}$, we get

$$u=\frac{1}{2\beta }+\frac{J}{2}\left(1-{\left(m^{\mathrm{eq}}\right)}^2\right)+\beta J{m}^{\mathrm{eq}}\frac{\mathrm{d}{m}^{\mathrm{eq}}}{\mathrm{d}\beta }+K-K\frac{\partial ln{\lambda }_{\max }(\beta J{m}^{\mathrm{eq}},\beta K)}{\partial \left(\beta K\right)},$$

(42)

with $m^{\mathrm{eq}}$ satisfying

$$m^{\mathrm{eq}}={\frac{\partial ln{\lambda }_{\max }\left(\beta Jz,\beta K\right)}{\partial \left(\beta Jz\right)}|}_{z=m^{\mathrm{eq}}}.$$

(43)

To proceed, we need to find ${\lambda }_{\max }(\beta Jz,\beta K)$. We consider separately the cases $J=0$ and $J\ne 0$.

3.2.1. $J=0$

In this case, it may be easily checked that the eigenvalues of $\mathcal{T}$ are given by $2\pi I_m\left(\beta K\right)$ with the corresponding eigenvector given by plane waves $exp\left(iq\theta \right)/\sqrt{2\pi }$ [14]. Using $I_0\left(x\right)>I_1\left(x\right)>I_2\left(x\right)\dots$, we conclude that ${\lambda }_{\max }(0,\beta K)=I_0\left(\beta K\right)$. Equation (43) then yields $m^{\mathrm{eq}}=0$, while Equation (42) gives

$$u=\frac{1}{2\beta }+K\left(1-\frac{I_1\left(\beta K\right)}{I_0\left(\beta K\right)}\right),$$

(44)

where we have used the result $\mathrm{d}{I}_0\left(x\right)/\mathrm{d}x=I_1\left(x\right)$. Equation (44) is the desired caloric curve of the model (18) within the canonical ensemble for $J=0,K\ne 0$.

3.2.2. $J\ne 0$

In this case, not knowing the analytic forms of the eigenvalues of $\mathcal{T}$, we resort to a numerical scheme to estimate the largest eigenvalue ${\lambda }_{\max }(\beta Jz,\beta K)$. To this end, we discretize the angles over the interval $[0,2\pi )$ as ${\theta }_j^{\left(a_j\right)}=a_j\Delta \theta$, with $a_j=1,2,\dots ,P$ and $\Delta \theta =2\pi /P$ for any large positive integer P (we choose $P=30$). The operator $\mathcal{T}(\theta ,{\theta }^{\prime })$ then takes the form of a matrix of size $P\times P$, whose largest eigenvalue may be estimated numerically by employing the so-called power method [19] (A FORTRAN90 library that implements the power method and is distributed under the GNU Lesser General Public License (GPL) is available at [20]). Noting that $\mathcal{T}(\theta ,{\theta }^{\prime })$ is a finite-dimensional real square matrix with positive entries, the application of the Perron–Frobenius theorem implies the existence of its largest eigenvalue that is real and non-degenerate. At given values of $T,K,J,z$, once ${\lambda }_{\max }(\beta Jz,\beta K)$ has been estimated numerically, we compute the free-energy function $\tilde{\varphi }(\beta ,z)$ as a function of z by using Equation (38). We then find numerically the value of z at which the computed free-energy function attains its minimum value. As discussed above, this minimizer is the equilibrium magnetization of the system at the given values of $T,K,J$. In order to obtain the caloric curve, one has to estimate numerically the derivative $\partial ln{\lambda }_{\max }(\beta J{m}^{\mathrm{eq}},\beta K)/\partial \left(\beta K\right)$, and then use Equation (42).

4. Results and Discussion

In this section, we discuss the results on equilibrium as well as relaxation properties of model (18) obtained by performing numerical integration of the Nosé–Hoover equations of motion (21). The numerical integration involved using a fourth-order Runge–Kutta method with timestep $\mathrm{d}t=0.01$.

4.1. Results in Equilibrium

Here, we discuss the Nosé–Hoover equilibrium properties for model (18). The initial condition corresponds to the ${\theta }_j$’s independently and uniformly distributed in $[0,2\pi )$ and the $p_j$’s independently sampled from a Gaussian distribution with zero mean and width equal to 0.5. The initial value of the parameter $\zeta$ is 2.0, while we have taken $\tau =0.01$. In Figure 1, we consider the case when only long-range interactions are present in system ($J=1.0,K=0.0$). Figure 1a shows for $T_{\mathrm{target}}=2.5$ that the average kinetic energy relaxes at long times to the value $T_{\mathrm{target}}/2$, as desired. Figure 1b shows for the same value of $T_{\mathrm{target}}$ that the average internal energy has the same value in the stationary state as the one in canonical equilibrium given by Equation (29); Figure 1c shows the single-particle momentum distribution $P\left(p\right)$ in the stationary state. We observe that $P\left(p\right)$ has the correct canonical-equilibrium form of a Gaussian distribution, which further corroborates the property of the Nosé–Hoover dynamics that the canonical distribution (16) is a stationary state of the dynamics. Figure 1d shows for a range of values of the temperature $T=T_{\mathrm{target}}$ that the caloric curve obtained within the Nosé–Hoover dynamics in equilibrium coincides with that within the canonical ensemble given by Equation (29). Figure 1a–c refer to the system size $N=128$, while Figure 1d refers to two system sizes, namely, $N=128$ and $N=512$. The aforementioned observed properties of the Nosé–Hoover dynamics have been checked to hold for (i) the case when only short-range interactions are present in the system (see Figure 2 that corresponds to $J=0.0,K=1.0$), in which case the caloric curve within the canonical ensemble is given by Equation (44), and (ii) when both long- and short-range interactions are present in the system (data not shown; see, however, Figure 3c).

4.2. Results out of Equilibrium

Here, we discuss the relaxation properties of the Nosé–Hoover dynamics for model (18). The initial condition corresponds to the so-called water-bag distribution that has both $\theta$ and p uniformly distributed over given intervals [10]. We consider ${\theta }_j$’s to be independently and uniformly distributed in $[0,2\pi )$ and the $p_j$’s to be independently and uniformly distributed in $[-\sqrt{1.5},\sqrt{1.5}]$. The initial value of the parameter $\zeta$ is 2.0, while we have taken $\tau =1.0$.

Let us start with a discussion of the results in Figure 4 that corresponds to the case when only long-range interactions are present in the system (18). In Figure 4a, we see that, for four different system sizes, the average kinetic energy density relaxes at long times to the target value $T_{\mathrm{target}}/2$ over a timescale that does not depend on the system size. A Gaussian distribution for the momentum, expected in canonical equilibrium, is characterized by a value 3 of the ratio $\langle p^4\rangle /{\langle p^2\rangle }^2$ (see Equation (17)). We see in Figure 4b that, in contrast to Figure 4a, this ratio, however, relaxes to the canonical equilibrium value over a time that depends on the system size, and which grows with increase of N. Figure 4c shows that the long-time magnetization value reached by the Nosé–Hoover dynamics coincides with the canonical equilibrium value for all system sizes. On the basis of these results, we conclude that, with only long-range interactions in system (18), only the second moment of the momentum distribution relaxes to its canonical equilibrium value over a size-independent timescale, while higher moments (and consequently, the whole of the momentum distribution) relax to their canonical equilibrium values over a time that grows with the system size. The latter fact is demonstrated in Figure 4d that shows for $N=512$ the time evolution of the single-particle momentum distribution.

The feature of a size-independent timescale for the relaxation of the average kinetic energy density to its canonical equilibrium value, observed in the case of purely long-range interactions in model (18), also holds on adding short-range interactions to the model and when the latter are the only interactions present in the system (see Figure 3a and Figure 5a). Moreover, in all cases, the long-time value of the magnetization matches with its canonical equilibrium value (see Figure 3c and Figure 5c). The most significant difference in the relaxation properties that is observed on adding short-range interactions may be inferred by comparing Figure 3b and Figure 4b: the very strong size-dependence observed in the relaxation of the ratio $\langle p^4\rangle /{\langle p^2\rangle }^2$ to its canonical equilibrium value gets substantially weakened on adding short-range interactions with coupling strength as low as $K=0.1$ compared to the value of the long-range coupling constant $J=1.0$. Similar inference may be drawn from a comparison of Figure 3d and Figure 4d. This observation has an immediate and an important implication: additional short-range interactions speed up the relaxation of the momentum distribution towards canonical equilibrium. The aforementioned system-size dependence vanishes on turning off long-range interactions, so that the only remnant interactions in the system are the short-range ones (see Figure 5b,d).

5. Conclusions

In this paper, we investigated the relaxation properties of the Nosé–Hoover dynamics of many-body interacting Hamiltonian systems, with an emphasis on the effect of inter-particle interactions. The dynamics aim to generate the canonical equilibrium distribution of a system at the desired temperature by employing time-reversible, deterministic dynamics. To pursue our study, we considered a representative model comprising N classical $XY$-spins occupying the sites of a one-dimensional periodic lattice. The spins interact with one another via both a long-range interaction, modelled as a mean–field interaction in which every spin interacts with every other, and a short-range one, modelled as a nearest-neighbor interaction in which every spin interacts with its left and right neighboring spins. We studied the Nosé–Hoover dynamics of the model through N-body integration of the corresponding equations of motion. Canonical equilibrium is characterized by a momentum distribution that is Gaussian. We found that the equilibrium properties of our model system evolving according to Nosé–Hoover dynamics are in excellent agreement with exact analytic results for the equilibrium properties derived within the canonical ensemble. Moreover, while starting from out-of-equilibrium initial conditions, the average kinetic energy of the system relaxes to its target value over a size-independent timescale. However, quite unexpectedly, we found that under the same conditions and with only long-range interactions present in the system, the momentum distribution relaxes to its Gaussian form in equilibrium over a scale that grows with N. The N-dependence gets weaker on adding short-range interactions, and vanishes when the latter are the only inter-particle interactions present in the system.

Viewed from the perspective of LRI systems, the slow relaxation observed within the Nosé–Hoover dynamics allows for drawing an analogy with a similar slow relaxation observed within the microcanonical dynamics of isolated LRI systems, a phenomenon that leads to the occurrence of nonequilibrium quasistationary states (QSSs) that have lifetimes diverging with the system size [10,21]. Within a kinetic theory approach, the QSSs are understood as stable, stationary solutions of the so-called Vlasov equation that governs the time evolution of the single-particle phase space distribution. The Vlasov equation is obtained as the first equation of the Bogoliubov–Born–Green–Yvon–Kirkwood (BBGKY) hierarchy by neglecting the correlation between particle trajectories, with corrections that decrease with an increase of N. For large but finite N, the eventual relaxation of QSSs towards equilibrium is understood as arising due to these finite-N corrections, the so-called collisional terms, to the Vlasov equation. In models in which the momentum variables are one-dimensional, it has been shown by analyzing the behavior of the dominant collisional term that Vlasov-stable phase-space distributions that are homogeneous in the coordinates evolve on times much larger than N, thereby leading for the distributions to characterize QSSs that have lifetimes diverging with N [8,11,22]. In light of the foregoing discussions, it is evidently pertinent and of immediate interest to invoke a kinetic theory approach and investigate in the context of the Nosé–Hoover dynamics of long-range systems whether additional short-range interactions play the role of collisional dynamics that speed up the relaxation of the system towards canonical equilibrium. Work in this direction is in progress and will be reported elsewhere.

The agreement reported in this paper in the value of the average kinetic energy computed in canonical equilibrium and within the Nosé–Hoover dynamics is reminiscent of a similar agreement in the large-system limit between ensemble and time averages predicted by Khinchin for the so-called sum-functions, that is, functions such as the kinetic energy that are sums of single-particle contributions [23]. The result was obtained for rarefied gases, which was later observed to also hold for systems with short-range interactions [24,25]. Our work hints at the validity of such a result even for long-range systems, as is evident from the agreement in the value of the average kinetic energy computed within the Nosé–Hoover dynamics and in canonical equilibrium (see Figure 4a). This point warrants a more detailed investigation that will be left for future studies.