# Accelerating, to some extent, the $p$-spin dynamics

Federico Ghimenti<EMAIL_ADDRESS>Frédéric van Wijland
Laboratoire Matière et Systèmes Complexes, Université de Paris and CNRS (UMR
7057), 10 rue Alice Domon et Léonie Duquet, 75013 Paris, France

###### Abstract

We consider a detailed-balance violating dynamics whose stationary state is a
prescribed Boltzmann distribution. Such dynamics have been shown to be faster
than any equilibrium counterpart. We quantify the gain in convergence speed
for a system whose energy landscape displays one, and then an infinite number
of, energy barriers. In the latter case, we work with the mean-field
disordered $p$-spin, and show that the convergence to equilibrium or to the
nonergodic phase is accelerated, both during the $\beta$ and
$\alpha$-relaxation stages. An interpretation in terms of trajectories in
phase space and of an accidental fluctuation-dissipation theorem is provided.

††preprint: APS/123-QED

##  I Markovian acceleration of the dynamics

Designing algorithms able to locate the minima of a cost function in a large
dimensional space is at the core of optimization problems lying well beyond
the realm of physics. Yet the methods used, e.g. in computer science or
applied mathematics, often borrow from physical approaches. Indeed the
Boltzmann distribution $P_{\text{\tiny B}}\sim\text{e}^{-\beta H}$ for a
system with energy $H$ concentrates, at low temperature $T=1/\beta$, in
regions where $H$ is the smallest, and one of the important questions is thus
whether one can efficiently sample the Boltzmann distribution in the low
temperature regime, thereby locating the minima of $H$ [1]. This is achieved
by endowing the degrees of freedom with dynamical evolution rules, often
inspired by Langevin dynamics.

The question of how to fine-tune these dynamical evolution rules to accelerate
the convergence to the steady state has grown into a field of its own, that is
discussed at length in the physics literature as it touches to such diverse
problems as the glass transition or finding the stable configuration of large
proteins [2]. It is of course abundantly discussed in the fields of
optimization [3] and machine learning [4] as well. For example, to accelerate
the Stochastic Gradient Method, Nesterov [5] introduced a modified version of
a momentum based (a.k.a. underdamped in the physics language) method. In the
physics literature too, many cleverly accelerated dynamics have been
introduced, such as the replica exchange [6], parallel tempering [7] or the
SWAP Monte Carlo method [8, 9, 10].

When designing alternative dynamics, it is also possible not to respect the
detailed balance condition, now a widely used procedure in Monte Carlo schemes
[11, 12], that leads to a faster convergence towards equilibrium [13]. A
particular subclass of algorithms is obtained by subjecting the system of
interest to a divergenceless force such that the stationary state remains the
Boltzmann distribution. This has been proposed by several authors over the
last three decades [14, 15, 16, 17, 18, 19, 20, 21, 22]. The terminology for a
concrete implementation refers to a ”lifting” procedure[23, 24, 25] , which
consists in introducing auxiliary degrees of freedom which drive the system
away from equilibrium (but the statistics of which preserves the Boltzmann
distribution for the degrees of freedom of interest). In the physics
literature this was for instance demonstrated on the mean-field Ising model by
[26, 27].

In this work we shall build on the version proposed by Ichiki and Ohzeki [18]
(which we hereafter call the IO dynamics). They start from two independent and
identical systems with a reversible dynamics whose equilibrium distribution
$P_{\text{\tiny B}}$ is known. They then add nonconservative solenoidal forces
that couple the two systems but chosen such that $P_{\text{\tiny B}}$ remains
the stationary distribution of the resulting dissipative dynamics for each
individual system. The strength $\gamma$ of the coupling between the two
copies of the system is a free parameter. As expected on general grounds, the
relaxation time to the stationary state is shorter in the presence of these
peculiar forces. They have illustrated their findings on several simple
examples. A key open question which we address in this work is to quantify the
extent to which convergence to the Boltzmann distribution is accelerated in
the presence of a single energy barrier, and in the presence of a large number
of barriers in the energy landscape of the system of interest.

To answer these questions, we begin by asking how the IO dynamics behaves in
the presence of a single energy barrier. This is a special case of the general
study conducted by Bouchet and Reygner [28] in which they determined a formal
expression for the modified Kramers escape time. In addition to specializing
the results of [28] we give a heuristic argument and we compare how the
optimal trajectory deviates from the standard equilibrium counterpart. Then we
implement the IO dynamics on a mean-field disordered $p$-spin system whose
energy landscape is famous for harboring a very large number of metastable
basins. We chose the disordered $p$-spin for the following reason: below a so-
called dynamical temperature $T_{d}$, its equilibrium dynamics freezes and the
resulting nonergodic behavior prevents the system from correctly sampling its
target equilibrium Boltzmann distribution, so that it can be viewed as the
most challenging system to test new ideas in optimization and sampling
acceleration. We find that the dynamics is indeed accelerated, at various
stages of the dynamics, though ergodicity breaking occurs at the same
dynamical temperature.

##  II The Ichiki-Ohzeki dynamics

Let us begin by considering a particle with position r that evolves according
to an overdamped Langevin dynamics in some energy landscape $V(\text{\bf r})$
at temperature $T$. The equation of motion reads

$\dot{\text{\bf r}}=-\bm{\nabla}_{\text{\bf r}}V(\text{\bf
r})+\sqrt{2T}\bm{\eta}$ (1)

where $\bm{\eta}$ is a Gaussian white noise with unit correlations
$\langle\bm{\eta}(t)\otimes\bm{\eta}(t^{\prime})\rangle=\mathds{1}\delta(t-t^{\prime})$.
The companion Fokker-Planck equation

$\partial_{t}P=-\bm{\nabla}_{\text{\bf r}}\cdot{\bf j},\,\,{\bf j}(\text{\bf
r},t)=-T\bm{\nabla}_{\text{\bf r}}P-\bm{\nabla}_{\text{\bf r}}VP$ (2)

has the Boltzmann distribution $P_{\text{\tiny B}}(\text{\bf
r})\propto\text{e}^{-\beta V(\text{\bf r})}$ as its stationary solution.
Because of the underlying detailed balance condition, this stationary solution
$P_{\text{\tiny B}}$ is actually an equilibrium one, as it makes the
probability current $\bf j$ vanish, and no entropy is produced in the
equilibrium steady-state.

Consider now two identical particles with positions $\text{\bf r}^{(1)}$ and
$\text{\bf r}^{(2)}$ evolving according to the coupled Langevin dynamics with
independent white Gaussian noises $\bm{\eta}^{(1)}$ and $\bm{\eta}^{(2)}$,

$\displaystyle\dot{\text{\bf r}}^{(1)}=-\bm{\nabla}_{\text{\bf
r}^{(1)}}V(\text{\bf r}^{(1)})+\gamma\bm{\nabla}_{\text{\bf
r}^{(2)}}V(\text{\bf r}^{(2)})+\sqrt{2T}\bm{\eta}^{(1)}$ (3)
$\displaystyle\dot{\text{\bf r}}^{(2)}=-\bm{\nabla}_{\text{\bf
r}^{(2)}}V(\text{\bf r}^{(2)})-\gamma\bm{\nabla}_{\text{\bf
r}^{(1)}}V(\text{\bf r}^{(1)})+\sqrt{2T}\bm{\eta}^{(2)}$ (4)

The parameter $\gamma$ introduces a skew symmetric coupling between the two
systems. Two independent overdamped Langevin dynamics for systems $1$ and $2$
are recovered when $\gamma=0$. Note that when $\gamma\neq 0$ the zero
temperature dynamics of a single system is not a gradient descent down its
potential landscape. Yet a gradient descent dynamics for the composite system
is recovered, since by denoting $U(\text{\bf r}^{(1)},\text{\bf
r}^{(2)})\equiv V(\text{\bf r}^{(1)})+V(\text{\bf r}^{(2)})$ we have

$\frac{\differential}{\differential t}U(\text{\bf r}^{(1)},\text{\bf
r}^{(2)})=-\sum_{a=1,2}\left(\bm{\nabla}_{\text{\bf r}^{(a)}}V(\text{\bf
r}^{(a)})\right)^{2}<0$ (5)

Therefore, should one of the subsystems get stuck in some local minimum, then
the other one is able to temporarily increase its energy. The latter system
then acts as a source of fluctuations for the stuck system. At nonzero
temperature, these fluctuations add to the thermal ones, allowing for a faster
exploration of the configuration space.

Let’s verify that the nonequilibrium stationary state measure is the Boltzmann
measure. The probability current in the Fokker-Planck equation for
$P(\text{\bf r}^{(1)},\text{\bf r}^{(2)},t)$ follows from Eq. (3). With
respect to the expression of Eq. (2) picks up an additional contribution for
each of the two subspaces, namely

$\begin{split}\delta{\bf j}^{(1)}=\gamma\bm{\nabla}_{\text{\bf
r}^{(1)}}V(\text{\bf r}^{(2)})P(\text{\bf r}^{(1)},\text{\bf r}^{(2)},t)\\\
\delta{\bf j}^{(2)}=-\gamma\bm{\nabla}_{\text{\bf r}^{(2)}}V(\text{\bf
r}^{(1)})P(\text{\bf r}^{(1)},\text{\bf r}^{(2)},t)\end{split}$ (6)

This additional current arises from the presence of the nonconservative forces
in Eq. (3), and it is divergence-free when $P(\text{\bf r}^{(1)},\text{\bf
r}^{(2)},t)=P_{\text{\tiny B}}(\text{\bf r}_{1})P_{\text{\tiny B}}(\text{\bf
r}_{2})$. We therefore know that this factored Boltzmann distribution is the
stationary state of the coupled process. However the dynamical rules in Eq.
(3) are accompanied by a positive entropy production rate $\dot{\Sigma}$ given
by

$\dot{\Sigma}=2\beta\gamma^{2}\langle(\bm{\nabla}V)^{2}\rangle_{\text{\tiny
B}}=2\gamma^{2}\langle\Delta V\rangle_{\text{\tiny B}}>0$ (7)

where $\langle\ldots\rangle_{\text{\tiny B}}$ refers to an average with
respect to the Boltzmann distribution, thus confirming the nonequilibrium
nature of the dynamics, and that $P_{\text{\tiny B}}$ is a nonequilibrium
steady-state distribution. This calculation tells us about the typical time
scale $1/\dot{\Sigma}$ it takes for the dynamics to display irreversible
features.

We now give a physicist’s argument explaining why the IO dynamics is faster
than the standard Langevin dynamics of Eq. (1).

###  II.1 A simple argument for the decrease of the relaxation time

A mathematical proof accounting for the decrease of the relaxation time for
arbitrary $\gamma$ was given in [14, 15] and it was rediscovered in the case
of discrete state Markov processes in [17]. However, for small $\gamma$, there
is a simple physicist’s argument to account for the increase in the relaxation
rate towards the steady state of the Fokker-Planck operator. To see this we
use the Darboux [29] and we rewrite the Fokker-Planck operator $\mathbb{W}$
associated to Eq. (3) in a rotated basis, $\tilde{\mathbb{W}}=P_{\text{\tiny
B}}^{1/2}\mathbb{W}P_{\text{\tiny B}}^{-1/2}$, which reads

$\tilde{\mathbb{W}}=\tilde{\mathbb{W}}_{0}^{(1)}+\tilde{\mathbb{W}}_{0}^{(2)}+\delta\tilde{\mathbb{W}}$
(8)

where

$\begin{split}\tilde{\mathbb{W}}_{0}^{(a)}&\equiv T\bm{\nabla}_{\text{\bf
r}^{(a)}}^{2}\\\ &-\frac{1}{4T}\left[(\bm{\nabla}_{\text{\bf
r}^{(a)}}V(\text{\bf r}^{(a)}))^{2}-2T\bm{\nabla}^{2}_{\text{\bf
r}^{(a)}}V(\text{\bf r}^{(a)})\right]\end{split}$ (9)

is the symmetrized Fokker-Planck operator for particle $a=1,2$ in the absence
of the coupling parameter $\gamma$ and

$\delta\tilde{\mathbb{W}}\equiv\gamma\left[\left(\bm{\nabla}_{\text{\bf
r}_{1}}V(\text{\bf r}_{1})\right)\cdot\bm{\nabla}_{\text{\bf
r}_{2}}-\left(\bm{\nabla}_{\text{\bf r}_{2}}V(\text{\bf
r}_{2})\right)\cdot\bm{\nabla}_{\text{\bf r}_{1}}\right]$ (10)

As expected, for $\gamma=0$ this is a Hermitian operator (which reflects the
detailed balance condition), but for $\gamma\neq 0$ a skew-Hermitian
perturbation is introduced in Eq. (8). We denote by $\ket{\Psi_{n}}$ the
eigenstates of the unperturbed Hermitian operator
$\tilde{\mathbb{W}}_{0}^{(1)}+\tilde{\mathbb{W}}_{0}^{(2)}$, and with $E_{n}$
the associated eigenvalue. We are interested in the change $\delta E$ to the
spectral gap $\Delta E=E_{1}-E_{0}$ of the operator $-\tilde{\mathbb{W}}$
between the ground state $\ket{\Psi_{0}}$ and the first excited state
$\ket{\Psi_{1}}$. We first observe that the ground state and its (zero) energy
are left unchanged by the perturbation. Therefore, to second order in
perturbation theory, we have

$\begin{split}\Delta
E&=-\gamma\matrixelement{\Psi_{1}}{\delta\tilde{\mathbb{W}}}{\Psi_{1}}\\\
&-\gamma^{2}\sum_{n>1}\frac{\absolutevalue{\matrixelement{\Psi_{1}}{\delta\tilde{\mathbb{W}}}{\psi_{n}}}^{2}}{E_{1}-E_{n}}\end{split}$
(11)

Due to the skew-Hermiticity of the perturbation the first term in Eq. (11) is
imaginary, and the second term is real and positive. To second order in
$\gamma$, the eigenvalues of $\tilde{\mathbb{W}}$ are shifted in the opposite
direction to what would occur in the case of a Hermitian perturbation. Only
the real part of the spectral gap is relevant in the computation of the
relaxation times, and since it is increased by the perturbation, the
relaxation time toward the steady state distribution is reduced.

This increase in the relaxation rate has been explored for simple systems in
[18, 19, 30], but for systems whose dynamics is controlled by rare—yet
important—events, the efficiency of the IO approach deserves deeper
investigation. As a first step, we revisit Kramers’ escape problem with the IO
dynamics Eq. (3) to understand qualitatively how it copes with a single energy
barrier.

##  III The barrier crossing problem

Consider a one-dimensional potential landscape $V(x)$ with a stable
equilibrium position at $x_{m}$ and ask about the time $\tau$ it takes for
particle $1$ to hop over a potential barrier $\Delta V$, passing through the
local maximum at $x_{M}$, starting from the bottom of the well. As we know
from Kramers’ analysis [31], for equilibrium dynamics, $\tau$ is given by the
Arrhenius formula $\tau\simeq\text{e}^{\beta\Delta V}$. Consider now a nonzero
$\gamma$ for two coupled particles. Both particles start from the bottom of
the potential well. We ask how long it takes particle 1 to hop over the
potential barrier irrespective of the motion of particle 2. We begin our
presentation with a heuristic analysis before connecting to a recent work by
Bouchet and Reygner [28].

###  III.1 Heuristic argument

Before proceeding with a detailed analysis, there is a simple way to
understand how a nonzero $\gamma$ speeds barrier crossing up. We are
interested in the rare event in which particle $1$ will hop over an energy
barrier while, at low temperature, the typical motion of particle 2 remains
confined within the bottom of the well. The position of particle 2 then
evolves according to

$\begin{split}\dot{x}^{(2)}(t)=&-V^{\prime}(x^{(2)}(t))-\gamma
V^{\prime}(x^{(1)}(t))+\sqrt{2T}\eta^{(2)}(t)\\\
\simeq&-k(x^{(2)}(t)-x_{m})-\gamma
V^{\prime}(x^{(1)}(t))+\sqrt{2T}\eta^{(2)}(t)\end{split}$ (12)

where $k_{m}\equiv V^{\prime\prime}(x_{m})>0$ is the stiffness of the well. We
can integrate this equation and substitute $x^{(2)}$ into the evolution
equation for $x^{(1)}$:

$\begin{split}\dot{x}^{(1)}(t)=&-\int_{0}^{t}\differential
sM_{R}(t-s)V^{\prime}(x^{(1)}(s))+\sqrt{2T}\eta(t)\\\ \end{split}$ (13)

where $\eta(t)$ is now a colored Gaussian noise with correlations
$\langle\eta(t)\eta(t^{\prime})\rangle=M_{C}(t-t^{\prime})$. The mobility
kernel $M_{R}(t)=\delta(t)+\gamma^{2}k_{m}\text{e}^{-k_{m}t}\theta(t)$ and the
correlation kernel $M_{C}(t)=\frac{1}{2}M_{R}(t)+\frac{1}{2}M_{R}(-t)$ are
related by Kubo’s fluctuation theorem for generalized Langevin equations.
While the two-dimensional process $(x^{(1)},x^{(2)})$ is a nonequilibrium one,
the effective process $x^{(1)}$, forgetting the information on $x^{(2)}$, is a
truly equilibrium process (at least when $x^{(2)}$ evolves in a quadratic
well). In the large $k_{m}$ limit, $M_{R}(t)\simeq(1+\gamma^{2})\delta(t)$ and
the mobility of the particle picks up a $\gamma^{2}$ contribution so that the
Kramers’ escape time is reduced accordingly: $\tau\simeq\text{e}^{\beta\Delta
V}/(1+\gamma^{2})$.

This heuristic approach suggests that the IO dynamics contributes to changes
in the barrier crossing time by rescaling the prefactor of the Arrhenius
formula. We will confirm this observation in the general case by looking at
the most probable path taken by the system during the barrier crossing
process, also known as the instanton trajectory.

To do so, it is convenient to cast the equations of motion for the two coupled
one-dimensional processes in a matrix form

$\mathbf{\dot{x}}=-\left(\mathds{1}+\bm{\Gamma}\right)\cdot\bm{\nabla}U(\mathbf{x})+\sqrt{2T}\bm{\eta}(t)$
(14)

with $\mathbf{x}\equiv\begin{bmatrix}x^{(1)}&x^{(2)}\end{bmatrix}^{T}$,
$U(\mathbf{x})\equiv V(x^{(1)})+V(x^{(2)})$ and

$\bm{\Gamma}\equiv\begin{bmatrix}0&-\gamma\\\ \gamma&0\end{bmatrix}$ (15)

This two-dimensional system presents a local minimum at
$\mathbf{x}_{m}\equiv\begin{bmatrix}x_{m}&x_{m}\end{bmatrix}^{T}$ and a saddle
point at $\mathbf{x}_{s}\equiv\begin{bmatrix}x_{M}&x_{m}\end{bmatrix}^{T}$.
The process in which particle $1$ overcomes the energetic barrier while
particle $2$ remains confined in the well corresponds, for the composite
system, to a trajectory from $\mathbf{x}_{m}$ to $\mathbf{x}_{s}$ and from
$\mathbf{x}_{s}$ to the minimum on the other side of the saddle. We now
proceed with a careful analysis of the instanton trajectory.

###  III.2 Instanton trajectory: no change in the exponential factor

The Kramers escape rate $\tau^{-1}$ can be expressed within a path-integral
formulation,

$\tau^{-1}=\int{\mathscr{D}}\mathbf{\hat{x}}{\mathscr{D}}\text{\bf
x}\text{e}^{-\frac{1}{T}S\left[\mathbf{x},\mathbf{\hat{x}}\right]}$ (16)

where the path-integral runs over paths connecting $\mathbf{x}_{m}$ to
$\mathbf{x}_{s}$, and where the Janssen-De Dominicis action functional reads

$S\left[\mathbf{x},\mathbf{\hat{x}}\right]=\int_{-\infty}^{+\infty}\differential
t\left[\mathbf{\dot{x}}\cdot\mathbf{\hat{x}}-H\left(\mathbf{x},\mathbf{\hat{x}}\right)\right]$
(17)

The response field $\mathbf{\hat{x}}$ expresses the effect of the noise on the
dynamics of $\mathbf{x}$. In the small noise $T\rightarrow 0$ limit, the path
that dominates the integral minimizes $S$, and the resulting saddle-point
trajectories follow a Hamiltonian-like dynamics (with respect to $H$ in Eq.
(17)) in a space of extended variables where $\mathbf{\hat{x}}$ stands for a
momentum. The escape time $\tau$ then satisfies the logarithmic equivalence

$\begin{cases}\tau\asymp\text{e}^{\frac{\phi}{T}}\\\
\phi\equiv\min_{\mathbf{x}(t)|\mathbf{x}(-\infty)=\mathbf{x}_{m}\text{,
}\mathbf{x}(+\infty)=\mathbf{x}_{s}}S\left[\mathbf{x},\mathbf{\hat{x}}\right]\end{cases}$
(18)

The explicit expression for the Hamiltonian $H(\mathbf{x},\mathbf{\hat{x}})$
is

$H(\mathbf{x},\mathbf{\hat{x}})\equiv-\mathbf{\hat{x}}\cdot\left(\mathds{1}+\bm{\Gamma}\right)\cdot\bm{\nabla}U(\mathbf{x})+\mathbf{\hat{x}}^{2}$
(19)

The minimum of $S$ for the trajectories of interest is found by looking for
the solutions of the Hamiltonian system induced by $H$ along the $H=0$ energy
manifold, so that we have the differential equations

$\begin{cases}\mathbf{\dot{x}}&=\frac{\partial
H}{\partial\mathbf{\hat{x}}}=-\left(\mathds{1}+\bm{\Gamma}\right)\cdot\bm{\nabla}U(\mathbf{x})+2\mathbf{\hat{x}}\\\
\mathbf{\dot{\hat{x}}}&=-\frac{\partial
H}{\partial\mathbf{x}}=\left(\mathds{1}+\mathbf{\Gamma}\right)\cdot\text{Hess
}U(\mathbf{x})\cdot\mathbf{\hat{x}}\end{cases}$ (20)

together with the condition $H(\mathbf{x},\mathbf{\hat{x}})=0$. Here
$\text{Hess }U(\mathbf{x})_{ab}\equiv V^{\prime\prime}(x^{(a)})\delta_{ab}$ is
the Hessian matrix of $U$. These coupled differential equations are equivalent
to

$\begin{cases}\mathbf{\dot{x}}&=\left(\mathds{1}-\mathbf{\Gamma}\right)\cdot\mathbf{\bm{\nabla}}U(\mathbf{x})\\\
\mathbf{\hat{x}}&=\mathbf{\bm{\nabla}}U(\mathbf{x})\end{cases}$ (21)

which, in principle, can be solved explicitly if a potential $V$ is specified.
We solve the instanton trajectory $\dot{\mathbf{x}}$ within perturbation
theory up to second order in $\gamma$ as a functional of the external
potential $V$, assuming the $\gamma=0$ trajectory is fully known. A linear
analysis of the trajectory close to the saddle point $\mathbf{x}_{s}$ shows
(see App. A) that the perturbation expansion must have the form

$\displaystyle x^{(1)}$
$\displaystyle=x^{(1)}_{0}+\gamma^{2}x^{(1)}_{2}+\mathcal{O}(\gamma^{4})$ (22)
$\displaystyle x^{(2)}$ $\displaystyle=x^{(2)}_{0}+\gamma
x^{(2)}_{1}+\mathcal{O}(\gamma^{3})$

with $\dot{x}^{(1)}_{0}=V^{\prime}\left(x^{(1)}_{0}(t)\right)$ and
$x^{(2)}_{0}=0$. By taking into account the boundary condition for the
perturbation $x(+\infty)^{(a)}_{i}=x(-\infty)^{(a)}_{i}=0$ for $a=1,2$, $i>0$
we obtain

$\displaystyle x^{(2)}_{1}$
$\displaystyle=-\text{e}^{k_{m}t}\int_{t}^{+\infty}\differential\tau\;\text{e}^{-k_{m}\tau}V^{\prime}\left(x^{(1)}_{0}(\tau)\right)$
(23) $\displaystyle x^{(1)}_{2}$
$\displaystyle=k_{m}\text{e}^{G(t)}\left[\int_{t}^{+\infty}\differential\tau\;\text{e}^{-G(\tau)}x^{(2)}_{1}(\tau)+C\right]$
(24)

where $G(t)\equiv\int\differential\tau V^{\prime\prime}(x^{(1)}_{0}(\tau))$
and $C$ is a constant determined by matching the asymptotic expression for
$x^{(1)}_{2}(t)$ for $t\rightarrow+\infty$ to the linearized expansion at
$\mathbf{x}_{s}$. We illustrate the role played by $\gamma$ in deforming the
equilibrium trajectory in Fig. 1. As done in previous studies of the escape
time out of equilibrium [32], we have chosen as an explicit potential $V$ the
well-known double well $V(x)\equiv\frac{1}{4}x^{4}-\frac{1}{2}x^{2}$. We see
that indeed system 2 finds it optimal to rise in its own energy landscape to
favor the crossing of system 1 over the barrier.

Figure 1: Instanton trajectory from Eq. (21) for $\gamma=0$ and its
perturbative expansion up to second order in $\gamma$ for the double well
potential $V(x)=\frac{1}{4}x^{4}-\frac{1}{2}x^{2}$.

Interestingly, the action for the instanton trajectory is unchanged with
respect to the equilibrium dynamics, namely

$S\left[\mathbf{x},\mathbf{\hat{x}}\right]=\int_{-\infty}^{+\infty}\differential
t\mathbf{\dot{x}}\cdot\mathbf{\bm{\nabla}}U(\mathbf{x})=\Delta U$ (25)

where $\Delta U=V(x_{M})-V(x_{m})$. That this result is identical to the one
obtained with the equilibrium dynamics at $\gamma=0$ supports the the
heuristic analysis of Sec. III.1. While the instanton trajectory is clearly
altered, the Arrhenius factor is left unchanged by the IO dynamics. In the
next section, we will investigate the role of $\gamma$ on the prefactor to the
Arrhenius exponential.

###  III.3 Full expression of $\tau$

Bouchet and Reygner [28] have derived an Eyring-Kramers formula for
nonequilibrium dynamics which do not necessarily sample the Boltzmann
distribution. Here we specialize their general result to the IO dynamics,
which allows us to arrive at a concrete $\gamma$-dependence of the escape
time. When applied to the IO dynamics the time $\tau$ reads

$\tau=\frac{2\pi}{\lambda_{+}}\sqrt{\frac{k_{M}}{k_{m}}}\exp\left(\frac{\Delta
V}{T}\right)$ (26)

where $k_{M}\equiv-V^{\prime\prime}({x}_{M})>0$, $k_{m}\equiv
V^{\prime\prime}(x_{m})>0$, and $\lambda_{+}$ is the modulus of the eigenvalue
related to the unstable direction of the dynamics linearized around the saddle
point $\mathbf{x}_{s}$:

$\lambda_{+}=\frac{1}{2}\left[k_{M}-k_{m}+\sqrt{\left(k_{M}+k_{m}\right)^{2}+4\gamma^{2}k_{M}k_{m}}\right]$
(27)

For $\gamma=0$ Eq. (26) yields the well-known Kramers’ time

$\tau=\frac{2\pi}{\sqrt{k_{M}k_{m}}}\exp\left(\frac{\Delta V}{T}\right)$ (28)

while for $\gamma\neq 0$ we have $\lambda_{+}>k_{M}$ and the time required for
particle $1$ to hop over the energetic barrier is reduced.

We realize a posteriori that the heuristic approach of Sec. III.1 was
physically sound. It can indeed be recovered by letting
$k_{m}\rightarrow+\infty$ in Eq. (27), yielding $\lambda_{+}\approx
k_{M}(1+\gamma^{2})$. This limiting case is actually the smallest achievable
value of $\tau$ at given $\gamma$, $k_{M}$ and $\Delta V$, since
$\partial_{k_{m}}\lambda_{+}$ is always positive.

The results of this section precisely quantify the gain in speed that $\gamma$
leads to in the crossing of a single steep energy barrier. A natural question
that arises next is whether any similar gain –to be quantified as a function
of $\gamma$— can be harvested in the presence of a large number of steep
energy barriers. We do not expect the IO dynamics to overcome the ergodicity
breaking that occurs in some systems where the energy barriers become
infinitely high, as is for instance the case in some mean-field disordered
system in the large system size limit. In what follows, our goal is to put the
IO dynamics to a severe test. We shall now investigate how it performs on a
specific $N$-body system known for displaying strong ergodicity breaking
behavior at low temperatures, the disordered $p$-spin model for $p\geq 3$.

##  IV $p$-spin with Ichiki-Ohzeki dynamics

The spherical $p$-spin is a disordered fully connected model of $N$
interacting spins $\sigma_{1},\ldots,\sigma_{N}$ fulfilling the constraint
$\sum_{i=1}^{N}\sigma_{i}^{2}=N$. The energy of a configuration is

${\mathscr{H}}[\\{\sigma_{i}\\}_{i=1,\ldots,N}]=-\sum_{i_{1}<\ldots<i_{p}}J_{i_{1}\ldots
i_{p}}\sigma_{i_{1}}\ldots\sigma_{i_{p}}$ (29)

where the $J_{i_{1}\ldots i_{p}}$ are independent quenched random couplings
with variance $\frac{2J^{2}}{p!N^{p-1}}$. To cut a long story short, we refer
the reader to the huge body of literature [33, 34, 35, 36, 37, 38] devoted to
the $p$-spin to understand the role it has played in shaping up our
understanding of disordered and glassy systems (including active ones [39]). A
summary of the properties relevant to the present discussion is that this
model undergoes a transition from a paramagnetic state at high temperatures to
a spin-glass phase at low temperatures $T<T_{c}$ [33]. From a static point of
view, its peculiarity is that there exists an intermediate temperature
$T_{d}>T_{c}$ such that in the temperature range $T_{c}<T<T_{d}$ an
exponentially large number of metastable states appears, without however
affecting the paramagnetic nature of the system. From a dynamical standpoint,
while the system relaxes above $T_{d}$, it fails to explore this complex
energy landcsape as soon as $T<T_{d}$ [34]. When endowed with purely
relaxational equilibrium dynamics, the spin-spin correlation function
$C(t,t^{\prime})=\frac{1}{N}\sum_{i}\langle\sigma_{i}(t)\sigma_{i}(t^{\prime})\rangle$,
in the regime with $t^{\prime}$ large and $t-t^{\prime}$ finite with time-
translation invariance, fails to relax to $0$ for $T<T_{d}$, reaching instead
a temperature-dependent plateau value $q$. The $p$-spin is thus the ideal
test-bench to quantify the acceleration of the IO dynamics.

###  IV.1 Mean field IO dynamics

In order to implement the IO dynamics we consider two coupled $p$-spin systems
(with the same realization of the disorder), with spins $\sigma_{i}^{(1)}$ and
$\sigma_{i}^{(2)}$ evolving according to

$\displaystyle\begin{split}\partial_{t}\sigma_{i}^{(1)}&=-\frac{\partial{\mathscr{H}}}{\partial\sigma_{i}^{(1)}}+\gamma\frac{\partial{\mathscr{H}}}{\partial\sigma_{i}^{(2)}}\\\
&-\mu^{(1)}(t)\sigma_{i}^{(1)}+\gamma\mu^{(2)}(t)\sigma_{i}^{(2)}+\sqrt{2T}\eta_{i}^{(1)}\end{split}$
(30)
$\displaystyle\begin{split}\partial_{t}\sigma_{i}^{(2)}&=-\frac{\partial{\mathscr{H}}}{\partial\sigma_{i}^{(2)}}-\gamma\frac{\partial{\mathscr{H}}}{\partial\sigma_{i}^{(1)}}\\\
&-\mu_{2}(t)\sigma_{i}^{(2)}-\gamma\mu^{(1)}(t)\sigma_{i}^{(1)}+\sqrt{2T}\eta_{i}^{(2)}\end{split}$
(31)

where the functions $\mu^{(a)}(t)$ are adjusted so that the spherical
constraint
$\sum_{i=1}^{N}\left\langle\left(\sigma_{i}^{(a)}(t)\right)^{2}\right\rangle=N$
is fulfilled throughout the dynamical evolution. There is a large body of
literature on the effect of a nonequilibrium drive on the $p$-spin (or on
other complex systems with a similar glassy behavior at low temperatures). In
these works it has consistently been found that in the presence of a
nonequilibrium drive (be it in a randomized non-relaxational dynamics [40], or
in a shear [41] or with an active force [39]) a lower temperature than the
equilibrium one is required in order to witness the onset of ergodicity
breaking. However, in those cases the nonequilibrium stationary-state that is
being sampled is not the Boltzmann distribution, and it is in general unknown.
We stress that, unlike in the latter series of works, the IO dynamics does
correctly sample, in its stationary state, the Boltzmann distribution
$P_{\text{\tiny B}}\sim\text{e}^{-\beta{\mathscr{H}}}$ for
${\mathscr{H}}={\mathscr{H}}[\\{\sigma_{i}^{(1)}\\}]+{\mathscr{H}}[\\{\sigma_{i}^{(2)}\\}]$
(with ${\mathscr{H}}$ given in Eq. (29)), despite being a nonequilibrium
dynamics. From an analytical standpoint, a remarkable feature of the standard
dynamics of the $p$-spin, which is preserved by Eq. (30) is that the mean-
field nature of the problem allows us to write the stochastic evolution of a
single spin in terms of unknown functions (the spin-spin correlation, the
response function that enters the memory kernel, and noise correlations) that
are determined self-consistently. Such a procedure can be applied to our case.
Note however that in the context of the IO dynamics, the spin-spin correlation
becomes a collection of four spin-spin correlations
$C_{ab}(t,t^{\prime})=\langle\sigma^{(a)}(t)\sigma^{(b)}(t^{\prime})\rangle$,
and the same holds for the response function $R_{ab}(t,t^{\prime})$ and the
integrated response
$F_{ab}(t,t^{\prime})\equiv-\int_{t^{\prime}}^{t}\differential\tau
R_{ab}(t,\tau)$ . In practice one is able to write an effective Langevin
equation for the dynamics of a single spin $\sigma^{(1)}$ in system $1$ and
its counterpart $\sigma^{(2)}$ in system $2$, from which self-consistent
equations of motion for the correlation and the integrated response can be
derived. Since our interest goes to the dynamics in the steady state we assume
time translational invariance to hold,
$C_{ab}(t,t^{\prime})=C_{ab}(t-t^{\prime})$ and
$F_{ab}(t,t^{\prime})=F_{ab}(t-t^{\prime})$, and by applying standard path
integral techniques we arrive at the following Langevin equation for
$\bm{\sigma}(t)\equiv\begin{bmatrix}\sigma^{(1)}&\sigma^{(2)})\end{bmatrix}^{T}$
and the dynamical equations for the correlation matrix $\mathbf{C}$ and the
integrated response matrix $\mathbf{F}$:

$\displaystyle\partial_{t}{\bm{\sigma}}(t)$
$\displaystyle=-\left(\mathds{1}+\bm{\Gamma}\right)\cdot{\bm{\mu}}(t)\cdot{\bm{\sigma}}(t)+\int_{0}^{t}\differential\tau\mathbf{M_{F}}(t-\tau)\cdot{\bm{\sigma}}(\tau)+\frac{1}{T}\left(\mathds{1}+\bm{\Gamma}\right)\cdot\mathbf{C}_{0}(p,t)\cdot{\bm{\sigma}}(0)+{\bm{\xi}}(t)$
(32) $\displaystyle\partial_{t}\mathbf{F}(t)$
$\displaystyle=-\mathds{1}-\left(\mathds{1}+\bm{\Gamma}\right)\cdot\bm{\mu}(t)\cdot\mathbf{F}(t)+\int_{t^{\prime}}^{t}d\tau\mathbf{M_{F}}(t-\tau)\cdot\mathbf{F}(\tau)$
(33) $\displaystyle\partial_{t}\mathbf{C}(t)$
$\displaystyle=-\left(\mathds{1}+\bm{\Gamma}\right)\cdot\bm{\mu}(t)\cdot\mathbf{C}(t)+\int_{0}^{t}d\tau\mathbf{M_{F}}(t-\tau)\cdot\mathbf{C}(\tau)+\frac{1}{T}\left(\mathds{1}+\bm{\Gamma}\right)\cdot\mathbf{C}_{0}(p,t)$
(34)
$\displaystyle\begin{split}\bm{\mu}(t)&=\text{Diag}\Biggl{[}T\mathds{1}_{2}+\int_{0}^{t}d\tau\mathbf{M}(t-\tau)\cdot\mathbf{C^{T}}(t-\tau)+\int_{0}^{t}d\tau\mathbf{D}(t-\tau)\cdot\partial_{\tau}\mathbf{F^{T}}(t-\tau)\\\
&+\frac{1}{T}\left(\mathds{1}+\bm{\Gamma}\right)\cdot\mathbf{C}_{0}(p,t)\cdot\mathbf{C^{T}}(t)\Biggr{]}\end{split}$
(35)

where the Gaussian noise has correlations
$\langle\bm{\xi}(t)\otimes\bm{\xi}(t^{\prime})^{T}\rangle=2T{\bf
1}_{2}\delta(t-t^{\prime})+{\bf D}(t,t^{\prime})$. The three kernels ${\bf
M_{F}}$, ${\bf C}_{0}$ and ${\bf D}$ that appear in Eqs. Eqs. (32, 33, 34) are
nonlinear functions of the elements of $\mathbf{C}$ and $\mathbf{F}$. The
overall structure of these equations is similar to the equilibrium setting,
though their derivation is somewhat more involved and we have deferred it to
App. B. These equations are solved with the boundary conditions
$\mathbf{C}(0)=\mathds{1}$, $\mathbf{F}(0)=\mathbf{0}$, while the explicit
expression for the diagonal matrix $\bm{\mu}(t)$ enforces the spherical
constraint.

In the case of standard equilibrium dynamics, the response and the
correlations are related by the fluctuation-dissipation theorem (FDT),
$\partial_{t}\mathbf{F}(t)=-\frac{1}{T}\partial_{t}\mathbf{C}(t)$ and it is
possible to obtain a closed equation for the spin-spin correlation function of
the system alone. By focusing on the long time limit of the dynamics one can
then find the temperature below which the correlation function displays a
nonzero plateau at infinite times when ergodicity is broken. Since the IO
dynamics is a nonequilibrium one, we cannot resort to the fluctuation
dissipation theorem, and analytical manipulations of Eqs. (32, 33, 34, 35) may
seem, at first sight, out of reach. In the following section, however, we will
exhibit a simple argument that suggests that for the specific correlations and
integrated response considered above, a relation formally similar to the FDT
might still survive in the nonequilibrium steady-state.

###  IV.2 Response and fluctuations

####  IV.2.1 Motivations

The entropy production is a standard way of measuring the distance to
equilibrium. Here we find, using Eq. (7), that

$\dot{\Sigma}=2\gamma^{2}N\mu$ (36)

so that the time scale at which nonequilibrium effects come into play goes as
$\frac{1}{\gamma^{2}}$. It is tempting to believe that the distance to
equilibrium can generically be viewed as being of order $O(\gamma^{2})$. This
may be true, but the connection between the response and the correlation
functions, in our particular systems, runs almost as deep as in equilibrium.

Before we address this question, it is useful, as a preliminary exercise, to
work out the example of a harmonic oscillator undergoing the IO dynamics. This
is a bit of an artificial example, as coarse-graining out one the two
oscillators leads to an effective equilibrium dynamics for the remaining one
(with memory, as discussed in Sec. III.1). Consider the two coupled equations
for $\text{\bf r}=(x_{1},x_{2})$

$\dot{\text{\bf r}}=-(\mathds{1}+{\mathbf{\Gamma}})k\text{\bf
r}+\sqrt{2T}\beta$ (37)

The evolution operator has a spectrum of relaxation rates
$-k(n_{1}+n_{2}+i\gamma(n_{1}-n_{2}))$ (with $n_{1},n_{2}\in\mathbb{N}$) and
we realize that harmonic oscillators are marginal with respect to the
acceleration of their dynamics: the relaxation rates only pick up a nonzero
imaginary part expressing an overall rotation in phase space, but no
modification of the real part means net acceleration. This system of two
coupled oscillators is however indeed out of equilibrium since
$\dot{\Sigma}=2k\gamma^{2}>0$. The correlation matrix
$\mathbf{C}(t,t^{\prime})=\langle\text{\bf r}(t)\otimes\text{\bf
r}(t^{\prime})\rangle$ and the response function
${\mathbf{R}}(t,t^{\prime})=\left.\frac{\delta\langle\text{\bf
r}(t)\rangle}{\delta{\bf f}(t^{\prime})}\right|_{\bf=\text{\bf 0}}$ to an
external perturbing force $\mathbf{f}$ are related, in the stationary state,
by

$\mathbf{R}(t,t^{\prime})(\mathds{1}+\bm{\Gamma})=-\frac{1}{T}\partial_{t}{\mathbf{C}}(t,t^{\prime})$
(38)

which can be cast in an integral form by introducing the integrated response
$\mathbf{F}(t,t^{\prime})\equiv-\int_{t^{\prime}}^{t}\differential\tau\mathbf{R}(t,\tau)$
as:

$\mathbf{F}(t,t^{\prime})(\mathds{1}+\bm{\Gamma})=\frac{1}{T}\left({\mathbf{C}}(t,t^{\prime})-\mathbf{C}(t^{\prime},t^{\prime})\right)$
(39)

The fluctuation-dissipation theorem is not fulfilled, not even in some
modified form, but it is true that the specific position-position correlation
function is directly related to the integrated response in a linear fashion.
There is no such generic connection for arbitrary observables. For lack of a
better name, we christen this connection between these specific response and
the correlations as an accidental FDT (aFDT). We shall now turn to a numerical
investigation of the connection between $\mathbf{F}$ and $\mathbf{C}$ for the
$p$-spin. We keep the a priori unrelated result of Eq. (38) under our belt for
the analysis of the numerics which we carry out in the next subsection.

####  IV.2.2 Numerics and aFDT

We integrate the dynamical equations Eqs. (33, 34, 35) for $p=3$ using a
finite difference method implemented in [42, 43, 44, 45] and detailed in App.
C. We check that our implementation of the IO dynamics samples the correct
Boltzmann distribution by measuring the energy per spin $E_{\infty}$ at long
times in the ergodic phase. This is a static quantity independent from
$\gamma$, for which a value of $E_{\infty}=-\frac{1}{2T}$ is predicted [33].
In the numerical integration, we computed the energy by exploiting the
relation

$\bm{\mu}(t\rightarrow+\infty)=\left(T-pE_{\infty}\right)\mathds{1}$ (40)

The agreement between the measured value of $E_{\infty}$ and its theoretical
prediction is shown in Fig. 2. The errors are of order $10^{-12}$ and they are
therefore negligible.

Figure 2: Energy per spin as a function of temperature $T$ obtained from
numerical integration (blue dots) at $\gamma=1$ and theoretical prediction
$E_{\infty}=-\frac{1}{2T}$ (black dashed line).

Having ascertained that our implementation of the IO dynamics correctly
samples a static quantity, we proceed by investigating the time behavior of
the correlation function of one subsystem. Our results are presented in Fig.
3. In the range of parameters explored, with $\gamma$ up to $2.4$, the IO
dynamics is accelerated by a non negligible amount compared to the standard
one. For $\gamma\approx 2.4$, the highest coupling achievable while preserving
the convergence of the integration scheme, the convergence to $0$ of the
correlation function is anticipated by almost an order of magnitude. For
$T<T_{d}$, the system freezes even at $\gamma\neq 0$, and the position of the
infinite time plateau is left unchanged.

Figure 3: Correlation function of system $1$, $C_{11}(t)$ for $\gamma=2.4$
(solid lines) and $\gamma=0$ (dashed lines) for various temperatures above and
below $T_{d}\approx 0.6124$. The dotted curves represent the cross
correlations $C_{12}(t)$ for $\gamma=2.4$. The shift between the dashed and
the solid line of a given color, which quantifies the acceleration of the
dynamics, somewhat increases as the temperature is decreased.

We quantify the speedup obtained by means of the IO dynamics by measuring the
relaxation time $\tau_{\alpha}(\gamma)$, defined here as the time necessary
for $C_{11}$ to decay to the value $1/e$ in the ergodic phase. Our results are
shown in Fig. 4. The ratio between $\tau_{\alpha}(\gamma)$ and the relaxation
time in the standard dynamics, $\tau_{\alpha}(0)$ appears to scale
approximatily as $(1+\gamma^{2})^{-1}$. This scaling the same as the one found
in Sec. III.1 when heuristically discussing the Kramers’ problem. We will
provide the reader with further interpretation for this phenomenon in the next
section.

Figure 4: Ratio between the relaxation time of $C_{11}$ for the IO dynamics,
$\tau_{\alpha}(\gamma)$, and the equilibrium dynamics $\tau_{\alpha}(0)$ for
temperatures close to $T_{d}$.

The departure of the IO dynamics from the equilibrium dynamics is accompanied
by an early growth in the cross correlations of the system. This can be seen
in Fig. 3, where we plot $C_{12}$ in dotted lines, and is presented in more
detail in Fig. 5, where a close up of the cross correlations for various
values of $\gamma$ is shown. The dynamics of these quantities is entirely
located in the region of $\beta$-relaxation of $C_{11}$. They exhibit a bump
at short times, reaching a point of maximum modulus and then decaying to zero.
This maximum grows linearly with $\gamma$ for $\gamma\lesssim 1$. The location
of the maximum moves left to shorter times as $\gamma$ grows, as expected from
Eq. (7).

Figure 5: Cross correlation $C_{12}$ (dotted lines) and $C_{21}$ (dashed
lines)as a function of time for $T=0.613$ at different values of $\gamma$.

Guided by the analysis of the harmonic oscillator in Sec. IV.2.1, we
investigate whether Eq. (38) holds for the $p$-spin, by studying the
relationship between each entry of the correlation matrix and the right hand
side of Eq. (39) and integrated response matrices. An example for a particular
entry and particular values of $T$ and $\gamma$ is shown in Fig. 6. There the
equality of Eq. (39) previously found for the harmonic oscillator is clearly
verified. Similar results are obtained for other entries and other values of
coupling and temperature. This is a consequence of the effective linear
dynamics of $\bm{\sigma}$ in Eq. (32).

Figure 6: Parametric plot of the left and right hand side of the first entry
of Eq. (39) from the integration of Eqs. (33, 34, 35) at $T=0.613$ and
$\gamma=2.4$.

The numerical investigation performed in this section raises a series of
questions to be answered. Can we show that the aFDT actually holds in the
$p$-spin IO dynamics? Can this result explain the absence of a shift in the
dynamical temperature $T_{d}$? Is it possible to justify the scaling of the
relaxation time of $C_{11}$ with $\gamma$? In the following section we will
build a theoretical understanding of our numerical findings.

####  IV.2.3 Interpretation of the aFDT

We establish the validity of the aFDT relation for the IO dynamics of the
$p-$spin. We can indeed show that if Eqs. (34, 35) and the aFDT relation Eq.
(38) are true then the remaining equation Eq. (33) must also hold. The details
of the proof are reported in App. D.

By using the aFDT relation we can reduce the sets of Eqs. (33, 34, 35) to a
single equation for the correlation matrix $\mathbf{C}(t)$:

$\begin{split}\partial_{t}\mathbf{C}(t)&=-T\left(\mathds{1}+\bm{\Gamma}\right)\cdot\mathbf{C}(t)\\\
&-\frac{1}{T}\left(\mathds{1}+\bm{\Gamma}\right)\int_{0}^{t}d\tau\mathbf{C}_{0}(p,t-\tau)\cdot\partial_{\tau}\mathbf{C}(\tau)\end{split}$
(41)

with $C_{0}(p,t)_{ij}\equiv\frac{p}{2}C_{ij}^{p-1}$. Apart from the matrix
prefactor $\left(\mathds{1}+\bm{\Gamma}\right)$ the structure of this equation
is the same as the one obtained in the standard equilibrium dynamics of the
$p$-spin [34, 36, 37]. We therefore study its long time behavior by assuming
that $\lim_{t\rightarrow+\infty}\mathbf{C}(t)=\mathbf{Q}$, with $\mathbf{Q}$ a
matrix of long time plateaus generalizing the nonergodic parameter. We
introduce the Laplace transform of a matrix $\mathbf{A}(t)$ as
$\mathbf{\tilde{A}}(z)\equiv\int_{0}^{+\infty}\differential
t\mathbf{A}(t)\text{e}^{-zt}$ and we observe that
$\mathbf{\tilde{C}}(z)\approx\frac{\mathbf{Q}}{z}$ and
$\mathbf{\tilde{C}}_{0}(p,z)\approx\frac{\mathbf{Q}_{0}(p)}{z}$ as
$z\rightarrow 0$, with $Q_{0}(p)_{ij}\equiv Q_{ij}^{p-1}$. Applying the
Laplace transform to both sides of Eq. (41) and keeping only the terms that
diverge as $z\to 0$ yields

$T\mathbf{Q}-\frac{1}{T}\mathbf{Q}_{0}(p)\cdot\left(\mathds{1}-\mathbf{Q}\right)=0$
(42)

We now see that the dependence on $\mathbf{\Gamma}$ has simply disappeared.
The form of the plateaus matrix can therefore not depend on $\gamma$. In
particular, it must match the solution obtained at $\gamma=0$. We can
therefore write $\mathbf{Q}=Q\mathds{1}$ and
$\mathbf{Q}_{0}(p)=Q^{p-1}\mathds{1}$. Eq. (42) then becomes the standard
equation for the long time plateau of the $p$-spin equilibrium dynamics

$\frac{2T^{2}}{p}=Q^{p-2}(1-Q)$ (43)

In analogy with the Kramers problem and the harmonic oscillator of Sec. III we
see that an increased mobility $\left(\mathds{1}+\bm{\Gamma}\right)$ arises in
the IO dynamics of the $p$-spin. This rescaling of the mobility can always be
absorbed into a rescaling of the time units. Below $T_{d}$, where the plateau
extends to infinite times, the dynamical ergodicity breaking cannot be avoided
while above $T_{d}$ an acceleration is achieved.

We support the above picture by carrying out a perturbation expansion to
second order in $\gamma$. In particular we will recover the scaling of the
relaxation time $\tau_{\alpha}$ shown in Fig. 4 starting from the dynamical
equation for the correlation matrix Eq. (41). We multiply both sides by
$\mathds{1}+\bm{\Gamma^{T}}$ and we expand $\mathbf{C}$ in powers of $\gamma$,
$\mathbf{C}\approx\delta_{0}\mathbf{C}+\gamma\delta_{1}\mathbf{C}+\gamma^{2}\delta_{2}\mathbf{C}$,
thereby obtaining:

$\begin{split}\partial_{t}\mathbf{C}(t)+\bm{\Gamma^{T}}\cdot\partial_{t}\delta_{1}\mathbf{C}(t)=&-(1+\gamma^{2})\Biggl{(}\mathbf{C}(t)\\\
&+\frac{1}{T}\int_{0}^{t}d\tau\mathbf{C}_{0}(p,t-\tau)\cdot\partial_{\tau}\mathbf{C}(\tau)\Biggr{)}\end{split}$
(44)

The second term on the left hand side contains $\delta_{1}\mathbf{C}$, the
first order perturbation of $\mathbf{C}$. It is a skew symmetric matrix and
therefore it contributes to the expansion of the cross correlations shown in
Fig. 5. Motivated by the numerical observation of Sec. IV.2.2 we assume that
this term essentially vanishes at times larger than the time scale at which
the $\beta$ relaxation takes place. Moreover, because of the skew symmetry of
the first order perturbation, the diagonal entries of Eq. (44) are not
affected by the short time dynamics of $\delta_{1}\mathbf{C}$. We therefore
obtain, to second order in $\gamma$:

$\begin{split}\partial_{t}C_{aa}(t)=&-(1+\gamma^{2})\Biggl{(}TC_{aa}(t)\\\
&+\frac{1}{T}\int_{0}^{t}d\tau
C_{aa,0}(p,t-\tau)\partial_{\tau}C_{aa}(\tau)\Biggr{)}\end{split}$ (45)

for $a=1,2$. This is the same equation as that of the equilibrium dynamics, up
to a rescaling of time $t\rightarrow t/(1+\gamma^{2})$, in agreement with the
findings of Fig. 4.

##  V Prospects

The nonequilibrium dynamics investigated here accelerates the relaxation of
the $p$-spin model without shifting the temperature at which a dynamical
ergodicity breaking occurs. A quantitative estimate of the speed up of the
dynamics is that an overall factor $1+\gamma^{2}$ in time is gained, for not
too large values of $\gamma$. We have come up with an explanation for this,
building on the analysis of the Kramers’ time for the IO dynamics. The
dynamics effectively yields an enhanced mobility, which accelerates the
dynamical evolution without affecting the effective height of the energetic
barrier. The presence of enhanced mobility explains both the preservation of
$T_{d}$ and the acceleration in the ergodic phase. We have provided evidence
for the scaling form of the acceleration with $\gamma$. At finite number $N$
of spins, it would interesting to probe the interplay of $\gamma$ with the
system size, as in principle, a sufficiently large $\gamma$ could overcome
energy barriers. To adjust $\gamma$ appropriately, one could rely on the
analysis of Stariolo and Cugliandolo [46, 47], where the relevant time-scales
are carefully sorted out. A related interesting question is how the instanton
trajectories are modified by a nonzero $\gamma$ in a finite size mean-field
model, expanding on the analysis of the Kramers problem made in the present
article. A relevant comparison could be provided by the study of Ros et al.
[48] where the instanton trajectory in the spherical $p$-spin model under
equilibrium overdamped Langevin dynamics is characterized.

While the numerical scheme we resort to prohibits the investigation of values
of $\gamma$ beyond $O(1)$, it is unclear at which value of $\gamma$ the
saturation of the speed up will occur. However, we speculate that such a
saturation is difficult to avoid. An interpretation of this effect that we
predict can be connected to the simple harmonic oscillator. Indeed, as
$\gamma$ is increased, one system is injected with a larger amount of energy
by the other one, but most of the time is spent swirling around in a circle a
finite distance from the bottom of the well, and this is exactly what would
happen in a harmonic oscillator, for which no acceleration is ever achieved.
We believe that coupling $n$ systems together would also lead to a marginal
gain as a similar saturation effect in terms of
$\gamma_{\text{eff}}\sim\sqrt{n}\gamma$ would occur for too large $n$. This,
in turn, raises the question of how to fine tune the optimal value of
$\gamma$. As our results suggest, it is possible that the best acceleration is
achieved by a nontrivial time-dependent protocol $\gamma(t)$ (in which
$\gamma(\to+\infty)\to\gamma_{\infty}$, the latter being possibly zero). This
question is touched upon in [16].

Beyond these formal but interesting questions, what we believe to be of true
physical interest is the application of such dynamics to molecular glass
formers. Unlike the $p$-spin, realistic glasses display only finite energy
barriers and the effect of $\gamma$ should be felt at all temperatures. How
this approach fares as compared to other ones such as SWAP [8, 9, 10] is, in
our opinion, and based on preliminary results [49], the most interesting open
question.

###### Acknowledgements.

We acknowledge very insightful exchanges with A. Altieri, T. Arnoulx de Pirey,
L. Berthier, L.F. Cugliandolo, J. Tailleur. We acknowledge the financial
support of the ANR grant THEMA.

##  Appendix A Instanton trajectory close to the saddle point

We study the instanton trajectory in the vicinity of the saddle point
$\mathbf{x}_{s}$. The first line of Eq. (21) reads

$\dot{\mathbf{x}}=\mathbf{A}\cdot\mathbf{x}$ (46)

with

$\mathbf{A}\equiv\begin{bmatrix}-k_{M}&\gamma k_{m}\\\ \gamma
k_{M}&k_{m}\end{bmatrix}$ (47)

This linear systems has one unstable and one stable direction along which the
dynamics takes the form $C_{\pm}e^{\lambda_{\pm}t}\mathbf{v}_{\pm}$,with
$\lambda_{\pm}$, $\mathbf{v}_{\pm}$ the eigenvalues and eigenvectors of
$\mathbf{A}$ and $C_{\pm}$ determined by the initial condition. Both the
eigenvalues and the eigenvectors depend on $\gamma$. The $t\rightarrow+\infty$
limit of the instanton dynamics in the first line of Eq. (46) corresponds to
the dynamics of system along its stable direction $\mathbf{v}_{-}$. By
expanding the latter around $\gamma=0$ we obtain

$e^{\lambda_{-}t}\mathbf{v}_{-}=\begin{bmatrix}-\text{e}^{-k_{M}t}+\gamma^{2}\text{e}^{-k_{M}t}\frac{k_{M}^{2}+2k_{m}t(k_{m}+k_{M})}{2(k_{m}+k_{M})^{2}}+{O}(\gamma^{4})\\\
\gamma\text{e}^{-k_{M}t}\frac{k_{M}}{k_{m}+k_{M}}+{O}(\gamma^{3})\end{bmatrix}$
(48)

which justifies the form of the perturbation expansion adopted in Eq. (22).

##  Appendix B Dynamical equations for the $p$-spin evolving under the IO
dynamics

We consider two copies of the spherical $p$-spin model evolving according to
the Ichiki-Ohzeki (IO) dynamics, for which the equations of motion read

$\displaystyle\partial_{t}\sigma_{i}^{(1)}=-\frac{\partial{\mathscr{H}}}{\partial\sigma_{i}^{(1)}}+\gamma\frac{\partial{\mathscr{H}}}{\partial\sigma_{i}^{(2)}}-\mu^{(1)}(t)\sigma_{i}^{(1)}+\gamma\mu^{(2)}(t)\sigma_{i}^{(2)}+\sqrt{2T}\xi_{i}^{(1)}$
(49)
$\displaystyle\partial_{t}\sigma_{i}^{(2)}=-\frac{\partial{\mathscr{H}}}{\partial\sigma_{i}^{(2)}}-\gamma\frac{\partial{\mathscr{H}}}{\partial\sigma_{i}^{(1)}}-\mu^{(2)}(t)\sigma_{i}^{(2)}-\gamma\mu^{(1)}(t)\sigma_{i}^{(1)}+\sqrt{2T}\xi_{i}^{(2)}$
(50)

where $\mu^{(a)}(t)$ are elastic constants that enforce the spherical
constraints
$\sum_{i=1}^{N}\left\langle\left(\sigma_{i}^{(a)}\right)^{2}(t)\right\rangle=N$
for $a=1,2$, $\xi_{i}^{(a)}$ are independent Gaussian white noises
$\left\langle\xi_{i}^{(a)}(t)\xi_{j}^{(b)}(t^{\prime})\right\rangle=\delta_{ab}\delta_{ij}\delta(t-t^{\prime})$
and the total Hamiltonian of the system $\mathscr{H}$ reads

$\mathscr{H}\left[\\{\sigma_{i}\\}^{(a)}_{i=1,\ldots,N,a=1,2}\right]=-\sum_{i_{1}<\ldots<i_{p}}J_{i_{1}\ldots
i_{p}}\sum_{a=1,2}\sigma^{(a)}_{i_{1}}\ldots\sigma^{(a)}_{i_{p}}$ (51)

Note that the couplings $J_{i_{1}\ldots i_{p}}$ have the same realization for
both systems, and they are independently distributed according to a Gaussian
distribution of mean $0$ and variance $\frac{p!}{2N^{p-1}}$.

We suppose that the systems start from an initial equilibrium condition at a
temperature $T\equiv\beta^{-1}$, which is kept constant throughout the
dynamical evolution. Our goal is to derive effective dynamical equations of
motion for the spins in the two system by averaging over the disorder. To
achieve this, we use standard path integral techniques as carried out in [50]
and further detailed in [36] or [37]. The starting point is the partition
function for the dynamics of the composite system in Eq. (49) at temperature
$T$. The initial condition for the two subsystems is given by two independent
realization of the Boltzmann distribution at temperature $T$. We have
therefore

$\displaystyle 1\equiv Z$
$\displaystyle=\int\mathcal{D}\bm{\sigma}_{0}\frac{1}{Z_{\text{eq}}}e^{-\beta\mathscr{H}[\sigma_{0}^{(1)},\sigma_{0}^{(2)}]}$
(52)
$\displaystyle\times\int\mathcal{D}\bm{\sigma}\mathcal{D}\bm{\hat{\sigma}}\exp[-\sum_{i=1}^{N}\int_{0}^{t}\differential\tau\hat{\sigma}_{i}^{(1)}\left(\partial_{\tau}\sigma^{(1)}_{i}+\frac{\partial{\mathscr{H}}}{\partial\sigma_{i}^{(1)}}-\gamma\frac{\partial{\mathscr{H}}}{\partial\sigma_{i}^{(2)}}+\mu^{(1)}\sigma_{i}^{(1)}-\gamma\mu^{(2)}\sigma_{i}^{(2)}\right)$
$\displaystyle+\hat{\sigma}_{i}^{(2)}\left(\partial_{\tau}\sigma^{(2)}_{i}+\frac{\partial{\mathscr{H}}}{\partial\sigma_{i}^{(2)}}+\gamma\frac{\partial{\mathscr{H}}}{\partial\sigma_{i}^{(1)}}+\mu^{(2)}\sigma_{i}^{(2)}+\gamma\mu^{(1)}\sigma_{i}^{(2)}\right)-T\left(\hat{\sigma}_{i}^{(1)}\right)^{2}-T\left(\hat{\sigma}_{i}^{(2)}\right)^{2}\Biggr{]}$

with $Z_{\text{eq}}$ the static partition function of the system,
$Z_{\text{eq}}\equiv\int\mathcal{D}\bm{\sigma}_{0}e^{-\beta\mathscr{H}[\bm{\sigma}_{0}]}$.
We are using a vector notation for the spin variables,
$\bm{\sigma}(t)=\left(\sigma^{(1)}(t),\sigma^{(2)}(t)\right)^{T}$, the
auxiliary field
$\bm{\hat{\sigma}}(t)=\left(\hat{\sigma}^{(1)}(t),\hat{\sigma}^{(2)}(t)\right)^{T}$
and the initial condition
$\bm{\sigma}_{0}(t)=\left(\sigma^{(1)}(0),\sigma^{(2)}(0)\right)^{T}$ As noted
in [51], since we are using thermalized initial conditions, we must resort to
the replica trick in order to average $Z$ over the disorder:

$\frac{1}{Z_{\text{eq}}}=\lim_{n\rightarrow 0}Z_{\text{eq}}^{n-1}$ (53)

and follow the dynamics of each replica. However, since the initial condition
belongs to the replica symmetric phase, the dynamics of the different replicas
are decoupled and it suffices to follow the behavior of any single one among
them. By using notation of the form
$\left(\sigma^{(a)}\cdot\sigma^{(b)}\right)(t,t^{\prime})\equiv\sum_{i=1}^{N}\sigma^{(a)}_{i}(t)\sigma^{(b)}_{i}(t^{\prime})$,
with $a$ and $b$ indices identifying the two systems, we can write in the
thermodynamic limit

$\displaystyle\overline{Z}$
$\displaystyle=\int\mathcal{D}\bm{\sigma}(0)\mathcal{D}\bm{\sigma}\mathcal{D}\bm{\hat{\sigma}}\exp\\{-\sum_{i=1}^{N}\Biggl{[}\int\differential\tau\hat{\sigma}_{i}^{(1)}\left(\partial_{\tau}\sigma_{i}^{(1)}+\mu^{(1)}\sigma_{i}^{(1)}-\gamma\mu^{(2)}\sigma_{i}^{(2)}\right)-T\left(\hat{\sigma}_{i}^{(1)}\right)^{2}$
(54)
$\displaystyle+\hat{\sigma}_{i}^{(2)}\left(\partial_{\tau}\sigma_{i}^{(2)}+\mu^{(2)}\sigma_{i}^{(2)}+\gamma\mu^{(1)}\sigma_{i}^{(1)}\right)-T\left(\hat{\sigma}_{i}^{(2)}\right)^{2}\Biggr{]}+\frac{1}{4N^{p-1}}\int\differential
t\int\differential
t^{\prime}I_{1}(t,t^{\prime})+\frac{1}{4TN^{p-1}}\int\differential
tI_{2}(t)\Biggr{\\}}$

With

$\displaystyle I_{1}(t,t^{\prime})$ $\displaystyle\equiv
p\left(\hat{\sigma}^{(1)}\cdot\hat{\sigma}^{(1)}\right)\left(\sigma^{(1)}\cdot\sigma^{(1)}\right)^{p-1}+p(p-1)\left(i\hat{\sigma}^{(1)}\cdot\sigma^{(1)}\right)\left(\sigma^{(1)}\cdot\hat{\sigma}^{(1)}\right)\left(\sigma^{(1)}\cdot\sigma^{(1)}\right)^{p-2}$
(55)
$\displaystyle-\gamma\left[p\left(\hat{\sigma}^{(1)}\cdot\hat{\sigma}^{(1)}\right)\left(\sigma^{(1)}\cdot\sigma^{(2)}\right)^{p-1}+p(p-1)\left(\hat{\sigma}^{(1)}\cdot\sigma^{(2)}\right)\left(\sigma^{(1)}\cdot\hat{\sigma}^{(1)}\right)\left(\sigma^{(1)}\cdot\sigma^{(2)}\right)^{p-2}\right]$
$\displaystyle+p\left(\hat{\sigma}^{(1)}\cdot\hat{\sigma}^{(2)}\right)\left(\sigma^{(1)}\cdot\sigma^{(2)}\right)^{p-1}+p(p-1)\left(\hat{\sigma}^{(1)}\cdot\sigma^{(2)}\right)\left(\sigma^{(1)}\cdot\hat{\sigma}^{(2)}\right)\left(\sigma^{(1)}\cdot\sigma^{(2)}\right)^{p-2}$
$\displaystyle+\gamma\left[p\left(\hat{\sigma}^{(1)}\cdot\hat{\sigma}^{(2)}\right)\left(\sigma^{(1)}\cdot\sigma^{(1)}\right)^{p-1}+p(p-1)\left(\hat{\sigma}^{(1)}\cdot\sigma^{(1)}\right)\left(\sigma^{(1)}\cdot\hat{\sigma}^{(2)}\right)\left(\sigma^{(1)}\cdot\sigma^{(1)}\right)^{p-2}\right]$
$\displaystyle-\gamma\left[p\left(\hat{\sigma}^{(1)}\cdot\hat{\sigma}^{(1)}\right)\left(\sigma^{(2)}\cdot\sigma^{(1)}\right)^{p-1}+p(p-1)\left(\hat{\sigma}^{(1)}\cdot\sigma^{(1)}\right)\left(\sigma^{(2)}\cdot\hat{\sigma}^{(1)}\right)\left(\sigma^{(2)}\cdot\sigma^{(1)}\right)^{p-2}\right]$
$\displaystyle+\gamma^{2}\left[p\left(\hat{\sigma}^{(1)}\cdot\hat{\sigma}^{(1)}\right)\left(\sigma^{(2)}\cdot\sigma^{(2)}\right)^{p-1}+p(p-1)\left(\hat{\sigma}^{(1)}\cdot\sigma^{(2)}\right)\left(\sigma^{(2)}\cdot\hat{\sigma}^{(1)}\right)\left(\sigma^{(2)}\cdot\sigma^{(2)}\right)^{p-2}\right]$
$\displaystyle-\gamma\left[p\left(\hat{\sigma}^{(1)}\cdot\hat{\sigma}^{(2)}\right)\left(\sigma^{(2)}\cdot\sigma^{(2)}\right)^{p-1}+p(p-1)\left(\hat{\sigma}^{(1)}\cdot\sigma^{(2)}\right)\left(\sigma^{(2)}\cdot\hat{\sigma}^{(2)}\right)\left(\sigma^{(2)}\cdot\sigma^{(2)}\right)^{p-2}\right]$
$\displaystyle-\gamma^{2}\left[p\left(\hat{\sigma}^{(1)}\cdot\hat{\sigma}^{(2)}\right)\left(\sigma^{(2)}\cdot\sigma^{(1)}\right)^{p-1}+p(p-1)\left(\hat{\sigma}^{(1)}\cdot\sigma^{(1)}\right)\left(\sigma^{(2)}\cdot\hat{\sigma}^{(2)}\right)\left(\sigma^{(2)}\cdot\sigma^{(1)}\right)^{p-2}\right]$
$\displaystyle+\left(1\leftrightarrow 2,\gamma\leftrightarrow-\gamma\right)$

and

$\displaystyle\frac{1}{2}I_{2}(t,0)$ $\displaystyle\equiv
p\left(\hat{\sigma}^{(1)}\cdot\sigma_{0}^{(1)}\right)\left(\sigma^{(1)}\cdot\sigma_{0}^{(1)}\right)^{p-1}-\gamma
p\left(\hat{\sigma}^{(1)}\cdot\sigma^{(1)}_{0}\right)\left(\sigma^{(2)}\cdot\sigma_{0}^{(1)}\right)^{p-1}$
(56)
$\displaystyle+p\left(\hat{\sigma}^{(1)}\cdot\sigma_{0}^{(2)}\right)\left(\sigma^{(1)}\cdot\sigma_{0}^{(2)}\right)^{p-1}-\gamma
p\left(\hat{\sigma}^{(1)}\cdot\sigma^{(2)}_{0}\right)\left(\sigma^{(2)}\cdot\sigma_{0}^{(2)}\right)^{p-1}$
$\displaystyle+\left(1\leftrightarrow 2,\gamma\leftrightarrow-\gamma\right)$

We now exploit the thermodynamic limit $N\rightarrow\infty$ through a saddle
point evaluation of the path integral. In order to do so, we introduce a set
of dynamical overlaps in the partition function:

$\displaystyle\overline{Z}$
$\displaystyle=\int\mathcal{D}\mathbf{Q}\Pi_{a,b=1,2}\delta\left(NQ_{1}^{(ab)}(t,t^{\prime})-\hat{\sigma}^{(a)}\cdot\hat{\sigma}^{(b)}\right)\delta\left(NQ_{2}^{(ab)}(t,t^{\prime})-\sigma^{(a)}\cdot\sigma^{(b)}\right)$
(57)
$\displaystyle\times\delta\left(NQ_{3}^{(ab)}(t,t^{\prime})-\hat{\sigma}^{(a)}\cdot\sigma^{(b)}\right)\delta\left(NQ_{4}^{(ab)}(t,t^{\prime})-\sigma^{(a)}\cdot\hat{\sigma}^{(b)}\right)$
$\displaystyle\times\delta\left(NQ_{5}^{(ab)}(t)-\hat{\sigma}^{(a)}\cdot\sigma_{0}^{(b)}\right)\times\delta\left(NQ_{6}^{(ab)}(t)-\sigma^{(a)}\cdot\sigma_{0}^{(b)}\right)\times\left(\text{r.h.s.
Eq. \eqref{eq:Zavg1}}\right)$
$\displaystyle=\int\mathcal{D}\mathbf{P}\int\mathcal{D}\mathbf{Q}\Pi_{a,b=1,2}\exp\left[iN\int\differential
t\differential
t^{\prime}\left(P_{1}^{(ab)}(t,t^{\prime})Q_{1}^{(ab)}-P_{1}^{(ab)}(t,t^{\prime})\hat{\sigma}^{(a)}\cdot\hat{\sigma}^{(b)}\right)\right]$
$\displaystyle\times\exp\left[iN\int\differential t\differential
t^{\prime}\left(P_{2}^{(ab)}Q_{2}^{(ab)}-P_{2}^{(ab)}\sigma^{(a)}\cdot\sigma^{(b)}\right)\right]\times\exp\left[iN\int\differential
t\differential
t^{\prime}\left(P_{3}^{(ab)}(t,t^{\prime})Q_{3}^{(ab)}-P_{3}^{(ab)}\hat{\sigma}^{(a)}\cdot\sigma^{(b)}\right)\right]$
$\displaystyle\times\exp\left[iN\int\differential t\differential
t^{\prime}\left(P_{4}^{(ab)}Q_{4}^{(ab)}-P_{4}^{(ab)}\sigma^{(a)}\cdot\hat{\sigma}^{(b)}\right)\right]\times\exp\left[iN\int\differential
t\left(P_{5}^{(ab)}Q_{5}^{(ab)}-P_{5}^{(ab)}\hat{\sigma}^{(a)}\cdot\sigma_{0}^{(b)}\right)\right]$
$\displaystyle\times\exp\left[iN\int\differential
t\left(P_{6}^{(ab)}Q_{6}^{(ab)}-P_{6}^{(ab)}\sigma^{(a)}\cdot\sigma_{0}^{(b)}\right)\right]\times\left(\text{r.h.s
of Eq. \eqref{eq:Zavg1}}\right)$

The dynamical overlaps have a natural interpretation if one introduces the
correlation and the response matrices $\mathbf{C}(t,t^{\prime})$ and
$\mathbf{R}(t,t^{\prime})$:

$\displaystyle C_{ab}(t,t^{\prime})$
$\displaystyle\equiv\frac{1}{N}\sum_{i}\langle\sigma_{i}^{(a)}(t)\sigma_{i}^{(b)}(t^{\prime})\rangle$
(58) $\displaystyle R_{ab}(t,t^{\prime})$
$\displaystyle\equiv\frac{1}{N}\sum_{i}\langle\sigma_{i}^{(a)}(t)\hat{\sigma}_{i}^{(b)}(t^{\prime})\rangle$
(59)

We therefore see that in the thermodynamic limit the following self consistent
relations hold:

$\begin{cases}Q_{2}^{(ab)}(t,t^{\prime})&=C_{ab}(t,t^{\prime})\\\
Q_{3}^{(ab)}(t,t^{\prime})&=R_{ab}(t^{\prime},t)\\\
Q_{4}^{(ab)}(t,t^{\prime})&=R_{ab}(t,t^{\prime})\\\
Q_{5}^{(ab)}(t)&=C_{ab}(t,0)\\\ Q_{6}^{(ab)}(t)&=R_{ab}(t,0)\end{cases}$ (60)

Note that $Q_{3}^{(ab)}(t,t^{\prime})=Q_{4}^{(ba)}(t^{\prime},t)$. Moreover,
we claim that, as in the standard case [52], $Q_{1}^{(ab)}=0$. The saddle
point equations are obtained by differentiating the argument of the
exponentials in Eq. (54) with respect to the set of dynamical overlaps
$Q_{i}^{(ab)}$ and they read:

$\begin{cases}iP_{1}^{(11)}&=\frac{p}{4}Q_{2}^{(11)p-1}-\gamma\frac{p}{4}\left(Q_{2}^{(12)p-1}+Q_{2}^{(21)p-1}\right)+\gamma^{2}\frac{p}{4}Q_{2}^{(22)p-1}\\\
iP_{2}^{(11)}&=0\\\
iP_{3}^{(11)}&=\frac{p}{4}(p-1)Q_{4}^{(11)}Q_{2}^{(11)p-2}+\gamma\frac{p}{4}(p-1)Q_{4}^{(12)}Q_{2}^{(11)p-2}-\gamma\frac{p}{4}(p-1)Q_{4}^{(21)}Q_{2}^{(21)p-2}-\gamma^{2}\frac{p}{4}(p-1)Q_{4}^{(11)}Q_{2}^{(21)p-2}\\\
iP_{4}^{(11)}&=\frac{p}{4}(p-1)Q_{3}^{(11)}Q_{2}^{(11)p-2}+\gamma\frac{p}{4}(p-1)Q_{3}^{(21)}Q_{2}^{(11)p-2}-\gamma\frac{p}{4}(p-1)Q_{3}^{(12)}Q_{2}^{(12)p-2}-\gamma^{2}\frac{p}{4}(p-1)Q_{3}^{(11)}Q_{2}^{(12)p-2}\\\
iP_{5}^{(11)}&=\frac{p}{2T}Q_{5}^{(11)}-\gamma\frac{p}{2T}Q_{5}^{(21)}\\\
iP_{6}^{(11)}&=0\\\ \\\
iP_{1}^{(12)}&=\frac{p}{4}Q_{2}^{(12)p-1}-\gamma\frac{p}{4}\left(Q_{2}^{(22)p-1}-Q_{2}^{(11)p-1}\right)-\gamma^{2}\frac{p}{4}Q_{2}^{(22)p-1}\\\
iP_{2}^{(12)}&=0\\\
iP_{3}^{(12)}&=\frac{p}{4}(p-1)Q_{4}^{(12)}Q_{2}^{(12)p-2}-\gamma\frac{p}{4}(p-1)Q_{4}^{(11)}Q_{2}^{(12)p-2}-\gamma\frac{p}{4}(p-1)Q_{4}^{(22)}Q_{2}^{(22)p-2}+\gamma^{2}\frac{p}{4}(p-1)Q_{4}^{(21)}Q_{2}^{(22)p-2}\\\
iP_{4}^{(12)}&=\frac{p}{4}(p-1)Q_{3}^{(12)}Q_{2}^{(12)p-2}+\gamma\frac{p}{4}(p-1)Q_{3}^{(11)}Q_{2}^{(11)p-2}+\gamma\frac{p}{4}(p-1)Q_{3}^{(22)}Q_{2}^{(12)p-2}+\gamma^{2}\frac{p}{4}(p-1)Q_{3}^{(21)}Q_{2}^{(11)p-2}\\\
iP_{5}^{(12)}&=\frac{p}{2T}Q_{5}^{(12)}-\gamma\frac{p}{2T}Q_{5}^{(22)}\\\
iP_{6}^{(12)}&=0\\\ \\\
iP_{1}^{(21)}&=\frac{p}{4}Q_{2}^{(21)p-1}-\gamma\frac{p}{4}\left(Q_{2}^{(11)p-1}-Q_{2}^{(22)p-1}\right)-\gamma^{2}\frac{p}{4}Q_{2}^{(12)p-1}\\\
iP_{2}^{(21)}&=0\\\
iP_{3}^{(21)}&=\frac{p}{4}(p-1)Q_{4}^{(21)}Q_{2}^{(21)p-2}+\gamma\frac{p}{4}(p-1)Q_{4}^{(11)}Q_{2}^{(11)p-2}+\gamma\frac{p}{4}(p-1)Q_{4}^{(22)}Q_{2}^{(21)p-2}+\gamma^{2}\frac{p}{4}(p-1)Q_{4}^{(12)}Q_{2}^{(11)p-2}\\\
iP_{4}^{(21)}&=\frac{p}{4}(p-1)Q_{3}^{(21)}Q_{2}^{(21)p-2}-\gamma\frac{p}{4}(p-1)Q_{3}^{(11)}Q_{2}^{(21)p-2}-\gamma\frac{p}{4}(p-1)Q_{3}^{(22)}Q_{2}^{(12)p-2}+\gamma^{2}\frac{p}{4}(p-1)Q_{3}^{(12)}Q_{2}^{(22)p-2}\\\
iP_{5}^{(21)}&=\frac{p}{2T}Q_{5}^{(21)}+\gamma\frac{p}{2T}Q_{5}^{(11)}\\\
iP_{6}^{(21)}&=0\\\ \\\
iP_{1}^{(22)}&=\frac{p}{4}Q_{2}^{(22)p-1}+\gamma\frac{p}{4}\left(Q_{2}^{(12)p-1}+Q_{2}^{(21)p-1}\right)+\gamma^{2}\frac{p}{4}Q_{2}^{(11)p-1}\\\
iP_{2}^{(22)}&=0\\\
iP_{3}^{(22)}&=\frac{p}{4}(p-1)Q_{4}^{(22)}Q_{2}^{(22)p-2}-\gamma\frac{p}{4}(p-1)Q_{4}^{(21)}Q_{2}^{(22)p-2}+\gamma\frac{p}{4}(p-1)Q_{4}^{(12)}Q_{2}^{(12)p-2}-\gamma^{2}\frac{p}{4}(p-1)Q_{4}^{(22)}Q_{2}^{(12)p-2}\\\
iP_{4}^{(22)}&=\frac{p}{4}(p-1)Q_{3}^{(22)}Q_{2}^{(22)p-2}-\gamma\frac{p}{4}(p-1)Q_{3}^{(12)}Q_{2}^{(22)p-2}+\gamma\frac{p}{4}(p-1)Q_{3}^{(21)}Q_{2}^{(21)p-2}-\gamma^{2}\frac{p}{4}(p-1)Q_{3}^{(22)}Q_{2}^{(21)p-2}\\\
iP_{5}^{(22)}&=\frac{p}{2T}Q_{5}^{(22)}+\gamma\frac{p}{2T}Q_{5}^{(12)}\\\
iP_{6}^{(22)}&=0\\\ \end{cases}$ (61)

Note that $iP_{2}=iP_{6}=0$ because $Q_{1}=0$ and because of causality, which
implies that $Q_{3}^{(ab)}(t,t^{\prime})Q_{4}^{(ab)}(t,t^{\prime})=0$.
Substitution of the saddle point expressions of Eq. (61) in Eq. (57), together
with Eq. (60) yields the effective dynamical equations for the spins of the
system. Using a vector notation,
$\bm{\sigma}(t)=\begin{bmatrix}\sigma^{(1)}&\sigma^{(2)}\end{bmatrix}^{T}$, we
obtain

$\partial_{t}\bm{\sigma}(t)=-\left(\mathds{1}+\bm{\Gamma}\right)\cdot\bm{\mu}(t)\cdot\bm{\sigma}(t)+\int_{0}^{t}d\tau\mathbf{M_{R}}(t,\tau)\cdot\bm{\sigma}(\tau)+\frac{1}{T}\left(\mathds{1}+\bm{\Gamma}\right)\cdot\mathbf{C}_{0}(p,t)\cdot\bm{\sigma}(0)+\bm{\xi}(t)$
(62)

where $\mu_{ab}(t)\equiv\delta_{ab}\mu^{(a)}(t)$,
$\bm{\Gamma}\equiv\begin{bmatrix}0&-\gamma\\\ \gamma&0\end{bmatrix}$. The
correlations of the noise $\bm{\xi}(t)$ are

$\left\langle\bm{\xi}(t)\bm{\xi}(t^{\prime})\right\rangle=2T\delta(t-t^{\prime})\mathds{1}_{2}+\mathbf{D}(t,t^{\prime})$
(63)

while the expression for the memory kernel $\mathbf{M_{R}}(t,t^{\prime})$, the
noise correlations $\mathbf{D}(t,t^{\prime})$ and the coupling with the
initial conditions $\mathbf{C}_{0}(t,0)$ are the following:

$\displaystyle\mathbf{M_{R}}(t,\tau)$
$\displaystyle=\frac{p(p-1)}{2}\begin{bmatrix}C_{11}^{p-2}\left(R_{11}+\gamma
R_{12}\right)-\gamma C_{21}^{p-2}\left(R_{21}+\gamma
R_{22}\right)&&&C_{12}^{p-2}\left(R_{12}-\gamma R_{11}\right)-\gamma
C_{22}^{p-2}\left(R_{22}-\gamma R_{21}\right)\\\ \\\
C_{21}^{p-2}\left(R_{21}+\gamma R_{22}\right)+\gamma
C_{11}^{p-2}\left(R_{11}+\gamma
R_{12}\right)&&&C_{22}^{p-2}\left(R_{22}-\gamma R_{21}\right)+\gamma
C_{12}^{p-2}\left(R_{12}-\gamma R_{11}\right)\end{bmatrix}$ (64)
$\displaystyle\mathbf{D}(t,t^{\prime})$
$\displaystyle\equiv\left(\mathds{1}+\bm{\Gamma}\right)\cdot\mathbf{C}_{0}(p,t)\cdot\left(\mathds{1}+\bm{\Gamma^{T}}\right)$
(65) $\displaystyle C_{0}(p,t,0)_{ij}$
$\displaystyle=\frac{p}{2}C_{ij}^{p-1}(t,0)$ (66)

From Eq. (62) by using the definition (58) and the equalities

$\displaystyle R_{ab}(t,t^{\prime})$
$\displaystyle=\left\langle\frac{\delta\sigma^{(a)}(t)}{\delta\xi^{(b)}(t^{\prime})}\right\rangle$
(67) $\displaystyle\langle\xi^{(a)}(t)\sigma^{(b)}(t^{\prime})\rangle$
$\displaystyle=2TR_{ab}(t^{\prime},t)+\int
d\tau\left[\mathbf{D}(t,\tau)\cdot\mathbf{R^{T}}(t^{\prime},\tau)\right]_{ab}$
(68)

we obtain

$\displaystyle\partial_{t}\mathbf{R}(t,t^{\prime})$
$\displaystyle=-\left(\mathds{1}+\bm{\Gamma}\right)\cdot\bm{\mu}(t)\cdot\mathbf{R}(t,t^{\prime})+\int_{t^{\prime}}^{t}d\tau\mathbf{M_{R}}(t,\tau)\cdot\mathbf{R}(\tau,t^{\prime})+\delta(t-t^{\prime})\mathds{1}_{2}$
(69) $\displaystyle\partial_{t}\mathbf{C}(t,t^{\prime})$
$\displaystyle=-\left(\mathds{1}+\bm{\Gamma}\right)\cdot\bm{\mu}(t)\cdot\mathbf{C}(t,t^{\prime})+\int_{0}^{t}d\tau\mathbf{M_{R}}(t,\tau)\cdot\mathbf{C}(\tau,t^{\prime})+2T\mathbf{R}(t^{\prime},t)$
(70)
$\displaystyle+\int_{0}^{t^{\prime}}d\tau\mathbf{D}(t,\tau)\cdot\mathbf{R^{T}}(t^{\prime},\tau)+\frac{1}{T}\left(\mathds{1}+\bm{\Gamma}\right)\cdot\mathbf{C}_{0}(p,t)\cdot\mathbf{C}(0,t^{\prime})$

By imposing the spherical constraint
$\sum_{i=1}^{N}\left\langle\left(\sigma_{i}^{(a)}\right)^{2}(t)\right\rangle=N$
for $a=1,2$ one can find a self-consistent expression for the restoring forces
$\bm{\mu}(t)$. If we denote the diagonal part of a matrix by
$\text{Diag}(\mathbf{A})_{ab}\equiv\delta_{ab}A_{ab}$ the spherical constraint
condition reads $\text{Diag}(\mathbf{C}(t,t))=1$. By differentiating with
respect to $t$ we get

$\lim_{t^{\prime}\rightarrow
t}\text{Diag}\left(\partial_{t}\mathbf{C}(t,t^{\prime})+\partial_{t^{\prime}}\mathbf{C}(t,t^{\prime})\right)=\lim_{t^{\prime}\rightarrow
t}\text{Diag}\left(\partial_{t}\mathbf{C}(t,t^{\prime})+\partial_{t^{\prime}}\mathbf{C^{T}}(t^{\prime},t)\right)=0$
(71)

which leads to

$\displaystyle\begin{split}\bm{\mu}(t)&=\text{Diag}\Biggl{[}T\mathds{1}_{2}+\int_{0}^{t}\differential\tau\mathbf{M}(t,\tau)\cdot\mathbf{C^{T}}(t,\tau)+\int_{0}^{t}d\tau\mathbf{D}(t,\tau)\cdot\mathbf{R^{T}}(t,\tau)\\\
&+\frac{1}{T}\left(\mathds{1}+\bm{\Gamma}\right)\cdot\mathbf{C}_{0}(p,t)\cdot\mathbf{C^{T}}(t,0)\Biggr{]}\end{split}$
(72)

Since the dynamics of the system starts within the steady state, we assume
that time translation invariance holds:
$\mathbf{C}(t,t^{\prime})\equiv\mathbf{C}(t-t^{\prime})$ and
$\mathbf{R}(t,t^{\prime})\equiv\mathbf{R}(t-t^{\prime})$ The equations of
motion for the response and correlation matrices therefore read, when taking
$t^{\prime}=0$

$\displaystyle\partial_{t}\mathbf{R}(t)$
$\displaystyle=-\left(\mathds{1}+\bm{\Gamma}\right)\cdot\bm{\mu}(t)\cdot\mathbf{R}(t)+\int_{t^{\prime}}^{t}\differential\tau\mathbf{M_{R}}(t-\tau)\cdot\mathbf{R}(\tau)+\delta(t)\mathds{1}_{2}$
(73) $\displaystyle\partial_{t}\mathbf{C}(t)$
$\displaystyle=-\left(\mathds{1}+\bm{\Gamma}\right)\cdot\bm{\mu}(t)\cdot\mathbf{C}(t)+\int_{0}^{t}\differential\tau\mathbf{M_{R}}(t-\tau)\cdot\mathbf{C}(\tau)+\frac{1}{T}\left(\mathds{1}+\bm{\Gamma}\right)\cdot\mathbf{C}_{0}(p,t)$
(74)
$\displaystyle\begin{split}\bm{\mu}(t)&=\text{Diag}\Biggl{[}T\mathds{1}_{2}+\int_{0}^{t}\differential\tau\mathbf{M}(t-\tau)\cdot\mathbf{C^{T}}(t-\tau)+\int_{0}^{t}\differential\tau\mathbf{D}(t-\tau)\cdot\mathbf{R^{T}}(t-\tau)\\\
&+\frac{1}{T}\left(\mathds{1}+\bm{\Gamma}\right)\cdot\mathbf{C}_{0}(p,t)\cdot\mathbf{C^{T}}(t)\Biggr{]}\end{split}$
(75)

For numerical purposes it is more convenient to work with the integrated
response matrix $\mathbf{F}(t,t^{\prime})$ [43], defined by

$F_{ab}(t,t^{\prime})\equiv-\int_{t^{\prime}}^{t}d\tau R_{ab}(t,\tau)$ (76)

since it has a smoother relaxation than $\mathbf{R}(t,t^{\prime})$. The
equations of motion in the time translation invariant regime are

$\displaystyle\partial_{t}\mathbf{F}(t)$
$\displaystyle=-\mathds{1}-\left(\mathds{1}+\bm{\Gamma}\right)\cdot\bm{\mu}(t)\cdot\mathbf{F}(t)+\int_{t^{\prime}}^{t}d\tau\mathbf{M_{F}}(t-\tau)\cdot\mathbf{F}(\tau)$
(77) $\displaystyle\partial_{t}\mathbf{C}(t)$
$\displaystyle=-\left(\mathds{1}+\bm{\Gamma}\right)\cdot\bm{\mu}(t)\cdot\mathbf{C}(t)+\int_{0}^{t}d\tau\mathbf{M_{F}}(t-\tau)\cdot\mathbf{C}(\tau)+\frac{1}{T}\left(\mathds{1}+\bm{\Gamma}\right)\cdot\mathbf{C}_{0}(p,t)$
(78)
$\displaystyle\begin{split}\bm{\mu}(t)&=\text{Diag}\Biggr{[}T\mathds{1}_{2}+\int_{0}^{t}d\tau\mathbf{M_{F}}(t-\tau)\cdot\mathbf{C^{T}}(t-\tau)+\int_{0}^{t}d\tau\mathbf{D}(t-\tau)\cdot\partial_{\tau}\mathbf{F^{T}}(t-\tau)\\\
&+\frac{1}{T}\left(\mathds{1}+\bm{\Gamma}\right)\cdot\mathbf{C}_{0}(p,t)\cdot\mathbf{C^{T}}(t)\Biggr{]}\end{split}$
(79)

and the memory kernel $\mathbf{M_{F}}(t)$ that reads

$\begin{split}\mathbf{M}_{\mathbf{F}}(t-\tau)&\equiv-\frac{p(p-1)}{2}\times\\\
&\begin{bmatrix}C_{11}^{p-2}\left(\partial_{\tau}F_{11}+\gamma\partial_{\tau}F_{12}\right)-\gamma
C_{21}^{p-2}\left(\partial_{\tau}F_{21}+\gamma\partial_{\tau}F_{22}\right)&&C_{12}^{p-2}\left(\partial_{\tau}F_{12}-\gamma\partial_{\tau}F_{11}\right)-\gamma
C_{22}^{p-2}\left(\partial_{\tau}F_{22}-\gamma\partial_{\tau}F_{21}\right)\\\
\\\
C_{21}^{p-2}\left(\partial_{\tau}F_{21}+\gamma\partial_{\tau}F_{22}\right)+\gamma
C_{11}^{p-2}\left(\partial_{\tau}F_{11}+\gamma\partial_{\tau}F_{12}\right)&&C_{22}^{p-2}\left(\partial_{\tau}F_{22}-\gamma\partial_{\tau}F_{21}\right)+\gamma
C_{12}^{p-2}\left(\partial_{\tau}F_{12}-\gamma\partial_{\tau}F_{11}\right)\end{bmatrix}\end{split}$
(80)

the companion initial conditions are $\mathbf{C}(0)\equiv\mathds{1}$,
$\mathbf{F}(0)=\mathbf{0}$. Together with Eq. (62), these equations correspond
to Eqs. (32, 33, 34, 35) presented in the main text.

In the following section we provide the reader with some details on the
algorithm implemented to integrate the dynamical equations.

##  Appendix C Numerical Integration of the dynamical equations

We present here the integration scheme adopted to numerically solve the
dynamical equations of motion, Eqs. (77, 78, 79). The basic idea for the
algorithm was first given in [42] in the context of the mode coupling
equations. It was then applied to the aging of the $p-$spin model in [43] and
generalised in [44] for the case of non stationary systems. A pedagogical
description of the algorithm can be found in [45].

To integrate Eqs. (77, 78, 79), we take a grid of $M$ points with a
discretisation step $\Delta t$, so that any function $f(t)$ can be discretised
as $f_{k}\equiv f(t_{k})$ with $t_{k}\equiv k\Delta t$. Suppose that we know
the first $n$ value of the correlation and the integrated response along the
discretized line and that we want to know their value at the step $n+1$. By
resorting to finite difference approximations, we can obtain a set of non-
linear equation for the variable
$\mathbf{V}_{i}\equiv\begin{bmatrix}\mathbf{C}_{i}&\mathbf{F}_{i}\end{bmatrix}$
of the form

$\mathbf{V}_{n+1}=\mathbf{NL}\left(\mathbf{V}_{0},\mathbf{V}_{1},...,\mathbf{V}_{n},\mathbf{V}_{n+1}\right)$
(81)

where $\mathbf{NL}$ is a generic nonlinear operator. This set of equation can
be solved iteratively starting from the initial condition
$\mathbf{V}_{n+1}=\mathbf{V}_{n}$. Once convergence is achieved, one can move
to the next point along the discretised line. This is the propagation step.

Once all the values of $\mathbf{C}$ and $\mathbf{F}$ have been obtained along
the grid, one performs a decimation and a rescaling to compress the
information obtained and reach longer times. Half of the points are discarded
by performing the substitution $\mathbf{V}_{2i}\rightarrow\mathbf{V}_{i}$ for
$i=0,...,M-1$ and the time step is rescaled by a factor two, $\Delta
t\rightarrow 2\Delta t$. The propagation step is then repeated starting from
the point $n+1=\frac{M}{2}$. Large times can be quickly achieved by iterating
these two steps.

The finite difference approximations that we adopt are the following. For the
left hand side of Eqs. (77, 78) the time derivative of a function at step
$k+1$, $\partial_{t}f_{k+1}$ can be approximated as

$\partial_{t}f_{k+1}=\frac{1}{2\Delta
t}\left(3f_{k+1}-4f_{k}+f_{k-1}\right)+O(\Delta t^{2})$ (82)

while the integrals on the right hand side have the general form

$\displaystyle I_{n+1}$
$\displaystyle\equiv\int_{0}^{t_{n+1}}\differential\tau
A(t-\tau)\partial_{\tau}B(t-\tau)C(\tau)$ (83) $\displaystyle J_{n+1}$
$\displaystyle\equiv\int_{0}^{t_{n+1}}\differential\tau
A(t-\tau)\partial_{\tau}B(t-\tau)C(\tau)$ (84)

and they can be approximated as

$\displaystyle I_{n+1}$
$\displaystyle=\sum_{j=0}^{n}\frac{1}{4}\left(A_{n+1-j}+A_{n-j}\right)\left(B_{n-j}-B_{n+1-j}\right)\left(C_{j}+C_{j+1}\right)+O(\Delta
t^{2})$ (85) $\displaystyle J_{n+1}$
$\displaystyle=\sum_{j=0}^{n}\frac{1}{4}\left(A_{n+1-j}+A_{n-j}\right)\left(B_{n-j}-B_{n+1-j}\right)\left(C_{n+1-j}+C_{n-j}\right)+O(\Delta
t^{2})$ (86)

For the results presented in the main text we have adopted a grid of $M=1024$
points and an initial timestep of $\Delta t=10^{-10}$.

##  Appendix D aFDT for the IO dynamics of the $p-$spin

Based on the results for the correlation and the response of an harmonic
oscillator we make the following ansatz concerning the form of the violation
of the FDT for the Twisted Dynamics of the $p-$spin:

$\frac{1}{T}\left[\mathbf{C}(t)-\mathbf{C}(0)\right]=\mathbf{F}(t)\cdot\left(\mathds{1}+\bm{\Gamma}\right)$
(88)

that is, an analogous form of the FDT relation holds if one looks at the
proper combination of response functions. We call this the aFDT relation. In
this section we will show that this relation is indeed satisfied by Eqs. (77,
78, 79). The strategy of the proof is the following: starting from Eqs. (78,
79) and the aFDT we will obtain Eq. (77). To do so, we first rewrite Eq. (88)
in a differential form:

$\frac{1}{T}\partial_{t}\mathbf{C}(t)=\partial_{t}\mathbf{F}(t)\cdot\left(\mathds{1}+\bm{\Gamma}\right)$
(89)

using Eq. (89) we can rewrite the memory kernel $\mathbf{M_{F}}$ in a simpler
form:

$\mathbf{M_{F}}(t-\tau)=\frac{1}{T}\left(\mathds{1}+\bm{\Gamma}\right)\cdot\partial_{\tau}\mathbf{C}_{0}(t-\tau)$
(90)

We now proceed to compute the left hand side of Eq. (89) using Eq. (74):

$\displaystyle\begin{split}\frac{1}{T}\mathbf{C}(t)&=-\frac{1}{T}\left(\mathds{1}+\bm{\Gamma}\right)\cdot\bm{\mu}(t)\cdot\mathbf{C}(t)+\frac{1}{T}\int_{0}^{t}d\tau\mathbf{M_{F}}(t-\tau)\cdot\mathbf{C}(\tau)+\frac{1}{T^{2}}\left(\mathds{1}+\bm{\Gamma}\right)\cdot\mathbf{C}_{0}(t)\\\
&=-\left(\mathds{1}+\bm{\Gamma}\right)\cdot\left(\frac{1}{T}\mathds{1}+\mathbf{F}(t)\cdot\left(\mathds{1}+\bm{\Gamma}\right)\right)+\frac{1}{T}\int_{0}^{t}d\tau\mathbf{M_{F}}(t-\tau)\cdot\mathbf{C}(\tau)+\frac{1}{T^{2}}\left(\mathds{1}+\bm{\Gamma}\right)\cdot\mathbf{C}_{0}(t)\end{split}$
(91)
$\displaystyle\begin{split}\bm{\mu}(t)&=\text{Diag}\Biggl{[}T\mathds{1}_{2}+\int_{0}^{t}d\tau\mathbf{M_{F}}(t-\tau)\cdot\mathbf{C^{T}}(t-\tau)+\int_{0}^{t}d\tau\mathbf{D}(t-\tau)\cdot\partial_{\tau}\mathbf{F^{T}}(t-\tau)\\\
&+\frac{1}{T}\left(\mathds{1}+\bm{\Gamma}\right)\cdot\mathbf{C}_{0}(p-1,t)\cdot\mathbf{C^{T}}(t)\Biggr{]}\end{split}$
(92)

to evaluate the time integral, we can use Eq. (90), integration by parts, the
aFDT relation and the boundary condition $\mathbf{C}(0)=\mathds{1}$,
$\mathbf{F}(0)=\mathbf{0}$ to obtain

$\displaystyle\frac{1}{T}\int_{0}^{t}d\tau\mathbf{M_{F}}(t-\tau)\cdot\mathbf{C}(\tau)$
$\displaystyle=\frac{1}{T^{2}}\left(\frac{p}{2}\mathds{1}-\mathbf{C}_{0}(t)\right)+\int_{0}^{t}d\tau\mathbf{M_{F}}(t-\tau)\cdot\mathbf{F}(\tau)\cdot\left(\mathds{1}+\bm{\Gamma}\right)$
(93)
$\displaystyle\int_{0}^{t}d\tau\mathbf{D}(t-\tau)\cdot\partial_{\tau}\mathbf{F^{T}}(t-\tau)$
$\displaystyle=\frac{1}{2T}\left(\mathds{1}+\bm{\Gamma}\right)\cdot\int_{0}^{t}d\tau\mathbf{C}_{0}\cdot\partial_{\tau}\mathbf{C^{T}}(\tau)$
(94)
$\displaystyle\int_{0}^{t}d\tau\mathbf{M_{F}}(t-\tau)\cdot\mathbf{C^{T}}(t-\tau)$
$\displaystyle=\frac{1}{T}\left[\frac{p}{2}\mathds{1}-\mathbf{C}_{0}(t)\cdot\mathbf{C^{T}}(t)-\frac{1}{T}\int_{0}^{t}d\tau\mathbf{C}_{0}(t-\tau)\cdot\partial_{\tau}\mathbf{C^{T}}(\tau)\right]$
(95)

So that $\bm{\mu}$ reads after simplifying and taking the diagonal part:

$\bm{\mu}(t)=\left(T+\frac{p}{2T}\right)\mathds{1}\equiv\bm{\mu}_{\infty}$
(96)

putting everything back in Eq. (91) and simplifying we get that

$\frac{1}{T}\partial_{t}\mathbf{C}(t)=\partial_{t}\mathbf{F}(t)\cdot\left(\mathds{1}+\bm{\Gamma}\right)$
(97)

with $\partial_{t}\mathbf{F}(t)$ given by Eq. (77), and the proof is
completed.

## References

  * Kirkpatrick  _et al._ [1983] S. Kirkpatrick, C. D. Gelatt Jr, and M. P. Vecchi, Optimization by simulated annealing, science 220, 671 (1983).
  * Mitsutake  _et al._ [2001] A. Mitsutake, Y. Sugita, and Y. Okamoto, Generalized-ensemble algorithms for molecular simulations of biopolymers, Peptide Science: Original Research on Biomolecules 60, 96 (2001).
  * Welling and Teh [2011] M. Welling and Y. W. Teh, Bayesian learning via stochastic gradient langevin dynamics, in  _Proceedings of the 28th international conference on machine learning (ICML-11)_ (Citeseer, 2011) pp. 681–688.
  * Andrieu  _et al._ [2003] C. Andrieu, N. De Freitas, A. Doucet, and M. I. Jordan, An introduction to mcmc for machine learning, Machine learning 50, 5 (2003).
  * Nesterov [1983] Y. E. Nesterov, A method for solving the convex programming problem with convergence rate $o(1/k^{2})$, in  _Dokl. akad. nauk Sssr_ , Vol. 269 (1983) pp. 543–547.
  * Swendsen and Wang [1986] R. H. Swendsen and J.-S. Wang, Replica monte carlo simulation of spin-glasses, Physical review letters 57, 2607 (1986).
  * Marinari and Parisi [1992] E. Marinari and G. Parisi, Simulated tempering: a new monte carlo scheme, EPL (Europhysics Letters) 19, 451 (1992).
  * Grigera and Parisi [2001] T. S. Grigera and G. Parisi, Fast monte carlo algorithm for supercooled soft spheres, Physical Review E 63, 045102 (2001).
  * Berthier  _et al._ [2016] L. Berthier, D. Coslovich, A. Ninarello, and M. Ozawa, Equilibrium sampling of hard spheres up to the jamming density and beyond, Physical review letters 116, 238002 (2016).
  * Ninarello  _et al._ [2017] A. Ninarello, L. Berthier, and D. Coslovich, Models and algorithms for the next generation of glass transition studies, Physical Review X 7, 021039 (2017).
  * Suwa and Todo [2010] H. Suwa and S. Todo, Markov chain monte carlo method without detailed balance, Physical review letters 105, 120603 (2010).
  * Dress and Krauth [1995] C. Dress and W. Krauth, Cluster algorithm for hard spheres and related systems, Journal of Physics A: Mathematical and General 28, L597 (1995).
  * Kapfer and Krauth [2017] S. C. Kapfer and W. Krauth, Irreversible local markov chains with rapid convergence towards equilibrium, Physical review letters 119, 240603 (2017).
  * Hwang  _et al._ [1993] C.-R. Hwang, S.-Y. Hwang-Ma, and S.-J. Sheu, Accelerating gaussian diffusions, The Annals of Applied Probability , 897 (1993).
  * Hwang  _et al._ [2005] C.-R. Hwang, S.-Y. Hwang-Ma, and S.-J. Sheu, Accelerating diffusions, The Annals of Applied Probability 15, 1433 (2005).
  * Lelievre  _et al._ [2013] T. Lelievre, F. Nier, and G. A. Pavliotis, Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion, Journal of Statistical Physics 152, 237 (2013).
  * Ichiki and Ohzeki [2013] A. Ichiki and M. Ohzeki, Violation of detailed balance accelerates relaxation, Physical Review E 88, 020101 (2013).
  * Ohzeki and Ichiki [2015a] M. Ohzeki and A. Ichiki, Langevin dynamics neglecting detailed balance condition, Physical Review E 92, 012105 (2015a).
  * Ohzeki and Ichiki [2015b] M. Ohzeki and A. Ichiki, Mathematical understanding of detailed balance condition violation and its application to langevin dynamics, Journal of Physics: Conference Series 638, 012003 (2015b).
  * Rey-Bellet and Spiliopoulos [2015a] L. Rey-Bellet and K. Spiliopoulos, Irreversible langevin samplers and variance reduction: a large deviations approach, Nonlinearity 28, 2081 (2015a).
  * Rey-Bellet and Spiliopoulos [2015b] L. Rey-Bellet and K. Spiliopoulos, Variance reduction for irreversible langevin samplers and diffusion on graphs, Electronic Communications in Probability 20, 1 (2015b).
  * Duncan  _et al._ [2016] A. B. Duncan, T. Lelievre, and G. A. Pavliotis, Variance reduction using nonreversible langevin samplers, Journal of statistical physics 163, 457 (2016).
  * Chen  _et al._ [1999] F. Chen, L. Lovász, and I. Pak, Lifting markov chains to speed up mixing, in  _Proceedings of the thirty-first annual ACM symposium on Theory of computing_ (1999) pp. 275–281.
  * Diaconis  _et al._ [2000] P. Diaconis, S. Holmes, and R. M. Neal, Analysis of a nonreversible markov chain sampler, Annals of Applied Probability , 726 (2000).
  * Apers  _et al._ [2017] S. Apers, A. Sarlette, and F. Ticozzi, When does memory speed-up mixing?, in  _2017 IEEE 56th Annual Conference on Decision and Control (CDC)_ (Ieee, 2017) pp. 4940–4945.
  * Turitsyn  _et al._ [2011] K. S. Turitsyn, M. Chertkov, and M. Vucelja, Irreversible monte carlo algorithms for efficient sampling, Physica D: Nonlinear Phenomena 240, 410 (2011).
  * Vucelja [2016] M. Vucelja, Lifting—a nonreversible markov chain monte carlo algorithm, American Journal of Physics 84, 958 (2016).
  * Bouchet and Reygner [2016] F. Bouchet and J. Reygner, Generalisation of the eyring–kramers transition rate formula to irreversible diffusion processes, in  _Annales Henri Poincaré_ , Vol. 17 (Springer, 2016) pp. 3499–3532.
  * Darboux [1946] G. Darboux,  _Leçons sur la théorie générale des surfaces et les applications géométriques du calcul infinitésimal: 4e partie: Déformation infiniment petite et représentation sphérique_ (Gauthier-Villars, 1946).
  * Ichiki and Ohzeki [2021] A. Ichiki and M. Ohzeki, Stochastic dynamics without detailed balance condition connecting simple gradient method and hamiltonian monte carlo, arXiv preprint arXiv:2107.08633 (2021).
  * Kramers [1940] H. A. Kramers, Brownian motion in a field of force and the diffusion model of chemical reactions, Physica 7, 284 (1940).
  * Bray and McKane [1989] A. Bray and A. McKane, Instanton calculation of the escape rate for activation over a potential barrier driven by colored noise, Physical review letters 62, 493 (1989).
  * Crisanti and Sommers [1992] A. Crisanti and H.-J. Sommers, The spherical p-spin interaction spin glass model: the statics, Zeitschrift für Physik B Condensed Matter 87, 341 (1992).
  * Crisanti  _et al._ [1993] A. Crisanti, H. Horner, and H.-J. Sommers, The sphericalp-spin interaction spin-glass model, Zeitschrift für Physik B Condensed Matter 92, 257 (1993).
  * Barrat  _et al._ [1996] A. Barrat, R. Burioni, and M. Mézard, Dynamics within metastable states in a mean-field spin glass, Journal of Physics A: Mathematical and General 29, L81 (1996).
  * Barrat [1997] A. Barrat, The p-spin spherical spin glass model, arXiv preprint cond-mat/9701031 (1997).
  * Castellani and Cavagna [2005] T. Castellani and A. Cavagna, Spin-glass theory for pedestrians, Journal of Statistical Mechanics: Theory and Experiment 2005, P05012 (2005).
  * De Dominicis and Giardina [2006] C. De Dominicis and I. Giardina,  _Random fields and spin glasses: a field theory approach_ (Cambridge University Press, 2006).
  * Berthier and Kurchan [2013] L. Berthier and J. Kurchan, Non-equilibrium glass transitions in driven and active matter, Nature Physics 9, 310 (2013).
  * Cugliandolo  _et al._ [1997] L. F. Cugliandolo, J. Kurchan, P. Le Doussal, and L. Peliti, Glassy behaviour in disordered systems with nonrelaxational dynamics, Physical review letters 78, 350 (1997).
  * Berthier  _et al._ [2000] L. Berthier, J.-L. Barrat, and J. Kurchan, A two-time-scale, two-temperature scenario for nonlinear rheology, Physical Review E 61, 5464 (2000).
  * Fuchs  _et al._ [1991] M. Fuchs, W. Gotze, I. Hofacker, and A. Latz, Comments on the alpha-peak shapes for relaxation in supercooled liquids, Journal of Physics: Condensed Matter 3, 5047 (1991).
  * Kim and Latz [2001] B. Kim and A. Latz, The dynamics of the spherical p-spin model: From microscopic to asymptotic, EPL (Europhysics Letters) 53, 660 (2001).
  * Berthier  _et al._ [2007] L. Berthier, G. Biroli, J.-P. Bouchaud, W. Kob, K. Miyazaki, and D. R. Reichman, Spontaneous and induced dynamic correlations in glass formers. ii. model calculations and comparison to numerical simulations, The Journal of chemical physics 126, 184504 (2007).
  * Folena [2020] G. Folena,  _The mixed p-spin model: selecting, following and losing states_ , Ph.D. thesis, Université Paris-Saclay; Università degli studi La Sapienza (Rome) (2020).
  * Stariolo and Cugliandolo [2019] D. A. Stariolo and L. F. Cugliandolo, Activated dynamics of the ising p-spin disordered model with finite number of variables, EPL (Europhysics Letters) 127, 16002 (2019).
  * Stariolo and Cugliandolo [2020] D. A. Stariolo and L. F. Cugliandolo, Barriers, trapping times, and overlaps between local minima in the dynamics of the disordered ising p-spin model, Physical Review E 102, 022126 (2020).
  * Ros  _et al._ [2021] V. Ros, G. Biroli, and C. Cammarota, Dynamical instantons and activated processes in mean-field glass models, SciPost Physics 10, 002 (2021).
  * Ghimenti [2022] F. Ghimenti, In preparation. (2022).
  * Cugliandolo and Kurchan [1993] L. F. Cugliandolo and J. Kurchan, Analytical solution of the off-equilibrium dynamics of a long-range spin-glass model, Physical Review Letters 71, 173 (1993).
  * Houghton  _et al._ [1983] A. Houghton, S. Jain, and A. Young, Role of initial conditions in the mean-field theory of spin-glass dynamics, Physical Review B 28, 2630 (1983).
  * Sompolinsky and Zippelius [1982] H. Sompolinsky and A. Zippelius, Relaxational dynamics of the edwards-anderson model and the mean-field theory of spin-glasses, Physical Review B 25, 6860 (1982).