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

Let us label the eigenvalues of J_ij in such a way that they are ordered: $\lambda_1 \leqslant \lambda_2 \leqslant \dots \leqslant \lambda_N$. We call their associated eigenvectors $\vec {v}_\mu$ with $\mu=1, \dots, N$ and we take them to be orthonormal, such that $v_\mu^2 = 1$. We consider the potential energy landscape

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H_J[\{s_i\},z] =H_{\rm pot}[\{s_i\}] + H_{\rm const}[\{s_i\},z] \nonumber \end{align} \tag{ 12 }$$

with z taken as a variable.

In the absence of a magnetic field, all eigenstates of the interaction matrix are stationary points of the potential energy hyper-surface,

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \left. \frac{\partial H_J}{\partial s_i}\right|_{\vec {s}^*} = - \sum_{j (\neq i)}^N J_{ij} s_j + z s_i {\large |}_{\vec {s}^*,z^*} = 0 \;\;\;\; \forall i , \;\;\;\; \Rightarrow \;\;\;\; \vec {s}^*=\pm \sqrt{N} \vec {v}_\mu , \; z^*=\lambda_\mu ,\;\;\; \forall \mu . \nonumber \end{align*}$$

These stationary points are associated to metastable states at zero temperature, apart from two of them that are the marginally stable ground states, and their number is linear in N, the number of spins. (The role of marginal stability in the physical behaviour of condensed matter systems was recently summarised in [24].) These statements are shown in the way described in the next paragraph.

The Hessian of the potential energy surface on each stationary point is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left. \frac{\partial H_J}{\partial s_i \partial s_j}\right|_{\vec {s}^*, z^*} = \left. - J_{ij} + z \delta_{ij}\right|_{\vec {s}^*, z^*} = - J_{ij} + \lambda_\mu \delta_{ij} . \nonumber \end{align} \tag{ 13 }$$

This matrix can be easily diagonalised and one finds $\newcommand{\e}{{\rm e}} D_{\nu\eta}=(-\lambda_{\nu}+\lambda_\eta)\delta_{\nu\eta}$. Thus, on the stationary point, $\vec {s}^*=\pm\sqrt{N} \vec {v}_\mu$, the Hessian has one vanishing eigenvalue (for $\nu=\mu$), $\mu-1$ positive eigenvalues (for $\nu <\mu$), and $N-\mu$ negative eigenvalues (for $\nu > \mu$). Positive (negative) eigenvalues of the Hessian indicate stable (unstable) directions. This implies that each saddle point labeled by μ has one marginally stable direction, $\mu-1$ stable directions and $N-\mu$ unstable directions. (In other words, the number of stable directions plus the marginally stable one is given by the index μ labelling the eigenvalue associated to the stationary state.) In conclusion, there are two maxima, $\vec {s}^* = \pm \sqrt{N} \vec {v}_1$, in general pair of saddles $\vec {s}^* = \pm \sqrt{N} \vec {v}_I$ with $I=\mu-1$ stable directions and N  −  I unstable ones, with I running with μ as $I=\mu-1$ and $\mu=3, \dots, N$, and finally two (marginally stable) minima, $\vec {s}^* = \pm \sqrt{N} \vec {v}_N$. In the large N limit the density of eigenvalues of the Hessian at each stationary point μ is a translated semi circle law [25].

The zero temperature energy of a generic configuration under no applied field is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H_J[\vec {s}^* = \vec {v}_\mu),z^*] = - \frac{1}{2} \sum_{ij} J_{ij} s_i s_j + \frac{z}{2} \left( \sum_i s_i^2 - N\right) = -\frac{1}{2} \sum_\mu (\lambda_\mu -z) s_\mu^2 - \frac{z}{2} N . \nonumber \end{align} \tag{ 14 }$$

At each stationary point ${\vec {s}}^* = \pm \sqrt{N}\vec {v}_\mu$ and $z^*=\lambda_\mu$ this energy is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H_\mu^* \equiv H_J[\vec {s}^* = \pm\sqrt{N} \vec {v}_\mu,z^*=\lambda_\mu] = -\frac{N}{2} \sum_{i} v_i^\mu \sum_j J_{ij} v_j^\mu + \frac{z^*}{2} \left(N \sum_i (v_i^\mu)^2 - N \right) = -\frac{1}{2} \lambda_\mu N . \nonumber \end{align} \tag{ 15 }$$

Here we used the notation $v_i^\mu$ to indicate the ith component of the μth eigenvector $\vec {v}_\mu$. The energy difference between the minima and the lowest saddles depends on the distribution of eigenvalues, a semi-circle law for the Gaussian distributed interaction matrices that we consider here.

The analysis of large dimensional random potential energy landscapes [27–30] is a research topic in itself with implications in condensed matter physics, notably in glass theory [24, 31], but also claimed to play a role in string theory [32, 33], evolution [34] or other fields. The p  =  2 spherical model provides a particularly simple case in which the potential energy landscape can be completely elucidated.

This special (almost) quadratic model allows for the complete evaluation of its free-energy density for a typical realisation of the disordered exchanges. The traditional derivation of the disorder averaged free-energy density can also be done using the replica method and a simple replica symmetric Ansatz solves this problem completely. We recall how the two methods [19, 35] work in this section.

The partition function reads

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle Z_J = \prod_{i=1}^N \int_{-\infty}^\infty {\rm d}s_i \; {\rm e}^{\frac{\beta}{2} \sum\limits_{i\neq j} J_{ij} s_i s_j} \; \frac{1}{2\pi{\rm i}} \int_{c-{\rm i}\infty}^{c+{\rm i}\infty} {\rm d}z \; {\rm e}^{-\frac{\beta z}{2} \left(\sum\limits_{i=1}^N s_i^2 -N \right)} \nonumber \end{align*}$$

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 in terms of the projection of the spin vector $\vec {s}$ on the eigenvectors of J_ij, $\newcommand{\e}{{\rm e}} s_\mu \equiv \vec {s} \cdot \vec {v}_\mu$. This operation can be done for any particular realisation of the interaction matrix. The new variables $s_\mu$ are also continuous and unbounded and the partition function can be recast as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Z_J = \prod_{\mu=1}^N \int_{-\infty}^\infty {\rm d}s_\mu \; \frac{1}{2\pi{\rm i}} \int_{c-{\rm i}\infty}^{c+{\rm i}\infty} {\rm d}z \; {\rm e}^{\sum\limits_{\mu=1}^N \beta (\lambda_\mu-z) s_\mu^2/2 + \beta z N/2 } . \label{eq:part-funct-p2-2} \nonumber \end{align} \tag{ 16 }$$

Assuming that one can exchange the quadratic integration over $s_\mu$ with the one over the Lagrange multiplier, and that c is such that the influence of eigenvalues $\lambda_\mu>c$ is negligible, one obtains

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle Z_J =\frac{1}{2\pi{\rm i}} \int_{c-{\rm i}\infty}^{c+{\rm i}\infty} {\rm d}z \; {\rm e}^{-N \left[ -\beta z/2 +(2N)^{-1} \sum\limits_\mu \ln [\beta (z-\lambda_\mu)/2] \right]} . \label{eq:part-funct-p2-3} \nonumber \end{align} \tag{ 17 }$$

In the saddle-point approximation valid for $N\to\infty$ the Lagrange multiplier is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 1 = \langle\langle \; k_{\rm B}T \left(z_{\rm sp} - \lambda_\mu \right)^{-1} \; \rangle\rangle \label{eq:z-sp-statics} \nonumber \end{align} \tag{ 18 }$$

and this equation determines the different phases in the model. We indicate with double brackets the sum over the eigenvalues of the matrix J_ij that in the limit $N\to\infty$ can be traded for an integration over its density:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:sum_int} \frac{1}{N} \sum_{\mu=1}^N g(\lambda_\mu) = \int {\rm d}\lambda_\mu \; \rho(\lambda_\mu)\, g(\lambda_\mu) \equiv \langle \langle \, g(\lambda_\mu) \, \rangle \rangle. \nonumber \end{align} \tag{ 19 }$$

The high temperature solution to equation (18)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:z_eq_pm_mf} z_{\rm sp} = z_{\rm eq} = T + \frac{J^2}{T} \nonumber \end{align} \tag{ 20 }$$

can be smoothly continued to lower temperatures until the critical point

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (k_{\rm B}T_{\rm c})^{-1} = \langle\langle \; (z_{\rm sp}-\lambda_\mu)^{-1} \; \rangle\rangle \label{eq:p2-Ts} \nonumber \end{align} \tag{ 21 }$$

is reached where z_sp touches the maximum eigenvalue of the J_ij matrix, and it sticks to it for all $T<T_{\rm c}=J$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z_{\rm sp}= z_{\rm eq} = \lambda_{\rm max} = 2J \;\;\;\;\;\;\;\; T\leqslant T_{\rm c} . \nonumber \end{align} \tag{ 22 }$$

T_c is the static critical temperature. (A magnetic field with a component on the largest eigenvalue, $\vec {h} \cdot \vec {v}_{\rm max} \neq 0$, acts as an ordering field and erases the phase transition.)

If one now checks whether the spherical constraint is satisfied by these saddle-point Lagrange multiplier values, one verifies that it is in the high temperature phase, but it is not in the low temperature phase, where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sum_{\mu=1}^N \, \langle s_\mu^2 \rangle = \frac{T}{T_{\rm c}} N . \nonumber \end{align} \tag{ 23 }$$

The way out is to give a macroscopic weight to the projection of the spin in the direction of the eigenvector that corresponds to the largest eigenvalue:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle s_N = m_0 \sqrt{N} + \delta s_N = \sqrt{\left(1-\frac{T}{T_{\rm c}} \right) N} + \delta s_N \nonumber \end{align} \tag{ 24 }$$

with $\langle \delta s_N\rangle =0$ so that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sum_{\mu=1}^N \, \langle s_\mu^2 \rangle = \langle s_N^2 \rangle + \sum^{N-1}_{\nu=1} \langle s_\nu^2 \rangle = \left(1-\frac{T}{T_{\rm c}} \right) \, N + \frac{T}{T_{\rm c}} N = N. \nonumber \end{align} \tag{ 25 }$$

The thermal average of the projection of the spin vector on each eigenvalue vanishes in the high temperature phase and reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle s_{\mu} \rangle &= \left\{\begin{array}{@{}ll@{}} [N (1- T/T_{\rm c})]^{\frac12} & \;\;\;\;\;\; \lambda_\mu = \lambda_{\rm max} , \nonumber \\ 0 & \;\;\;\;\;\; \lambda_\mu < \lambda_{\rm max} , \end{array} \right. \nonumber \end{align} \tag{ 26 }$$

below the phase transition (once we have chosen one of the ergodic components with the spontaneous symmetry breaking of the $\vec {s} \to - \vec {s}$ invariance). The configuration condenses onto the eigenvector associated to the largest eigenvalue of the exchange matrix that carries a weight proportional to $\sqrt{N}$. Going back to the original spin basis, the mean magnetisation per site is zero at all temperatures but the thermal average of the square of the local magnetisation, that defines the Edwards–Anderson parameter, is not when $T<T_{\rm c}$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle m^2_i \rangle = 1-T/T_{\rm c} \;\;\;\; \Rightarrow \;\;\;\; q_{\rm EA} \equiv [\langle m^2_i \rangle]_J = 1-T/T_{\rm c} \;\;\;\; {\rm{with}} \;\;\;\; T_{\rm c}=J . \label{eq:p2-qea} \nonumber \end{align} \tag{ 27 }$$

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 distribution has long tails, as for the model is defined on a sparse random graph [36–38], the energy density diverges and the behaviour is more subtle [39, 40].

The disorder averaged free-energy density can also be computed using the replica trick [41] and a replica symmetric Ansatz. This Ansatz corresponds to an overlap matrix between replicas $\newcommand{\e}{{\rm e}} Q_{ab} = \delta_{ab} + q_{\rm EA} \epsilon_{ab}$ with $\newcommand{\e}{{\rm e}} \epsilon_{ab} =1$ for $a\neq b$ and $\newcommand{\e}{{\rm e}} \epsilon_{ab} =0$ for a  =  b. When $N\to\infty$ the saddle point equations fixing the parameter q_EA yield 0 above T_c and a marginally stable solution with $q_{\rm EA}=1-T/T_{\rm c}$ and identical physical properties to the ones discussed above below T_c.

The equilibrium energy is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {e}^{\rm pot}_{\rm eq} &= \left\{\begin{array}{@{}ll@{}} - \frac{J^2}{2T} \left[ 1- \left(1-\frac{T}{J}\right)^2 \right] = \frac{1}{2} (k_{\rm B}T - \lambda_{\rm max}) & \qquad \qquad T<T_{\rm c} , \nonumber \\ -\frac{J^2}{2T} & \qquad \qquad T>T_{\rm c} . \end{array} \right. \label{eq:epot-equil} \nonumber \end{align} \tag{ 28 }$$

We added here a superscript ${\rm pot}$ since in the modified model that we will study in this paper the total energy will also have a kinetic energy contribution. The entropy diverges at low temperatures as $\ln T$, just as for the classical ideal gas, as usual in classical continuous spin models.

In the $N\to\infty$ limit, the only relevant correlation and linear response functions are

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \qquad \qquad C_{ab}(t_1,t_2) = N^{-1} \sum_{i=1}^N \; [\langle s^a_i(t_1) s^b_i(t_2)\rangle] , \nonumber \end{align} \tag{ 62 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \qquad \qquad C_{ab}(t_1,0) = N^{-1} \sum_{i=1}^N \; [\langle s^a_i(t_1) s^b_i(0)\rangle] , \nonumber \end{align} \tag{ 63 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \qquad \qquad R_{ab}(t_1,t_2) = \left. N^{-1} \frac{\delta \;\;\;}{\delta h_b(t_2)} \sum_{i=1}^N \; [\langle {s^a}^{(h)}_{i}(t_1)\rangle] \right|_{h=0} , \nonumber \end{align} \tag{ 64 }$$

for $t_1, \ t_2 > 0$, where the infinitesimal perturbation h is linearly coupled to the spin $H\mapsto H - h \sum_{i=1}^N s_i$ at time 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_1>t_2$. The square brackets denote here and everywhere in the paper the average over quenched disorder. The angular brackets indicate the average over thermal noise if the system is coupled to an environment, and over the initial conditions of the dynamics sampled, say, with a probability distribution. When the coupling to the bath is set to zero, as we do in this paper, the last average is the only one remaining in the angular brackets operation. The meaning of the indices $a, b$ is given in the next paragraph.

The dynamical equations starting from a random state are well-known and can be found in [20, 71–73]. They are usually derived from the dynamical Martin–Siggia–Rose–Janssen–deDominicis generating function. The method has been modified to take into account the effect of equilibrium initial conditions in [74] and it was applied to the relaxational p spin model in [75, 76]. The average over disorder now becomes non-trivial and needs the use of the replica trick. The scripts $a, b$ indicate then the replica indices $a, b = 1, \dots, N$. For initial conditions in equilibrium the replica structure is replica symmetric (see section 2.2 and [19]), with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C_{aa}(0,0) = Q_{aa}=1 \qquad {\rm{and}} \qquad C_{a \, b\neq a}(0,0) = Q_{a \, b\neq a}=q_{\rm in} \nonumber \end{align} \tag{ 65 }$$

and $q_{\rm in}=0$ in the paramagnetic state while $q_{\rm in}\neq 0$ in the condensed phase. This structure has an effect on the equation for the time-dependent correlation function that will keep the initial replica structure. There will be two kinds of correlations with the initial condition

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C_1(t,0) \equiv C_{a=1 \, b=1}(t,0) \qquad {\rm{and}}\qquad C_{b\neq 1}(t,0) \equiv C_{a=1 \, b\neq1}(t,0) , \nonumber \end{align} \tag{ 66 }$$

where we singled out the replica labeled one and, with the definitions, we shortened the notation from two to one subscripts as there should not be room for confusion here. Since there is no reason to think that the replicas that are not the one numbered one behave differently, we follow the dynamics of

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C_2(t_1,t_2) \nonumber \end{align} \tag{ 67 }$$

as a representative of this group. The interpretation of the correlations $C_1(t_1, t_2)$ and $C_2(t_1, t_2)$ can be given in terms of real replicas. The former is the self-correlation between the configuration of one replica of the system $\{s_i\}(t_2)$ and the same replica evolved until a later time t₁, $\{s_i\}(t_1)$. For this reason, we will eliminate the subscript ₁ and call $C_1(t_1, t_2) \mapsto C(t_1, t_2)$ in most places hereafter. The latter is the correlation between replica labeled 2, let us say $\{\sigma_i\}(t_2)$, at time t₂ and the singled out replica one evolved until time t₁ and represented by $\{s_i\}(t_1)$. Although we could write an evolution equation for the two-time $C_2(t_1, t_2)$, actually we do not need it since only $C_2(t_1, 0)$, the correlation with replica 2 evaluated at the initial time, intervenes in the other equations. This is the only correlation between different replicas that will appear in the calculations.

In the $N\to\infty$ limit one derives the dynamical equations that read

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left(m\partial_{t_{1}}^{2}+z(t_{1})\right)R(t_{1},t_{2}) = J^2 \int_{t_{2}}^{t_{1}}{\rm d}t'R(t_{1},t') R(t',t_{2}) +\delta(t_1-t_2), \nonumber \end{align} \tag{ 68 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left(m\partial_{t_{1}}^{2}+z(t_{1})\right)C(t_{1},t_{2}) =&~ J^2\int_{0}^{t_{1}}{\rm d}t'R(t_{1},t') C(t',t_{2}) \nonumber \\ & + J^2\int_{0}^{t_{2}}{\rm d}t'R(t_{2},t')C(t_{1},t') +\frac{J J_0}{T'} \sum_a C_a(t_{1},0) C_a(t_2,0) , \label{eq:dyn-eqs-C-0} \nonumber \end{align} \tag{ 69 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left(m\partial_{t_{1}}^{2}+z(t_{1})\right)C_a(t_{1},0) = J^2\int_{0}^{t_{1}}{\rm d}t' \; R(t_{1},t') C_a(t',0) +\frac{J J_0}{T'} \sum_b C_b(t_{1},0) Q_{ab} , \label{eq:dyn-eqs-C1-0} \nonumber \end{align} \tag{ 70 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z(t_{1}) =& -m\partial_{t_{1}}^{2}C(t_{1},t_2)\vert_{t_{2}\rightarrow t_{1}^{-}} + 2J^2\int_{0}^{t_{1}}{\rm d}t'R(t_{1},t')C(t_{1},t') \nonumber \\ & + \frac{J J_0}{T'} \sum_a (C_a(t_{1},0))^2 \label{eq:dyn-eqs-z-0} . \nonumber \end{align} \tag{ 71 }$$

The border conditions are

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}ll@{}} &C_1(t_1,0)= C(t_1,0) \qquad {\rm{implying}} \qquad C_1(0,0)= C(0,0) =1 , \nonumber \\ & C_2(0,0) = Q_{1,2} \qquad \qquad {\rm{implying}} \qquad C_{2}(0,0) = q_{\rm in} . \end{array} \nonumber \end{align} \tag{ 72 }$$

Note that the initial condition is not the same for all $C_b(t_1, 0)$. It is equal to 1 for b  =  1 and equal to q_in for $b\neq 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 correspond to p  =  2 and they have the extra ingredient of the influence of equilibrium initial conditions with a non-trivial replica structure, allowing for condensed initial states in proper thermal equilibrium.

The sums over the replica indices appearing in equations (69)–(71) can be readily computed in the $n\to 0$ limit; they read

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{J J_0}{T'} \sum_a C_a(t_{1},0) C_a(t_2,0) &= \frac{J J_0}{T'} \left[ C_1(t_1,0) C_1(t_2,0) + (n-1) C_2(t_1,0) C_2(t_2,0) \right] \nonumber \\ &\mathop{\to}\limits_{n\to 0} \frac{J J_0}{T'} \left[ C_1(t_1,0) C_1(t_2,0) - C_2(t_1,0) C_2(t_2,0) \right] , \label{sum1} \nonumber \end{align} \tag{ 73 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{J J_0}{T'} \sum_b C_b(t_{1},0) Q_{1b} &= \frac{J J_0}{T'} \left[ C_1(t_1,0) \, 1 + (n-1) C_{b(\neq 1)}(t_1,0) q_{\rm in} \right] \nonumber \\ &\mathop{\to}\limits_{n\to 0} \frac{J J_0}{T'} \left[ C_a(t_1,0) - C_{b(\neq 1)}(t_1,0) q_{\rm in} \right] , \label{sum2} \nonumber \end{align} \tag{ 74 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{J J_0}{T'} \sum_b C_b(t_{1},0) Q_{(a\neq 1)b} &= \frac{J J_0}{T'} \left[ C_{1}(t_1,0) q_{\rm in} + C_2(t_1,0) \, 1 + (n-2) C_{b(\neq (1,2))}(t_1,0) q_{\rm in} \right] \nonumber \\ &\mathop{\to}\limits_{n\to 0} \frac{J J_0}{T'} \left[ q_{\rm in} C_1(t_1,0) + (1-2q_{\rm in}) C_2(t_1,0) \right]. \label{sum3} \nonumber \end{align} \tag{ 75 }$$

Consequently, the terms induced in the equation for $C(t_1, t_2)$ and $C_2(t_1, 0)$ are different.

With inertia and no coupled bath, the equal-time conditions are

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C(t_1,t_1) &= 1 , \nonumber \\ R(t_1,t_1) &= 0 , \nonumber \\ \partial_{t_1}C(t_1,t_2)\vert_{t_2\rightarrow t_1^{-}}=\partial_{t_1}C(t_1,t_2)\vert_{t_2\rightarrow t_1^{+}} &= 0 , \nonumber \\ \partial_{t_1}C_2(t_1,t_2)\vert_{t_2\rightarrow t_1^{-}}=\partial_{t_1}C_2(t_1,t_2)\vert_{t_2\rightarrow t_1^{+}} &= 0 , \nonumber \\ \partial_{t_1}C_2(t_1,0)\vert_{t_1\rightarrow 0^+} &= 0 , \nonumber \\ \partial_{t_1}R(t_1,t_2)\vert_{t_2\rightarrow t_1^{-}} &= \frac{1}{m} , \nonumber \\ R(t_1,t_2)\vert_{t_2\rightarrow t_1^{+}} &= 0 , \nonumber \end{align} \tag{ 76 }$$

for all times $t_1, t_2$ larger than or equal to 0⁺ , when the dynamics start.

We found convenient to numerically integrate the integro-differential equations using an expression of the Lagrange multiplier that trades the second-time derivative of the correlation function into the total conserved energy after the quench. Following the same steps explained in [7] we deduce

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z(t_{1})=2e(t_1)+ 4J^2 \int_{0}^{t_{1}}{\rm d}t' \; R(t_{1},t')C(t_{1},t')+ \frac{2JJ_0}{T'} \ \left[ C^2(t_{1},0) -C_2^2(t_{1},0)\right] , \label{eq:z-Nicolas} \nonumber \end{align} \tag{ 77 }$$

where we used equation (73) evaluated at $t_1=t_2$. It seems that we have simply traded z(t₁) by e(t₁). Indeed, taking advantage of the fact that for an isolated system $e(t_1)=e_f$, a constant, the numerical solution of the evolution equations now becomes easier since it does not involve the second time derivative of the correlation function. In practice, in the numerical algorithm we fix the total energy e_f to its post-quench value derived in section 4.3, and we then integrate the set of coupled integro-differential equations with a standard Runge–Kutta method.

In short, the set of equations that fully determine the evolution of the system from an initial condition in canonical Boltzmann equilibrium at any temperature $T'$ are

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left(m\partial_{t_{1}}^{2}+z(t_{1})\right)R(t_{1},t_{2}) = \delta(t_1-t_2) + J^2 \int_{t_{2}}^{t_{1}}{\rm d}t'R(t_{1},t') R(t',t_{2}) , \label{eq:dyn-eqs-R} \nonumber \end{align} \tag{ 78 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left(m\partial_{t_{1}}^{2}+z(t_{1})\right)C(t_{1},t_{2}) =&~ J^2\int_{0}^{t_{1}}{\rm d}t'R(t_{1},t') C(t',t_{2}) + J^2\int_{0}^{t_{2}}{\rm d}t'R(t_{2},t')C(t_{1},t') \; \nonumber \\ & +\frac{J J_0}{T'} \; [C(t_{1},0) C(t_2,0) - C_2(t_1,0) C_2(t_2,0)] , \label{eq:dyn-eqs-C} \nonumber \end{align} \tag{ 79 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left(m\partial_{t_{1}}^{2}+z(t_{1})\right)C_2(t_{1},0) = J^2\int_{0}^{t_{1}}{\rm d}t' \; R(t_{1},t') C_2(t',0) + \frac{J J_0}{T'} \left[ q_{\rm in} C(t_1,0) + (1-2q_{\rm in}) C_2(t_1,0) \right] , \label{eq:dyn-eqs-C1} \nonumber \end{align} \tag{ 80 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z(t_{1}) =&~ -m\partial_{t_{1}}^{2}C(t_{1},t_2)\vert_{t_{2}\rightarrow t_{1}^{-}} + 2J^2\int_{0}^{t_{1}}{\rm d}t'R(t_{1},t')C(t_{1},t') \nonumber \\ & + \frac{J J_0}{T'} \; [C^2(t_{1},0) - C^2_2(t_1,0)] \label{eq:dyn-eqs-z} . \nonumber \end{align} \tag{ 81 }$$

High and low temperature initial states are distinguished by $q_{\rm in}=0$ for $T'>J_0$, and $q_{\rm in} =1-T'/J_0$ for $T'<J_0$, respectively. The equation for C(t₁,0) is just the one for $C(t_1, t_2)$ evaluated at t₂  =  0 so we do not write it explicitly.

The equation for C₂(t,0) admits the solution $C_2(t_1, 0)=q_{\rm in} = {\rm{cst}}$ in the case in which no quench is performed. Indeed, setting $C_2(t_1, 0)=q_{\rm in}$ in equation (80) one has

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z(t_1) q_{\rm in} = J^2 q_{\rm in} \int_0^{t_1} {\rm d}t' \, R(t_1,t') + \frac{J^2}{T'} [ q_{\rm in} C(t_1,0) + (1-2 q_{\rm in}) q_{\rm in}]. \nonumber \end{align} \tag{ 82 }$$

This equation has the solution $q_{\rm in}=0$, the one of the paramagnetic phase, and a non-vanishing $q_{\rm in}\neq 0$ solution relevant in the ordered phase. Let us now focus on the case $q_{\rm in}\neq 0$. Using FDT, a property of equilibrium, the integral can be performed, the contribution from $t'=0$ cancels the first term in the square brackets, and the one from $t'=t_1$ combines with the second term in the square bracket; the ensuing equation simplifies to read

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z(t_1) = \frac{2J^2}{T'} (1 - q_{\rm in}) = z_f. \nonumber \end{align} \tag{ 83 }$$

Therefore, z(t₁) is also a constant. The equation for z(t₁), using FDT, becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z(t_1) &= T' + \frac{J^2}{T'} [C^2(t_1,t_1) - C^2(t_1,0)] + \frac{J^2}{T'} [C^2(t_1,0) - q_{\rm in}^2] \nonumber \\ &= T' + \frac{J^2}{T'} (1 - q_{\rm in}^2) \ = \ z_f . \nonumber \end{align} \tag{ 84 }$$

The two equations yield the low-temperature values z_f  =  2J and $q_{\rm in}=1-T'/J$.

In equilibrium C₂(t,0) and z(t₁) remain constant and equal to their initial values, q_in and z_f.

Knowing that z(t₁) remains constant in equilibrium, one can easily analyse the response equation in Fourier space. The equation that determines its dynamical evolution is transformed into

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat R(\omega) = {\frac{1}{-m \omega^2 + z_f - \Sigma(\omega)} } . \label{Rstat} \nonumber \end{align} \tag{ 85 }$$

For this model $\Sigma(\omega)=J^2 \, \hat R(\omega)$ and then

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat R(\omega) = \frac{1}{2 J^2} \left[ (-m\omega^2 + z_f) \pm \sqrt{(-m\omega^2 + z_f)^2 - 4 J^2} \right]. \label{eq:response_fourier_transform-text} \nonumber \end{align} \tag{ 86 }$$

$\hat R(\omega)$, and also $R(t)$, are independent of temperature for $T'<T_{\rm c}=J$ while they depend on temperature through z_f for $T'>T_{\rm c}=J$.

The terms under the square root can be more conveniently written as functions of the special values

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle m\omega_{\pm}^2 = z_f\pm 2J \label{eq:omega_pm} \nonumber \end{align} \tag{ 87 }$$

and the linear response is recast as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat R(\omega) = \frac{1}{2 J^2} \left[ -m\omega^2 + z_f \pm m \sqrt{(\omega_-^2-\omega^2) (\omega_+^2-\omega^2)} \right]. \label{eq:response_fourier_transform-text1} \nonumber \end{align} \tag{ 88 }$$

The imaginary and real parts of $\hat R(\omega)$ are then

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm{Im}} \hat R(\omega) &=& \frac{1}{2J^2} \left\{\begin{array}{@{}ll@{}} \ -m \sqrt{(\omega^2-\omega^2_-)(\omega_+^2-\omega^2) } \qquad & \qquad \qquad \;\; {\rm{for}} \;\; \qquad \omega_-\leqslant \omega \leqslant \omega_+ \nonumber \\ 0 \qquad \qquad & \qquad \qquad \;\; {\rm{otherwise}} \end{array} \right. \nonumber \end{align} \tag{ 89 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm{Re}} \hat R(\omega) &=& \frac{1}{2J^2} \left\{\begin{array}{@{}ll@{}} z_f - m\omega^2 - m \sqrt{(\omega^2_- -\omega^2)(\omega_+^2-\omega^2)} & \qquad {\rm{for}} \;\; \omega \leqslant \omega_- \nonumber \\ z_f - m\omega^2 & \qquad {\rm{for}} \;\; \omega_-\leqslant \omega \leqslant \omega_+ \nonumber \\ z_f - m\omega^2+ m \sqrt{(\omega^2 - \omega^2_-)(\omega^2 - \omega_+^2)} & \qquad {\rm{for}} \;\; \omega_+ \leqslant \omega \end{array} \right. \nonumber \end{align} \tag{ 90 }$$

(note the unusual choice of sign for the imaginary part that we adopted.) In terms of the physical parameters, $\hat R(\omega)$ is real for $|-m\omega^2+z_f|> 2J$. In the low temperature phase, since z_f  =  2J, this implies $\omega_-=0$ and the imaginary part of the linear response is gapless. One can easily check that $|\hat R(\omega)|^2 = 1/J^2$ in the interval $\omega_-\leqslant \omega \leqslant \omega_+$. Away from this interval the modulus of the linear response is a complicated function of the frequency.

The zero frequency linear response

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat R(\omega=0) = \int_0^\infty {\rm d}\tau \, R(\tau) = {\frac{1}{2 J^2}} \left[ z_f \pm \sqrt{z_f^2 - 4 J^2} \right] , \label{eq:Romega0} \nonumber \end{align} \tag{ 91 }$$

with τ a time delay, must be a real quantity and this form is a manifestation of the condition $z_f \geqslant 2J$. The static susceptibility, or zero frequency response, is then given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi_{\rm st} \equiv \hat R(\omega=0)= \left\{\begin{array}{@{}ll@{}} 1/J & \qquad {\rm{for}} \quad T' < T_{\rm c}=J \nonumber \\ 1/T' & \qquad {\rm{for}} \quad T' > T_{\rm c}=J \end{array} \right. \nonumber \end{align} \tag{ 92 }$$

with the result in the second line being ensured by the choice of minus sign in front of the square root. The frequency-dependent linear response can then be transformed back to real-time and thus get its full time-evolution.

In the lower limit of the spectrum the imaginary part of the linear response goes as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm{Im}} \hat R(\omega_{-} + \epsilon) \sim \left\{\begin{array}{@{}ll@{}} {- \frac{1}{J} \sqrt{\frac{m}{J}}} \; \ \epsilon \qquad & \quad {\rm{at}} \quad T<T_{\rm c} \nonumber \\ {- \frac{(2 m J \omega_{-})^{1/2}}{J^2} } \sqrt{\epsilon} & \quad {\rm{at}} \quad T>T_{\rm c} \end{array} \right. \quad & \qquad {\rm{as}} \qquad \epsilon \rightarrow 0^{+} \nonumber \end{align} \tag{ 93 }$$

while in the upper limit it vanishes as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm{Im}} \hat R(\omega_{+} - \epsilon) \sim {- \frac{(2 m J \omega_{+})^{1/2}}{J^2} } \sqrt{\epsilon} \qquad {\rm{as}} \qquad \epsilon \rightarrow 0^{+} \nonumber \end{align} \tag{ 94 }$$

with the corresponding $\omega_{\pm}$ at $T>T_{\rm c}$ or $T<T_{\rm c}$.

The FDT in the frequency domain ${\rm{Im}} \hat R(\omega) = - \beta \omega \hat C(\omega)$ implies that $\hat C(\omega)$ should vanish in the same frequency intervals in which the linear response is real. In the $\omega=0$ case the linear response is real and the FDT as written above only imposes that $\hat C(\omega=0)=\int_{-\infty}^\infty {\rm d}t \, C(t)$ must be finite. One has to bear in mind that in the cases in which the correlation function approaches a non-vanishing constant q asymptotically the Fourier transform to be computed is the one of $C_{\rm st}(t_1-t_2) = C(t_1-t_2)-q$ with respect to $t_1-t_2$ and the integral should then yield $C_{\rm st}(\omega=0)=\int_{-\infty}^\infty {\rm d}t \, C_{\rm st}(t)= 1-q$; more details are given in appendix A.

Consistently with the FDT constraints discussed in the previous paragraph, the time-dependent equation for $C_{\rm st}(t_1-t_2) = C(t_1-t_2)-q_{\rm in}$ can be treated following the steps explained in [9] and appendix A.1.3 (that do not assume FDT) and one finds

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C_{\rm st}(\omega)=J^2 C_{\rm st}(\omega) |\hat R(\omega)|^2 . \nonumber \end{align} \tag{ 95 }$$

This equation has two solutions, either $C_{\rm st}(\omega)=0$ or $|\hat R(\omega)|^2 = 1/J^2$. The latter holds for frequencies in the interval $\omega \in [\omega_-, \omega_+]$. Outside of this interval $C_{\rm st}(\omega)$ must vanish.

By taking the derivative of equation (79) with respect to t₂ one readily checks that it equals equation (78) if the FDT between R and C is satisfied for all times and $\partial_{t_2}C_2(t_2, 0)=0$ implying

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C_2(t_2,0)=q_{\rm in} . \nonumber \end{align} \tag{ 96 }$$

This last condition is a property of equilibrium as we have already discussed.

Having established that C₂ is a constant, equation (80) enforces $q_{\rm in}=0$, the high temperature initial value, or

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z_f = \frac{J^2}{T'} [1-C(t_1,0)] + \frac{J^2}{T'} [C(t_1,0) + (1-2q_{\rm in})] = 2J , \nonumber \end{align} \tag{ 97 }$$

the low temperature Lagrange multiplier.

Concerning the correlation function C(t₁,0), we write it as $C(t_1, 0)=C_{\rm st}(t_1, 0)+q_0$ allowing for a non-vanishing asymptotic value, q₀, and taking $C_{\rm st}(t_1, 0)$ such that it vanishes in the long t₁ limit. The equation (79) evaluated at t₂  =  0 is then rewritten as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left(m\partial_{t_{1}}^{2}+z_f\right)[C_{\rm st}(t_{1},0) + q_0] = J^2\int_{0}^{t_{1}}{\rm d}t' \; R(t_{1},t') [C_{\rm st}(t',0)+q_0] +\frac{J^2}{T'} \; [C_{\rm st}(t_{1},0) +q_0 - q_{\rm in}^2]. \label{eq:dyn-eqs-Cbis} \nonumber \end{align} \tag{ 98 }$$

This equation has three terms that do not depend on $C_{\rm st}(t_1, 0)$ explicitly

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z_f q_0 - J^2 q_0 \lim_{t_1\to\infty} \int_0^{t_1} {\rm d}t' \; R(t_1,t') - \frac{J^2}{T'} (q_0 - q_{\rm in}^2) \nonumber \end{align} \tag{ 99 }$$

and their sum should vanish in the long t₁ limit. It trivially does for high temperature initial states since $q_0=q_{\rm in}=0$ and, in equilibrium at low temperatures, we can assume $q_0=q_{\rm in}$, use FDT, and find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z_f =\frac{2J^2}{T'} (1-q_{\rm in}) =2J \nonumber \end{align} \tag{ 100 }$$

confirming the assumption $q_0=q_{\rm in}$. The remaining equation fixes the time-dependence of C(t₁,0).

Let us assume that the limiting value of the Lagrange multiplier is a constant

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \lim_{t_1\to\infty} z(t_1) = z_f. \nonumber \end{align} \tag{ 111 }$$

The limit of the correlation function can be zero if the system decorrelates completely at sufficiently long time-delays or non-zero if it remains within a confined state. We therefore call

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle q = \lim_{t_1-t_2\to\infty} \lim_{t_2\to \infty} C(t_1,t_2) \nonumber \end{align} \tag{ 112 }$$

the asymptotic value of the full two time correlation, or its stationary part in possible ageing cases. Similarly,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle q_0 = \lim_{t_1\to \infty} C(t_1,0), \nonumber \end{align} \tag{ 113 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle q_2 = \lim_{t_1\to\infty}C_2(t_1,0). \nonumber \end{align} \tag{ 114 }$$

The equation for the linear response function does not depend explicitly on the pre-quench parameters, it only does on the post-quench mass m and interaction strength J. (An implicit dependence on the initial state is not excluded, since the value taken by the Lagrange multiplier may depend on it.) The analysis that we developed for the constant energy dynamics applies to the sudden quench case too. The solution of the response equation in Fourier space yields

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat R(\omega)=\frac{1}{2 J^2} \left[ (-m\omega^2 + z_f) \pm \sqrt{(-m\omega^2 + z_f)^2 - 4 J^2} \right] \nonumber \end{align} \tag{ 115 }$$

with the zero frequency value

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat R(\omega=0) = {\frac{1}{2 J^2}} \left(z_f - \sqrt{z_f^2 - 4 J^2} \right). \label{eq:Romega00} \nonumber \end{align} \tag{ 116 }$$

From the numerical solution of the full equations that we will present in section 6 we infer that for m  =  m₀ and $J\neq J_0$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z_f = \left\{\begin{array}{@{}ll@{}} 2J \qquad & \qquad x>y , \nonumber \\ T' + J^2/T' \qquad & \qquad x<y , \end{array} \right. \label{eq:prediction-zf} \nonumber \end{align} \tag{ 117 }$$

with, we recall, x  =  J/J₀ and $y=T'/J_0$. Replacing in equation (116) one notices that the minus sign has to be selected for x  <  y at low frequency and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi_{\rm st} \equiv \hat R(\omega=0) = \left\{\begin{array}{@{}ll@{}} 1/T' & \qquad x<y , \nonumber \\ 1/J & \qquad x>y . \end{array} \right. \label{eq:prediction-Romega} \nonumber \end{align} \tag{ 118 }$$

After the quench Im$\hat R(\omega)$ is non-zero in a finite interval of frequencies $[\omega_-, \omega_+]$ with $m\omega^2_{\pm} = (z_f\pm 2J)$ as in equation (87) and z_f taking the values in equation (117). Therefore, Im$\hat R(\omega)$ 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 equation (117). We have verified them with the complete numerical solution of the $N\to\infty$ dynamic equations, see section 6, on all sectors of the phase diagram. We can now Fourier back to real time to get the full time dependence of the linear response function. In the numerical section we will compare this functional form, named 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,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \overline z = 2 \overline e_{\rm kin} - 2 \overline e_{\rm pot} , \qquad \qquad e_{\rm tot}(0^+) = \overline e_{\rm kin} + \overline e_{\rm pot} , \nonumber \end{align} \tag{ 119 }$$

where the overlines represent a long-time average defined in equation (54).

Using these relations one derives the parameter dependence of the kinetic and potential energies in the four relevant regions of the phase diagram parametrised by x  =  J/J₀ and $y=T'/J_0$.

• (I)  
  y  >  1 and x  <  y ($q_{\rm in}=0$, $z_f=T'+J^2/T'$, $\chi_{\rm st}=1/T'$)

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}rcl@{}} 4{\overline e_{\rm kin}} & = & {- \frac{JJ_0}{T'} + 2T' + \frac{J^2}{T'} \ = \ J_0 \left(- \frac{x}{y} + 2y + \frac{x^2}{y}\right) } \nonumber \\ 4{\overline e_{\rm pot}} & = & {- \frac{JJ_0}{T'} - \frac{J^2}{T'} \ = \ J_0 \left(- \frac{x}{y} - \frac{x^2}{y}\right) } \label{eq:I} \end{array} \nonumber \end{align} \tag{ 120 }$$
• (II)  
  y  >  1 and x  >  y ($q_{\rm in}=0$, z_f  =  2J, $\chi_{\rm st}=1/J$)

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}rcl@{}} 4{\overline e_{\rm kin}} & = & {- \frac{JJ_0}{T'} + T' + 2J \ = \ J_0 \left(-\frac{x}{y} + y + 2 x\right) } \nonumber \\ 4{\overline e_{\rm pot}} & = & {- \frac{JJ_0}{T'} + T' - 2J \ = \ J_0 \left(- \frac{x}{y} + y - 2 x\right) } \end{array} \nonumber \end{align} \tag{ 121 }$$
• (III)  
  y  <  1 and x  >  y ($q_{\rm in}\neq 0$, z_f  =  2J, $\chi_{\rm st}=1/J$)

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}rcl@{}} 4{\overline e_{\rm kin}} & = & {T' + \frac{T'J}{J_0} \ = \ J_0 \left(y + yx \right) } \nonumber \\ 4{\overline e_{\rm pot}} & = & {T' + \frac{T'J}{J_0} -4J \ = \ J_0 \left(y + yx - 4x\right) } \end{array} \label{eq:asymptotic-kin-pot-energy-III} \nonumber \end{align} \tag{ 122 }$$
• (IV)  
  y  <  1 and x  <  y ($q_{\rm in} \neq 0$, $z_f = T'+J^2/T'$, $\chi_{\rm st}=1/T'$)

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}rcl@{}} 4{\overline e_{\rm kin}} & = & {\frac{JT'}{J_0} - 2J + 2 T' + \frac{J^2}{T'} \ = \ J_0 \left(xy - 2 x + 2 y + \frac{x^2}{y} \right) } \nonumber \\ 4{\overline e_{\rm pot}} & = & {\frac{JT'}{J_0} - 2J - \frac{J^2}{T'} \ = \ J_0 \left(xy - 2 x - \frac{x^2}{y} \right) } \label{eq:IV} \end{array} \nonumber \end{align} \tag{ 123 }$$

The minimum potential energy density, $e_{\rm pot}=-J$ is realised for $T'=0$ in Sector III.

We stress that we have found these results without using FDT and they can therefore hold out of thermal equilibrium. We will investigate later which other conditions impose the use of FDT, a strong Gibbs–Boltzmann equilibrium condition.

We can identify a kinetic temperature from the kinetic energy densities derived in the previous subsection by simply imposing $T_{\rm kin} = 2\overline{e}_{\rm kin}$. This operation leads to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}rclrcl@{}} {T^{({\rm I})}_{\rm kin}} &=& {- \frac{JJ_0}{2T'} + T' + \frac{J^2}{2T'}} , \qquad \qquad \qquad \; T^{({\rm II})}_{\rm kin} &=& {- \frac{JJ_0}{2T'} + \frac{T'}{2} + J } , \nonumber \\ {T^{({\rm III})}_{\rm kin}} &=& {\frac{T'}{2} \left(1+\frac{J}{J_0}\right) } , \qquad \qquad \qquad \qquad {T^{({\rm IV})}_{\rm kin}} &=& {\frac{JT'}{2J_0} - J + T' + \frac{J^2}{2T'} } , \end{array} \label{eq:kin-temps} \nonumber \end{align} \tag{ 124 }$$

in the four Sectors of the phase diagram distinguished in the previous subsection.

In appendix A.3.1 we explain how we can exploit the conservation of the total energy, under the assumption that the asymptotic kinetic and potential energies take the form of Gibbs–Boltzmann equilibrium paramagnetic and condensed equilibrium phases at a single temperature T_f to fix its value. This means that we require $2\overline e_{\rm kin} = T_f$ and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \overline{e}_{\rm pot} &= \left\{\begin{array}{@{}lll@{}} {-\frac{J^2}{2T_f}} \qquad \qquad & & {\rm{paramagnetic~or}} \nonumber \\ {-\frac{J^2}{2T_f} (1-q^2)} \qquad \qquad & & {\rm{condensed~or~two-step}} , \end{array} \right. \nonumber \end{align} \tag{ 125 }$$

(the 'two-step' refers to the ageing Ansatz that will be discussed in section 4.4.7) with q  =  1  −  T_f/J. For x  >  y we find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{T_f}{J_0} &= \left\{\begin{array}{@{}ll@{}} {\frac{J^2}{J_0} \; \left[ J -\frac{T'}{2} - \frac{JJ_0}{2T'} \right]^{-1} } \qquad \qquad & {\rm{for}} \;\; 1\leqslant y \;\; {\rm{(IIa)}} , \nonumber \\ {\frac{J}{J_0} +\frac{T'}{2J_0} - \frac{J}{2T'} =x +\frac{y}{2} - \frac{x}{2y} } \qquad \qquad & {\rm{for}} \;\; 1\leqslant y \;\; {\rm{(IIb)}} , \nonumber \\ {\frac{T'}{2J_0} \left(1+\frac{J}{J_0}\right) = \frac{y}{2} \left(1+x \right) } \qquad \qquad & {\rm{for}} \;\; y\leqslant 1 \;\; {\rm{(III)}} . \end{array} \right. \label{eq:Tf-asymptotic-values} \nonumber \end{align} \tag{ 126 }$$

The roman numbers between parenthesis refer to the cases listed in equations (120)–(123) and to the four sectors in the phase diagram in figure 5. The difference between (IIa) and (IIb) is that in the first case we used a paramagnetic potential energy and in the second case a condensed one at T_f. The values for x  <  y are given in appendix A.3.1. One can check that in the cases (IIb) and (III) $T_f=T_{\rm kin}$, the kinetic temperatures given in equation (124). Instead, final and kinetic temperatures do not coincide in (IIa) nor in (I) and (IV). This fact excludes the possibility of reaching 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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 0 \leqslant T_f \leqslant J \nonumber \end{align} \tag{ 127 }$$

serve to find special curves on the phase diagram. The first bound is natural since we do not want to have a negative kinetic energy. The second one ensures that $q\geqslant 0$. The implications of these bounds, that are examined in appendix A.3.2, are

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y &\leqslant \left\{\begin{array}{@{}ll@{}} \sqrt{x} & \qquad \qquad {\rm{for}} \;\; y\geqslant 1 , \nonumber \\ {\frac{2x}{x+1}} & \qquad \qquad {\rm{for}} \;\; y\leqslant 1 . \end{array} \right. \nonumber \end{align} \tag{ 128 }$$

They mean that an asymptotic state with a single temperature T_f, or a two-step regime with temperature T_f for $C:1 \to q$ and $T_{\rm ef\,\!f}\to\infty$ for $C:q\to 0$ (like the one explained in section 4.4.7) could only exist below the piecewise curve $y(x)$. One can simply check that the piece for $y\leqslant 1$ lies in Sector IV, since y  <  x, and it is therefore irrelevant given that we have already shown that there cannot be a single temperature scenario in this Sector. The limit then moves to x  =  y for y  <  1. Parameters on the curve $y=\sqrt{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\to\infty$ limit seems possible for y  <  1 and x  >  y (III), and below the curve $y=\sqrt{x}$ for y  >  1 in II. It is instructive to work out the limiting values of $T_f/J_0$ and q  =  1  −  T_f/J on the borders of the region y  <  1 and y  <  x (III) of the phase diagram, and the no-quench case x  =  1:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{T_f}{J_0} \ = \ \left\{\begin{array}{@{}ll@{}} 0 & \qquad y=0 \nonumber \\ {\frac{x(x+1)}{2} }& \qquad x =y \nonumber \\ {\frac{x+y}{2} } & \qquad y=1 \nonumber \\ y & \qquad x=1 \end{array} \right. \qquad & \qquad q \ = \ \left\{\begin{array}{@{}ll@{}} 1 & \qquad y=0 \nonumber \\ {\frac{1-x}{2} }& \qquad x =y \nonumber \\ {\frac{x-1}{2x} } & \qquad y=1 \nonumber \\ 1-y & \qquad x=1. \end{array} \right. \qquad \quad {\rm{(III)}} \nonumber \end{align} \tag{ 129 }$$

These values match at x  =  y  =  1. q is larger than zero on the lines x  =  y and y  =  1 that mark the end of what we call Sector III. Moreover, the approximate asymptotic analysis of the mode dynamics of the finite N system will lead to $T_\mu=T_f$ in Sector III, see equation (186).

On the limiting curve in Sector II, $y=\sqrt{x}$ for y  >  1, T_f  =  J and q  =  0.

As one could have intuitively expected, $T_f>T'$ for energy injection quenches (x  <  1) and $T_f<T'$ for energy extraction quenches (x  >  1).

We now add one assumption to the analysis: that the FDT, at the single temperature T_f, relates the linear response to the correlation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle R(t_1-t_2) = - \frac{1}{T_f} \; \frac{dC(t_1-t_2)}{d(t_1-t_2)} \ \theta(t_1-t_2)\label{eq:fdt}. \nonumber \end{align} \tag{ 130 }$$

Static susceptibility.

The values of the zero frequency linear response

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat R(\omega = 0) = \int_0^\infty {\rm d}\tau \; R(\tau) = \frac{1}{T_f} (1-q) , \nonumber \end{align} \tag{ 131 }$$

where one must recall that τ is a time-difference, force

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat R(\omega = 0) = \left\{\begin{array}{@{}lllll@{}} 1/T' \quad & {\rm{for}} \qquad x<y \qquad & \Rightarrow \qquad T_f=T' \quad & {\rm{if}} \quad q=0 \quad {\rm{(I)}} \; \& \; {\rm{(IV)}} \nonumber \\ 1/J \quad & {\rm{for}}\qquad x>y \qquad & \Rightarrow \qquad T_f=J \quad & {\rm{if}} \quad q=0 \quad {\rm{(IIa)}} \nonumber \\ 1/J \quad & {\rm{for}}\qquad x>y \qquad & \Rightarrow \qquad T_f=J/(1-q)\quad & {\rm{if}} \quad q\neq 0 \quad {\rm{(IIb)}} \; \& \; {\rm{(III)}}. \end{array} \right. \nonumber \end{align} \tag{ 132 }$$

The result in the last line, valid for (IIb) and (III), does not put any additional constraint on T_f. Instead, the conditions in the first two lines are incompatible with the expressions imposed by the energetic considerations. They corroborate the impossibility of having a single T_f in (I) and (IV) or with q  =  0 in (II).

Limits of validity of the single temperature scenario.

One can also look for a two step solution with similar characteristics to the ageing one found for dissipative dynamics [44] and summarised in section 2.3. Asking for the relation between correlation and linear response in equation (130) to hold in a stationary regime of relaxation in which C decays from 1 to q, and that the effective temperature T_eff characterising the second regime of decay from q to 0 is infinity, as in the dissipative case [44], one recovers

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}lll@{}} z_f = 2 J , \qquad \chi_{\rm st} = {\frac{1}{J}} , \qquad & {q = 1- \frac{T_f}{J}} \quad \qquad \;\; {\rm{and}} \qquad & q_0 = q_2=0 \end{array} \nonumber \end{align} \tag{ 133 }$$

with the same T_f and $q\neq 0$ as in equation (126). This kind of ageing Ansatz takes the form of the one realised in the dissipative dynamics [44] and in the p  =  3 isolated dynamics [7], but this regime is not realised by the dynamics of the p  =  2 isolated model as the numerical results of the following sections show.

From the analysis of the equation ruling the two-time correlation function, assuming stationarity and hence a dependence on time-difference only even after a quench, one deduces (the details of the derivation are given in appendix A.1.3, see also [9] for a general treatment)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat C_{\rm st}(\omega) = J^2 \hat C_{\rm st}(\omega) |\hat R(\omega)|^2 \label{eq:SD-corr} \nonumber \end{align} \tag{ 134 }$$

where $\hat C_{\rm st}(\omega)$ is the Fourier transform of the time-varying part (subtracting the possible non-vanishing asymptotic value q). Note that the relation between the correlation and the linear response is the same as the one that we derived in the constant energy no-quench problem. It implies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat C_{\rm st}(\omega) \neq 0 \;\; {\rm{and}} \;\; |\hat R(\omega)|^2 = 1/J^2 \qquad {\rm{or}} \qquad \hat C_{\rm st}(\omega) = 0 , \nonumber \end{align} \tag{ 135 }$$

independently of the control parameters. Below we check numerically that these relations are satisfied in various quenches. In particular, from the analytic form of $\hat R(\omega)$ one can easily see that $|\hat R(\omega)|^2 = 1/J^2$ in the frequency interval in which the linear response is complex.

Importantly enough, we cannot yield an explicit analytic form of $\hat C(\omega)$ since it is factorised on both sides of the identity (134). We are forced to go back to the full set of dynamic equations and solve them numerically to get insight on the behaviour of $C(t_1, t_2)$.

We will see that a state with a single T_f characterising the fluctuation-dissipation relation is reached numerically in the following two cases only:

• (1)  
  the dynamics are run at constant parameters (no quench), x  =  1;
• (2)  
  the special relation $y=\sqrt{x}$ between parameters holds and y  >  1 (within sector II).

In all other cases no equilibrium results à la Gibbs–Boltzmann are found for the global observables (correlation function, linear response function, kinetic and potential energies) but a different statistical description, of a generalised kind, should be adopted. In particular, in Sector III where the conditions derived from the energy conservation and static susceptibility allow for a single temperature scenario, the full solution of the complete set of questions will prove that this is not realised. A detailed explanation is given in the numerical section of the paper and in the analysis of the N  −  1 integrals of motion presented in section 7.

We recall that the $p\geqslant 3$ strongly interacting case behaves very differently [7]. On the one hand, equilibrium towards a proper paramagnetic state, and within confining metastable states, were reached in two sectors of its dynamic phase diagram. On the other hand, an ageing asymptotic state in a tuned regime of parameters was also found for more than two spin interactions in the potential energy. In the p  =  2 integrable model we do not find an ageing asymptotic state. Moreover, Gibbs–Boltzmann equilibrium is achieved in the two very particular cases listed above only.

We begin with the kinetic energy right before the quench,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {e}^{(N)}_{{\rm kin}}(0^-)=\frac{m}{2N}\sum_{\mu=1}^{N}\langle \, \dot s^{2}_{\mu} \, \rangle_{{\rm eq}}=\frac{T'}{2} . \nonumber \end{align} \tag{ 172 }$$

It coincides with the infinite-N mean-field result.

Next we analyze the pre-quench potential energy,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:e_pot_ic_finite_n} {e}^{(N)}_{\rm pot}(0^-) = - \frac{1}{2N} \sum^{N}_{\mu=1} \lambda^{(0)}_\mu \langle \, s_\mu^2 \, \rangle_{\rm eq} = -\frac{T'}{2N}\left(\sum^{N}_{\mu=1}\frac{\lambda^{(0)}_{\mu}-z^{(N)}_{{\rm eq}}}{z^{(N)}_{{\rm eq}}-\lambda^{(0)}_{\mu}} + z^{(N)}_{{\rm eq}}\sum^{N}_{\mu=1}\frac{1}{z^{(N)}_{{\rm eq}}-\lambda^{(0)}_{\mu}} \right) = \frac{T^{\prime}}{2}-\frac{z^{(N)}_{{\rm eq}}}{2} , \nonumber \end{align} \tag{ 173 }$$

where we used equation (166). In the equilibrium condensed phase we can rely on $z^{(N)}_{{\rm eq}}=\lambda^{(0)}_{{\rm max}}+T^{\prime}/(N q^{(N)}_{{\rm in}})$ to obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:e_pot_finite_n} {e}^{(N)}_{{\rm pot}}(0^-)=-\frac{\lambda^{(0)}_{{\rm max}}}{2}+\frac{T^{\prime}}{2}-\frac{T^{\prime}}{2Nq^{(N)}_{{\rm in}}} . \nonumber \end{align} \tag{ 174 }$$

The $N\to\infty$ result is $e_{{\rm pot}}(0^-)=-J_0+T'/2$, consistent with equation (174), since $\lim_{N\rightarrow\infty}\lambda^{(N)}_{{\rm max}}=2J_0$.

In the paramagnetic phase $q^{(N)}_{{\rm in}}\ll 1$ and the third term in equation (174) induces important corrections. In this case, using equation (169) we can write

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {e}^{(N)}_{{\rm pot}}(0^-)\simeq -\frac{J^2_0}{2T^{\prime}}-\left(\frac{\lambda^{(0)}_{{\rm max}}}{2}-J_0\right)-\frac{s(T^{\prime})}{2N} \nonumber \end{align} \tag{ 175 }$$

and one readily recovers the $N\to\infty$ limit $e_{{\rm pot}}(0^-)=-J_0^2/(2T^{\prime})$.

Now we will compute the values of the kinetic and potential energy, and the Lagrange multiplier after an interaction quench

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \lambda^{(0)}_{\mu}\rightarrow\lambda_{\mu}=\frac{J}{J_0}\lambda^{(0)}_{\mu}. \nonumber \end{align} \tag{ 176 }$$

The kinetic energy is not affected by the quench in the interaction and, just as in the $N\to\infty$ limit (see section 4.3), we have that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {e}^{(N)}_{{\rm kin}}(0^+)={e}^{(N)}_{{\rm kin}}(0^-)=\frac{T'}{2}. \nonumber \end{align} \tag{ 177 }$$

For the potential energy it is enough to note that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {e}^{(N)}_{{\rm pot}}(0^+)=-\frac{1}{2N}\sum^{N}_{\mu=1}\lambda_{\mu}\langle \, s^2_{\mu} \, \rangle_{{\rm eq}} = -\frac{1}{2N}\frac{J}{J_0}\sum^{N}_{\mu=1}\lambda^{(0)}_{\mu}\langle \, s^2_{\mu} \, \rangle_{{\rm eq}} = \frac{J}{J_0}{e}^{(N)}_{{\rm pot}}(0^-). \nonumber \end{align} \tag{ 178 }$$

Using that $z^{(N)}(0^+)=2({e}^{(N)}_{{\rm kin}}(0^+)-{e}^{(N)}_{{\rm pot}}(0^+)) = T'-2J/J_0 \, e_{\rm pot}^{(N)}(0^-)$ it is now easy to find the initial value of the Lagrange multiplier. When the initial conditions are taken from the condensed phase, $q^{(N)}_{{\rm in}}=\mathcal{O}(1)$, and we can write

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z^{(N)}(0^+)=\lambda_{{\rm max}}+T^{\prime}\left(1-\frac{J}{J_0} \right) + \frac{JT^{\prime}}{NJ_0q^{(N)}_{{\rm in}}}. \nonumber \end{align} \tag{ 179 }$$

For initial states in the paramagnetic phase

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z^{(N)}(0^+)\simeq \frac{JJ_0}{T^{\prime}} + T^{\prime}+\left(\lambda_{{\rm max}}-2J\right)+\frac{J}{J_0}\frac{s(T^{\prime})}{N}. \nonumber \end{align} \tag{ 180 }$$

In this section we use parameters in the paramagnetic (PM) equilibrium phase. We fix $J_0=m_0=m=1$, we equilibrate the initial conditions at $T^{\prime}=1.25$ and we evolve with parameters J  =  J₀  =  1 and m  =  m₀  =  1 (no quench). In figures 7 and 8 we show the numerical results for infinite N and finite N, respectively. The system should remain in Gibbs–Boltzmann equilibrium in the PM phase in both cases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 7 we analyse the results of the numerical integration of the $N\to\infty$ Schwinger–Dyson equations for $T'=1.25$ and all other parameters set to 1. The figure shows the time evolution of the correlation function (figure 7(a)), the fluctuation dissipation parametric plot (figure 7(b)), the response function, $R(t_1, t_2)$ together with $-(1/T^{\prime}) \, \partial_{t_1}C(t_1, t_2)$ (figure 7(c)), the deviation of the Lagrange multiplier from its equilibrium value (figure 7(d)), the potential and kinetic contributions to the energy density (figure 7(e)), and two studies of the fluctuation-dissipation theorem in the frequency domain (figures 7(f) and (g)).

The correlation with the initial condition, C(t₁,0), and the ones between two different times, $C(t_1, t_2)$, are identical, 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 figure 8(a) which were obtained by integrating Newton equations for finite N.

The fluctuation-dissipation relation [68] is satisfied with the temperature of the initial condition, that is the same as the one of the final state. This fact can be proven in general for Newtonian evolution of initial configurations drawn from Gibbs–Boltzmann equilibrium [85]. Indeed in figure 7(b) we display the parametric plot $\chi(t_1, t_2)=\int_{t_2}^{t_1} {\rm d}t' \, R(t_1, t')$ versus $C(t_1, t_2)$ for two waiting times t₂ and, with a dashed line, the equilibrium result $-1/T' \; C$ finding perfect agreement within our numerical accuracy. As a further confirmation of the validity of the fluctuation-dissipation theorem, we show the response function, $R(t_1, t_2)$, and $-(1/T^{\prime}) \, \partial_{t_1}C(t_1, t_2)$, both plotted against $t_1-t_2$, for two different values of t₂ in figure 7(c). The two quantities coincide almost perfectly. In the same panel we checked that the response function coincides with the one derived by taking the inverse of the theoretical Fourier transform of the stationary asymptotic response, given by equation (86). The theoretical curve is indicated as $R_{{\rm st}}$.

The Lagrange multiplier and the potential and kinetic energies remain constant throughout the evolution of the system and equal to their predicted values, apart from small deviations due to the numerical errors introduced by the numerical integration scheme, see panels (d) and (f) of figure 7, and appendix C. The plot showing $z(t)$ proves that the relative error in this quantity is at most of order 10⁻⁷. All these results are compatible with Gibbs–Boltzmann equilibrium in the paramagnetic phase. We do not show the time evolution of the off-diagonal correlation with the initial configuration, C₂(t,0), since it is identically zero at all times.

In figure 7(f) we show the Fourier transforms of the correlation and response functions, $\hat{C}(\omega, t_2)$ and $\hat{R}(\omega, t_2)$ respectively (the transform is performed on the variable $\tau=t_1-t_2$ with t₂ fixed), for two different values of t₂. Note that we are showing only the real part of $\hat{C}(\omega, t_2)$, since we implicitly assume that $C(t_2+\tau, t_2)=C(t_2-\tau, t_2)$. The black solid lines represent the theoretical prediction for the real and imaginary parts of the Fourier transform of the response function in the stationary regime, $\hat{R}_{{\rm st}}(\omega)$, given by equation (86), the inverse Fourier transform of which is plotted in figure 7(c). In panel (g) of figure 7, the ratio $-{\rm Im}{\hat{R}(\omega)}/(\omega \hat{C}(\omega))$ together with the prediction $1/T^{\prime}$ from FDT indicated by a horizontal dashed line are shown. Note the deviation from the flat result at the right edge of the frequency spectrum. This is due to the fact that the ratio approaches zero over zero and the numerical error incurred for those large frequencies is much amplified. At the left end of the spectrum, the more interesting low frequency regime, the oscillations are only present for the t  =  0 curve.

In figure 8 we show results obtained by solving the dynamics of each mode in a finite N system with the method explained in section 5.8. In figure 8(a) we see that the correlations with the initial condition quickly relax to 0, as expected in the PM phase. They do with a weak size-dependence in the long time-delay tails. We only show the correlation with the initial configuration since we have checked that the time-delayed one is stationary. Also included in this panel is the same correlation function computed using the Schwinger–Dyson equation valid in the $N\to\infty$ limit. We see perfect agreement with the finite N results at short times and small deviations at longer times. In figure 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) of figure 8 where we plot it subtracting $z_{\rm eq}^{(N)}$, calculated as the solution to equation (166) (a non-linear equation) that in the $N\rightarrow\infty$ limit yields $z_{\rm eq}^{(N)}\rightarrow T'+ J^2/T'$. The very weak (oscillatory) deviation from $z_{\rm eq}^{(N)}$ decreases with the size of the time-step used in the numerical solution of the dynamic equations (in the $N\to\infty$ case we have a similar effect, see appendix C).

The mode-by-mode analysis of the finite N dynamics is performed in panels (d) and (e) of figure 8. These panels display the time-delay dependence of the correlation function of the first and last mode. In figure 8(f) we display the mode temperatures $T_{\mu}(t)$ at the initial time and after a long time evolution. The mode temperatures coincide with the expected equilibrium value, except for the largest modes, where there is a very small deviation. These variations represent small numerical errors due to the finite time-step discretisation used numerically, and are hard to improve algorithmically unless by using a still smaller integration step. We have also checked (not shown) that the mode correlations $C_{\mu}(t_1, t_2)$ and the mode response function $R_{\mu}(t_1, t_2)$ satisfy the fluctuation-dissipation relation with a temperature given by $T^{\prime}$ for all modes.

We now turn to the constant energy dynamics in the condensed, low temperature equilibrium phase.

In figure 9 we show the mode dynamics for initial conditions in equilibrium at $T'=0.5$. From figure 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 figure 9(c) we show the Lagrange multiplier, which should be constant in equilibrium. In the numerical solution, the Lagrange multiplier exhibits oscillations around the initial value. Their amplitude decreases consistently with the integration step δ, implying that for $\delta\rightarrow0$ we recover the expected constant behaviour. The two-time global correlation $C(t_1, t_2)$ is stationary for all system sizes (not shown), so we focus on the particular case with t₂  =  0. We can see from figure 9(a) that, at variance with the paramagnetic case, the dynamics of the correlation function C(t₁,0) has a strong dependence on the system size. After a fast decay from the initial value, the correlation shows a plateau, the lifetime of which increases with system size, approaching asymptotically the value predicted by the $N\to\infty$ treatment that is shown with a (red) horizontal line. The source of this size dependence is the behaviour of the largest mode $\mu=N$. In figure 9(f) we show the time dependence of the largest mode correlation function C_N(t,0) for different system sizes. We observe that its oscillation frequency decreases as we increase the system size. A similar finite size effect is seen in the dynamics of the next-to-largest mode in panel (g). Since the largest modes dominate the long-time dynamics, this effect causes the size dependence of the plateau lifetime. For $N\rightarrow \infty$ the oscillation frequency of the Nth mode goes to zero, allowing for the presence of an infinite plateau, see figure 11. The modes lying in the middle and other end of the spectrum have almost no size dependence, as shown in (d).

Zoom In Zoom Out Reset image size

Figure 10 investigates the mode kinetic and potential energies and the mode temperatures that can be extracted from them. In figure 10(c) we show the mode temperatures $T_{\mu}(t)$ at two measuring times, as a function of the mode index, what we will call temperature spectrum later. We observe deviations from the expected behaviour $T^{\rm tot}_{\mu}=T'\;\forall\mu$ only close to $\mu=N$. To gain insight into the mode temperature deviations, in panels (a) and (b) we show the time evolution of the kinetic, potential and total energies for the first and last modes. For the lowest modes, the kinetic and potential energies are constant, equal, and their sum is identical to $T'$, consistently with equilibrium dynamics. However, for the modes that are closer to $\mu=N$ the energies weakly oscillate with an amplitude that decreases with the integration step and should vanish for $\delta \to 0$. This implies that in such limit, all mode temperatures should be equal to the equilibrium temperature, even for modes close to the right edge of the spectrum.

Zoom In Zoom Out Reset image size

We now turn to the analysis of the dynamics in the $N\to\infty$ limit. In figure 11 we show the results obtained through the numerical integration of the Schwinger–Dyson equations for the two-time correlation and response functions, and the same choice of parameters, that is $T^{\prime}=0.5$ and J  =  1. Stationarity is satisfied as well as the FDT with the initial temperature, see (b) and (c). The main difference with the case in figure 7 is that the correlation functions, both with the initial condition and with the configuration at a waiting-time t₂, relax to a non-vanishing value (a). Within numerical accuracy we observe $q_0 = q \simeq 1-T^{\prime}/J = 0.5$ and this value as well as the potential and kinetic energies (not shown) are consistent with equilibrium at $T'$. Also in this case, we checked that the response $R(t_2+\tau, t_2)$, for $t_2 \geqslant 0$, coincides with the (numerical) inverse Fourier transform of the theoretical prediction given by equation (86) for the Fourier transform of the stationary asymptotic R, see panel (c).

Zoom In Zoom Out Reset image size

In panel (d) we show the time-evolution of the off-diagonal correlation, C₂(t,0). In the case of equilibrium dynamics, C₂(t,0) should be constant and equal to $q_{{\rm in}}$. As one can see, the value of C₂(t,0) obtained by numerical integration is not exactly a constant function, but it approaches a constant in the long time limit which differs from $q=q_{{\rm in}}=0.5$ by a very small amount. This deviation is only due to the approximations introduced by the numerical integration scheme. The same can be said about the behaviour of the numerical $z(t)$, see panel (e). Its value oscillates around the expected equilibrium value $z_{{\rm eq}}=2 J = 2$ at $T^{\prime}=0.5$, with oscillations amplitude of order 10⁻⁷ 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 figure 11(f). Again, note that we are showing only the real part of $\hat{C}(\omega, t_2)$, since we implicitly assume that $C(t_2+\tau, t_2)=C(t_2-\tau, t_2)$. The black solid lines represent the theoretical prediction for the real and imaginary parts of the Fourier transform of the response function in the stationary regime, $\hat{R}_{{\rm st}}(\omega)$, given by equation (86), the inverse transform of which is plotted in (c). In panel (g) we display the ratio $-{\rm Im}{\hat{R}(\omega)}/(\omega \hat{C}(\omega))$ together with the prediction $1/T^{\prime}$ from FDT indicated by a dashed horizontal line. In all presentations we find good agreement with the validity of FDT with the proviso that in the plot in (g) the high frequency regime is contaminated by the numerical error, and the low frequency regime by the fact that we can perform the Fourier transform on a finite time window only, and this causes the dependence on t₂ shown in the plot.

In this subsection we summarise the dynamics of a system prepared in equilibrium in the paramagnetic state y  >  1 and quenched to a value of J such that y  >  x. This is what we called Sector I in the phase diagram. Two cases need to be further distinguished within this Sector. For x  <  1 energy is injected in the quench and this problem is treated in figures 12 and 13. For y  >  x  >  1, instead, a small amount of energy is extracted from the system and we explain the difference that this implies in figure 14 and the discussion around it.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The first observation is that we confirm that, in both cases x  >  1 and x  <  1, $\hat R(\omega=0)=1/T'$ and $z_f = T' + J^2/T'$. The latter claim can be verified in figure 12(d) for x  >  1 and $N\to\infty$, for example. Next we check that the dynamics approach a stationary asymptotic state in which the global one-time quantities are constant, as seen for example in panels (d) and (e) in the same figure, and the two-time correlation and linear response depend on the time delay only, see panels (a) and (c). The last question concerns the fluctuation dissipation relation between linear response and correlation function. We used the parametric plot (b) to demonstrate that the evolution does not approach a state characterised by a single temperature but, instead, $\chi(C)$ is curved and not even single valued. Having said this, the comparison of the time-dependent $R(t_1, t_2)$ and $-\partial_{t_1} C(t_1, t_2)$ in panel (c) could have fooled us into the belief that the FDT holds with a single temperature. Indeed, the difference between the two functions is not visible in this scale. In panel (f) we show the Fourier transform of the correlation function and linear response function. The Fourier transform of the linear response function is consistent with the analytic prediction presented in section 4.4.2, in particular, the imaginary part is non-vanishing only in the interval $[ \omega_{-}, \omega_{+}]$ (save some small oscillations around zero due to numerical error outside this interval), with $m \omega^2_{\pm} = z_f \pm 2 J$. The fluctuation-dissipation ratio is shown in (g) and we will come back to its analysis below, where we present the mode dependent results for finite N.

We next study the evolution under the same parameters in a single system with N  =  1024 modes. In the first four panels in figure 13 we display the time dependence of the kinetic and potential energies, $e_\mu^{\rm kin}(t)$ and $\newcommand{\e}{{\rm e}} \epsilon_\mu^{\rm pot}(t)$. These functions oscillate as a function of time with amplitude that slowly decreases in the cases $\mu=1, \mu=N-1, \mu =N$, while it slowly increases for $\mu=2$, in this time window. The total energy of each mode shown with a solid red line displays a small downward drift towards, one may expect, a constant value. The Lagrange multiplier shows a residual time variation around a value that is slightly above the $N\to\infty$ prediction but is very close to its finite N modified value $T' + \lambda_{\rm max}^2/(4T')$ (shown with a solid horizontal line). We then determine the mode temperatures from the equipartition of the kinetic and potential energies averaged over a sufficiently long time window. All modes are in equilibrium with themselves in the sense that the potential, kinetic and total mode temperatures coincide except for small deviations present in the largest modes. Panel (f) shows the non-equilibrium temperature profile with the largest modes having higher temperature than the others. T_f, the temperature obtained assuming a paramagnetic final state in equilibrium at a single temperature, see equation (A.21), is clearly different from the mode temperatures and is not the average of them either (not shown), confirming that under these quenches the system does not equilibrate to the paramagnetic state. Finally, panel (g) displays the comparison between the mode inverse temperatures and the frequency dependence of the fluctuation-dissipation inverse temperature. The agreement is very good (except at the edge of the spectrum where both numerator and denominator vanish and it is very hard to control the limiting behaviour even in the equilibrium case, see figure 7). The (yellow) continuous line is the approximate theoretical prediction in equation (186) that captures the numerical behaviour rather accurately. We recall that it was derived assuming that the Lagrange multiplier takes its asymptotic constant value immediately after the quench, at time t  =  0⁺ , and that the ensuing dynamics is the one of independent harmonic oscillators with mode-dependent frequencies.

As already mentioned, parameters in Sector I also permit quenches with energy extraction, realised by x  >  1. We repeated the analysis of the $N\to\infty$ equations for parameters in this regime, choosing, for example, $T'=1.2 \, J_0$ and $J= 1.1 \, J_0$. We do not present the data here as there are no fundamental variations with respect to what we have already discussed for the energy injection case. Having said so, the mode temperatures do present an interesting difference that we discuss with the support of figure 14. Also in this case, the temperatures of the modes are approximately the same for the low lying modes while a mode dependence appears close to the edge of the spectrum. However, the temperatures of the largest modes are, in this case, lower than the temperatures of the lower modes, see panel (a) and its inset. We ascribe this feature to the fact that the quench extracts energy from the system. In panel (b) we confront the mode inverse temperatures to the ones extracted from the fluctuation dissipation ratio in the frequency domain and, once again, the agreement between numerical curves for $N\to\infty$ and finite N data is very good. Moreover, the data are also in good agreement with the outcome of the assumption $z(0^+) =z_f$ that leads to equation (186), shown with a yellow solid line.

For these parameter values the Lagrange multiplier approaches z_f  =  2J. We confirmed this claim with the study of several cases in this Sector (see figure 3). Concerning the behaviour of the other global observables, energies, correlation and linear response, and the mode-dependent ones, we differentiate three cases lying in Sector II: $y > \sqrt{x}$, $y=\sqrt{x}$ and $y<\sqrt{x}$, all with y  >  1 and y  <  x.

$y>\sqrt{x}$. An example of the decay of the two-time global correlation function can be seen in panel (j) in figure 6. It is stationary and it rapidly approaches zero with progressively decaying oscillations around this value. The linear response and correlation function are not related by an FDT with a single temperature (not shown) and in this respect there are no important differences regarding what we have just shown for energy extraction processes in Sector I. For these reasons we chose not to show more data for these parameters.

$y=\sqrt{x}$. A prediction from the asymptotic analysis of the Schwinger–Dyson equations, the fact that on the curve $y=\sqrt{x}$ and $y\geqslant 1$ FDT is satisfied with T_f  =  J, is made explicit in figure 15 where the parametric plot $T_f \chi(t_1-t_2)$ against $C(t_1-t_2)$ is constructed for four pairs of $(x, y=\sqrt{x})$ given in the key. The dotted line is the FDT prediction with T_f  =  J. The agreement between numerics and analytics is very good. Additional support on the fact that the dynamics after the quench occur as in equilibrium at T_f is given in panel (b) in the same figure where quenched and equilibrium correlations are indistinguishable. The latter are obtained by drawing the initial conditions drawn from equilibrium at $T'=T_f=J$ and running the code with the same parameter J. Coincidence of a similar quality (not shown) is found for the other three sets of parameters used in (a).

Zoom In Zoom Out Reset image size

$y<\sqrt{x}$. We chose to show data for $T'=1.25 \, J_0$ and $J=2 \, J_0$, parameters that lie in the region $y<\sqrt{x}$ of the phase diagram where the naive analysis of the $N\to\infty$ equations allows for ageing behaviour. The $N\to\infty$ and finite N data are shown in figure 16. The Lagrange multiplier and kinetic and potential energies approach their expected asymptotic values with an oscillatory behaviour and smoothly decaying amplitude (not shown). The two-time correlation is stationary and shows no ageing. The decorrelation from the initial configuration and from a later configuration at time t₂ behave very similarly. The time-delayed $C(t_1, t_2)$ shows a rapid decay towards a value close to 0.1 around which it oscillates once and next decays towards zero with damped oscillations (a). The linear response function shows a similar effect in the sense of having a fast variation at short time differences and a slower one later (c). The value of R obtained from the numerical integration agrees very well with $R_{{\rm st}}$ from equation (86). The most interesting results concern the comparison between the linear response and the correlation function in the parametric plot in panel (b). The very short time delay behaviour, for C close to 1 and χ close to 0 seems to be described by the slope dictated by T_f, the value of the temperature deduced from an ageing like asymptotic scenario. However, the curves deviate from this straight line for smaller C and larger χ. When the correlation reaches a value close to 0.1 corresponding to its first oscillation, the parametric plot approaches a flat form that ensures the limit $\chi_{\rm st}=1/J$. This second behaviour is reminiscent of what happens in the ageing relaxation of the same model [44].

Zoom In Zoom Out Reset image size

The asymptotic analysis of the $N\to\infty$ equations allow for an ageing solution with a diverging effective temperature for correlation values below the plateau at q, in this region of the parameter space. (In the p  =  3 model, for similar sets of parameters an ageing solution though with a finite T_eff is realised [7].) The numerical solution of the complete set of Schwinger–Dyson equations demonstrates a behaviour with some features that are similar to the approximate asymptotic solution, but no signature of ageing. Indeed, the parametric plot could be interpreted as formed by two pieces, one in which the form is non-trivial close to C  =  1 and another one that is, on average, flat, separated by the correlation value at which the first oscillation occurs. A flat piece in the parametric plot means that the integrated response, and all other functions such as the potential energy, do not get contributions from these time scales and it corresponds to an infinite [44] effective temperature [51, 52]. In practice, the parametric plot is not locally flat for small values of C but, as we can see from the mode analysis in figure 16(e), the temperatures of the corresponding modes do take very large values.

In figure 17 we start discussing the behaviour of the $N\to\infty$ model in Sector III. We use $T'=0.5$ and J  =  0.8, a quench that injects a small amount of energy from the system. Panel (a) proves that the two time correlation approaches a non-vanishing value asymptotically. The horizontal dashed line is at $q=1-T_f/J\simeq 0.44$, the theoretical value derived from the assumption of equilibration à la Gibbs–Boltzmann. Further information about the decay of the correlation functions is given in (c) where we show $C(t, 0)$ and the off-diagonal correlation with the initial configuration C₂(t,0) against time. C₂ starts at $q_{{\rm in}}=0.5$ and decreases monotonically. $C(t, 0)$ quickly decays from 1 with superimposed oscillations. Both curves should reach $q_2=q_0\neq q$ asymptotically and the data are compatible with this claim.

Zoom In Zoom Out Reset image size

Panel (a) in figure 15 shows the value of the asymptotic potential energy as a function of J/J₀ for three different values of $T'/J_0$. The agreement between $e_{\rm pot}^{\,f}$ and $T'/4 + T'J(4J_0)-4J$, the parameter dependence found from the energy conservation and the asymptotic value of z is extremely good.

For quenches with $y = T^{\prime}/J_0 < 1$ and $x>x_{\rm c}(y)=y$, the asymptotic values of $C(t_1, t_2)$ and $C(t_1, 0)=C_2(t_1, 0)$, q and q₀, respectively, do not vanish. However, it is not always easy to extract their functional dependence on x and y, especially for parameters that get away from the shallow quenches $x\simeq 1$. Figure 15(b) shows q and q₀ against J/J₀ for three values of $T'/J_0$, together with the naive single temperature prediction q  =  1  −  T_f/J with $2T_f=2T_{\rm kin}= T'(1+J/J_0)$ drawn with black solid lines. The agreement is good close to x  =  1 though deviations are clear close to the critical line $x_{\rm c}(y)$ and for large energy extraction deep inside this parameter sector. Note that the prediction of equilibration à la Gibbs–Boltzmann is such that $q\neq 0$ at $x_{\rm 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  <  1 (injection) and q₀  <  q for x  >  1 (extraction). We will come back to the behaviour of the correlation function below.

Zoom In Zoom Out Reset image size

Panel (b) in figure 17 shows the parametric $\chi(C)$ curve where one sees that the t  =  0 (red) data are very different from the ones for long time delays t₂. The data for long t₂ suggest that FDT has established at temperature T_f. Another test of FDT, now in the time domain, is shown in (d) with the comparison between the linear response and the time derivative of the correlation. The other panels show the asymptotic values of the Lagrange multiplier (e) and kinetic and potential energies (f). These yield further support to the asymptotic value of the q parameter and T_f estimated analytically under the Gibbs–Boltzmann assumption, since they demonstrate perfect agreement with the asymptotic contributions to the total energy. Nevertheless, we will revisit this claim below when studying deeper quenches in the same sector.

The companion curves for finite N are in figure 19. First of all, panels (a)–(d) display the time dependence of the $\mu=1, \, N/2, \ N-1, \ N$ mode energies in a system with N  =  1024. While the modes $\mu=1, \ N/2$ show the usual oscillatory behaviour of a harmonic oscillator, the largest modes $\mu=N-1$ and $\mu=N$ are clearly out of equilibrium. The Lagrange multiplier is approximately constant and equal to the largest eigenvalue, within numerical accuracy. The spectrum of mode temperatures is plotted in (f) with a zoom over the largest modes in its inset. The profile is an equilibrium one, with $T_\mu$ being independent of μ, apart from the deviations close to the edge of the spectrum. Finally, (g) shows a comparison between the inverse mode temperatures computed numerically for the finite N system and the inverse temperature, $\frac{1}{T_f}$, extracted from the assumption of independent harmonic oscillators (given by equation 186). Higher modes (low frequency) have temperatures slightly above the temperature T_f while lower modes (high frequencies) have temperatures slightly below T_f, that is shown with an horizontal yellow line. This fact is another warning concerning the claim of complete equilibration to a Gibbs–Boltzmann measure.

Zoom In Zoom Out Reset image size

Indeed, although the data for the shallow quench that we have just described might have suggested equilibration to a Gibbs–Boltzmann probability distribution, the detailed comparison of the full time-delay dependence of the correlation function after a quench and in equilibrium (no quench) at parameters J and T_f that are the ones that the equilibrium measure would have, prove that such a steady state is not reached by the dynamics. This statement is proven in figure 20 where we display the self-correlations stemming from the two procedures for three choices of quenches: to the critical line x  =  y (a case that will be further studied in the next subsubsection), the shallow quench, and a deep quench.

Zoom In Zoom Out Reset image size

We note that the comparison of the linear response functions after a quench and in equilibrium yield identical results. It is the correlation function the one that deviates from Gibbs–Boltzmann equilibrium.

Figure 20(a) has already shown non-equilibrium dynamics on the critical line y  =  x for y  <  1. We give here further evidence for this fact. In figure 21 we show results for $T'=0.5 < T^0_{\rm c}$ and J  =  0.5, that is to say, point c on the phase diagram in figure 5. This quench injects a relatively small amount of energy into the system, $\Delta e\simeq 0.375$ and takes the parameters to be on the critical line y  =  x.

Zoom In Zoom Out Reset image size

The self correlation is shown in panel (a) of figure 21. Stationarity is clear for short time delays $t_1-t_2$, for t₂  >  0 and there is some remanent waiting time dependence at these short t₂. This double behaviour is reminiscent of what is seen in the relaxational dynamics where a sharp separation of time-scales exists. The data suggest that $q_0=\lim_{t\rightarrow \infty} C(t, 0)$ is different from $q=\lim_{t_1-t_2\rightarrow \infty} \lim_{t_2\to\infty}C(t_1, t_2)$, although it seems hard to determine these values from the numerics with good precision. In the event of a two step relaxation of $C(t_1, t_2)$ with an approach to a plateau at q and a further decay from it to zero, the plateau should be at q  =  1  −  T_f/J  =  0.25, shown with a dashed horizontal line. The data still lie below this value and we infer that they might not converge to a plateau but simply decay to zero.

Further information about the decay of the correlation functions is given in (c) where we show $C(t, 0)$ and the off-diagonal correlation with the initial configuration C₂(t,0) against time. C₂ starts at $q_{{\rm in}}=0.5$ and decreases monotonically. $C(t, 0)$ quickly decays from 1 with superimposed oscillations. Both curves seem to join and slowly and monotonically decay to zero.

The parametric plot of the linear response function, $\chi(t_1, t_{2})$, against the correlation, $C(t_{1}, t_2)$, for fixed t₂ and using $t_1-t_2$ as a parameter, for three different values of the waiting time t₂, is shown in panel (b). The $\chi(C)$ curve for t  =  0 does not have any special form. The curves for late t₂ are close to the straight line 1/T_f for time-delays that correspond to the first oscillations of the correlation and linear response, they then oscillate, and for longer time-delays the parametric construction becomes flat, with χ approaching, for $C\to 0$, the expected static susceptibility 1/J  =  2. Again, this second flat regime is reminiscent of the one found for the relaxation stochastic dynamics at and below the critical temperature [86, 87].

The fluctuation-dissipation theorem relating χ to C in a stationary regime with target temperature T_f  =  0.375 obtained from equation (A.19), is indicated as a straight dashed line. The presentation in panel (d) confirms the agreement between linear response and time variation of the time-delayed correlation function dictated by the FDT at temperature T_f for $C\geqslant 0.2$, say. However, it clearly breaks for smaller values of C. As for the other quenches, we checked that the response $R(t_2+\tau, t_2)$ coincides with the one derived by anti-transforming the theoretical Fourier transform of the stationary asymptotic response, given by equation (86), $R_{{\rm st}}$.

In (e) we provide a close look at the time evolution of $z(t)$. As one can see, $z(t)$ approaches the predicted value z_f  =  2J. From the time evolution of the energy (f), we observe that the total energy is constant in time, as it should be, while the kinetic and potential energies quickly approach their asymptotic values, which coincide within numerical accuracy with the ones predicted in section 4.4.3. These values could be mistaken for ${e}^{\,f}_{\rm kin}=T_f/2$ and ${e}^{\,f}_{\rm pot}=-J^2(1-q^2)/(2T_f)$, with q  =  1  −  T_f/J  =  0.25, the Gibbs–Boltzmann equilibrium predictions, though, the system is not in equilibrium.

In figure 22 we show results for $T'=0.5 < T^0_{\rm c}$ and J  =  0.25. This quench injects a large amount of energy into the system, $\Delta e=0.5625$, which is sufficient to take it out of the initial condensed state. Had the system reached an equilibrium paramagnetic state asymptotically after the quench its temperature would be $T_f \simeq 0.320$, so that $T_f/J \simeq 1.281$, from Equation (A.20). This, however, is not consistent within the $N\to\infty$ analysis and, accordingly, it is not realised dynamically.

Zoom In Zoom Out Reset image size

The self correlations shown in (a) present large oscillations with weakly decaying amplitude around the expected asymptotic value $\lim_{t_1-t_2\rightarrow \infty} \lim_{t_2\to\infty} C(t_1, t_2) = 0$ for all values of t₂. $C(t_1, t_2)$ satisfies stationarity for sufficiently late t₂, as one can see from the plot where the curves for different values of t₂ coincide, within numerical accuracy. The parametric plot of the susceptibility χ against the correlation C, shown in (b), is rather complex and very far from linear. The FDT relating χ to C in a stationary regime with putative target temperature $T_f \simeq 0.320$ is indicated as a straight dashed line. This behaviour is confirmed by the time evolution of the stationary response function, $R(t_1, t_2)$, see panel (d). $R(t_1, t_2)$ does not coincide with $-\frac{1}{T_f} \partial_{t_1}C(t_1, t_2)$, with T_f the target temperature. In panel (c) we also show the time evolution of the off-diagonal correlation with the initial configuration, C₂(t,0). In contrast to the previous cases, C₂(t,0) is oscillating in time around the value q₂  =  0. Panel (e) demonstrates that $z(t)$ is far from $z_f=T_f + J^2/T_f$ relative to an equilibrium paramagnetic state but oscillates around $T'+J^2/T'$, consistently with our claims so far. From panel (f) we observe that the kinetic and the potential energies relax to the predicted asymptotic values specified in section 4.4.3. The Fourier transforms of the correlation and response functions for two different values of t₂, as indicated in the key, are shown in panel (g). The black solid lines represent the theoretical prediction for the real and imaginary part of the Fourier transform of the response function in the stationary regime, $\hat{R}_{{\rm st}}(\omega)$, given by equation (86), with parameters m  =  1, J  =  0.25 and z_f  =  0.625. They are in excellent agreement with the numerical results. In panel (h) we show the ratio $-{\rm Im}{\hat{R}(\omega)}/(\omega\hat{C}(\omega))$. The dashed line is  −1/T_f and is clearly off the data. The different curves were computed for various waiting times t₂ given in the key. In figure 23 we compare the FDR (averaged over different waiting times to get rid of the undesired oscillations) to the mode temperatures of the finite N system. The correspondence between the two ways of extracting the temperatures is good. The yellow line is the approximate prediction in equation (186) stemming from the independent harmonic oscillator approaximation that for this kind of quench is quite far from the numerical results though it captures the qualitative features.

Zoom In Zoom Out Reset image size

The analysis of the constants of motion should shed light on the 'distance' from complete equilibration to a Gibbs–Boltzmann probability density, especially in Sector III where shallow quenches in the $N\to\infty$ Schwinger–Dyson formalism suggested proximity from this description. We here compare the actual values of the $\langle I_\mu\rangle$ to the 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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:imuf} \overline{\langle I_{\mu\neq N} \rangle}_{T_f} = \frac{T_f}{z_f-\lambda_\mu} + \frac{T^2_fJ_0}{J} \ \frac{1}{z_f-\lambda_\mu} \ S_\mu^{(1)} + T_f \ \frac{q_f}{\lambda_N-\lambda_\mu} + \frac{{T_f}^2 J^2_0}{J^2} \; S^{(2)}_\mu , \nonumber \end{align} \tag{ 214 }$$

where we allowed for the condensation of the largest mode. We wish to compare this expression to the one in equation (210). After some lengthy calculations, using the $N\to\infty$ values

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle q_{\rm in} = 1 - \frac{T'}{J_0} , \qquad z_f = 2J = \frac{J}{J_0} z(0^-) , \qquad T_f = \frac{T'}{2} \left(1 + \frac{J}{J_0}\right) , \qquad q_f =1-\frac{T_f}{J} , \nonumber \end{align} \tag{ 215 }$$

we find that the difference $\Delta I_{\mu\neq N} = \langle I_{\mu\neq N}(0^+) \rangle - \overline{\langle I_{\mu\neq N}\rangle}_{T_f}$ is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta I_{\mu\neq N} = - \frac{{T'}^2 }{2} \left(\frac{J_0}{J}-1\right)^2 \ \frac{1}{\lambda_N^{(0)}-\lambda^{(0)}_\mu} \ S_\mu^{(1)} \nonumber \end{align} \tag{ 216 }$$

a finite, non-zero, value for all $\mu\neq N$.

The Nth mode.

The Nth mode has a different scaling with system size. We have already computed $\langle I_N(0^+)\rangle$ in equation (212) and we want to compare it to what it should read in equilibrium at T_f:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \overline{\langle I_N\rangle}_f = q_f \left(1+ \frac{T_f J_0}{J} \ S_N^{(1)} \right) N +\frac{T^2_f J^2_0}{J^2} \ S_N^{(2)} . \nonumber \end{align} \tag{ 217 }$$

The difference between the two, $\newcommand{\e}{{\rm e}} \Delta I_{N} \equiv \langle I_{N}(0^+)\rangle - \overline{\langle I_{N}\rangle}_{T_f}$, is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta I_{N} = \left[ q_{\rm in} - q_f + (q_{\rm in} T' - q_f T_f) \frac{J_0}{J} \ S_N^{(1)} \right] N + \left({T'}^2 - T_f^2 \frac{J_0}{J} \right) \frac{J_0}{J} \ S_N^{(2)} . \nonumber \end{align} \tag{ 218 }$$

The first term diverges linearly with N. We have already argued that the second one is $O(1)$, see the discussion after equation (212). 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'$ in the large N limit, and $S_N^{(1)} \to -1/J_0$,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:slope_in} \frac{\Delta I_{N}}{N} \to - \frac{{T'}^2}{4J_0^2} \left(\frac{J_0}{J} - 1\right)^2 . \nonumber \end{align} \tag{ 219 }$$

This form is validated by the numerical data in figure 25 for three pairs of $T'$ and J, all in Sector III.

Zoom In Zoom Out Reset image size

Quite clearly, the differences between the actual $\langle I_\mu\rangle$ and the ones in equilibrium at a temperature T_f are proportional to J₀/J  −  1, a factor that vanishes for J₀  =  J. Otherwise, the differences are finite for $\mu\neq N$ and proportional to N for the last mode, making the mode-by-mode difference non-vanishing for $J\neq J_0$. This fact confirms, then, the lack of equilibration to a Gibbs–Boltzmann measure with a single T_f.

The special line ${y=\sqrt{x}}$.

On the curve $y=\sqrt{x}$ the asymptotic analysis for $N\rightarrow\infty$ predicts thermalization at temperature T_f  =  J. The numerical analysis of the Schwinger–Dyson equations confirms this prediction as the correlation and linear response are linked by FDT and the time-dependence of the correlation function after an instantaneous quench is identical (within numerical accuracy) to the one found in equilibrium at this temperature. We shall briefly analyse the behaviour of the constants of motion in this case.

On the special line $y=\sqrt{x}$, ${T'}^2 = J J_0$ and equation (213) yields the following expression for the averaged integrals of motion

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle I_{\mu}(0^+)\rangle=\frac{1}{z(0^-)/J_0-\overline{\lambda}_{\mu}}\left(\sqrt{\frac{J}{J_0}}+1+2S^{(1n)}_{\mu} \right) , \nonumber \end{align} \tag{ 220 }$$

where $\overline{\lambda}_{\mu}=\lambda^{(0)}_{\mu}/J_0$ are the normalised eigenvalues that vary in the interval $(-2, 2)$. From equation (214) with q_f  =  0, z_f  =  2J valid in the $N\rightarrow\infty$ limit, and using equation (208) the constants of motion in an equilibrium state at temperature T_f are

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \overline{\langle I_{\mu}\rangle}_{T_f}=\frac{1}{2-\overline{\lambda}_{\mu}}\left(2+2S^{(1n)}_{\mu} \right) . \nonumber \end{align} \tag{ 221 }$$

We observe that the two sets of values are very similar if J is close to J₀. Numerically, we find that, in the special case $T^{\prime}=1.25 \ J_0$ and $J=1.5625 \ J_0$, the $\Delta I_{\mu}$ is of order of 10⁻². We conclude that even in this case, in which the asymptotic analysis predicts that the global, mode-averaged quantities, behave as in equilibrium, a prediction that seems to be confirmed by the numerics, the constants of motion are not exactly the same in the initial state and in the thermal state the system would reach in case of thermalisation.

In order to properly interpret these results, it is important to keep in mind that, strictly speaking, the conserved energy dynamics of an isolated (finite size) system should keep memory of the initial conditions, even if the system is non-integrable. In our problem, we see this information encoded in the $I_\mu$s. More so, not even in the $N\to\infty$ limit this memory is erased as the $\Delta I_{\mu}$s remain finite.

We have seen in section 5.7 that the $N\to\infty$ system decouples into independent harmonic oscillators in the asymptotic long time limit (taken after $N\to\infty$) since $z(t) \to z_f$. A natural idea is to check whether we can identify the integrals of motion $\langle I_\mu(0^+)\rangle$ with the ones that an ensemble of harmonic oscillators with spring constants $m \omega_\mu^2 = z_f-\lambda_\mu$ in equilibrium at the temperatures $T_\mu$ would have.

On the one hand, we note that $\langle s_\mu^2(t)\rangle$ and $\langle\,p_\mu^2(t)\rangle$ are not constant for harmonic oscillators but their time averages are, so we evaluate $\overline{\langle I_\mu \rangle}$ finding, in this framework,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \overline{\langle I_\mu \rangle}_{\rm osc} = \frac{T_\mu}{z_f-\lambda_\mu} + \frac{T_\mu}{N} \sum_{\nu(\neq \mu)} \left[ \frac{T_\nu}{z_f-\lambda_\mu} + \frac{T_\nu}{z_f-\lambda_\nu} \right] \frac{1}{\lambda_\nu-\lambda_\mu} . \nonumber \end{align} \tag{ 222 }$$

In the $N\to\infty$ limit the value of z_f is expected to be $T'+J^2/T'$ for x  <  y and 2J for y  <  x. Imposing the identity between the $\overline{\langle I_\mu \rangle}_{\rm osc}$ and the $\langle I_\mu(0^+)\rangle$ given in equations (209) and (211) should yield a better estimate of the temperatures $T_\mu$ than the one explained in section 5.7. We leave this analysis aside for the moment.

On the other hand, once the oscillators have decoupled and the Lagrange multiplier stabilised, the mode total energies, or the time-averaged kinetic and potential energies are also constants of motion. In the next subsection we investigate to what extent the $\langle I_\mu\rangle$ are proportional to $e_\mu^{\rm tot}= 2\overline{e_\mu^{\rm kin}}$.

Figure 26 shows the spectrum of integrals of motion $J\langle I_{\mu}\rangle$ (red data points) together with the mode kinetic temperatures $T^{{\rm kin}}_{\mu}=2\overline{e_\mu^{\rm kin}}$ (that, we have checked, are in agreement within numerical accuracy with the potential ones $T^{{\rm pot}}_{\mu}$ 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\langle I_\mu \rangle$ are noisier than the ones for $T_\mu^{\rm 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 equation (124) and we see that the mode temperatures are very close to it again away from the edge of the spectrum and satisfying the constraint $N^{-1} \sum_\mu T_\mu^{\rm kin}=T_{\rm kin}$.

Zoom In Zoom Out Reset image size

The parameters in Sector II, panel (b), are on the curve $y=\sqrt{x}$ on which the data for the global correlation and linear response show equilibrium at a single temperature T_f  =  J. The blue data points for the mode kinetic energies are precisely on this value. The integrals of motion scatter, however, quite a lot around this value. Importantly enough, the difference $\Delta I_{\mu}$ is very small in this case.

In sector III, panel (c), the data for $J \langle I_\mu\rangle$ tend to be relatively flat in the bulk of the spectrum and very close to $T_{\rm kin}=T_f=T'(1+J/J_0)/2$. This result can be derived from equation (209) taking advantage of a simple rearrangement of the two sums, $S_\mu^{(2)} = (S_\mu^{(1)}-S_N^{(1)})/(z(0^-)-\lambda_\mu^{(0)})$ valid for $\mu\neq N$. In particular, taking $\lambda_\mu^{(0)}$ right at the middle of the spectrum, $S_\mu^{(1)}=0$ by symmetry and $J\langle I_0\rangle = T_f$ for $N\to\infty$. The incipient divergence close to the right edge of the spectrum, with the deviation from this constant value, is also clear. By using a maximal value of 8 in the vertical axes we have explicitly left aside the points for $\mu\simeq N$ that take much larger values. For instance, $J \langle I_N \rangle\simeq 180$. The approximate mode temperatures in equation (186) are all identical to $T_f=T_{\rm kin}$ in this Sector, consistently with the numerical data away from the edge. The fact that the system is not in proper Gibbs–Boltzmann equilibrium is due to the fact that the higher lying integrals of motion do not comply with this temperature.

The data for $J \langle I_\mu\rangle$ in Sector III look very similar to the ones found in Sector IV, see panel (d). Indeed, the data almost coincide far from the edge, since they are both determined by equation (209) that has a weak dependence on J (the only parameter that takes a different value in panels (c) and (d)). Close to the edge, the $J \langle I_\mu\rangle$ differ since in III (c) there is scaling with N while in IV (d) there is not. The mode kinetic energies are mode independent and identical to $T'(1+J/J_0)/2$ contrary to what the approximation in the last line of equation (186) tells. This means that for these quenches the assumption in equation (183) fails.

We reckon that the kinetic temperature $T_{\rm kin} =2 \overline{e}_{\rm kin}$ should be equal to the sum of the mode kinetic temperatures $T_{\rm kin} = 1/N \, \sum_{\mu=1}^N T_\mu^{\rm kin}$. We have checked numerically this property.

We leave a more detailed analysis of the comparison between the two and their use to build a GGE for a future publication.

Let us assume that the limiting value of the Lagrange multiplier is a constant

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \lim_{t_1\to\infty} z(t_1) = z_f . \nonumber \end{align} \tag{ A.1 }$$

Recalling the definitions given in the main part of the paper, the limits of the correlation functions are

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle q = \lim_{t_1-t_2\to\infty} \lim_{t_2\to \infty} C(t_1,t_2) , \qquad q_0 = \lim_{t_1\to \infty} C(t_1,0) , \qquad q_2 = \lim_{t_1\to\infty}C_2(t_1,0) . \nonumber \end{align} \tag{ A.2 }$$

The linear response was analysed in the main body of the paper and, independently of the quench parameters it is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat R(\omega) = \frac{1}{2 J^2} \left[ (-m\omega^2 + z_f) \pm \sqrt{(-m\omega^2 + z_f)^2 - 4 J^2} \right] \label{eq:response_fourier_transform} \nonumber \end{align} \tag{ A.3 }$$

in the frequency domain, where the Fourier transform has been computed with respect to the time difference $t_1-t_2$. This result only assumes a long-time limit in which z_f is time-independent.

We could estimate the asymptotic value of z(t₁) taking the long t₁ limit of equation (81); however, without the use of FDT we cannot compute the integral involved.

We are tempted to propose

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle q_0 = q_2. \nonumber \end{align} \tag{ A.4 }$$

The interpretation of this result in terms of the evolution of real replicas is that the asymptotic value of the self-correlation between times t₁ and 0, q₀, is the same as the asymptotic value of the correlation between two replicas evaluated at the same times t₁ and 0 with t₁ diverging.

Allowing for the two-time correlation function not to decay to zero but to a finite value q, we separate this contribution explicitly and we write $C(t_1, t_2) = q + C_{\rm st}(t_1-t_2)$ with $\lim_{t_1-t_2\to\infty} C_{\rm st}(t_1-t_2)=0$. We also use $\lim_{t_1\to\infty} C_2(t_1, 0) =q_2$ leaving the possibility of there being a different asymptotic value for this quantity open.

The dynamic equation now becomes an equation on $\tilde C$ that reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left(m \partial^2_{t_1-t_2 } + z_f \right) (q + C_{\rm st}(t_1-t_2)) =&~ \frac{JJ_0}{T'} (q_0^2 - q_2^2) + 2 J^2 q \; \lim_{t_1\to\infty} \int_0^{t_1} {\rm d}t' \, R(t_1-t') \nonumber \\ & + J^2 \int_0^{t_1} {\rm d}s \; R(t_1-t_2+s) C_{\rm st}(-s) + J^2 \int_0^{t_2} {\rm d}s \; C_{\rm st}(t_1-t_2+s) R(s) . \label{Cst} \nonumber \end{align} \tag{ A.5 }$$

Using the causal properties of the linear response, one can extend the lower limit of the last integral to $-\infty$ and safely take the upper limit to infinity since we are interested in the long time limit $t_2\to\infty$ while keeping $t_1-t_2$ fixed. The Fourier transform of this term with respect to $t_1-t_2$ yields $R(-\omega) C_{\rm st}(\omega)$ that using $R(-\omega) = R^*(\omega)$ simply becomes $R^*(\omega) C_{\rm st}(\omega)$. Proceeding in a similar way with the first integral one finds that it equals $\hat R(\omega) C_{\rm st}(\omega)$. Thus,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left(- m \omega^2 +z_f - J^2 \hat R(\omega) \right) C_{\rm st}(\omega) =& \left(-z_f q + \frac{JJ_0}{T'} (q_0^2 - q_2^2) + 2 J^2 \, q \lim_{t_1\to\infty} \int_0^{t_1} {\rm d}t' \, R(t_1-t') \right) \delta(\omega) \nonumber \\ & + J^2 \, C_{\rm st}(\omega) R^*(\omega) . \nonumber \end{align} \tag{ A.6 }$$

The factor between parenthesis in the first term on the right-hand-side can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {-z_f q + \frac{JJ_0}{T'} (q_0^2 - q_2^2) + 2 J^2 \, q \lim_{t_1\to\infty} \int_0^{t_1} {\rm d}t' \, R(t_1-t') =0 . } \label{cond2} \nonumber \end{align} \tag{ A.7 }$$

If we now assume $q_0=q_2$, this equation imposes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm{if}} \;\; q\neq 0 \qquad {\rm{then}} \;\; z_f = 2J \;\; {\rm{and}} \;\; \lim_{t_1\to\infty} \int_0^{t_1} {\rm d}t' \, R(t_1-t') = 1/J, \nonumber \end{align} \tag{ A.8 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm{otherwise}} \;\; q=0. \nonumber \end{align} \tag{ A.9 }$$

The remaining equation, using equation (85) to replace the parenthesis on the left-hand-side by $1/\hat R(\omega)$ is recast as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J^2 C_{\rm st}(\omega) |\hat R(\omega)|^2 = C_{\rm st}(\omega) . \nonumber \end{align} \tag{ A.10 }$$

At each frequency this equation has two possible solutions

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C_{\rm st}(\omega) \neq 0 \;\;{\rm{and}}\;\; |\hat R(\omega)|^2=1/J^2 \qquad {\rm{or}}\qquad C_{\rm st}(\omega) = 0 . \nonumber \end{align} \tag{ A.11 }$$

In the frequency domain in which the linear response (A.3) is complex, that is Im$\hat R(\omega)\neq 0$, one can easily check that it satisfies $|\hat R(\omega)|^2=1/J^2$. 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 1/J². For instance, at zero frequency in cases in which x  >  y, $\hat R(\omega=0)=1/J$, verifying that its square is 1/J². However, for x  <  y, $\hat R(\omega=0)=1/T'$ and its square is different from 1/J². In these cases, the Fourier transform of the decaying part of the correlation, $\hat C(\omega)$ vanishes. We have tested this statement numerically and several examples can be seen in the main part of the paper.

The asymptotic limit of the equation for C(t₁,0) implies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z_f q_0 = J^2 q_0 \lim_{t_1\to\infty} \int_0^{t_1} {\rm d}t' \, R(t_1-t') + \frac{JJ_0}{T'} (q_0 -q_{{\rm in}} q_2) \label{eq:q0} \nonumber \end{align} \tag{ A.12 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle q_{{\rm in}} &= \left\{\begin{array}{@{}ll@{}} 1 - \frac{T'}{J_0} & \qquad T<T^0_s \nonumber \\ 0 & \qquad T>T^0_s \end{array} \right. \nonumber \end{align} \tag{ A.13 }$$

and $T^0_s$ the critical temperature of the initial potential energy equation (A.12) is the equivalent of equation (82) in [7].

This equation admits the solution $q_0=q_2=0$ or

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z_f = J^2 \lim_{t_1\to\infty} \int_0^{t_1} {\rm d}t' \, R(t_1-t') + \frac{JJ_0}{T'} (1-q_{\rm in}) \label{eq:zfnueva0} \nonumber \end{align} \tag{ A.14 }$$

if we assume $q_0=q_2$. Recalling that $q_{\rm in} =1 -T'/J_0$, the remaining equation becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z_f = J^2 \lim_{t_1\to\infty} \int_0^{t_1} {\rm d}t' \, R(t_1-t') + J \label{eq:zfnueva} \nonumber \end{align} \tag{ A.15 }$$

and is consistent for z_f  =  2J and $\lim_{t_1\to\infty} \int_0^{t_1} {\rm d}t' \, R(t_1-t') =1/J$.

For condensed (first line in the left) and paramagnetic (second line in the left) initial conditions going to condensed (first line in the right) and paramagnetic-like (second line in the right) states the energy conservation law reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{m}{m_0} \frac{T'}{2} - \frac{J J_0}{2T'} \left[ \begin{array}{@{}c@{}} 1 - \left(1-\frac{T'}{J_0}\right)^2 \nonumber \\ 1 \end{array} \right] = \frac{T_f}{2} - \frac{J^2}{2T_f} \left[ \begin{array}{@{}c@{}} 1 - \left(1-\frac{T_f}{J}\right)^2 \nonumber \\ 1 \end{array} \right] - \frac{J J_0}{2T'} (q_0^2 - q_2^2) . \nonumber \end{align} \tag{ A.18 }$$

We have already determined $q_0=q_2$ so the last term vanishes identically. From this condition we find the following expressions for T_f depending on the kind of quench performed:

• From condensed to condensed
  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_f = \frac{J}{2} \left(\frac{m J_0}{m_0 J} + 1 \right) \frac{T'}{J_0} \qquad \qquad {\rm{(III)}} . \label{eq:Tf_from_cond_to_cond} \nonumber \end{align} \tag{ A.19 }$$
• From condensed to paramagnetic
  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_f = \frac{J}{2} \frac{J_0}{T'} \left\{\frac{mJ_0}{m_0 J} \left(\frac{T'}{J_0} \right)^2 - (1-q_{\rm in}^2) + \sqrt{\left[ \frac{mJ_0}{m_0 J} \left(\frac{T'}{J_0} \right)^2 - (1-q_{\rm in}^2) \right]^2 + \left(\frac{2T'}{J_0}\right)^2 } \right\} \;\; \qquad {\rm{(IV)}} . \label{eq:Tf_from_cond_to_para} \nonumber \end{align} \tag{ A.20 }$$
• From paramagnetic to paramagnetic
  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_f = \frac{J}{2} \frac{J_0}{T'} \left\{\frac{mJ_0}{m_0 J} \left(\frac{T'}{J_0} \right)^2 - 1 + \sqrt{\left[ \frac{mJ_0}{m_0 J} \left(\frac{T'}{J_0} \right)^2 - 1 \right]^2 + \left(\frac{2T'}{J_0}\right)^2 } \right\} \qquad {\rm{(I)~and~(IIa)}} . \label{eq:Tf_from_para_to_para} \nonumber \end{align} \tag{ A.21 }$$
• From paramagnetic to ageing
  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_f = J \left[ 1 + \frac{m J_0}{m_0 J} \frac{T'}{2J_0} - \frac{J_0}{2T'} \right] \qquad {\rm{(IIb)}}. \label{eq:Tf_from_para_to_ageing} \nonumber \end{align} \tag{ A.22 }$$

We have chosen to single out in these expression the parameters $T'/J_0$ and $(m J_0)/(m_0 J)$ that we have used in [7] to characterise the dynamic phase diagram of the $p\geqslant 3$ model. In the two cases (condensed to condensed and paramagnetic to paramagnetic) in which the no quench limit can be taken we naturally recover $T_f=T'$. In the ageing case T_f should be the temperature of the stationary regime. As we will show from the numerical solution to the full dynamic equations, and the mode by mode analysis, not all these cases are realised.

Below we prove analytically that a state with a single T_f characterising the fluctuation-dissipation relation can be realised for particular choices of the parameters only. Moreover, these conditions are not restrictive enough, as the numerical solution of the full equations show that FDT is not realised for quenches allowed by them.

Let us consider the cases $y\geqslant 1$ and $y\leqslant 1$ separately.

The case ${y \geqslant 1}$: Quench from a paramagnetic equilibrium state.

The double condition $0\leqslant T_f \leqslant J$ translates into

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle -1 \leqslant \frac{1}{2} \frac{mJ_0}{m_0J} \frac{T'}{J_0} - \frac{1}{2} \frac{J_0}{T'} \leqslant 0 \nonumber \end{align} \tag{ A.23 }$$

that in terms of x and y reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle -1 \leqslant \frac{y}{2x}- \frac{1}{2y} \leqslant 0 \nonumber \end{align} \tag{ A.24 }$$

and yields

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 0 \leqslant y^2-x + 2 xy \qquad {\rm{and}}\qquad y^2\leqslant x \label{eq:upper-bound0} . \nonumber \end{align} \tag{ A.25 }$$

The first condition in equation (A.25) is satisfied for

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 0 \leqslant (y-y_+) (y-y_-) \nonumber \end{align} \tag{ A.26 }$$

with y₊ and y₋ the two roots of the quadratic equation, that is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y_{\pm} = -x \pm \sqrt{x^2 + x} , \nonumber \end{align} \tag{ A.27 }$$

a positive and a negative value. One should then have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y \geqslant y_+ \qquad \qquad {\rm{or}} \qquad \qquad y \leqslant y_- . \nonumber \end{align} \tag{ A.28 }$$

As $y\geqslant 1$, $y\leqslant y_-$ is not possible. Instead, y  >  y₊ is trivially satisfied. The second condition in equation (A.25) is the only restrictive one and reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y \leqslant \sqrt{x}\equiv f(x) . \nonumber \end{align} \tag{ A.29 }$$

This upper bound is the dashed line representing $f(x)$ in the phase diagram for y  >  1. Note that this curve has no finite limit.

The case ${y \leqslant 1}$: Quench from a condensed equilibrium state.

Using $q_{\rm in} =1 -T'/J_0$, one has

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{T_f}{J} = \frac{1}{2} \frac{J_0}{J} \frac{T'}{J_0} - \frac{1}{2} \frac{T'}{J_0} \nonumber \end{align} \tag{ A.30 }$$

and the condition to use reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 0 \leqslant \frac{y}{2} \left(\frac{1}{x} + 1 \right) \leqslant 1. \nonumber \end{align} \tag{ A.31 }$$

The lower bound is trivially satisfied while the upper bound implies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y \leqslant \frac{2x}{1+x} \equiv g(x) \qquad \qquad {\rm{for}} \;\; y\leqslant 1 . \nonumber \end{align} \tag{ A.32 }$$

The dashed line for y  <  1 in the phase diagram represents $g(x)$.

In the condensed case the energy conservation implies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_f = J \; \frac{y(1+x)}{2x} , \qquad q=1-\frac{y}{2} - \frac{y}{2x} \qquad {\rm{and}}\qquad \overline e_{\rm pot} = -J \left(1-\frac{x}{4y}- \frac{y}{4} \right) . \label{eq:Tf-condensed} \nonumber \end{align} \tag{ A.33 }$$

We note that q vanishes on the curve $y=2x/(1+x)$, that it equals one for y  =  0 and that, for fixed x  >  1 it increases from $(x-1)/(2x)$ at y  =  1 to 1 at y  =  0. Moreover, q is different from zero on the lines y  =  x and y  =  1.

We now reason as follows. We know, from the numerical analysis of the Schwinger–Dyson set of equations, that $\chi_{\rm st} = \lim_{t\to\infty} \int_0^t {\rm d}t' R(t-t') = 1/T'$ for x  <  y and $\chi_{\rm st} = \int_0^t R(t) = 1/J$ for x  >  y. If we evaluate the integral of the linear response using FDT at a single temperature we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi_{\rm st} = \lim_{t\to\infty} \int_0^t {\rm d}t' R(t-t') = \lim_{t\to\infty} \int_0^t {\rm d}t' \frac{1}{T_f} \frac{\partial}{\partial t'} C_{\rm st}(t-t') = \left\{\begin{array}{@{}ll@{}} {\frac{1}{T_f}} & \qquad y>1 \nonumber \\ {\frac{1}{T_f} (1-q) }= \frac{1}{J} & \qquad y<1. \end{array} \right. \nonumber \end{align} \tag{ 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 from the conservation of the energy is consistent with these conditions. We list below the conclusions drawn in the four parameter Sectors of the phase diagram.

y  >  1 and x  <  y.

y  >  1 and x  >  y.

y  <  1 and x  >  y.

y  <  1 and x  <  y.