(Non-equilibrium) thermodynamics of integrable models: The Generalized Gibbs Ensemble description of the classical Neumann model

Interest in the long-time dynamics of quantum isolated systems has continuously grown since the celebrated quantum Newton's cradle experiment [1], which proved that a quenched one-dimensional Bose gas does not reach standard thermal equilibrium. Soon after, a Generalized Gibbs Ensemble (GGE) was proposed to describe typical observables in the steady state of systems with an extensive number of conserved quantities, say $I_\mu$ with $\mu=1,\dots, N$, and N the number of degrees of freedom [2,3]. The pertinence of such density matrix was studied in a myriad of different one-dimensional quantum cases where the necessity to focus on quasi-local observables and constants of motion was stressed [4–7].

Although most studies of quenches of isolated systems have focused on quantum systems, non-ergodic dynamics are not specifically quantum: classical integrable systems [8–10] are not expected to reach equilibrium as dictated by conventional statistical mechanics either. One can then ask whether a GGE description could apply to their long-term evolution as well and, if so, under which conditions. Yuzbashyan argued that the Generalized Microcanonical Ensemble (GME), in which the values of all independent constants of motion are fixed, is exact for classical integrable systems [11]. More explicitly, this means long-time averages of the Newtonian dynamics, $\overline{A} = \lim_{\tau\to\infty} \lim_{t_{\rm st} \gg t_0}\tau^{-1} \int_{t_{\rm st}}^{t_{\rm st} + \tau} \textrm{d}t' \, A(t')$, and statistical averages calculated with the flat GME measure, $\langle A \rangle_{\textrm{GGE}} =\sum_{\textrm{conf}} A \;\prod_{\mu=1}^N \delta(I_\mu - I_\mu(0))/ \sum_{\rm conf} \prod_{\mu=1}^N \delta(I_\mu - I_\mu(0))$, of typical observables coincide (the sums run over all possible configurations). Here, $I_\mu(0)$ are the values of the constants of motion at the initial time t = 0, and constrain the manifold in phase space in which the dynamics occurs. However, this does not ensure that a canonical GGE, of the form $e^{-\sum_\mu \gamma_\mu I_\mu}/Z_{\textrm{GGE}}$, could be derived from such a GME, especially in long-range interacting systems in which the notion of a subsystem is not straightforward and the additivity of the conserved quantities is not justified [12,13]. It is therefore of paramount importance to explicitly construct the GGE of a classical integrable interacting model (beyond the independent harmonic-oscillator case) and put to the test its main statement, that in the stationary limit ¹ the long-time average, $\overline{A}$, and the phase space average, $\langle A \rangle_{\textrm{GGE}} =\sum_{\textrm{conf}} A e^{-\sum_\mu \gamma_\mu I_\mu(\textrm{conf})}/Z_{\textrm{GGE}}$, coincide (for any not explicitly time-dependent and non-pathological observable A). The parameters $\gamma_\mu$ are Lagrange multipliers fixed by requiring that the phase space averages of the N constants of motion, $\langle I_\mu\rangle_{\textrm{GGE}}$, equal their values at the initial conditions, $I_\mu(0)$. For an early discussion of the GGE for a classical system see [14], and for an approach based on generalized hydrodynamics see [15,16].

Our goal here is to exhibit one such non-trivial classical model, the Neumann Model. We used a mixed analytic-numerical treatment to prove that in the thermodynamic limit, N → ∞, taken before the long-time limit, $t\gg t_{\textrm{st}}$, with $t_{\textrm{st}}$ a characteristic time-scale, the model reaches a stationary state which satisfies the extended ergodic hypothesis with a GGE measure in which the $I_\mu$ are integrals of motion in involution (with quartic dependencies on the phase space variables) [17]. In so doing, we elucidated the dynamic phase diagram and we evidenced condensation phenomena and macroscopic fluctuations that should be of importance, as we explain, in quenches of Bose-Einstein condensates.

The Neumann Model (NM) is a non-trivial classical integrable system [18] that describes the motion of a particle on a sphere embedded in an N-dimensional space, $S_{N-1}$, under fully anisotropic harmonic forces. The Hamiltonian is

$$\begin{equation} H_{\textrm{quad}} = H_{\mathrm{kin}}+H_{\mathrm{pot}}=\frac{1}{2m} \sum_\mu p_\mu^2 - \frac{1}{2} \sum_\mu \lambda_\mu s_\mu^2, \end{equation} \tag{ 1 }$$

with $s_\mu$, $\mu=1, \ldots, N$, the coordinates of the position vector, $p_\mu$ the corresponding momentum components, m the mass, and $-\lambda_\mu$ the spring constants. The primary and secondary spherical constraints are

$$\begin{eqnarray} &&C_1 \equiv \sum_\mu s_\mu^2 = N, \qquad C_2 \equiv \sum_\mu s_\mu p_\mu = 0. \end{eqnarray} \tag{ 2 }$$

The equations of motion, subject to the constraints (2) can be derived with the Poisson-Dirac method and read

$$\begin{equation} \dot p_\mu = (\lambda_\mu - z)s_\mu. \end{equation} \tag{ 3 }$$

The "Lagrange multiplier" z is given by

$$\begin{equation} z=\frac{1}{N}\sum_{\mu}( p^2_{\mu}/m + \lambda_{\mu}s^2_{\mu} ), \end{equation} \tag{ 4 }$$

makes the μ-modes interact, and ensures the validity of C₁ and C₂. For any initial condition satisfying these constraints, the dynamics conserve them, the quadratic Hamiltonian, $H_{\textrm{quad}}$, as well as the N Uhlenbeck integrals of motion in involution [19–22],

$$\begin{equation} I_\mu = s_\mu^2 + \frac{1}{mN} \sum_{\nu(\neq\mu)} \frac{s_\mu^2 p_\nu^2 + s_\nu^2 p_\mu^2 -2 s_\mu p_\mu s_\nu p_\nu}{\lambda_\nu-\lambda_\mu}. \end{equation} \tag{ 5 }$$

The latter verify $\sum_\mu I_\mu = C_1$, and $\sum_\mu \lambda_\mu I_\mu = - 2 H_{\textrm{pot}} - 2H_{\textrm{kin}} \, C_1/N + 1/(mN) \, C_2^2$, which equals $-2H_{\textrm{quad}}$ thanks to the constraints in eq. (2).

We are interested in developing a statistical description of the NM dynamics. This can make sense only in the limit N → ∞ taken before any long-time limit. In this setting one can expect the fluctuations of z to be suppressed, and

$$\begin{equation} z(t)\mapsto \langle z(t)\rangle_{i.c.}, \end{equation} \tag{ 6 }$$

where we made the time dependency of z explicit. The angular brackets represent an average over any distribution of initial conditions satisfying the constraints C₁ and C₂. The replacement in eq. (6) is a kind of mean-field approximation that softens the strict conditions on C₁ and C₂ and imposes them on average, $\langle C_1\rangle_{i.c.} = N$ and $\langle C_2\rangle_{i.c.} = 0$. In this way, we introduce a variation of the original model that we call the Soft Neumann Model (SNM). With the replacement of the phase space function z by its average, the equations of motion (3) are reduced to those of a system of harmonic oscillators with frequencies

$$\begin{equation} \omega^2_{\mu}(t)=\frac{\langle z(t)\rangle_{i.c.}-\lambda_{\mu}}{m}. \end{equation} \tag{ 7 }$$

As already said, this model satisfies both constraints only on average over initial conditions. Similarly, it has no strictly conserved quantities but

$$\begin{equation} \langle H_{\textrm{quad}}\rangle_{i.c.} = \textrm{const} \quad \textrm{and} \quad \langle I_\mu\rangle_{i.c.} = \textrm{const}\quad \forall\mu, \end{equation} \tag{ 8 }$$

are also conserved on average. The conditions under which the NM and SNM are equivalent will be analyzed below.

Quadratic potential energies combined with a global spherical constraint as the one in eq. (1) are common in statistical physics. Depending on the choice of the spring constants $\lambda_\mu$ one finds, e.g., the celebrated spherical ferromagnet [23,24] or the so-called p = 2 disordered spherical model [25,26]. Problems of particles embedded in large-dimensional spherical spaces and subject to random potentials can also be of this kind. For convenience, and to make a closer connection with the physics of disordered systems, we order the λ's such that $\lambda_1 < \lambda_2 < \cdots< \lambda_N$ and in the large N numerical applications we take them to be distributed by a Wigner semi-circle law on the interval $[-2J,2J]$. In this way, they can be thought of as the eigenvalues of a two-body interaction matrix with zero-mean Gaussian-distributed entries that couple the coordinates in a different basis (e.g., real spins with a global spherical constraint). The fact that they take values within a real interval with an edge ensures that the total energy is bounded from below. We expect to find similar results for other choices of the distribution of λ's with an edge in the N → ∞ limit.

In most studies of quantum quenches, the initial condition is taken to be the ground state of a Hamiltonian which is suddenly modified. However, equilibrium finite-temperature initial states [27–30] are more relevant to describe, for instance, experiments in ultracold Bose gases [31]. Along these lines, we draw the initial conditions from a proper Gibbs-Boltzmann equilibrium measure

$$\begin{equation} \rho_{i.c.}=\frac{1}{Z(T')} \, \exp\left[-\beta' H_{\textrm{quad}}- \frac{\beta' z_{\textrm{eq}}}{2} \, \Big(\sum_\mu s_\mu^2-N\Big)\right], \end{equation} \tag{ 9 }$$

where $z_{\textrm{eq}}$ is the equilibrium value of the Lagrange multiplier enforcing the spherical constraint (also on average) at inverse temperature $\beta'=1/T'$ with $k_B=1$. $Z(T')$ is the canonical partition function, and $H_{\textrm{quad}}$ is given in eq. (1) with spring constants $\lambda_\mu^{(0)}$ in the interval $[-2J_0, 2J_0]$. The equilibrium $z_{\textrm{eq}}$ takes the values

$$\begin{eqnarray} z_{\textrm{eq}} = \left\{ \begin{array}{l@{\qquad}l} T' + J_0^2/T', & \text{for} ~~T'\geq J_0,\\ 2J_0, & \text{for}~~T'\leq J_0 . \end{array}\right. \end{eqnarray} \tag{ 10 }$$

Depending on $T'/J_0$ being larger or smaller than one, the initial conditions belong to an extended phase in which the variances of all modes are ${\mathcal O}(1)$, or to a condensed phase in which the averaged squared N-th mode, $\langle s^2_N\rangle_{i.c.}$, scales as ${\mathcal O}(N)$ [25]. Two scenarii for the condensation phenomenon are possible: a mixed two-pure-state measure with the possibility of symmetry breaking induced by a vanishing pinning field, or a Gaussian measure centered at zero with diverging dispersion [24,32,33]. In the magnetic interpretation, $T'=J_0$ is a critical point between a disordered and a magnetically ordered phase. The analogy with Bose-Einstein Condensation (BEC) was already reckoned in [25] with $s_N^{2}$ playing the role of the ground state density. The transition, in this case, is towards the condensed phase.

We drive the system out of equilibrium by performing a sudden interaction quench in which we rescale all spring constants, $\lambda_\mu^{(0)} \mapsto \lambda_\mu$, with the same factor $J/J_0$ that controls the amount of energy injected ($J/J_0<1$) or extracted ($J/J_0<1$). This procedure mimics the quenches performed in isolated quantum systems [4–6]. Right after the instantaneous quench, the initial kinetic energy of all modes is ${\mathcal O}(1)$ and the averaged Uhlenbeck constants are ${\mathcal O}(1)$ for $T'>J_0$ while $\langle I_N\rangle_{i.c.} = {\mathcal O}(N)$ for $T'<J_0$. Each $\langle I_\mu\rangle_{i.c.}$ is a function of $\lambda_\mu/J$ and the adimensional parameters $T'/J_0$ and $J/J_0$ that can be easily calculated:

$$\begin{equation} \langle I_\mu\rangle_{i.c.} = \frac{T'^2}{J_0 J} \frac{(J_0+J) J_0/T' - \lambda_\mu^{(0)}}{z_{\rm eq}(T',J_0) - \lambda_\mu^{(0)}} \end{equation} \tag{ 11 }$$

for $T' \geq J_0$ or $T'\leq J_0$ and $\mu \neq N$, and

$$\begin{equation} \langle I_N\rangle_{i.c.} = (1-T'/J_0) (1-T'/J) N + {\mathcal O}(1) \end{equation} \tag{ 12 }$$

for $T'\leq J_0$.

Insight into the long-time dynamics of the SNM was gained in [14,34]. In these papers we studied the Schwinger-Dyson equations that couple the global two-time correlation, C, and linear response, R, averaged over the initial measure $\rho_{i.c.}$ and, also, the harmonic spring constants (quenched randomness), in the strict N → ∞ limit. This approach bears resemblance with dynamic mean theory [35]. The (replica) method used to impose the thermal initial conditions ensures symmetry breaking for $T'<J_0$. Four phases were identified in the $(J/J_0, \, T'/J_0)$ phase diagram (energy injection/initial condition characteristics) as deduced from $\chi_\infty = \lim_{t\to\infty} \int_0^t \textrm{d}t' \, R(t,t')$, which equals $1/J$ for $T'<J$ (II, III) and $1/T'$ for $T'>J$ (I, IV), and $q_0 = \lim_{t\to\infty} C(t,0)$, which takes a non-zero value for $T'<J_0$ and $T'<J$ (III), see fig. 1. The asymptotic value of the Lagrange multiplier is strictly larger than $\lambda_N$ for $T'>J$, whereas it locks to $\lambda_N=2J$ for $T'<J$ implying that the potential on the N-th mode flattens and the gap of the effective Hamiltonian closes for t → ∞ after N → ∞. Noteworthy, all these observables approach constant limits algebraically with superimposed oscillations [34].

Zoom In Zoom Out Reset image size

In this letter we work with a fixed (and typical) realization of the $\lambda_\mu$. On the one hand, we solve the coordinate dynamics for finite but large N and, ideally, long times with a modification of the semi-analytic phase-Ansatz method used in [36] to study the O(N) field theory, and adapted in [14] to the present case. With this method we compute the time averages $\overline{\langle s_\mu^2}\rangle_{i.c.}$ and $\overline{\langle p_\mu^2 }\rangle_{i.c.}$ (controlling the deviations from the ideal limit t → ∞ after N → ∞). The particle's behaviour in each Sector of the phase diagram of fig. 1 is sketched in fig. 3. In Sector I the particle is delocalised on the full sphere initially and it remains like this subsequently. In Sector II, the initial extended state with negligible projection on the N-th mode develops into configurations with a sub-extensive projection on this particular direction. In Sector III, the motion remains condensed, similarly to how it was initially. Finally, in Sector IV the N-th mode acquires extensive kinetic energy. More quantitative details are given later, when we compare dynamic and GGE averages.

On the other hand, we calculate the GGE partition sum

$$\begin{equation} Z_{\textrm{GGE}} = \int {\mathcal D}s \, {\mathcal D} p\, \textrm{d}z_{\textrm{GGE}} \, e^{-\sum_\mu \gamma_\mu I_\mu - \frac{z_{\textrm{GGE}}}{2} (\sum_\mu s_\mu^2 -N)}, \end{equation} \tag{ 13 }$$

with ${\mathcal D} s = \prod_\mu \textrm{d}s_\mu$, ${\mathcal D} p = \prod_\mu \textrm{d}p_\mu$ and $z_{\textrm{GGE}}$ the Lagrange multiplier that imposes the spherical constraint (which in this formulation could be reabsorbed in the definition of $\gamma_\mu$ thanks to $\sum_\mu I_\mu = C_1$). The standard Gibbs-Boltzmann equilibrium partition sum (relevant to describe the case $J=J_0$ and any $T'$) is recovered by setting $\gamma_\mu=-\beta'\lambda_\mu/2$ and $z_{\textrm{GGE}} = \beta' z_{\textrm{eq}}$. We evaluate the averages $\langle s^2_\mu\rangle_{\textrm{GGE}}$ and $\langle p^2_\mu\rangle_{\textrm{GGE}}$ that we compare to the dynamic ones $\overline{\langle s_\mu^2}\rangle_{i.c.}$ and $\overline{\langle p_\mu^2 }\rangle_{i.c.}$. We analyze the fluctuations of the constraints C_{1, 2} (dynamically and with the GGE) and from their scaling we determine in which cases the SNM is equivalent to the proper NM.

The partition sum $Z_{\textrm{GGE}}$ is a non-trivial object since the $I_\mu$ are quartic functions of the phase space variables, see eq. (5). Still, we managed to calculate it by adapting methods that are common in the treatment of disordered systems and random matrices [17]. Firstly, we used auxiliary variables to decouple the quartic terms. Secondly, for N → ∞, we transformed $\lambda_\mu$ into a continuous variable λ, all $N^{-1} \sum_\mu A_\mu$ into $\int \textrm{d}\lambda \, \rho(\lambda) A(\lambda)$ for any $A(\lambda)$, and $\sum_{\nu (\neq \mu)} \frac{A_\nu}{\lambda_\mu - \lambda_\nu} \mapsto \dashint \textrm{d}\lambda' \, {\rho(\lambda')} \, \frac{A(\lambda')}{\lambda-\lambda'}$ with $\dashint$ the Cauchy principal value. In both integrals $\rho(\lambda)$ is the density of spring constants taken to be of Wigner form. In some cases we separated the contribution of the N-th mode which may be macroscopic and scale differently from the ones in the bulk. Thirdly, we evaluated $Z_{\textrm{GGE}}$ by saddle-point. Then, we showed that the harmonic Ansatz

$$\begin{eqnarray} \begin{array}{rcl} \langle s^2(\lambda)\rangle_{\textrm{GGE}} &=& T(\lambda)/(\tilde z -\lambda),\\ \langle p^2(\lambda)\rangle_{\textrm{GGE}} &=& m T(\lambda) \end{array} \end{eqnarray} \tag{ 14 }$$

solves the saddle-point equations. Finally, we exploit the conditions $\langle I(\lambda)\rangle_{\textrm{GGE}} = \langle I(\lambda) \rangle_{i.c.}$, with $\langle I(\lambda)\rangle_{\textrm{GGE}} = - \partial \ln Z_{\textrm{GGE}} /\partial \gamma_\mu$ evaluated at the saddle point to close the system of equations. In the absence of initial condition condensation, $T'>J_0$, all Uhlenbeck constants are ${\mathcal O}(1)$ and their GGE averages should equal the ones over initial conditions given in eq. (11)

$$\begin{eqnarray} &&\langle I(\lambda)\rangle_{i.c.} =\frac{2T(\lambda)}{\tilde z -\lambda} \left[1-\dashint \textrm{d}\lambda' \frac{\rho(\lambda')\,T(\lambda')}{\lambda-\lambda'}\right]. \end{eqnarray} \tag{ 15 }$$

When the initial state is condensed, $T'<J_0$, eq. (15) applies to all λ with the exception of $\lambda_N$, for which

$$\begin{equation} \frac{ \langle I_N\rangle_{i.c.} } {2\langle s_N^2 \rangle_{\textrm{GGE}}} =\left[1-\dashint \textrm{d}\lambda' \frac{\rho(\lambda') T(\lambda')}{\lambda_N-\lambda'} - \frac{\langle s_N^2 \rangle_{\textrm{GGE}}}{2N}\right] \end{equation} \tag{ 16 }$$

plus o(1) corrections, and $\langle I_N\rangle_{i.c.}$ given in eq. (12). Together with the constraint $\langle C_1\rangle_{\textrm{GGE}} =N$, these are the central equations that allow us to solve the problem. Their numerical solution yields the spectrum of mode temperatures, $T(\lambda)$, $\tilde z$ and $\langle s^2_N\rangle_{\textrm{GGE}}$, and with them we deduce the expectation value of any observable. A selected number of results are shown in fig. 2 where we compare the GGE averages to the dynamic ones for parameters in Sectors I and IV of the phase diagram displayed in fig. 1. We collect dynamic data for N = 100, 1024 and GGE data for N = 100 (with a more naive saddle point evaluation) and N → ∞. The agreement is very good. The rather small extent of finite size effects in the bulk can also be appreciated in the figure (the double logarithmic scale enhances the appearance of the deviations, which are actually restricted to the neighborhood of the edge in (c) and (d)). In the insets in (a) and (c) the spectra of the Lagrange multipliers $\gamma_\mu$ for finite N are shown, which can be compared to the one of $T(\lambda)$. Results of similar quality are obtained in Sectors II and III (not shown).

Zoom In Zoom Out Reset image size

The dynamics in each Sector and the equivalence between the NM and the SNM can be rationalized according to the scaling properties of the last mode and the fluctuations of the constraints

$$\begin{eqnarray} &&\Delta C_i \equiv \langle (C_i - \langle C_i \rangle)^2 \rangle\qquad \textrm{for}\qquad i=1,2, \end{eqnarray} \tag{ 17 }$$

which can be studied both dynamically and with the GGE. When the scaling of these fluctuations is ${\mathcal O}(N^2)$ the SNM is not equivalent to the NM.

In Sector I, $\langle s^2_{\mu}\rangle_{\textrm{GGE}}$ and $\langle p^2_{\mu}\rangle_{\textrm{GGE}}$ are $\mathcal{O}(1)$ for all μ, including $\mu=N$. In a sense, this is the simplest possible generalization of the Boltzmann equilibrium extended phase. In Sector II, we have numerical evidence for $\langle s^2_{N} \rangle_{\textrm{GGE}}$ scaling as $N^{1/2}$, while $\langle s^2_{\mu\neq N}\rangle_{\textrm{GGE}}$ and $\langle p^2_{\mu}\rangle_{\textrm{GGE}}$ should be $o(N^{1/2})$. This is a quasi-condensed phase in which the weight of the last mode is large but not extensive. Since there is no condensation, the energy conserving dynamics in the extended and quasi-condensed phases explore the full sphere in the course of time as sketched in figs. 3(a), (b) with a red dot and the red sphere, respectively. Moreover, $\Delta C_i = o(N^2)$ and the NM and SNM models are equivalent in Sectors I and II.

Zoom In Zoom Out Reset image size

In Sectors III and IV, where the initial conditions are drawn from the Boltzmann measure of the SNM at $T'<J_0$. As explained above, they can be of two kinds: i) $s_N(t=0) \propto \pm \sqrt{N}$ with negligible fluctuations, or ii) $s_N(t=0)$ Gaussian distributed, centered at zero with $\sqrt{N}$ fluctuations [24,32]. In both cases $\langle s_N^2\rangle_{i.c.} \propto N$, but the dynamics are different and have to be discussed separately.

In case i), Sector III is a properly s-condensed phase with $\langle s^2_{N}\rangle_{\textrm{GGE}}$ scaling as N, while $\langle s^2_{\mu\neq N}\rangle_{\textrm{GGE}} = o(N)$ and $\langle p^2_{\mu}\rangle_{\textrm{GGE}}=\mathcal{O}(1)$. The system precesses around one of the two states with $|s_N| = {\mathcal O}(N^{1/2})$, the one selected by the symmetry-broken initial conditions, and comparably negligible projection on all other directions, see the symmetrically placed green dots and green trajectory in figs. 3(a), (b), respectively. The constraints C₁ and C₂ are strictly satisfied up to sub-extensive corrections and the NM and SNM are equivalent. Remarkably, in Sector IV both $\langle s^2_{N}\rangle_{\textrm{GGE}}$ and $\langle p^2_{N}\rangle_{\textrm{GGE}}$ scale as N, and the N-th mode captures ${\mathcal O}(N)$ kinetic energy. We call this Sector an s, p-condensed phase. The last mode is in a superposition of states associated to each initial condition. At any instant t, the configurations are distributed on an ellipse in the plane (s_N , p_N ) with axes ${\mathcal O}(N^{1/2})$, as in the closed motion of a harmonic oscillator, see the violet ellipse and cylinder in figs. 3(a), (b), respectively. The average over trajectories implies, in particular, that the limit correlation q₀ vanishes. The constraints C_{1, 2} are only verified on average over the initial conditions and the SNM and NM are not equivalent. We note that $\Delta C_{1,2}$ are averages of a quartic functions of the phase variables; had we evaluated only quadratic functions of s_N we would have not noticed the inequivalence between the two models. Quite surprisingly, the averaged dynamics cannot be boiled down to the ones of a typical trajectory with its own z(t).

In case ii), the initial conditions imply $\Delta C_1 = {\mathcal O}(N^2)$ at all times due to the large fluctuations of the last mode. One can show that, in Sector III, $\Delta C_2 = o(N^2)$ at all times. In this situation, due to the large fluctuations in C₁, zero-mean initial conditions are appropriate for the soft model but not for the strictly spherical one. In practice, in the SNM we average over spherical trajectories with different radius determined by the initial condition. In Sector IV, due to the condensation of p_N , the dynamics do not preserve the scaling properties of C₂ either. In other words, the fluctuations of the secondary constraint, which vanish in the initial condition, get macroscopically amplified by the dynamics. In conclusion, we average over trajectories that no longer move on the sphere. In this Sector, the fluctuations of all the quantities that are conserved on average, $H_{\textrm{quad}}$, C_{1, 2} and I_N , condense, which implies that the dynamics do not conserve the quadratic energy, are not restricted to a sphere and are not strictly integrable. The behaviours in Sectors III and IV are represented in figs. 3(c), (d), with the same colour code as the one we used before.

Contrary to the quantum mechanical subtleties [37,38], the notion of classical integrability is clear [8–10]. The dynamics should be ergodic on the portion of phase space compatible with the constants of motion [11]. Still, the fact that a canonical GGE could describe the time-averages of generic observables in a classical interacting integrable system is not obvious. We modified the celebrated Neumann model by imposing the spherical constraint on average over the initial conditions and we were then able to solve it in the thermodynamic limit. We thus provided an explicit example in which identities between temporal and statistical averages, for all kinds of thermal initial conditions (on average) and observables not correlated with the constants of motion and post-quench parameters, can be demonstrated.

Importantly enough, for condensed initial states, $\langle s^2_N\rangle_{i.c.}$ and $\langle I_N\rangle_{i.c.}$ are macroscopic and stay so after the quench. In these cases, we distinguished symmetry-broken initial conditions and symmetric ones with zero mean and condensed fluctuations. Quadratic observables are insensitive to the changes that the latter induce but quartic ones are not. For symmetry-broken initial conditions, the SNM behaves just as the NM in the phase in which only the coordinate is condensed but it loses its equivalence with the NM in the phase in which not only the coordinate but also the momentum condenses. For initial states with macroscopic fluctuations, although the dynamics conserve on average an extensive number of phase space functions, the individual trajectories do not. Energy conservation is violated in the condensed Sectors of the phase diagram and the SNM and NM are not equivalent. Interestingly enough, given the similarity between the phase transitions and condensation in this model and in BEC [25,33] we may expect similar phenomena in quenches of thermal initial states of the latter.

In conclusion, we provided a new class of classical integrable systems in which the stationary dynamics is captured by the GGE measure. Notice that the integrals of motion in the NM and SNM are non-local, in the sense that the $I_\mu$ are given by sums involving all coordinates and associated momenta. These models are, therefore, different from the one-dimensional cases, with (quasi)-local constants of motion, that are usually dealt with in equilibration studies of integrable systems.

We thank J.-B. Zuber for very helpful discussions.