Statistical mechanics of systems with long-range interactions and negative absolute temperature

TL;DR

This paper investigates a Hamiltonian model with long-range interactions, revealing negative temperature regimes where ensemble equivalence breaks down, supported by analytical calculations and molecular dynamics simulations.

Contribution

It demonstrates the existence of negative-temperature regimes and ensemble nonequivalence in a long-range interacting system through combined analytical and numerical methods.

Findings

Negative temperature regime identified

Ensemble nonequivalence confirmed

Fluctuation-Dissipation Theorem validated at negative temperatures

Abstract

A Hamiltonian model living in a bounded phase space and with long-range interactions is studied. It is shown, by analytical computations, that there exists an energy interval in which the microcanonical entropy is a decreasing convex function of the total energy, meaning that ensemble equivalence is violated in a negative-temperature regime. The equilibrium properties of the model are then investigated by molecular dynamics simulations: first, the caloric curve is reconstructed for the microcanonical ensemble and compared to the analytical prediction, and a generalized Maxwell-Boltzmann distribution for the momenta is observed; then, the nonequivalence between the microcanonical and canonical descriptions is explicitly shown. Moreover, the validity of Fluctuation-Dissipation Theorem is verified through a numerical study, also at negative temperature and in the region where the two ensembles are nonequivalent.

Equations

$$S(E)=\ln {\int }_{\Omega }d\mathbf{X}\delta (H(\mathbf{X})-E)$$

$$H(\mathbf{X})$$

$$v(m)$$

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

$${\rho }_j(p|E)=\frac{1}{\omega (E)}{\int }_{\Omega }d\mathbf{X}\delta (H(\mathbf{X})-E)\delta (p-p_j)$$

$$\omega (E)={\int }_{\Omega }d\mathbf{X}\delta (H(\mathbf{X})-E).$$

$${\rho }_j(p|E)$$

$$={\omega }^{\prime }(E-k(p))/\omega (E)$$

$$=\exp \left[S^{\prime }(E-k(p))-S(E)\right]$$

$${\rho }_j(p|E)\propto \exp \left[-\beta k(p)\right]$$

$$⟨\cos p_j⟩=-I_1(\beta )/I_0(\beta )$$

$$⟨A(\mathbf{X})⟩\simeq \frac{1}{\mathcal{T}}{\int }_0^{\mathcal{T}}A(\mathbf{X}(t))dt$$

$$X_j\to X_j^{\prime }=X_j+\delta X_j$$

$$\delta A(t)=A({\mathbf{X}}^{\prime }(t))-A(\mathbf{X}(t)),$$

$$\langle\delta A\rangle(t)=-\sum_{j}\left\langle A(\mathbf{X}(t))\frac{\partial\ln\rho(\mathbf{X})}{\partial X_{j}}\bigg{|}_{\begin{subarray}{c}t=0\end{subarray}}\right\rangle\delta X_{j}(0)$$

$$\frac{⟨\delta A⟩(t)}{\delta p_j(0)}=\beta ⟨\sin (p_j(t))\cdot \sin (p_j(0))⟩$$

$$H(\mathbf{X})=\bar{H}({\mu }_1(\mathbf{X}),...,{\mu }_n(\mathbf{X})).$$

$$S(E)=\ln \int \prod _{j}d{\bar{\mu }}_j\delta (\bar{H}({\bar{\mu }}_1,...,{\bar{\mu }}_n)-E)\exp [N\bar{s}({\bar{\mu }}_1,...,{\bar{\mu }}_n)]$$

$$\bar{s}({\bar{\mu }}_1,...,{\bar{\mu }}_n)=\frac{1}{N}\ln {\int }_{\Omega }d\mathbf{X}\delta ({\mu }_j(\mathbf{X})-{\bar{\mu }}_j).$$

$$\bar{s}=\underset{{\lambda }_1,...,{\lambda }_n}{\inf }\left\{\sum _j{\lambda }_j{\bar{\mu }}_j+\frac{\ln Z({\lambda }_1,...,{\lambda }_n)}{N}\right\}$$

$$Z({\lambda }_1,...,{\lambda }_n)={\int }_{\Omega }{e}^{-N{\sum }_j{\lambda }_j{\mu }_j(\mathbf{X})}d\mathbf{X}.$$

$$\kappa =\frac{1}{N}\sum _{i=1}^N\left(1-\cos (p_i)\right)$$

$$\frac{H}{N}=\kappa -\frac{J}{2}(m_x^2+m_y^2)-\frac{K}{4}(m_x^4+2m_x^2{m}_y^2+m_y^4)$$

$$\bar{s}=\underset{{\lambda }_{\kappa }}{\inf }$$

$$+\underset{{\lambda }_x,{\lambda }_y}{\inf }\{{\lambda }_x{m}_x+{\lambda }_y{m}_y+\ln I_0\left(\sqrt{{\lambda }_x^2+{\lambda }_y^2}\right)\}$$

$$=(\kappa$$

$$+\ln (4{\pi }^2)-m{B}_{inv}(m)+\ln I_0(B_{inv}(m))$$

$$\bar{H}(\kappa ,m_x,m_y)=E.$$

$$B_{inv}(\tilde{m})=(J\tilde{m}+K{\tilde{m}}^3)B_{inv}\left(1-\frac{E}{N}-\frac{J}{2}{\tilde{m}}^2-\frac{K}{4}{\tilde{m}}^4\right).$$

$$S(E)=N\bar{s}(\tilde{\kappa },\tilde{m}).$$

$$F(\beta )=\underset{\kappa ,m}{\inf }\left\{H(\kappa ,m)-\frac{N}{\beta }\bar{s}(\kappa ,m)\right\}.$$

$$\dot{\mathbf{X}}=\left(\begin{array}{c}\dot{q}\\ \dot{p}\end{array}\right)=\left(\begin{array}{c}{\partial }_{p}H\\ -{\partial }_{q}H-D\beta {\partial }_{p}H+\sqrt{2D}\xi (t)\end{array}\right)$$

$$q(t_0+h)=q_0+{\partial }_{p}K(p_0)\frac{h}{2}+O(h^2)$$

Full text

I Introduction

It is well known that in physical systems with long-range interactions the equivalence of statistical ensembles can fail Dauxois et al. (2002); Campa et al. (2009), meaning that there exist equilibrium states described by the microcanonical probability density function (p.d.f.) that do not correspond to any state described by the canonical one: in other words, the average of a macroscopic observable can give different results, at the same temperature $T$, if the system is isolated or closed. This phenomenon can be observed in a large variety of physical contexts including self-gravitating systems, perfect fluids in two dimensions and spin models with mean field interactions Padmanabhan (1990); Lynden-Bell et al. (1968); Lynden-Bell (1999); Chavanis (2002, 2006); Ramírez-Hernández et al. (2008); Hovhannisyan et al. (2017).

For long-range interacting systems, inequivalence of statistical ensembles is due to the lack of additivity of the total energy $E$, which can result in a change of concavity for the entropy $S(E)$ (i.e. in negative specific heat): since the Legendre-Fenchel transform that relates the free energy $F(T)$ to $S(E)$ is not invertible in this case, there is no one-to-one correspondence between the microcanonical and the canonical description Touchette (2009); Patelli and Ruffo (2014). In particular, it can be shown that the caloric curve $T$ vs $E$ in the canonical ensemble can be obtained from the microcanonical analogue through the Maxwell construction, i.e. by replacing the free energy with its convex envelope. The situation is somehow reminiscent of that arising in van der Waals equation for non ideal gases, obtained by a mean field approach, where the pressure is a non-monotonic function of the volume: Maxwell construction was indeed introduced in this context, in order to describe the coexistence of two different phases in the “unphysical” region at negative compressibility Ma (1985); Huang (1988). Let us stress, however, that in this case there is no real ensemble inequivalence: as discussed in Lebowitz and Penrose (1966), the negative compressibility is a mere consequence of the homogeneity assumption, and it does not appear if this hypothesis is relaxed; the equation of state can be derived through a rigorous limit procedure on the interaction potential (the so-called van der Waals scaling), and no negative-compressibility region is found. Ensemble inequivalence arises instead if the long-range interaction part of the potential is mean-field Kiessling and Percus (1995).

The aim of this paper is to investigate the statistical features of a system with ensemble inequivalence in a negative temperature regime: in particular we show that the momenta, at equilibrium, follow a generalized Maxwell-Boltzmann distribution, and that a Fluctuation-Dissipation theorem holds.

Let us recall that a system with Hamiltonian $H(\mathbf{X})$, where $\mathbf{X}$ is a point in the phase space $\Omega$, is said to have a negative absolute temperature if its Boltzmann entropy

[TABLE]

is a decreasing function of the energy $E$ in a certain energy interval (here and in the following, we set the Boltzmann constant $k_{B}$ equal to 1): this essentially means that the volume of the accessible phase-space region shrinks when increasing the energy of the system Ramsey (1956). It is easy to understand that Hamiltonians in which the kinetic energy has the usual form, quadratic in the momenta, cannot achieve negative temperature: for the ideal gas $S(E)\propto\ln E$, and additional degrees of freedom cannot change the sign of $dS(E)/dE$.

The possibility of negative absolute temperature states in Hamiltonian systems is well established. An important example is already described in the Onsager’s work on the statistical hydrodynamics of point vortices Onsager (1949); this seminal contribution was the starting point for the investigation of the dynamical and statistical properties of hydrodynamics systems Marchioro and Pulvirenti (1994). These models are interesting both for applications (e.g. in plasma physics) and for their statistical mechanics properties. Let us note, for instance, that an example of a system with ensemble inequivalence at negative temperature is already studied in Smith and O’Neil (1990): here the considered system is a guiding-center model for a plasma in a cylindric domain, equivalent to a 2-D point vortices system; it is shown that the equilibrium solution of the mean-field Vlasov equation, in a high-energy range at negative temperature, presents ensemble inequivalence. Rigorous results on the statistical mechanics of point vortices in bounded domains have been proven Kiessling ; Caglioti et al. (1992, 1995), relating the inequivalence of statistical ensembles to the non-uniqueness of the solution, at high energy and for some kind of bounded domains in two dimensions, of the mean-field equation describing the system. Recently, the problem of equivalence of statistical ensembles for systems of point voritices on a spherical surface has been also addressed Kiessling and Wang (2012).

The relevance of the presence of negative temperature has been investigated in another important Hamiltonian model, the Discretized Nonlinear Schrödinger Equation: their emergence in out-of-equilibrium conditions have been recently investigated Iubini et al. (2012, 2013, 2017). Experimental evidence of negative temperature equilibrium states, on the other hand, dates back to the studies on nuclear spins by Purcell, Pound and Ramsey in the 50’s Purcell and Pound (1951); Ramsey (1956). In more recent times, experiments have shown that negative temperature is also present in systems of cold atoms Braun et al. (2013).

The structure of this paper is the following: in Section II we describe the Hamiltonian model we are interested in, and define the mechanical observable that we will use to measure the temperature of the system in numerical simulations; in Section III we present our results at equilibrium, showing that it is possible to have ensemble inequivalence at negative temperature; Section IV is devoted to the study of response at negative temperature (in particular, the validity of the Fluctuation-Dissipation Theorem is verified); in Section V we briefly sketch our conclusions. Appendix A contains the calculations that lead to our analytical results, while Appendix B illustrate the details of the algorithm we have used for our simulations.

II The model

Let us consider a Hamiltonian system consisting of $N$ degrees of freedom described by conjugated coordinates $\{p_{i},\theta_{i}\}$, $i=1,...,N$. Positions $\{\theta_{i}\}$ and momenta $\{p_{i}\}$ are both angular variables, meaning that they are constrained in the interval $(-\pi,\pi]$ with periodic boundary conditions. The total Hamiltonian reads:

[TABLE]

where $m$ is the modulus of the magnetization vector defined as

[TABLE]

$J$ and $K$ are parameters that can assume, in general, both positive and negative values; the additive constant in Eqn. (2b) is actually unessential for the dynamics: in the following we will choose it in such a way that the minimal energy achievable by the system is zero.

First, let us comment on the unusual dependence of Hamiltonian (2a) on the momenta. As discussed in the Introduction, negative temperature equilibrium states cannot be observed in mechanical systems ruled by quadratic kinetic energy; a more promising class of models is obtained by considering “modified” kinetic terms, defined as periodic functions of the momenta: in this way the phase-space is bounded and we can expect the presence of an energy range with decreasing entropy. This choice has been considered in previous works about negative temperature Cerino et al. (2015); Baldovin et al. (2017, 2018), and its consistence has been checked in different contexts. Let us just notice here that in the small energy limit such terms reduce to the usual quadratic ones.

The potential term defined by Eqn. (2b) is the one that characterize the so-called “Generalized Hamiltonian Mean Field” (GHMF) model de Buyl et al. (2005) describing the mean field interactions between magnetized rotators. The GHMF is an extended version of the Hamiltonian Mean Field model Antoni and Ruffo (1995) that includes also a quartic dependence on the magnetization; among other interesting properties, this system is paradigmatic for the study of inequivalence between canonical ad microcanonical ensembles and has been extensively studied in past years Bouchet et al. (2008); Campa et al. (2009); Baldovin (2018).

Let us note that the bound phase space (due to the specific shape of the “kinetic” terms) and the long-range interactions are two features that our model shares with two dimensional vortex systems in bounded domains Onsager (1949); Smith and O’Neil (1990); Marchioro and Pulvirenti (1994). Therefore it is quite natural to wonder about the possibility of ensemble inequivalence in the negative temperature regime.

Hamiltonian $H$ depends on the angular positions through the magnetization modulus $m$, which is maximal ($m=1$) when all the “modified” rotators are parallel and vanishes if their angular positions are homogeneously distributed in $(-\pi,\pi]$; the $N$ factor in front of the potential is needed in order to ensure the extensivity of the system (Kac’s prescription). Let us stress that for the system (2a) it is possible to write down an analytical computation of the equilibrium properties using large deviation approaches Patelli and Ruffo (2014): the procedure is quite similar to the one presented in Campa et al. (2009), and it is explicitly carried out in Appendix A.

In order to compare the theoretical prediction for the equilibrium states of the system to the results of numerical simulations, as a preliminary step we need to define a proper mechanical observable for the inverse temperature $\beta$. In Hamiltonian systems with the usual kinetic energy $\sum_{i}p_{i}^{2}/2m$ one would exploit the equipartition theorem Ma (1985); Huang (1988) and measure $T=1/\beta$ as the (temporal) average $\langle p_{i}^{2}/m\rangle$ – assuming that ergodic hypothesis holds. Since model (2a) is not quadratic in the momenta we cannot follow this approach.

The microcanonical p.d.f. for the momentum of the $j$-th particle is given by

[TABLE]

where $\mathbf{X}$ is a point in the phase space $\Omega$ and $\omega(E)$ represents the density of states (d.o.s.)

[TABLE]

P.d.f. (4) can be evaluated as discussed in Cerino et al. (2015): let us split the phase space into the product $\Omega=\Gamma_{j}\times\Omega^{\prime}$, where $\Gamma_{j}\equiv(-\pi,\pi]$ is the interval to whom variable $p_{j}$ belongs; we get

[TABLE]

where $k(p)=1-\cos(p)$ is the “kinetic” contribution of the single particle, so that $H^{\prime}(\mathbf{X}^{\prime})=H(\mathbf{X})-k(p_{j})$ does not depend on $p_{j}$. Here $\omega^{\prime}(E)$ and $S^{\prime}(E)$ are the d.o.s. and the microcanonical entropy for the Hamiltonian $H^{\prime}(\mathbf{X}^{\prime})$, and we have used the definition of entropy (1). Expanding equation (6) in $k(p)$ ($E$ is of the order of $N$ and $k(p)$ is a bounded function, so that $|k(p)|\ll E$) we finally get

[TABLE]

(assuming that $dS^{\prime}/dE=dS/dE\equiv\beta$ in the thermodynamical limit). In order to measure the temperature of the system we can study the observable

[TABLE]

where $I_{n}(x)$ is the $n-$th modified Bessel function of the first kind. It can be easily proved that $B(x)=I_{1}(x)/I_{0}(x)$ is an odd function, thus it can be inverted in order to find $\beta$.

Let us also stress that, since the r.h.s. of the above equation can be positive, $\beta$ can assume, in principle, also negative values.

III Equilibrium properties

Molecular dynamics simulations have been performed in order to study the functional dependence of the inverse temperature $\beta$ on the energy $E$ of the system. We used a second-order symplectic Velocity Verlet-like integrator (see Appendix B, case $D=0$), choosing the integration step $\Delta t$ so as to observe total-energy relative fluctuations $\Delta E/E\simeq O(10^{-5})$. Since we are interested in equilibrium properties, the total integration time $\mathcal{T}$ has always been chosen to be much longer than the typical characteristic times of the system (which can be shown to be $O(10)$: see e.g. the plots in Section IV). In order to fix initial conditions with the desired value of the total energy, some care has to be devoted to the fact that our “kinetic” terms are bounded; first, one has to choose the angular positions in such a way that $0<E-Nv(m)<2N$, e.g. by repeated uniform extractions over a suitable interval; once the above constraint is satisfied, momenta can be randomly extracted in such a way that the sum of the kinetic terms amounts to the residual energy. We also include a small external potential $v_{ext}(\mathbf{X})=\cos(\theta_{1})$, acting only on the first particle, in order to break the angular symmetry of the magnetization $\mathbf{m}$. All interesting observables are computed as temporal averages:

[TABLE]

In Fig. 1 we show the theoretical caloric curve $\beta(E/N)$ for a choice of $J$ and $K$ that leads to inequivalence of statistical ensembles (namely, to the existence of microcanonical equilibrium states with negative specific heat that have no canonical counterpart). In this case both $J$ and $K$ are positive, so that the interacting potential has a minimum for $m=1$ (“ferromagnetic” limit).

Let us notice that the values of $\beta$ measured in numerical simulations through the observable (8) show a very good agreement with the analytical prediction.

The figure also shows the difference between the system we are considering and the corresponding GHMF model with the same parameters. At low energies the behavior of the two systems is quite similar, due to the fact that the kinetic terms are equal up to order $O(p_{i}^{3})$ when energies are small. The scenario changes for large values of $E/N$: in this case, the “bounded phase-space” version of the GHMF model can reach equilibrium states with $\beta\equiv dS/dE<0$, while the original GHMF model behaves like an ideal gas whose (positive) $\beta$ asymptotically tends to zero (infinite temperature).

Fig. 2 shows a different example, in which both $J$ and $K$ are chosen to be negative; in this case $v(m)$ is minimum at $m=0$ (“antiferromagnetic” limit).

As can be seen in the figure, in a certain energy range negative microcanonical specific heat and negative values of $\beta$ coexist – meaning that microcanonical and canonical ensembles are inequivalent for some negative temperatures. The microcanonical specific heat of the corresponding GHMF model is, instead, always positive, and the statistical ensembles are always equivalent. In the insets of Fig. 2 we also show the single-particle momentum distribution in a couple of cases (one at positive and one at negative $\beta$): the empirical distribution and the analytical expression for the p.d.f. given by Eqn. (7) are in excellent agreement.

Let us now focus on the energy range in which statistical ensembles are inequivalent in this case. Fig. 3 shows the detail of this region. In panel (a), the microcanonical caloric curve is compared to the results of numerical simulations. The dashed line is the Maxwell’s construction, which individuates the transition temperature $1/\beta^{*}$ in the canonical ensemble Touchette (2009): all the microcanonical equilibrium states with energies corresponding to the Maxwell’s construction (red circles in the figure) are metastable or unstable Campa et al. (2009) and have no equivalent counterpart in the canonical ensemble. The microcanonical caloric curve also shows a first order phase transition, which can be understood by looking at panel (b): the systems goes from a regime with $m=0$ to a magnetized phase with $m>0$.

Ensemble inequivalence can be checked through numerical simulations. We start from an initial condition whose energy corresponds to an unstable microcanonical state, and we let the system evolve through a molecular dynamics simulation for a long time $\mathcal{T}_{0}$. So far the system is isolated, and we can measure the inverse temperature $\beta$ of the corresponding microcanonical state by studying the above discussed observable. The final configuration is then taken as the initial condition for two different evolution protocols: on the one hand we let the system evolve further trough a symplectic deterministic integrator, as before, for an additional time $\mathcal{T}_{1}$; on the other hand we can perform a stochastic evolution of the system, for the same time, in order to mimic the presence of a thermal reservoir and let the system reach a canonical equilibrium state, in which $\beta$ is fixed to the value measured before. The details about the integrator are discussed in Appendix B. We fix the diffusivity constant $D$ that appears in the generalized Langevin Equation (26) in such a way that $1/|\beta|D$, the typical time-scale of the stochastic dynamics, is comparable to the characteristic times of the system.

The results of our simulations are shown in Fig.4 for two choices of the initial energy that lead to unstable microcanonical states, namely $E/N=1.84$ (panels (a) and (b)) and $E/N=1.94$ (panels (c) and (d)). In the left part of the figure, the single-particle momentum distributions, microcanonical and canonical, are displayed: in both cases the plots are almost identical, meaning that our simulation protocol is able to mimic a canonical ensemble with the same temperature of the starting isolated state. On the right panels we show the corresponding distributions for the angular positions: in both cases we can see a clear difference between the magnetization of the original microcanonical state and that of the canonical state at the same temperature.

IV Response theory

In this section we study the response of the system to a small perturbation. In particular we want to check the validity of the well-known Fluctuation-Dissipation Relation (FDR) (see Ref.Marconi et al. (2008)) that describes the temporal evolution of a mechanical observable, on average, after the system has been perturbed.

Suppose that at time $t=0$ the microscopic state of the system $\mathbf{X}(0)$ is instantaneously changed into a new configuration $\mathbf{X}^{\prime}(0)$ by modifying the canonical coordinate (position or momentum) $X_{j}$:

[TABLE]

where $|\delta X_{j}|$ is small on the typical scales of the system. We can follow the evolution of both systems and measure the mechanical observable $A(\mathbf{X})$ in both cases. Defining

[TABLE]

the FDR ensures that

[TABLE]

where $\rho(\mathbf{X})$ is the equilibrium p.d.f. and $\langle\cdot\rangle$ represents the corresponding average.

In our case, if we choose $A(\mathbf{X})=\sin p_{j}$, Eqn. (12) can be rewritten as:

[TABLE]

where a perturbation on the $j$-th momentum has been considered.

The above relation can be simply checked by numerical simulations. The strategy consists in preparing the system in a certain configuration $\mathbf{X}$, perturbing the momentum of a randomly chosen particle by a quantity $\delta p(0)$ in order to obtain a new configuration $\mathbf{X}^{\prime}$ and integrating the dynamics of both systems for a time $\tau$. During the evolution one can measure $A(\mathbf{X})$ and $A(\mathbf{X}^{\prime})$. The procedure is then repeated $M\gg 1$ times, starting from the final configuration of the previous evolution of $\mathbf{X}$: since the integrator is symplectic, in this way we are sure to consider initial conditions at the same energy, and we can average over all the measured evolutions.

In Fig. 5 we show the results of our numerical simulations for two different choices of $E$. The agreement between the measured values of $\delta A(t)$ and Eqn. (13) is clear.

Let us notice that the FDR can be used, in general, to give an operative definition of the inverse temperature $\beta$ of the system through the measure of responses and autocorrelation functions Marconi et al. (2008): the presence of negative temperature and ensemble inequivalence does not hinder this possibility.

V Conclusions

In this paper we studied a Hamiltonian system with long-range interactions and modified ”kinetic” terms (i.e., the usual $p^{2}/2$ has been replaced by $1-\cos p$ for all particles). In such a model analytical computations and numerical simulations show the existence of an energy range with absolute negative temperature. Comparing the results obtained with microcanonical and canonical ensembles we show that such ensembles are non-equivalent. Even in presence of such non standard features, the validity of Fluctuation-Dissipation Theorem is verified through a numerical study of the response of the system to a small perturbation.

In spite of the fact that there has been a certain confusion in the literature Romero-Rochín (2013); Dunkel and Hilbert (2014); Hilbert et al. (2014), systems with negative temperature do not show any kind of peculiarity and there is no evidence of behaviors in disagreement with the basic features established in statical mechanics; this has been clearly discussed in several recent papers Frenkel and Warren (2015); Buonsante et al. (2016); Puglisi et al. (2017); Swendsen (2018). We briefly list a series of facts showing that systems with negative absolute temperature do not show any anomaly Puglisi et al. (2017):

1. a)

the statistical properties of the energy fluctuations are ruled by the inverse absolute temperature $\beta$ defined by the microcanonical entropy; 2. b)

the p.d.f. of the momentum is given by a generalized Maxwell-Boltzmann distribution where $\beta$ appears; 3. c)

$\beta$ is a measurable quantity, without the appearance of any inconsistency; 4. d)

in presence of slow variables (e.g. a “heavy intruder” in a chain) it is possible to introduce, in a consistent way, a Langevin equation Baldovin et al. (2018); 5. e)

possibility of non equivalence of the statistical ensembles Smith and O’Neil (1990); Caglioti et al. (1995); 6. f)

validity of the Fluctuation-Dissipation Theorem.

Acknowledgements.

We thank A. Campa for helpful discussions and useful comments.

Appendix A Large deviation analysis

In this Appendix we sketch the large-deviation strategy used to determine the equilibrium properties of model (2a). Our argument follows quite closely the usual derivation that can be done for the GHMF model (see, e.g. Campa et al. (2009)), with an ad-hoc handling of the “kinetic” terms. Let us just recall that large deviation approaches can be exploited to study the thermodynamic behavior of a system composed of $N$ particles if $n\ll N$ mean quantities $\mu_{j}(\mathbf{X})=\sum_{i=1}^{N}g_{j}(q_{i},p_{i})/N$, $j=1,...,n$, of the phase-space exist such that, for a suitable choice of $\bar{H}$, one has

[TABLE]

In such a case the following relation holds:

[TABLE]

where $\bar{s}$ is the so-called “entropy of macrostates”:

[TABLE]

Furthermore, large deviation theory ensures that the entropy of macrostates can be computed, in the $N\gg 1$ limit, as

[TABLE]

with

[TABLE]

In the present case, defining

[TABLE]

we can write Eqn. (2a) (apart from unessential constants) as

[TABLE]

where $m_{x}$ and $m_{y}$ are the components of the magnetization vector $\mathbf{m}$ (see Eqn. (3)). The entropy of macrostates can be shown to be equal to:

[TABLE]

where $I_{n}(x)$ is the n-th Modified Bessel function of the first kind and $B_{inv}(x)$ is the inverse function of $I_{1}(x)/I_{0}(x)$. Let us notice that $\bar{s}$ depends actually only on $\kappa$ and $m=\sqrt{m_{x}^{2}+m_{y}^{2}}$. Now we can evaluate the microcanonical entropy $S(E)$ by using its relation with $\bar{s}(\kappa,m)$. The r.h.s. of Eqn. (15) in the $N\gg 1$ limit can be estimated through a constrained extremal problem, namely as the supremum of $\bar{s}(\kappa,m)$ with the condition

[TABLE]

It can be shown that such problem is fulfilled by the value $\tilde{m}$ of the magnetization that verifies the following condition:

[TABLE]

Once $\tilde{m}$ is known (Eqn. (23) can be solved by numerical methods) the corresponding value of $\kappa$ can be found from Eqn. (22) and we finally have

[TABLE]

In a similar way we can derive the free energy of the system:

[TABLE]

Appendix B Numerical integrator for the Langevin dynamics

In this Appendix we briefly discuss the stochastic Verlet-like integration scheme used to simulate Hamiltonian systems with generalized “kinetic” terms, subjected to a thermal bath at fixed $\beta$. We are therefore interested in integration schemes for the “generalized” Langevin equation

[TABLE]

where $\xi(t)$ is a white noise with unitary variance and $D$ is a positive constant Baldovin et al. (2018).

Since in the limit $D=0$ the dynamics reduces to that of an isolated Hamiltonian system, we would like to find an algorithm that is symplectic in this limit. This problem has been addressed in Melchionna (2007) for the case of quadratic kinetic energy: in the following we will apply the same reasoning to a wider class of Hamiltonian systems.

It can be shown Honeycutt (1992) that the integration of Eqn. (26) over a time-step $h$ leads to errors of order $h^{3/2}$ if the deterministic and the stochastic parts are evolved independently. In order to improve the stability of the algorithm one can alternate the integrations of positions and momenta, as in the usual Position Verlet. In the present case, the integration of the positions can be trivially written as

[TABLE]

where $(q_{0},p_{0})$ is the state for the considered degree of freedom at time $t_{0}$ and $K(p)$ is the kinetic term.

On the other hand, in order to evolve the momenta we have to deal with the following partial differential equation:

[TABLE]

where $F_{0}=-\partial_{q}V(t_{0})$). If $K(p)=p^{2}/2m$, such equation can be solved exactly; in the general case we can use the approximation

[TABLE]

with an error of order $h^{2}$ on $p(h)$. Substituting the above expression into Eqn. (28) we get an equation of the form

[TABLE]

where A and B are coefficients that are constant during the integration step. Solving the above equation, one gets the following Position Verlet-like integration scheme:

[TABLE]

where

[TABLE]

An equivalent Velocity Verlet-like algorithm can be found by inverting the integration order. Let us notice that in the deterministic limit $D=0$ such second-order integrators are also symplectic.