Quenched dynamics of classical isolated systems: the spherical spin model with two-body random interactions or the Neumann integrable model
Leticia F. Cugliandolo, Gustavo S. Lozano, Nicolas Nessi, Marco Picco, Alessandro Tartaglia
QQuenched dynamics of classical isolated systems:the spherical spin model with two-body random interactions orthe Neumann integrable model
Leticia F. Cugliandolo , Gustavo S. Lozano , Nicolás Nessi , Marco Picco and Alessandro Tartaglia Laboratoire de Physique Théorique et Hautes Energies, UMR 7589,Sorbonne Universités et CNRS, 4 place Jussieu, 75252 Paris Cedex 05, France Departamento de Física,FCEYN Universidad de Buenos Aires & IFIBA CONICET,Pabellón 1 Ciudad Universitaria, 1428 Buenos Aires, ArgentinaNovember 12, 2018
Abstract
We study the Hamiltonian dynamics of the spherical spin model with fully-connected two-body interactionsdrawn from a zero-mean Gaussian probability distribution. In the statistical physics framework, the potentialenergy is of the so-called p = 2 spherical disordered kind, closely linked to the O ( N ) scalar field theory. Mostimportantly for our setting, the energy conserving dynamics are equivalent to the ones of the Neumann integrablesystem. We take initial conditions from the Boltzmann equilibrium measure at a temperature that can beabove or below the static phase transition, typical of a disordered (paramagnetic) or of an ordered (disguisedferromagnetic) equilibrium phase. We subsequently evolve the configurations with Newton dynamics dictatedby a different Hamiltonian, obtained from an instantaneous global rescaling of the elements in the interactionrandom matrix. In the limit of infinitely many degrees of freedom, N → ∞ , we identify three dynamical phasesdepending on the parameters that characterise the initial state and the final Hamiltonian. We obtain the global dynamical observables (energy density, self correlation function, linear response function, static susceptibility,etc.) 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 notat) thermal equilibrium à la Gibbs-Boltzmann. Surprisingly enough, in the N → ∞ limit and for a particularrelation between parameters the global observables comply Gibbs-Boltzmann equilibrium. We next set theanalysis of the system with finite number of degrees of freedom in terms of N non-linearly coupled modes. Theseare the projections of the vector spin (or particle’s position on the sphere) on the eigenvectors of the interactionmatrix, the most relevant being those linked to the eigenvalues at the edge of the spectrum. We argue that in asystem with infinite size the modes decouple at long times. We evaluate the mode temperatures and we relatethem to the frequency-dependent effective temperature measured with the fluctuation-dissipation relation in thefrequency domain, similarly to what was recently proposed for quantum integrable cases. Finally, we analysethe N − integrals of motion, notably, their scaling with N , and we use them to show that the system is outof equilibrium in all phases, even for parameters that show an apparent Gibbs-Boltzmann behaviour of globalobservables. We elaborate on the role played by these constants of motion in the post-quench dynamics and webriefly discuss the possible description of the asymptotic dynamics in terms of a Generalised Gibbs Ensemble. a r X i v : . [ c ond - m a t . s t a t - m ec h ] D ec ontents p = 2 spin model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2 The statics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.2.1 The potential energy landscape . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2.2 The free-energy density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.3 Relaxation dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.4 Relation with the O ( N ) λφ model in the large N limit . . . . . . . . . . . . . . . . . . . . . . . 11 Ansatz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 244.4.8 The correlation function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 244.4.9 Numerical results preview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 N . . . . . . . . . . . . . . . . . . . . . . . . . 285.6 Energy variation at the quench . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 315.6.1 Pre-quench energies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 315.6.2 Post-quench energies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 315.7 Independent harmonic oscillators in the asymptotic limit . . . . . . . . . . . . . . . . . . . . . . . 325.8 Exact solution of the mode dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 Numerical results 34 N → ∞ limit 67 A.1 Stationary dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67A.1.1 The asymptotic values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67A.1.2 The parameters q and q . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67A.1.3 The two-time correlation function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67A.1.4 The correlation with the initial condition . . . . . . . . . . . . . . . . . . . . . . . . . . . 68A.2 The off-diagonal correlations in the no quench problem . . . . . . . . . . . . . . . . . . . . . . . . 69A.3 Energy conservation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69A.3.1 Final temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69A.3.2 Limits of validity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70A.3.3 Consistency with FDT at T f . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 B The harmonic oscillator 71
B.1 Equilibrium initial conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71B.2 Potential energy quench . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72B.3 Fluctuation dissipation relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73B.4 The Generalized Gibbs Ensemble . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
C Discrete time version 74 Introduction
In the past decade, atomic physics experiments have been able to test the global coherent quantum dynamicsof interacting systems. This achievement has boosted research on the dynamics and possible equilibration ofmany-body isolated systems [1]. Some of the quantum instances realised in the laboratory are low dimensionaland considered to be integrable. Therefore they are not able to act as a bath on themselves and questions on howto describe their asymptotic dynamics pose naturally. With the aim of describing their asymptotic states, theGeneralised Gibbs Ensemble (GGE), an extension of the canonical Gibbs-Boltzmann density operator that aimsto 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 isolatedinteracting 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 modelsbelong to the class of p spin spherical disordered models with, however, properties that render their constantenergy dynamics very different, as we will show here.The spherical fully-connected p -spin disordered models are paradigms in the mean-field description of glassyphysics. They are solvable models (in the thermodynamic limit) that successfully mimic the physics of (fragile)glasses for p ≥ and domain growth for p = 2 . The connection between the p = 2 model, in its classical andquantum formulations, with coarsening phenomena is made via its relation to the celebrated O ( N ) λφ modelin 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 sphericalsystems is vast. Numerous aspects of their behaviour are very well characterised, even analytically (see, e.g. , thereview articles [9, 10, 11]).In Ref. [7] we studied the Hamiltonian dynamics of the p = 3 spherical disordered spin model. By adding akinetic 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 − dimensional sphere under theeffect of a complex quenched random potential [12, 13, 14]. Here we will focus on the Newtonian dynamics ofthe 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 ofclassical mechanics [15], the constrained motion of a classical particle on S N − under a harmonic potential, fora special choice of the spring constants. The only difference is that in the p = 2 model one imposes the sphericalconstraint 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 − conserved charges in involution [16, 17, 18]. With these at hand, wewill be able to study the statistical properties in depth and construct a candidate GGE to describe the long-timedynamics.We perform instantaneous quenches towards a post-quench disordered potential that keeps memory of thepre-quench one, mimicking in this way the “quantum quench” procedure in a classical setting. The changein the potential energy landscape induces finite injection or extraction of energy density in the sample. Thesubsequent dynamics conserve the total energy. We sample the initial conditions from canonical equilibrium ata tuned temperature, choosing in this way initial configurations typical of a paramagnetic equilibrium state athigh temperature or a condensed, ferromagnetic-like, state at low temperature. The control parameters in thedynamic phase diagram that we will establish are the amount of energy injected or extracted and the initialtemperature 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 inwhich the infinite size system remains confined (condensed, in the statistical physics language) and cases in whichit does not. In none of them a Gibbs-Boltzmann equilibrium measure is reached, contrary to what happens inthe strongly interacting p = 3 case. The role played by the N − integrals of motion on the lack of equilibrationof 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 obtainedthem can go directly to Sec. 8 where we sum up our findings and we present a thorough comparison betweenthe 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 odel (static, metastable and relaxational dynamic properties) studied from the statistical mechanics point ofview. In the following Section we explain the relation between the disordered spin system and the integrableNeumann model of classical mechanics. We also explain in this Section the statistical description of the long-term dynamics of integrable systems provided by the Generalised Gibbs Ensemble proposal. In Sec. 4 we presentour analytic results for the dynamics of the model in the N → ∞ limit and in Sec. 5 we go further and we setthe 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 − integrals of motion inthe 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 ofthe paper. 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 simplestpossible magnetic model with quenched random interactions. This model, as we will explain below, shares manypoints in common with the O ( N ) model of ferromagnetism when treated in the infinite N limit. Its staticproperties have been derived with a direct calculation and using the replica trick. Its relaxational dynamics arealso analytically solvable. The reader familiar with this model can jump over this Section and go directly to thenext one where the relation with Neumann’s model and integrability are discussed. p = 2 spin model The p = 2 spin model is a system with two-spin interactions mediated by quenched random couplings J ij .The potential energy is H pot [ { s i } ] = − N (cid:88) i (cid:54) = j J ij s i s j . (1)The coupling exchanges are independent identically distributed random variables taken from a Gaussian distri-bution with average and variance [ J ij ] = 0 , [ J ij ] = J N . (2)The parameter J characterises the width of the Gaussian distribution. In its standard spin-glass setting the spinsare Ising variables and the model is the Sherrington-Kirkpatrick spin-glass. We will, instead, use continuousvariables, −√ N ≤ s i ≤ √ N with i = 1 , . . . , N , globally forced to satisfy (on average) a spherical constraint, (cid:80) Ni =1 s i = N , with N the total number of spins [19]. Such spherical constraint is imposed by adding a term H constr = z (cid:32) N (cid:88) i =1 s i − N (cid:33) (3)to the Hamiltonian, with z a Lagrange multiplier. The spins thus defined do not have intrinsic dynamics. Instatistical physics applications their temporal evolution is given by the coupling to a thermal bath, via a MonteCarlo 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 forglasses, with potential energy H pot = − (cid:80) i ...i p J i ...i p s i . . . s i p and p integer. Even more generally, the form(1) is one instance of a generic random potential V ( { s i } ) with zero mean and correlations [12, 13, 14] [ V ( { s i } ) V ( { s (cid:48) i } )] = − N V ( | (cid:126)s − (cid:126)s (cid:48) | /N ) (4)with V ( | (cid:126)s − (cid:126)s (cid:48) | /N ) = − J ( (cid:126)s · (cid:126)s (cid:48) /N ) .Similarly to what was done in [7] in the study of the p ≥ Hamiltonian dynamics, the model can be endowedwith conservative dynamics by changing the “spin” interpretation into a “particle” one. In this way, a kineticenergy [21, 22] H kin [ { ˙ s i } ] = m N (cid:88) i =1 ( ˙ s i ) , (5) an be added to the potential energy. The total energy of the Hamiltonian spherical p -spin model is then H syst = H kin + H pot + H constr . (6)This model represents a particle constrained to move on an N -dimensional hyper-sphere with radius √ N . Theposition of the particle is given by an N-dimensional vector (cid:126)s = ( s , . . . , s N ) and its velocity by another N -dimensional vector ˙ (cid:126)s = ( ˙ s , . . . , ˙ s N ) . The N coordinates s i are globally constrained to lie, as a vector, on thehypersphere with radius √ N . The velocity vector ˙ (cid:126)s is, on average, perpendicular to (cid:126)s , due to the sphericalconstraint. The parameter m is the mass of the particle.The generic set of N equations of motion for the isolated system is m ¨ s i ( t ) + z ( t ) s i ( t ) = (cid:88) j ( (cid:54) = i ) J ij s j ( t ) , (7) i = 1 , . . . , n , where the Lagrange multiplier needs to be time-dependent to enforce the spherical constraint inthe course of time.The initial condition will be taken to be { s i , ˙ s i } ≡ { s i (0) , ˙ s 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 differenttemperatures T (cid:48) .From Eq. (7) one derives an identity between the energy density and the Lagrange multiplier. By multiplyingthe equation by s i ( t ) and taking t → t lim t → t − m∂ t C ( t , t ) + z ( t ) = − e pot ( t ) . (8)The first term can be rewritten as m lim t → t − ∂ t C ( t , t ) = − m lim t → t − (cid:80) Ni =1 ˙ s i ( t ) ˙ s i ( t ) = − m (cid:80) Ni =1 ( ˙ s i ( t )) .Therefore z ( t ) = − e pot ( t ) + 2 e kin ( t ) . (9)The Lagrange multiplier takes the form of an action density, as a difference between kinetic and potential energydensities. Using now the conservation of the total energy, e f = e pot ( t ) + e kin ( t ) , z ( t ) = 2 e f − e pot ( t ) = − e f + 4 e kin ( t ) . (10)The p = 2 model belongs to a different universality class from the one of the p ≥ cases, in the sense that itsfree-energy landscape and relaxation dynamics are much simpler. It is, indeed, a model that resembles stronglythe large N , O ( N ) model for ferromagnetism. A hint on the simpler properties of its potential energy landscapeis given by the fact that the equations derived for general p simplify considerably for p = 2 . For example, thestatic and dynamic transitions occur at the same temperature T c = T d , and the number of metastable states isdrastically reduced. We recall these properties in the rest of this Section. The static properties of the p = 2 spherical model were elucidated in [19]. The trick is to project the spinvector (cid:126)s onto the basis of eigenvectors of the interaction matrix. One calls λ µ and (cid:126)v µ the µ -th eigenvalue andeigenvector of the matrix J ij , and s µ = (cid:126)s · (cid:126)v µ the projection of (cid:126)s on the eigenvector (cid:126)v µ . In terms of the latterthe Hamiltonian is not only quadratic but also diagonal. The extrema of the potential energy landscape andthe partition function can then be easily computed. In the thermodynamic limit, N → ∞ , the eigenvalues aredistributed according to the Wigner-Dyson semi-circle form [23] ρ ( λ µ ) = 12 πJ (cid:113) (2 J ) − λ µ θ (2 J − | λ µ | ) . (11)For finite N the distance between the largest and next to largest eigenvalues is order N − / /N / = N − / . .2.1 The potential energy landscape Let us label the eigenvalues of J ij in such a way that they are ordered: λ ≤ λ ≤ · · · ≤ λ N . We calltheir associated eigenvectors (cid:126)v µ with µ = 1 , . . . , N and we take them to be orthonormal, such that v µ = 1 . Weconsider the potential energy landscape H J ( { s i } , z ) with z taken as a variable.In the absence of a magnetic field, all eigenstates of the interaction matrix are stationary points of thepotential energy hyper-surface, ∂H J ∂s i (cid:12)(cid:12)(cid:12)(cid:12) (cid:126)s ∗ = − N (cid:88) j ( (cid:54) = i ) J ij s j + zs i | (cid:126)s ∗ ,z ∗ = 0 ∀ i , ⇒ (cid:126)s ∗ = ±√ N(cid:126)v µ , z ∗ = λ µ , ∀ µ . These stationary points are metastable states at zero temperature, apart from two of them that are the marginallystable ground states, and their number is linear in N , the number of spins. (The role of marginal stability in thephysical behaviour of condensed matter systems was recently summarised in [24].) These statements are shownin the way described in the next paragraph.The Hessian of the potential energy surface on each stationary point is ∂H J ∂s i ∂s j (cid:12)(cid:12)(cid:12)(cid:12) (cid:126)s ∗ ,z ∗ = − J ij + zδ ij | (cid:126)s ∗ ,z ∗ = − J ij + λ µ δ ij . (12)This matrix can be easily diagonalised and one finds D νη = ( − λ ν + λ η ) δ νη . Thus, on the stationary point, (cid:126)s ∗ = ±√ N(cid:126)v µ , the Hessian has one vanishing eigenvalue (for ν = µ ), µ − 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, µ − stabledirections and N − µ unstable directions. (In other words, the number of stable directions plus the marginallystable one is given by the index µ labelling the eigenvalue associated to the stationary state.) In conclusion,there are two maxima, (cid:126)s ∗ = ±√ N(cid:126)v , in general two saddles (cid:126)s ∗ = ±√ N(cid:126)v I with I = µ − stable directions and N − I unstable ones, with I running with µ as I = µ − and µ = 3 , . . . , N , and finally two (marginally stable)minima, (cid:126)s ∗ = ±√ N(cid:126)v N . 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 H J = − (cid:88) ij J ij s i s j + z (cid:32)(cid:88) i s i − N (cid:33) = − (cid:88) µ ( λ µ − z ) s µ − z N . (13)At each stationary point (cid:126)s ∗ = ±√ N(cid:126)v µ and z ∗ = λ µ this energy is H ∗ µ ≡ H J ( (cid:126)s ∗ = (cid:126)v µ ) = − (cid:88) i v µi (cid:88) j J ij v µj + z ∗ (cid:32)(cid:88) i ( v µi ) − N (cid:33) = − λ µ N . (14)Here we used the notation v µi to indicate the i th component of the µ th eigenvector (cid:126)v µ . The energy differencebetween the minima and the lowest saddles depends on the distribution of eigenvalues, a semi-circle law for theGaussian 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 ∂H J ∂s i (cid:12)(cid:12)(cid:12)(cid:12) (cid:126)s ∗ = − N (cid:88) j ( (cid:54) = i ) J ij s j + zs i − h i | (cid:126)s ∗ ,z ∗ = 0 , ∀ i , ⇒ s ∗ i = ( z ∗ − J ) − ij h j and z ∗ is fixed by imposing the spherical constraint on (cid:126)s ∗ . One then finds two solutions for the Lagrangemultiplier that lie outside the interval of variation of the eigenvalues of the J ij matrix: | z ∗ | > λ N . The stabilityanalysis shows that the stationary points are just one fully stable minimum and a fully unstable maximum. Theelimination 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 initself with implications in condensed matter physics, notably in glass theory [31, 24], but also claimed to play arole in string theory [32, 33], evolution [34] or other fields. The p = 2 spherical model provides a particularlysimple case in which the potential energy landscape can be completely elucidated. .2.2 The free-energy density This special (almost) quadratic model allows for the complete evaluation of its free-energy density for atypical realisation of the disordered exchanges. The traditional derivation of the disorder averaged free-energydensity can also be done using the replica method and a simple replica symmetric
Ansatz solves this problemcompletely. We recall how the two methods [19, 35] work in this Section.The partition function reads Z J = N (cid:89) i =1 (cid:90) ∞−∞ ds i e β (cid:80) i (cid:54) = j J ij s i s j πi (cid:90) c + i ∞ c − i ∞ dz e − βz ( (cid:80) Ni =1 s i − N ) where c is a real constant to be fixed below.It is convenient to diagonalise the matrix J ij with an orthogonal transformation and write the exponent interms of the projection of the spin vector (cid:126)s on the eigenvectors of J ij , s µ ≡ (cid:126)s · (cid:126)v µ . This operation can be done forany particular realisation of the interaction matrix. The new variables s µ are also continuous and unboundedand the partition function can be recast as Z J = N (cid:89) µ =1 (cid:90) ∞−∞ ds µ πi (cid:90) c + i ∞ c − i ∞ dz e (cid:80) Nµ =1 β ( λ µ − z ) s µ / βzN/ . (15)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 Z J = 12 πi (cid:90) c + i ∞ c − i ∞ dz e − N [ − βz/ N ) − (cid:80) µ ln[ β ( z − λ µ ) / ] . (16)In the saddle-point approximation valid for N → ∞ the Lagrange multiplier is given by (cid:104)(cid:104) k B T ( z sp − λ µ ) − (cid:105)(cid:105) (17)and this equation determines the different phases in the model. We indicate with double brackets the sum overthe eigenvalues of the matrix J ij that in the limit N → ∞ can be traded for an integration over its density: N N (cid:88) µ =1 g ( λ µ ) = (cid:90) dλ µ ρ ( λ µ ) g ( λ µ ) ≡ (cid:104)(cid:104) g ( λ µ ) (cid:105)(cid:105) . (18)Let us discuss the problem in the absence of a magnetic field. The high temperature solution to Eq. (17) z sp = z eq = T + J T (19)can be smoothly continued to lower temperatures until the critical point ( k B T c ) − = (cid:104)(cid:104) ( z sp − λ µ ) − (cid:105)(cid:105) (20)is reached where z sp touches the maximum eigenvalue of the J ij matrix, and it sticks to it for all T < T c = J : z sp = z eq = λ max = 2 J T ≤ T c . (21) T c is the static critical temperature. (A magnetic field with a component on the largest eigenvalue, (cid:126)h · (cid:126)v max (cid:54) = 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 multipliervalues, one verifies that it is in the high temperature phase, but it is not in the low temperature phase, where N (cid:88) µ =1 (cid:104) s µ (cid:105) = TT c N . (22)The way out is to give a macroscopic weight to the projection of the spin in the direction of the eigenvector thatcorresponds to the largest eigenvalue: s N = m √ N + δs N = (cid:115)(cid:18) − TT c (cid:19) N + δs N (23) ith (cid:104) δs N (cid:105) = 0 so that (cid:104) s N (cid:105) + N − (cid:88) ν =1 (cid:104) s ν (cid:105) = (cid:18) − TT c (cid:19) N + TT c N = N . (24)The thermal average of the projection of the spin vector on each eigenvalue vanishes in the high temperaturephase and reads (cid:104) s µ (cid:105) = (cid:40) [ N (1 − T /T c )] λ µ = λ max , λ µ < λ max , (25)below the phase transition (once we have chosen one of the ergodic components with the spontaneous symmetrybreaking of the (cid:126)s → − (cid:126)s invariance). The configuration condenses onto the eigenvector associated to the largesteigenvalue of the exchange matrix that carries a weight proportional to √ N . Going back to the original spinbasis, the mean magnetisation per site is zero at all temperatures but the thermal average of the square of thelocal magnetisation, that defines the Edwards-Anderson parameter, is not when T < T c : (cid:104) m i (cid:105) = 1 − T /T c ⇒ q EA ≡ [ (cid:104) m i (cid:105) ] J = 1 − T /T c with T c = J . (26)The order parameter q EA vanishes at T c and the static transition is of second order.The condensation phenomenon occurs for any distribution of exchanges with a finite support. If the distri-bution has long tails, as when the model is defined on a sparse random graph [36, 37, 38], the energy densitydiverges 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 replicasymmetric Ansatz . This
Ansatz corresponds to an overlap matrix between replicas Q ab = δ ab + q EA (cid:15) ab with (cid:15) ab = 1 for a (cid:54) = b and (cid:15) ab = 0 for a = b . When N → ∞ the saddle point equations fixing the parameter q EA yield above T c and a marginally stable solution with q EA = 1 − T /T c and identical physical properties to theones discussed above below T c .The equilibrium energy is given by e poteq = (cid:40) − J T (cid:104) − (cid:0) − TJ (cid:1) (cid:105) = ( k B T − λ max ) T < T c , − J T T > T c . (27)We added a superscript pot since in the modified model that we will study in this paper the total energy will alsohave a kinetic energy contribution. The entropy diverges at low temperatures as ln T , just as for the classicalideal gas, as usual in classical continuous spin models. The over-damped relaxation dynamics of the spherical p = 2 spin model (coupled to a Markovian bath) werestudied in [42, 43, 44, 26, 45]. One of the settings considered in these papers evolve a completely random initialcondition, { s i } , that corresponds (formally) to an infinite temperature initial state. The system is then subjectto an instantaneous temperature quench by changing the temperature of the bath to a final value T < + ∞ .Initial conditions drawn from equilibrium at temperature T (cid:48) < T c , and evolving in contact with a bath at thesame temperature, T = T (cid:48) , were considered in [44] and it was shown in this paper that equilibrium at the sametemperature is maintained ever after. A quench of the dissipative system from equilibrium at T (cid:48) < T c to anothersubcritical temperature T < T c was also studied in [44] and it was there shown that equilibrium at the targettemperature is achieved.The coupling to the bath is modeled with a stochastic equation of Langevin kind. This equation can beexactly solved in the basis of eigenvectors of the interaction matrix s µ ( t ) = s µ (0) e − λ µ t − (cid:82) t dt (cid:48) z ( t (cid:48) ) + (cid:90) t dt (cid:48) e − λ µ ( t − t (cid:48) ) − (cid:82) tt (cid:48) dt (cid:48) z ( t (cid:48)(cid:48) ) [ ξ µ ( t (cid:48) ) + h µ ( t (cid:48) )] (28)and the Lagrange multiplier z ( t ) can be fixed by imposing the spherical constraint C ( t, t ) = 1 N (cid:88) µ s µ ( t ) = 1 (29) t all times. This yields a self-consistent equation for z ( t ) . Notably, the asymptotic solution depends on thechoice of initial state, as we expose below. The applied field h µ is used to compute the linear response function R ( t , t ) = 1 N N (cid:88) i =1 δ (cid:104) s i ( t ) (cid:105) h δh i ( t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h =0 = 1 N N (cid:88) µ =1 δ (cid:104) s µ ( t ) (cid:105) h δh µ ( t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h =0 . (30)For quenches of initial conditions drawn from equilibrium at T (cid:48) → ∞ , and evolution in contact with a bathat temperature T > T c , the dynamics quickly approach equilibrium at the new temperature. The Lagrange mul-tiplier quickly converges to z eq = T + J /T . The correlation and linear response are invariant under translationsof 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 (cid:48) → ∞ , and evolution in contact with a bath attemperature T < T c , 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, ( t − t ) /t (cid:28) , and one non-stationaryfor long time differences ( t − t ) /t (cid:29) . The detailed time-dependence in the two regimes can be extractedusing the procedure sketched above. It yields the behaviour of the self correlation and linear response that scalein 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 largesteigenvalue 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 reachedin finite times with respect to N .For low temperature quenches from random initial conditions, the Lagrange multiplier approaches z f = 2 J as a power law, z ( t ) − z f (cid:39) − / (4 t ) . The slow approach to the asymptotic value is determinant to allow for thenon-stationary slow relaxation. The global correlation and linear response are computed from the spin solution s µ ( t ) and they can be cast as [44] C ( t , t ) = C st ( t − t ) + C ag ( t , t ) , (31) R ( t , t ) = R st ( t − t ) + R ag ( t , t ) , (32)with the stationary and a non-stationary terms linked by the FDT at the temperature of the bath R st ( t − t ) = − T dC st ( t − t ) d ( t − t ) , (33)and a modified FDT at an effective temperature T eff [51, 52] selected by the dynamics, R ag ( t − t ) = 1 T eff ( t , t ) ∂C ag ( t , t ) ∂t , (34)always with t ≥ t . In the asymptotic limit, the two terms added to form C and R evolve in different regimesin the sense that when one changes the other one is constant and vice versa. The limiting values of the twocontributions to the correlation function are C st (0) = 1 − q , lim t − t →∞ C st ( t − t ) = 0 , (35) lim t → t − C ag ( t , t ) = q , lim t (cid:29) t C ag ( t , t ) = 0 , (36)with the parameter q being equal to q = 1 − TJ . (37)This is the correct expression of the Edwards-Anderson parameter for the equilibrium low temperature solution,see Eq. (26), and once again one finds T c = J from q = 0 . The complete solution of the Langevin equationsallows one to deduce the exact scaling forms of the stationary and ageing contributions to the correlation andlinear response. These are C ag ( t , t ) = f C (cid:18) t t (cid:19) R ag ( t , t ) = t − / f R (cid:18) t t (cid:19) (38)with f C ( x ) and f R ( 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 ≥ cases [20] it is not constant but grows with time. ore precisely, it scales as T eff ( t , t ) (cid:39) t / f T ( t /t ) and it diverges asymptotically as t / . This implies, inparticular, that the ageing regime does not contribute to the asymptotic potential energy that, after a quenchto T < T c , reads e potasymp = − J (cid:20) T (1 − q ) (cid:21) = e poteq (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 (cid:48) < T c → T < T c , for example taking initial conditions in equilibriumat zero temperature, s µ (0) = √ Nδ µ,N and evolving them at T < T c , the Lagrange multiplier approaches z f = 2 J 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 J ij , can also be imple-mented in the case in which there is inertia and the differential equation has a second order time derivative. Thedynamics are recast into the ones of harmonic oscillators coupled by a self-consistent time-dependent Lagrangemultiplier. We will use this formulation in Sec. 5. Although a full analytical solution is not possible with thesecond time derivative, a performant numerical algorithm will allow us to monitor the evolution of the differentmodes. O ( N ) λφ model in the large N limit The λφ scalar field theory in d dimensions is defined by the Hamiltonian H = (cid:90) d d x (cid:20)
12 ( ∇ φ ) + r φ + λ φ (cid:21) , (40)where r and λ are two parameters. This model is the Ginzburg-Landau free energy for the local order parameterof the paramagnetic-ferromagnetic transition controlled by the parameter r going from r > to r < . Whenthe field is upgraded to a vector with N components and the limit N → ∞ is taken the quartic term, firstconveniently normalised by N , can be approximated by λ N φ (cid:55)→ λ N (cid:104) φ (cid:105) φ . The quantity z ( t ) = (cid:104) φ (cid:105) is notexpected to fluctuate and plays the role of the Lagrange multiplier in the spherical disordered model. Oncethis approximation made, the model becomes quadratic in the field and its statics and relaxation dynamics canbe easily studied. The only difficulty lies in imposing the self-consistent constraint that determines z ( t ) . Thecondensation phenomenon that we discussed in the disordered model is also present in the field theory and itcorresponds to a condensation on the zero wave vector mode. In the dynamic problem this corresponds to theprogressive approach to the homogeneous field configuration [9].The conserved energy dynamics of the λφ model, especially after sudden quenches, has been studied by anumber of groups. Details on the behaviour of the scalar problem, as well as a review of general equilibrationand pre-equilibration issues can be found in J. Berges’ Les Houches Lecture Notes [53], see also [54]. Thedynamics of the large N limit of the O ( N ) model was analysed in [55, 56, 57, 58, 59, 60]. More recent works userenormalisation group techniques to study the short time dynamics [61, 62] at the dynamic phase transition. In this Section we explain the relation between the Hamiltonian p = 2 disordered model and the integrablemodel of Neumann [15]. We start by recalling some basic properties of classical integrable systems in the senseof 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 formalismto analyse certain aspects of the long time dynamics of the system, and we set the stage for a future study ofthe eventual approach to a GGE ensemble. In classical mechanics, systems are said to be Liouville integrable if there exist sufficiently many well-behavedfirst 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 ore concrete terms, an integrable dynamical system consists of a N -dimensional phase space Γ together with N independent functions O , . . . , O N : Γ → R , such that the mutual Poisson brackets vanish, { O j , O l } = 0 for all j, l . (41)We will assume henceforth that the O i do not depend explicitly on time and that dO i /dt = 0 is equivalent to { H, O i } = 0 . Conventionally, the first function O is the Hamiltonian itself and the first constant of motion isthe energy. All other O i with i (cid:54) = 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 N − N = N in whichall configurations share the initial values of all the conserved quantities O i ( t ) = O i (0) . Under these conditionsHamilton’s equations of motion are solvable by performing a canonical transformation into action-angle variables ( I i , φ i ) with i = 1 , . . . , N such that the Hamiltonian is rewritten as ˜ H ( I ) and I i ( t ) = I i (0) , φ i ( t ) = φ i (0) + t ∂ ˜ H ( I ) ∂I i = φ k (0) + t ω i ( I ) . (42)The action functions I i are conserved quantities and we collected them in I in the dependence of the frequencies ω i and the Hamiltonian ˜ H . The remaining evolution is given by N circular motions with constant angularvelocities. Both deciding whether a system is integrable and finding the canonical transformation that leads tothe pairs ( I i , φ i ) are in practice very difficult questions. Whenever the system is integrable, and one knows theaction-angle pairs, the statement in Eq. (42) is part of the Liouville-Arnold theorem [65]. The model proposed by Neumann in 1850 describes the dynamics of a particle constrained to move on the N − dimensional sphere under the effect of harmonic forces [15]. The Hamiltonian is H = 14 N (cid:88) k (cid:54) = l L kl + 12 (cid:88) k a k x k (43)where the L kl are the elements of an angular momentum anti-symmetric matrix √ m L kl = x k p l − p k x l , (44)and p k and x k are phase space variables with canonical Poisson brackets { x k , p l } = δ kl . The global sphericalconstraint N (cid:88) k =1 x k = N (45)ensures that the motion takes place on S N − . Using the fact that L kk = 0 to rewrite the double sum in the firstterm 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 (cid:80) k (cid:54) = l L kl = m (cid:80) k,l L kl = 2 (cid:80) k x k (cid:80) l p l − (cid:80) k x k p k (cid:80) l x l p l . Imposing next the spherical constraint, thatalso implies (cid:80) k x k p k = 0 , the sum simply becomes (cid:88) k (cid:54) = l L kl = 2 Nm (cid:88) k p k . (46)We note that we added a factor /N in ifront of the kinetic energy in Eq. (43) in order to ensure that the twoterms 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 latteris written in the basis of eigenvectors of the interaction matrix J ij .The N − integrals of motion of this problem were constructed by K. Uhlenbeck [16] and more recentlyrederived by Babelon & Talon [18] with a separation of variables method. In a notation that is convenient forour application they read I k = x k + 1 N (cid:88) l ( (cid:54) = k ) L kl a k − a l = x k + 1 mN (cid:88) l ( (cid:54) = k ) x k p l + x l p k − x k p l x l p k a k − a l (47) In the sense that the gradients (cid:126) ∇ O i are linearly independent vectors on a tangent space to any point in Γ nd satisfy (cid:80) k I k = N and (cid:80) k a k I k = H . In the definition of our Hamiltonian and equations of motion weused a convention such that a k (cid:55)→ − λ µ (note the minus sign). After a trivial translation to the variables of the p = 2 spherical model we then have I µ = s µ + 1 mN (cid:88) ν ( (cid:54) = µ ) s µ p ν + s ν p µ − s µ p ν s ν p µ λ ν − λ µ , (48) (cid:80) µ I µ = 1 and (cid:80) µ λ µ I µ = − H . Let (cid:126)X = ( x , p , . . . , x N , p N ) be a generic point in phase space. The fact that the microcanonical measure ρ GME ( (cid:126)X ) = c − N (cid:89) j =1 δ ( I j ( (cid:126)X ) − i j ) , (49)with c = (cid:82) d (cid:126)X (cid:81) Nj =1 δ ( I j ( (cid:126)X ) − i j ) be sampled asymptotically is ensured by the Liouville-Arnold theorem [65], ifthe frequencies of the periodic motion on the torus are independent, that is, (cid:126)k · (cid:126)ω = 0 for (cid:126)k = ( k , . . . , k N ) withinteger k j has the unique solution (cid:126)k = 0 . One can call this ensemble the Generalized Microcanonical Ensemble.In principle, the Generalized Canonical Ensemble, commonly called Generalized Gibbs Ensemble (GGE), cannow be constructed from the Generalized Microcanonical Ensemble following the usual steps. The idea is tolook for the joint probability distribution of N extensive (as for the Hamiltonian in the usual case) constants ofmotion of a subsystem P ( i , . . . , i N ) di . . . di N . As in cases with just one conserved quantity, it is convenient tointerpret P as a probability over position and momenta variables, and write P GGE ( (cid:126)X ) = Z − ( λ , . . . , λ N ) e − (cid:80) j ζ j I j ( (cid:126)X ) . (50)This form can be derived under the same kind of assumptions used in the derivation of the canonical measure fromthe 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 ( { i , . . . , i N } ) = g ( { i (1)1 , . . . , i (1) N } ) g ( { i (2)1 , . . . , i (2) N } ) , (ii) additivity ofthe conserved quantities i (2) j = i j − i (1) j , (iii) small system 1 ( i (1) j (cid:28) i (2) j ) , (iv) constant inverse ‘temperatures’ ζ j ≡ ∂ i j S ( { i , . . . , i N } ) = k B ∂ i j ln g ( { i , . . . , i N } ) . An inspiring discussion along these lines appeared in [66].The conditions just listed imply a locality requirement on the i j s, otherwise (ii) and (iii) would be violated. Thisis similar to the requirement of having short-range interactions to derive the equivalence between the canonicaland microcanonical ensembles in standard statistical mechanics.In quenching procedures, the parameters ζ j should be determined by requiring that the expectation valueof each conserved quantity I j calculated on ρ GGE matches the (conserved) initial value I j (0 + ) (right after thequench): I j (0 + ) = (cid:90) d (cid:126)X I j ( (cid:126)X ) P GGE ( (cid:126)X ) . (51)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 / , with λ µ theeigenvalues of the random interaction matrix, the GGE measure would be P GGE ( (cid:126)X ) = Z − e − β f ( − (cid:80) µ λ µ I µ = Z − e − β f H = P GB ( (cid:126)X ) , (52)the Gibbs-Boltzmann one. Take now a generic function of the phase space variables A ( (cid:126)X ) that does not depend explicitly on time andis not conserved. Birkhoff’s theorem [67] states that its long-time average exists and reaches a constant, A ≡ lim τ →∞ τ (cid:90) t st + τt st dt (cid:48) A ( (cid:126)X ( t (cid:48) )) = cst (53) or τ sufficiently long and t st 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 begiven by the averages over the statistical measure P GGE : A = (cid:90) D (cid:126)X A ( (cid:126)X ) P GGE ( (cid:126)X ) . (54)Which are the observables for which this result should hold is an interesting question that needs to be answeredwith care.The GGE proposal [2, 3] and most, if not all, of its discussion appeared in the treatment and study of quantumisolated systems and, especially, the dynamics following an instantaneous quench performed as a sudden changein a parameter of the system’s Hamiltonian. A series of review articles are [4, 5, 6]. The main motivation for ourresearch project is to ask similar questions in the classical context, with the aim of disentangling the quantumaspects from the bare consequences of isolation and integrability. The fluctuation-dissipation theorem (FDT) [68] is a model independent equilibrium relation between thetime-delayed linear response of a chosen observable and its companion correlation function. In Gibbs-Boltzmannequilibrium this relation is independent of the specific system and observable and it only involves the inversetemperature of the system. For classical systems it admits simple expressions in the time and frequency domains R AB ( t − t ) = β∂ t C AB ( t − t ) θ ( t − t ) or Im ˆ R AB ( ω ) = − βω ˆ C AB ( ω ) . (55)Out of canonical equilibrium, the fluctuation-dissipation relations (FDR) between the linear response andthe correlation function have been used to quantify the departure from equilibrium [20]. Indeed, the possiblytime and observable dependent parameter that replaces β in far form equilibrium systems yields an effectivetemperature 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 domainIm ˆ R AB ( ω ) ω ˆ C AB ( ω ) = − β AB eff ( ω ) (56)that concretely defines the frequency dependent, and also possibly observable dependent, inverse effective tem-perature β AB eff .It was shown in [69, 70] that the Lagrange multipliers ζ j of the GGE, seen as inverse temperatures β j , of anumber of isolated integrable quantum systems which reach a stationary state can be read from the FDR’s ofproperly chosen observables β j = β eff ( ω j ) for all j . (57)In this paper we will show that this statement also applies to the classical integrable system that we analyse. We now enter the heart of our study and we consider the dynamics of the isolated system after different kindsof 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: lim t →∞ lim N →∞ (58) We start by giving a short description of steps that lead to the dynamic equations that couple linear responseand correlation function and fully characterise the evolution of the model in the N → ∞ limit. When dealing with the numerical data we used a Fourier transform convention such that R ( ω, t ) = (cid:80) k R ( t k + t , t ) e iωt k and C ( ω, t ) = (cid:80) k ( C ( t k + t , t ) − q ) cos( ωt k ) with t k the discrete times on which we collect the data points. The Fourier transform ofthe correlation is computed for t k > only, taking advantage of the long-time stationarity property C st ( − t ) = C st ( t ) . For this reason,there is no factor in the left-hand-side of the FDT in the frequency domain. .1.1 The Schwinger-Dyson equations In the N → ∞ limit, the only relevant correlation and linear response functions are C ab ( t , t ) = N − N (cid:88) i =1 [ (cid:104) s ai ( t ) s bi ( t ) (cid:105) ] , (59) C ab ( t ,
0) = N − N (cid:88) i =1 [ (cid:104) s ai ( t ) s bi (0) (cid:105) ] , (60) R ab ( t , t ) = N − δδh b ( t ) N (cid:88) i =1 [ (cid:104) s a ( h ) i ( t ) (cid:105) ] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h =0 , (61)for t , t > , where the infinitesimal perturbation h is linearly coupled to the spin H (cid:55)→ H − h (cid:80) Ni =1 s i attime t 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 t > t . The square brackets denote hereand everywhere in the paper the average over quenched disorder. The angular brackets indicate the averageover thermal noise if the system is coupled to an environment, and over the initial conditions of the dynamicssampled, say, with a probability distribution. When the coupling to the bath is set to zero, as we do in thispaper, 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 itwas applied to the relaxational p spin model in [75, 76]. The average over disorder now becomes non-trivial andneeds the use of the replica trick. The scripts a, b indicate then the replica indices a, b = 1 , . . . , N . For initialconditions in equilibrium the replica structure is replica symmetric (see Sec. 2.2 and [19]), with Q aa = 1 and Q a (cid:54) = b = q in (62)and q in in the paramagnetic state while q in (cid:54) = 0 in the condensed phase. This structure has an effect on theequation for the time-dependent correlation function that will keep the initial replica structure. There will betwo kinds of correlations with the initial condition C ( t, and C a (cid:54) =1 ( t, , (63)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 C ( t , t ) (64)as a representative of this group. The interpretation of the correlations C ( t , t ) and C ( t , can be given interms of real replicas. The former is the self-correlation between the configuration of one replica of the system { s i } ( t ) and the same replica evolved until a later time t , { s i } ( t ) . For this reason, we will eliminate thesubscript and call C ( t , t ) (cid:55)→ C ( t , t ) in most places hereafter. The latter is the correlation between onereplica of the system { σ i } (0) at the initial time and another replica of the system evolved until time t andrepresented by { s i } ( t ) . Although we could also write an evolution equation for the two-time C ( t , t ) we donot need it here since only C ( t , intervenes in the other equations. n the N → ∞ limit one derives the dynamical equations that read (cid:0) m∂ t + z ( t ) (cid:1) R ( t , t ) = J (cid:90) t t dt (cid:48) R ( t , t (cid:48) ) R ( t (cid:48) , t ) + δ ( t − t ) , (65) (cid:0) m∂ t + z ( t ) (cid:1) C ( t , t ) = J (cid:90) t dt (cid:48) R ( t , t (cid:48) ) C ( t (cid:48) , t )+ J (cid:90) t dt (cid:48) R ( t , t (cid:48) ) C ( t , t (cid:48) ) + JJ T (cid:48) (cid:88) a C a ( t , C a ( t , , (66) (cid:0) m∂ t + z ( t ) (cid:1) C a ( t ,
0) = J (cid:90) t dt (cid:48) R ( t , t (cid:48) ) C a ( t (cid:48) ,
0) + JJ T (cid:48) (cid:88) b C b ( t , Q ab , (67) z ( t ) = − m∂ t C ( t , t ) | t → t − + 2 J (cid:90) t dt (cid:48) R ( t , t (cid:48) ) C ( t , t (cid:48) )+ JJ T (cid:48) (cid:88) a ( C a ( t , . (68)The border conditions are C ( t ,
0) = C ( t , implying C (0 ,
0) = C (0 ,
0) = 1 ,C (0 ,
0) = Q , implying C (0 ,
0) = q in . (69)Note that the initial condition is not the same for all C b ( t , . It is equal to 1 for b = 1 and equal to q in for b (cid:54) = 1 . One can check that these equations coincide with the ones in [75, 76] when inertia is neglected, p = 2 and J = J , and a coupling to a bath is introduced. With respect to the equations studied in [7], they correspondto p = 2 and they have the extra ingredient of the influence of equilibrium initial conditions with a non-trivialreplica 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 → limit; they read JJ T (cid:48) (cid:88) a C a ( t , C a ( t ,
0) = JJ T (cid:48) [ C ( t , C ( t ,
0) + ( n − C ( t , C ( t , → n → JJ T (cid:48) [ C ( t , C ( t , − C ( t , C ( t , , (70) JJ T (cid:48) (cid:88) b C b ( t , Q b = JJ T (cid:48) (cid:2) C ( t ,
0) 1 + ( n − C b ( (cid:54) =1) ( t , q in (cid:3) → n → JJ T (cid:48) (cid:2) C a ( t , − C b ( (cid:54) =1) ( t , q in (cid:3) , (71) JJ T (cid:48) (cid:88) b C b ( t , Q ( a (cid:54) =1) b = JJ T (cid:48) (cid:2) C ( t , q in + C ( t ,
0) 1 + ( n − C b ( (cid:54) =(1 , ( t , q in (cid:3) → n → JJ T (cid:48) [ q in C ( t ,
0) + (1 − q in ) C ( t , . (72)Consequently, the terms induced in the equation for C ( t , t ) and C ( t , are different.With inertia and no coupled bath, the equal-time conditions are C ( t , t ) = 1 ,R ( t , t ) = 0 ,∂ t C ( t , t ) | t → t − = ∂ t C ( t , t ) | t → t +1 = 0 ,∂ t C ( t , t ) | t → t − = ∂ t C ( t , t ) | t → t +1 = 0 , (73) ∂ t C ( t , | t → + = 0 ,∂ t R ( t , t ) | t → t − = 1 m ,R ( t , t ) | t → t +1 = 0 , for all times t , t larger than or equal to + , when the dynamics start.We found convenient to numerically integrate the integro-differential equations using an expression of theLagrange multiplier that trades the second-time derivative of the correlation function into the total conserved nergy after the quench. Following the same steps explained in [7] we deduce z ( t ) = 2 e ( t ) + 4 J (cid:90) t dt (cid:48) R ( t , t (cid:48) ) C ( t , t (cid:48) ) + 2 JJ T (cid:48) (cid:2) C ( t , − C ( t , (cid:3) , (74)where we used Eq. (70) evaluated at t = t . It seems that we have simply traded z ( t ) by e ( t ) . Indeed, takingadvantage of the fact that for an isolated system e ( t ) = e f , a constant, the numerical solution of the evolutionequations becomes now easier since it does not involve the second time derivative of the correlation function. Inpractice, in the numerical algorithm we fix the total energy e f to its post-quench value derived in Sec. 4.3, andwe 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 incanonical Boltzmann equilibrium at any temperature T (cid:48) are (cid:0) m∂ t + z ( t ) (cid:1) R ( t , t ) = δ ( t − t ) + J (cid:90) t t dt (cid:48) R ( t , t (cid:48) ) R ( t (cid:48) , t ) , (75) (cid:0) m∂ t + z ( t ) (cid:1) C ( t , t ) = J (cid:90) t dt (cid:48) R ( t , t (cid:48) ) C ( t (cid:48) , t ) + J (cid:90) t dt (cid:48) R ( t , t (cid:48) ) C ( t , t (cid:48) )+ JJ T (cid:48) [ C ( t , C ( t , − C ( t , C ( t , , (76) (cid:0) m∂ t + z ( t ) (cid:1) C ( t ,
0) = J (cid:90) t dt (cid:48) R ( t , t (cid:48) ) C ( t (cid:48) ,
0) + JJ T (cid:48) [ q in C ( t ,
0) + (1 − q in ) C ( t , , (77) z ( t ) = − m∂ t C ( t , t ) | t → t − + 2 J (cid:90) t dt (cid:48) R ( t , t (cid:48) ) C ( t , t (cid:48) )+ JJ T (cid:48) [ C ( t , − C ( t , . (78)High and low temperature initial states are distinguished by q in = 0 for T (cid:48) > J , and q in = 1 − T (cid:48) /J for T (cid:48) < J ,respectively. The equation for C ( t , is just the one for C ( t , t ) evaluated at t = 0 so we do not write itexplicitly. In order to ensure constant energy dynamics we set J = J and m = m in this Subsection. We verify thatthe equations consistently conserve the equilibrium conditions. Moreover, we derive a number of properties ofthe linear response function that will be useful in the analysis of the instantaneous quenches as well. The equation for C ( t, admits the solution C ( t ,
0) = q in = cst in the case in which no quench is performed.Indeed, setting C ( t ,
0) = q in in Eq. (77) one has z ( t ) q in = J q in (cid:90) t dt (cid:48) R ( t , t (cid:48) ) + J T (cid:48) [ q in C ( t ,
0) + (1 − q in ) q in ] . (79)This equation has the solution q in = 0 , the one of the paramagnetic phase, and a non-vanishing q in (cid:54) = 0 solutionrelevant in the ordered phase. Let us now focus on the case q in (cid:54) = 0 . Using FDT, a property of equilibrium, theintegral can be performed, the contribution from t (cid:48) = 0 cancels the first term in the square brackets, and theone from t (cid:48) = t combines with the second term in the square bracket; the ensuing equation simplifies to read z ( t ) = 2 J T (cid:48) (1 − q in ) = z f . (80)Therefore, z ( t ) is also a constant. The equation for z ( t ) , using FDT, becomes z ( t ) = T (cid:48) + J T (cid:48) [ C ( t , t ) − C ( t , J T (cid:48) [ C ( t , − q ]= T (cid:48) + J T (cid:48) (1 − q ) = z f . (81)The two equations yield the low-temperature values z f = 2 J and q in = 1 − T (cid:48) /J .In equilibrium C ( t, and z ( t ) remain constant and equal to their initial values, q in and z f . .2.2 The linear response in the frequency domain Knowing that z ( t ) remains constant in equilibrium, one can easily analyse the response equation in Fourierspace. The equation that determines its dynamical evolution is transformed into ˆ R ( ω ) = 1 − mω + z f − Σ( ω ) . (82)For this model Σ( ω ) = J ˆ R ( ω ) and then ˆ R ( ω ) = 12 J (cid:104) ( − mω + z f ) ± (cid:112) ( − mω + z f ) − J (cid:105) . (83) ˆ R ( ω ) , and also R ( t ) , are independent of temperature for T (cid:48) < T c = J while they depend on temperature through z f for T (cid:48) > T c = J .The terms under the square root can be more conveniently written as functions of the special values mω ± = z f ± J (84)and the linear response is recast as ˆ R ( ω ) = 12 J (cid:20) − mω + z f ± m (cid:113) ( ω − − ω )( ω − ω ) (cid:21) . (85)The imaginary and real parts of ˆ R ( ω ) are thenIm ˆ R ( ω ) = 12 J (cid:40) − m (cid:113) ( ω − ω − )( ω − ω ) for ω − ≤ ω ≤ ω + otherwise (86)Re ˆ R ( ω ) = 12 J z f − mω − m (cid:113) ( ω − − ω )( ω − ω ) for ω ≤ ω − z f − mω for ω − ≤ ω ≤ ω + z f − mω + m (cid:113) ( ω − ω − )( ω − ω ) for ω + ≤ ω (87)(note the unusual choice of sign for the imaginary part that we adopted.) In terms of the physical parameters, ˆ R ( ω ) is real for | − mω + z f | > J . In the low temperature phase, since z f = 2 J , this implies ω − = 0 andthe imaginary part of the linear response is gapless. One can easily check that | ˆ R ( ω ) | = 1 /J in the interval ω − ≤ ω ≤ ω + . Away from this interval the modulus of the linear response is a complicated function of thefrequency.The zero frequency linear response ˆ R ( ω = 0) = (cid:90) ∞ dτ R ( τ ) = 12 J (cid:104) z f ± (cid:113) z f − J (cid:105) , (88)with τ a time delay, must be a real quantity and this form is a manifestation of the condition z f ≥ J . Thestatic susceptibility is then given by χ st = (cid:40) /J for T (cid:48) < T c = J /T (cid:48) for T (cid:48) > T c = J (89)with the result in the second line being ensured by the choice of minus sign in front of the square root. Thefrequency-dependent linear response can then be transformed back to real-time and thus get its full time-evolution.In the lower limit of the spectrum the imaginary part of the linear response goes asIm ˆ R ( ω − + (cid:15) ) ∼ − J (cid:114) mJ (cid:15) at T < T c − (2 mJω − ) / J √ (cid:15) at T > T c as (cid:15) → + (90)while in the upper limit it vanishes asIm ˆ R ( ω + − (cid:15) ) ∼ − (2 mJω + ) / J √ (cid:15) as (cid:15) → + (91)with the corresponding ω ± at T > T c or T < T c . .2.3 The correlation functions The FDT in the frequency domain Im ˆ R ( ω ) = − βω ˆ C ( ω ) implies that ˆ C ( ω ) should vanish in the same frequencyintervals in which the linear response is real. In the ω = 0 case the linear response is real and the FDT as writtenabove only imposes that ˆ C ( ω = 0) = (cid:82) ∞−∞ dt C ( t ) must be finite. One has to bear in mind that in the cases inwhich the correlation function approaches a non-vanishing constant q asymptotically the Fourier transform tobe computed is the one of C st ( t − t ) = C ( t − t ) − q with respect to t − t and the integral should then yield C st ( ω = 0) = (cid:82) ∞−∞ dt C st ( 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 C st ( t − t ) = C ( t − t ) − q in can be treated following the steps explained in [9] and App. A.1.3 (that do notassume FDT) and one finds C st ( ω ) = J C st ( ω ) | ˆ R ( ω ) | . (92)This equation has two solutions, either C st ( ω ) = 0 or | ˆ R ( ω ) | = 1 /J . The latter holds for frequencies in theinterval ω ∈ [ ω − , ω + ] . Outside of this interval C st ( ω ) must vanish.By taking the derivative of Eq. (76) with respect to t one readily checks that it equals Eq. (75) if the FDTbetween R and C is satisfied for all times and ∂ t C ( t ,
0) = 0 implying C ( t ,
0) = q in . (93)This last condition is a property of equilibrium as we have already discussed.Having established that C is a constant, Eq. (77) enforces q in = 0 , the high temperature initial value, or z f = J T (cid:48) [1 − C ( t , J T (cid:48) [ C ( t ,
0) + (1 − q in )] = 2 J , (94)the low temperature Lagrange multiplier.Concerning the correlation function C ( t , , we write it as C ( t ,
0) = C st ( t ,
0) + q allowing for a non-vanishing asymptotic value, q , and taking C st ( t , such that it vanishes in the long t limit. The equation (76)evaluated at t = 0 is then rewritten as (cid:0) m∂ t + z f (cid:1) [ C st ( t ,
0) + q ] = J (cid:90) t dt (cid:48) R ( t , t (cid:48) )[ C st ( t (cid:48) ,
0) + q ] + J T (cid:48) [ C st ( t ,
0) + q − q ] . (95)This equation has three terms that do not depend on C st ( t , explicitly z f q − J q lim t →∞ (cid:90) t dt (cid:48) R ( t , t (cid:48) ) − J T (cid:48) ( q − q ) (96)and their sum should vanish in the long t limit. It trivially does for high temperature initial states since q = q in = 0 and, in equilibrium at low temperatures, we can assume q = q in , use FDT, and find z f = 2 J T (cid:48) (1 − q in ) = 2 J (97)confirming the assumption q = q in . The remaining equation fixes the time-dependence of C ( t , . For the sake of completeness, we compute the energy variation due to a simultaneous change of the mass m → m and the variance of the interaction strengths J → J . In the applications and numerical tests we willfocus on the latter changes only.The kinetic energy density before the quench is e kin (0 − ) = 1 N N (cid:88) i =1 m s i (0 − )) = T (cid:48) , (98)the last equality being due to the fact that we take equilibrium initial conditions at temperature T (cid:48) . Thepotential energy density before the quench depends on the system being paramagnetic or condensed initially: e pot (0 − ) = − J T (cid:48) (cid:2) − q ( T (cid:48) /J ) (cid:3) = − J T (cid:48) (cid:34) − (cid:18) − T (cid:48) J (cid:19) (cid:35) condensed e pot (0 − ) = − J T (cid:48) paramagnetic (99) he kinetic energy density right after the quench is e kin (0 + ) = 1 N N (cid:88) i =1 m s i (0 + )) (100)and, since the velocities do not change in the infinitesimal interval taking from − to + , e kin (0 + ) = 1 N N (cid:88) i =1 m s i (0 + )) = 1 N N (cid:88) i =1 m s i (0 − )) = mm T (cid:48) . (101)The post-quench potential energy density can be estimated from the relation between the Lagrange multiplierand the energy e pot (0 + ) = − z (0 + ) + e kin (0 + ) (102)and the equation for z ( t ) (see Sec. 4.1) z ( t ) = − m lim t → t − ∂ t C ( t , t ) + JJ T (cid:48) [ C ( t , − C ( t , . (103)Thus e pot (0 + ) = − JJ T (cid:48) lim t → + [ C ( t , − C ( t , . (104)Assuming that the spin configuration did not change between the infinitesimal time step going from 0 to + , e pot (0 + ) = − JJ T (cid:48) (cid:34) − (cid:18) − T (cid:48) J (cid:19) (cid:35) condensed e pot (0 + ) = − JJ T (cid:48) paramagnetic (105)All the values just derived imply the changes in the total energy ∆ e tot = ∆ e kin + ∆ e pot = (cid:18) mm − (cid:19) T (cid:48) − J ( J − J )2 T (cid:48) (cid:34) − (cid:18) − T (cid:48) J (cid:19) (cid:35) , ∆ e tot = ∆ e kin + ∆ e pot = (cid:18) mm − (cid:19) T (cid:48) − J ( J − J )2 T (cid:48) , (106)respectively.We will concentrate on potential energy quenches only, and we will trace the phase diagram using theparameters y = T (cid:48) J x = JJ . (107) x > corresponds to energy extraction and x < to energy injection. We will show that the parameter spaceis split in several sectors displaying fundamentally different 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 notalways satisfied by the actual evolution found with the numerical integration, we deduce some properties of theLagrange multiplier, linear response and correlation function, in the long time limit. In this Section we state theassumptions, we summarise the results, and we leave most details of how these are derived to App. A.Consider the system in equilibrium at T (cid:48) with parameters J , m 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 theLagrange multiplier, reach a constant. Later we will further suppose that (in most cases) the correlation functionbecomes, itself, invariant under time-translations. Finally, we will explore in which circumstances a fluctuation-dissipation theorem (FDT) can establish with respect to a temperature T f for all time-delays, in stationarycases, 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 thenext Section, do not apply to all quenches. Still, we find useful to explore their consequences and derive fromthem a set of relations between the control parameters for which special behaviour arises, that we will later putto the numerical test. .4.1 Asymptotic values in a steady state Let us assume that the limiting value of the Lagrange multiplier is a constant lim t →∞ z ( t ) = z f . (108)The limit of the correlation function can be zero if the system decorrelates completely at sufficiently longtime-delays or non-zero if it remains within a confined state. We therefore call q = lim t − t →∞ lim t →∞ C ( t , t ) (109)the asymptotic value of the full two time correlation, or its stationary part in possible ageing cases, after thequench. Similarly, q = lim t →∞ C ( t , , (110) q = lim t →∞ C ( t , . (111) The equation for the linear response function does not depend explicitly on the pre-quench parameters, itdoes only on the post-quench mass m and interaction strength J . (An implicit dependence on the initial stateis not excluded, since the value taken by the Lagrange multiplier may depend on it.) The analysis that wedeveloped for the constant energy dynamics applies to the sudden quench case too. The solution of the responseequation in Fourier space yields ˆ R ( ω ) = 12 J (cid:104) ( − mω + z f ) ± (cid:112) ( − mω + z f ) − J (cid:105) . (112)with the zero frequency value ˆ R ( ω = 0) = 12 J (cid:16) z f − (cid:113) z f − J (cid:17) . (113)From the numerical solution of the full equations that we will present in Sec. 6 we infer that for m = m and J (cid:54) = J z f = (cid:40) J x > y ,T (cid:48) + J /T (cid:48) x < y , (114)with, we recall, x = J/J and y = T (cid:48) /J . Replacing in Eq. (113) one notices that the minus sign has to beselected for x < y at low frequency and χ st ≡ ˆ R ( ω = 0) = (cid:40) /T (cid:48) x < y , /J x > y , (115)where we called χ st , as a static susceptibility, the zero frequency response. After the quench Im ˆ R ( ω ) is non-zeroin a finite interval of frequencies [ ω − , ω + ] with mω ± = ( z f ± J ) as in Eq. (84) and z f taking the values inEq. (114). Therefore, Im ˆ R ( ω ) 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 time-independent 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 getthe full time dependence of the linear response function. In the numerical Section we will compare this functionalform, named R st , to the outcome of the full integration of the dynamic equations. From the relation between z and the energies, the conservation of energy, and Birkhoff’s theorem, z = 2 e kin − e pot , e tot (0 + ) = e kin + e pot , (116) here the overlines represent a long-time average defined in Eq. (53). Note that the Lagrange multiplier takesthe 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 fourrelevant regions of the phase diagram parametrised by x = J/J and y = T (cid:48) /J .(I) y > and x < y ( q in = 0 , z f = T (cid:48) + J /T (cid:48) , χ st = 1 /T (cid:48) ) e kin = − JJ T (cid:48) + 2 T (cid:48) + J T (cid:48) = J (cid:18) − xy + 2 y + x y (cid:19) e pot = − JJ T (cid:48) − J T (cid:48) = J (cid:18) − xy − x y (cid:19) (117)(II) y > and x > y ( q in = 0 , z f = 2 J , χ st = 1 /J ) e kin = − JJ T (cid:48) + T (cid:48) + 2 J = J (cid:18) − xy + y + 2 x (cid:19) e pot = − JJ T (cid:48) + T (cid:48) − J = J (cid:18) − xy + y − x (cid:19) (118)(III) y < and x > y ( q in (cid:54) = 0 , z f = 2 J , χ st = 1 /J ) e kin = T (cid:48) + T (cid:48) JJ = J ( y + yx )4 e pot = T (cid:48) + T (cid:48) JJ − J = J ( y + yx − x ) (119)(IV) y < and x < y ( q in (cid:54) = 0 , z f = T (cid:48) + J /T (cid:48) , χ st = 1 /T (cid:48) ) e kin = JT (cid:48) J − J + 2 T (cid:48) + J T (cid:48) = J (cid:18) xy − x + 2 y + x y (cid:19) e pot = JT (cid:48) J − J − J T (cid:48) = J (cid:18) xy − x − x y (cid:19) (120)The minimum potential energy density, e pot = − J is realised for T (cid:48) = 0 in Sector III.We stress that we have found these results without using FDT and they can therefore hold out of thermalequilibrium. We will investigate later which other conditions impose the use of FDT, a strong Gibbs-Boltzmannequilibrium condition. We can identify a kinetic temperature from the kinetic energy densities derived in the previous Subsectionby simply imposing T kin = 2 e kin . This operation leads to T ( I )kin = − JJ T (cid:48) + T (cid:48) + J T (cid:48) , T ( II )kin = − JJ T (cid:48) + T (cid:48) J ,T ( III )kin = T (cid:48) (cid:18) JJ (cid:19) , T ( IV )kin = JT (cid:48) J − J + T (cid:48) + J T (cid:48) , (121)in the four Sectors of the phase diagram distinguished in the previous Subsection. In App. A.3.1 we explain how we can exploit the conservation of the total energy, under the assumption thatthe asymptotic kinetic and potential energies take the form of Gibbs-Boltzmann equilibrium paramagnetic andcondensed equilibrium phases at a single temperature T f to fix its value. This means that we require e kin = T f and e pot = − J T f paramagnetic or − J T f (1 − q ) condensed or two-step , (122) ith q = 1 − T f /J . For x > y we find T f J = J J (cid:20) J − T (cid:48) − JJ T (cid:48) (cid:21) − for ≤ y (IIa) ,JJ + T (cid:48) J − J T (cid:48) = x + y − x y for ≤ y (IIb) ,T (cid:48) J (cid:18) JJ (cid:19) = y x ) for y ≤ (III) . (123)The roman numbers between parenthesis refer to the cases listed in Eqs. (117)-(120) and to the four sectors in thephase diagram in Fig. 5. The difference between (IIa) and (IIb) is that in the first case we used a paramagneticpotential energy and in the second case a condensed one at T f . The values for x < y are given in App. A.3.1.One can check that in the cases (IIb) and (III) T f = T kin , the kinetic temperatures given in Eq. (121). Instead,final and kinetic temperatures do not coincide in (IIa) nor in (I) and (IV). This fact excludes the possibility ofreaching Gibbs-Boltzmann equilibrium in (I) and (IV), that is, for x < y and it leaves the possibility open in(II) and (III) at the price of considering a potential energy density with a non-vanishing q = 1 − T f /J . Limits of validity
The two bounds ≤ T f ≤ J (124)serve to find special curves on the phase diagram. The first bound is natural since we do not want to have anegative kinetic energy. The second one ensures that q ≥ . The implications of these bounds, that are examinedin App. A.3.2, are y ≤ √ x for y ≥ , xx + 1 for y ≤ . (125)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 → , could only exist below the piecewise curve y ( x ) . One can simply checkthat the piece for y ≤ lies in Sector IV, since y < x , and it is therefore irrelevant given that we have alreadyshown that there cannot be a single temperature scenario in this Sector. The limit then moves to x = y for y < . 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 possiblefor y < and x > y (III), and below the curve y = √ x for y > in II. It is instructive to work out the limitingvalues of T f /J and q = 1 − T f /J on the borders of the region y < and y < x (III) of the phase diagram, andthe no-quench case x = 1 : T f J = y = 0 x ( x + 1)2 x = yx + y y = 1 y x = 1 q = y = 01 − x x = yx − x y = 11 − y x = 1 (III) (126)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 whatwe call Sector III. Moreover, the approximate asymptotic analysis of the mode dynamics of the finite N systemwill lead to T µ = T f in Sector III, see Eq. (183).On the limiting curve in Sector II, y = √ x for y > , T f = J and q = 0 .As one could have intuitively expected, T f > T (cid:48) for energy injection quenches ( x < ) and T f < T (cid:48) for energyextraction quenches ( x > ). .4.6 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 linearresponse to the correlation R ( t − t ) = − T f dC ( t − t ) d ( t − t ) θ ( t − t ) . (127) Static susceptibility
The values of the zero frequency linear response ˆ R ( ω = 0) = (cid:90) ∞ dτ R ( τ ) = 1 T f (1 − q ) , (128)where one must recall that τ is a time-difference, force ˆ R ( ω = 0) = /T (cid:48) for x < y ⇒ T f = T (cid:48) if q = 0 (I) & (IV) /J for x > y ⇒ T f = J if q = 0 (IIa) /J for x > y ⇒ T f = J/ (1 − q ) if q (cid:54) = 0 (IIb) & (III) (129)The result in the last line, valid for (IIb) and (III), does not put any additional constraint on T f . Instead, the otherconditions 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 < (sector III) or below the special curve y = √ x , for y > (lying inside sector II). Whether this is realisedor not needs to be investigated numerically. Ansatz
One can also look for a two step solution with similar characteristics to the ageing one found for dissipativedynamics [44] and summarised in Sec. 2.3. Asking for the relation between correlation and linear response inEq. (127) to hold in a stationary regime of relaxation in which C decays from to q , and that the effectivetemperature T eff characterising the second regime of decay from q to diverges, one recovers z f = 2 J , χ st = 1 J , q = 1 − T f J and q = q = 0 (130)with the same T f and q (cid:54) = 0 as in Eq. (123). From the analysis of the equation ruling the two-time correlation function, assuming stationarity and hencea dependence on time-difference only, one deduces (the details of the derivation are given in App. A.1.3, seealso [9] for a general treatment) ˆ C st ( ω ) = J ˆ C st ( ω ) | ˆ R ( ω ) | (131)where ˆ C st ( ω ) is the Fourier transform of the time-varying part (subtracting the possible non-vanishing asymptoticvalue q ). Note that the relation between the correlation and the linear response is the same as the one that wederived in the constant energy no-quench problem. It implies ˆ C st ( ω ) (cid:54) = 0 and | ˆ R ( ω ) | = 1 /J or ˆ C st ( ω ) = 0 , (132)independently of the control parameters. Below we check numerically that these relations are satisfied in variousquenches. In particular, from the analytic form of ˆ R ( ω ) one can easily see that | ˆ R ( ω ) | = 1 /J in the frequencyinterval in which the linear response is complex.Importantly enough, we cannot yield an explicit analytic form of ˆ C ( ω ) since it is factorised on both sides ofthe identity (131). We are forced to go back to the full set of dynamic equations and solve them numerically toget insight on the behaviour of C ( t , t ) . .4.9 Numerical results preview We will see that a state with a single T f characterising the fluctuation-dissipation relation is reached numer-ically 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 > (within sector II).In all other cases no equilibrium results à la Gibbs-Boltzmann are found for the global observables (correlationfunctions, linear response functions, kinetic and potential energies) but a different statistical description, of ageneralised kind, should be adopted. In particular, in Sector III where the conditions derived from the energyconservation and static susceptibility allowed for a single temperature scenario, the full solution of the completeset of questions will prove that this is not realised. A detailed explanation is given in the numerical section ofthe paper and in the analysis of the N − integrals of motion presented in Sec. 7.We recall that the p ≥ strongly interacting case behaves very differently [7]. On the one hand, equilibriumtowards a proper paramagnetic state, and within confining metastable states, were reached in two sectors of itsdynamic phase diagram. On the other hand, an ageing asymptotic state in a tuned regime of parameters wasalso found for more than two spin interactions in the potential energy. In the p = 2 integrable model we do notfind an ageing asymptotic state. Moreover, Gibbs-Boltzmann equilibrium is achieved in the two very particularcases listed above only. In this Section we describe how the finite system size dynamics can be solved by using a convenient basis inwhich the evolution becomes the one of harmonic oscillators coupled only through the Lagrange multiplier. Weshow that these oscillators decouple under the assumption z ( t ) → z f allowing for a simple approximate solutionof the problem that can, however, be relevant for N → ∞ only. We then explain a way to numerically solve thedynamics for finite N . Take a system with finite N . The post-quench matrix J ij has µ = 1 , ..., N eigenmodes with eigenvalues λ µ .If we denote s µ ( t ) = (cid:126)s ( t ) · (cid:126)v µ (133)the projection of the spin vector in the direction of the µ -th eigenvector of J ij , the N rotated equations of motionread m ¨ s µ ( t ) + ( z ( t ) − λ µ ) s µ ( t ) = 0 . (134)This set of equations has to be complemented with the initial conditions s µ (0) and ˙ s µ (0) . They are very similarto the equations for a parametric oscillator, the difference being that, in our case, the time-dependent frequencydepends on the variables via the Lagrange multiplier. Furthermore, they are identical to the equations of theNeumann’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 (cid:126)s using (cid:126)s ( t ) = (cid:80) µ s µ ( t ) (cid:126)v µ . In particular, the correlation function is given by C J ( t , t ) = 1 N (cid:88) µ (cid:104) s µ ( t ) s µ ( t ) (cid:105) , (135)where the subscript J means that the result depends, in principle, on the interaction matrix chosen, and theangular brackets represent an average over initial conditions. One could then perform the disorder average oranalyse the self-averageness properties of the correlation in different time regimes. The J dependence shoulddisappear in the N → ∞ limit.Since we are interested in an uniform interaction quench, it is easy to see that λ (0) µ = J J λ µ . (136) .2 Behaviour under stationary conditions Let us assume that the system reaches stationarity and that the Lagrange multiplier approaches a constant lim t →∞ z ( t ) = z f . (137)In order to simplify the notation, in the rest of this section we will measure time with respect to a referencetime t st for which the stationary regime for z ( t ) has already been established. (We insist upon the fact that thisassumption can only hold for N → ∞ .)The equation of motion of each mode becomes m ¨ s µ ( t ) + ( z f − λ µ ) s µ ( t ) = 0 (138)and can be thought of as Newton’s equation for the mode Hamiltonian H µ = 12 m ˙ s µ + 12 mω µ s µ (139)with ω µ ≡ ( z f − λ µ ) /m . This equation has three types of solutions depending on the sign of ω µ : ω µ > s µ ( t ) = s µ ( t st ) cos [ ω µ ( t − t st )] + ˙ s µ ( t st ) ωµ sin [ ω µ ( t − t st )] ,ω µ = 0 s µ ( t ) = s µ ( t st ) + ˙ s µ ( t st )( t − t st ) ,ω µ < s µ ( t ) = (cid:16) s µ ( t st ) − ˙ s µ ( t st ) | ω µ | (cid:17) e −| ω µ | ( t − t st ) + (cid:16) s µ ( t st ) + ˙ s µ ( t st ) | ω µ | (cid:17) e | ω µ | ( t − t st ) , (140)that is to say, oscillatory solutions with constant amplitude in the first case, diffusive behaviour in the interme-diate case and exponentially diverging solutions in the last case. We insist upon the fact that the initial timehere is taken to be t st 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, z eq > λ N and such a fully oscillating behaviour is expected.In equilibrium in the low temperature condensed phase z eq = λ max = λ N for N → ∞ and the µ = N modeshould grow linearly in time while all other modes should oscillate with frequency ω µ = (cid:112) ( z f − λ µ ) /m . Theamplitude of each mode is determined by the initial conditions, that are actually matching conditions at time t st , when stationarity is reached in this case. (Recall that λ N − is at distance N − / from λ N [23]. This meansthat, under the assumption z f → λ N , z f − λ N − = λ N − λ N − (cid:39) N − / and there will be almost diffusive modesclose 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 greaterthan λ N , although very close to it. The diffusive behaviour of the N th mode (in the N → ∞ limit) would beobtained as the limit of zero frequency of a (finite N ) oscillating N th mode. At variance with the N → ∞ approach, the finite size study allows to access the details of the dynamics ofeach mode. In this section we define some mode-observables that will provide valuable information. Of particularinterest are the mode energies, which can be defined as (cid:15) kin µ ( t ) = m (cid:104) ˙ s µ ( t ) (cid:105) ,(cid:15) pot µ ( t ) = 12 ( z ( t ) − λ µ ) (cid:104) s µ ( t ) (cid:105) ,(cid:15) tot µ ( t ) = e kin µ ( t ) + e pot µ ( t ) . (141)Note that in the analysis of the N → ∞ model the potential energy density is e pot = − / (2 N ) (cid:80) µ λ µ (cid:104) s µ (cid:105) without the term proportional to z ( t ) . For this reason we use here the different symbol (cid:15) pot µ for the modepotential energies that include the term proportional to z ( t ) . The values of these energies at t = 0 − are givenby the fact that all modes are in equilibrium at the same temperature: (cid:15) tot µ (0 − ) = 2 (cid:15) kin µ (0 − ) = 2 (cid:15) pot µ (0 − ) = T (cid:48) . (142) mmediately after the quench, they are (cid:15) kin µ (0 + ) = m (cid:104) ˙ s µ (0 + ) (cid:105) , (cid:15) pot µ (0 + ) = T (cid:48) z (0 + ) − λ µ z (0 − ) − λ (0) µ . (143)In order to study the eventual thermalisation of the system, we can define an effective time dependent modetemperature through the total mode energy T µ ( t ) ≡ (cid:15) tot µ ( t ) (144)based on the fact that the modes are (quasi) decoupled. Whenever the system enters a stationary regime inwhich z ( t ) is constant, see Section 5.2, the mode temperatures T µ are independent of time, since the systembehaves as a collection of non-coupled harmonic oscillators. We can also define mode temperatures using thekinetic and potential mode energies that oscillate around their mean if we take their time-average [67] e kin , pot µ = lim (cid:28) τ τ (cid:90) t st + τt st dt (cid:48) e kin , pot µ ( t (cid:48) ) (145)and propose T kin , pot µ ≡ (cid:15) kin , pot µ . (146)In the stationary regime, as shown in Section 5.2, the different mode temperatures should verify T kin µ = T pot µ . (147)Other useful observables are the time-delayed mode correlation functions C µ ( t , t ) = (cid:104) s µ ( t ) s µ ( t ) (cid:105) (148)and the mode linear response functions R µ ( t , t ) = δ (cid:104) s µ ( t ) (cid:105) h δh µ ( t ) (cid:12)(cid:12)(cid:12)(cid:12) h =0 (149)that is defined and measured as follows. If we add an external field h µ linearly coupled to each mode s µ , theequations of motion are modified into m ¨ s µ + ( z ( t ) − λ µ ) s µ − h µ = 0 (150)and its solution reads s (inhom) µ ( t ) = s (hom) µ ( t ) + s Pµ ( t ) , (151)where s (hom) µ ( t ) is the solution to the Newton equation without the external field and s Pµ ( t ) is a particular solutionof the inhomogeneous problem with initial condition s Pµ (0) = 0 and ˙ s Pµ (0) = 0 . The linear response function R µ ( t , t ) ≡ δs µ ( t ) /δh µ ( t ) | h =0 of the mode µ can be defined equivalently through s Pµ ( t ) = (cid:90) t dt R µ ( t , t ) h µ ( t ) . (152)In practice, to measure the linear response function numerically we apply a small external field localised intime h µ ( t ) = h δ ( t − t ) and we solve the inhomogeneous problem to obtain s (inhom) µ ( t ) . Using Eq. (152) weobtain the linear response as R µ ( t , t ) = s (inhom) µ ( t ) − s (hom) µ ( t ) h , (153)where s (hom) µ ( t ) must have been calculated independently. In thermal equilibrium the linear response andcorrelation function are related by the fluctuation dissipation relation, R µ ( t , t ) = 1 T µ ∂ t C µ ( t , t ) θ ( t − t ) . (154)Whether the time-evolving correlation and linear response satisfy this relation, whether the mode temperaturesare the same as the ones obtained from the energy characteristics of the modes and, finally, whether they alltake the same value, are issues that we will explore. .4 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 (cid:15) kin µ ( t ) = 12 m ˙ s µ ( t ) , (cid:15) pot µ ( t ) = 12 mω µ s µ ( t ) . (155)Clearly, neither the kinetic nor the potential mode energies are constant, but, in the steady state limit thesum of the two, that is to say the total mode energy, is (cid:15) kin µ ( t ) + (cid:15) pot µ ( t ) = (cid:15) kin µ ( t st ) + (cid:15) pot µ ( t st ) ≡ (cid:15) tot µ , (156)with t st the time at which the steady state is established and t > t st .As expected from Birkhoff’s theorem [67], the long-time averages, say taken after t st , should be constantsand one can expect them to be equal to half the total energy (cid:15) kin µ = 12 (cid:15) tot µ = 12 m (cid:20) ω µ s µ ( t st ) + 12 ˙ s µ ( t st ) (cid:21) ,(cid:15) pot µ = 12 (cid:15) tot µ = 12 mω µ (cid:20) s µ ( t st ) + 12 ˙ s µ ( t st ) ω µ (cid:21) . (157)(In practice, the average over a few periods is enough to obtain the constant value.) If one now associates atemperature to these values, arguing equipartition of quadratic degrees of freedom, one has T µ = 2 (cid:15) kin µ = 2 (cid:15) pot µ = (cid:15) tot µ . (158)The mode temperatures depend on the averages at the end of the transient, and the mode frequency Ω f µ thatitself depends on the asymptotic limit of the Lagrange multiplier z f and the eigenvalue λ µ .In the argument above we implicitly assumed that ω µ does not vanish. The case µ = N is tricky. If onenaively sets ω µ to zero from the outset (cid:15) pot N = ω µ (cid:104) s N (cid:105) apparently vanishes. The correct way of treating thelargest mode is to remember that the projection on the largest mode condenses and that (cid:104) s N (cid:105) is proportionalto N . This will ensure that lim N (cid:29) (cid:104) s N (cid:105) ∝ N , in such a way that lim N →∞ ω µ (cid:104) s N (cid:105) = 2 (cid:15) kin µ , similarly to whathappens in equilibrium, where (cid:104) s N (cid:105) = q in N and the Lagrange multiplier is such that ( z f − λ N ) q in N = T (cid:48) .We will see in the next Sections that, in some cases, the scenario described in this Section is actually realisedby the dynamics. Which are the quenches in which such a behaviour is observed will be determined by thecomplete solution of Newton’s equations with the methods that we will now describe. N In this section we address the calculation of equilibrium averages at finite N in order to provide suitableinitial 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, (cid:126)s (0) = ( s , . . . , s N ) and ˙ (cid:126)s (0) = ( ˙ s , . . . , ˙ s N ) , mimicking an initial thermal state at finite temperature, be it T (cid:48) > T c = J or T (cid:48) < T c = J , for a given realisation of the N × N interaction matrix. Averages over these initial states of theinteresting observables should then be computed. This method is computationally expensive as a large numberof initial state should be considered to get smooth and reliable results. Instead, the numerical method thatwe will explain in Sec. 5.8 is such that only the averages (cid:104) s µ (cid:105) eq and (cid:104) ˙ s µ (cid:105) eq are needed as input for the initialconditions. 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 µ , s µ } at temperature T (cid:48) , for agiven realization of disorder, is P GB ( { p µ , s µ } ) = 1 Z exp (cid:34) − T (cid:48) (cid:88) µ p µ m − T (cid:48) (cid:88) µ ( z ( N )eq − λ (0) µ ) s µ (cid:35) , (159)with Z the partition function. The statistical averages are computed as integrals over this measure. The integralsover p µ range from −∞ to ∞ . The quadratic averages of the velocities are thus simply given by (cid:104) ˙ s µ (cid:105) eq = T (cid:48) m ∀ µ, N , (160) ust as for the infinite N case, and the initial conditions will be (cid:104) ˙ s µ (0 + ) (cid:105) = (cid:104) ˙ s µ (cid:105) eq .As long as the equilibrium value of the Lagrange multiplier be strictly larger than the maximum eigenvalue, z ( N )eq > λ (0)max , (161)the weight of the coordinates s µ are well-defined independent Gaussian factors. We will see that the self-consistentsolution complies with this bound. Relying on the spherical constraint being imposed by the Lagrange multiplier,we extend the s µ integrals to ±∞ and (cid:104) s µ (cid:105) eq = T (cid:48) z ( N )eq − λ (0) µ ∀ µ, N . (162)The difference between the two equilibrium phases will be codified in the value of z ( N )eq , which can be obtainedas the solution of the spherical constraint equation N (cid:88) µ =1 (cid:104) s µ (cid:105) eq = N (cid:88) µ =1 T (cid:48) z ( N )eq − λ (0) µ = N . (163)We solved this equation numerically to determine z ( N )eq and we found that the solution turns out to be alwaysgreater than λ (0)max , for any value of the temperature and finite N . In Fig. 1 (a) we show z ( N )eq as a functionof temperature for three values of N and a single realisation of the random matrix in each case. At hightemperatures all the curves collapse (on the scale of the figure) on the paramagnetic curve z eq = T (cid:48) + J /T (cid:48) ,irrespective of the system size. At low temperatures (inset), z ( N )eq 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 (cid:104) s µ (0 + ) (cid:105) = (cid:104) s µ (cid:105) eq for the mode dynamics. . . . . . . .
05 0 . . z ( N ) e q T N=2565121024 N → ∞ PM z ( N ) e q − λ ( ) m a x T . .
01 0 0 .
002 0 . z ( N ) e q − λ ( ) m a x /NT = 0 . J . J . J . J . J Figure 1:
Equilibrium dynamics of the finite N system. (a) Equilibrium Lagrange multiplier at finite N . Solutionto Eq. (163) as a function of temperature using one particular realisation of disorder for each size. The non-monotonic N dependence of the plateau is of the order of magnitude of the variation with N of the largest eigenvalue. Inset: differencebetween the Lagrange multiplier and the maximum eigenvalue of the interaction matrix in the condensed region as a functionof temperature. The trend is now monotonic in N . (b) System size scaling of the Lagrange multiplier in the condensedphase. Difference between the Lagrange multiplier, as obtained from the solution to Eq. (163), and the maximum eigenvalueof the interaction matrix as a function of /N for different system sizes, using one particular realization of disorder for eachsize. The dashed lines are T (cid:48) / ( Nq in ) , with q in = 1 − T (cid:48) /J the value of the self-overlap in the N → ∞ limit.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 /N . The straight dashed lines have slope T (cid:48) /q in ,where q in = 1 − T (cid:48) /J is the N → ∞ value of the self-overlap. For temperatures sufficiently below the transition, he finite size data, obtained for one particular realisation of the random matrix J ij , follow the infinite sizeresults for all system sizes analysed. For temperatures close to the transition, there appear deviations for thesmallest system sizes (largest /N ). In conclusion, we find that for large system sizes or temperatures not tooclose to the transition, the solution to Eq. (163) behaves as z ( N )eq (cid:39) λ (0)max + T (cid:48) Nq in T (cid:48) < J . (164)Based on this, we define a finite size version of the equilibrium self-overlap q ( N )in ≡ T (cid:48) N ( z ( N )eq − λ (0) N ) , (165)which is finite if the highest mode is macroscopically populated. For N → ∞ , q in = 1 − T (cid:48) /J for T (cid:48) < J andzero for T (cid:48) > J . In Fig. 2 (a) we show q ( N )in as a function of temperature. We can observe the convergence ofthe finite size results towards the N → ∞ predictions as the system size is increased. . . . . . q ( N ) i n T N=2565121024 N → ∞ . . .
75 0 0 .
002 0 . z ( N ) e q − λ ( ) m a x /NT = J . J . J . J . J Figure 2:
Equilibrium dynamics of the finite N system. (a) Equilibrium q in at finite N . Self overlap as a functionof temperature for different system sizes, using one particular realization of disorder for each size. In this plot we check theleading finite order of q in and its dependence on T (cid:48) as − T (cid:48) /J far from the transition. Close to the transition there arefinite N corrections. (b) System size scaling of the Lagrange multiplier in the paramagnetic phase. Difference between theLagrange multiplier, as obtained from the solution of Eq. (163), and the largest eigenvalue of the interaction matrix as afunction of /N for different system sizes, using one particular realisation of disorder for each size. The non-vanishing valueat /N (cid:28) corresponds to z ( N →∞ )eq − J .Finally, we investigate the finite N corrections to z ( N )eq in the paramagnetic phase, T (cid:48) > J . We find that alinear scaling in /N also applies here, but the value of z ( N )eq − λ (0)max at N → ∞ does not vanish and it is givenby T (cid:48) + J /T (cid:48) − J . Then, in the paramagnetic phase we find z ( N )eq − λ (0)max (cid:39) z ( N →∞ )eq − J + s ( T (cid:48) ) N T (cid:48) > J (166)where s ( T (cid:48) ) = s is the slope of the dashed lines that we obtained from a fit and turns out to be independent oftemperature (all the dashed curves are parallel straight lines).Using the definition in Eq. (165), we can express the Lagrange multiplier as z ( N )eq = λ (0)max + T (cid:48) / ( Nq ( N )in ) , andwe can verify that the potential energy of the highest mode (note that we included the term proportional to theLagrange multiplier and we therefore compute (cid:15) pot N instead of e pot N ) (cid:15) pot N = 12 ( z ( N )eq − λ (0) N ) (cid:104) s N (cid:105) eq = T (cid:48) , (167)assumes the correct value in equilibrium, i.e., the one consistent with the equipartition theorem, if (cid:104) s N (cid:105) eq = q ( N )in N . (168) .6 Energy variation at the quench We here compute the finite- N equilibrium values of the Lagrange multiplier, kinetic and potential modeenergies using the finite size averages for (cid:104) s µ (cid:105) eq and (cid:104) ˙ s µ (cid:105) eq proposed in the previous Section, and we comparethem 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 ofthe observables at t = 0 + . We begin with the kinetic energy right before the quench, e ( N )kin (0 − ) = m N N (cid:88) µ =1 (cid:104) ˙ s µ (cid:105) eq = T (cid:48) . (169)It coincides with the infinite- N mean-field result.Next we analyze the pre-quench potential energy, e ( N )pot (0 − ) = − N N (cid:88) µ =1 λ (0) µ (cid:104) s µ (cid:105) eq = − T (cid:48) N (cid:32) N (cid:88) µ =1 λ (0) µ − z ( N )eq z ( N )eq − λ (0) µ + z ( N )eq N (cid:88) µ =1 z ( N )eq − λ (0) µ (cid:33) = T (cid:48) − z ( N )eq , (170)where we used Eq. (163). In the equilibrium condensed phase we can rely on z ( N )eq = λ (0)max + T (cid:48) / ( Nq ( N )in ) toobtain e ( N )pot (0 − ) = − λ (0)max T (cid:48) − T (cid:48) Nq ( N )in . (171)The N → ∞ result is e pot (0 − ) = − J + T (cid:48) / , consistent with Eq. (171), since lim N →∞ λ ( N )max = 2 J .In the paramagnetic phase q ( N )in (cid:28) and the third term in Eq. (171) induces important corrections. In thiscase, using Eq. (166) we can write e ( N )pot (0 − ) (cid:39) − J T (cid:48) − (cid:32) λ (0)max − J (cid:33) − s ( T (cid:48) )2 N (172)and one readily recovers the N → ∞ limit e pot (0 − ) = − J / (2 T (cid:48) ) . Now we will compute the values of the kinetic and potential energy, and the Lagrange multiplier after aninteraction quench λ (0) µ → λ µ = JJ λ (0) µ . (173)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 )kin (0 + ) = e ( N )kin (0 − ) = T (cid:48) . (174)For the potential energy it is enough to note that e ( N )pot (0 + ) = − N N (cid:88) µ =1 λ µ (cid:104) s µ (cid:105) eq = − N JJ N (cid:88) µ =1 λ (0) µ (cid:104) s µ (cid:105) eq = JJ e ( N )pot (0 − ) . (175)Using that z ( N ) (0 + ) = 2( e ( N )kin (0 + ) − e ( N )pot (0 + )) = T (cid:48) − J/J e ( N )pot (0 − ) it is now easy to find the initial valueof the Lagrange multiplier. When the initial conditions are taken from the condensed phase, q ( N )in = O (1) , andwe can write z ( N ) (0 + ) = λ max + T (cid:48) (cid:18) − JJ (cid:19) + JT (cid:48) NJ q ( N )in . (176)For initial states in the paramagnetic phase z ( N ) (0 + ) (cid:39) JJ T (cid:48) + T (cid:48) + ( λ max − J ) + JJ s ( T (cid:48) ) N . (177) .7 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 ofour non trivial problem. In equilibrium at time t = 0 − the initial frequencies of the modes are ω µ = [ z eq ( T (cid:48) , J ) − λ (0) µ ] /m . (178)In the asymptotic limit after the quench we identify the frequencies with ω µ = [ z f − λ µ ] /m (179)where we assumed that the Lagrange multiplier reached the constant z f = lim t →∞ z ( t ) .The analysis of the harmonic oscillator does not need any long-time assumption to set its spring constant, orfrequency, to a constant value. In our problem, the dynamics may approach the ones of independent harmonicoscillators with constant spring constants only asymptotically. During the transient evolution the mode energiesvary. In reality, we do not know the values they take at the end of the transient regime. We can make a roughapproximation in which we assume that z f is reached instantaneously after the quench, z (0 + ) = z f , so that wecan use (cid:104) p µ (0 + ) (cid:105) /m ≈ mω µ (cid:104) s µ (0 + ) (cid:105) ≈ T (cid:48) (180)instead of the unknown values at the end of the stationary regime.Under these assumptions the final mode temperatures are T fµ = T (cid:48) (cid:32) ω µ ω µ + 1 (cid:33) = T (cid:48) (cid:32) z f − λ µ z eq ( T (cid:48) , J ) − λ (0) µ + 1 (cid:33) , (181)see App. B for the details of the derivation. It is convenient to replace the post-quench eigenvalues λ µ by theirexpression in terms of the pre-quench ones and the quench parameter x = J/J , λ µ = xλ (0) µ . We can thendistinguish the four cases (I)-(IV) depending on the values of y = 2 T (cid:48) λ (0) N and x = λ N λ (0) N = JJ . (182)They are T fµ T (cid:48) = T (cid:48) + λ N T (cid:48) − λ µ T (cid:48) + ( λ (0) N ) T (cid:48) − λ (0) µ + 1 for y > x and y > (I) λ N − λ µ T (cid:48) + ( λ (0) N ) T (cid:48) − λ (0) µ + 1 for x > y and y > (II) λ N − λ µ λ (0) N − λ (0) µ + 1 for x > y and y < (III) T (cid:48) + λ N T (cid:48) − λ µ λ (0) N − λ (0) µ + 1 for y > x and y < (IV) (183)Several comments are in order. The expression for x > y and y < (sector III) is the same as the one that wederived from the analysis of the N → ∞ Schwinger-Dyson equations, see Eq. (123), simply T fµ = T (cid:48) ( x + 1) / T f = T (III)kin . The no-quench case x = 1 in realised in Sectors I and III and one rapidly checks that T fµ = T (cid:48) inboth cases. On the curve y = √ x the mode temperature do not take the same value. We will argue later thatthe approximation used in the Section yields a qualitatively erroneous result in this case. Continuity betweensectors I and IV on the one side, and II and III on the other, are ensured setting y = 1 that is to say T (cid:48) = λ (0) N .Finally, continuity across the dynamic transition at y = x or T (cid:48) dyn = λ N (184)is also verified. e would like to know which is the condition satisfied by z f under this approximation. In order to obtainsuch equation, we first note that the time-dependent spherical constraint imposes that N (cid:88) µ =1 (cid:104) s µ ( t ) (cid:105) = N . (185)In particular, this implies that N (cid:88) µ =1 (cid:104) s µ ( t ) (cid:105) = N . (186)Inserting the approximation in Eq. (180) in the time-averaged spherical constraint, we find an equation for z ( N ) f T (cid:48) m N (cid:88) µ =1 ( ω µ ) − (cid:34)(cid:18) ω µ ω µ (cid:19) + 1 (cid:35) = T (cid:48) N (cid:88) µ =1 (cid:34) z ( N )eq − λ (0) µ + 1 z ( N ) f − λ µ (cid:35) = N . (187)Since z ( N )eq is chosen in such a way that N (cid:88) µ =1 T (cid:48) z ( N )eq − λ (0) µ = N , (188)we find that the equation for z ( N ) f simplifies to N (cid:88) µ =1 T (cid:48) z ( N ) f − λ µ = N . (189)In other words, under this approximation, z ( N ) f is the equilibrium Lagrange multiplier for a system in equilibriumat temperature T (cid:48) with variance of the disorder distribution equal to J . In the N → ∞ limit z f = 2 J if T (cid:48) < J and z f = T (cid:48) + J /T (cid:48) if T (cid:48) > J .We will put these predictions to the test in Sec. 6 using the numerical solution to the finite N evolutionequations with the numerical method that we describe in the next Subsection. In various regions of the phasediagram these a priori approximate forms are in strikingly good agreement with the numerical data. In othersthey are not and we discuss why this is so. One possible approach to solve the dynamics of each mode starting from canonical equilibrium initial con-ditions is to take a large ensemble of initial configurations drawn from the Gibbs-Boltzmann distribution, nu-merically integrate the Newton equations Eq. (134) for each initial condition, and then calculate the observablesaveraging over the trajectories corresponding to the different initial states. Such approach is feasible but com-putationally very demanding. In this Section we describe a more convenient method to solve the dynamics foreach mode that uses heavily the tools developed to treat a paradigmatic problem in classical mechanics, the oneof parametric oscillators [77, 78, 79].In order to solve Eq. (134) we propose an amplitude-phase
Ansatz [77, 78, 79, 80] s µ ( t ) = A (cid:112) Ω µ ( t ) exp (cid:20) − i (cid:90) t dt (cid:48) Ω µ ( t (cid:48) ) (cid:21) . (190)Inserting this Ansatz in the µ th mode Newton equation, we obtain an equation for the mode and time dependentauxiliary function Ω µ ( t ) ,
12 ¨Ω µ ( t )Ω µ ( t ) − (cid:18) ˙Ω µ ( t )Ω µ ( t ) (cid:19) + Ω µ ( t ) = ω µ ( t ) , (191)where ω µ ( t ) ≡ ( z ( t ) − λ µ ) /m . The last equation has to be complemented by the initial values ω µ (0) and ˙Ω µ (0) .If we choose ˙Ω µ (0) = 0 , (192)we find that the projection of the spin configuration is s µ ( t ) = s µ (0) (cid:115) Ω µ (0)Ω µ ( t ) cos (cid:18)(cid:90) t dt (cid:48) Ω µ ( t (cid:48) ) (cid:19) + ˙ s µ (0) (cid:112) Ω µ ( t )Ω µ (0) sin (cid:18)(cid:90) t dt (cid:48) Ω µ ( t (cid:48) ) (cid:19) , (193) hich is reminiscent of the general solution of the harmonic oscillator problem, see Eq. (140), here with atime-dependent “frequency” Ω µ ( t ) .We still have to specify the initial condition for Ω µ (0) . A possible choice is Ω µ (0) = z (0) − λ µ (194)that enforces ¨Ω µ (0) = 0 [80]. However, this choice is consistent with real Ω µ ( t ) only if z (0) − λ µ ≥ , which isverified uniquely for J ≤ J , i.e., uniquely for energy injection. An initial condition ensuring real and positive Ω µ ( t ) for all µ for any quench is Ω µ (0) = λ N − λ µ . (195)We choose this initial condition for the numerical calculations.In order to solve for Ω µ ( t ) we consider the equal-times mode correlation function C µ ( t, t ) = (cid:104) s µ ( t ) (cid:105) = (cid:104) s µ (0) (cid:105) Ω µ (0)Ω µ ( t ) cos (cid:18)(cid:90) t dt (cid:48) Ω µ ( t (cid:48) ) (cid:19) + (cid:104) ˙ s µ (0) (cid:105) Ω µ (0)Ω µ ( t ) sin (cid:18)(cid:90) t dt (cid:48) Ω µ ( t (cid:48) ) (cid:19) + (cid:104) s µ (0) ˙ s µ (0) (cid:105) Ω µ ( t ) sin (cid:18)(cid:90) t dt (cid:48) Ω µ ( t (cid:48) ) (cid:19) cos (cid:18)(cid:90) t dt (cid:48) Ω µ ( t (cid:48) ) (cid:19) , (196)in terms of which we write the potential energy as e pot ( t ) = − N (cid:88) µ λ µ C µ ( t, t ) . (197)Replacing this equation in z ( t ) = 2 e f − e pot ( t ) , we find an expression of the Lagrange multiplier as a functionof the mode correlations at equal times z ( t ) = 2 e f + 2 N (cid:88) µ λ µ C µ ( t, t ) . (198)Finally, we note that the system conformed by Eqs. (191), (196) and (198) is closed and allows to find thetime evolution of the Lagrange multiplier and the auxiliary functions Ω µ ( t ) . This set of equations is amenableto numerical integration. Once we obtain Ω µ ( t ) , the most interesting observables can be calculated using thegeneral solution in the form given in Eq. (190). The advantage of this method is that we do not need to drawinitial states { s µ (0) , ˙ s µ (0) } but we only have to specify the initial averages (cid:104) s µ (0) (cid:105) and (cid:104) ˙ s µ (0) (cid:105) that we will taketo be the ones enforced by equilibrium at T (cid:48) , that is to say, the forms given in Eqs. (160) and (162). This Section summarises what we found numerically by solving the N → ∞ Schwinger-Dyson equations thatcouple the global correlation and linear response C and R (see Sec. 4), and the finite N ones acting on the modeprojections 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 wherebythe dynamics take place. We chose to start this Section with the summary of the dynamical phase diagram andthe behaviour of the quantities that determine it. Later, we give further details on the dynamics at constantenergy (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 statecharacterised 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 provethe lack of Gibbs-Boltzmann equilibrium but also because they lead to the evaluation of the mode temperaturesthat will play a role in the statistical description of the steady states. .1 The phase diagram The phase diagram is determined by the asymptotic behaviour of the zero frequency linear response orsusceptibility, χ st = ˆ R ( ω = 0) , and the asymptotic value of the Lagrange multiplier. We determine their valuesthrough the variation of the parameter J/J for fixed T (cid:48) /J . In the phase diagram presented in Fig. 5 and theensuing discussion we call y = T (cid:48) /J the vertical axis and x = J/J the horizontal one. The former determinesthe initial state and the latter the kind of quench performed with injection of energy for x < and extraction ofenergy for x > .We study separately parameters in four Sectors of the phase diagram, although the final results will allowus to distinguish three different phases. The Sectors are indicated with Roman numbers and the phases withdifferent colours or shades in Fig. 5. We also mark the line x = 1 (equilibrium dynamics) and the curve y = √ x with y > where special dynamics are found.We recall that dynamic phase transitions have been found in the quench dynamics of quantum isolatedsystems, see, e.g. [81, 82, 56, 57, 58, 83, 84, 61]. Here, and in [7], we see dynamic phase transitions arise in theNewtonian 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 (cid:48) ˆ R ( ω = 0) against J/T (cid:48) and we see the change in behaviour from ˆ R ( ω = 0) = 1 /T (cid:48) to ˆ R ( ω = 0) = 1 /J at x c ( y ) = y , thatis T (cid:48) = J . . . . . . χ s t · T ′ J/T ′ T ′ J = 0 . . . . Figure 3:
The zero frequency linear response , computed from the Schwinger-Dyson N → ∞ equations, for severalchoices of initial conditions given in the key, with both y < (condensed) and y > (paramagnetic) cases, together with theanalytic prediction in Eq. (115) plotted with a dashed line.The change in ˆ R ( ω = 0) is accompanied by a change in the asymptotic value of z as a function of the quenchparameter x = J/J . This fact is confirmed numerically in Fig. 4 where data for N → ∞ and N finite are shownin 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 polynomialfunction f ( x ) = a x + b x + c . We obtained very good results with a (cid:39) /y , b (cid:39) and c (cid:39) y (the precision of thefit is really very high in terms of reduced χ ). These results strongly suggest the following functional dependenceof z f on the parameters x and y , z f ( x, y ) = J x y + y for x ≤ y , x for x ≥ y . (199)To get a visual confirmation of this argument, in Fig. 4 (a) we plotted the functions f ( x ) = x /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 theprediction is very good. z f J/J N → + ∞ (a) T ′ /J = 0 . . . . . . . . z f J/J N = 1024 T /J = 0 . . . . . . Figure 4:
Estimated asymptotic value of z ( t ) as function of J/J , for different values of T (cid:48) /J , as indicated in the keys.(a) N → ∞ results. The dashed curved lines are functions of the form f ( x, y ) = y + x /y for x < y where x = J/J and y = T (cid:48) /J . We also show the diagonal x to let the reader see the crossover between the two regimes. (b) Finite size systemwith N = 1024 . The straight dashed line is Jλ max , the other curves are T (cid:48) + λ / (4 T (cid:48) ) , the finite size version of the infinite N fits. The analysis of the finite N data was done along the same lines, see Fig. 4 (b), with the difference thatthe data for x < y were fitted by T (cid:48) + λ / (4 T (cid:48) ) and the ones for x > y with Jλ max finding again very goodagreement. (We found an appreciable deviation in the fit for x > y had we used J instead of Jλ max . Regardingthe results for x < y we could have used T (cid:48) + J /T (cid:48) with a similar quality for the fit.) Remarkably, the functionaldependence proposed in Eq. (199) is the one predicted by the independent harmonic oscillators approximationin Eq. (189).By using the change in the dependence of ˆ R ( ω = 0) and z f on the quench parameter x = J/J as a criteriumto track the dynamical phase transition, we obtained the numerical estimates of the critical curve x c ( y ) . InFig. 5 we show the data for x c ( y ) (with error bars) for some values of the control parameter y . The data stronglysuggest that there is a linear relation between the critical value x c and the parameter y for any y ; in short, weconfirm x c ( y ) = y .Concerning the long-time behaviour of C ( t , t ) , it is useful to distinguish the cases y < and y > , that is,quenches that start from equilibrium in the condensed phase from quenches that start from equilibrium in theparamagnetic phase.We observe the following trends: • For x < x c ( y ) , C ( t , t ) tends to be stationary, though within the time scales of the numerics it has notreached this limit yet when y is too small. In most instances, C ( t , t ) oscillates around , exceptions beingthe critical quench and the case with both y > and x > where zero is approached asymptotically frombelow. The time average of C computed on intervals far from the initial transitory regime vanishes in allcases suggesting an effective q = 0 . • For x ≥ x c ( y ) , C ( t , t ) is rapidly stationary and one very clearly identifies the asymptotic constant q = lim t (cid:29) t t − t (cid:29) t C ( t , t ) . For y < it is different from zero while for y > it equals zero. The asymptotic q = lim t (cid:29) t C ( t , 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 severalchoices 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 thisSection we will support these claims with the detailed study of all relevant observables. . .
52 0 0 . . T ′ J JJ fa c hb g id e j χ st = T ′ q = 0 χ st = J q > χ st = J q = 0 I IIIIIIV
Figure 5:
The dynamic phase diagram.
The parameter x = J/J controls the energy injection/extraction, with x < corresponding to energy injection ( ∆ e > ), while x > to energy extraction ( ∆ e < ). The parameter y = T (cid:48) /J representsthe pre-quench equilibrium temperature. The dotted lines are the functions f ( x ) = √ x for x > and g ( x ) = 2 x/ (1 + x ) for x ≤ , respectively, discussed in the text. The three phases that are characterised by different behaviour of z f , χ st and thelong times limit of C ( t , t ) are highlighted with different colors. The red data points equipped with error bars indicate thenumerical estimate of x c for several values of y = T (cid:48) /J (see the main text for more details). The crosses indicate the casesshown in Fig. 6 with the corresponding labels. y = 0 . x = 0 . y = 0 . x = 0 . y = 0 . x = 0 . y = 1 . x = 0 . y = 1 . x = 1 . − . − . . . t − t (a) t = 0153045 t − t (b) t − t (c) t − t (d) t − t (e) y = 0 . x = 0 . y = 0 . x = 1 . y = 0 . x = 1 . y = 1 . x = 1 . y = 1 . x = 1 . − . − . . . t − t (f) t = 0153045 t − t (g) t − t (h) t − t (i) t − t (j) Figure 6:
The time-delayed correlation function C ( t , t ) from the numerical integration of the Schwinger-Dysonequations, for different values of the parameters y = T (cid:48) /J and x = J/J , as indicated above the plots. Crosses accompaniedby the labels (a)-(j) are marked at the corresponding location in the phase diagram in Fig. 5. In the first row x ≤ y while inthe second x > y . The first three panels in both rows correspond to y < and the last two to y > . The third panel in thefirst row is on the critical line x = y . We first checked that for J = J and m = m , that is to say ∆ e = 0 , and for all initial conditions, the systemhas stationary evolution and the total energy as well as other conserved quantities are indeed conserved. As the ynamics of generic Hamiltonian systems is hard to control numerically, we included this analysis to validateour algorithms. Moreover, this study allowed us to know which is the order of the numerical error incurred into.A time discretisation step δ = 0 . in the integration of the N → ∞ Schwinger-Dyson equations was sufficientto assure numerical convergence of our results. In the integration of the finite N problem we found a weakdependence on δ but the value δ = 0 . gave acceptable results.We studied the equilibrium dynamics for initial paramagnetic configurations ( T (cid:48) > J ) and condensed ones( T (cid:48) < J ). We present and discuss the results in two Subsections. As already said, they serve as benchmarks forthe more interesting quenching cases that we put forward later. In this Section we use parameters in the paramagnetic (PM) equilibrium phase. We fix J = m = m = 1 ,we equilibrate the initial conditions at T (cid:48) = 1 . and we evolve with parameters J = J = 1 and m = m = 1 (no quench). In Figs. 7 and 8 we show the numerical results for finite N and infinite N , respectively. The systemshould 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 (cid:48) = 1 . and all other parameters set to . The figure shows the time evolution of the correlation function (a), thefluctuation dissipation parametric plot (b), the response function, R ( t , t ) together with − (1 /T (cid:48) ) ∂ t C ( t , t ) (c),the deviation of the Lagrange multiplier from its equilibrium value (d), the potential and kinetic contributionsto 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 ( t , , and the ones between two different times, C ( t , t ) , areidentical, meaning that time translation invariance is satisfied. All the correlations relax to zero q = q = 0 .Moreover, the curves coincide approximately with the ones in Fig. 8 (a) which were obtained by integratingNewton equations for finite N .The fluctuation-dissipation relation [68] is satisfied with the temperature of the initial condition, that isthe same as the one of the final state. This fact can be proven in general for Newtonian evolution of initialconfigurations drawn from Gibbs-Boltzmann equilibrium [85]. Indeed in panel (b) we display the parametricplot χ ( t , t ) = (cid:82) t t dt (cid:48) R ( t , t (cid:48) ) vs C ( t , t ) for two waiting times t and, with a dashed line, the equilibrium result − /T (cid:48) C finding perfect agreement within our numerical accuracy. As a further confirmation of the validity ofthe fluctuation-dissipation theorem, we show the response function, R ( t , t ) , and − (1 /T (cid:48) ) ∂ t C ( t , t ) , bothplotted against t − t , for two different values of t 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 inverseof the theoretical Fourier transform of the stationary asymptotic response, given by Eq. (83). The theoreticalcurve is indicated as R st .The Lagrange multiplier and the potential and kinetic energies remain constant throughout the evolution ofthe system and equal to their predicted values, apart from small deviations due to the numerical errors introducedby the numerical integration scheme, see panels (d) and (f), and App. C. The plot showing z ( t ) proves that therelative error in this quantity is at most of order − . All these results are compatible with Gibbs-Boltzmannequilibrium in the paramagnetic phase. We do not show the time evolution of the off-diagonal correlation withthe initial configuration, C ( t, , since it is identically zero at all times.In panel (f) we show the Fourier transforms of the correlation and response functions, ˆ C ( ω, t ) and ˆ R ( ω, t ) respectively (the transform is performed on the variable τ = t − t with t fixed), for two different values of t .Note that we are showing only the real part of ˆ C ( ω, t ) , since we implicitly assume that C ( t + τ, t ) = C ( t − τ, t ) .The black solid lines represent the theoretical prediction for the real and imaginary parts of the Fourier transformof the response function in the stationary regime, ˆ R st ( ω ) , given by Eq. (83), the inverse Fourier transform ofwhich is plotted in (c). In panel (g), the ratio − Im ˆ R ( ω ) / ( ω ˆ C ( ω )) together with the prediction /T (cid:48) from FDTindicated by a horizontal dashed line are shown. Note the deviation from the flat result at the right edge of thefrequency spectrum. This is due to the fact that the ratio approaches zero over zero and the numerical errorincurred for those large frequencies is much amplified. At the left end of the spectrum, the more interesting lowfrequency regime, the oscillations are only present for the t = 0 curve.In Fig. 8 we show results obtained by solving the dynamics of each mode in a finite N system with the . − . . . . .
81 0 20 40 60 80(a) 00 . . . − . . . . . C ( t , t ) t − t t = 07 . . . χ C t = 022 . T ′ = 1 . t − t R ( t , t ) , t = 022 . − T ′ ∂ t C ( t , t ) , t = 022 . R st − · − − · − · − − . − . . . . z ( t ) − z e q t e ( t ) t − − . .
51 0 2 4 6 8 10(f) 0 . .
751 0 . . ω ˆ C , t = 0Re ˆ R Im ˆ R ˆ C , t = 7 . R Im ˆ R − I m ˆ R / ( ω ˆ C ) ω t = 07 . . /T ′ Figure 7:
Constant energy dynamics of the N → ∞ system in the PM phase . ∆ e = 0 is ensured by J = J and m = m . T (cid:48) = 1 . > T c = 1 and the initial condition is in the paramagnetic phase. (a) Dynamics of the correlation functionfor various choices of the waiting time given in the key. (b) Linear-response vs. correlation parametric plot for two valuesof the waiting time t . The dashed line shows the FDT with the initial temperature. (c) The response function, R ( t , t ) , − (1 /T (cid:48) ) ∂ t C ( t , t ) , and R st , the (numerical) inverse Fourier transform of the theoretical prediction given by Eq. (83) withparameters m = 1 , J = 1 and z f = z eq = T (cid:48) + J /T (cid:48) = 2 . , against t − t , for two values of t . (d) The difference betweenthe numerical Lagrange multiplier, z ( t ) , and the expected value at equilibrium, z eq . (e) Time evolution of the kinetic energydensity, e kin (red line), the potential energy density, e pot (blue line) and the total one, e f (green line). (f) Fourier transformsof the correlation (real part) and the response, for two values of t . The black solid lines represent the theoretical predictionsfor ˆ R ( ω ) , given by Eq. (83), the inverse Fourier transform of which is plotted in (c). We recall that we chose to use aconvention such that the imaginary part of ˆ R is negative. (g) The ratio − Im ˆ R ( ω ) / ( ω ˆ C ( ω )) together with /T (cid:48) indicated bya dashed horizontal line.method explained in Sec. 5.8. In panel (a) we see that the correlations with the initial condition quickly relaxto , as expected in the PM phase. They do with a weak size-dependence in the long time-delay tails. We onlyshow the correlation with the initial configuration since we have checked that the time-delayed one is stationary. . .
51 0 10 20 30(a) − . . − . . .
001 0 5 10(c) − − . . . µ = 1 − − . . . µ = N . . . . .
26 0 0 . C ( t , ) t N = 2565121024 N → ∞ e ( t ) te kin ( t ) e pot ( t ) e ( t ) z ( t ) − z ( N ) e q tδ = 0 . . . C µ ( t , ) / C µ ( , ) tN = 2565121024 C µ ( t , ) / C µ ( , ) tN = 2565121024 T µ ( t ) µ/Nt = 030 Figure 8:
Constant energy dynamics of a finite N = 1024 system in the PM phase . T (cid:48) = 1 . and all otherparameters are set to . (a) The correlation function between a configuration at time t and the initial condition for differentsystem sizes (grey curves) and the one for N → ∞ (red curve). (b) The kinetic, potential and total energies as a functionof time. (c) Dynamics of the Lagrange multiplier referred to the finite N equilibrium value z ( N )eq , for different discretisationsteps used in the numerical code. (d) and (e) Mode correlation for different system sizes for µ = 1 and µ = N , respectively.(f) Mode temperatures at different times. In (a), (b), (c)-(e) δ = 0 . , the discretisation step that we adopt in all furtherstudies.Also included in this panel is the same correlation function computed using the Schwinger-Dyson equation validin the N → ∞ limit. We see perfect agreement with the finite N results at short times and small deviations atlonger 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 ( N )eq , calculated as thesolution to Eq. (163) (a non-linear equation) that in the N → ∞ limit yields z ( N )eq → T (cid:48) + J /T (cid:48) . The very weak(oscillatory) deviation from z ( N )eq decreases with the size of the time-step used in the numerical solution of thedynamic equations (in the N → ∞ case we have a similar effect, see App. C).The mode-by-mode analysis of the finite N dynamics is performed in Fig. 8 (d) and (e). Two panels displaythe time-delay dependence of the correlation function of the first and last mode. In Fig. 8 (f) we display the modetemperatures T µ ( t ) at the initial time and after a long time evolution. The mode temperatures coincide with theexpected equilibrium value, except for the largest modes, where there is a very small deviation. These variationsrepresent small numerical errors due to the finite time-step discretisation used numerically, and are hard toimprove algorithmically unless by using a still smaller integration step. We have also checked (not shown) thatthe mode correlations C µ ( t , t ) and the mode response function R µ ( t , t ) satisfy the fluctuation-dissipationrelation with a temperature given by T (cid:48) for all modes. We now turn to the constant energy dynamics in the condensed, low temperature equilibrium phase. .
51 0 10 20 30(a) − . − . − . .
25 0 5 10(b) 00 .
001 0 5 10 15 20(c) − − . . . µ = 1 − − . . . µ = N − − . .
51 0 10 20 30(f) µ = N C ( t , ) tq in N → ∞ N = 2565121024 e ( t ) te kin ( t ) e pot ( t ) e ( t ) z ( t ) − z ( N ) e q tδ = 0 . . . C µ ( t , ) / C µ ( , ) t N=2565121024 C µ ( t , ) / C µ ( , ) t N=2565121024 C µ ( t , ) / C µ ( , ) t N=2565121024
Figure 9:
Constant energy dynamics in the condensed phase.
Results from the integration of the Newton equationsfor the individual modes at T (cid:48) = 0 . for a single disorder realisation. (a) Dynamics of the global correlation function with theinitial condition for different system sizes. (b) Evolution of the energetic contributions and the total energy. (c) Dynamicsof the Lagrange multiplier for different system sizes referred to the equilibrium z ( N )eq (that is slightly larger than λ max , seethe text). (d)-(f) Mode correlation function C µ ( t , t ) , normalised by their values at equal times, for µ = 1 , N − , N anddifferent system sizes indicated in the key. The mode correlations are stationary, so we only show results for t = 0 .In Fig. 9 we show the mode dynamics for initial conditions in equilibrium at T (cid:48) = 0 . . 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, whichshould be constant in equilibrium. In the numerical solution, the Lagrange multiplier exhibits oscillations aroundthe initial value. Their amplitude decreases consistently with the integration step δ , implying that for δ → werecover the expected constant behaviour. The two-time global correlation C ( t , t ) is stationary for all systemsizes (not shown), so we focus on the particular case with t = 0 . We can see from Fig. 9 (a) that, at variancewith the paramagnetic case, the dynamics of the correlation function C ( t , has a strong dependence on thesystem size. After a fast decay from the initial value, the correlation shows a plateau, the lifetime of whichincreases with system size, approaching asymptotically the value predicted by the N → ∞ treatment that isshown 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 C N ( t, for differentsystem sizes. We observe that its oscillation frequency decreases as we increase the system size. A similar finitesize effect is seen in the dynamics of the next-to-largest mode in panel (g). Since the largest modes dominate thelong-time dynamics, this effect causes the size dependence of the plateau lifetime. For N → ∞ the oscillationfrequency of the N -th mode goes to zero, allowing for the presence of an infinite plateau, see Fig. 11. The modeslying 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 beextracted from them. In Fig. 10 (c) we show the mode temperatures T µ ( t ) at two measuring times, as a functionof the mode index, what we will call temperature spectrum later. We observe deviations from the expected . . .
75 0 10 20 30 40(a) µ = 1 00 . . .
751 0 10 20(b) µ = N . .
505 0 0 . e µ ( t ) t KineticPotentialTotal T e µ ( t ) t KineticPotentialTotal δ = 0 . . . T T µ ( t ) µ/N t = 020 Figure 10:
Constant energy dynamics in the condensed phase. (a), (b) Mode energies for T (cid:48) = 0 . in equilibriumfor different mode indices in a system with N = 1024 . (c) Mode temperatures at different times. The total energy (red) isidentical to T (cid:48) (grey) and they both equal . . The potential and kinetic energies are equal to T (cid:48) / . .behaviour T tot µ = T (cid:48) ∀ µ 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 (cid:48) ,consistently with equilibrium dynamics. However, for the modes that are closer to µ = N the energies weaklyoscillate with an amplitude that decreases with the integration step and should vanish for δ → . This impliesthat in such limit, all mode temperatures should be equal to the equilibrium temperature, even for modes closeto 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 obtainedthrough the numerical integration of the Schwinger-Dyson equations for the two-time correlation and responsefunctions, and the same choice of parameters, that is T (cid:48) = 0 . and J = 1 . Stationarity is satisfied as well asthe FDT with the initial temperature, see (b) and (c). The main difference with the case in Fig. 7 is that thecorrelation functions, both with the initial condition and with the configuration at a waiting-time t , relax to anon-vanishing value (a). Within numerical accuracy we observe q = q (cid:39) − T (cid:48) /J = 0 . and this value as wellas the potential and kinetic energies (not shown) are consistent with equilibrium at T (cid:48) . Also in this case, wechecked that the response R ( t + τ, t ) , for t ≥ , coincides with the (numerical) inverse Fourier transform ofthe 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, C ( t, . In the case of equilibriumdynamics, C ( t, should be constant and equal to q in . As one can see, the value of C ( t, obtained by numericalintegration is not exactly a constant function, but it approaches a constant in the long time limit which differsfrom q = q in = 0 . by a very small amount. This deviation is only due to the approximations introduced by thenumerical 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 z eq = 2 J = 2 at T (cid:48) = 0 . , with oscillations amplitudeof order − for short times and decreasing with time.We next show the Fourier transforms of the correlation and response functions, for two different values of t in Fig. 7 (f). Again, note that we are showing only the real part of ˆ C ( ω, t ) , since we implicitly assume that C ( t + τ, t ) = C ( t − τ, t ) . The black solid lines represent the theoretical prediction for the real and imaginaryparts of the Fourier transform of the response function in the stationary regime, ˆ R st ( ω ) , given by Eq. (83), theinverse transform of which is plotted in (c). In panel (g) we display the ratio − Im ˆ R ( ω ) / ( ω ˆ C ( ω )) together withthe prediction /T (cid:48) from FDT indicated by a dashed horizontal line. In all presentations we find good agreementwith the validity of FDT with the proviso that in the plot in (g) the high frequency regime is contaminated bythe numerical error, and the low frequency regime by the fact that we can perform the Fourier transform on afinite time window only, and this causes the dependence on t shown in the plot. . . .
81 0 5 10 15 20 25 30(a) 00 . . . . . . . . . . . C ( t , t ) t − t t = 03 . . χ C t = 11 . T ′ = 0 . t − t R ( t , t ) , t = 015 − T ′ ∂ t C ( t , t ) , t = 015 R st − · − · − · − · − · − − · − − · − · − · − C ( t , ) − q t z ( t ) − z e q t − − . .
51 0 2 4 6 8 10(f) 1 . . . . ω ˆ C , t = 0Re ˆ R Im ˆ R ˆ C , t = 7 . R Im ˆ R − I m ˆ R / ( ω ˆ C ) ω t = 07 . . /T ′ Figure 11:
Constant energy dynamics of the N → ∞ system in the condensed phase. ∆ e = 0 , ensured by J = J . T (cid:48) = 0 . < T c = 1 . (a) Stationary dynamics of the two-time correlation function. The asymptotic limit is q = 0 . , shownas a dotted horizontal line. (b) Linear-response vs. correlation parametric plot. The dashed line shows the FDT with theinitial temperature. (c) The linear response function, R ( t , t ) , and the quantity − (1 /T (cid:48) ) ∂ t C ( t , t ) , both plotted against t − t , for two values of t . The curve indicated by R st is the (numerical) inverse Fourier transform of the theoreticalprediction given by Eq. (83) with parameters m = 1 , J = 1 and z f = z eq = 2 . (d) Difference between C ( t, and q . (e)The difference between the time-dependent Lagrange multiplier, z ( t ) , and the expected asymptotic value z eq = 2 J . (f) TheFourier transforms of the correlation and linear response functions, for two different values of t indicated in the key. Theblack solid lines represent the theoretical prediction for the real and imaginary parts of the Fourier transform of the responsefunction in the stationary regime, ˆ R st ( ω ) , given by Eq. (83), the inverse transform of which is plotted in (c). In panel (g),the ratio − Im ˆ R ( ω ) / ( ω ˆ C ( ω )) together with the prediction /T (cid:48) from FDT plotted with a dashed horizontal line. .3 Instantaneous quenches We shall now vary the initial temperature T (cid:48) and the coupling J of the Hamiltonian that drives the timeevolution to consider specific instantaneous quenching processes that inject or extract energy. The aim is toillustrate the analytical results of the previous Section and put them to the test. As we have already announcedthe structure of the phase diagram, we will consider representative quenches in the four sectors, labeled I ( y > and y > x ), II ( y > and y < x ), III ( y < and y < x ) and IV ( y < and y > x ).For each of the quenches performed we investigate if the system reaches a stationary state by checkingwhether: • The one-time quantities z ( t ) and e f tot = lim t → + ∞ e tot ( t ) approach constants that we measure and compareto the ones predicted in Sec. 4.4. • The kinetic and potential energy average over time to constant e kin and e pot that we also compare to theones predicted in Sec. 4.4. • The two-time correlation and linear response become stationary after a sufficiently long time, that is C ( t , t ) ∼ C st ( t − t ) and R ( t , t ) ∼ R st ( t − t ) for t large. • The long-time limits q = lim t → + ∞ C ( t, and q = lim t → + ∞ C ( t, are equal or differ.We also evaluate whether: • The susceptibility χ ( t , t ) = (cid:82) t t d t R ( t , t ) and C ( t , t ) satisfy the FDT relation χ ( t , t ) = (1 − C ( t , t )) /T f , with a single temperature T f . • The response function R ( t , t ) coincides with − /T f ∂ t C ( t , t ) , for different values of t , and with the(numerical) inverse Fourier transform of the analytic form given in Eq. (83).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 andpotential mode energies and we compare it to the global fluctuation dissipation ratio in the frequency domain ineach type of quench. This is motivated by the suggestion in [69, 70] (for quenches of isolated quantum integrablesystems) that these should be related. In this Subsection we summarise the dynamics of a system prepared in equilibrium in the paramagnetic state y > and quenched to a value of J such that y > x . This is what we called Sector I in the phase diagram. Twocases need to be further distinguished within this Sector. For x < energy is injected in the quench and thisproblem is treated in Figs. 12 and 13. For y > x > , instead, a small amount of energy is extracted from thesystem 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 > and x < , ˆ R ( ω = 0) = 1 /T (cid:48) and z f = T (cid:48) + J /T (cid:48) . The latter claim can be verified in Fig. 12 (d) for x > and N → ∞ , for example. Next we checkthat 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 responsedepend on the time delay only, see panels (a) and (c). The last question concerns the fluctuation dissipationrelation between linear response and correlation function. We used the parametric plot (b) to demonstrate thatthe evolution does not approach a state characterised by a single temperature but, instead, χ ( C ) is curved andnot even single valued. Having said this, the comparison of the time-dependent R ( t , t ) and − ∂ t C ( t , t ) inpanel (c) could have fooled us into the belief that the FDT holds with a single temperature. Indeed, the differencebetween the two functions is not visible in this scale. The Fourier analysis in (f) demonstrates that the frequencydependence of the real and imaginary parts of the linear response are highly non-trivial, though respecting thelimits in the frequency interval where it does not vanish derived in Sec. 4.4.2. The fluctuation-dissipation ratiois shown in (g) and we will come back to its analysis below, where we present the mode dependent results forfinite N .We next study the evolution under the same parameters in a single system with N = 1024 modes. Inthe first four panels in Fig. 13 we display the time dependence of the kinetic and potential energies, e kin µ ( t ) . − . . . (a) . . . . − . − . . . (b) − . . . ( ) C ( t , t ) t − t t = 015306090 χ C t = 03060 T f = 1 . t − t R ( t ,t ) ,t = 3060 − T f ∂ t C ( t ,t ) ,t = 3060 R st . . . .
625 0 20 40 60 80 100 (d) − . . . (e) z ( t ) t e ( t ) t − − (f) . . .
751 0 .
75 1 1 .
25 1 . (g) ω ˆ C , t = 30Re ˆ R Im ˆ R ˆ C , t = 60Re ˆ R Im ˆ R − I m ˆ R / ( ω ˆ C ) ω t = 30456075 T f Figure 12:
Sector I. Energy injection on a paramagnetic initial state. T (cid:48) = 1 . J , J = 0 . J and ∆ e = 0 . .(a) The time-delayed correlation function satisfies stationarity for t (cid:29) . The horizontal dashed line is at q = 0 . (b) Theparametric plot of the linear response function, χ ( t + τ, t ) , against the correlation, C ( t + τ, t ) , for two waiting times t .The black dashed line shows the FDT relation with T f (cid:39) . from Eq. (A.21). The numerical data exclude this possibilityfor the full range of C . In (c) R ( t , t ) and − /T f ∂ t C ( t , t ) as a function of t − t . In this scale we see no differencebetween the two functions; the presentation is misleading, since there is a difference, see (b). (d) Time evolution of theLagrange multiplier, z ( t ) , along with the constants T f + J /T f (cid:39) . (below) and T (cid:48) + J /T (cid:48) = 1 . (above) represented bydashed horizontal lines. (e) From top to bottom: the kinetic energy, the total energy and the potential energy densities invery good agreement with their expected values. (f) The Fourier transforms of the correlation and response functions withrespect to time delay. The black solid lines represent the real and imaginary parts of ˆ R st ( ω ) , given by Eq. (83). In panel (g),the ratio − Im ˆ R ( ω ) / ( ω ˆ C ( ω )) together with /T f (cid:39) . that is off the data but not very far away from the constant at thelimit ω → ω + c (recall that the downward trend close ω + c is due to numerical inaccuracy). We note that /T kin (cid:39) . is stillfarther away from the high frequency limit. Results for the same parameters and finite N are shown in Fig. 13.and (cid:15) pot µ ( t ) . These functions oscillate as a function of time with amplitude that slowly decreases in the cases µ = 1 , µ = N − , µ = N , while it slowly increases for µ = 2 , in this time window. The total energy of each modeshown with a solid red line displays a small downward drift towards, one may expect, a constant value. TheLagrange multiplier shows a residual time variation around a value that is slightly above the N → ∞ prediction . . µ = 1 0 . . µ = N/ µ = N − µ = N . . e µ ( t ) t KineticPotentialTotal e µ ( t ) t KineticPotentialTotal e µ ( t ) t KineticPotentialTotal e µ ( t ) t KineticPotentialTotal z ( t ) t z f . . . .
98 0 .
99 1 T µ µ/NT pot µ T kin µ T µ ( t = 40)Analytical T µ µ/N . . .
751 0 .
75 1 1 .
25 1 . β ( ω ) ω [ T µ ( ω µ )] − , N = 1024 − Im[ ˆ R ( ω )] / ( ω ˆ C ( ω )), N → ∞ Analytical
Figure 13:
Sector I. Energy injection on a paramagnetic initial state.
A system with N = 1024 and parameters T (cid:48) = 1 . J and J = 0 . J , leading to ∆ e = 0 . . (a)-(d) The time-dependent energy of four selected modes. (e) TheLagrange multiplier compared to the N → ∞ prediction z f = T (cid:48) + λ / (4 T (cid:48) ) . There is a small deviation most likelydue to numerical error. (f) Mode temperatures extracted from the use of equipartition, T kin , pot µ = 2 e kin , pot µ . Inset: Detailof the behaviour of the largest modes. (g) Comparison between the modes inverse temperatures and the inverse effectivetemperature from the fluctuation dissipation relation of the N → ∞ equations in the frequency domain (both measured inunits of J ). In (f) and (g) the yellow curves are given by the approximate prediction for T µ in Eq. (183). The frequencyinterval where Im ˆ R ( ω ) (cid:54) = 0 is [ ω − , ω + ] with ω − (cid:39) . and ω + (cid:39) . in this case.but is very close to its finite N modified value T (cid:48) + λ / (4 T (cid:48) ) (shown with a solid horizontal line). We thendetermine the mode temperatures from the equipartition of the kinetic and potential energies averaged over asufficiently 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, seeEq. (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 thefluctuation-dissipation inverse temperature. The agreement is very good (except at the edge of the spectrumwhere both numerator and denominator vanish and it is very hard to control the limiting behaviour even inthe equilibrium case, see Fig. 7). The (yellow) continuous line is the approximate theoretical prediction inEq. (183) that captures the numerical behaviour rather accurately. We recall that it was derived assuming thatthe Lagrange multiplier takes its asymptotic constant value immediately after the quench, at time t = 0 + , and hat 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 > . We repeated the analysis of the N → ∞ equations for parameters in this regime, choosing, for example, T (cid:48) = 1 . J and J = 1 . J . We do not present the data here as there are no fundamental variations with respectto what we have already discussed for the energy injection case. Having said so, the mode temperatures dopresent an interesting difference that we discuss with the support of Fig. 14. Also in this case, the temperaturesof the modes are approximately the same for the low lying modes while a mode dependence appears close tothe edge of the spectrum. However, the temperatures of the largest modes are, in this case, lower than thetemperatures of the lower modes, see panel (a) and its inset. We ascribe this feature to the fact that the quenchextracts energy from the system. In panel (b) we confront the mode inverse temperatures to the ones extractedfrom the fluctuation dissipation ratio in the frequency domain and, once again, the agreement between numericalcurves for N → ∞ and finite N data is very good. Moreover, the data are also in good agreement with theoutcome of the assumption z (0 + ) = z f that leads to Eq. (183), shown with a yellow solid line. . . . . . . . . . .
97 0 .
98 0 .
99 1 T µ µ/NT pot µ T kin µ T µ ( t = 40)Analytical T µ µ/N . .
25 0 . . . β ( ω ) ω [ T µ ( ω µ )] − , N = 1024 − Im[ ˆ R ( ω )] / ( ω ˆ C ( ω )), N → ∞ Analytical
Figure 14:
Sector I. Energy extraction from a paramagnetic initial state. T (cid:48) = 1 . J and J = 1 . J such that y > x > . (a) Non-equilibrium temperature profile with the temperatures of the largest modes being smaller than the rest forenergy extraction. Inset: Zoom over the behaviour of the largest modes. (b) Comparison of the inverse mode temperatureswith the frequency-dependent effective inverse temperature of the global fluctuation dissipation relation. Also shown with ayellow solid line is the approximate analytic prediction in Eq. (183). The frequency interval where Im ˆ R ( ω ) (cid:54) = 0 is [ ω − , ω + ] with ω − = 0 and ω + (cid:39) . in this case. For these parameter values the Lagrange multiplier approaches z f = 2 J . We confirmed this claim with thestudy 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 inSector II: y > √ x , y = √ x and y < √ x , all with y > and y < x . y > √ x An example of the decay of the two-time global correlation function can be seen in panel (j) in Fig. 6. It isstationary and it rapidly approaches zero with progressively decaying oscillations around this value. The linearresponse and correlation function are not related by an FDT with a single temperature (not shown) and in thisrespect there are no important differences regarding what we have just shown for energy extraction processes inSector I. For these reasons we chose not to show more data for these parameters. = √ x A prediction from the asymptotic analysis of the Schwinger-Dyson equations, the fact that on the curve y = √ x and y ≥ FDT is satisfied with T f = J , is made explicit in Fig. 18 where the parametric plot T f χ ( t − t ) against C ( t − t ) is constructed for four pairs of ( x, y = √ x ) given in the key. The dotted line is theFDT prediction with T f = J . The agreement between numerics and analytics is very good. Additional supporton the fact that the dynamics after the quench occur as in equilibrium at T f is given in panel (b) in the samefigure where quenched and equilibrium correlations are indistinguishable. The latter are obtained by drawingthe initial conditions drawn from equilibrium at T (cid:48) = 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). . . . . − .
25 0 0 .
25 0 . .
75 1 (a) χ · T f C T ′ J = 1 . . . − . . . . . .
01 0 . (b) C ( t , t ) t − t t = 1530 T ′ = J = J = 3 . Figure 15:
Sector II. Energy extraction from a paramagnetic initial state for parameters lying on the special curve y = √ x of the phase diagram. (a) Check of FDT at T f = J for four sets of parameters on this curve. (b) Comparison betweenthe time-delayed correlation after a quench with T (cid:48) = 1 . and J = 3 . , and the equilibrium (no quench) correlation at T f = J = 3 . . The agreement is perfect. y < √ x We chose to show data for T (cid:48) = 1 . J and J = 2 J , parameters that lie in the region y < √ x of the phasediagram 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 expectedasymptotic values with an oscillatory behaviour and smoothly decaying amplitude (not shown). The two-timecorrelation is stationary and shows no ageing. The decorrelation from the initial configuration and from a laterconfiguration at time t behave very similarly. The time-delayed C ( t , t ) shows a rapid decay towards a valueclose to . around which it oscillates once and next decays towards zero with damped oscillations (a). Thelinear response function shows a similar effect in the sense of having a fast variation at short time differences anda slower one later (c). The value of R obtained from the numerical integration agrees very well with R st fromEq. (83). The most interesting results concern the comparison between the linear response and the correlationfunction in the parametric plot in panel (b). The very short time delay behaviour, for C close to and χ closeto seems to be described by the slope dictated by T f , the value of the temperature deduced from an ageinglike asymptotic scenario. However, the curves deviate from this straight line for smaller C and larger χ . Whenthe correlation reaches a value close to . corresponding to its first oscillation, the parametric plot approaches aflat form that ensures the limit χ st = 1 /J . This second behaviour is reminiscent of what happens in the ageingrelaxation of the same model [44].The asymptotic analysis of the N → ∞ equations allow for an ageing solution with a diverging effectivetemperature 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 numericalsolution of the complete set of Schwinger-Dyson equations demonstrates a behaviour with some features that are . . . .
81 0 10 20 30 40 50 (a) . . . . . . . (b) − . . . . ( ) C ( t , t ) t − t t = 0306090120 χ C t = 03060 T f = 1 . t − t R ( t , t ) , t = 3060 − T f ∂ t C ( t , t ) , t = 3060 R st . .
97 0 .
98 0 .
99 1 T µ µ/NT pot µ T kin µ T µ ( t = 40)Analytical T µ µ/N . . . . . . . . β ( ω ) ω [ T µ ( ω µ )] − , N = 10241 /T f Analytical
Figure 16:
Sector II. Energy extraction from a paramagnetic initial state ( y < √ x ) T (cid:48) = 1 . > T c , and post-quench coupling J = 2 leading to ∆ e = − . . (a) The correlation function, C ( t , t ) against t − t for different t . C ( t , t ) vanishes as t − t → + ∞ . (b) The parametric plot χ ( τ, t ) against C ( t + τ, t ) , for two different values of the waiting time t . The black dashed line is the FDT at T f (cid:39) . J , see Eq. (A.22). (c) R and − /T f ∂ t C as a function of t − t ,with the same T f as in (b). R st is reproduced from Eq. (83) with m = 1 , J = 2 and z f = 4 . Again, in this representationthe difference with − /T f ∂ t C is not visible, see panel (b) for a better understanding. (d) Almost all modes in the bulk ofthe spectrum have the same temperature, and it is very close to T f shown with an horizontal dashed line. Inset: Detail ofthe largest modes. (e) The inverse mode temperatures (red data) and the outcome of the harmonic oscillator approximation(yellow curve) in Eq. (183).similar to the approximate asymptotic solution, but no signature of ageing. Indeed, the parametric plot couldbe interpreted as formed by two pieces, one in which the form is non-trivial close to C = 1 and another one thatis, on average, flat, separated by the correlation value at which the first oscillation occurs. A flat piece in theparametric plot means that the integrated response, and all other functions such as the potential energy, do notget contributions from these time scales and it corresponds to an infinite [44] effective temperature [51, 52]. Inpractice, the parametric plot is not locally flat for small values of C but, as we can see from the mode analysisin Fig. 16 (e), the temperatures of the corresponding modes do take very large values. . . .
81 0 5 10 15 20 25 30 C ( t , t ) t − t (a)0 . . . t (c) 00 . . . . . . . χ C (b) − . . . . . t − t (d) t = 03 . . . C ( t, C ( t, q = 0 . t = 07 . . T f = 0 . R ( t , t ) , t = 3 . . − T f ∂ t C ( t , t ) , t = 3 . . R st . . . . − . − . − . . z ( t ) t e ( t ) t Figure 17:
Sector III. Energy injection, shallow quench from condensed to condensed. T (cid:48) = 0 . < T c and J = 0 . . The small energy injection, ∆ e = 0 . , is not sufficient to drive the system out of a condensed state. (a) Dynamicsof the correlation function. The horizontal dashed line is at q = 1 − T f /J (cid:39) . , with T f from Eq. (A.19). (b) χ ( t , t ) against C ( t , t ) for fixed t and using t − t as a parameter. The black dashed line shows the FDT with T f = 0 . . In(d) R , R st from Eq. (83), and − /T f ∂ t C as a function of t − t . (e) Time evolution of the Lagrange multiplier, z ( t ) , and z f = 2 J = 1 . with a dashed line. (f) From top to bottom: the kinetic energy (in good agreement with e f kin = T f / ), thetotal energy (constant in time) and the potential energy (in good agreement with e f pot = − J / (2 T f )(1 − q ) ). In Fig. 17 we start discussing the behaviour of the N → ∞ model in Sector III. We use T (cid:48) = 0 . and J = 0 . ,a quench that injects a small amount of energy from the system. Panel (a) proves that the two time correlationapproaches a non-vanishing value asymptotically. The horizontal dashed line is at q = 1 − T f /J (cid:39) . , thetheoretical value derived from the assumption of equilibration à la Gibbs-Boltzmann. Further information aboutthe decay of the correlation functions is given in (c) where we show C ( t, and the off-diagonal correlation withthe initial configuration C ( t, against time. C starts at q in = 0 . and decreases monotonically. C ( t, quicklydecays from with superimposed oscillations. Both curves should reach q = q (cid:54) = q asymptotically and thedata are compatible with this claim.Panel (a) in Fig. 18 shows the value of the asymptotic potential energy as a function of J/J for the samethree values of T (cid:48) /J . The agreement between e f pot and T (cid:48) / T (cid:48) J (4 J ) − J , the parameter dependence foundfrom the energy conservation and the asymptotic value of z is extremely good.For quenches with y = T (cid:48) /J < and x > x c ( y ) = y , the asymptotic values of C ( t , t ) and C ( t ,
0) = ( t , , q and q , respectively, do not vanish. However, it is not always easy to extract their functionaldependence on x and y , especially for parameters that get away from the shallow quenches x (cid:39) . Figure 18 (b)shows q and q against J/J for three values of T (cid:48) /J , together with the naive single temperature prediction q = 1 − T f /J with T f = 2 T kin = T (cid:48) (1 + J/J ) drawn with black solid lines. The agreement is good close to x = 1 though deviations are clear close to the critical line x c ( y ) and for large energy extraction deep inside thisparameter sector. Note that the prediction of equilibration à la Gibbs-Boltzmann is such that q (cid:54) = 0 at x c ( 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 q . One clearly finds that q > q for x < (injection) and q < q for x > (extraction). We will come back tothe behaviour of the correlation function below. − . − − . − − .
50 0 0 . . . (a) e f p o t J/J T ′ /J = 0 . . . . . . .
81 0 0 . . (b) J/J qq Figure 18:
Sector III, shallow and deep quenches from condensed to condensed.
The three datasets correspondto T (cid:48) /J = 0 . (red), . (green), . (blue). (a) The final potential energy as a function of J/J . The dotted lines are T (cid:48) / T (cid:48) J (4 J ) − J . (b) The estimated asymptotic value of the correlation function, q = lim ( t − t ) → + ∞ ,t (cid:29) C ( t , t ) (simple points), and q = lim t → + ∞ C ( t, (circles), against J/J for the three values of T (cid:48) /J < (increasing from top tobottom). Data are equipped with error bars. The solid lines are the equilibrium predictions for q . Notice the deviationsclose to x c = y and for x far from . The long time limits of C ( t , t ) and C ( t , are ordered according to q > q for J < J and q < q for J > J .Panel (b) in Fig. 17 shows the parametric χ ( C ) curve where one sees that the t = 0 (red) data are very differentfrom the ones for long t . The data for long t suggest that FDT has established at temperature T f . Another testof FDT, now in the time domain, is shown in (d) with the comparison between the linear response and the timederivative of the correlation. The other panels show the asymptotic values of the Lagrange multiplier (e) andkinetic 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 withthe asymptotic contributions to the total energy. Nevertheless, we will revisit this claim below when studyingdeeper quenches in the same sector.The companion curves for finite N are in Fig. 19. First of all, panels (a)-(d) display the time dependence ofthe µ = 1 , N/ , N − , N mode energies in a system with N = 1024 . While the modes µ = 1 , N/ show theusual oscillatory behaviour of a harmonic oscillator, the largest modes µ = N − and µ = N are clearly outof equilibrium. The Lagrange multiplier is approximately constant and equal to the largest eigenvalue, withinnumerical accuracy. The spectrum of mode temperatures is plotted in (f) with a zoom over the largest modesin its inset. The profile is an equilibrium one, with T µ being independent of µ , apart from the deviations closeto 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 resultscoincide (except at high frequencies where the result is biased by the numerical limitations in the computationof the Fourier transform). Higher modes (low frequency) have temperatures slightly below the temperature T f . . µ = 1 0 . . µ = N/ − . . . .
75 0 10 20 30(c) µ = N − − − . . .
52 0 10 20 30(d) µ = N . .
75 0 10 20 30(e) e µ ( t ) t KineticPotentialTotal e µ ( t ) t KineticPotentialTotal e µ ( t ) t KineticPotentialTotal e µ ( t ) t KineticPotentialTotal z ( t ) t λ max . . . . . . . . .
97 0 .
98 0 .
99 1 T µ µ/NT pot µ T kin µ T µ ( t = 40) T f T µ µ/N . . . . . . β ( ω ) ω [ T µ ( ω µ )] − , N = 1024Analytical Figure 19:
Sector III, condensed initial conditions and small energy change. T (cid:48) = 0 . and J = 0 . . (a)-(d) Modeenergies in a system with N = 1024 . The largest modes are neatly out of equilibrium. (e) The Lagrange multiplier. (f) Modetemperatures, with a zoom over the largest modes in the inset. Close to equilibrium like profile apart from the deviationsclose to the edge of the spectrum. (g) Comparison of the inverse mode temperatures with the ones of independent harmonicoscillators.while lower modes (high frequencies) have temperatures slightly above T f , that is shown with an horizontalyellow line. This fact is another warning concerning the claim of complete equilibration to a Gibbs-Boltzmannmeasure.Indeed, although the data for the shallow quench that we have just described might have suggested equilibra-tion to a Gibbs-Boltzmann probability distribution, the detailed comparison of the full time-delay dependence ofthe correlation function after a quench and in equilibrium (no quench) at parameters J and T f that are the onesthat the equilibrium measure would have, prove that such a steady state is not reached by the dynamics. Thisstatement is proven in Fig. 20 where we display the self-correlations stemming from the two procedures for threechoices 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 identicalresults. 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 < . We givehere further evidence for this fact. In Fig. 21 we show results for T (cid:48) = 0 . < T c and J = 0 . , that is to say, . . . . .
01 0 . (a) C ( t , t ) t − t t = 3060eq. at T f . . . .
01 0 . (b) t − t t = 3060eq. at T f . . . . . .
01 0 . (c) t − t t = 3060eq. at T f Figure 20:
Lack of Gibbs-Boltzmann equilibration in Sector III.
The time-delay correlation function after a quenchfrom J to J in Sector III (red and green curves) and in equilibrium at parameters J and T f , the single temperature thatwould correspond to Gibbs-Boltzmann equilibrium (black curves). The three panels are for J = 0 . (quench to the criticalline), J = 0 . (quench close to x = 1 ) and J = 2 (deep quench). The dotted horizontal line is q = 1 − T f /J .point c on the phase diagram in Fig. 5. This quench injects a relatively small amount of energy into the system, ∆ e (cid:39) . and takes the parameters to be on the critical line y = x .The self correlation is shown in panel (a) of Fig. 21. Stationarity is clear for short time delays t − t , for t > and there is some remanent waiting time dependence at these short t . This double behaviour is reminiscent ofwhat is seen in the relaxational dynamics where a sharp separation of time-scales exists. The data suggest that q = lim t →∞ C ( t, is different from q = lim t − t →∞ lim t →∞ C ( t , t ) , although it seems hard to determinethese values from the numerics with good precision. In the event of a two step relaxation of C ( t , t ) with anapproach to a plateau at q and a further decay from it to zero, the plateau should be at q = 1 − T f /J = 0 . ,shown with a dashed horizontal line. The data still lie below this value and we infer that they might not convergeto 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, and theoff-diagonal correlation with the initial configuration C ( t, against time. C starts at q in = 0 . and decreasesmonotonically. C ( t, quickly decays from with superimposed oscillations. Both curves seem to join and slowlyand monotonically decay to zero.The parametric plot of the linear response function, χ ( t , t ) , against the correlation, C ( t , t ) , for fixed t and using t − t as a parameter, for three different values of the waiting time t , is shown in panel (b). The χ ( C ) curve for t = 0 does not have any special form. The curves for late t are close to the straight line /T f fortime-delays that correspond to the first oscillations of the correlation and linear response, they then oscillate, andfor longer time-delays the parametric construction becomes flat, with χ approaching, for C → , the expectedstatic susceptibility /J = 2 . Again, this second flat regime is reminiscent of the one found for the relaxationstochastic 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 . obtained from Eq. (A.19), is indicated as a straight dashed line. The presentation in panel (d) confirms theagreement between linear response and time variation of the time-delayed correlation function dictated by theFDT at temperature T f for C ≥ . , say. However, it clearly breaks for smaller values of C . As for the otherquenches, we checked that the response R ( t + τ, t ) coincides with the one derived by anti-transforming thetheoretical Fourier transform of the stationary asymptotic response, given by Eq. (83), R st .In (e) we provide a close look at the time evolution of z ( t ) . As one can see, z ( t ) approaches the predictedvalue z f = 2 J . 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 coincidewithin numerical accuracy with the ones predicted in Sec. 4.4.3. These values could be mistaken for e f kin = T f / and e f pot = − J (1 − q ) / (2 T f ) , with q = 1 − T f /J = 0 . , the Gibbs-Boltzmann equilibrium predictions, though,the system is not in equilibrium. .
51 0 20 40 60 80 100 120 140 C ( t , t ) t − t (a)0 . . . . . t (c) 012 0 0 . . . . χ C (b)00 . . . t − t (d) t = 07 . . . . . C ( t, C ( t, q = 0 . t = 037 . . T f = 0 . R ( t , t ) , t = 37 . . − T f ∂ t C ( t , t ) , t = 37 . . R st . . . − . − . . z ( t ) t e ( t ) t Figure 21:
Energy injection from condensed to the critical line y = x . T (cid:48) = 0 . < T c , J = 0 . and ∆ e = 0 . .(a) The correlation function and q = 1 − T f /J = 0 . (horizontal dotted line). (b) χ ( t , t ) against C ( t , t ) , for fixed t and using t − t as a parameter. The black dashed line is the FDT with T f = 0 . . (d) The curve indicated by R st is the(numerical) inverse Fourier transform of the theoretical prediction given by Eq. (83). (e) Time evolution of the Lagrangemultiplier, z ( t ) , and J = 1 with a dashed horizontal line. (f) From top to bottom: the kinetic energy (in good agreementwith e f kin = T f / (cid:39) . ), the total energy (constant in time with value e f (cid:39) − . ) and the potential energy (in apparentagreement with e f pot = − J T f (1 − q ) (cid:39) − . , though see the text for a revision of this claim). In Fig. 22 we show results for T (cid:48) = 0 . < T c and J = 0 . . This quench injects a large amount of energy intothe system, ∆ e = 0 . , which is sufficient to take it out of the initial condensed state. Had the system reachedan equilibrium paramagnetic state asymptotically after the quench its temperature would be T f (cid:39) . , so that T f /J (cid:39) . , 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 theexpected asymptotic value lim t − t →∞ lim t →∞ C ( t , t ) = 0 for all values of t . C ( t , t ) satisfies stationarity forsufficiently late t , as one can see from the plot where the curves for t coincide, within numerical accuracy. Theparametric plot of the susceptibility χ against the correlation C , shown in (b), is rather complex and very far fromlinear. The FDT relating χ to C in a stationary regime with putative target temperature T f (cid:39) . is indicatedas a straight dashed line. This behaviour is confirmed by the time evolution of the stationary response function, R ( t , t ) , see panel (d). R ( t , t ) does not coincide with − T f ∂ t C ( t , t ) , with T f the target temperature. In . . . (a) − . . ( ) − . . . (b) − . . . (d) C ( t , t ) t − t t = 0154575 t C ( t, C ( t, χ C t = 0607590 T f = 0 . t − t R ( t , t ) , t = 7590 − T f ∂ t C ( t , t ) , t = 7590 R st . . . . . . (e) − . − . . . . (f) z ( t ) t e ( t ) t − − . .
75 1(h) ω ˆ C , t = 75Re ˆ R Im ˆ R ˆ C , t = 105Re ˆ R Im ˆ R − I m ˆ R / ( ω ˆ C ) ω t = 45751051 /T f Figure 22:
Sector IV. Large energy injection on a condensed state. T (cid:48) = 0 . < T c , J = 0 . and ∆ e = 0 . . (a)Dynamics of the correlation function. The horizontal line is at q = 0 . (b) χ ( τ + t , t ) against C ( τ + t , t ) for four waitingtimes t specified in the key. The black dotted line is the FDT at T f (cid:39) . . (c) Comparison between C ( t, and C ( t, . C starts at q in = 0 . and oscillates around , even though the amplitude of oscillations is slowly decreasing with time. (d) R , R st and − /T f ∂ t C as functions of t − t . (e) The Lagrange multiplier z ( t ) , T f + J /T f (dotted line) and z f = T (cid:48) + J /T (cid:48) (dashed line). (f) From top to bottom: kinetic, total (constant), and potential energy densities. (g) Fourier transforms ofthe correlation and response functions for two t indicated in the key. The black solid lines are the theoretical predictionsfor the real and imaginary part of the Fourier transform of the linear response function, given by Eq. (83), with parameters m = 1 , J = 0 . and z f = 0 . . (h) The ratio − Im ˆ R ( ω ) / ( ω ˆ C ( ω )) where the various curves correspond to different waitingtimes t . The dashed line is at /T f (cid:39) . and /T kin = 2 . would be even below it.panel (c) we also show the time evolution of the off-diagonal correlation with the initial configuration, C ( t, . Incontrast to the previous cases, C ( t, is oscillating in time around the value q = 0 . Panel (e) demonstrates that z ( t ) is far from z f = T f + J /T f relative to an equilibrium paramagnetic state but oscillates around T (cid:48) + J /T (cid:48) , onsistently with our claims so far. From panel (f) we observe that the kinetic and the potential energies relaxto the predicted asymptotic values specified in Sec. 4.4.3. The Fourier transforms of the correlation and responsefunctions for two different values of t , as indicated in the key, are shown in panel (g). The black solid linesrepresent the theoretical prediction for the real and imaginary part of the Fourier transform of the responsefunction in the stationary regime, ˆ R st ( ω ) , given by Eq. (83), with parameters m = 1 , J = 0 . and z f = 0 . .They are in excellent agreement with the numerical results. In panel (h) we show the ratio − Im ˆ R ( ω ) / ( ω ˆ C ( ω )) .The dashed line is − /T f and is clearly off the data. The different curves were computed for various waitingtimes t given in the key. In Fig. 23 we compare the FDR (averaged over different waiting times to get rid ofthe undesired oscillations) to the mode temperatures of the finite N system. The correspondence between thetwo 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 farfrom the numerical results though it captures the qualitative features. . . . β ( ω ) ω β f [ T µ ( ω µ )] − , N = 1024 − Im[ ˆ R ( ω )] / ( ω ˆ C ( ω )), N → ∞ Analytical
Figure 23:
Sector IV. Large energy injection on a condensed state.
Same parameters as in Fig. 22. Comparisonbetween the mode inverse temperature for a finite size system (red), the inverse temperature for the fluctuation dissipationratio in the infinite size limit (black curve) and the analytical expression for the mode inverse temperatures for independentharmonic oscillators (yellow curve). All measured in units of J . The frequency interval in which Im ˆ R ( ω ) is non-zero is ( ω − , ω + ) = (1 / (2 √ , / (2 √ (cid:39) (0 . , . . In Sec. 3 we recalled the relation between the p = 2 spherical disordered model and the classical integrablemodel introduced by Neumann. In this Section we will present some results concerning the behaviour of theintegrals of motion. A key issue we address here is how these influence the statistical properties in the steadystate. In Sec. 2.2.1 we studied the potential energy landscape of the p = 2 disordered model and we found that ithas N extrema corresponding to (cid:126)s = ±√ N(cid:126)v µ with (cid:126)v µ the N eigenvectors of the interaction matrix J ij . Thesedirections turn out to be the extrema of the integrals of motion landscape as well [88]. This claim is proveneasily. Indeed, ∂I µ ∂p η = 0 , (200) here we labeled the constants of motion defined in Eq. (48) with µ , the eigenvalue index, implies the followingtwo conditions s µ p η − s µ p µ s η for µ (cid:54) = η , mN (cid:88) ν ( (cid:54) = η ) s ν p η − p ν s ν s η λ ν − λ η for µ = η . It is clear that the first relation makes the second one hold identically. Moreover, replaced in the definition ofthe I µ s one finds − λ µ I µ = − λ µ s µ (201)that has to be extremised under the global spherical constraint on the s µ . This is just the analysis of the potentialenergy landscape that we performed in Sec. 2.2.1, leading to (cid:126)s ∗ = ±√ N(cid:126)v µ and z ∗ = λ µ for all the saddles in thelandscape. In our setting we use random initial conditions and we average over them. We will therefore focus on theaveraged integrals of motion, (cid:104) I µ (cid:105) = (cid:104) s µ (cid:105) + 1 mN (cid:88) ν ( (cid:54) = µ ) (cid:104) s µ p ν (cid:105) + (cid:104) s ν p µ (cid:105) − (cid:104) s µ p ν s ν p µ (cid:105) λ ν − λ µ . (202)Right after the instantaneous quench the initial values (cid:104) s µ (0 + ) (cid:105) and (cid:104) p µ (0 + ) (cid:105) are the ones right before thequench, (cid:104) s µ (0 − ) (cid:105) and (cid:104) p µ (0 − ) (cid:105) . Under the quench the eigenvalues transform as λ µ = J/J λ (0) µ . We cantherefore compute the (cid:104) I µ (0 + ) (cid:105) using these values. Owing to the fact that the initial conditions are drawnfrom an equilibrium probability density, we have (cid:104) s µ (0 + ) p ν (0 + ) (cid:105) = (cid:104) s µ (0 + ) (cid:105)(cid:104) p ν (0 + ) (cid:105) and (cid:104) p µ (0 + ) p ν (0 + ) (cid:105) = (cid:104) s µ (0 + ) s ν (0 + ) (cid:105) = (cid:104) s µ (0 + ) p ν (0 + ) (cid:105) = 0 for all µ (cid:54) = ν , and (cid:104) s µ (0 + ) p µ (0 + ) (cid:105) = 0 for all µ . The constants are then (cid:104) I µ (0 + ) (cid:105) = (cid:104) s µ (0 + ) (cid:105) + 1 mN (cid:88) ν ( (cid:54) = µ ) (cid:104) s µ (0 + ) (cid:105)(cid:104) p ν (0 + ) (cid:105) + (cid:104) s ν (0 + ) (cid:105)(cid:104) p µ (0 + ) (cid:105) λ ν − λ µ . (203) y < : condensed initial states. In the cases in which y < the initial state is condensed and the integrals of motion of the modes µ (cid:54) = N and µ = N scale very differently with N . For the modes in the bulk (cid:104) I µ (cid:54) = N (0 + ) (cid:105) = T (cid:48) z (0 − ) − λ (0) µ + T (cid:48) N (cid:104) s N (0 + ) (cid:105) λ N − λ µ + T (cid:48) N (cid:88) ν ( (cid:54) = µ,N ) z (0 − ) − λ (0) ν λ ν − λ µ + T (cid:48) N (cid:88) ν ( (cid:54) = µ ) z (0 − ) − λ (0) µ λ ν − λ µ . (204)This expression involves two sums that will appear again and again in the rest of this Section, and depend onlyon the pre-quench parameters J and T (cid:48) , S (1) µ ≡ N (cid:88) ν ( (cid:54) = µ ) λ (0) ν − λ (0) µ = 1 J N (cid:88) ν ( (cid:54) = µ ) λ ( n ) ν − λ ( n ) µ ≡ J S (1 n ) µ ,S (2) µ ≡ N (cid:88) ν ( (cid:54) = µ,N ) z (0 − ) − λ (0) ν λ (0) ν − λ (0) µ . For µ (cid:54) = N the second sum reads S (2) µ (cid:54) = N = 1 z (0 − ) − λ (0) µ (cid:18) J + S (1) µ (cid:19) = 1 z (0 − ) − λ (0) µ J (cid:16) S (1 n ) µ (cid:17) . (205)The superscript n means that the eigenvalues have been rescaled by J in such a way that they vary in theinterval ( − , and S (1 n ) µ is just a number. Using these definitions we find (cid:104) I µ (cid:54) = N (0 + ) (cid:105) = T (cid:48) z (0 − ) − λ (0) µ + J J (cid:34) T (cid:48) q in λ (0) N − λ (0) µ + T (cid:48) S (2) µ + T (cid:48) z (0 − ) − λ (0) µ S (1) µ (cid:35) . (206) his expression has a rather simple dependence on J that only appears as a prefactor in front of the sum ofthree terms within the square brackets. (This fact justifies the similarity of the bulk numerical values of (cid:104) I µ (cid:105) in,for example, panels (c) and (d) in Fig. 26.) In the large N limit z (0 − ) → λ (0) N and (cid:104) I µ (cid:54) = N (0 + ) (cid:105) = T (cid:48) λ (0) N − λ (0) µ (cid:20) J J + 2 T (cid:48) J S (1 n ) µ (cid:21) . (207) . . − −
200 0 2000 4000(b) h I N i / N NT = 0 . J = 0 . T = 0 . J = 1 . T = 0 . J = 1 . h I N i NT = 1 . J = 0 . T = 1 . J = 0 . T = 1 . J = 1 . Figure 24:
Uhlenbeck’s integral of motion (cid:104) I N (cid:105) for the largest mode, µ = N , as a function of system size for (a) threequenches with T (cid:48) /J < , and (b) three quenches with T (cid:48) /J > . We observe a clear linear scaling with system size in thecase of a condensed initial state, while the result is order 1 for a PM initial state. In panel (a) the dashed lines are the slopespredicted by Eq. (209).For the largest mode µ = N , if y < , we obtain (cid:104) I N (0 + ) (cid:105) = q in (cid:18) T (cid:48) J J S (1) N (cid:19) N + T (cid:48) J J S (2) N . (208)In the limit N → ∞ we can use S (1) N → − /J and the known form of q in to find (cid:104) I N (0 + ) (cid:105) (cid:55)→ (cid:18) − T (cid:48) J (cid:19) (cid:18) − T (cid:48) J (cid:19) N + T (cid:48) J J S (2) N ( J ) . (209)The parameter dependence of the slope in the right-hand-side seen as a function of N is verified numerically inFig. 24 (a). The good match between this form and the numerics indicates that S (2) N must be negligible withrespect to the first term that is O ( N ) ; indeed, we have checked numerically that S (2) N is O (1) . y < : paramagnetic initial states. If, instead, y > , there is no condensed mode and the integrals of motion are (cid:104) I µ (0 + ) (cid:105) = T (cid:48) z (0 − ) − λ (0) µ (cid:18) T (cid:48) J J S (1) µ (cid:19) + T (cid:48) J J S (2) µ = J z (0 − ) − λ (0) µ (cid:34) T (cid:48) J + T (cid:48) JJ + 2 T (cid:48) JJ S (1 n ) µ (cid:35) (210)and all of order . The particular case µ = N is displayed in the panel (b) of Fig. 24, showing no dependenceon N for N > ∼ , as expected. Summary.
Summarising, Fig. 24 displays the integral of motion associated to the mode at the edge of the spectrum for y < (a) and y > (b). All data points were obtained using a single realisation of the random matrix. For the ondensed initial conditions we clearly see the linear scaling with N . For the non-condensed ones the variationsshow a weak N dependence for small N plus a possible variation due to the fluctuations in the realisation of theeigenvalues. 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 areconserved (cid:104) I µ (0 + ) (cid:105) = (cid:104) I µ ( t ) (cid:105) for all µ , and that they satisfy the two constraints (cid:80) µ I µ (0 + ) = (cid:80) µ I µ ( t ) = N and (cid:80) µ λ µ I µ (0 + ) = (cid:80) µ λ µ I µ ( t ) = − e tot N , with e tot the total energy density. The analysis of the constants of motion should shed light on the “distance” from complete equilibration to aGibbs-Boltzmann probability density, especially in Sector III where shallow quenches in the N → ∞ Schwinger-Dyson formalism suggested proximity from this description. We here compare the actual values of the (cid:104) I µ (cid:105) tothe ones a system in Gibbs-Boltzmann equilibrium at a single temperature T f would have. The modes in the bulk
In a final Gibbs-Boltzmann equilibrium state the integrals of motion should read (cid:104) I µ (cid:54) = N (cid:105) T f = T f z f − λ µ + T f J J z f − λ µ S (1) µ + T f q f λ N − λ µ + T f J J S (2) µ , (211)where we allowed for the condensation of the largest mode. We wish to compare this expression to the one inEq. (207). After some lengthy calculations, using the N → ∞ values q in = 1 − T (cid:48) J , z f = 2 J = JJ z (0 − ) , T f = T (cid:48) (cid:18) JJ (cid:19) , q f = 1 − T f J , (212)we find that the difference ∆ I µ (cid:54) = N = (cid:104) I µ (cid:54) = N (0 + ) (cid:105) − (cid:104) I µ (cid:54) = N (cid:105) T f is ∆ I µ (cid:54) = N = − T (cid:48) (cid:18) J J − (cid:19) λ (0) N − λ (0) µ S (1) µ (213)a finite, non-zero, value for all µ (cid:54) = N . − −
200 0 2000 4000 ∆ I N NT = 0 . J = 0 . T = 0 . J = 1 . T = 0 . J = 1 . Figure 25:
Scaling of ∆ I N in sector III . The data points at five values of N are shown with points that are joined bycoloured straight lines. The dashed black lines are the predictions of Eq. (216). The N th mode he N th mode has a different scaling with system size. We have already computed (cid:104) I N (0 + ) (cid:105) in Eq. (209)and we want to compare it to what it should read in equilibrium at T f : (cid:104) I N (cid:105) f = q f (cid:18) T f J J S (1) N (cid:19) N + T f J J S (2) N . (214)The difference between the two, ∆ I N ≡ (cid:104) I N (0 + ) (cid:105) − (cid:104) I N (cid:105) T f , is ∆ I N = (cid:20) q in − q f + ( q in T (cid:48) − q f T f ) J J S (1) N (cid:21) N + (cid:18) T (cid:48) − T f J J (cid:19) J J S (2) N . (215)The first term diverges linearly with N . We have already argued that the second one is O (1) , see the discussionafter Eq. (209). Therefore, we focus on the slope of the difference, seen as a function of N . After replacing q f , q in and T f by their expressions in terms of J , J and T (cid:48) in the large N limit, and S (1) N → − /J , ∆ I N N → − T (cid:48) J (cid:18) J J − (cid:19) . (216)This form is validated by the numerical data in Fig. 25 for three pairs of T (cid:48) and J , all in Sector III.Quite clearly, the differences between the actual (cid:104) I µ (cid:105) and the ones in equilibrium at a temperature T f areproportional to J /J − , a factor that vanishes for J = J . Otherwise, the differences are finite for µ (cid:54) = N andproportional to N for the last mode, making the mode-by-mode difference non-vanishing for J (cid:54) = J . This factconfirms, then, the lack of equilibration to a Gibbs-Boltzmann measure with a single T f . The special line y = √ x 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 linearresponse are linked by FDT and the time-dependence of the correlation function after an instantaneous quenchis identical (within numerical accuracy) to the one found in equilibrium at this temperature. We shall brieflyanalyse the behaviour of the constants of motion in this case.On the special line y = √ x , T (cid:48) = JJ and Eq. (210) yields the following expression for the averaged integralsof motion (cid:104) I µ (0 + ) (cid:105) = 1 z (0 − ) /J − λ µ (cid:32)(cid:114) JJ + 1 + 2 S (1 n ) µ (cid:33) , (217)where λ µ = λ (0) µ /J are the normalised eigenvalues that vary in the interval ( − , . From Eq. (211) with q f = 0 , z f = 2 J valid in the N → ∞ limit, and using Eq. (205) the constants of motion in an equilibrium state attemperature T f are (cid:104) I µ (cid:105) T f = 12 − λ µ (cid:16) S (1 n ) µ (cid:17) . (218)We observe that the two sets of values are very similar if J is close to J . Numerically, we find that, in thespecial case T (cid:48) = 1 . J and J = 1 . J , the ∆ I µ is of order of − . We conclude that even in this case,in which the asymptotic analysis predicts that the global, mode-averaged quantities, behave as in equilibrium, aprediction that seems to be confirmed by the numerics, the constants of motion are not exactly the same in theinitial 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, theconserved energy dynamics of an isolated (finite size) system should keep memory of the initial conditions, evenif the system is non-integrable. In our problem, we see this information encoded in the I µ s. More so, not evenin the N → ∞ limit this memory is erased as the ∆ I µ s remain finite. We have seen in Sec. 5.7 that the N → ∞ system decouples into independent harmonic oscillators in theasymptotic long time limit (taken after N → ∞ ) since z ( t ) → z f . A natural idea is to check whether we canidentify the integrals of motion (cid:104) I µ (0 + ) (cid:105) with the ones that an ensemble of harmonic oscillators with springconstants mω µ = z f − λ µ in equilibrium at the temperatures T µ would have.
48 0 0 . . . . J h I µ i , T k i n µ T = 1 . J = 0 .
5, sec. I h I µ (0 + ) i T kin µ T kin T = 1 . J = 1 . h I µ (0 + ) i T kin µ T kin J h I µ i , T k i n µ µ/NT = 0 . J = 0 .
8, sec. III h I µ (0 + ) i T kin µ T kin µ/NT = 0 . J = 0 .
25, sec. IV h I µ (0 + ) i T kin µ T kin Figure 26:
Comparison between the averaged Uhlenbeck’s integrals of motion J (cid:104) I µ (cid:105) (red data points) and the mode tem-peratures T kin µ = 2 e kin µ (blue datapoints) of each mode in the four sectors of the phase diagram (we have checked that (cid:15) pot µ yield equivalent results). We used the same vertical scale in all plots, leaving aside the values of (cid:104) I µ (cid:105) close to the edge in (c)for which, in particular, (cid:104) I N (cid:105) (cid:39) . There is no such divergence at the edge of the spectrum in the other panels. The blackhorizontal lines represent the global kinetic temperatures T kin in the three sectors. They are located at T kin (cid:39) . (a), T kin (cid:39) . (b), T kin (cid:39) . (c) and T kin (cid:39) . (d), according to Eq. (121). We have also checked that N − (cid:80) µ T kin µ = T kin .On the one hand, we note that (cid:104) s µ ( t ) (cid:105) and (cid:104) p µ ( t ) (cid:105) are not constant for harmonic oscillators but their timeaverages are, so we evaluate (cid:104) I µ (cid:105) finding, in this framework, (cid:104) I µ (cid:105) osc = T µ z f − λ µ + T µ N (cid:88) ν ( (cid:54) = µ ) (cid:20) T ν z f − λ µ + T ν z f − λ ν (cid:21) λ ν − λ µ . (219)In the N → ∞ limit the value of z f is expected to be T (cid:48) + J /T (cid:48) for x < y and J for y < x . Imposing theidentity between the (cid:104) I µ (cid:105) osc and the (cid:104) I µ (0 + ) (cid:105) given in Eqs. (206) and (208) should yield a better estimate of thetemperatures 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 totalenergies, or the time-averaged kinetic and potential energies are also constants of motion. In the next Subsectionwe investigate to what extent the (cid:104) I µ (cid:105) are proportional to e tot µ = 2 e kin µ . .2.3 The integrals of motion and the mode temperatures Figure 26 shows the spectrum of integrals of motion J (cid:104) I µ (cid:105) (red data points) together with the mode kinetictemperatures T kin µ = 2 e kin µ (that, we have checked, are in agreement within numerical accuracy with the potentialones 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 (cid:104) I µ (cid:105) are noisier than theones 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 seethat the mode temperatures are very close to it again away from the edge of the spectrum and satisfying theconstraint N − (cid:80) µ T kin µ = T kin .The parameters in Sector II, panel (b), are on the curve y = √ x on which the data for the global correlationand linear response show equilibrium at a single temperature T f = J . The blue data points for the mode kineticenergies 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.In sector III, panel (c), the data for J (cid:104) I µ (cid:105) tend to be relatively flat in the bulk of the spectrum and veryclose to T kin = T f = T (cid:48) (1 + J/J ) / . This result can be derived from Eq. (206) taking advantage of a simplerearrangement of the two sums, S (2) µ = ( S (1) µ − S (1) N ) / ( z (0 − ) − λ (0) µ ) valid for µ (cid:54) = N . In particular, taking λ (0) µ right at the middle of the spectrum, S (1) µ = 0 by symmetry and J (cid:104) I (cid:105) = T f for N → ∞ . The incipientdivergence close to the right edge of the spectrum, with the deviation from this constant value, is also clear. Byusing a maximal value of 8 in the vertical axes we have explicitly left aside the points for µ (cid:39) N that take muchlarger values. For instance, J (cid:104) I N (cid:105) (cid:39) . The approximate mode temperatures in Eq. (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 isnot in proper Gibbs-Boltzmann equilibrium is due to the fact that the higher lying integrals of motion do notcomply with this temperature.The data for J (cid:104) I µ (cid:105) in Sector III look very similar to the ones found in Sector IV, see panel (d). Indeed, thedata almost coincide far from the edge, since they are both determined by Eq. (206) that has a weak dependenceon J (the only parameter that takes a different value in panels (c) and (d)). Close to the edge, the J (cid:104) I µ (cid:105) differ since in III (c) there is scaling with N while in IV (d) there is not. The mode kinetic energies are modeindependent and identical to T (cid:48) (1 + J/J ) / 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 = 2 e kin should be equal to the sum of the mode kinetic temper-atures T kin = 1 /N (cid:80) 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 afuture publication. The fact that the integrable system reaches a state that is very close to thermal equilibrium in sector III, thesector with the lowest total energy in the phase diagram, allows us to infer some properties of the phase spacestructure of the model. Only a system for which it is possible to visit all configurations with the same energyin the course of its dynamical evolution is capable of reaching thermal equilibrium. An integrable model cannotachieve this goal since it is bound to wander in a region of the phase space compatible with the values that theintegrals of motion (IOM) take on the initial configuration. The dynamics is constrained inside the phase spaceregion composed by configurations which have the same values of the I µ s ∀ µ . Such regions can be labeled withthe values of the I µ s (which also define the energy of the group of configurations since − (cid:80) λ µ I µ = 2 H ), and weshall call them iso-IOMs-regions in the following.A close-to-thermalised dynamics in an integrable system should be indicative of a substantial overlap betweenthe constant energy manifold and the equal IOMs region in phase space. This claim is, however, highly non-trivial since the equal energy manifold is N − dimensional while the iso-IOMs-region has only N − 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 . . − . − . . . − . − . . . − . − . . . . . − . − . . . σ µ e N = 20481024512 σ N − / N e N = 20481024512 σ N − / N e N = 20481024512 σ N / N e N = 20481024512 Figure 27:
Variance σ µ ( e ) for different system sizes. (a) Behaviour of the average over the first three fourths of thespectrum. (b), (c) σ µ /N for modes near the edge of the spectrum. (d) σ N /N . Note the logarithmic scale in the vertical axisin panels (a) and (d).the constant energy manifolds have a substantial overlap with any iso-IOMs-region that include a configurationwith the given energy. In order to test this guess, we studied the following quantity: σ µ ( e ) = (cid:90) N (cid:89) i =1 ds i dp i δ ( H [ s i , p i ] /N − e )( I µ [ s i , p i ] − (cid:104) I µ [ s i , p i ] (cid:105) ) , (220)where the average is a microcanonical one given by (cid:104) I µ [ s i , p i ] (cid:105) = (cid:90) N (cid:89) i =1 ds i dp i δ ( H [ s i , p i ] /N − e ) I µ [ s i , p i ] . (221)The quantity σ µ ( e ) measures how large are the fluctuations in the value of a given IOM I µ in the set ofconfigurations with the same energy e . According to the discussion at the beginning of this Section, if, for agiven 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 acanonical one, introducing a Lagrange multiplier β and a measure exp( − βH ) /Z , fixing the average energy densityof the ensemble. For large N , the fluctuations of the energy average are small and we get a good approximationto the microcanonical mean. The advantage of the canonical measure is that the (cid:104) I µ (cid:105) s are expressed in terms ofcanonical averages of the same kind as those which we were using in the previous Sections to describe the initial tate of the dynamics. Moreover, once z is fixed, the Hamiltonian is quadratic, and this allows one to expressthe higher order average (cid:104) I µ (cid:105) , which includes products of , and phase space variables, in terms of quadraticaverages (cid:104) s i (cid:105) and (cid:104) p i (cid:105) . A straightforward numerical calculation of σ µ ( e ) is then possible. We show numericalresults in Fig. 27. In panel (a) we plot σ µ , the average of σ µ ( e ) in the first three fourths of the spectrum σ µ = N/ (cid:88) µ =1 σ µ ( e ) . (222)We can clearly observe that it is very small for low energies and that it increases by several orders of magnitudeas we increase the energy density of the system. Such is the behaviour of σ µ ( e ) for the great majority of themodes. In panels (b) and (c) we show the behaviour of σ µ ( e ) for µ close to N , the edge of the spectrum. Weobserve that σ µ ( e ) exhibits a maximum at low energies. Finally, in panel (d), σ µ = N ( e ) is plotted. At very lowenergies, σ N ( e ) is very large and scales with N . As we increase energy, σ µ ( e ) decreases abruptly until it reachesa 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 smallfor sufficiently low energies. This means that the low energy configurations of the model have very similar valuesof the I µ s, at least for µ far from the edge of the spectrum. For modes close or at the edge of the spectrumfluctuations can be very large. These results suggest that the iso-IOMs-regions which lie at low energies havevery 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 thespectrum. In fact, in sector III, the sector with the lowest energies in the phase diagram, we have verified inprevious Section that the initial state and the final thermal state at temperature T f which partially describes thelong-time dynamics have very similar values of the IOMs far from the edge of the spectrum, with discrepanciesonly near the edge of the spectrum. This paper continues our study of the conserved energy dynamics following sudden quenches in classicaldisordered isolated models ruled by Newton dynamics.In the p = 3 strongly interacting case [7] all quenches reach an asymptotic regime in which a single (ordouble) temperature dynamical regime establishes. The systems either equilibrate to a paramagnetic state witha proper temperature, they remain confined in a metastable state with restricted Gibbs-Boltzmann equilibriumat a single temperature, or age indefinitely after a quench to the threshold with the dynamics being characterisedby one temperature at short time delays and another one at long time delays, similarly to what happens in therelaxational case [20, 72, 71, 73, 9]. The two temperatures T f and T eff depend on the pre and post quenchparameters 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). T s J . T d J . .
64 0 . .
75 1 1 .
25 1 . T ! J JJ T s J . T d J . .
64 0 . .
75 1 1 .
25 1 . χ st = T f χ st = − qT f + qT eff χ st = − qT f . .
52 0 0 . . T ′ J JJ χ st = T ′ q = 0 χ st = J q > χ st = J q = 0 Figure 28:
Dynamic phase diagrams. (a) p = 3 . (b) p = 2 . n 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 itsrelation 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 limitof the two-time correlation function. The phase diagram is shown in Fig. 28 (b), and it can be comparedto 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 diagramdepending on the initial state (being condensed or not) and the final value of the static susceptibility. Wereduced these four sectors to three phases, indicated with different colours in the Figure, in which thedynamic behaviour is different. Basically, they are distinguished by two “order parameters”, χ st and q , thestatic 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 nosingle temperature characterising the values taken by different observables in the long time limits, not evenafter being averaged over long time intervals. • There is one case, quenches from the condensed equilibrium state in which energy is extracted or injectedin 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 observablesexhibits the fact that the modes are not in Gibbs-Boltzmann equilibrium. Moreover, for deep quenches inthe same sector III one clearly sees that standard thermal equilibration is not reached. • Another special case is provided by quenches with T (cid:48) > J and T (cid:48) = JJ . On this special curve the globalobservables 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 aformalism that allows one to show that in the infinite size and long time limit (taken after the former one) themodes decouple and become independent harmonic oscillators. Once this regime is reached, the mode energiescan 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 avariation from sector to sector of the phase diagram. We compared these forms with the numerical measurementsand 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 testthe recently proposed relation between them and the frequency dependent effective temperature stemming fromthe fluctuation-dissipation relation of the spin observable [69, 70]. We found excellent agreement between thetwo 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 andstudied having only a few degrees of freedom. Here, as we are interested in searching for a statistical descriptionof the post-quench dynamics, we dealt with the limit of large, and even diverging, number of degrees of freedom.The N − integrals motion of the Neumann model have been identified by K. Uhlenbeck [16, 18]. After atrivial 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 alsoscale with N . The distance between their values and the ones they would have taken in equilibrium at a singletemperature 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 (cid:48) = JJ in sector or phase II, the integrals of motion are notidentical to the ones expected in equilibrium. This proves that not even in this case the system is able to fullyequilibrate.The N − integrals of motion could be used to build a putative Generalised Gibbs Ensemble, or they maybe a guideline to choose the ones with good scaling properties. We will investigate this problem in a sequel tothis publication. cknowledgements We are deeply indebted to O. Babelon for letting us know that we were dealing with Neumann’s model, and forilluminating discussions. LFC thanks M. Serbyn for useful discussions. We acknowledge financial support fromECOS-Sud A14E01 and PICS 506691 (CNRS-CONICET Argentina). LFC is a member of Institut Universitairede France. 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 thecorrelation C and linear response R functions derived in the N → ∞ limit. Based on a number of hypothesesthat 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 (cid:48) with parameters J , m 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 aconstant. This assumption does not apply to certain quenches of the isolated system. Still, we investigate theconsequences of the stationary assumption.
A.1.1 The asymptotic values
Let us assume that the limiting value of the Lagrange multiplier is a constant lim t →∞ z ( t ) = z f . (A.1)Recalling the definitions given in the main part of the paper, the limits of the correlation functions are q = lim t − t →∞ lim t →∞ C ( t , t ) , q = lim t →∞ C ( t , , q = lim t →∞ C ( t , . (A.2)The linear response was analysed in the main body of the paper and, independently of the quench parametersit is given by ˆ R ( ω ) = 12 J (cid:104) ( − mω + z f ) ± (cid:112) ( − mω + z f ) − J (cid:105) . (A.3)in the frequency domain, where the Fourier transform has been computed with respect to the time difference t − t . This result only assumes a long-time limit in which z f is time-independent. A.1.2 The parameters q and q We could estimate the asymptotic value of z ( t ) taking the long t limit of Eq. (78); however, without theuse of FDT we cannot compute the integral involved.We are tempted to propose q = q . (A.4)The interpretation of this result in terms of the evolution of real replicas is that the asymptotic value of theself-correlation between times t and , q , is the same as the asymptotic value of the correlation between tworeplicas evaluated at the same times t and with t 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 thiscontribution explicitly and we write C ( t , t ) = q + C st ( t − t ) with lim t − t →∞ C st ( t − t ) = 0 . We also use lim t →∞ C ( t ,
0) = q leaving the possibility of there being a different asymptotic value for this quantity open.The dynamic equation now becomes an equation on ˜ C that reads (cid:0) m∂ t − t + z f (cid:1) ( q + C st ( t − t )) = JJ T (cid:48) ( q − q ) + 2 J q lim t →∞ (cid:90) t dt (cid:48) R ( t − t (cid:48) )+ J (cid:90) t ds R ( t − t + s ) C st ( − s ) + J (cid:90) t ds C st ( t − t + s ) R ( s ) . (A.5)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 t → ∞ while keeping t − t fixed. The Fourier transform of this term with respect to t − t yields R ( − ω ) C st ( ω ) that using R ( − ω ) = ∗ ( ω ) simply becomes R ∗ ( ω ) C st ( ω ) . Proceeding in a similar way with the first integral one finds that it equals ˆ R ( ω ) C st ( ω ) . Thus, (cid:16) − mω + z f − J ˆ R ( ω ) (cid:17) C st ( ω ) = (cid:18) − z f q + JJ T (cid:48) ( q − q ) + 2 J q lim t →∞ (cid:90) t dt (cid:48) R ( t − t (cid:48) ) (cid:19) δ ( ω )+ J C st ( ω ) R ∗ ( ω ) . (A.6)The factor between parenthesis in the first term on the right-hand-side can be written as − z f q + JJ T (cid:48) ( q − q ) + 2 J q lim t →∞ (cid:90) t dt (cid:48) R ( t − t (cid:48) ) = 0 . (A.7)If we now assume q = q , this equation imposesif q (cid:54) = 0 then z f = 2 J and lim t →∞ (cid:90) t dt (cid:48) R ( t − t (cid:48) ) = 1 /J , (A.8)otherwise q = 0 . (A.9)The remaining equation, using Eq. (82) to replace the parenthesis on the left-hand-side by / ˆ R ( ω ) is recast as J C st ( ω ) | ˆ R ( ω ) | = C st ( ω ) . (A.10)At each frequency this equation has two possible solutions C st ( ω ) (cid:54) = 0 and | ˆ R ( ω ) | = 1 /J or C st ( ω ) = 0 . (A.11)In the frequency domain in which the linear response (A.3) is complex, that is Im ˆ R ( ω ) (cid:54) = 0 , one can easilycheck that it satisfies | ˆ R ( ω ) | = 1 /J . This holds for any set of parameters x, y .In the cases in which the linear response is real, it is not always true that its square equals /J . For instance,at zero frequency in cases in which x > y , ˆ R ( ω = 0) = 1 /J , verifying that its square is /J . However, for x < y , ˆ R ( ω = 0) = 1 /T (cid:48) and its square is different from /J . In these cases, the Fourier transform of the decayingpart of the correlation, ˆ C ( ω ) vanishes. We have tested this statement numerically and several examples can beseen in the main part of the paper. A.1.4 The correlation with the initial condition
The asymptotic limit of the equation for C ( t , implies z f q = J q lim t →∞ (cid:90) t dt (cid:48) R ( t − t (cid:48) ) + JJ T (cid:48) ( q − q in q ) (A.12)where q in = (cid:40) − T (cid:48) J T < T s T > T s (A.13)and T s the critical temperature of the initial potential energy Equation (A.12) is the equivalent of Eq. (82)in [7].This equation admits the solution q = q = 0 or z f = J lim t →∞ (cid:90) t dt (cid:48) R ( t − t (cid:48) ) + JJ T (cid:48) (1 − q in ) (A.14)if we assume q = q . Recalling that q in = 1 − T (cid:48) /J , the remaining equation becomes z f = J lim t →∞ (cid:90) t dt (cid:48) R ( t − t (cid:48) ) + J (A.15)and is consistent for z f = 2 J and lim t →∞ (cid:82) t dt (cid:48) R ( t − t (cid:48) ) = 1 /J . .2 The off-diagonal correlations in the no quench problem In the no quench case we can check whether C ( t ,
0) = q = q in for all t is a solution of the correspondingequation. Plugging this form in the evolution equation for C we find that, at t = 0 , either q in = 0 (theparamagnetic case) or z f (0) = 2 J /T (cid:48) (1 − q in ) and after replacing q in by its expression as a function of T (cid:48) /J thecorrect equilibrium z f = 2 J is recovered. However, for t = δ and later times, C ( t , cannot remain constant forinitial conditions with q in (cid:54) = 0 , due to the non-trivial contribution of the term J (cid:82) t dt (cid:48) R ( t , t (cid:48) ) q in = J /m q in δ .If we assume that C ( t , approaches a constant q possibly different from q in , the equation for C ( t , inthe long t limit reads z f q = J q lim t →∞ (cid:90) t dt (cid:48) R ( t , t (cid:48) ) + JJ T (cid:48) [ q q in + (1 − q in ) q ] (A.16)We have already shown q = q . Therefore, this equation has solution q = 0 or it becomes Eq. (A.15).In the no quench case we can use FDT and obtain z f q = J T (cid:48) q (1 − q ) + JJ T (cid:48) [ q q in + (1 − q in ) q ] (A.17)This is a new equation with no parallel in [7] since in this reference we only studied the dynamics starting fromparamagnetic-like initial states. A.3 Energy conservation
In this Section we impose that the kinetic and potential energies of the asymptotic final state correspond tothe ones of an equilibrium paramagnetic of condensed state at temperature T f . This means that e pot = − J /T f in the former case and e pot = − J /T f (1 − q ) with q = 1 − T f /J in the latter. We then derive, from theconservation of the total energy, the temperature T f as a function of the control parameters T (cid:48) , J , J or, inother words, x and y . Here, we will see the internal limits of validity of these assumptions. Elsewhere we willdetermine where it is realised numerically. A.3.1 Final temperature
For condensed (first line in the left) and paramagnetic (second line in the left) initial conditions going tocondensed (first line in the right) and paramagnetic-like (second line in the right) states the energy conservationlaw reads mm T (cid:48) − JJ T (cid:48) − (cid:16) − T (cid:48) J (cid:17) = T f − J T f − (cid:16) − T f J (cid:17) − JJ T (cid:48) ( q − q ) . (A.18)We have already determined q = q so the last term vanishes identically. From this condition we find thefollowing expressions for T f depending on the kind of quench performed: • From condensed to condensed T f = J (cid:18) mJ m J + 1 (cid:19) T (cid:48) J (III) . (A.19) • From condensed to paramagnetic T f = J J T (cid:48) mJ m J (cid:18) T (cid:48) J (cid:19) − (1 − q ) + (cid:118)(cid:117)(cid:117)(cid:116)(cid:34) mJ m J (cid:18) T (cid:48) J (cid:19) − (1 − q ) (cid:35) + (cid:18) T (cid:48) J (cid:19) (IV) . (A.20) • From paramagnetic to paramagnetic T f = J J T (cid:48) mJ m J (cid:18) T (cid:48) J (cid:19) − (cid:118)(cid:117)(cid:117)(cid:116)(cid:34) mJ m J (cid:18) T (cid:48) J (cid:19) − (cid:35) + (cid:18) T (cid:48) J (cid:19) (I) and (IIa) . (A.21) From paramagnetic to ageing T f = J (cid:20) mJ m J T (cid:48) J − J T (cid:48) (cid:21) (IIb) . (A.22)We have chosen to single out in these expression the parameters T (cid:48) /J and ( mJ ) / ( m J ) that we have usedin [7] to characterise the dynamic phase diagram of the p ≥ model. In the two cases (condensed to condensedand paramagnetic to paramagnetic) in which the no quench limit can be taken we naturally recover T f = T (cid:48) .In the ageing case T f should be the temperature of the stationary regime. As we will show from the numericalsolution 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 relationcan be realised for particular choices of the parameters only. Moreover, these conditions are not restrictiveenough, as the numerical solution of the full equations show that FDT is not realised for quenches allowed bythem. A.3.2 Limits of validity
Let us consider the cases y ≥ and y ≤ separately. The case y ≥ : Quench from a paramagnetic equilibrium state The double condition ≤ T f ≤ J translates into − ≤ mJ m J T (cid:48) J − J T (cid:48) ≤ (A.23)that in terms of x and y reads − ≤ y x − y ≤ (A.24)and yields ≤ y − x + 2 xy and y ≤ x . (A.25)The first condition in Eq. (A.25) is satisfied for ≤ ( y − y + )( y − y − ) (A.26)with y + and y − the two roots of the quadratic equation, that is y ± = − x ± (cid:112) x + x , (A.27)a positive and a negative value. One should then have y ≥ y + or y ≤ y − . (A.28)As y ≥ , y ≤ y − is not possible. Instead, y > y + is trivially satisfied. The second condition in Eq. (A.25) is theonly restrictive one and reads y ≤ √ x ≡ f ( x ) . (A.29)This upper bound is the dashed line representing f ( x ) in the phase diagram for y > . Note that this curve hasno finite limit. The case y ≤ : Quench from a condensed equilibrium state Using q in = 1 − T (cid:48) /J , one has T f J = 12 J J T (cid:48) J − T (cid:48) J (A.30)and the condition to use reads ≤ y (cid:18) x + 1 (cid:19) ≤ (A.31)The lower bound is trivially satisfied while the upper bound implies y ≤ x x ≡ g ( x ) for y ≤ . (A.32) he dashed line for y < in the phase diagram represents g ( x ) .In the condensed case the energy conservation implies T f = J y (1 + x )2 x , q = 1 − y − y x and e pot = − J (cid:18) − x y − y (cid:19) . (A.33)We note that q vanishes on the curve y = 2 x/ (1 + x ) , that it equals one for y = 0 and that, for fixed x > itincreases from ( x − / (2 x ) at y = 1 to 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 = lim t →∞ (cid:82) t dt (cid:48) R ( t − t (cid:48) ) = 1 /T (cid:48) for x < y and χ st = (cid:82) t R ( t ) = 1 /J for x > y . If we evaluate the integralof the linear response using FDT at a single temperature we obtain χ st = lim t →∞ (cid:90) t dt (cid:48) R ( t − t (cid:48) ) = lim t →∞ (cid:90) t dt (cid:48) T f ∂∂t (cid:48) C st ( t − t (cid:48) ) = T f y > T f (1 − q ) = 1 J y < (A.34)Using these results we can check whether there is a set of x, y for which T f derived in the previous Section fromthe conservation of the energy is consistent with these conditions. We list below the conclusions drawn in thefour parameter Sectors of the phase diagram. y > and x < y (Sector I) T f can equal T (cid:48) only for x = 1 that implies no-quench. y > and x > y (Sector II)Only on the curve y = √ x and y > FDT holds at T f . The validity of FDT at T f for parameters lying onthe curve y = √ x is verified numerically in Fig. 18. For all other values of x, y in this Sector FDT cannot besatisfied. y < and x > y (Sector III)In this Sector we do not find any contradiction. This reasoning suggests that the dynamics may satisfy FDTin this Sector. The no-quench case x = 1 for y < is obviously included. y < 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 harmonicpotential V ( x ) . The dynamics of this problem is given by the familiar equation m ¨ x + mω x = 0 , (B.1)with initial conditions x (0) = x and p (0) = p . B.1 Equilibrium initial conditions
Let us take initial conditions in canonical equilibrium within a harmonic potential V ( x ) . The probabilitydistribution of the initial conditions is P ( x , p ) = Z − e − β (cid:48) H = Z − e − β (cid:48) [ p m + V ( x )] (B.2) ith β (cid:48) = 1 / ( k B T (cid:48) ) the inverse temperature, using the same notation adopted in the main part of the paperfor the initial temperature T (cid:48) . The averaged kinetic energy of the ensemble of initial states sampled with thisprobability distribution function is m (cid:104) p (cid:105) = k B T (cid:48) , (B.3)the equipartition of the kinetic energy. We recall that angular brackets indicate average over the initial conditionssampled as above. The averaged total energy is (cid:104) H (cid:105) = − ∂∂β (cid:48) ln Z ( β (cid:48) ) . (B.4) B.2 Potential energy quench
Make now a quench in the potential that corresponds to V (cid:55)→ V and do it so quickly that the phase spacevariables do not change and remain p , x . By performing this abrupt change one injects or extracts a finiteamount of energy, ∆ E = H ( x , p ) − H ( x , p ) = V ( x ) − V ( x ) . (B.5)The energy surface on which the dynamics will take place is the one of the post-quench energy E (0 + ) = p / (2 m ) + V ( x ) .We now focalise on V being a harmonic potential. The Newton evolution of each initial configuration is x ( t ) = x cos ωt + p mω sin ωt ,p ( t ) = − mωx sin ωt + p cos ωt . (B.6)Let us call y = x ( t ) and z = p ( t ) the position and momentum at a time t . The probability density of y, z attime t is P ( y, z, t ) = (cid:90) dx (cid:90) dp P ( x , p ) δ ( y − x cos ωt − p mω sin ωt ) × δ ( z + mωx sin ωt − p cos ωt ) . We use the second δ function to integrate over p , P ( y, z, t ) = (cid:90) dx P (cid:16) x , z cos ωt + mωx tan ωt (cid:17) ωt × δ ( y − x cos ωt − z + mx ω sin ωtmω cos ωt sin ωt ) . The remaining δ function implies y − zmω tan ωt − x (cos ωt + tan ωt sin ωt ) = y − zmω tan ωt − x ωt = 0 (B.7)and we use it to integrate over x . Indeed, replacing x = y cos ωt − zmω sin ωt and taking care of the Jacobianone finds P ( y, z, t ) = P (cid:16) y cos ωt − zmω sin ωt, z cos ωt + mωy sin ωt (cid:17) . (B.8)In the case in which no quench is performed the initial potential has to be, then, harmonic V ( u ) = mω u / and the equation above implies ln Z + ln P ( y, z, t ) = − β i m ( z cos ωt + mωy sin ωt ) − β i mω (cid:16) y cos ωt − zmω sin ωt (cid:17) = − β i m z − β i mω y . The equilibrium distribution is conserved by the dynamics, as it should.Imagine now that the quench corresponds to a change in the spring parameter of the quadratic potential ω (cid:55)→ ω . The equation for P ( y, z ) implies ln Z + ln P ( y, z, t ) = − β i m ( z cos ωt + mωy sin ωt ) − β i mω (cid:16) y cos ωt − zmω sin ωt (cid:17) . xpanding the squares and collecting terms ln Z + ln P ( y, z, t ) = − β i (cid:18) (cos ωt + ω ω sin ωt ) z m + mω (sin ωt + ω ω cos ωt ) y +2 zymω (1 + ω ω ) cos ωt sin ωt (cid:19) . Although the measure P is still Gaussian it does not have the same covariance as the initial P .The averages and the variances of the position and momentum can be computed directly from the solutionsto the equations of motion. The averages vanish and for the variances one finds σ x ( t ) = (cid:104) x ( t ) (cid:105) = (cid:104) x (cid:105) cos ωt + (cid:104) p (cid:105) m ω sin ωt ,σ p ( t ) = (cid:104) p ( t ) (cid:105) = (cid:104) x (cid:105) m ω sin ωt + (cid:104) p (cid:105) cos ωt . (B.9)Replacing now the averages of the initial values (cid:104) p (cid:105) /m = mω (cid:104) x (cid:105) = T (cid:48) mω σ x ( t ) = mω (cid:104) x ( t ) (cid:105) = T (cid:48) (cid:18) ω ω cos ωt + sin ωt (cid:19) , m σ p ( t ) = 1 m (cid:104) p ( t ) (cid:105) = T (cid:48) (cid:18) cos ωt + ω ω sin ωt (cid:19) . (B.10)One readily verifies that, as expected, the averaged total energy is conserved (cid:104) E ( t ) (cid:105) = mω σ x ( t ) + 1 m σ p ( t ) = T (cid:48) (cid:18) ω ω (cid:19) for t > (B.11)since each trajectory does conserve its initial energy. The averaged total energy is, however, different from theone right before the quench, T (cid:48) = (cid:104) e tot ( t = 0 − ) (cid:105) (cid:54) = (cid:104) e tot ( t = 0 + ) (cid:105) = T (cid:48) (cid:0) ω /ω (cid:1) .Time-independent values of the variances are found from the average over a long time window: mω σ x ( t ) = T (cid:48) (cid:18) ω ω + 1 (cid:19) , m σ p ( t ) = T (cid:48) (cid:18) ω ω + 1 (cid:19) , x ( t ) p ( t ) = 0 (B.12)and from these one can identify the final temperature T f = T (cid:48) (cid:18) ω ω + 1 (cid:19) . (B.13) B.3 Fluctuation dissipation relations
The linear response of the harmonic oscillator defined in (B.1) to an infinitesimal perturbation modifying itsenergy as H (cid:55)→ H − hx , and therefore adding a linear term in h to Eq. (B.1) instantaneously at time t , is R ( t , t ) = δx ( t ) δh ( t ) (cid:12)(cid:12)(cid:12)(cid:12) h =0 = sin ω ( t − t ) mω θ ( t − t ) (B.14)while the product of the unperturbed position at the times t and t , after a long-time average over the referencetime t , is C ( t , t ) = x ( t ) x ( t ) = 12 (cid:18) x + p m ω (cid:19) cos ω ( t − t ) . (B.15)The Fourier transform with respect to the time delay, t − t → ν , of these two expressions yieldI m ˜ R ( ν ) = π mω δ ( ν − ω ) , (B.16) ν ˜ C ( ν ) = πνmω (cid:18) mω x + p m (cid:19) δ ( ν − ω ) . (B.17)Focusing on the frequency ν where both are non-zero due to the Dirac deltas, the ratio between the two yieldsthe ratio between the prefactors β FDR ( ν ) = 2 I m ˜ R ( ν ) ν ˜ C ( ν ) = 1 mω x + p m = 1 e tot (0 + ) for ν = ω , (B.18)where e tot is the total energy of the harmonic oscillator. .4 The Generalized Gibbs Ensemble The harmonic oscillator dynamics conserves its total energy e tot and the GGE density function is p GGE ( x, p ) = Z − e − β GGE H ( x,p ) (B.19)with Z = (cid:82) dxdp e − β GGE H ( x,p ) the normalisation constant or partition function. The inverse temperature of theGGE ensemble is determined by the condition (cid:104) H (cid:105) GGE = (cid:90) dxdp p GGE ( x, p ) H ( x, p ) = e tot (0 + ) (B.20)and with a simple calculation one finds β GGE = 1 e tot (0 + ) . (B.21)Therefore, β FDR ( ν = ω ) = β GGE . (B.22)The generalisation of this result to an ensemble of N harmonic oscillators is straightforward: β FDR ( ν = ω µ ) = β GGE µ for all µ . (B.23) C Discrete time version
Let us now consider the discrete time version of the C ( t , and z ( t ) equations and look for the discretisationneeded to recover the results C ( t ,
0) = q in and z ( t ) = 2 J at low temperatures.First of all, we use the Taylor expansions R ( δ,
0) = δ/m and C ( δ,
0) = 1 − δ T (cid:48) /m , for δ → .The discretized version of Eq. (77) evaluated at t = nδ is given by mδ [ C (( n + 1) δ, − C ( nδ,
0) + C (( n − δ, − z ( nδ ) C ( nδ,
0) + JJ T (cid:48) [ q in C ( nδ,
0) + (1 − q in ) C ( nδ, J δ n (cid:88) k =0 w ( n ) k R ( nδ, kδ ) C ( kδ, (C.1)where we approximated the integral by the discrete sum (cid:90) nδ dt (cid:48) R ( nδ, t (cid:48) ) C ( t (cid:48) ,
0) = δ n (cid:88) k =0 w ( n ) k R ( nδ, kδ ) C ( kδ, (C.2)with w ( n ) k coefficients the value pf which depend on the particular approximation used. The most commonmethod of approximation is the one given by the Newton-Cotes formulas, which, for example, gives for n = 1 , w (1)0 = and w (1)1 = (Trapezoidal rule), for n = 2 , w (2)0 = , w (2)1 = and w (2)2 = (Simpson’s rule), for n = 3 , w (3)0 = , w (3)1 = , w (3)2 = and w (3)3 = (Simpson’s rule), and so on.Consider now the equilibrium dynamics, that is, set J = J and assume stationarity conditions for C and RC ( t , t ) = C st ( t − t ) ,R ( t , t ) = R st ( t − t ) . (C.3)Moreover, suppose that z ( kδ ) = z = 2 J T (cid:48) (1 − q in ) and C ( kδ,
0) = q in for k = 0 , , ..., n .We want C ( t, to be a constant, so we have to enforce C (( n + 1) δ,
0) = q in . Using these assumptions,Eq. (C.1) can be rewritten as C (( n + 1) δ,
0) = q in + q in δ m (cid:34) − z + 2 J T (cid:48) (1 − q in ) − J T (cid:48) (1 − C st ( nδ )) + J δ n (cid:88) k =0 w ( n ) k R st (( n − k ) δ ) (cid:35) = q in + q in δ J m (cid:34) − T (cid:48) (1 − C st ( nδ )) + δ n (cid:88) k =0 w ( n ) k R st (( n − k ) δ ) (cid:35) = q in + q in δ J m (cid:34) − T (cid:48) (1 − C st ( nδ )) + δ n (cid:88) k =0 u ( n ) k R st ( kδ ) (cid:35) , (C.4) here we used the notation u ( n ) k = w ( n ) n − k (in the case of the coefficients of the Newton-Cotes formulas, one has w ( n ) k = w ( n ) n − k , that is, the coefficients appearing in symmetric positions in the sum are equal). Notice that, if q in = 0 , C ( t, is trivially constant, being always identical to zero independently of the choice of the coefficients w ( n ) k . To enforce C (( n + 1) δ,
0) = q in (cid:54) = 0 we need to satisfy the following equation δ n (cid:88) k =0 u ( n ) k R st ( kδ ) = 1 T (cid:48) (1 − C st ( nδ )) (C.5)that ensures that the terms between square brackets cancel identically. Take the case n = 1 , for example.Eq. (C.5) is equivalent to the condition C ( δ,
0) = C st ( δ ) = 1 − δ T (cid:48) u (1)1 R ( δ, . (C.6)The Taylor expansions C st ( δ ) = 1 − δ T (cid:48) /m and R st ( δ ) δ/m satisfy this equation when we use u (1)1 = 1 .Take now the generic n case and suppose that the discrete version of the FDT holds, namely R st ( nδ ) = − T (cid:48) δ [ C st ( nδ ) − C st (( n − δ )] . (C.7)Equation (C.5) reduces to − T (cid:48) n (cid:88) k =1 u ( n ) k [ C st ( kδ ) − C st (( k − δ )] = 1 T (cid:48) (1 − C st ( nδ )) (C.8)that can be rewritten as ( u ( n ) n − C st ( nδ ) + 1 − u ( n )1 + n − (cid:88) k =1 (cid:104) u ( n ) k − u ( n ) k +1 (cid:105) C st ( kδ ) = 0 (C.9)where we used R st (0) = 0 and C st (0) = 1 .As one can see, wether or not the left-hand-side indeed vanishes, implying C (( n + 1) δ,
0) = C ( nδ,
0) = . . . = C (0 , via Eq. (C.4), depends on the particular choice of the coefficients w ( n ) k used to approximate theintegrals. A simple way to satisfy Eq. (C.9) is to use w ( n ) k = 1 , for any n, k . We adopt, therefore, this rule toexpress the sums that represent the integrals.Let us now investigate what does this discretisation rule implies for the discrete time evolution of the Lagrangemultiplier z . The discrete version of Eq. (74) evaluated at t = nδ , and rewritten in a way that makes z eq appearis z ( nδ ) = 2 e f + 2 J T (cid:48) (1 − q ) + 2 J T (cid:48) (cid:2) C st ( nδ ) − (cid:3) + 4 J δ n (cid:88) k =0 w ( n ) k R st (( n − k ) δ ) C st (( n − k ) δ )= z eq + 2 J T (cid:48) (cid:2) C st ( nδ ) − (cid:3) + 4 J δ n (cid:88) k =0 u ( n ) k R st ( kδ ) C st ( kδ ) (C.10)where z eq = 2 e f + J T (cid:48) (1 − q ) = T (cid:48) + J T (cid:48) (1 − q ) is the equilibrium value of the Lagrange multiplier (for any q in ), u ( n ) k = w ( n ) n − k , and we have assumed stationarity for C and R .To enforce z ( nδ ) = z eq , we need to satisfy the following equation T (cid:48) (cid:2) C st ( nδ ) − (cid:3) + δ n (cid:88) k =0 u ( n ) k R st ( kδ ) C st ( kδ ) = 0 . (C.11)If one supposes that the discrete version of the FDT, Eq. (C.7), holds then Eq. (C.11) reduces to T (cid:48) (cid:2) C st ( nδ ) − (cid:3) − T (cid:48) n (cid:88) k =1 u ( n ) k [ C st ( kδ ) − C st (( k − δ )] C st ( kδ ) = 0 (C.12)that can be rewritten as (1 − u ( n ) n ) C st ( nδ ) − u ( n )1 C st ( δ ) − n − (cid:88) k =1 (cid:104) u ( n ) k C st ( kδ ) − u ( n ) k +1 C st (( k + 1) δ ) (cid:105) C st ( kδ ) = 0 . (C.13) f we choose u ( n ) k = 1 for any n, k , we obtain C st ( nδ ) + 1 − C st ( δ ) + n − (cid:88) k =1 (cid:8) C st ( kδ ) − C st (( k + 1) δ ) + [ C st (( k + 1) δ ) − C st ( kδ )] (cid:9) = 0 (C.14)that simplifies to n − (cid:88) k =0 [ C st (( k + 1) δ ) − C st ( kδ )] = 0 . (C.15) eferences [1] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. , 885 (2008).[2] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Phys. Rev. Lett. , 050405 (2007).[3] M. Rigol, V. Dunjko, and M. Olshanii, Nature , 854 (2008).[4] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. , 863 (2011).[5] P. Calabrese, J. Stat. Mech. , 064001 (2016).[6] C. Gogolin and J. Eisert, Rep. Prog. Phys. , 056001 (2016).[7] L. F. Cugliandolo, G. S. Lozano, and N. Nessi, J. Stat. Mech , P083301 (2017).[8] T. Scaffidi and E. Altman„ Semiclassical Theory of Many-Body Quantum Chaos and its Bound ,arXiv:1711.04768.[9] L. F. Cugliandolo, Slow Relaxations and nonequilibrium dynamics in condensed matter, In: Barrat JL.,Feigelman M., Kurchan J., Dalibard J. (eds), Les Houches Lecture Notes, arXiv:cond-mat/0210312 (2002).[10] A. Cavagna, Phys. Rep. , 51 (2009).[11] L. Berthier and G. Biroli, Rev. Mod. Phys. , 587 (2011).[12] A. Engel, Nucl. Phys. B [FS] , 617 (1993).[13] S. Franz and M. Mézard, EPL , 209 (1994).[14] L. F. Cugliandolo and P. Le Doussal, Phys. Rev. E , 152 (1996).[15] C. Neumann, Crelle Journal , 46 (1850).[16] K. K. Uhlenbeck, Spinger Lecture Notes in Mathematics , 146 (1982).[17] J. Avan and M. Talon, Int. J. Mod. Phys. A , 4477 (1990).[18] O. Babelon and M. Talon, Nucl. Phys. B , 321 (1992).[19] J. M. Kosterlitz, D. J. Thouless, and R. C. Jones, Phys. Rev. Lett , 1217 (1976).[20] L. F. Cugliandolo and J. Kurchan, Phys. Rev. Lett. , 173 (1993).[21] L. F. Cugliandolo and G. S. Lozano, Phys. Rev. Lett. , 4979 (1998).[22] L. F. Cugliandolo and G. S. Lozano, Phys. Rev. B , 915 (1999).[23] M. L. Mehta, Random matrices , Elsevier, Amsterdam, 2004.[24] M. Mueller and M. Wyart, Ann. Rev. Cond. Matt. Phys. , 177 (2015).[25] L. Laloux and J. Kurchan, J. Phys. A: Math. Gen. , 1929 (1996).[26] L. F. Cugliandolo and D. S. Dean, J. Phys. A: Math. Gen. , L453 (1995).[27] R. A. Adler, The Geometry of Random Fields , J. Wiley and Sons, New York, 1981.[28] A. J. Bray and D. S. Dean, Phys. Rev. Lett. , 150201 (2007).[29] Y. V. Fyodorov, Markov Proc. and Rel. Fields , 483 (2015).[30] G. Ben Arous and L. Sagun and V. Ugur Guney and Yann Le Cun, Explorations on high dimensionallandscapes , in International Conference on learning representations, ICLR, 2015.[31] D. J. Wales,
Energy Landscapes: With Applications to Clusters, Biomolecules and Glasses , CambridgeUniversity Press, Cambridge, U.K., 2004.[32] L. Susskind, hep-th/0302219.[33] M. R. Douglas, B. Shiffman, and S. Zelditch, Commun. Math. Phys. , 325 (2004).[34] Y. V. Fyodorov and B. A. Khoruzhenko, Proc. Nat. Ac. Sc. , 6827 (2016).[35] K. H. Fischer and J. A. Hertz,
Spin glasses , Cambridge University Press, Cambridge, 1991.[36] S. F. Edwards and R. C. Jones, J. Phys. A: Math. Gen. , 1595 (1976).
37] G. Semerjian and L. F. Cugliandolo, J. Phys. A: Math. Gen. , 4837 (2002).[38] F. Slanina, Phys. Rev. E , 011118 (2011).[39] G. Semerjian and L. F. Cugliandolo, Europhys. Lett. , 247 (2003).[40] R. Kuehn, Acta Physica Polonica B , 1653 (2015).[41] M. Mézard, G. Parisi, and M. A. Virasoro, Spin Glass Theory and Beyond: An Introduction to the ReplicaMethod and Its Applications , World Scientific, Singapore, 1986.[42] P. Shukla and S. Singh, J. Phys. C , L81 (1981).[43] S. Ciuchi and F. di Pasquale, Nucl. Phys. B [FS] , 31 (1988).[44] L. F. Cugliandolo and D. S. Dean, J. Phys. A: Math. Gen. , 4213 (1995).[45] Y. V. Fyodorov, A. Perret, and G. Schehr, J. Stat. Mech. , P11017 (2015).[46] A. Dembo, A. Guionnet, and C. Mazza, J. Stat. Phys. , 781 (2007).[47] A. J. Bray, Adv. Phys. , 357 (1994).[48] A. Onuki, Phase transition dynamics , Cambridge University Press, 2004.[49] S. Puri and V. Wadhawan, editors,
Kinetics of phase transitions , Taylor and Francis Group, 2009.[50] F. Corberi and P. Politi, Comptes Rendus de Physique , 255 (2015).[51] L. F. Cugliandolo, J. Kurchan, and L. Peliti, Phys. Rev. E , 3898 (1997).[52] L. F. Cugliandolo, J. Phys. A , 483001 (2011).[53] J. Berges, in Strongly interacting quantum systems out of equilibrium, T. Giamarchi, A. Millis, O. Parcolletand L. F. Cugliandolo (eds), Les Houches Lecture Notes, arXiv:1503.02907 (2015).[54] D. Boyanovsky, C. Destri, and H. J. de Vega, Phys. Rev. D , 045003 (2004).[55] D. Boyanovsky, H. J. de Vega, R. Holman, and J. Salgado, Phys. Rev. D , 125009 (1999).[56] B. Sciolla and G. Biroli, Phys. Rev. Lett. , 220401 (2010).[57] B. Sciolla and G. Biroli, J. Stat. Mech. , P11003 (2011).[58] B. Sciolla and G. Biroli, Phys. Rev. B , 201110 (2013).[59] A. Chandran, A. Nanduri, S. S. Gubser, and S. L. Sondhi, Phys. Rev. B , 024306 (2013).[60] A. Maraga, A. Chiocchetta, A. Mitra, and A. Gambassi, Phys. Rev. E , 042151 (2015).[61] A. Chiocchetta, M. Tavora, A. Gambassi, and A. Mitra, Phys. Rev. B , 134311 (2016).[62] A. Chiocchetta, A. Gambassi, S. Diehl, and J. Marino, Phys. Rev. Lett. , 135701 (2017).[63] O. Babelon, D. Bernard, and M. Talon, Introduction to Classical Integrable Systems , Cambridge UniversityPress, 2009.[64] M. Dunajski, Integrable systems, Cambridge University Lectures (2012).[65] V. I. Arnold,
Mathematical Methods of Classical Mechanics , Springer-Verlag, Berlin, 1978.[66] E. Yuzbashyan, Annals of Physics , 288 (2016).[67] A. Khinchin,
Mathematical foundations of statistical mechanics , Dover, New York, 1949.[68] R. Kubo, Rep. Prog. Phys , 255 (1966).[69] L. Foini, A. Gambassi, R. Konik, and L. F. Cugliandolo, Phys. Rev. E , 052116 (2017).[70] J. de Nardis et al., SciPost , 023 (2017).[71] A. Barrat, The p-spin spherical spin glass model, cond-mat/9701031, 1997.[72] L. F. Cugliandolo and J. Kurchan, Phil. Mag. B , 501 (1995).[73] T. Castellani and A. Cavagna, J. Stat. Mech. , P05012 (2005).[74] A. Houghton, S. Jain, and A. P. Young, Phys. Rev. B , 2630 (1983).
75] S. Franz and G. Parisi, J. Phys. I France , 1401 (1995).[76] A. Barrat, R. Burioni, and M. Mézard, J. Phys. A: Math. Gen. , L81 (1996).[77] V. Ermakov, Univ. Izv. Kiev, Series III , 1 (1880).[78] W. Milne, Phys. Rev. , 863 (1930).[79] E. Pinney, Proc. Am. Math. Soc. , 681 (1950).[80] S. Sotiriadis and J. Cardy, Phys. Rev. B , 134305 (2010).[81] M. Eckstein, M. Kollar, and P. Werner, Phys. Rev. Lett. , 056403 (2009).[82] M. Schiró and M. Fabrizio, Phys. Rev. Lett. , 076401 (2010).[83] N. Tsuji, M. Eckstein, and P. Werner, Phys. Rev. Lett. , 136404 (2013).[84] N. Tsuji and P. Werner, Phys. Rev. B , 165115 (2013).[85] C. Aron, G. Biroli, and L. F. Cugliandolo, J. Stat. Mech. , P11018 (2010).[86] P. Calabrese and A. Gambassi, J. Phys. A , R133 (2005).[87] F. Corberi, E. Lippiello, and M. Zannetti, J. Stat. Mech. , P07002 (2007).[88] O. Babelon, private communication., R133 (2005).[87] F. Corberi, E. Lippiello, and M. Zannetti, J. Stat. Mech. , P07002 (2007).[88] O. Babelon, private communication.