Quenched dynamics of classical isolated systems: the spherical spin model with two-body random interactions or the Neumann integrable model

We study the Hamiltonian dynamics of the spherical spin model with fully-connected two-body interactions drawn from a Gaussian probability distribution. In the statistical physics framework, the potential energy is of the so-called $p=2$ spherical disordered kind. Most importantly for our setting, the energy conserving dynamics are equivalent to the ones of the Neumann integrable system. We take initial conditions in thermal equilibrium and we subsequently evolve the configurations with Newton dynamics dictated by a different Hamiltonian. We identify three dynamical phases depending on the parameters that characterise the initial state and the final Hamiltonian. We obtain the {\it global} dynamical observables with numerical and analytic methods and we show that, in most cases, they are out of thermal equilibrium. We note, however, that for shallow quenches from the condensed phase the dynamics are close to (though not at) thermal equilibrium. Surprisingly enough, for a particular relation between parameters the global observables comply Gibbs-Boltzmann equilibrium. We next set the analysis of the system with finite number of degrees of freedom in terms of $N$ non-linearly coupled modes. We evaluate the mode temperatures and we relate them to the frequency-dependent effective temperature measured with the fluctuation-dissipation relation in the frequency domain, similarly to what was recently proposed for quantum integrable cases. Finally, we analyse the $N-1$ integrals of motion and we use them to show that the system is out of equilibrium in all phases, even for parameters that show an apparent Gibbs-Boltzmann behaviour of global observables. We elaborate on the role played by these constants of motion in the post-quench dynamics and we briefly discuss the possible description of the asymptotic dynamics in terms of a Generalised Gibbs Ensemble.


Introduction
In the past decade, atomic physics experiments have been able to test the global coherent quantum dynamics of interacting systems. This achievement has boosted research on the dynamics and possible equilibration of many-body isolated systems [1]. Some of the quantum instances realised in the laboratory are low dimensional and considered to be integrable. Therefore they are not able to act as a bath on themselves and questions on how to describe their asymptotic dynamics pose naturally. With the aim of describing their asymptotic states, the Generalised Gibbs Ensemble (GGE), an extension of the canonical Gibbs-Boltzmann density operator that aims to include the effect of all relevant conserved charges, was proposed [2,3] (see, e.g., the review articles [4,5,6]).
Similar equilibration problems can arise in classical isolated systems. A first study of the dynamics of isolated interacting mean field disordered models appeared in [7]. We continue developing this project and we analyse, in this paper, the quench dynamics of a classical integrable system with (weak) frozen randomness. Both models belong to the class of p spin spherical disordered models with, however, properties that render their constant energy dynamics very different, as we will show here.
The spherical fully-connected p-spin disordered models are paradigms in the mean-field description of glassy physics. They are solvable models (in the thermodynamic limit) that successfully mimic the physics of (fragile) glasses for p ≥ 3 and domain growth for p = 2. The connection between the p = 2 model, in its classical and quantum formulations, with coarsening phenomena is made via its relation to the celebrated O(N ) λφ 4 model in the infinite N limit. Furthermore, the model has recently appeared in a semiclassical study of the Sachdev-Ye-Kitaev model [8]. The literature on the static, metastable and stochastic dynamics of the p spin spherical systems is vast. Numerous aspects of their behaviour are very well characterised, even analytically (see, e.g., the review articles [9,10,11]).
In Ref. [7] we studied the Hamiltonian dynamics of the p = 3 spherical disordered spin model. By adding a kinetic term to the standard potential energy we induced energy conserving dynamics to the real valued spins. In this setting, the dynamics correspond to the motion of a particle on an N − 1 dimensional sphere under the effect of a complex quenched random potential [12,13,14]. Here we will focus on the Newtonian dynamics of the particle under conservative forces arising from a quadratic potential, the p = 2 case.
The Hamiltonian p = 2 disordered model turns out to be equivalent to the Neumann integrable system of classical mechanics [15], the constrained motion of a classical particle on SN−1 under a harmonic potential, for a special choice of the spring constants. The only difference is that in the p = 2 model one imposes the spherical constraint on average while in Neumann's model one does strictly, on each trajectory. This difference, however, should not be important in the N → ∞ limit. We will exploit results from the Integrable Systems literature, notably the exact expressions of the N − 1 conserved charges in involution [16,17,18]. With these at hand, we will be able to study the statistical properties in depth and construct a candidate GGE to describe the long-time dynamics.
We perform instantaneous quenches towards a post-quench disordered potential that keeps memory of the pre-quench one, mimicking in this way the "quantum quench" procedure in a classical setting. The change in the potential energy landscape induces finite injection or extraction of energy density in the sample. The subsequent dynamics conserve the total energy. We sample the initial conditions from canonical equilibrium at a tuned temperature, choosing in this way initial configurations typical of a paramagnetic equilibrium state at high temperature or a condensed, ferromagnetic-like, state at low temperature. The control parameters in the dynamic phase diagram that we will establish are the amount of energy injected or extracted and the initial temperature of the system, both measured with respect to the same energy scale.
The dynamic evolution in the different phases of the phase diagram will be pretty different, with cases in which the infinite size system remains confined (condensed, in the statistical physics language) and cases in which it does not. In none of them a Gibbs-Boltzmann equilibrium measure is reached, contrary to what happens in the strongly interacting p = 3 case. The role played by the N − 1 integrals of motion on the lack of equilibration of the infinite size system will be discussed.
The reader just interested in a summary of our results and not so much in the way in which we obtained them can go directly to Sec. 8 where we sum up our findings and we present a thorough comparison between the dynamics of the isolated p = 2 and p = 3 cases.
The paper is organised as follows. In Sec. 2 we recall the main features of the p = 2 spherical disordered model (static, metastable and relaxational dynamic properties) studied from the statistical mechanics point of view. In the following Section we explain the relation between the disordered spin system and the integrable Neumann model of classical mechanics. We also explain in this Section the statistical description of the longterm dynamics of integrable systems provided by the Generalised Gibbs Ensemble proposal. In Sec. 4 we present our analytic results for the dynamics of the model in the N → ∞ limit and in Sec. 5 we go further and we set the analysis of the evolution of the finite N case. Section 6 is devoted to the numerical study of the N → ∞ and finite N dynamic equations. In Sec. 7 we investigate the behaviour of the N − 1 integrals of motion in the various sectors of the phase diagram and we discuss their influence against the equilibration of the system. Finally, Sec. 8 presents our conclusions. Several Appendices complement the presentation in the main part of the paper.

Background
This Section presents a short account of the equilibrium properties and relaxation dynamics of the spherical p = 2 disordered model, first introduced and studied by Kosterlitz, Thouless and Jones [19] as the simplest possible magnetic model with quenched random interactions. This model, as we will explain below, shares many points in common with the O(N ) model of ferromagnetism when treated in the infinite N limit. Its static properties have been derived with a direct calculation and using the replica trick. Its relaxational dynamics are also analytically solvable. The reader familiar with this model can jump over this Section and go directly to the next one where the relation with Neumann's model and integrability are discussed.

The Hamiltonian spherical p = 2 spin model
The p = 2 spin model is a system with two-spin interactions mediated by quenched random couplings Jij. The potential energy is The coupling exchanges are independent identically distributed random variables taken from a Gaussian distribution with average and variance The parameter J characterises the width of the Gaussian distribution. In its standard spin-glass setting the spins are Ising variables and the model is the Sherrington-Kirkpatrick spin-glass. We will, instead, use continuous variables, − √ N ≤ si ≤ √ N with i = 1, . . . , N , globally forced to satisfy (on average) a spherical constraint, N i=1 s 2 i = N , with N the total number of spins [19]. Such spherical constraint is imposed by adding a term to the Hamiltonian, with z a Lagrange multiplier. The spins thus defined do not have intrinsic dynamics. In statistical physics applications their temporal evolution is given by the coupling to a thermal bath, via a Monte Carlo rule or a Langevin equation [20]. The quadratic model is a particular case of the family of p-spin models, the celebrated mean-field model for glasses, with potential energy Hpot = − i 1 ...ip Ji 1 ...ip si 1 . . . si p and p integer. Even more generally, the form (1) is one instance of a generic random potential V ({si}) with zero mean and correlations [12,13,14] [ with V(| s − s |/N ) = − J 2 2 ( s · s /N ) 2 . Similarly to what was done in [7] in the study of the p ≥ 3 Hamiltonian dynamics, the model can be endowed with conservative dynamics by changing the "spin" interpretation into a "particle" one. In this way, a kinetic energy [21,22] can be added to the potential energy. The total energy of the Hamiltonian spherical p-spin model is then Hsyst = H kin + Hpot + Hconstr .
This model represents a particle constrained to move on an N -dimensional hyper-sphere with radius √ N . The position of the particle is given by an N-dimensional vector s = (s1, . . . , sN ) and its velocity by another Ndimensional vector˙ s = (ṡ1, . . . ,ṡN ). The N coordinates si are globally constrained to lie, as a vector, on the hypersphere with radius √ N . The velocity vector˙ s is, on average, perpendicular to s, due to the spherical constraint. The parameter m is the mass of the particle.
The generic set of N equations of motion for the isolated system is i = 1, . . . , n, where the Lagrange multiplier needs to be time-dependent to enforce the spherical constraint in the course of time.
The initial condition will be taken to be {s 0 i ,ṡ 0 i } ≡ {si(0),ṡi(0)} and chosen in ways that we specify below. We will be interested in using equilibrium initial states drawn from a Gibbs-Boltzmann measure at different temperatures T .
From Eq. (7) one derives an identity between the energy density and the Lagrange multiplier. By multiplying the equation by si(t2) and taking t2 → t1 The first term can be rewritten as m lim t 2 →t − The Lagrange multiplier takes the form of an action density, as a difference between kinetic and potential energy densities. Using now the conservation of the total energy, e f = epot(t1) + e kin (t1), The p = 2 model belongs to a different universality class from the one of the p ≥ 3 cases, in the sense that its free-energy landscape and relaxation dynamics are much simpler. It is, indeed, a model that resembles strongly the large N , O(N ) model for ferromagnetism. A hint on the simpler properties of its potential energy landscape is given by the fact that the equations derived for general p simplify considerably for p = 2. For example, the static and dynamic transitions occur at the same temperature Tc = T d , and the number of metastable states is drastically reduced. We recall these properties in the rest of this Section.

The statics
The static properties of the p = 2 spherical model were elucidated in [19]. The trick is to project the spin vector s onto the basis of eigenvectors of the interaction matrix. One calls λµ and vµ the µ-th eigenvalue and eigenvector of the matrix Jij, and sµ = s · vµ the projection of s on the eigenvector vµ. In terms of the latter the Hamiltonian is not only quadratic but also diagonal. The extrema of the potential energy landscape and the partition function can then be easily computed. In the thermodynamic limit, N → ∞, the eigenvalues are distributed according to the Wigner-Dyson semi-circle form [23] ρ(λµ) = 1 2πJ 2 (2J) 2 − λ 2 µ θ(2J − |λµ|) . (11) For finite N the distance between the largest and next to largest eigenvalues is order N −1/6 /N 1/2 = N −2/3 .

The potential energy landscape
Let us label the eigenvalues of Jij in such a way that they are ordered: λ1 ≤ λ2 ≤ · · · ≤ λN . We call their associated eigenvectors vµ with µ = 1, . . . , N and we take them to be orthonormal, such that v 2 µ = 1. We consider the potential energy landscape HJ ({si}, z) with z taken as a variable.
These stationary points are metastable states at zero temperature, apart from two of them that are the marginally stable ground states, and their number is linear in N , the number of spins. (The role of marginal stability in the physical behaviour of condensed matter systems was recently summarised in [24].) These statements are shown in the way described in the next paragraph. The Hessian of the potential energy surface on each stationary point is ∂HJ ∂si∂sj s * ,z * = −Jij + zδij| s * ,z * = −Jij + λµδij .
This matrix can be easily diagonalised and one finds Dνη = (−λν + λη)δνη. Thus, on the stationary point, s * = ± √ N vµ, the Hessian has one vanishing eigenvalue (for ν = µ), µ − 1 positive eigenvalues (for ν < µ), and N − µ negative eigenvalues (for ν > µ). Positive (negative) eigenvalues of the Hessian indicate stable (unstable) directions. This implies that each saddle point labeled by µ has one marginally stable direction, µ − 1 stable directions and N − µ unstable directions. (In other words, the number of stable directions plus the marginally stable one is given by the index µ labelling the eigenvalue associated to the stationary state.) In conclusion, there are two maxima, s * = ± √ N v1, in general two saddles s * = ± √ N vI with I = µ − 1 stable directions and N − I unstable ones, with I running with µ as I = µ − 1 and µ = 3, . . . , N , and finally two (marginally stable) minima, s * = ± √ N vN . In the large N limit the density of eigenvalues of the Hessian at each metastable state µ is a translated semi circle law [25].
The zero temperature energy of a generic configuration under no applied field is At each stationary point s * = ± √ N vµ and z * = λµ this energy is Here we used the notation v µ i to indicate the ith component of the µth eigenvector vµ. The energy difference between the minima and the lowest saddles depends on the distribution of eigenvalues, a semi-circle law for the Gaussian distributed interaction matrices that we consider here.
A magnetic field reduces the number of stationary points from a macroscopic number to just two. Indeed, the stationary state equation now reads Jijsj + zsi − hi| s * ,z * = 0 , ∀i , ⇒ s * i = (z * − J) −1 ij hj and z * is fixed by imposing the spherical constraint on s * . One then finds two solutions for the Lagrange multiplier that lie outside the interval of variation of the eigenvalues of the Jij matrix: |z * | > λN . The stability analysis shows that the stationary points are just one fully stable minimum and a fully unstable maximum. The elimination of saddle-points by an external field has important consequences on the dynamics of the system [26].
In this paper we do not apply any external field. The analysis of large dimensional random potential energy landscapes [27,28,29,30] is a research topic in itself with implications in condensed matter physics, notably in glass theory [31,24], but also claimed to play a role in string theory [32,33], evolution [34] or other fields. The p = 2 spherical model provides a particularly simple case in which the potential energy landscape can be completely elucidated.

The free-energy density
This special (almost) quadratic model allows for the complete evaluation of its free-energy density for a typical realisation of the disordered exchanges. The traditional derivation of the disorder averaged free-energy density can also be done using the replica method and a simple replica symmetric Ansatz solves this problem completely. We recall how the two methods [19,35] work in this Section.
The partition function reads where c is a real constant to be fixed below.
It is convenient to diagonalise the matrix Jij with an orthogonal transformation and write the exponent in terms of the projection of the spin vector s on the eigenvectors of Jij, sµ ≡ s · vµ. This operation can be done for any particular realisation of the interaction matrix. The new variables sµ are also continuous and unbounded and the partition function can be recast as Assuming that one can exchange the quadratic integration over sµ with the one over the Lagrange multiplier, and that c is such that the influence of eigenvalues λµ > c is negligible, one obtains In the saddle-point approximation valid for N → ∞ the Lagrange multiplier is given by (17) and this equation determines the different phases in the model. We indicate with double brackets the sum over the eigenvalues of the matrix Jij that in the limit N → ∞ can be traded for an integration over its density: Let us discuss the problem in the absence of a magnetic field. The high temperature solution to Eq. (17) can be smoothly continued to lower temperatures until the critical point is reached where zsp touches the maximum eigenvalue of the Jij matrix, and it sticks to it for all T < Tc = J: Tc is the static critical temperature. (A magnetic field with a component on the largest eigenvalue, h · vmax = 0, acts as an ordering field and erases the phase transition.) If one now checks whether the spherical constraint is satisfied by these saddle-point Lagrange multiplier values, one verifies that it is in the high temperature phase, but it is not in the low temperature phase, where N µ=1 The way out is to give a macroscopic weight to the projection of the spin in the direction of the eigenvector that corresponds to the largest eigenvalue: with δsN = 0 so that The thermal average of the projection of the spin vector on each eigenvalue vanishes in the high temperature phase and reads below the phase transition (once we have chosen one of the ergodic components with the spontaneous symmetry breaking of the s → − s invariance). The configuration condenses onto the eigenvector associated to the largest eigenvalue of the exchange matrix that carries a weight proportional to √ N . Going back to the original spin basis, the mean magnetisation per site is zero at all temperatures but the thermal average of the square of the local magnetisation, that defines the Edwards-Anderson parameter, is not when T < Tc: The order parameter qEA vanishes at Tc and the static transition is of second order. The condensation phenomenon occurs for any distribution of exchanges with a finite support. If the distribution has long tails, as when the model is defined on a sparse random graph [36,37,38], the energy density diverges and the behaviour is more subtle [39,40].
The disorder averaged free-energy density can also be computed using the replica trick [41] and a replica symmetric Ansatz. This Ansatz corresponds to an overlap matrix between replicas Q ab = δ ab + qEA ab with ab = 1 for a = b and ab = 0 for a = b. When N → ∞ the saddle point equations fixing the parameter qEA yield 0 above Tc and a marginally stable solution with qEA = 1 − T /Tc and identical physical properties to the ones discussed above below Tc.
The equilibrium energy is given by We added a superscript pot since in the modified model that we will study in this paper the total energy will also have a kinetic energy contribution. The entropy diverges at low temperatures as ln T , just as for the classical ideal gas, as usual in classical continuous spin models.

Relaxation dynamics
The over-damped relaxation dynamics of the spherical p = 2 spin model (coupled to a Markovian bath) were studied in [42,43,44,26,45]. One of the settings considered in these papers evolve a completely random initial condition, {s 0 i }, that corresponds (formally) to an infinite temperature initial state. The system is then subject to an instantaneous temperature quench by changing the temperature of the bath to a final value T < +∞. Initial conditions drawn from equilibrium at temperature T < Tc, and evolving in contact with a bath at the same temperature, T = T , were considered in [44] and it was shown in this paper that equilibrium at the same temperature is maintained ever after. A quench of the dissipative system from equilibrium at T < Tc to another subcritical temperature T < Tc was also studied in [44] and it was there shown that equilibrium at the target temperature is achieved.
The coupling to the bath is modeled with a stochastic equation of Langevin kind. This equation can be exactly solved in the basis of eigenvectors of the interaction matrix and the Lagrange multiplier z(t) can be fixed by imposing the spherical constraint at all times. This yields a self-consistent equation for z(t). Notably, the asymptotic solution depends on the choice of initial state, as we expose below. The applied field hµ is used to compute the linear response function For quenches of initial conditions drawn from equilibrium at T → ∞, and evolution in contact with a bath at temperature T > Tc, the dynamics quickly approach equilibrium at the new temperature. The Lagrange multiplier quickly converges to zeq = T + J 2 /T . The correlation and linear response are invariant under translations of time and they are related by the fluctuation dissipation theorem [44,46], see Eq. (33) below.
For quenches of initial conditions drawn from equilibrium at T → ∞, and evolution in contact with a bath at temperature T < Tc, the correlation and linear response functions behave as in coarsening systems [47,48,49,50], decaying in two time regimes, one stationary for short time differences, (t1 − t2)/t2 1, and one non-stationary for long time differences(t1 − t2)/t2 1. The detailed time-dependence in the two regimes can be extracted using the procedure sketched above. It yields the behaviour of the self correlation and linear response that scale in the same way as these do in the O(N ) model of ferromagnetism studied in the large N limit, see Sec. 2.4. The progressive condensation of the spin "vector" in the direction of the eigenvector corresponding to the largest eigenvalue of the interaction matrix is the equivalent of the ordering process in the O(N ) model, that is to say, the condensation on the zero wave-vector mode. Complete alignment with an overlap of order √ N is not reached in finite times with respect to N .
For low temperature quenches from random initial conditions, the Lagrange multiplier approaches z f = 2J as a power law, z(t) − z f −3/(4t). The slow approach to the asymptotic value is determinant to allow for the non-stationary slow relaxation. The global correlation and linear response are computed from the spin solution sµ(t) and they can be cast as [44] with the stationary and a non-stationary terms linked by the FDT at the temperature of the bath and a modified FDT at an effective temperature T eff [51,52] selected by the dynamics, always with t1 ≥ t2. In the asymptotic limit, the two terms added to form C and R evolve in different regimes in the sense that when one changes the other one is constant and vice versa. The limiting values of the two contributions to the correlation function are with the parameter q being equal to This is the correct expression of the Edwards-Anderson parameter for the equilibrium low temperature solution, see Eq. (26), and once again one finds Tc = J from q = 0. The complete solution of the Langevin equations allows one to deduce the exact scaling forms of the stationary and ageing contributions to the correlation and linear response. These are with fC (x) and fR(x) known analytically. The behaviour of the effective temperature is special in the p = 2 model in the sense that contrary to what happens in the p ≥ 3 cases [20] it is not constant but grows with time.
More precisely, it scales as T eff (t1, t2) t 1/2 fT (t1/t2) and it diverges asymptotically as t 1/2 1 . This implies, in particular, that the ageing regime does not contribute to the asymptotic potential energy that, after a quench to T < Tc, reads (39) and is identical to the equilibrium one, see Eq. (27), once q = 1 − T /J is used. For quenches within the ordered phase, T < Tc → T < Tc, for example taking initial conditions in equilibrium at zero temperature, sµ(0) = √ N δµ,N and evolving them at T < Tc, the Lagrange multiplier approaches z f = 2J faster than any power law and the system rapidly equilibrates to the after quench conditions [44].
The same technique, based on the projection of the spin vector on the eigenvectors of Jij, can also be implemented in the case in which there is inertia and the differential equation has a second order time derivative. The dynamics are recast into the ones of harmonic oscillators coupled by a self-consistent time-dependent Lagrange multiplier. We will use this formulation in Sec. 5. Although a full analytical solution is not possible with the second time derivative, a performant numerical algorithm will allow us to monitor the evolution of the different modes.

Relation with the O(N ) λφ 4 model in the large N limit
The λφ 4 scalar field theory in d dimensions is defined by the Hamiltonian where r and λ are two parameters. This model is the Ginzburg-Landau free energy for the local order parameter of the paramagnetic-ferromagnetic transition controlled by the parameter r going from r > 0 to r < 0. When the field is upgraded to a vector with N components and the limit N → ∞ is taken the quartic term, first conveniently normalised by N , can be approximated by λ The quantity z(t) = φ 2 is not expected to fluctuate and plays the role of the Lagrange multiplier in the spherical disordered model. Once this approximation made, the model becomes quadratic in the field and its statics and relaxation dynamics can be easily studied. The only difficulty lies in imposing the self-consistent constraint that determines z(t). The condensation phenomenon that we discussed in the disordered model is also present in the field theory and it corresponds to a condensation on the zero wave vector mode. In the dynamic problem this corresponds to the progressive approach to the homogeneous field configuration [9].
The conserved energy dynamics of the λφ 4 model, especially after sudden quenches, has been studied by a number of groups. Details on the behaviour of the scalar problem, as well as a review of general equilibration and pre-equilibration issues can be found in J. Berges' Les Houches Lecture Notes [53], see also [54]. The dynamics of the large N limit of the O(N ) model was analysed in [55,56,57,58,59,60]. More recent works use renormalisation group techniques to study the short time dynamics [61,62] at the dynamic phase transition.

Neumann's model, integrability and equilibration
In this Section we explain the relation between the Hamiltonian p = 2 disordered model and the integrable model of Neumann [15]. We start by recalling some basic properties of classical integrable systems in the sense of Liouville [63,64]. We then recall the definition of Neumann's model and we compare it to the p = 2 one. Finally, we explain the ideas behind the Generalised Gibbs Ensemble. Later, in Sec. 7, we will use this formalism to analyse certain aspects of the long time dynamics of the system, and we set the stage for a future study of the eventual approach to a GGE ensemble.

Integrable systems
In classical mechanics, systems are said to be Liouville integrable if there exist sufficiently many well-behaved first integrals or constant of motions in involution such that the problem can then be solved by quadratures [63,64], in other words, the solution can be reduced to a finite number of algebraic operations and integrations. In more concrete terms, an integrable dynamical system consists of a 2N -dimensional phase space Γ together with N independent functions 1 O1, . . . , ON : Γ → R, such that the mutual Poisson brackets vanish, We will assume henceforth that the Oi do not depend explicitly on time and that dOi/dt = 0 is equivalent to {H, Oi} = 0. Conventionally, the first function O1 is the Hamiltonian itself and the first constant of motion is the energy. All other Oi with i = 1 are also constants of motion since their Poisson bracket with H vanishes. The dynamics of the system can then be seen as the motion in a manifold of dimension 2N − N = N in which all configurations share the initial values of all the conserved quantities Oi(t) = Oi(0). Under these conditions Hamilton's equations of motion are solvable by performing a canonical transformation into action-angle variables (Ii, φi) with i = 1, . . . , N such that the Hamiltonian is rewritten asH(I) and The action functions Ii are conserved quantities and we collected them in I in the dependence of the frequencies ωi and the HamiltonianH. The remaining evolution is given by N circular motions with constant angular velocities. Both deciding whether a system is integrable and finding the canonical transformation that leads to the pairs (Ii, φi) are in practice very difficult questions. Whenever the system is integrable, and one knows the action-angle pairs, the statement in Eq. (42) is part of the Liouville-Arnold theorem [65].

Neumann's model and its integrals of motion
The model proposed by Neumann in 1850 describes the dynamics of a particle constrained to move on the N − 1 dimensional sphere under the effect of harmonic forces [15]. The Hamiltonian is where the L kl are the elements of an angular momentum anti-symmetric matrix and p k and x k are phase space variables with canonical Poisson brackets {x k , p l } = δ kl . The global spherical constraint N k=1 ensures that the motion takes place on SN−1. Using the fact that L kk = 0 to rewrite the double sum in the first term in H as an unconstrained sum, and replacing L kl by its explicit form in terms of x k and p k , one derives m k =l L 2 kl = m k,l L 2 kl = 2 k x 2 k l p 2 l − 2 k x k p k l x l p l . Imposing next the spherical constraint, that also implies k x k p k = 0, the sum simply becomes k =l We note that we added a factor 1/N in ifront of the kinetic energy in Eq. (43) in order to ensure that the two terms in N be extensive and the thermodynamic limit non-trivial. It is quite clear that Neumann's model is therefore identical to the Hamiltonian p = 2 model once the latter is written in the basis of eigenvectors of the interaction matrix Jij.
The N − 1 integrals of motion of this problem were constructed by K. Uhlenbeck [16] and more recently rederived by Babelon & Talon [18] with a separation of variables method. In a notation that is convenient for our application they read and satisfy k I k = N and 1 2 k a k I k = H. In the definition of our Hamiltonian and equations of motion we used a convention such that a k → −λµ (note the minus sign). After a trivial translation to the variables of the p = 2 spherical model we then have µ Iµ = 1 and µ λµIµ = −2H.

Statistical measures for integrable systems
Let X = (x1, p1, . . . , xN , pN ) be a generic point in phase space. The fact that the microcanonical measure with c = d X N j=1 δ(Ij( X) − ij) be sampled asymptotically is ensured by the Liouville-Arnold theorem [65], if the frequencies of the periodic motion on the torus are independent, that is, k · ω = 0 for k = (k1, . . . , kN ) with integer kj has the unique solution k = 0. One can call this ensemble the Generalized Microcanonical Ensemble.
In principle, the Generalized Canonical Ensemble, commonly called Generalized Gibbs Ensemble (GGE), can now be constructed from the Generalized Microcanonical Ensemble following the usual steps. The idea is to look for the joint probability distribution of N extensive (as for the Hamiltonian in the usual case) constants of motion of a subsystem P (i1, . . . , iN )di1 . . . diN . As in cases with just one conserved quantity, it is convenient to interpret P as a probability over position and momenta variables, and write This form can be derived under the same kind of assumptions used in the derivation of the canonical measure from the microcanonical one, that is (i) independence of the chosen subsystem with respect to the rest, in other words, the factorisation of the density of states g({i1, . . . , iN }) = g({i 1 , . . . , i N }), (ii) additivity of the conserved quantities i j ), (iv) constant inverse 'temperatures' ζj ≡ ∂i j S2({i1, . . . , iN }) = kB∂i j ln g2({i1, . . . , iN }). An inspiring discussion along these lines appeared in [66]. The conditions just listed imply a locality requirement on the ijs, otherwise (ii) and (iii) would be violated. This is similar to the requirement of having short-range interactions to derive the equivalence between the canonical and microcanonical ensembles in standard statistical mechanics.
In quenching procedures, the parameters ζj should be determined by requiring that the expectation value of each conserved quantity Ij calculated on ρGGE matches the (conserved) initial value Ij(0 + ) (right after the quench): The ζj are then the Lagrange multipliers that enforce this set of N constraints. In the p = 2 or Neumann model a set of conserved quantities in involution are the Iµ defined in Eq. (48). We will study them in Sec. 7. We note that if the Lagrange multipliers became, under some special conditions −λµβ f /2, with λµ the eigenvalues of the random interaction matrix, the GGE measure would be the Gibbs-Boltzmann one.

Averages in the long time limit
Take now a generic function of the phase space variables A( X) that does not depend explicitly on time and is not conserved. Birkhoff's theorem [67] states that its long-time average exists and reaches a constant, for τ sufficiently long and tst a reference transient time. We will use this fact at various points in our study.
The claim of equilibration to a Generalised Gibbs Ensemble is that the long time averages should also be given by the averages over the statistical measure PGGE: Which are the observables for which this result should hold is an interesting question that needs to be answered with care. The GGE proposal [2,3] and most, if not all, of its discussion appeared in the treatment and study of quantum isolated systems and, especially, the dynamics following an instantaneous quench performed as a sudden change in a parameter of the system's Hamiltonian. A series of review articles are [4,5,6]. The main motivation for our research project is to ask similar questions in the classical context, with the aim of disentangling the quantum aspects from the bare consequences of isolation and integrability.

The GGE temperatures and the fluctuation dissipation theorem
The fluctuation-dissipation theorem (FDT) [68] is a model independent equilibrium relation between the time-delayed linear response of a chosen observable and its companion correlation function. In Gibbs-Boltzmann equilibrium this relation is independent of the specific system and observable and it only involves the inverse temperature of the system. For classical systems it admits simple expressions in the time and frequency domains 2 Out of canonical equilibrium, the fluctuation-dissipation relations (FDR) between the linear response and the correlation function have been used to quantify the departure from equilibrium [20]. Indeed, the possibly time and observable dependent parameter that replaces β in far form equilibrium systems yields an effective temperature that in certain cases with slow dynamics admits the interpretation of a proper temperature [51,52]. Specially useful for our purposes is the fluctuation dissipation relation (FDR) in the frequency domain that concretely defines the frequency dependent, and also possibly observable dependent, inverse effective temperature β AB eff . It was shown in [69,70] that the Lagrange multipliers ζj of the GGE, seen as inverse temperatures βj, of a number of isolated integrable quantum systems which reach a stationary state can be read from the FDR's of properly chosen observables βj = β eff (ωj) for all j .
In this paper we will show that this statement also applies to the classical integrable system that we analyse.

Analytic results for the dynamics of the infinite size system
We now enter the heart of our study and we consider the dynamics of the isolated system after different kinds of quenches. In this Section we use an analytic treatment of the global dynamics in the thermodynamic limit. Long time regimes will be considered only after the diverging number of degrees of freedom:

Dynamical equations
We start by giving a short description of steps that lead to the dynamic equations that couple linear response and correlation function and fully characterise the evolution of the model in the N → ∞ limit. 2 When dealing with the numerical data we used a Fourier transform convention such that R(ω, t 2 ) = k R(t k + t 2 , t 2 )e iωt k and C(ω, t 2 ) = k (C(t k + t 2 , t 2 ) − q) cos(ωt k ) with t k the discrete times on which we collect the data points. The Fourier transform of the correlation is computed for t k > 0 only, taking advantage of the long-time stationarity property Cst(−t) = Cst(t). For this reason, there is no factor 2 in the left-hand-side of the FDT in the frequency domain.

The Schwinger-Dyson equations
In the N → ∞ limit, the only relevant correlation and linear response functions are for t1, t2 > 0, where the infinitesimal perturbation h is linearly coupled to the spin H → H − h N i=1 si at time t2 and the upperscript (h) indicates that the configuration is measured after having applied the field h. Since causality is respected, the linear response is non-zero only for t1 > t2. The square brackets denote here and everywhere in the paper the average over quenched disorder. The angular brackets indicate the average over thermal noise if the system is coupled to an environment, and over the initial conditions of the dynamics sampled, say, with a probability distribution. When the coupling to the bath is set to zero, as we do in this paper, the last average is the only one remaining in the angular brackets operation. The meaning of the indices a, b is given in the next paragraph.
The dynamical equations starting from a random state are well-known and can be found in Refs. [71,20,72,73]. They are usually derived from the dynamical Martin-Siggia-Rose-Janssen-deDominicis generating function. The method has been modified to take into account the effect of equilibrium initial conditions in [74] and it was applied to the relaxational p spin model in [75,76]. The average over disorder now becomes non-trivial and needs the use of the replica trick. The scripts a, b indicate then the replica indices a, b = 1, . . . , N . For initial conditions in equilibrium the replica structure is replica symmetric (see Sec.

and [19]), with
and Q a =b = qin (62) and qin in the paramagnetic state while qin = 0 in the condensed phase. This structure has an effect on the equation for the time-dependent correlation function that will keep the initial replica structure. There will be two kinds of correlations with the initial condition where we singled out the replica a = 1. Since there is no reason to think that the replicas that are not a = 1 behave differently, we follow the dynamics of C2(t1, t2) as a representative of this group. The interpretation of the correlations C1(t1, t2) and C2(t1, 0) can be given in terms of real replicas. The former is the self-correlation between the configuration of one replica of the system {si}(t2) and the same replica evolved until a later time t1, {si}(t1). For this reason, we will eliminate the subscript 1 and call C1(t1, t2) → C(t1, t2) in most places hereafter. The latter is the correlation between one replica of the system {σi}(0) at the initial time 0 and another replica of the system evolved until time t1 and represented by {si}(t1). Although we could also write an evolution equation for the two-time C2(t1, t2) we do not need it here since only C2(t1, 0) intervenes in the other equations.
In the N → ∞ limit one derives the dynamical equations that read The border conditions are Note that the initial condition is not the same for all C b (t1, 0). It is equal to 1 for b = 1 and equal to qin for b = 1. One can check that these equations coincide with the ones in [75,76] when inertia is neglected, p = 2 and J = J0, and a coupling to a bath is introduced. With respect to the equations studied in [7], they correspond to p = 2 and they have the extra ingredient of the influence of equilibrium initial conditions with a non-trivial replica structure, allowing for condensed initial states in proper thermal equilibrium. The sums over the replica indices appearing in Eqs. (66), (67) and (68) can be readily computed in the n → 0 limit; they read Consequently, the terms induced in the equation for C(t1, t2) and C2(t1, 0) are different. With inertia and no coupled bath, the equal-time conditions are for all times t1, t2 larger than or equal to 0 + , when the dynamics start. We found convenient to numerically integrate the integro-differential equations using an expression of the Lagrange multiplier that trades the second-time derivative of the correlation function into the total conserved energy after the quench. Following the same steps explained in [7] we deduce where we used Eq. (70) evaluated at t1 = t2. It seems that we have simply traded z(t1) by e(t1). Indeed, taking advantage of the fact that for an isolated system e(t1) = e f , a constant, the numerical solution of the evolution equations becomes now easier since it does not involve the second time derivative of the correlation function. In practice, in the numerical algorithm we fix the total energy e f to its post-quench value derived in Sec. 4.3, and we then integrate the set of coupled integro-differential equations with a standard Runge-Kutta method.
In short, the set of equations that fully determine the evolution of the system from an initial condition in canonical Boltzmann equilibrium at any temperature T are High and low temperature initial states are distinguished by qin = 0 for T > J0, and qin = 1 − T /J0 for T < J0, respectively. The equation for C(t1, 0) is just the one for C(t1, t2) evaluated at t2 = 0 so we do not write it explicitly.

Constant energy dynamics
In order to ensure constant energy dynamics we set J = J0 and m = m0 in this Subsection. We verify that the equations consistently conserve the equilibrium conditions. Moreover, we derive a number of properties of the linear response function that will be useful in the analysis of the instantaneous quenches as well.

Consistency with equilibrium parameters
The equation for C2(t, 0) admits the solution C2(t1, 0) = qin = cst in the case in which no quench is performed. Indeed, setting C2(t1, 0) = qin in Eq. (77) one has z(t1)qin = J 2 qin This equation has the solution qin = 0, the one of the paramagnetic phase, and a non-vanishing qin = 0 solution relevant in the ordered phase. Let us now focus on the case qin = 0. Using FDT, a property of equilibrium, the integral can be performed, the contribution from t = 0 cancels the first term in the square brackets, and the one from t = t1 combines with the second term in the square bracket; the ensuing equation simplifies to read Therefore, z(t1) is also a constant. The equation for z(t1), using FDT, becomes The two equations yield the low-temperature values z f = 2J and qin = 1 − T /J. In equilibrium C2(t, 0) and z(t1) remain constant and equal to their initial values, qin and z f .

The linear response in the frequency domain
Knowing that z(t1) remains constant in equilibrium, one can easily analyse the response equation in Fourier space. The equation that determines its dynamical evolution is transformed intô For this model Σ(ω) = J 2R (ω) and then R(ω), and also R(t), are independent of temperature for T < Tc = J while they depend on temperature through z f for T > Tc = J.
The terms under the square root can be more conveniently written as functions of the special values and the linear response is recast aŝ The imaginary and real parts ofR(ω) are then (note the unusual choice of sign for the imaginary part that we adopted.) In terms of the physical parameters, R(ω) is real for | − mω 2 + z f | > 2J. In the low temperature phase, since z f = 2J, this implies ω− = 0 and the imaginary part of the linear response is gapless. One can easily check that |R(ω)| 2 = 1/J 2 in the interval ω− ≤ ω ≤ ω+. Away from this interval the modulus of the linear response is a complicated function of the frequency. The zero frequency linear responsê with τ a time delay, must be a real quantity and this form is a manifestation of the condition z f ≥ 2J. The static susceptibility is then given by with the result in the second line being ensured by the choice of minus sign in front of the square root. The frequency-dependent linear response can then be transformed back to real-time and thus get its full timeevolution.
In the lower limit of the spectrum the imaginary part of the linear response goes as while in the upper limit it vanishes as with the corresponding ω± at T > Tc or T < Tc.

The correlation functions
The FDT in the frequency domain ImR(ω) = −βωĈ(ω) implies thatĈ(ω) should vanish in the same frequency intervals in which the linear response is real. In the ω = 0 case the linear response is real and the FDT as written above only imposes thatĈ(ω = 0) = ∞ −∞ dt C(t) must be finite. One has to bear in mind that in the cases in which the correlation function approaches a non-vanishing constant q asymptotically the Fourier transform to be computed is the one of Cst(t1 − t2) = C(t1 − t2) − q with respect to t1 − t2 and the integral should then yield Cst(ω = 0) = ∞ −∞ dt Cst(t) = 1 − q; more details are given in App. A. Consistently with the FDT constraints discussed in the previous paragraph, the time-dependent equation for Cst(t1 − t2) = C(t1 − t2) − qin can be treated following the steps explained in [9] and App. A.1.3 (that do not assume FDT) and one finds This equation has two solutions, either Cst(ω) = 0 or |R(ω)| 2 = 1/J 2 . The latter holds for frequencies in the interval ω ∈ [ω−, ω+]. Outside of this interval Cst(ω) must vanish. By taking the derivative of Eq. (76) with respect to t2 one readily checks that it equals Eq. (75) if the FDT between R and C is satisfied for all times and ∂t 2 C2(t2, 0) = 0 implying This last condition is a property of equilibrium as we have already discussed.
Having established that C2 is a constant, Eq. (77) enforces qin = 0, the high temperature initial value, or the low temperature Lagrange multiplier. Concerning the correlation function C(t1, 0), we write it as C(t1, 0) = Cst(t1, 0) + q0 allowing for a nonvanishing asymptotic value, q0, and taking Cst(t1, 0) such that it vanishes in the long t1 limit. The equation (76) evaluated at t2 = 0 is then rewritten as This equation has three terms that do not depend on Cst(t1, 0) explicitly and their sum should vanish in the long t1 limit. It trivially does for high temperature initial states since q0 = qin = 0 and, in equilibrium at low temperatures, we can assume q0 = qin, use FDT, and find confirming the assumption q0 = qin. The remaining equation fixes the time-dependence of C(t1, 0).

The energy before and after a quench
For the sake of completeness, we compute the energy variation due to a simultaneous change of the mass m0 → m and the variance of the interaction strengths J0 → J. In the applications and numerical tests we will focus on the latter changes only.
The kinetic energy density before the quench is the last equality being due to the fact that we take equilibrium initial conditions at temperature T . The potential energy density before the quench depends on the system being paramagnetic or condensed initially: The kinetic energy density right after the quench is and, since the velocities do not change in the infinitesimal interval taking from 0 − to 0 + , The post-quench potential energy density can be estimated from the relation between the Lagrange multiplier and the energy and the equation for z(t1) (see Sec. 4.1) Thus Assuming that the spin configuration did not change between the infinitesimal time step going from 0 to 0 + , All the values just derived imply the changes in the total energy respectively. We will concentrate on potential energy quenches only, and we will trace the phase diagram using the parameters x > 1 corresponds to energy extraction and x < 1 to energy injection. We will show that the parameter space is split in several sectors displaying fundamentally different dynamics.

Asymptotic analysis of the quench dynamics
We now study the full set of equations (75)-(78), derived in the N → ∞ limit, that couple the correlation C and linear response R functions. Using a number of hypotheses that we carefully list below, and that are not always satisfied by the actual evolution found with the numerical integration, we deduce some properties of the Lagrange multiplier, linear response and correlation function, in the long time limit. In this Section we state the assumptions, we summarise the results, and we leave most details of how these are derived to App. A.
Consider the system in equilibrium at T with parameters J0, m0 and let it evolve in isolation with parameters J, m. We will assume that the dynamics approach a long times limit in which one-time quantities, such as the Lagrange multiplier, reach a constant. Later we will further suppose that (in most cases) the correlation function becomes, itself, invariant under time-translations. Finally, we will explore in which circumstances a fluctuationdissipation theorem (FDT) can establish with respect to a temperature T f for all time-delays, in stationary cases, or for correlation values that are in the stationary regime, when we look for ageing solutions.
These assumptions are not obvious and, as we will show analytically in some cases and numerically in the next Section, do not apply to all quenches. Still, we find useful to explore their consequences and derive from them a set of relations between the control parameters for which special behaviour arises, that we will later put to the numerical test.

Asymptotic values in a steady state
Let us assume that the limiting value of the Lagrange multiplier is a constant The limit of the correlation function can be zero if the system decorrelates completely at sufficiently long time-delays or non-zero if it remains within a confined state. We therefore call the asymptotic value of the full two time correlation, or its stationary part in possible ageing cases, after the quench. Similarly, q0 = lim

The linear response function
The equation for the linear response function does not depend explicitly on the pre-quench parameters, it does only on the post-quench mass m and interaction strength J. (An implicit dependence on the initial state is not excluded, since the value taken by the Lagrange multiplier may depend on it.) The analysis that we developed for the constant energy dynamics applies to the sudden quench case too. The solution of the response equation in Fourier space yieldŝ with the zero frequency valueR From the numerical solution of the full equations that we will present in Sec. 6 we infer that for m = m0 and J = J0 with, we recall, x = J/J0 and y = T /J0. Replacing in Eq. (113) one notices that the minus sign has to be selected for x < y at low frequency and where we called χst, as a static susceptibility, the zero frequency response. After the quench ImR(ω) is non-zero in a finite interval of frequencies [ω−, ω+] with mω 2 ± = (z f ± 2J) as in Eq. (84) and z f taking the values in Eq. (114). Therefore, ImR(ω) is gapless for x > y and it is gapped for x < y.
These results are exact and do not assume anything apart from a long-time limit in which z f is timeindependent and given by Eq. (114). We have verified them with the complete numerical solution of the N → ∞ dynamic equations, see Sec. 6, on all sectors of the phase diagram. We can now Fourier back to real time to get the full time dependence of the linear response function. In the numerical Section we will compare this functional form, named Rst, to the outcome of the full integration of the dynamic equations.

The asymptotic kinetic and potential energies
From the relation between z and the energies, the conservation of energy, and Birkhoff's theorem, where the overlines represent a long-time average defined in Eq. (53). Note that the Lagrange multiplier takes the form of an action density, as a difference between kinetic and potential energy densities. Using these relations one derives the parameter dependence of the kinetic and potential energies in the four relevant regions of the phase diagram parametrised by x = J/J0 and y = T /J0.
The minimum potential energy density, epot = −J is realised for T = 0 in Sector III. We stress that we have found these results without using FDT and they can therefore hold out of thermal equilibrium. We will investigate later which other conditions impose the use of FDT, a strong Gibbs-Boltzmann equilibrium condition.

Kinetic temperature from the kinetic energy density
We can identify a kinetic temperature from the kinetic energy densities derived in the previous Subsection by simply imposing T kin = 2e kin . This operation leads to in the four Sectors of the phase diagram distinguished in the previous Subsection.

Final temperature under thermal equilibrium assumption
In App. A.3.1 we explain how we can exploit the conservation of the total energy, under the assumption that the asymptotic kinetic and potential energies take the form of Gibbs-Boltzmann equilibrium paramagnetic and condensed equilibrium phases at a single temperature T f to fix its value. This means that we require 2e kin = T f and with q = 1 − T f /J. For x > y we find The roman numbers between parenthesis refer to the cases listed in Eqs. (117)-(120) and to the four sectors in the phase diagram in Fig. 5. The difference between (IIa) and (IIb) is that in the first case we used a paramagnetic potential energy and in the second case a condensed one at

Limits of validity
The two bounds serve to find special curves on the phase diagram. The first bound is natural since we do not want to have a negative kinetic energy. The second one ensures that q ≥ 0. The implications of these bounds, that are examined They mean that an asymptotic state with a single temperature T f , or a double regime with temperature T f for C : 1 → q and T eff → ∞ for C : q → 0, could only exist below the piecewise curve y(x). One can simply check that the piece for y ≤ 1 lies in Sector IV, since y < x, and it is therefore irrelevant given that we have already shown that there cannot be a single temperature scenario in this Sector. The limit then moves to x = y for y < 1. Parameters on the special curve y = √ x will play a special role, as we will show below.

Particular values
For the moment, a single temperature scenario for the global observables in the N → ∞ limit seems possible for y < 1 and x > y (III), and below the curve y = √ x for y > 1 in II. It is instructive to work out the limiting values of T f /J0 and q = 1 − T f /J on the borders of the region y < 1 and y < x (III) of the phase diagram, and the no-quench case x = 1: These values match at x = y = 1. q is larger than zero on the lines x = y and y = 1 that mark the end of what we call Sector III. Moreover, the approximate asymptotic analysis of the mode dynamics of the finite N system will lead to Tµ = T f in Sector III, see Eq. (183). On the limiting curve in Sector II, y = √ x for y > 1, T f = J and q = 0. As one could have intuitively expected, T f > T for energy injection quenches (x < 1) and T f < T for energy extraction quenches (x > 1).

Results under FDT at a single temperature
We now add one assumption to the analysis: that the FDT, at the single temperature T f , relates the linear response to the correlation

Static susceptibility
The values of the zero frequency linear responsê where one must recall that τ is a time-difference, forcê The result in the last line, valid for (IIb) and (III), does not put any additional constraint on T f . Instead, the other conditions in the first two lines are incompatible with the expressions imposed by the energetic considerations. They corroborate the impossibility of having a single T f in (I) and (IV) or with q = 0 in (II).

Limits of validity of the single temperature scenario
As explained above and in App. A.3.3, the consistency between the static susceptibility values and the T f derived from conservation of energy impose that FDT with a single temperature may only hold for x > y and y < 1 (sector III) or below the special curve y = √ x, for y > 1 (lying inside sector II). Whether this is realised or not needs to be investigated numerically.

Two step (possibly ageing) Ansatz
One can also look for a two step solution with similar characteristics to the ageing one found for dissipative dynamics [44] and summarised in Sec. 2.3. Asking for the relation between correlation and linear response in Eq. (127) to hold in a stationary regime of relaxation in which C decays from 1 to q, and that the effective temperature T eff characterising the second regime of decay from q to 0 diverges, one recovers with the same T f and q = 0 as in Eq. (123).

The correlation function
From the analysis of the equation ruling the two-time correlation function, assuming stationarity and hence a dependence on time-difference only, one deduces (the details of the derivation are given in App. A.1.3, see also [9] for a general treatment)Ĉ whereĈst(ω) is the Fourier transform of the time-varying part (subtracting the possible non-vanishing asymptotic value q). Note that the relation between the correlation and the linear response is the same as the one that we derived in the constant energy no-quench problem. It implieŝ independently of the control parameters. Below we check numerically that these relations are satisfied in various quenches. In particular, from the analytic form ofR(ω) one can easily see that |R(ω)| 2 = 1/J 2 in the frequency interval in which the linear response is complex. Importantly enough, we cannot yield an explicit analytic form ofĈ(ω) since it is factorised on both sides of the identity (131). We are forced to go back to the full set of dynamic equations and solve them numerically to get insight on the behaviour of C(t1, t2).

Numerical results preview
We will see that a state with a single T f characterising the fluctuation-dissipation relation is reached numerically in the following two cases only: (1) the dynamics are run at constant parameters (no quench), x = 1; (2) the special relation y = √ x between parameters holds and y > 1 (within sector II). In all other cases no equilibrium results à la Gibbs-Boltzmann are found for the global observables (correlation functions, linear response functions, kinetic and potential energies) but a different statistical description, of a generalised kind, should be adopted. In particular, in Sector III where the conditions derived from the energy conservation and static susceptibility allowed for a single temperature scenario, the full solution of the complete set of questions will prove that this is not realised. A detailed explanation is given in the numerical section of the paper and in the analysis of the N − 1 integrals of motion presented in Sec. 7.
We recall that the p ≥ 3 strongly interacting case behaves very differently [7]. On the one hand, equilibrium towards a proper paramagnetic state, and within confining metastable states, were reached in two sectors of its dynamic phase diagram. On the other hand, an ageing asymptotic state in a tuned regime of parameters was also found for more than two spin interactions in the potential energy. In the p = 2 integrable model we do not find an ageing asymptotic state. Moreover, Gibbs-Boltzmann equilibrium is achieved in the two very particular cases listed above only.

Analytic results for the dynamics of the finite size system
In this Section we describe how the finite system size dynamics can be solved by using a convenient basis in which the evolution becomes the one of harmonic oscillators coupled only through the Lagrange multiplier. We show that these oscillators decouple under the assumption z(t) → z f allowing for a simple approximate solution of the problem that can, however, be relevant for N → ∞ only. We then explain a way to numerically solve the dynamics for finite N .

Newton equations in a rotated basis
Take a system with finite N . The post-quench matrix Jij has µ = 1, ..., N eigenmodes with eigenvalues λµ. If we denote the projection of the spin vector in the direction of the µ-th eigenvector of Jij, the N rotated equations of motion read This set of equations has to be complemented with the initial conditions sµ(0) andṡµ(0). They are very similar to the equations for a parametric oscillator, the difference being that, in our case, the time-dependent frequency depends on the variables via the Lagrange multiplier. Furthermore, they are identical to the equations of the Neumann's integrable classical system [15], see Sec. 3.2.
Once the equations of motion for the sµ are solved, we can recover the trajectories for s using s(t) = µ sµ(t) vµ. In particular, the correlation function is given by where the subscript J means that the result depends, in principle, on the interaction matrix chosen, and the angular brackets represent an average over initial conditions. One could then perform the disorder average or analyse the self-averageness properties of the correlation in different time regimes. The J dependence should disappear in the N → ∞ limit. Since we are interested in an uniform interaction quench, it is easy to see that

Behaviour under stationary conditions
Let us assume that the system reaches stationarity and that the Lagrange multiplier approaches a constant In order to simplify the notation, in the rest of this section we will measure time with respect to a reference time tst for which the stationary regime for z(t) has already been established. (We insist upon the fact that this assumption can only hold for N → ∞.) The equation of motion of each mode becomes and can be thought of as Newton's equation for the mode Hamiltonian with ω 2 µ ≡ (z f − λµ)/m. This equation has three types of solutions depending on the sign of ω 2 µ : that is to say, oscillatory solutions with constant amplitude in the first case, diffusive behaviour in the intermediate case and exponentially diverging solutions in the last case. We insist upon the fact that the initial time here is taken to be tst the time needed to reach the stationary state.
If the Lagrange multiplier approaches, then, a value that is larger than λN , all modes oscillate indefinitely. In Gibbs-Boltzmann equilibrium in the PM phase, zeq > λN and such a fully oscillating behaviour is expected. In equilibrium in the low temperature condensed phase zeq = λmax = λN for N → ∞ and the µ = N mode should grow linearly in time while all other modes should oscillate with frequency ωµ = (z f − λµ)/m. The amplitude of each mode is determined by the initial conditions, that are actually matching conditions at time tst, when stationarity is reached in this case. (Recall that λN−1 is at distance N −2/3 from λN [23]. This means that, under the assumption z f → λN , z f − λN−1 = λN − λN−1 N −2/3 and there will be almost diffusive modes close to the largest one in the large N limit.) However, the simulations at finite N show that for finite N , z f is always greater than λN and all modes are oscillatory. For "condensed-type" dynamics z f will still be greater than λN , although very close to it. The diffusive behaviour of the N th mode (in the N → ∞ limit) would be obtained as the limit of zero frequency of a (finite N ) oscillating N th mode.

Mode observables
At variance with the N → ∞ approach, the finite size study allows to access the details of the dynamics of each mode. In this section we define some mode-observables that will provide valuable information. Of particular interest are the mode energies, which can be defined as kin Note that in the analysis of the N → ∞ model the potential energy density is epot = −1/(2N ) µ λµ s 2 µ without the term proportional to z(t). For this reason we use here the different symbol pot µ for the mode potential energies that include the term proportional to z(t). The values of these energies at t = 0 − are given by the fact that all modes are in equilibrium at the same temperature: Immediately after the quench, they are In order to study the eventual thermalisation of the system, we can define an effective time dependent mode temperature through the total mode energy based on the fact that the modes are (quasi) decoupled. Whenever the system enters a stationary regime in which z(t) is constant, see Section 5.2, the mode temperatures Tµ are independent of time, since the system behaves as a collection of non-coupled harmonic oscillators. We can also define mode temperatures using the kinetic and potential mode energies that oscillate around their mean if we take their time-average [67] e kin,pot and propose T kin,pot In the stationary regime, as shown in Section 5.2, the different mode temperatures should verify Other useful observables are the time-delayed mode correlation functions and the mode linear response functions that is defined and measured as follows. If we add an external field hµ linearly coupled to each mode sµ, the equations of motion are modified into and its solution reads where s (hom) µ (t) is the solution to the Newton equation without the external field and s P µ (t) is a particular solution of the inhomogeneous problem with initial condition s P µ (0) = 0 andṡ P µ (0) = 0. The linear response function Rµ(t1, t2) ≡ δsµ(t1)/δhµ(t2)| h=0 of the mode µ can be defined equivalently through In practice, to measure the linear response function numerically we apply a small external field localised in time hµ(t) = h0δ(t − t2) and we solve the inhomogeneous problem to obtain s (inhom) µ (t). Using Eq. (152) we obtain the linear response as where s (hom) µ (t1) must have been calculated independently. In thermal equilibrium the linear response and correlation function are related by the fluctuation dissipation relation, Whether the time-evolving correlation and linear response satisfy this relation, whether the mode temperatures are the same as the ones obtained from the energy characteristics of the modes and, finally, whether they all take the same value, are issues that we will explore.

Kinetic and potential mode energies in the stationary state
As already mentioned, in a steady state, z(t) → z f , the modes kinetic and potential energies are Clearly, neither the kinetic nor the potential mode energies are constant, but, in the steady state limit the sum of the two, that is to say the total mode energy, is with tst the time at which the steady state is established and t > tst.
As expected from Birkhoff's theorem [67], the long-time averages, say taken after tst, should be constants and one can expect them to be equal to half the total energy (In practice, the average over a few periods is enough to obtain the constant value.) If one now associates a temperature to these values, arguing equipartition of quadratic degrees of freedom, one has The mode temperatures depend on the averages at the end of the transient, and the mode frequency Ω 2 fµ that itself depends on the asymptotic limit of the Lagrange multiplier z f and the eigenvalue λµ.
In the argument above we implicitly assumed that ω 2 µ does not vanish. The case µ = N is tricky. If one naively sets ω 2 µ to zero from the outset 2 pot N = ω 2 µ s 2 N apparently vanishes. The correct way of treating the largest mode is to remember that the projection on the largest mode condenses and that s 2 N is proportional to N . This will ensure that limN 1 s 2 N ∝ N , in such a way that limN→∞ ω 2 µ s 2 N = 2 kin µ , similarly to what happens in equilibrium, where s 2 N = qinN and the Lagrange multiplier is such that (z f − λN )qinN = T . We will see in the next Sections that, in some cases, the scenario described in this Section is actually realised by the dynamics. Which are the quenches in which such a behaviour is observed will be determined by the complete solution of Newton's equations with the methods that we will now describe.

Initial conditions: equilibrium averages with finite N
In this section we address the calculation of equilibrium averages at finite N in order to provide suitable initial conditions for the numerical integration of the mode dynamics explained in Sec. 5.8.
If we were to naively integrate the mode equations, we would need to draw initial vectors, s(0) = (s1, . . . , sN ) and˙ s(0) = (ṡ1, . . . ,ṡN ), mimicking an initial thermal state at finite temperature, be it T > Tc = J0 or T < Tc = J0, for a given realisation of the N × N interaction matrix. Averages over these initial states of the interesting observables should then be computed. This method is computationally expensive as a large number of initial state should be considered to get smooth and reliable results. Instead, the numerical method that we will explain in Sec. 5.8 is such that only the averages s 2 µ eq and ṡ 2 µ eq are needed as input for the initial conditions. We then focus on determining these averages in a finite size system in equilibrium.
The canonical equilibrium probability density of the configuration {pµ = mṡµ, sµ} at temperature T , for a given realization of disorder, is with Z the partition function. The statistical averages are computed as integrals over this measure. The integrals over pµ range from −∞ to ∞. The quadratic averages of the velocities are thus simply given by just as for the infinite N case, and the initial conditions will be ṡ 2 µ (0 + ) = ṡ 2 µ eq.
As long as the equilibrium value of the Lagrange multiplier be strictly larger than the maximum eigenvalue, the weight of the coordinates sµ are well-defined independent Gaussian factors. We will see that the self-consistent solution complies with this bound. Relying on the spherical constraint being imposed by the Lagrange multiplier, we extend the sµ integrals to ±∞ and The difference between the two equilibrium phases will be codified in the value of z (N ) eq , which can be obtained as the solution of the spherical constraint equation We solved this equation numerically to determine z (N ) eq and we found that the solution turns out to be always greater than λ is always larger than λ (0) max and, as expected, the difference between them decreases with system size.
Once the finite size Lagrange multiplier is obtained, we replace it in Eq. (162) to obtain the initial conditions s 2 µ (0 + ) = s 2 µ eq for the mode dynamics.  To gain insight into the scaling with the system size, in Fig. 1 (b) we plot the difference between z (N ) eq and λ (0) max for temperatures in the condensed phase as a function of 1/N . The straight dashed lines have slope T /qin, where qin = 1 − T /J0 is the N → ∞ value of the self-overlap. For temperatures sufficiently below the transition, the finite size data, obtained for one particular realisation of the random matrix Jij, follow the infinite size results for all system sizes analysed. For temperatures close to the transition, there appear deviations for the smallest system sizes (largest 1/N ). In conclusion, we find that for large system sizes or temperatures not too close to the transition, the solution to Eq. (163) behaves as Based on this, we define a finite size version of the equilibrium self-overlap which is finite if the highest mode is macroscopically populated. For N → ∞, qin = 1 − T /J0 for T < J0 and zero for T > J0. In Fig. 2 (a) we show q (N ) in as a function of temperature. We can observe the convergence of the finite size results towards the N → ∞ predictions as the system size is increased.  Finally, we investigate the finite N corrections to z (N ) eq in the paramagnetic phase, T > J0. We find that a linear scaling in 1/N also applies here, but the value of z where s(T ) = s is the slope of the dashed lines that we obtained from a fit and turns out to be independent of temperature (all the dashed curves are parallel straight lines). Using the definition in Eq. (165), we can express the Lagrange multiplier as z in ), and we can verify that the potential energy of the highest mode (note that we included the term proportional to the Lagrange multiplier and we therefore compute pot N instead of e pot N ) assumes the correct value in equilibrium, i.e., the one consistent with the equipartition theorem, if

Energy variation at the quench
We here compute the finite-N equilibrium values of the Lagrange multiplier, kinetic and potential mode energies using the finite size averages for s 2 µ eq and ṡ 2 µ eq proposed in the previous Section, and we compare them with the equilibrium N → ∞ results obtained in Sec. 2.2.1. Analogously to the N → ∞ results in Sec. 4.3, these equilibrium values define the initial condition for the interaction quench and, therefore, set the values of the observables at t = 0 + .

Pre-quench energies
We begin with the kinetic energy right before the quench, e (N ) It coincides with the infinite-N mean-field result.
Next we analyze the pre-quench potential energy, where we used Eq. (163). In the equilibrium condensed phase we can rely on z The and one readily recovers the N → ∞ limit epot(0 − ) = −J 2 0 /(2T ).

Post-quench energies
Now we will compute the values of the kinetic and potential energy, and the Lagrange multiplier after an interaction quench The kinetic energy is not affected by the quench in the interaction and, just as in the N → ∞ limit (see Sec. 4.3), we have that e (N ) For the potential energy it is enough to note that e (N ) Using that z (N ) (0 + ) = 2(e For initial states in the paramagnetic phase

Independent harmonic oscillators in the asymptotic limit
We now use the results in App. B concerning the quench dynamics of a harmonic oscillator in the context of our non trivial problem. In equilibrium at time t = 0 − the initial frequencies of the modes are In the asymptotic limit after the quench we identify the frequencies with where we assumed that the Lagrange multiplier reached the constant z f = limt→∞ z(t).
The analysis of the harmonic oscillator does not need any long-time assumption to set its spring constant, or frequency, to a constant value. In our problem, the dynamics may approach the ones of independent harmonic oscillators with constant spring constants only asymptotically. During the transient evolution the mode energies vary. In reality, we do not know the values they take at the end of the transient regime. We can make a rough approximation in which we assume that z f is reached instantaneously after the quench, z(0 + ) = z f , so that we can use instead of the unknown values at the end of the stationary regime. Under these assumptions the final mode temperatures are They are for y > x and y > 1 (I) for y > x and y < 1 (IV) Several comments are in order. The expression for x > y and y < 1 (sector III) is the same as the one that we derived from the analysis of the N → ∞ Schwinger-Dyson equations, see Eq. (123), simply T f kin . The no-quench case x = 1 in realised in Sectors I and III and one rapidly checks that T f µ = T in both cases. On the curve y = √ x the mode temperature do not take the same value. We will argue later that the approximation used in the Section yields a qualitatively erroneous result in this case. Continuity between sectors I and IV on the one side, and II and III on the other, are ensured setting y = 1 that is to say T = λ is also verified.
We would like to know which is the condition satisfied by z f under this approximation. In order to obtain such equation, we first note that the time-dependent spherical constraint imposes that N µ=1 s 2 µ (t) = N . (185) In particular, this implies that Inserting the approximation in Eq. (180) in the time-averaged spherical constraint, we find an equation for z Since z is chosen in such a way that we find that the equation for z In other words, under this approximation, z (N ) f is the equilibrium Lagrange multiplier for a system in equilibrium at temperature T with variance of the disorder distribution equal to J. In the We will put these predictions to the test in Sec. 6 using the numerical solution to the finite N evolution equations with the numerical method that we describe in the next Subsection. In various regions of the phase diagram these a priori approximate forms are in strikingly good agreement with the numerical data. In others they are not and we discuss why this is so.

Exact solution of the mode dynamics
One possible approach to solve the dynamics of each mode starting from canonical equilibrium initial conditions is to take a large ensemble of initial configurations drawn from the Gibbs-Boltzmann distribution, numerically integrate the Newton equations Eq. (134) for each initial condition, and then calculate the observables averaging over the trajectories corresponding to the different initial states. Such approach is feasible but computationally very demanding. In this Section we describe a more convenient method to solve the dynamics for each mode that uses heavily the tools developed to treat a paradigmatic problem in classical mechanics, the one of parametric oscillators [77,78,79].
In order to solve Eq. (134) we propose an amplitude-phase Ansatz [77,78,79,80] Inserting this Ansatz in the µth mode Newton equation, we obtain an equation for the mode and time dependent auxiliary function Ωµ(t), 1 2Ω where ω 2 µ (t) ≡ (z(t) − λµ)/m. The last equation has to be complemented by the initial values ωµ(0) andΩµ(0). If we chooseΩ we find that the projection of the spin configuration is which is reminiscent of the general solution of the harmonic oscillator problem, see Eq. (140), here with a time-dependent "frequency" Ωµ(t).
We still have to specify the initial condition for Ωµ(0). A possible choice is that enforcesΩµ(0) = 0 [80]. However, this choice is consistent with real Ωµ(t) only if z(0) − λµ ≥ 0, which is verified uniquely for J ≤ J0, i.e., uniquely for energy injection. An initial condition ensuring real and positive Ωµ(t) for all µ for any quench is We choose this initial condition for the numerical calculations.
In order to solve for Ωµ(t) we consider the equal-times mode correlation function in terms of which we write the potential energy as Replacing this equation in z(t) = 2e f − 4epot(t), we find an expression of the Lagrange multiplier as a function of the mode correlations at equal times Finally, we note that the system conformed by Eqs. (191), (196) and (198) is closed and allows to find the time evolution of the Lagrange multiplier and the auxiliary functions Ωµ(t). This set of equations is amenable to numerical integration. Once we obtain Ωµ(t), the most interesting observables can be calculated using the general solution in the form given in Eq. (190). The advantage of this method is that we do not need to draw initial states {sµ(0),ṡµ(0)} but we only have to specify the initial averages s 2 µ (0) and ṡ 2 µ (0) that we will take to be the ones enforced by equilibrium at T , that is to say, the forms given in Eqs. (160) and (162).

Numerical results
This Section summarises what we found numerically by solving the N → ∞ Schwinger-Dyson equations that couple the global correlation and linear response C and R (see Sec. 4), and the finite N ones acting on the mode projections sµ (see Sec. 5). Some general considerations about the numerical algorithm used to integrate the N → ∞ equations are given in App. C.
The finite N results are consistent with the infinite N ones and help us understanding the mechanism whereby the dynamics take place. We chose to start this Section with the summary of the dynamical phase diagram and the behaviour of the quantities that determine it. Later, we give further details on the dynamics at constant energy (no quench) and in each of the dynamic phases identified after sudden quenches.
We signal here that we will make a special effort to show, in each case considered, that an asymptotic state characterised by the single temperature T f that the naive asymptotic analysis of the dynamic equations predicts is not attained. The investigations that lead to this conclusion are very instructive not only because they prove the lack of Gibbs-Boltzmann equilibrium but also because they lead to the evaluation of the mode temperatures that will play a role in the statistical description of the steady states.

The phase diagram
The phase diagram is determined by the asymptotic behaviour of the zero frequency linear response or susceptibility, χst =R(ω = 0), and the asymptotic value of the Lagrange multiplier. We determine their values through the variation of the parameter J/J0 for fixed T /J0. In the phase diagram presented in Fig. 5 and the ensuing discussion we call y = T /J0 the vertical axis and x = J/J0 the horizontal one. The former determines the initial state and the latter the kind of quench performed with injection of energy for x < 1 and extraction of energy for x > 1.
We study separately parameters in four Sectors of the phase diagram, although the final results will allow us to distinguish three different phases. The Sectors are indicated with Roman numbers and the phases with different colours or shades in Fig. 5. We also mark the line x = 1 (equilibrium dynamics) and the curve y = √ x with y > 1 where special dynamics are found.
We recall that dynamic phase transitions have been found in the quench dynamics of quantum isolated systems, see, e.g. [81,82,56,57,58,83,84,61]. Here, and in [7], we see dynamic phase transitions arise in the Newtonian dynamics of isolated classical interacting systems.
In Fig. 3 (a) we check the prediction (115) for the zero frequency linear response function. We plot T R (ω = 0) against J/T and we see the change in behaviour fromR(ω = 0) = 1/T toR(ω = 0) = 1/J at xc(y) = y, that is T = J. The change inR(ω = 0) is accompanied by a change in the asymptotic value of z as a function of the quench parameter x = J/J0. This fact is confirmed numerically in Fig. 4 where data for N → ∞ and N finite are shown in panels (a) and (b), respectively. For x < y, the numerically estimated z f (x) for fixed y in the case N → ∞ (a) were fitted with the polynomial function f (x) = a x 2 + b x + c. We obtained very good results with a 1/y, b 0 and c y (the precision of the fit is really very high in terms of reduced χ 2 ). These results strongly suggest the following functional dependence of z f on the parameters x and y, To get a visual confirmation of this argument, in Fig. 4 (a) we plotted the functions f (x) = x 2 /y + y for x < y, one for each one of the values of y that the numerical data refer to. The agreement between the data and the prediction is very good. The analysis of the finite N data was done along the same lines, see Fig. 4 (b), with the difference that the data for x < y were fitted by T + λ 2 max /(4T ) and the ones for x > y with Jλmax finding again very good agreement. (We found an appreciable deviation in the fit for x > y had we used 2J instead of Jλmax. Regarding the results for x < y we could have used T + J 2 /T with a similar quality for the fit.) Remarkably, the functional dependence proposed in Eq. (199) is the one predicted by the independent harmonic oscillators approximation in Eq. (189).
By using the change in the dependence ofR(ω = 0) and z f on the quench parameter x = J/J0 as a criterium to track the dynamical phase transition, we obtained the numerical estimates of the critical curve xc(y). In Fig. 5 we show the data for xc(y) (with error bars) for some values of the control parameter y. The data strongly suggest that there is a linear relation between the critical value xc and the parameter y for any y; in short, we confirm xc(y) = y.
Concerning the long-time behaviour of C(t1, t2), it is useful to distinguish the cases y < 1 and y > 1, that is, quenches that start from equilibrium in the condensed phase from quenches that start from equilibrium in the paramagnetic phase.
We observe the following trends: • For x < xc(y), C(t1, t2) tends to be stationary, though within the time scales of the numerics it has not reached this limit yet when y is too small. In most instances, C(t1, t2) oscillates around 0, exceptions being the critical quench and the case with both y > 1 and x > 1 where zero is approached asymptotically from below. The time average of C computed on intervals far from the initial transitory regime vanishes in all cases suggesting an effective q = 0.
• For x ≥ xc(y), C(t1, t2) is rapidly stationary and one very clearly identifies the asymptotic constant q = lim t 2 t 0 t 1 −t 2 t 0 C(t1, t2). For y < 1 it is different from zero while for y > 1 it equals zero. The asymptotic q0 = limt 2 t 0 C(t1, 0) is different from q in the cases in which both are non-vanishing.
These facts can be appreciated in Fig. 6 where we display the decay of the correlation function for several choices of the parameters in different regions of the phase diagram, marked with crosses in Fig. 5.
Having announced the main features of the dynamics after different types of quenches, in the rest of this Section we will support these claims with the detailed study of all relevant observables. x for x > 1 and g(x) = 2x/(1 + x) for x ≤ 1, respectively, discussed in the text. The three phases that are characterised by different behaviour of z f , χst and the long times limit of C(t1, t2) are highlighted with different colors. The red data points equipped with error bars indicate the numerical estimate of xc for several values of y = T /J0 (see the main text for more details). The crosses indicate the cases shown in Fig. 6 with the corresponding labels.   Fig. 5. In the first row x ≤ y while in the second x > y. The first three panels in both rows correspond to y < 1 and the last two to y > 1. The third panel in the first row is on the critical line x = y.

Constant energy dynamics
We first checked that for J = J0 and m = m0, that is to say ∆e = 0, and for all initial conditions, the system has stationary evolution and the total energy as well as other conserved quantities are indeed conserved. As the dynamics of generic Hamiltonian systems is hard to control numerically, we included this analysis to validate our algorithms. Moreover, this study allowed us to know which is the order of the numerical error incurred into. A time discretisation step δ = 0.001 in the integration of the N → ∞ Schwinger-Dyson equations was sufficient to assure numerical convergence of our results. In the integration of the finite N problem we found a weak dependence on δ but the value δ = 0.001 gave acceptable results.
We studied the equilibrium dynamics for initial paramagnetic configurations (T > J0) and condensed ones (T < J0). We present and discuss the results in two Subsections. As already said, they serve as benchmarks for the more interesting quenching cases that we put forward later.

Dynamics in the paramagnetic phase
In this Section we use parameters in the paramagnetic (PM) equilibrium phase. We fix J0 = m0 = m = 1, we equilibrate the initial conditions at T = 1.25 and we evolve with parameters J = J0 = 1 and m = m0 = 1 (no quench). In Figs. 7 and 8 we show the numerical results for finite N and infinite N , respectively. The system should remain in Gibbs-Boltzmann equilibrium in the PM phase in both cases.
In Fig. 7 we analyse the results of the numerical integration of the N → ∞ Schwinger-Dyson equations for T = 1.25 and all other parameters set to 1. The figure shows the time evolution of the correlation function (a), the fluctuation dissipation parametric plot (b), the response function, R(t1, t2) together with −(1/T ) ∂t 1 C(t1, t2) (c), the deviation of the Lagrange multiplier from its equilibrium value (d), the potential and kinetic contributions to the energy density (e), and two studies of the fluctuation-dissipation theorem in the frequency domain (f) and (g).
The correlation with the initial condition, C(t1, 0), and the ones between two different times, C(t1, t2), are identical, meaning that time translation invariance is satisfied. All the correlations relax to zero q = q0 = 0. Moreover, the curves coincide approximately with the ones in Fig. 8 (a) which were obtained by integrating Newton equations for finite N .
The fluctuation-dissipation relation [68] is satisfied with the temperature of the initial condition, that is the same as the one of the final state. This fact can be proven in general for Newtonian evolution of initial configurations drawn from Gibbs-Boltzmann equilibrium [85]. Indeed in panel (b) we display the parametric plot χ(t1, t2) = t 1 t 2 dt R(t1, t ) vs C(t1, t2) for two waiting times t2 and, with a dashed line, the equilibrium result −1/T C finding perfect agreement within our numerical accuracy. As a further confirmation of the validity of the fluctuation-dissipation theorem, we show the response function, R(t1, t2), and −(1/T ) ∂t 1 C(t1, t2), both plotted against t1 − t2, for two different values of t2 in panel (c). The two quantities coincide almost perfectly. In the same panel we checked that the response function coincides with the one derived by taking the inverse of the theoretical Fourier transform of the stationary asymptotic response, given by Eq. (83). The theoretical curve is indicated as Rst.
The Lagrange multiplier and the potential and kinetic energies remain constant throughout the evolution of the system and equal to their predicted values, apart from small deviations due to the numerical errors introduced by the numerical integration scheme, see panels (d) and (f), and App. C. The plot showing z(t) proves that the relative error in this quantity is at most of order 10 −7 . All these results are compatible with Gibbs-Boltzmann equilibrium in the paramagnetic phase. We do not show the time evolution of the off-diagonal correlation with the initial configuration, C2(t, 0), since it is identically zero at all times.
In panel (f) we show the Fourier transforms of the correlation and response functions,Ĉ(ω, t2) andR(ω, t2) respectively (the transform is performed on the variable τ = t1 − t2 with t2 fixed), for two different values of t2. Note that we are showing only the real part ofĈ(ω, t2), since we implicitly assume that C(t2+τ, t2) = C(t2−τ, t2). The black solid lines represent the theoretical prediction for the real and imaginary parts of the Fourier transform of the response function in the stationary regime,Rst(ω), given by Eq. (83), the inverse Fourier transform of which is plotted in (c). In panel (g), the ratio −ImR(ω)/(ωĈ(ω)) together with the prediction 1/T from FDT indicated by a horizontal dashed line are shown. Note the deviation from the flat result at the right edge of the frequency spectrum. This is due to the fact that the ratio approaches zero over zero and the numerical error incurred for those large frequencies is much amplified. At the left end of the spectrum, the more interesting low frequency regime, the oscillations are only present for the t = 0 curve.
In Fig. 8  . We recall that we chose to use a convention such that the imaginary part ofR is negative. (g) The ratio −ImR(ω)/(ωĈ(ω)) together with 1/T indicated by a dashed horizontal line. method explained in Sec. 5.8. In panel (a) we see that the correlations with the initial condition quickly relax to 0, as expected in the PM phase. They do with a weak size-dependence in the long time-delay tails. We only show the correlation with the initial configuration since we have checked that the time-delayed one is stationary. Also included in this panel is the same correlation function computed using the Schwinger-Dyson equation valid in the N → ∞ limit. We see perfect agreement with the finite N results at short times and small deviations at longer times. In Fig. 8 (b) we observe that the global kinetic, potential and total energy densities are constant, as expected. The Lagrange multiplier is studied in panel (c) where we plot it subtracting z The mode-by-mode analysis of the finite N dynamics is performed in Fig. 8 (d) and (e). Two panels display the time-delay dependence of the correlation function of the first and last mode. In Fig. 8 (f) we display the mode temperatures Tµ(t) at the initial time and after a long time evolution. The mode temperatures coincide with the expected equilibrium value, except for the largest modes, where there is a very small deviation. These variations represent small numerical errors due to the finite time-step discretisation used numerically, and are hard to improve algorithmically unless by using a still smaller integration step. We have also checked (not shown) that the mode correlations Cµ(t1, t2) and the mode response function Rµ(t1, t2) satisfy the fluctuation-dissipation relation with a temperature given by T for all modes.

Dynamics in the condensed phase
We now turn to the constant energy dynamics in the condensed, low temperature equilibrium phase. In Fig. 9 we show the mode dynamics for initial conditions in equilibrium at T = 0.5. From Fig. 9 (b) we notice that the total energy is conserved and that the kinetic and potential contributions are also constant, consistent with thermal equilibrium in the isolated system. In Fig. 9 (c) we show the Lagrange multiplier, which should be constant in equilibrium. In the numerical solution, the Lagrange multiplier exhibits oscillations around the initial value. Their amplitude decreases consistently with the integration step δ, implying that for δ → 0 we recover the expected constant behaviour. The two-time global correlation C(t1, t2) is stationary for all system sizes (not shown), so we focus on the particular case with t2 = 0. We can see from Fig. 9 (a) that, at variance with the paramagnetic case, the dynamics of the correlation function C(t1, 0) has a strong dependence on the system size. After a fast decay from the initial value, the correlation shows a plateau, the lifetime of which increases with system size, approaching asymptotically the value predicted by the N → ∞ treatment that is shown with a (red) horizontal line. The source of this size dependence is the behaviour of the largest mode µ = N . In Fig. 9 (f) we show the time dependence of the largest mode correlation function CN (t, 0) for different system sizes. We observe that its oscillation frequency decreases as we increase the system size. A similar finite size effect is seen in the dynamics of the next-to-largest mode in panel (g). Since the largest modes dominate the long-time dynamics, this effect causes the size dependence of the plateau lifetime. For N → ∞ the oscillation frequency of the N -th mode goes to zero, allowing for the presence of an infinite plateau, see Fig. 11. The modes lying in the middle and other end of the spectrum have almost no size dependence, as shown in (d). Figure 10 investigates the mode kinetic and potential energies and the mode temperatures that can be extracted from them. In Fig. 10 (c) we show the mode temperatures Tµ(t) at two measuring times, as a function of the mode index, what we will call temperature spectrum later. We observe deviations from the expected behaviour T tot µ = T ∀µ only close to µ = N . To gain insight into the mode temperature deviations, in panels (a) and (b) we show the time evolution of the kinetic, potential and total energies for the first and last modes. For the lowest modes, the kinetic and potential energies are constant, equal, and their sum is identical to T , consistently with equilibrium dynamics. However, for the modes that are closer to µ = N the energies weakly oscillate with an amplitude that decreases with the integration step and should vanish for δ → 0. This implies that in such limit, all mode temperatures should be equal to the equilibrium temperature, even for modes close to the right edge of the spectrum.
We now turn to the analysis of the dynamics in the N → ∞ limit. In Fig. 11 we show the results obtained through the numerical integration of the Schwinger-Dyson equations for the two-time correlation and response functions, and the same choice of parameters, that is T = 0.5 and J = 1. Stationarity is satisfied as well as the FDT with the initial temperature, see (b) and (c). The main difference with the case in Fig. 7 is that the correlation functions, both with the initial condition and with the configuration at a waiting-time t2, relax to a non-vanishing value (a). Within numerical accuracy we observe q0 = q 1 − T /J = 0.5 and this value as well as the potential and kinetic energies (not shown) are consistent with equilibrium at T . Also in this case, we checked that the response R(t2 + τ, t2), for t2 ≥ 0, coincides with the (numerical) inverse Fourier transform of the theoretical prediction given by Eq. (83) for the Fourier transform of the stationary asymptotic R, see panel (c).
In panel (d) we show the time-evolution of the off-diagonal correlation, C2(t, 0). In the case of equilibrium dynamics, C2(t, 0) should be constant and equal to qin. As one can see, the value of C2(t, 0) obtained by numerical integration is not exactly a constant function, but it approaches a constant in the long time limit which differs from q = qin = 0.5 by a very small amount. This deviation is only due to the approximations introduced by the numerical integration scheme. The same can be said about the behaviour of the numerical z(t), see panel (e). Its value oscillates around the expected equilibrium value zeq = 2J = 2 at T = 0.5, with oscillations amplitude of order 10 −7 for short times and decreasing with time.
We next show the Fourier transforms of the correlation and response functions, for two different values of t2 in Fig. 7 (f). Again, note that we are showing only the real part ofĈ(ω, t2), since we implicitly assume that C(t2 + τ, t2) = C(t2 − τ, t2). The black solid lines represent the theoretical prediction for the real and imaginary parts of the Fourier transform of the response function in the stationary regime,Rst(ω), given by Eq. (83), the inverse transform of which is plotted in (c). In panel (g) we display the ratio −ImR(ω)/(ωĈ(ω)) together with the prediction 1/T from FDT indicated by a dashed horizontal line. In all presentations we find good agreement with the validity of FDT with the proviso that in the plot in (g) the high frequency regime is contaminated by the numerical error, and the low frequency regime by the fact that we can perform the Fourier transform on a finite time window only, and this causes the dependence on t2 shown in the plot. (c) 25 15 χ C t 2 = 11.25 15 T ′ = 0.50  83), the inverse transform of which is plotted in (c). In panel (g), the ratio −ImR(ω)/(ωĈ(ω)) together with the prediction 1/T from FDT plotted with a dashed horizontal line.

Instantaneous quenches
We shall now vary the initial temperature T and the coupling J of the Hamiltonian that drives the time evolution to consider specific instantaneous quenching processes that inject or extract energy. The aim is to illustrate the analytical results of the previous Section and put them to the test. As we have already announced the structure of the phase diagram, we will consider representative quenches in the four sectors, labeled I (y > 1 and y > x), II (y > 1 and y < x), III (y < 1 and y < x) and IV (y < 1 and y > x).
For each of the quenches performed we investigate if the system reaches a stationary state by checking whether: • The one-time quantities z(t) and e f tot = limt→+∞ etot(t) approach constants that we measure and compare to the ones predicted in Sec. 4.4.
• The kinetic and potential energy average over time to constant e kin and epot that we also compare to the ones predicted in Sec. 4.4.
In all cases an asymptotic stationary regime is reached but it is not characterised by a single temperature. Consequently, we study the spectrum of mode temperatures as defined from the (time-averaged) kinetic and potential mode energies and we compare it to the global fluctuation dissipation ratio in the frequency domain in each type of quench. This is motivated by the suggestion in [69,70] (for quenches of isolated quantum integrable systems) that these should be related.

Sector I: a paramagnetic initial condition
In this Subsection we summarise the dynamics of a system prepared in equilibrium in the paramagnetic state y > 1 and quenched to a value of J such that y > x. This is what we called Sector I in the phase diagram. Two cases need to be further distinguished within this Sector. For x < 1 energy is injected in the quench and this problem is treated in Figs. 12 and 13. For y > x > 1, instead, a small amount of energy is extracted from the system and we explain the difference that this implies in Fig. 14 and the discussion around it.
The first observation is that we confirm that, in both cases x > 1 and x < 1,R(ω = 0) = 1/T and z f = T +J 2 /T . The latter claim can be verified in Fig. 12 (d) for x > 1 and N → ∞, for example. Next we check that the dynamics approach a stationary asymptotic state in which the global one-time quantities are constant, as seen for example in panels (d) and (e) in the same figure, and the two-time correlation and linear response depend on the time delay only, see panels (a) and (c). The last question concerns the fluctuation dissipation relation between linear response and correlation function. We used the parametric plot (b) to demonstrate that the evolution does not approach a state characterised by a single temperature but, instead, χ(C) is curved and not even single valued. Having said this, the comparison of the time-dependent R(t1, t2) and −∂t 1 C(t1, t2) in panel (c) could have fooled us into the belief that the FDT holds with a single temperature. Indeed, the difference between the two functions is not visible in this scale. The Fourier analysis in (f) demonstrates that the frequency dependence of the real and imaginary parts of the linear response are highly non-trivial, though respecting the limits in the frequency interval where it does not vanish derived in Sec. 4.4.2. The fluctuation-dissipation ratio is shown in (g) and we will come back to its analysis below, where we present the mode dependent results for finite N .
We next study the evolution under the same parameters in a single system with N = 1024 modes. In the first four panels in Fig. 13 we display the time dependence of the kinetic and potential energies, e kin µ (t)  In panel (g), the ratio −ImR(ω)/(ωĈ(ω)) together with 1/T f 0.92 that is off the data but not very far away from the constant at the limit ω → ω + c (recall that the downward trend close ω + c is due to numerical inaccuracy). We note that 1/T kin 0.87 is still farther away from the high frequency limit. Results for the same parameters and finite N are shown in Fig. 13. and pot µ (t). These functions oscillate as a function of time with amplitude that slowly decreases in the cases µ = 1, µ = N − 1, µ = N , while it slowly increases for µ = 2, in this time window. The total energy of each mode shown with a solid red line displays a small downward drift towards, one may expect, a constant value. The Lagrange multiplier shows a residual time variation around a value that is slightly above the N → ∞ prediction but is very close to its finite N modified value T + λ 2 max /(4T ) (shown with a solid horizontal line). We then determine the mode temperatures from the equipartition of the kinetic and potential energies averaged over a sufficiently long time window. All modes are in equilibrium with themselves in the sense that the potential, kinetic and total mode temperatures coincide except for small deviations present in the largest modes. Panel (f) shows the non-equilibrium temperature profile with the largest modes having higher temperature than the others. T f , the temperature obtained assuming a paramagnetic final state in equilibrium at a single temperature, see Eq. (A.21), is clearly different from the mode temperatures and is not the average of them either (not shown), confirming that under these quenches the system does not equilibrate to the paramagnetic state. Finally, panel (g) displays the comparison between the mode inverse temperatures and the frequency dependence of the fluctuation-dissipation inverse temperature. The agreement is very good (except at the edge of the spectrum where both numerator and denominator vanish and it is very hard to control the limiting behaviour even in the equilibrium case, see Fig. 7). The (yellow) continuous line is the approximate theoretical prediction in Eq. (183) that captures the numerical behaviour rather accurately. We recall that it was derived assuming that the Lagrange multiplier takes its asymptotic constant value immediately after the quench, at time t = 0 + , and that the ensuing dynamics is the one of independent harmonic oscillators with mode-dependent frequencies.
As already mentioned, parameters in Sector I also permit quenches with energy extraction, realised by x > 1. We repeated the analysis of the N → ∞ equations for parameters in this regime, choosing, for example, T = 1.2 J0 and J = 1.1 J0. We do not present the data here as there are no fundamental variations with respect to what we have already discussed for the energy injection case. Having said so, the mode temperatures do present an interesting difference that we discuss with the support of Fig. 14. Also in this case, the temperatures of the modes are approximately the same for the low lying modes while a mode dependence appears close to the edge of the spectrum. However, the temperatures of the largest modes are, in this case, lower than the temperatures of the lower modes, see panel (a) and its inset. We ascribe this feature to the fact that the quench extracts energy from the system. In panel (b) we confront the mode inverse temperatures to the ones extracted from the fluctuation dissipation ratio in the frequency domain and, once again, the agreement between numerical curves for N → ∞ and finite N data is very good. Moreover, the data are also in good agreement with the outcome of the assumption z(0 + ) = z f that leads to Eq. (183), shown with a yellow solid line. with ω− = 0 and ω+ 2.10 in this case.

Sector II: large energy extraction from a paramagnetic initial state
For these parameter values the Lagrange multiplier approaches z f = 2J. We confirmed this claim with the study of several cases in this Sector (see Fig. 3). Concerning the behaviour of the other global observables, energies, correlation and linear response, and the mode-dependent ones, we differentiate three cases lying in Sector II: y > √ x, y = √ x and y < √ x, all with y > 1 and y < x.
An example of the decay of the two-time global correlation function can be seen in panel (j) in Fig. 6. It is stationary and it rapidly approaches zero with progressively decaying oscillations around this value. The linear response and correlation function are not related by an FDT with a single temperature (not shown) and in this respect there are no important differences regarding what we have just shown for energy extraction processes in Sector I. For these reasons we chose not to show more data for these parameters.
A prediction from the asymptotic analysis of the Schwinger-Dyson equations, the fact that on the curve y = √ x and y ≥ 1 FDT is satisfied with T f = J, is made explicit in Fig. 18 where the parametric plot T f χ(t1 − t2) against C(t1 − t2) is constructed for four pairs of (x, y = √ x) given in the key. The dotted line is the FDT prediction with T f = J. The agreement between numerics and analytics is very good. Additional support on the fact that the dynamics after the quench occur as in equilibrium at T f is given in panel (b) in the same figure where quenched and equilibrium correlations are indistinguishable. The latter are obtained by drawing the initial conditions drawn from equilibrium at T = T f = J and running the code with the same parameter J. Coincidence of a similar quality (not shown) is found for the other three sets of parameters used in (a).
We chose to show data for T = 1.25 J0 and J = 2 J0, parameters that lie in the region y < √ x of the phase diagram where the naive analysis of the N → ∞ equations allows for ageing behaviour. The N → ∞ and finite N data are shown in Fig. 16. The Lagrange multiplier and kinetic and potential energies approach their expected asymptotic values with an oscillatory behaviour and smoothly decaying amplitude (not shown). The two-time correlation is stationary and shows no ageing. The decorrelation from the initial configuration and from a later configuration at time t2 behave very similarly. The time-delayed C(t1, t2) shows a rapid decay towards a value close to 0.1 around which it oscillates once and next decays towards zero with damped oscillations (a). The linear response function shows a similar effect in the sense of having a fast variation at short time differences and a slower one later (c). The value of R obtained from the numerical integration agrees very well with Rst from Eq. (83). The most interesting results concern the comparison between the linear response and the correlation function in the parametric plot in panel (b). The very short time delay behaviour, for C close to 1 and χ close to 0 seems to be described by the slope dictated by T f , the value of the temperature deduced from an ageing like asymptotic scenario. However, the curves deviate from this straight line for smaller C and larger χ. When the correlation reaches a value close to 0.1 corresponding to its first oscillation, the parametric plot approaches a flat form that ensures the limit χst = 1/J. This second behaviour is reminiscent of what happens in the ageing relaxation of the same model [44].
The asymptotic analysis of the N → ∞ equations allow for an ageing solution with a diverging effective temperature for correlation values below the plateau at q, in this region of the parameter space. (In the p = 3 model, for similar sets of parameters an ageing solution though with a finite T eff is realised [7].) The numerical solution of the complete set of Schwinger-Dyson equations demonstrates a behaviour with some features that are similar to the approximate asymptotic solution, but no signature of ageing. Indeed, the parametric plot could be interpreted as formed by two pieces, one in which the form is non-trivial close to C = 1 and another one that is, on average, flat, separated by the correlation value at which the first oscillation occurs. A flat piece in the parametric plot means that the integrated response, and all other functions such as the potential energy, do not get contributions from these time scales and it corresponds to an infinite [44] effective temperature [51,52]. In practice, the parametric plot is not locally flat for small values of C but, as we can see from the mode analysis in Fig. 16 (e), the temperatures of the corresponding modes do take very large values.

Sector III: initial and final condensed states
In Fig. 17 we start discussing the behaviour of the N → ∞ model in Sector III. We use T = 0.5 and J = 0.8, a quench that injects a small amount of energy from the system. Panel (a) proves that the two time correlation approaches a non-vanishing value asymptotically. The horizontal dashed line is at q = 1 − T f /J 0.44, the theoretical value derived from the assumption of equilibration à la Gibbs-Boltzmann. Further information about the decay of the correlation functions is given in (c) where we show C(t, 0) and the off-diagonal correlation with the initial configuration C2(t, 0) against time. C2 starts at qin = 0.5 and decreases monotonically. C(t, 0) quickly decays from 1 with superimposed oscillations. Both curves should reach q2 = q0 = q asymptotically and the data are compatible with this claim.
Panel (a) in Fig. 18 shows the value of the asymptotic potential energy as a function of J/J0 for the same three values of T /J0. The agreement between e f pot and T /4 + T J(4J0) − 4J, the parameter dependence found from the energy conservation and the asymptotic value of z is extremely good.
For quenches with y = T /J0 < 1 and x > xc(y) = y, the asymptotic values of C(t1, t2) and C(t1, 0) = C2(t1, 0), q and q0, respectively, do not vanish. However, it is not always easy to extract their functional dependence on x and y, especially for parameters that get away from the shallow quenches x 1. Figure 18 (b) shows q and q0 against J/J0 for three values of T /J0, together with the naive single temperature prediction q = 1 − T f /J with 2T f = 2T kin = T (1 + J/J0) drawn with black solid lines. The agreement is good close to x = 1 though deviations are clear close to the critical line xc(y) and for large energy extraction deep inside this parameter sector. Note that the prediction of equilibration à la Gibbs-Boltzmann is such that q = 0 at xc(y) and this fact is not clear at all from the data (we have extracted q from a very short plateau in the correlation, that continues to decrease possibly to zero, in these cases). Also shown in this plot is the asymptotic value of q0. One clearly finds that q0 > q for x < 1 (injection) and q0 < q for x > 1 (extraction). We will come back to the behaviour of the correlation function below.  ) and C(t1, 0) are ordered according to q > q0 for J < J0 and q < q0 for J > J0. Fig. 17 shows the parametric χ(C) curve where one sees that the t = 0 (red) data are very different from the ones for long t2. The data for long t2 suggest that FDT has established at temperature T f . Another test of FDT, now in the time domain, is shown in (d) with the comparison between the linear response and the time derivative of the correlation. The other panels show the asymptotic values of the Lagrange multiplier (e) and kinetic and potential energies (f). These yield further support to the asymptotic value of the q parameter and T f estimated analytically under the Gibbs-Boltzmann assumption, since they demonstrate perfect agreement with the asymptotic contributions to the total energy. Nevertheless, we will revisit this claim below when studying deeper quenches in the same sector.

Panel (b) in
The companion curves for finite N are in Fig. 19. First of all, panels (a)-(d) display the time dependence of the µ = 1, N/2, N − 1, N mode energies in a system with N = 1024. While the modes µ = 1, N/2 show the usual oscillatory behaviour of a harmonic oscillator, the largest modes µ = N − 1 and µ = N are clearly out of equilibrium. The Lagrange multiplier is approximately constant and equal to the largest eigenvalue, within numerical accuracy. The spectrum of mode temperatures is plotted in (f) with a zoom over the largest modes in its inset. The profile is an equilibrium one, with Tµ being independent of µ, apart from the deviations close to the edge of the spectrum. Finally, (g) shows a comparison between the inverse mode temperatures and the N → ∞ frequency dependent effective temperature extracted from the FDR. N → ∞ and finite N results coincide (except at high frequencies where the result is biased by the numerical limitations in the computation of the Fourier transform). Higher modes (low frequency) have temperatures slightly below the temperature T f while lower modes (high frequencies) have temperatures slightly above T f , that is shown with an horizontal yellow line. This fact is another warning concerning the claim of complete equilibration to a Gibbs-Boltzmann measure. Indeed, although the data for the shallow quench that we have just described might have suggested equilibration to a Gibbs-Boltzmann probability distribution, the detailed comparison of the full time-delay dependence of the correlation function after a quench and in equilibrium (no quench) at parameters J and T f that are the ones that the equilibrium measure would have, prove that such a steady state is not reached by the dynamics. This statement is proven in Fig. 20 where we display the self-correlations stemming from the two procedures for three choices of quenches: to the critical line x = y (a case that will be further studied in the next Subsubsection), the shallow quench, and a deep quench.
We note that the comparison of the linear response functions after a quench and in equilibrium yield identical results. It is the correlation function the one that deviates from Gibbs-Boltzmann equilibrium. Figure 20 (a) has already shown non-equilibrium dynamics on the critical line y = x for y < 1. We give here further evidence for this fact. In Fig. 21 we show results for T = 0.5 < T 0 c and J = 0.5, that is to say,  Fig. 5. This quench injects a relatively small amount of energy into the system, ∆e 0.375 and takes the parameters to be on the critical line y = x.

Transition between Sectors III & IV
The self correlation is shown in panel (a) of Fig. 21. Stationarity is clear for short time delays t1 −t2, for t2 > 0 and there is some remanent waiting time dependence at these short t2. This double behaviour is reminiscent of what is seen in the relaxational dynamics where a sharp separation of time-scales exists. The data suggest that q0 = limt→∞ C(t, 0) is different from q = limt 1 −t 2 →∞ limt 2 →∞ C(t1, t2), although it seems hard to determine these values from the numerics with good precision. In the event of a two step relaxation of C(t1, t2) with an approach to a plateau at q and a further decay from it to zero, the plateau should be at q = 1 − T f /J = 0.25, shown with a dashed horizontal line. The data still lie below this value and we infer that they might not converge to a plateau but simply decay to zero.
Further information about the decay of the correlation functions is given in (c) where we show C(t, 0) and the off-diagonal correlation with the initial configuration C2(t, 0) against time. C2 starts at qin = 0.5 and decreases monotonically. C(t, 0) quickly decays from 1 with superimposed oscillations. Both curves seem to join and slowly and monotonically decay to zero.
The parametric plot of the linear response function, χ(t1, t2), against the correlation, C(t1, t2), for fixed t2 and using t1 − t2 as a parameter, for three different values of the waiting time t2, is shown in panel (b). The χ(C) curve for t = 0 does not have any special form. The curves for late t2 are close to the straight line 1/T f for time-delays that correspond to the first oscillations of the correlation and linear response, they then oscillate, and for longer time-delays the parametric construction becomes flat, with χ approaching, for C → 0, the expected static susceptibility 1/J = 2. Again, this second flat regime is reminiscent of the one found for the relaxation stochastic dynamics at and below the critical temperature [86,87].
The fluctuation-dissipation theorem relating χ to C in a stationary regime with target temperature T f = 0.375 obtained from Eq. (A.19), is indicated as a straight dashed line. The presentation in panel (d) confirms the agreement between linear response and time variation of the time-delayed correlation function dictated by the FDT at temperature T f for C ≥ 0.2, say. However, it clearly breaks for smaller values of C. As for the other quenches, we checked that the response R(t2 + τ, t2) coincides with the one derived by anti-transforming the theoretical Fourier transform of the stationary asymptotic response, given by Eq. (83), Rst.
In (e) we provide a close look at the time evolution of z(t). As one can see, z(t) approaches the predicted value z f = 2J. From the time evolution of the energy (f), we observe that the total energy is constant in time, as it should be, while the kinetic and potential energies quickly approach their asymptotic values, which coincide within numerical accuracy with the ones predicted in Sec

Sector IV: large energy injection on a condensed state
In Fig. 22 we show results for T = 0.5 < T 0 c and J = 0.25. This quench injects a large amount of energy into the system, ∆e = 0.5625, which is sufficient to take it out of the initial condensed state. Had the system reached an equilibrium paramagnetic state asymptotically after the quench its temperature would be T f 0.320, so that T f /J 1.281, from Eq. (A.20). This, however, is not consistent within the N → ∞ analysis and, accordingly, it is not realised dynamically.
The self correlations shown in (a) present large oscillations with weakly decaying amplitude around the expected asymptotic value limt 1 −t 2 →∞ limt 2 →∞ C(t1, t2) = 0 for all values of t2. C(t1, t2) satisfies stationarity for sufficiently late t2, as one can see from the plot where the curves for t2 coincide, within numerical accuracy. The parametric plot of the susceptibility χ against the correlation C, shown in (b), is rather complex and very far from linear. The FDT relating χ to C in a stationary regime with putative target temperature T f 0.320 is indicated as a straight dashed line. This behaviour is confirmed by the time evolution of the stationary response function, R(t1, t2), see panel (d). R(t1, t2) does not coincide with − 1 T f ∂t 1 C(t1, t2), with T f the target temperature. In  The dashed line is −1/T f and is clearly off the data. The different curves were computed for various waiting times t2 given in the key. In Fig. 23 we compare the FDR (averaged over different waiting times to get rid of the undesired oscillations) to the mode temperatures of the finite N system. The correspondence between the two ways of extracting the temperatures is good. The yellow line is the approximate prediction in Eq. (183) stemming from the independent harmonic oscillator approaximation that for this kind of quench is quite far from the numerical results though it captures the qualitative features.  Figure 23: Sector IV. Large energy injection on a condensed state. Same parameters as in Fig. 22. Comparison between the mode inverse temperature for a finite size system (red), the inverse temperature for the fluctuation dissipation ratio in the infinite size limit (black curve) and the analytical expression for the mode inverse temperatures for independent harmonic oscillators (yellow curve). All measured in units of J0. The frequency interval in which ImR(ω) is non-zero is (ω−, ω+) = (1/(2 √ 2), 3/(2 √ 2)) (0.35, 1.06).

Integrals of motion
In Sec. 3 we recalled the relation between the p = 2 spherical disordered model and the classical integrable model introduced by Neumann. In this Section we will present some results concerning the behaviour of the integrals of motion. A key issue we address here is how these influence the statistical properties in the steady state.

The integrals of motion landscape
In Sec. 2.2.1 we studied the potential energy landscape of the p = 2 disordered model and we found that it has N extrema corresponding to s = ± √ N vµ with vµ the N eigenvectors of the interaction matrix Jij. These directions turn out to be the extrema of the integrals of motion landscape as well [88]. This claim is proven easily. Indeed, where we labeled the constants of motion defined in Eq. (48) with µ, the eigenvalue index, implies the following two conditions It is clear that the first relation makes the second one hold identically. Moreover, replaced in the definition of the Iµs one finds that has to be extremised under the global spherical constraint on the sµ. This is just the analysis of the potential energy landscape that we performed in Sec. 2.2.1, leading to s * = ± √ N vµ and z * = λµ for all the saddles in the landscape.
In the cases in which y < 1 the initial state is condensed and the integrals of motion of the modes µ = N and µ = N scale very differently with N . For the modes in the bulk This expression involves two sums that will appear again and again in the rest of this Section, and depend only on the pre-quench parameters J0 and T , For µ = N the second sum reads The superscript n means that the eigenvalues have been rescaled by J0 in such a way that they vary in the interval (−2, 2) and S (1n) µ is just a number. Using these definitions we find This expression has a rather simple dependence on J that only appears as a prefactor in front of the sum of three terms within the square brackets. (This fact justifies the similarity of the bulk numerical values of Iµ in, for example, panels (c) and (d) in Fig. 26.) In the large N limit z(0 − ) → λ (0) N and For the largest mode µ = N , if y < 1, we obtain In the limit N → ∞ we can use S N → −1/J0 and the known form of qin to find The parameter dependence of the slope in the right-hand-side seen as a function of N is verified numerically in Fig. 24 (a). The good match between this form and the numerics indicates that S N must be negligible with respect to the first term that is O(N ); indeed, we have checked numerically that S (2) N is O(1). y < 1: paramagnetic initial states.
If, instead, y > 1, there is no condensed mode and the integrals of motion are and all of order 1. The particular case µ = N is displayed in the panel (b) of Fig. 24, showing no dependence on N for N > ∼ 2000, as expected.

Summary.
Summarising, Fig. 24 displays the integral of motion associated to the mode at the edge of the spectrum for y < 1 (a) and y > 1 (b). All data points were obtained using a single realisation of the random matrix. For the condensed initial conditions we clearly see the linear scaling with N . For the non-condensed ones the variations show a weak N dependence for small N plus a possible variation due to the fluctuations in the realisation of the eigenvalues. The result is distinctly finite in this case.
As for the time-evolution of these averaged quantities, we have verified (not shown) that each of them are conserved Iµ(0 + ) = Iµ(t) for all µ, and that they satisfy the two constraints µ Iµ(0 + ) = µ Iµ(t) = N and The N th mode The N th mode has a different scaling with system size. We have already computed IN (0 + ) in Eq. (209) and we want to compare it to what it should read in equilibrium at T f : The difference between the two, The first term diverges linearly with N . We have already argued that the second one is O(1), see the discussion after Eq. (209). Therefore, we focus on the slope of the difference, seen as a function of N . After replacing q f , qin and T f by their expressions in terms of J, J0 and T in the large N limit, and S (1) This form is validated by the numerical data in Fig. 25 for three pairs of T and J, all in Sector III. Quite clearly, the differences between the actual Iµ and the ones in equilibrium at a temperature T f are proportional to J0/J − 1, a factor that vanishes for J0 = J. Otherwise, the differences are finite for µ = N and proportional to N for the last mode, making the mode-by-mode difference non-vanishing for J = J0. This fact confirms, then, the lack of equilibration to a Gibbs-Boltzmann measure with a single T f .
On the curve y = √ x the asymptotic analysis for N → ∞ predicts thermalization at temperature T f = J. The numerical analysis of the Schwinger-Dyson equations confirms this prediction as the correlation and linear response are linked by FDT and the time-dependence of the correlation function after an instantaneous quench is identical (within numerical accuracy) to the one found in equilibrium at this temperature. We shall briefly analyse the behaviour of the constants of motion in this case.
On the special line y = √ x, T 2 = JJ0 and Eq. (210) yields the following expression for the averaged integrals of motion where λµ = λ (0) µ /J0 are the normalised eigenvalues that vary in the interval (−2, 2). From Eq. (211) with q f = 0, z f = 2J valid in the N → ∞ limit, and using Eq. (205) the constants of motion in an equilibrium state at temperature T f are We observe that the two sets of values are very similar if J is close to J0. Numerically, we find that, in the special case T = 1.25 J0 and J = 1.5625 J0, the ∆Iµ is of order of 10 −2 . We conclude that even in this case, in which the asymptotic analysis predicts that the global, mode-averaged quantities, behave as in equilibrium, a prediction that seems to be confirmed by the numerics, the constants of motion are not exactly the same in the initial state and in the thermal state the system would reach in case of thermalisation.
In order to properly interpret these results, it is important to keep in mind that, strictly speaking, the conserved energy dynamics of an isolated (finite size) system should keep memory of the initial conditions, even if the system is non-integrable. In our problem, we see this information encoded in the Iµs. More so, not even in the N → ∞ limit this memory is erased as the ∆Iµs remain finite.

Independent harmonic oscillators
We have seen in Sec. 5.7 that the N → ∞ system decouples into independent harmonic oscillators in the asymptotic long time limit (taken after N → ∞) since z(t) → z f . A natural idea is to check whether we can identify the integrals of motion Iµ(0 + ) with the ones that an ensemble of harmonic oscillators with spring constants mω 2 µ = z f − λµ in equilibrium at the temperatures Tµ would have. On the one hand, we note that s 2 µ (t) and p 2 µ (t) are not constant for harmonic oscillators but their time averages are, so we evaluate Iµ finding, in this framework, In the N → ∞ limit the value of z f is expected to be T + J 2 /T for x < y and 2J for y < x. Imposing the identity between the Iµ osc and the Iµ(0 + ) given in Eqs. (206) and (208) should yield a better estimate of the temperatures Tµ than the one explained in Sec. 5.7. We leave this analysis aside for the moment.
On the other hand, once the oscillators have decoupled and the Lagrange multiplier stabilised, the mode total energies, or the time-averaged kinetic and potential energies are also constants of motion. In the next Subsection we investigate to what extent the Iµ are proportional to e tot µ = 2e kin µ . Figure 26 shows the spectrum of integrals of motion J Iµ (red data points) together with the mode kinetic temperatures T kin µ = 2e kin µ (that, we have checked, are in agreement within numerical accuracy with the potential ones T pot µ away from the edge of the spectrum), for parameters in the four sectors of the dynamic phase diagram. We show the data using the same vertical scale in all cases. Although the data for J Iµ are noisier than the ones for T kin µ , the two are in fairly good agreement in the bulk of the spectrum, far from the edge, in all panels. In the same figures we include the values of the global kinetic temperature, T kin defined in Eq. (121) and we see that the mode temperatures are very close to it again away from the edge of the spectrum and satisfying the constraint N −1 µ T kin µ = T kin . The parameters in Sector II, panel (b), are on the curve y = √ x on which the data for the global correlation and linear response show equilibrium at a single temperature T f = J. The blue data points for the mode kinetic energies are precisely on this value. The integrals of motion scatter, however, quite a lot around this value. Importantly enough, the difference ∆Iµ is very small in this case.

The integrals of motion and the mode temperatures
In sector III, panel (c), the data for J Iµ tend to be relatively flat in the bulk of the spectrum and very close to T kin = T f = T (1 + J/J0)/2. This result can be derived from Eq. (206) taking advantage of a simple rearrangement of the two sums, S  183) are all identical to T f = T kin in this Sector, consistently with the numerical data away from the edge. The fact that the system is not in proper Gibbs-Boltzmann equilibrium is due to the fact that the higher lying integrals of motion do not comply with this temperature.
The data for J Iµ in Sector III look very similar to the ones found in Sector IV, see panel (d). Indeed, the data almost coincide far from the edge, since they are both determined by Eq. (206) that has a weak dependence on J (the only parameter that takes a different value in panels (c) and (d)). Close to the edge, the J Iµ differ since in III (c) there is scaling with N while in IV (d) there is not. The mode kinetic energies are mode independent and identical to T (1 + J/J0)/2 contrary to what the approximation in the last line of Eq. (183) tells. This means that for these quenches the assumption in Eq. (180) fails.
We reckon that the kinetic temperature T kin = 2e kin should be equal to the sum of the mode kinetic temperatures T kin = 1/N N µ=1 T kin µ . We have checked numerically this property. We leave a more detailed analysis of the comparison between the two and their use to build a GGE for a future publication.

Fluctuations of the integrals of motion in the equal energy hypersurface
The fact that the integrable system reaches a state that is very close to thermal equilibrium in sector III, the sector with the lowest total energy in the phase diagram, allows us to infer some properties of the phase space structure of the model. Only a system for which it is possible to visit all configurations with the same energy in the course of its dynamical evolution is capable of reaching thermal equilibrium. An integrable model cannot achieve this goal since it is bound to wander in a region of the phase space compatible with the values that the integrals of motion (IOM) take on the initial configuration. The dynamics is constrained inside the phase space region composed by configurations which have the same values of the Iµs ∀µ. Such regions can be labeled with the values of the Iµs (which also define the energy of the group of configurations since − λµIµ = 2H), and we shall call them iso-IOMs-regions in the following.
A close-to-thermalised dynamics in an integrable system should be indicative of a substantial overlap between the constant energy manifold and the equal IOMs region in phase space. This claim is, however, highly nontrivial since the equal energy manifold is 2N − 1 dimensional while the iso-IOMs-region has only 2N − N = N dimensions. In the large N limit, the equal energy manifold is huge with respect to the equal IOMs one (for N = 2 the constant energy region is a volume and the equal IOMs-configuration is a surface).
Our hypothesis for the dynamics of the Neumann model with parameters close to x = 1 in phase III, is that the constant energy manifolds have a substantial overlap with any iso-IOMs-region that include a configuration with the given energy. In order to test this guess, we studied the following quantity: where the average is a microcanonical one given by The quantity σ 2 µ (e) measures how large are the fluctuations in the value of a given IOM Iµ in the set of configurations with the same energy e. According to the discussion at the beginning of this Section, if, for a given energy, we observe a small value in σµ(e), this indicates a tendency of the integrable system to thermalise. In order to perform the averages over equal energy configurations, we replace the microcanonical average by a canonical one, introducing a Lagrange multiplier β and a measure exp(−βH)/Z, fixing the average energy density of the ensemble. For large N , the fluctuations of the energy average are small and we get a good approximation to the microcanonical mean. The advantage of the canonical measure is that the Iµ s are expressed in terms of canonical averages of the same kind as those which we were using in the previous Sections to describe the initial state of the dynamics. Moreover, once z is fixed, the Hamiltonian is quadratic, and this allows one to express the higher order average I 2 µ , which includes products of 4, 6 and 8 phase space variables, in terms of quadratic averages s 2 i and p 2 i . A straightforward numerical calculation of σ 2 µ (e) is then possible. We show numerical results in Fig. 27. In panel (a) we plot σ 2 µ , the average of σ 2 µ (e) in the first three fourths of the spectrum We can clearly observe that it is very small for low energies and that it increases by several orders of magnitude as we increase the energy density of the system. Such is the behaviour of σ 2 µ (e) for the great majority of the modes. In panels (b) and (c) we show the behaviour of σ 2 µ (e) for µ close to N , the edge of the spectrum. We observe that σ 2 µ (e) exhibits a maximum at low energies. Finally, in panel (d), σ 2 µ=N (e) is plotted. At very low energies, σ 2 N (e) is very large and scales with N . As we increase energy, σ 2 µ (e) decreases abruptly until it reaches a value that is inversely proportional to N .
Summarising, we observe that, far from the edge of the spectrum, the fluctuations in the IOMs are very small for sufficiently low energies. This means that the low energy configurations of the model have very similar values of the Iµs, at least for µ far from the edge of the spectrum. For modes close or at the edge of the spectrum fluctuations can be very large. These results suggest that the iso-IOMs-regions which lie at low energies have very similar values for the Iµs for µ far from N , and only differ in the values of the Iµs close to the edge of the spectrum. In fact, in sector III, the sector with the lowest energies in the phase diagram, we have verified in previous Section that the initial state and the final thermal state at temperature T f which partially describes the long-time dynamics have very similar values of the IOMs far from the edge of the spectrum, with discrepancies only near the edge of the spectrum.

Conclusions
This paper continues our study of the conserved energy dynamics following sudden quenches in classical disordered isolated models ruled by Newton dynamics.
In the p = 3 strongly interacting case [7] all quenches reach an asymptotic regime in which a single (or double) temperature dynamical regime establishes. The systems either equilibrate to a paramagnetic state with a proper temperature, they remain confined in a metastable state with restricted Gibbs-Boltzmann equilibrium at a single temperature, or age indefinitely after a quench to the threshold with the dynamics being characterised by one temperature at short time delays and another one at long time delays, similarly to what happens in the relaxational case [20,72,71,73,9]. The two temperatures T f and T eff depend on the pre and post quench parameters in ways that were determined in [7]. For the sake of comparison, the dynamic phase diagram of the p = 3 isolated model is reproduced in Fig. 28 (a).  In the p = 2 model the conserved energy dynamics are quite different both from the relaxational ones [42,43,44,26,45] and the isolated p = 3 interacting case [7]. This is due to its integrability, made explicit by its relation to the Neumann integrable model. The main results in this paper are the following.
• We identified the dynamic phase diagram according to the asymptotic behaviour of the static susceptibility, the Lagrange multiplier imposing the spherical constraint (an action density), and the long time-delay limit of the two-time correlation function. The phase diagram is shown in Fig. 28 (b), and it can be compared to the one of the p = 3 isolated case shown on its left side.
• In the analysis of the Schwinger-Dyson equations we distinguished four sectors in the phase diagram depending on the initial state (being condensed or not) and the final value of the static susceptibility. We reduced these four sectors to three phases, indicated with different colours in the Figure, in which the dynamic behaviour is different. Basically, they are distinguished by two "order parameters", χst and q, the static susceptibility and the asymptotic value of the two-time correlation.
• In none of the phases the system equilibrates to a Gibbs-Boltzmann measure. Accordingly, there is no single temperature characterising the values taken by different observables in the long time limits, not even after being averaged over long time intervals.
• There is one case, quenches from the condensed equilibrium state in which energy is extracted or injected in small amounts, in which the dynamics of the global observables, those averaged over all modes, are close to the ones at a single temperature T f . However, a closer look into the mode dynamics and its observables exhibits the fact that the modes are not in Gibbs-Boltzmann equilibrium. Moreover, for deep quenches in the same sector III one clearly sees that standard thermal equilibration is not reached.
• Another special case is provided by quenches with T > J0 and T 2 = JJ0. On this special curve the global observables satisfy thermal equilibrium properties at T f = J.
Much has been learnt from the evolution of the system with finite number of degrees of freedom using a formalism that allows one to show that in the infinite size and long time limit (taken after the former one) the modes decouple and become independent harmonic oscillators. Once this regime is reached, the mode energies can be associated to mode temperatures, via standard equipartition, and a temperature spectrum obtained. A naive approximation to determine their dependence on the control parameters was explained and leads to a variation from sector to sector of the phase diagram. We compared these forms with the numerical measurements and the agreement is quite good in all cases.
Having obtained the total energies of the modes, and from them the mode temperatures, we put to the test the recently proposed relation between them and the frequency dependent effective temperature stemming from the fluctuation-dissipation relation of the spin observable [69,70]. We found excellent agreement between the two in all phases of the phase diagram.
The p = 2 spherical system turns out to be equivalent to the Neumann classical mechanics integrable model. We stress the fact that in the field of classical integrable systems, the model of Neumann was usually defined and studied having only a few degrees of freedom. Here, as we are interested in searching for a statistical description of the post-quench dynamics, we dealt with the limit of large, and even diverging, number of degrees of freedom.
The N − 1 integrals motion of the Neumann model have been identified by K. Uhlenbeck [16,18]. After a trivial extension that allows us to deal with the large N limit, we studied their scaling properties with system size. In cases in which the initial state is condensed, the integrals of motion associated to the edge of the spectrum also scale with N . The distance between their values and the ones they would have taken in equilibrium at a single temperature T f gave us a rough measure of distance from Gibbs-Boltzmann equilibrium. Importantly enough, in the particular case in which the global correlation and linear response behave as in thermal equilibrium at T f = J, that is to say, parameters on the curve T 2 = JJ0 in sector or phase II, the integrals of motion are not identical to the ones expected in equilibrium. This proves that not even in this case the system is able to fully equilibrate. The N − 1 integrals of motion could be used to build a putative Generalised Gibbs Ensemble, or they may be a guideline to choose the ones with good scaling properties. We will investigate this problem in a sequel to this publication.
A Asymptotic analysis in the N → ∞ limit In this Appendix we give additional details on the study of the full set of equations (75)-(78) that couple the correlation C and linear response R functions derived in the N → ∞ limit. Based on a number of hypotheses that we carefully list below, we analyse the behaviour of the model in the long times limit.

A.1 Stationary dynamics
Consider the system in equilibrium at T with parameters J0, m0 and let it evolve in isolation with parameters J, m. We will assume that the dynamics approach a steady state in which one-time quantities approach a constant. This assumption does not apply to certain quenches of the isolated system. Still, we investigate the consequences of the stationary assumption.

A.1.1 The asymptotic values
Let us assume that the limiting value of the Lagrange multiplier is a constant The linear response was analysed in the main body of the paper and, independently of the quench parameters it is given byR in the frequency domain, where the Fourier transform has been computed with respect to the time difference t1 − t2. This result only assumes a long-time limit in which z f is time-independent.
A.1.2 The parameters q 0 and q 2 We could estimate the asymptotic value of z(t1) taking the long t1 limit of Eq. (78); however, without the use of FDT we cannot compute the integral involved.
We are tempted to propose q0 = q2 . (A.4) The interpretation of this result in terms of the evolution of real replicas is that the asymptotic value of the self-correlation between times t1 and 0, q0, is the same as the asymptotic value of the correlation between two replicas evaluated at the same times t1 and 0 with t1 diverging.

A.1.3 The two-time correlation function
Allowing for the two-time correlation function not to decay to zero but to a finite value q, we separate this contribution explicitly and we write C(t1, t2) = q + Cst(t1 − t2) with limt 1 −t 2 →∞ Cst(t1 − t2) = 0. We also use limt 1 →∞ C2(t1, 0) = q2 leaving the possibility of there being a different asymptotic value for this quantity open.
The dynamic equation now becomes an equation onC that reads Using the causal properties of the linear response, one can extend the lower limit of the last integral to −∞ and safely take the upper limit to infinity since we are interested in the long time limit t2 → ∞ while keeping t1 − t2 fixed. The Fourier transform of this term with respect to t1 − t2 yields R(−ω)Cst(ω) that using R(−ω) = R * (ω) simply becomes R * (ω)Cst(ω). Proceeding in a similar way with the first integral one finds that it equalŝ R(ω)Cst(ω). Thus, The factor between parenthesis in the first term on the right-hand-side can be written as If we now assume q0 = q2, this equation imposes if q = 0 then z f = 2J and lim In the no quench case we can check whether C2(t1, 0) = q2 = qin for all t1 is a solution of the corresponding equation. Plugging this form in the evolution equation for C2 we find that, at t1 = 0, either qin = 0 (the paramagnetic case) or z f (0) = 2J 2 /T (1 − qin) and after replacing qin by its expression as a function of T /J the correct equilibrium z f = 2J is recovered. However, for t1 = δ and later times, C2(t1, 0) cannot remain constant for initial conditions with qin = 0, due to the non-trivial contribution of the term J 2 t 1 0 dt R(t1, t )qin = J 2 /m qin δ. If we assume that C2(t1, 0) approaches a constant q2 possibly different from qin, the equation for C2(t1, 0) in the long t1 limit reads • From paramagnetic to ageing We have chosen to single out in these expression the parameters T /J0 and (mJ0)/(m0J) that we have used in [7] to characterise the dynamic phase diagram of the p ≥ 3 model. In the two cases (condensed to condensed and paramagnetic to paramagnetic) in which the no quench limit can be taken we naturally recover T f = T . In the ageing case T f should be the temperature of the stationary regime. As we will show from the numerical solution to the full dynamic equations, and the mode by mode analysis, not all these cases are realised. Below we prove analytically that a state with a single T f characterising the fluctuation-dissipation relation can be realised for particular choices of the parameters only. Moreover, these conditions are not restrictive enough, as the numerical solution of the full equations show that FDT is not realised for quenches allowed by them.

A.3.2 Limits of validity
Let us consider the cases y ≥ 1 and y ≤ 1 separately.
The case y ≥ 1: Quench from a paramagnetic equilibrium state The double condition 0 ≤ T f ≤ J translates into This upper bound is the dashed line representing f (x) in the phase diagram for y > 1. Note that this curve has no finite limit.
The case y ≤ 1: Quench from a condensed equilibrium state The lower bound is trivially satisfied while the upper bound implies The dashed line for y < 1 in the phase diagram represents g(x).
In the condensed case the energy conservation implies We note that q vanishes on the curve y = 2x/(1 + x), that it equals one for y = 0 and that, for fixed x > 1 it increases from (x − 1)/(2x) at y = 1 to 1 at y = 0. Moreover, q is different from zero on the lines y = x and y = 1.

A.3.3 Consistency with FDT at T f
We now reason as follows. We know, from the numerical analysis of the Schwinger-Dyson set of equations, that χst = limt→∞ t 0 dt R(t − t ) = 1/T for x < y and χst = t 0 R(t) = 1/J for x > y. If we evaluate the integral of the linear response using FDT at a single temperature we obtain Using these results we can check whether there is a set of x, y for which T f derived in the previous Section from the conservation of the energy is consistent with these conditions. We list below the conclusions drawn in the four parameter Sectors of the phase diagram.
y > 1 and x < y (Sector I) T f can equal T only for x = 1 that implies no-quench.
y > 1 and x > y (Sector II) Only on the curve y = √ x and y > 1 FDT holds at T f . The validity of FDT at T f for parameters lying on the curve y = √ x is verified numerically in Fig. 18. For all other values of x, y in this Sector FDT cannot be satisfied.
y < 1 and x > y (Sector III) In this Sector we do not find any contradiction. This reasoning suggests that the dynamics may satisfy FDT in this Sector. The no-quench case x = 1 for y < 1 is obviously included. y < 1 and x < y (Sector IV) There is no solution and FDT at a single temperature is excluded.

B The harmonic oscillator
Take, as an example, a system constituted of a single point-like particle with mass m, under a harmonic potential V (x). The dynamics of this problem is given by the familiar equation

B.1 Equilibrium initial conditions
Let us take initial conditions in canonical equilibrium within a harmonic potential V0(x). The probability distribution of the initial conditions is with β = 1/(kBT ) the inverse temperature, using the same notation adopted in the main part of the paper for the initial temperature T . The averaged kinetic energy of the ensemble of initial states sampled with this probability distribution function is 1 2m the equipartition of the kinetic energy. We recall that angular brackets indicate average over the initial conditions sampled as above. The averaged total energy is

B.2 Potential energy quench
Make now a quench in the potential that corresponds to V0 → V and do it so quickly that the phase space variables do not change and remain p0, x0. By performing this abrupt change one injects or extracts a finite amount of energy, ∆E = H(x0, p0) − H0(x0, p0) = V (x0) − V0(x0) . (B.5) The energy surface on which the dynamics will take place is the one of the post-quench energy E(0 + ) = p 2 0 /(2m) + V (x0). We now focalise on V being a harmonic potential. The Newton evolution of each initial configuration is x(t) = x0 cos ωt + p0 mω sin ωt , p(t) = −mωx0 sin ωt + p0 cos ωt .
Although the measure P is still Gaussian it does not have the same covariance as the initial P0.
The averages and the variances of the position and momentum can be computed directly from the solutions to the equations of motion. The averages vanish and for the variances one finds σ 2 x (t) = x 2 (t) = x 2 0 cos 2 ωt + p 2 0 1 m 2 ω 2 sin 2 ωt , σ 2 p (t) = p 2 (t) = x 2 0 m 2 ω 2 sin 2 ωt + p 2 0 cos 2 ωt . (B.10) One readily verifies that, as expected, the averaged total energy is conserved for t > 0 (B.11) since each trajectory does conserve its initial energy. The averaged total energy is, however, different from the one right before the quench, T = etot(t = 0 − ) = etot(t = 0 + ) = T 1 + ω 2 /ω 2 0 . Time-independent values of the variances are found from the average over a long time window:

B.3 Fluctuation dissipation relations
The linear response of the harmonic oscillator defined in (B.1) to an infinitesimal perturbation modifying its energy as H → H − hx, and therefore adding a linear term in h to Eq. (B.1) instantaneously at time t2, is R(t1, t2) = δx(t1) δh(t2) h=0 = sin ω(t1 − t2) mω θ(t1 − t2) (B.14) while the product of the unperturbed position at the times t1 and t2, after a long-time average over the reference time t2, is where etot is the total energy of the harmonic oscillator.

B.4 The Generalized Gibbs Ensemble
The harmonic oscillator dynamics conserves its total energy etot and the GGE density function is