Federico Ghimenti

Frédéric Van Wijland

Accelerating, to some extent, the p-spin dynamics

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 β and α-relaxation stages. An interpretation in terms of trajectories in phase space and of an accidental fluctuation-dissipation theorem is provided.

I.

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 B ∼ e -βH for a system with energy H concentrates, at low temperature T = 1/β, 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 [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]. 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 [START_REF] Mitsutake | Generalizedensemble algorithms for molecular simulations of biopolymers[END_REF]. It is of course abundantly discussed in the fields of optimization [START_REF] Welling | Bayesian learning via stochastic gradient langevin dynamics[END_REF] and machine learning [START_REF] Andrieu | An introduction to mcmc for machine learning[END_REF] as well. For example, to accelerate the Stochastic Gradient Method, Nesterov [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate o(1/k 2 )[END_REF] 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 [START_REF] Swendsen | Replica monte carlo simulation of spin-glasses[END_REF], parallel tempering [START_REF] Marinari | Simulated tempering: a new monte carlo scheme[END_REF] or the SWAP Monte Carlo method [START_REF] Grigera | Fast monte carlo algorithm for supercooled soft spheres[END_REF][START_REF] Berthier | Equilibrium sampling of hard spheres up to the jamming density and beyond[END_REF][START_REF] Ninarello | Models and algorithms for the next generation of glass transition studies[END_REF].

When designing alternative dynamics, it is also possible not to respect the detailed balance condition, now a widely used procedure in Monte Carlo schemes [START_REF] Suwa | Markov chain monte carlo method without detailed balance[END_REF][START_REF] Dress | Cluster algorithm for hard spheres and related systems[END_REF], that leads to a faster convergence towards equilibrium [START_REF] Kapfer | Irreversible local markov chains with rapid convergence towards equilibrium[END_REF]. 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 [START_REF] Hwang | Accelerating gaussian diffusions[END_REF][START_REF] Hwang | Accelerating diffusions[END_REF][START_REF] Lelievre | Optimal nonreversible linear drift for the convergence to equilibrium of a diffusion[END_REF][START_REF] Ichiki | Violation of detailed balance accelerates relaxation[END_REF][START_REF] Ohzeki | Langevin dynamics neglecting detailed balance condition[END_REF][START_REF] Ohzeki | Mathematical understanding of detailed balance condition violation and its application to langevin dynamics[END_REF][START_REF] Rey-Bellet | Irreversible langevin samplers and variance reduction: a large deviations approach[END_REF][START_REF] Rey-Bellet | Variance reduction for irreversible langevin samplers and diffusion on graphs[END_REF][START_REF] Duncan | Variance reduction using nonreversible langevin samplers[END_REF]. The termi- * federico.ghimenti@u-paris.fr nology for a concrete implementation refers to a "lifting" procedure [START_REF] Chen | Lifting markov chains to speed up mixing[END_REF][START_REF] Diaconis | Analysis of a nonreversible markov chain sampler[END_REF][START_REF] Apers | When does memory speed-up mixing?[END_REF] , 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 [START_REF] Turitsyn | Irreversible monte carlo algorithms for efficient sampling[END_REF][START_REF] Vucelja | Lifting-a nonreversible markov chain monte carlo algorithm[END_REF].

In this work we shall build on the version proposed by Ichiki and Ohzeki [START_REF] Ohzeki | Langevin dynamics neglecting detailed balance condition[END_REF] (which we hereafter call the IO dynamics). They start from two independent and identical systems with a reversible dynamics whose equilibrium distribution P B is known. They then add nonconservative solenoidal forces that couple the two systems but chosen such that P B remains the stationary distribution of the resulting dissipative dynamics for each individual system. The strength γ 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 [START_REF] Bouchet | Generalisation of the eyringkramers transition rate formula to irreversible diffusion processes[END_REF] in which they determined a formal expression for the modified Kramers escape time. In addition to specializing the results of [START_REF] Bouchet | Generalisation of the eyringkramers transition rate formula to irreversible diffusion processes[END_REF] 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 (r) at temperature T . The equation of motion reads

ṙ = -∇ r V (r) + √ 2T η (1)
where η is a Gaussian white noise with unit correlations η(t)⊗η(t ) = 1δ(t-t ). The companion Fokker-Planck equation

∂ t P = -∇ r • j, j(r, t) = -T ∇ r P -∇ r V P (2)
has the Boltzmann distribution P B (r) ∝ e -βV (r) as its stationary solution. Because of the underlying detailed balance condition, this stationary solution P B is actually an equilibrium one, as it makes the probability current j vanish, and no entropy is produced in the equilibrium steady-state. Consider now two identical particles with positions r (1) and r (2) evolving according to the coupled Langevin dynamics with independent white Gaussian noises η (1) and η (2) , ṙ(1) = -∇ r (1) V (r (1) ) + γ∇ r (2) V (r (2) ) +

√

2T η (1) (3)

ṙ(2) = -∇ r (2) V (r (2) ) -γ∇ r (1) V (r (1) ) + √ 2T η (2) (4)
The parameter γ introduces a skew symmetric coupling between the two systems. Two independent overdamped Langevin dynamics for systems 1 and 2 are recovered when γ = 0. Note that when γ = 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 (r (1) , r (2) ) ≡ V (r (1) ) + V (r (2) ) we have

d dt U (r (1) , r (2) ) = - a=1,2 ∇ r (a) V (r (a) ) 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 (r (1) , r (2) , t) follows from Eq. [START_REF] Welling | Bayesian learning via stochastic gradient langevin dynamics[END_REF]. With respect to the expression of Eq. ( 2) picks up an additional contribution for each of the two subspaces, namely δj (1) = γ∇ r (1) V (r (2) )P (r (1) , r (2) , t) δj (2) = -γ∇ r (2) V (r (1) )P (r (1) , r (2) , t) [START_REF] Swendsen | Replica monte carlo simulation of spin-glasses[END_REF] This additional current arises from the presence of the nonconservative forces in Eq. ( 3), and it is divergencefree when P (r (1) , r (2) , t) = P B (r 1 )P B (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 Σ given by

Σ = 2βγ 2 (∇V ) 2 B = 2γ 2 ∆V B > 0 (7)
where . . . B refers to an average with respect to the Boltzmann distribution, thus confirming the nonequilibrium nature of the dynamics, and that P B is a nonequilibrium steady-state distribution. This calculation tells us about the typical time scale 1/ Σ 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).

A. A simple argument for the decrease of the relaxation time

A mathematical proof accounting for the decrease of the relaxation time for arbitrary γ was given in [START_REF] Hwang | Accelerating gaussian diffusions[END_REF][START_REF] Hwang | Accelerating diffusions[END_REF] and it was rediscovered in the case of discrete state Markov processes in [START_REF] Ichiki | Violation of detailed balance accelerates relaxation[END_REF]. However, for small γ, 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 [START_REF] 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[END_REF] and we rewrite the Fokker-Planck operator W associated to Eq. (3) in a rotated basis, W = P

1/2 B WP -1/2 B , which reads W = W(1) 0 + W(2) 0 + δ W (8)
where

W(a) 0 ≡ T ∇ 2 r (a) - 1 4T (∇ r (a) V (r (a) )) 2 -2T ∇ 2 r (a) V (r (a) ) (9)
is the symmetrized Fokker-Planck operator for particle a = 1, 2 in the absence of the coupling parameter γ and

δ W ≡ γ [(∇ r1 V (r 1 )) • ∇ r2 -(∇ r2 V (r 2 )) • ∇ r1 ] (10)
As expected, for γ = 0 this is a Hermitian operator (which reflects the detailed balance condition), but for γ = 0 a skew-Hermitian perturbation is introduced in Eq. [START_REF] Grigera | Fast monte carlo algorithm for supercooled soft spheres[END_REF]. We denote by |Ψ n the eigenstates of the unperturbed Hermitian operator W(1) 0 + W(2) 0 , and with E n the associated eigenvalue. We are interested in the change δE to the spectral gap ∆E = E 1 -E 0 of the operator -W between the ground state |Ψ 0 and the first excited state |Ψ 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

∆E = -γ Ψ 1 |δ W|Ψ 1 -γ 2 n>1 Ψ 1 |δ W|ψ n 2 E 1 -E n (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 γ, the eigenvalues of 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 [START_REF] Ohzeki | Langevin dynamics neglecting detailed balance condition[END_REF][START_REF] Ohzeki | Mathematical understanding of detailed balance condition violation and its application to langevin dynamics[END_REF][START_REF] Ichiki | Stochastic dynamics without detailed balance condition connecting simple gradient method and hamiltonian monte carlo[END_REF], 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 τ it takes for particle 1 to hop over a potential barrier ∆V , passing through the local maximum at x M , starting from the bottom of the well. As we know from Kramers' analysis [START_REF] Kramers | Brownian motion in a field of force and the diffusion model of chemical reactions[END_REF], for equilibrium dynamics, τ is given by the Arrhenius formula τ e β∆V . Consider now a nonzero γ 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 [START_REF] Bouchet | Generalisation of the eyringkramers transition rate formula to irreversible diffusion processes[END_REF].

A. Heuristic argument

Before proceeding with a detailed analysis, there is a simple way to understand how a nonzero γ 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

ẋ(2) (t) = -V (x (2) (t)) -γV (x (1) (t)) + √ 2T η (2) (t) -k(x (2) (t) -x m ) -γV (x (1) (t)) + √ 2T η (2) (t) (12)
where k m ≡ V (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) :

ẋ(1) (t) = - t 0 dsM R (t -s)V (x (1) (s)) + √ 2T η(t) (13)
where η(t) is now a colored Gaussian noise with correlations η(t)η(t 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)

) = M C (t -t ). The mobility ker- nel M R (t) = δ(t) + γ 2 k m e -kmt θ(t) and the correla- tion kernel M C (t) = 1 2 M R (t) + 1 2 M R (-
(1 + γ 2 )δ(t) and the mobility of the particle picks up a γ 2 contribution so that the Kramers' escape time is reduced accordingly: τ e β∆V /(1 + γ 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

ẋ = -(1 + Γ) • ∇U (x) + √ 2T η(t) (14)
with x ≡ x (1) x (2) T , U (x) ≡ V (x (1) ) + V (x (2) ) and

Γ ≡ 0 -γ γ 0 (15)
This two-dimensional system presents a local minimum at x m ≡ x m x m T and a saddle point at x s ≡

x M x m 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 x m to x s and from x s to the minimum on the other side of the saddle. We now proceed with a careful analysis of the instanton trajectory.

B. Instanton trajectory: no change in the exponential factor

The Kramers escape rate τ -1 can be expressed within a path-integral formulation,

τ -1 = D xDxe -1 T S[x,x] (16)
where the path-integral runs over paths connecting x m to x s , and where the Janssen-De Dominicis action functional reads

S [x, x] = +∞ -∞ dt [ ẋ • x -H (x, x)] (17)
The response field x expresses the effect of the noise on the dynamics of x. In the small noise T → 0 limit, the path that dominates the integral minimizes S, and the resulting saddle-point trajectories follow a Hamiltonianlike dynamics (with respect to H in Eq. ( 17)) in a space of extended variables where x stands for a momentum. The escape time τ then satisfies the logarithmic equivalence

τ e φ T φ ≡ min x(t)|x(-∞)=xm, x(+∞)=xs S [x, x] (18)
The explicit expression for the Hamiltonian H(x, x) is

H(x, x) ≡ -x • (1 + Γ) • ∇U (x) + x2 (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

ẋ = ∂H ∂x = -(1 + Γ) • ∇U (x) + 2x ẋ = -∂H ∂x = (1 + Γ) • Hess U (x) • x (20)
together with the condition H(x, x) = 0. Here Hess U (x) ab ≡ V (x (a) )δ ab is the Hessian matrix of U . These coupled differential equations are equivalent to

ẋ = (1 -Γ) • ∇U (x) x = ∇U (x) (21)
which, in principle, can be solved explicitly if a potential V is specified. We solve the instanton trajectory ẋ within perturbation theory up to second order in γ as a functional of the external potential V , assuming the γ = 0 trajectory is fully known. A linear analysis of the trajectory close to the saddle point x s shows (see App. A) that the perturbation expansion must have the form

x (1) = x (1) 0 + γ 2 x (1) 2 + O(γ 4 ) x (2) = x (2) 0 + γx (2) 1 + O(γ 3 ) (22) with ẋ(1) 0 = V x (1)
0 (t) and x

(2) 0 = 0. By taking into account the boundary condition for the perturbation x(+∞)

(a) i = x(-∞) (a) i = 0 for a = 1, 2, i > 0 we obtain x (2) 1 = -e kmt +∞ t dτ e -kmτ V x (1) 0 (τ ) (23)
x

(1) 2 = k m e G(t) +∞ t dτ e -G(τ ) x (2) 1 (τ ) + C (24)
FIG. 1. Instanton trajectory from Eq. ( 21) for γ = 0 and its perturbative expansion up to second order in γ for the double well potential

V (x) = 1 4 x 4 -1 2 x 2 .
where

G(t) ≡ dτ V (x (1)
0 (τ )) and C is a constant determined by matching the asymptotic expression for x (1) 2 (t) for t → +∞ to the linearized expansion at x s . We illustrate the role played by γ in deforming the equilibrium trajectory in Fig. 1. As done in previous studies of the escape time out of equilibrium [START_REF] Bray | Instanton calculation of the escape rate for activation over a potential barrier driven by colored noise[END_REF], we have chosen as an explicit potential V the well-known double well

V (x) ≡ 1 4 x 4 -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. Interestingly, the action for the instanton trajectory is unchanged with respect to the equilibrium dynamics, namely

S [x, x] = +∞ -∞ dt ẋ • ∇U (x) = ∆U ( 25
)
where ∆U = V (x M ) -V (x m ). That this result is identical to the one obtained with the equilibrium dynamics at γ = 0 supports the the heuristic analysis of Sec. III A.

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 γ on the prefactor to the Arrhenius exponential.

C. Full expression of τ

Bouchet and Reygner [START_REF] Bouchet | Generalisation of the eyringkramers transition rate formula to irreversible diffusion processes[END_REF] 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 γ-dependence of the escape time. When applied to the IO dynamics the time τ reads

τ = 2π λ + k M k m exp ∆V T (26)
where k M ≡ -V (x M ) > 0, k m ≡ V (x m ) > 0, and λ + is the modulus of the eigenvalue related to the unstable direction of the dynamics linearized around the saddle point x s :

λ + = 1 2 k M -k m + (k M + k m ) 2 + 4γ 2 k M k m (27)
For γ = 0 Eq. ( 26) yields the well-known Kramers' time

τ = 2π √ k M k m exp ∆V T (28)
while for γ = 0 we have λ + > 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 A was physically sound. It can indeed be recovered by letting k m → +∞ in Eq. ( 27), yielding λ + ≈ k M (1 + γ 2 ). This limiting case is actually the smallest achievable value of τ at given γ, k M and ∆V , since ∂ km λ + is always positive.

The results of this section precisely quantify the gain in speed that γ 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 γ-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 meanfield 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 ≥ 3.

IV. p-SPIN WITH ICHIKI-OHZEKI DYNAMICS

The spherical p-spin is a disordered fully connected model of N interacting spins σ 1 , . . . , σ N fulfilling the constraint

N i=1 σ 2 i = N . The energy of a configuration is H [{σ i } i=1,...,N ] = - i1<...<ip J i1...ip σ i1 . . . σ ip ( 29
)
where the J i1...ip are independent quenched random couplings with variance 2J 2 p!N p-1 . To cut a long story short, we refer the reader to the huge body of literature [START_REF] Crisanti | The spherical p-spin interaction spin glass model: the statics[END_REF][START_REF] Crisanti | The sphericalp-spin interaction spin-glass model[END_REF][START_REF] Barrat | Dynamics within metastable states in a mean-field spin glass[END_REF][START_REF] Barrat | The p-spin spherical spin glass model[END_REF][START_REF] Castellani | Spin-glass theory for pedestrians[END_REF][START_REF] Dominicis | Random fields and spin glasses: a field theory approach[END_REF] 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 [START_REF] Berthier | Non-equilibrium glass transitions in driven and active matter[END_REF]). 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 [START_REF] Crisanti | The spherical p-spin interaction spin glass model: the statics[END_REF]. 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 [START_REF] Crisanti | The sphericalp-spin interaction spin-glass model[END_REF]. When endowed with purely relaxational equilibrium dynamics, the spin-spin correlation function C(t, t ) = 1 N i σ i (t)σ i (t ) , in the regime with t large and t -t finite with time-translation invariance, fails to relax to 0 for T < T d , reaching instead a temperaturedependent plateau value q. The p-spin is thus the ideal test-bench to quantify the acceleration of the IO dynamics.

A. 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 σ (1) i and σ

(2) i evolving according to

∂ t σ (1) i = - ∂H ∂σ (1) i + γ ∂H ∂σ (2) i -µ (1) (t)σ (1) i + γµ (2) (t)σ (2) i + √ 2T η (1) i ( 30
)
∂ t σ (2) i = - ∂H ∂σ (2) i -γ ∂H ∂σ (1) i -µ 2 (t)σ (2)
i -γµ (1) (t)σ

(1) i + √ 2T η (2) i (31)
where the functions µ (a) (t) are adjusted so that the spherical constraint

N i=1 σ (a) i (t) 2
= 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 [START_REF] Cugliandolo | Glassy behaviour in disordered systems with nonrelaxational dynamics[END_REF], or in a shear [START_REF] Berthier | A two-timescale, two-temperature scenario for nonlinear rheology[END_REF] or with an active force [START_REF] Berthier | Non-equilibrium glass transitions in driven and active matter[END_REF]) 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 B ∼ e -βH for H = H [{σ (1) i }]+H [{σ (2)
i }] (with 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 ) = σ (a) (t)σ (b) (t ) , and the same holds for the response function R ab (t, t ) and the integrated response F ab (t, t ) ≡ -t t dτ R ab (t, τ ) . In practice one is able to write an effective Langevin equation for the dynamics of a single spin σ (1) in system 1 and its counterpart σ (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 ) = C ab (t -t ) and F ab (t, t ) = F ab (t -t ), and by applying standard path integral techniques we arrive at the following Langevin equation for σ(t) ≡ σ (1) σ (2) )

T and the dynamical equations for the correlation matrix C and the integrated response matrix F:

∂ t σ(t) = -(1 + Γ) • µ(t) • σ(t) + t 0 dτ M F (t -τ ) • σ(τ ) + 1 T (1 + Γ) • C 0 (p, t) • σ(0) + ξ(t) (32)
∂ t F(t) = -1 -(1 + Γ) • µ(t) • F(t) + t t dτ M F (t -τ ) • F(τ ) ( 33
)
∂ t C(t) = -(1 + Γ) • µ(t) • C(t) + t 0 dτ M F (t -τ ) • C(τ ) + 1 T (1 + Γ) • C 0 (p, t) (34)
µ(t) = Diag T 1 2 + t 0 dτ M(t -τ ) • C T (t -τ ) + t 0 dτ D(t -τ ) • ∂ τ F T (t -τ ) + 1 T (1 + Γ) • C 0 (p, t) • C T (t) (35)
where the Gaussian noise has correlations ξ(t) ⊗ ξ(t

) T = 2T 1 2 δ(t -t ) + D(t, t ).
The three kernels M F , C 0 and D that appear in Eqs. Eqs. [START_REF] Bray | Instanton calculation of the escape rate for activation over a potential barrier driven by colored noise[END_REF][START_REF] Crisanti | The spherical p-spin interaction spin glass model: the statics[END_REF][START_REF] Crisanti | The sphericalp-spin interaction spin-glass model[END_REF] are nonlinear functions of the elements of C and 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 C(0) = 1, F(0) = 0, while the explicit expression for the diagonal matrix µ(t) enforces the spherical constraint.

In the case of standard equilibrium dynamics, the response and the correlations are related by the fluctuationdissipation theorem (FDT),

∂ t F(t) = -1 T ∂ t 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. [START_REF] Bray | Instanton calculation of the escape rate for activation over a potential barrier driven by colored noise[END_REF][START_REF] Crisanti | The spherical p-spin interaction spin glass model: the statics[END_REF][START_REF] Crisanti | The sphericalp-spin interaction spin-glass model[END_REF][START_REF] Barrat | Dynamics within metastable states in a mean-field spin glass[END_REF] 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.

B. Response and fluctuations

Motivations

The entropy production is a standard way of measuring the distance to equilibrium. Here we find, using Eq. [START_REF] Marinari | Simulated tempering: a new monte carlo scheme[END_REF], that Σ = 2γ 2 N µ [START_REF] Barrat | The p-spin spherical spin glass model[END_REF] so that the time scale at which nonequilibrium effects come into play goes as 1 γ 2 . It is tempting to believe that the distance to equilibrium can generically be viewed as being of order O(γ 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 A). Consider the two coupled equations for r = (x 1 , x 2 )

ṙ = -(1 + Γ)kr + √ 2T β (37)
The evolution operator has a spectrum of relaxation rates -k(n 1 + n 2 + iγ(n 1 -n 2 )) (with n 1 , n 2 ∈ N) and we realize that harmonic oscillators are marginal with re-spect 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 Σ = 2kγ 2 > 0. The correlation matrix C(t, t ) = r(t) ⊗ r(t ) and the response function

R(t, t ) = δ r(t) δf (t ) =0
to an external perturbing force f are related, in the stationary state, by

R(t, t )(1 + Γ) = - 1 T ∂ t C(t, t ) (38)
which can be cast in an integral form by introducing the integrated response F(t, t ) ≡ -t t dτ R(t, τ ) as:

F(t, t )(1 + Γ) = 1 T (C(t, t ) -C(t , t )) (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 F and 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.

Numerics and aFDT

We integrate the dynamical equations Eqs. [START_REF] Crisanti | The spherical p-spin interaction spin glass model: the statics[END_REF][START_REF] Crisanti | The sphericalp-spin interaction spin-glass model[END_REF][START_REF] Barrat | Dynamics within metastable states in a mean-field spin glass[END_REF]) for p = 3 using a finite difference method implemented in [START_REF] Fuchs | Comments on the alpha-peak shapes for relaxation in supercooled liquids[END_REF][START_REF] Kim | The dynamics of the spherical pspin model: From microscopic to asymptotic[END_REF][START_REF] Berthier | Spontaneous and induced dynamic correlations in glass formers. ii. model calculations and comparison to numerical simulations[END_REF][START_REF] Folena | The mixed p-spin model: selecting, following and losing states[END_REF] 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 ∞ at long times in the ergodic phase. This is a static quantity independent from γ, for which a value of

E ∞ = -1
2T is predicted [START_REF] Crisanti | The spherical p-spin interaction spin glass model: the statics[END_REF]. In the numerical integration, we computed the energy by exploiting the relation

µ(t → +∞) = (T -pE ∞ ) 1 (40)
The agreement between the measured value of E ∞ and its theoretical prediction is shown in Fig. 2. The errors are of order 10 -12 and they are therefore negligible.

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 γ up to 2.4, the IO dynamics is accelerated by a non negligible amount compared to the standard one. For γ ≈ 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 γ = 0, and the position of the infinite time plateau is left unchanged.

We quantify the speedup obtained by means of the IO dynamics by measuring the relaxation time τ α (γ), 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 τ α (γ) and the relaxation time in the standard dynamics, τ α (0) appears to scale approximatily as (1 + γ 2 ) -1 . This scaling the same as the one found in Sec. III A when heuristically discussing the Kramers' problem. We will provide the reader with further interpretation for this phenomenon in the next section.

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 cor- relations for various values of is shown. The dynamics of these quantities is entirely located in the region of β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 γ for γ 1.

The location of the maximum moves left to shorter times as γ grows, as expected from Eq. ( 7).

Guided by the analysis of the harmonic oscillator in Sec. IV B 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 γ 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 σ in Eq. [START_REF] Bray | Instanton calculation of the escape rate for activation over a potential barrier driven by colored noise[END_REF].

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 γ? In the fol- 39) from the integration of Eqs. [START_REF] Crisanti | The spherical p-spin interaction spin glass model: the statics[END_REF][START_REF] Crisanti | The sphericalp-spin interaction spin-glass model[END_REF][START_REF] Barrat | Dynamics within metastable states in a mean-field spin glass[END_REF] at T = 0.613 and γ = 2.4.

lowing section we will build a theoretical understanding of our numerical findings.

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. [START_REF] Crisanti | The sphericalp-spin interaction spin-glass model[END_REF][START_REF] Barrat | Dynamics within metastable states in a mean-field spin glass[END_REF] 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. [START_REF] Crisanti | The spherical p-spin interaction spin glass model: the statics[END_REF][START_REF] Crisanti | The sphericalp-spin interaction spin-glass model[END_REF][START_REF] Barrat | Dynamics within metastable states in a mean-field spin glass[END_REF] to a single equation for the correlation matrix C(t):

∂ t C(t) = -T (1 + Γ) • C(t) - 1 T (1 + Γ) t 0 dτ C 0 (p, t -τ ) • ∂ τ C(τ ) ( 41
)
with C 0 (p, t) ij ≡ p 2 C p-1 ij . Apart from the matrix prefactor (1 + Γ) the structure of this equation is the same as the one obtained in the standard equilibrium dynamics of the p-spin [START_REF] Crisanti | The sphericalp-spin interaction spin-glass model[END_REF][START_REF] Barrat | The p-spin spherical spin glass model[END_REF][START_REF] Castellani | Spin-glass theory for pedestrians[END_REF]. We therefore study its long time behavior by assuming that lim t→+∞ C(t) = Q, with Q a matrix of long time plateaus generalizing the nonergodic parameter. We introduce the Laplace transform of a matrix A(t) as Ã(z) ≡ +∞ 0 dtA(t)e -zt and we observe that C(z)

≈ Q z and C0 (p, z) ≈ Q0(p) z as z → 0, with Q 0 (p) ij ≡ Q p-1 ij .
Applying the Laplace transform to both sides of Eq. ( 41) and keeping only the terms that diverge as z → 0 yields

T Q - 1 T Q 0 (p) • (1 -Q) = 0 (42)
We now see that the dependence on Γ has simply disappeared. The form of the plateaus matrix can therefore not depend on γ. In particular, it must match the solution obtained at γ = 0. We can therefore write Q = Q1 and Q 0 (p) = Q p-1 1. Eq. ( 42) then becomes the standard equation for the long time plateau of the p-spin equilibrium dynamics

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 (1 + Γ) 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 γ. In particular we will recover the scaling of the relaxation time τ α shown in Fig. 4 starting from the dynamical equation for the correlation matrix Eq. ( 41). We multiply both sides by 1 + Γ T and we expand C in powers of γ,

C ≈ δ 0 C + γδ 1 C + γ 2 δ 2 C, thereby obtaining: ∂ t C(t) + Γ T • ∂ t δ 1 C(t) = -(1 + γ 2 ) C(t) + 1 T t 0 dτ C 0 (p, t -τ ) • ∂ τ C(τ ) (44)
The second term on the left hand side contains δ 1 C, the first order perturbation of 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 B 2 we assume that this term essentially vanishes at times larger than the time scale at which the β 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 δ 1 C. We therefore obtain, to second order in γ:

∂ t C aa (t) = -(1 + γ 2 ) T C aa (t) + 1 T t 0 dτ C aa,0 (p, t -τ )∂ τ C aa (τ ) (45)
for a = 1, 2. This is the same equation as that of the equilibrium dynamics, up to a rescaling of time t → t/(1+γ 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 + γ 2 in time is gained, for not too large values of γ. 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 γ. At finite number N of spins, it would interesting to probe the interplay of γ with the system size, as in principle, a sufficiently large γ could overcome energy barriers. To adjust γ appropriately, one could rely on the analysis of Stariolo and Cugliandolo [START_REF] Stariolo | Activated dynamics of the ising p-spin disordered model with finite number of variables[END_REF][START_REF] Stariolo | Barriers, trapping times, and overlaps between local minima in the dynamics of the disordered ising p-spin model[END_REF], where the relevant time-scales are carefully sorted out. A related interesting question is how the instanton trajectories are modified by a nonzero γ 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. [START_REF] Ros | Dynamical instantons and activated processes in mean-field glass models[END_REF] 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 γ beyond O(1), it is unclear at which value of γ 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 γ 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 γ eff ∼ √ nγ would occur for too large n. This, in turn, raises the question of how to fine tune the optimal value of γ. As our results suggest, it is possible that the best acceleration is achieved by a nontrivial time-dependent protocol γ(t) (in which γ(→ +∞) → γ ∞ , the latter being possibly zero). This question is touched upon in [START_REF] Lelievre | Optimal nonreversible linear drift for the convergence to equilibrium of a diffusion[END_REF].

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 γ should be felt at all temperatures. How this approach fares as compared to other ones such as SWAP [START_REF] Grigera | Fast monte carlo algorithm for supercooled soft spheres[END_REF][START_REF] Berthier | Equilibrium sampling of hard spheres up to the jamming density and beyond[END_REF][START_REF] Ninarello | Models and algorithms for the next generation of glass transition studies[END_REF] is, in our opinion, and based on preliminary results [49], the most interesting open question.

+ γ ∂H ∂σ (2) i -µ (1) (t)σ (1)
i + γµ (2) (t)σ

(2) i + √ 2T ξ (1) i (B1)
∂ t σ (2) i = - ∂H ∂σ (2) i -γ ∂H ∂σ (1) i -µ (2) (t)σ (2)
i -γµ (1) (t)σ

(1)

i + √ 2T ξ (2) i (B2)
where µ (a) (t) are elastic constants that enforce the spherical constraints Note that the couplings J i1...ip have the same realization for both systems, and they are independently distributed according to a Gaussian distribution of mean 0 and variance p! 2N p-1 . We suppose that the systems start from an initial equilibrium condition at a temperature T ≡ β -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 [START_REF] Cugliandolo | Analytical solution of the off-equilibrium dynamics of a long-range spin-glass model[END_REF] and further detailed in [START_REF] Barrat | The p-spin spherical spin glass model[END_REF] or [START_REF] Castellani | Spin-glass theory for pedestrians[END_REF]. The starting point is the partition function for the dynamics of the composite system in Eq. (B1) 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

1 ≡ Z = Dσ 0 1 Z eq e -βH [σ (1) 0 ,σ (2) 0 ]
× DσD σ exp -

N i=1 t 0 dτ σ(1) i ∂ τ σ (1) i + ∂H ∂σ (1) i -γ ∂H ∂σ (2) i + µ (1) σ (1)
i -γµ (2) σ

(2) i + σ(2) i ∂ τ σ (2) i + ∂H ∂σ (2) i + γ ∂H ∂σ (1) i + µ (2) σ (2)
i + γµ (1) σ

(2) i -T σ(1) i 2 -T σ(2) i 2 (B4)
with Z eq the static partition function of the system, Z eq ≡ Dσ 0 e -βH [σ0] . We are using a vector notation for the spin variables, σ(t) = σ (1) (t), σ (2) (t) T , the auxiliary field σ(t) = σ(1) (t), σ(2) (t) T and the initial condition σ 0 (t) = σ (1) (0), σ (2) (0) T As noted in [START_REF] Houghton | Role of initial conditions in the mean-field theory of spin-glass dynamics[END_REF], since we are using thermalized initial conditions, we must resort to the replica trick in order to average Z over the disorder:

1 Z eq = lim n→0 Z n-1 eq (B5)
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 σ

(a) • σ (b) (t, t ) ≡ N i=1 σ (a) i (t)σ (b)
i (t ), with a and b indices identifying the two systems, we can write in the thermodynamic limit

Z = Dσ(0)DσD σ exp - N i=1 dτ σ(1) i ∂ τ σ (1)
i µ (1) σ

(1)

i -γµ (2) σ (2) i -T σ(1) i 2 + σ(2) i ∂ τ σ (2) i + µ (2) σ (2)
i + γµ (1) σ

(1) i -T σ(2) i 2 + 1 4N p-1 dt dt I 1 (t, t ) + 1 4T N p-1 dtI 2 (t) (B6) With I 1 (t, t ) ≡ p σ(1) • σ(1)
σ (1) • σ (1) p-1 + p(p -1) iσ (1) • σ (1) σ (1) • σ(1) σ (1) • σ (1) p-2

-γ p σ(1) • σ(1) σ (1) • σ (2) p-1 + p(p -1) σ(1) • σ (2) σ (1) • σ(1) σ (1) • σ (2) p-2 + p σ(1) • σ(2) σ (1) • σ (2) p-1 + p(p -1) σ(1) • σ (2) σ (1) • σ(2) σ (1) • σ (2) p-2 + γ p σ(1) • σ(2) σ (1) • σ (1) p-1 + p(p -1) σ(1) • σ (1) σ (1) • σ(2) σ (1) • σ (1) p-2 -γ p σ(1) • σ(1) σ (2) • σ (1) p-1 + p(p -1) σ(1) • σ (1) σ (2) • σ(1) σ (2) • σ (1) p-2 + γ 2 p σ(1) • σ(1) σ (2) • σ (2) p-1 + p(p -1) σ(1) • σ (2) σ (2) • σ(1) σ (2) • σ (2) p-2 -γ p σ(1) • σ(2) σ (2) • σ (2) p-1 + p(p -1) σ(1) • σ (2) σ (2) • σ(2) σ (2) • σ (2) p-2 -γ 2 p σ(1) • σ(2) σ (2) • σ (1) p-1 + p(p -1) σ(1) • σ (1) σ (2) • σ(2) σ (2) • σ (1) p-2 + (1 ↔ 2, γ ↔ -γ) (B7) and 1 2 I 2 (t, 0) ≡ p σ(1) • σ (1) 0 σ (1) • σ (1) 0 p-1 -γp σ(1) • σ (1) 0 σ (2) • σ (1) 0 p-1 + p σ(1) • σ (2) 0 σ (1) • σ (2) 0 p-1 -γp σ(1) • σ (2) 0 σ (2) • σ (2) 0 p-1 + (1 ↔ 2, γ ↔ -γ) (B8)
We now exploit the thermodynamic limit N → ∞ 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:

Z = DQΠ a,b=1,2 δ N Q (ab) 1 (t, t ) -σ(a) • σ(b) δ N Q (ab) 2 (t, t ) -σ (a) • σ (b) × δ N Q (ab) 3 (t, t ) -σ(a) • σ (b) δ N Q (ab) 4 (t, t ) -σ (a) • σ(b) × δ N Q (ab) 5 (t) -σ(a) • σ (b) 0 × δ N Q (ab) 6 (t) -σ (a) • σ (b) 0 × (r.h.s. Eq. (B6)) = DP DQΠ a,b=1,2 exp iN dtdt P (ab) 1 (t, t )Q (ab) 1 -P (ab) 1 (t, t )σ (a) • σ(b) × exp iN dtdt P (ab) 2 Q (ab) 2 -P (ab) 2 σ (a) • σ (b) × exp iN dtdt P (ab) 3 (t, t )Q (ab) 3 -P (ab) 3 σ(a) • σ (b) × exp iN dtdt P (ab) 4 Q (ab) 4 -P (ab) 4 σ (a) • σ(b) × exp iN dt P (ab) 5 Q (ab) 5 -P (ab) 5 σ(a) • σ (b) 0 × exp iN dt P (ab) 6 Q (ab) 6 -P (ab) 6 σ (a) • σ (b) 0 × (r.h.s of Eq. (B6)) (B9)
The dynamical overlaps have a natural interpretation if one introduces the correlation and the response matrices C(t, t ) and R(t, t ):

C ab (t, t ) ≡ 1 N i σ (a) i (t)σ (b) i (t ) (B10) R ab (t, t ) ≡ 1 N i σ (a) i (t)σ (b) i (t ) (B11)
We therefore see that in the thermodynamic limit the following self consistent relations hold:

               Q (ab) 2 (t, t ) = C ab (t, t ) Q (ab) 3 (t, t ) = R ab (t , t) Q (ab) 4 (t, t ) = R ab (t, t ) Q (ab) 5 (t) = C ab (t, 0) Q (ab) 6 (t) = R ab (t, 0) (B12) Note that Q (ab) 3 (t, t ) = Q (ba) 4
(t , t). Moreover, we claim that, as in the standard case [START_REF] Sompolinsky | Relaxational dynamics of the edwards-anderson model and the mean-field theory of spin-glasses[END_REF], Q

= 0. The saddle point equations are obtained by differentiating the argument of the exponentials in Eq. (B6) with respect to the set

= p 2T Q (11) 5 -γ p 2T Q (21) 5 iP (11) 6 = 0 iP (12) 1 = p 4 Q (12)p-1 2 -γ p 4 Q (22)p-1 2 -Q (11)p-1 2 -γ 2 p 4 Q (22)p-1 2 iP (12) 2 = 0 iP (12) 3 = p 4 (p -1)Q (12) 4 Q (12)p-2 2 -γ p 4 (p -1)Q (11) 4 Q (12)p-2 2 -γ p 4 (p -1)Q (22) 4 Q (22)p-2 2 + γ 2 p 4 (p -1)Q (21) 4 Q (22)p-2 2 iP (12) 4 = p 4 (p -1)Q (12) 3 Q (12)p-2 2 + γ p 4 (p -1)Q (11) 3 Q (11)p-2 2 + γ p 4 (p -1)Q (22) 3 Q (12)p-2 2 + γ 2 p 4 (p -1)Q (21) 3 Q (11)p-2 2 iP (12) 5 = p 2T Q (12) 5 -γ p 2T Q (22) 5 iP (12) 6 = 0 iP (21) 1 = p 4 Q (21)p-1 2 -γ p 4 Q (11)p-1 2 -Q (22)p-1 2 -γ 2 p 4 Q (12)p-1 2 iP (21) 2 = 0 iP (21) 3 = p 4 (p -1)Q (21) 4 Q
(21)p-2 2 + γ p 4 (p -1)Q (11) 4 Q
(11)p-2 2 + γ p 4 (p -1)Q (22) 4 Q (21)p-2 2 + γ 2 p 4 (p -1)Q (12) 4 Q (11)p-2 2 iP (21) 4
= p 4 (p -1)Q

(21)p-2 2 -γ p 4 (p -1)Q (21) 3 Q
(21)p-2 2 -γ p 4 (p -1)Q (11) 3 Q
(12)p-2 2 + γ 2 p 4 (p -1)Q (22) 3 Q
(22)p-2 2 iP (21) 5 = p 2T Q (21) 5 + γ p 2T Q (12) 3 Q
= 0

iP (22) 1 = p 4 Q (22)p-1 2 + γ p 4 Q (12)p-1 2 + Q (21)p-1 2 + γ 2 p 4 Q (11)p-1 2 iP (22) 2 = 0 iP (22) 3 = p 4 (p -1)Q (22) 4 Q (22)p-2 2 -γ p 4 (p -1)Q (21) 4 Q (22)p-2 2 + γ p 4 (p -1)Q (12) 4 Q (12)p-2 2 -γ 2 p 4 (p -1)Q (22) 4 Q (12)p-2 2 iP (22) 4 = p 4 (p -1)Q (22) 3 Q (22)p-2 2 -γ p 4 (p -1)Q (12) 3 Q (22)p-2 2 + γ p 4 (p -1)Q (21) 3 Q (21)p-2 2 -γ 2 p 4 (p -1)Q (22) 3 Q (21)p-2 2 iP (22) 5 = p 2T Q (22) 5 + γ p 2T Q (12) 5 iP (22) 6
= 0 (B13) Note that iP 2 = iP 6 = 0 because Q 1 = 0 and because of causality, which implies that Q (t, t ) = 0. Substitution of the saddle point expressions of Eq. (B13) in Eq. (B9), together with Eq. (B12) yields the effective dynamical equations for the spins of the system. Using a vector notation, σ(t) = σ (1) σ (2) T , we obtain Since the dynamics of the system starts within the steady state, we assume that time translation invariance holds: C(t, t ) ≡ C(t-t ) and R(t, t ) ≡ R(t-t ) The equations of motion for the response and correlation matrices therefore read, when taking t = 0

∂ t σ(t) = -(1 + Γ) • µ(t) • σ(t) + t 0 dτ M R (t, τ ) • σ(τ ) + 1 T (1 + Γ) • C 0 (p, t) • σ(0) + ξ(t) (
∂ t R(t) = -(1 + Γ) • µ(t) • R(t) + t t dτ M R (t -τ ) • R(τ ) + δ(t)1 2 (B25) ∂ t C(t) = -(1 + Γ) • µ(t) • C(t) + t 0 dτ M R (t -τ ) • C(τ ) + 1 T (1 + Γ) • C 0 (p, t) (B26) µ(t) = Diag T 1 2 + t 0 dτ M(t -τ ) • C T (t -τ ) + t 0 dτ D(t -τ ) • R T (t -τ ) + 1 T (1 + Γ) • C 0 (p, t) • C T (t) (B27)
For numerical purposes it is more convenient to work with the integrated response matrix F(t, t ) [START_REF] Kim | The dynamics of the spherical pspin model: From microscopic to asymptotic[END_REF], defined by the companion initial conditions are C(0) ≡ 1, F(0) = 0. Together with Eq. (B14), these equations correspond to Eqs. [START_REF] Bray | Instanton calculation of the escape rate for activation over a potential barrier driven by colored noise[END_REF][START_REF] Crisanti | The spherical p-spin interaction spin glass model: the statics[END_REF][START_REF] Crisanti | The sphericalp-spin interaction spin-glass model[END_REF][START_REF] Barrat | Dynamics within metastable states in a mean-field spin glass[END_REF] presented in the main text.

F
In the following section we provide the reader with some details on the algorithm implemented to integrate the dynamical equations. putting everything back in Eq. (D4) and simplifying we get that

1 T ∂ t C(t) = ∂ t F(t) • (1 + Γ) (D10)
with ∂ t F(t) given by Eq. (B29), and the proof is completed.

6 FIG. 3 .

 63 FIG.3. Correlation function of system 1, C11(t) for γ = 2.4 (solid lines) and γ = 0 (dashed lines) for various temperatures above and below T d ≈ 0.6124. The dotted curves represent the cross correlations C12(t) for γ = 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.

FIG. 4 .

 4 FIG. 4. Ratio between the relaxation time of C11 for the IO dynamics, τα(γ), and the equilibrium dynamics τα(0) for temperatures close to T d .

FIG. 5 .

 5 FIG. 5. Cross correlation C12 (dotted lines) and C21 (dashed lines)as a function of time for T = 0.613 at different values of γ.

12 FIG. 6 .

 126 FIG. 6. Parametric plot of the left and right hand side of the first entry of Eq. (39) from the integration of Eqs.[START_REF] Crisanti | The spherical p-spin interaction spin glass model: the statics[END_REF][START_REF] Crisanti | The sphericalp-spin interaction spin-glass model[END_REF][START_REF] Barrat | Dynamics within metastable states in a mean-field spin glass[END_REF] at T = 0.613 and γ = 2.4.

2 (

 2 t) = N for a = 1, 2, ξ (a) i are independent Gaussian white noises ξ

j 2 σ

 2 (t ) = δ ab δ ij δ(t -t ) and the total Hamiltonian of the system H readsH {σ i } (a) i=1,...,N,a=1,2 = -i1<...<ip J i1...ip a=1,

  ) and the coupling with the initial conditions C 0 (t, 0) are the following: 11 (R 11 + γR 12 ) -γC p-2 21 (R 21 + γR 22 ) C p-2 12 (R 12 -γR 11 ) -γC p-2 22 (R 22 -γR 21 ) C p-2 21 (R 21 + γR 22 ) + γC p-2 11 (R 11 + γR 12 ) C p-2 22 (R 22 -γR 21 ) + γC p-2 12 (R 12 -γR 11 )D(t, t ) ≡ (1 + Γ) • C 0 (p, t) • 1 + Γ T (t)σ (b) (t ) = 2T R ab (t , t) + dτ D(t, τ ) • R T (t , τ ) ab

	M R (t, τ ) =	p(p -1) 2	 C p-2 	 
					(B16)
					(B17)
	C 0 (p, t, 0) ij =	p 2	C p-1 ij (t, 0)	(B18)
	From Eq. (B14) by using the definition (B10) and the equalities
					R ab (t, t ) =	δσ (a) (t) δξ (b) (t )	(B19)
				ξ (a) (B20)
	we obtain
			∂ N i=1	(a) i σ
					B14)
	where µ ab (t) ≡ δ ab µ (a) (t), Γ ≡	0 -γ γ 0	. The correlations of the noise ξ(t) are

ξ(t)ξ(t ) = 2T δ(t -t )1 2 + D(t, t ) (B15)

while the expression for the memory kernel M R (t, t ), the noise correlations D(t,

t t R(t, t ) = -(1 + Γ) • µ(t) • R(t, t ) + t t dτ M R (t, τ ) • R(τ, t ) + δ(t -t )1 2 (B21) ∂ t C(t, t ) = -(1 + Γ) • µ(t) • C(t, t ) + t 0 dτ M R (t, τ ) • C(τ, t ) + 2T R(t , t) (B22) + t 0 dτ D(t, τ ) • R T (t , τ ) + 1 T (1 + Γ) • C 0 (p, t) • C(0, t )

By imposing the spherical constraint

2 (t) = N for a = 1,

2 one can find a self-consistent expression for the restoring forces µ(t). If we denote the diagonal part of a matrix by Diag(A) ab ≡ δ ab A ab the spherical constraint condition reads Diag(C(t, t)) = 1. By differentiating with respect to t we get lim t →t Diag (∂ t C(t, t ) + ∂ t C(t, t )) = lim t →t Diag ∂ t C(t, t ) + ∂ t C T (t , t) = 0 (B23) which leads to µ(t) = Diag T 1 2 + t 0 dτ M(t, τ ) • C T (t, τ ) + t 0 dτ D(t, τ ) • R T (t, τ ) + 1 T (1 + Γ) • C 0 (p, t) • C T (t, 0) (B24)

  since it has a smoother relaxation than R(t, t ). The equations of motion in the time translation invariant regime are∂ t F(t) = -1 -(1 + Γ) • µ(t) • F(t) + ∂ t C(t) = -(1 + Γ) • µ(t) • C(t) + 11 (∂ τ F 11 + γ∂ τ F 12 ) -γC p-2 21 (∂ τ F 21 + γ∂ τ F 22 ) C p-2 12 (∂ τ F 12 -γ∂ τ F 11 ) -γC p-2 22 (∂ τ F 22 -γ∂ τ F 21 ) C p-2 21 (∂ τ F 21 + γ∂ τ F 22 ) + γC p-2 11 (∂ τ F 11 + γ∂ τ F 12 ) C p-2 22 (∂ τ F 22 -γ∂ τ F 21 ) + γC p-2 12 (∂ τ F 12 -γ∂ τ F 11 )

							t
							dτ M F (t -τ ) • F(τ )	(B29)
							t
					0	t	dτ M F (t -τ ) • C(τ ) +	1 T	(1 + Γ) • C 0 (p, t)	(B30)
							t
		µ(t) = Diag T 1 2 +		dτ D(t -τ ) • ∂ τ F T (t -τ )
							0	(B31)
			+	1 T	(1 + Γ) • C 0 (p, t) • C T (t)
	and the memory kernel M F (t) that reads
	M F (t -τ ) ≡ -	p(p -1) 2	×
	 C p-2					
							
							(B32)

ab (t, t ) ≡ -t t dτ R ab (t, τ ) (B28) t 0 dτ M F (t -τ ) • C T (t -τ ) +

ACKNOWLEDGMENTS

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 x s . The first line of Eq. ( 21) reads

This linear systems has one unstable and one stable direction along which the dynamics takes the form C ± e λ±t v ± ,with λ ± , v ± the eigenvalues and eigenvectors of A and C ± determined by the initial condition. Both the eigenvalues and the eigenvectors depend on γ. The t → +∞ limit of the instanton dynamics in the first line of Eq. (A1) corresponds to the dynamics of system along its stable direction v -. By expanding the latter around γ = 0 we obtain

which justifies the form of the perturbation expansion adopted in Eq. [START_REF] Duncan | Variance reduction using nonreversible langevin samplers[END_REF].

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

and they read:

Appendix C: Numerical Integration of the dynamical equations

We present here the integration scheme adopted to numerically solve the dynamical equations of motion, Eqs. (B29, B30, B31). The basic idea for the algorithm was first given in [START_REF] Fuchs | Comments on the alpha-peak shapes for relaxation in supercooled liquids[END_REF] in the context of the mode coupling equations. It was then applied to the aging of the p-spin model in [START_REF] Kim | The dynamics of the spherical pspin model: From microscopic to asymptotic[END_REF] and generalised in [START_REF] Berthier | Spontaneous and induced dynamic correlations in glass formers. ii. model calculations and comparison to numerical simulations[END_REF] for the case of non stationary systems. A pedagogical description of the algorithm can be found in [START_REF] Folena | The mixed p-spin model: selecting, following and losing states[END_REF].

To integrate Eqs. (B29, B30, B31), we take a grid of M points with a discretisation step ∆t, so that any function f (t) can be discretised as f k ≡ f (t k ) with t k ≡ k∆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 V i ≡ C i F i of the form

where NL is a generic nonlinear operator. This set of equation can be solved iteratively starting from the initial condition V n+1 = 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 C and 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 V 2i → V i for i = 0, ..., M -1 and the time step is rescaled by a factor two, ∆t → 2∆t. The propagation step is then repeated starting from the point n + 1 = 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. (B29, B30) the time derivative of a function at step k + 1, ∂ t f k+1 can be approximated as

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

and they can be approximated as

For the results presented in the main text we have adopted a grid of M = 1024 points and an initial timestep of ∆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:

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. (B29, B30, B31). The strategy of the proof is the following: starting from Eqs. (B30, B31) and the aFDT we will obtain Eq. (B29). To do so, we first rewrite Eq. (D1) in a differential form:

using Eq. (D2) we can rewrite the memory kernel M F in a simpler form:

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

to evaluate the time integral, we can use Eq. (D3), integration by parts, the aFDT relation and the boundary condition

So that µ reads after simplifying and taking the diagonal part:


S Sy

A Fabbri

Isaline Gravaud

D Seyedi

Evaluation of the CO 2 leakage risk along the abandoned wells in the French context

Keywords: classification, leakage risk, wellbore integrity, CO 2 storage

 one, to estimate the leakage risk along the abandoned wells in the French context. The selected criteria of this classification are the wellbore type (with/without casing that intercepts the aquifer where the storage will take place), the abandonment date (before/after 2000 which corresponds to a significant improvement in the French wellbore abandonment regulation), the top level of primary cementation, the deviation of the wellbore, and the drilling date (before/after 1980). The input data necessary to establish the classification are mostly available. When a data is missing, we have considered either: 1) a most prejudicial option assuming that the maximum value is considered; 2) or less prejudicial option where the criterion does not affect the calculated risk category. The results show a great variability of the risk when one moves from the most prejudicial option to the less prejudicial one. This emphasizes the primary importance of the management of the lack of data in this kind of study. Thus, without any additional information, the evaluation of the CO 2 leakage risk along the abandoned wells should at least exhibit the two extreme values that respectively correspond to the less and the most prejudicial options.

Introduction

A preliminary step of any CO 2 storage project consists in selecting the most appropriate location to inject the CO 2 . This step is fundamental as it will strongly influence the technical and economical evaluation of the project, and thus can affect its success.

As previously discussed by Grataloup [START_REF] Grataloup | A site selection methodology for CO 2 underground storage in deep saline aquifers : case of the Paris Basin[END_REF], from the several targets that influence this selection, the safety of the storage is particularly important. Consequently, for this last point, the scanning of the abandoned wellbore distribution that can be affected by the potential storage is necessary. Indeed, the injected CO 2 could be in contact with abandoned wellbores penetrating the aquifers, and depending on the quality of the cementing job, or on the type of abandonment procedure or on the external solicitation that have stressed the well, they can either be considered as a sustainable barrier to CO 2 leakage or as a preferential leakage pathway to upper fresh water aquifers or to the atmosphere (e.g. Wertz et al. [START_REF] Wertz | CO2 -Thermodynamic Model in a Leaking Well[END_REF], Nordbotten et al. [START_REF] Nordbotten | Model for CO 2 leakage including geological layers and multiple leaky wells[END_REF]). However, in these preliminary studies, it should not be economically feasible to quantify precisely (with geophysical tools for example) the leakage potential of each wellbore near each scanned storage location. Actually, in order to select the optimal storage location, it should be sufficient to have an estimation of the leakage potential through a classification based on data that can be easily found in official reports.

To do so, we relied on the classification of Watson et al. [START_REF] Watson | Evaluation of the potential for gas and CO 2 leakage along wellbores[END_REF] on the risk of CO 2 leakage (shallow leakage criteria) in the Alberta sedimentary basin. Then we have applied this classification to the French data base "deep wells" of Vernoux et al. [START_REF] Vernoux | Inventaire et évaluation des risques engendrés par les forages profonds sur les nappes d'eau souterraine du bassin Seine-Normandie[END_REF], which collects information from drilling and abandonment reports of 3482 wells within the Seine-Normandie basin.

This audit showed that the criteria are available in this data base with the exception of two criteria that are not available for all wells. To calculate the risk of CO 2 leakage from each well, numerical values reflecting the influence on the potential of wellbore to leak, are assigned to each criterion, and where a data is missing in the data base, two options are considered. In the first option, the maximum values are assigned to criteria, whereas in the second one minimum values are allocated to them. So the leakage risk, which is regarded as the product of risk criteria, is calculated for each well, and four risk levels are defined. A geographic information system (GIS) is then used to visualize the wells and their corresponding risk level.

Identification of the risk criteria

In order to build the classification, the first step is to identify the criteria that drive the risk of leakage. According to the studies from Khodadadi [START_REF] Khodadadi | Investigation of gas migration in Khangiran Wells[END_REF] and Wojtanowics et al. [START_REF] Wojtanowics | Gas Flow in Wells After Cementing[END_REF], the operator, the well location, the material used (casing/cement type, etc.), the well geometry, etc. seem to be the key items. However, the main problem is the quantification of the impact of these criteria on the risk of leakage. On that point, an important work was made by Watson et al. [4;[START_REF] Watson | Identification of Wells with high CO 2 -Leackage Potential in Mature Oil Fields Developed for CO 2 -Enhanced Oil Recovery[END_REF], where the impact of these criteria is evaluated form the analysis of the leakage occurrence from Gas Migration and Surface Casing Vent Flow Measurements on more than 500 wells in the Alberta sedimentary basin (Canada). In the following, this classification system is adapted to the French context. The Dogger formation is considered as the target reservoir and the following criteria are taken into account:

Well type (cased well or not): We consider that a well is cased if it has a cased interval which intercepts the reservoir (Dogger). Indeed, the casing intercepting the reservoir is interpreted by Watson et al. [START_REF] Watson | Evaluation of the potential for gas and CO 2 leakage along wellbores[END_REF] as additional CO 2 leakage sources (through the interfaces), increasing occurrences of SCVF (Surface Casing Vent Flow) and GM (Gas Migration) Abandonment date: before/after 2000, which corresponds to a significant improvement in the French wellbore abandonment regulation Cement to surface: this criterion is assigned to a well if it has at least one cemented interval that reaches the surface (z = 0) Well deviation: A well is deviated if its total vertical depth (TVD) is less than its total depth (TD) Spud date: this criterion that permits to take into account the age of the well, is compared to 1980 corresponding to the inflexion point of the distribution curve of the total number of wells between 1954 and 2003 [START_REF] Vernoux | Inventaire et évaluation des risques engendrés par les forages profonds sur les nappes d'eau souterraine du bassin Seine-Normandie[END_REF] To test the applicability of this classification, the database developed by Vernoux et al. [START_REF] Vernoux | Inventaire et évaluation des risques engendrés par les forages profonds sur les nappes d'eau souterraine du bassin Seine-Normandie[END_REF] under an agreement between BRGM and the French Water agency is used. This database is composed of 2590 oil wells, 119 geothermal wells, 435 natural gas wells and 338 water wells. The main table of this data base gathers the main characteristics of wells, namely, their identification number, type, cementing, casing, and plugging characteristics, their Lambert coordinates, depth, and dates (end of drilling, and abandonment). It is also related to other tables that complement the description of wells. Additional information about the fluids within the wells, the cement types, geological levels crossed by wells, cementing and casing descriptions, may be found in these secondary tables of this data base.

Within these wells, we need to select the abandoned wells "likely" to intercept a CO 2 storage site. In this study, we considered the aquifer of the Dogger (Jurassic), which is one of the two deep aquifers traditionally regarded as a target for CO 2 storage in the Paris Basin.

All wells reaching at least the Middle Jurassic are supposed to reach or impact the Dogger aquifer. However, only areas where the top of the Middle Jurassic is deeper than 800m have been considered. The limit of 800 m corresponds to the depth at which the temperature and pressure conditions keep the CO 2 in a supercritical state, guaranteeing the injectivity and storage capacity.

There are 216 wells in the data base that reach the Middle Jurassic and therefore expected to intercept the Dogger aquifer, which represents about 7% of the total number of deep wells. Table 1 summarizes the availability of information within the database to estimate if a risk criterion is reached or not.

Calculation of the risk of CO 2 shallow leakage

The calculation of the risk of surface leakage associated with each well supposed to intercept the Dogger aquifer is made by assigning a weighting coefficient to each criterion. The overall risk is obtained by multiplying them. There is no universal rule to estimate these values, but they are classically estimated through experience feedback. Without any further information, the values determined by [START_REF] Watson | Identification of Wells with high CO 2 -Leackage Potential in Mature Oil Fields Developed for CO 2 -Enhanced Oil Recovery[END_REF] and reported in Table 2 are used in the present study. The maximum value of the overall risk that may be associated to a well is 900, while the minimum value is 1. The mean overall risk through the 216 wells, taking into account the lack of information, is equal to 133.6, while the standard deviation is 60. Moreover we have seen that wellbore reports and database cannot always give us all the information we need to evaluate the value of all criteria. Then, when information is missing, we have decided to consider the two following extreme cases:

The MIN case: Whatever the criterion, we assign the minimum value (1) in the case of a lack of information. This is the case considered by [START_REF] Watson | Identification of Wells with high CO 2 -Leackage Potential in Mature Oil Fields Developed for CO 2 -Enhanced Oil Recovery[END_REF] or the MAX case : Whatever the criterion, we assign the maximum value associated to the criterion in the case of a lack of information

Risk levels

Using the values of table 2, a risk tree of 32 branches has been determined (five criteria, two possibilities for each criterion and therefore 2 5 = 32 branches, cf. Figure 1). Some branches are overlapped. The risk tree shows four risk levels:

Elevated risk: the CO 2 surface leakage risk of the well is greater than or equal to 200 and less than or equal to 900. In the risk tree, the wells being at this risk level are those with casing, abandoned before 2000 and without cementing to surface Mean risk: the CO 2 surface leakage risk of the well is greater than or equal to 40 and less than or equal to 180. In the risk tree, the wells being at this risk level are those with casing, abandoned before 2000 and with cementing to surface, those with casing, abandoned after 2000 and without cementing to surface, and those without casing, with spud date before 1980, abandoned before 2000 and without cementing to surface Low risk: the CO 2 surface leakage risk of the well is greater than or equal to 24 and less than or equal to 37.5. In the risk tree, the wells being at this risk level are those with casing, abandoned after 2000 and with cementing to surface, and those without casing, with spud date after 1980, abandoned before 2000 and without cementing to surface Very low risk: CO 2 surface leakage risk of the well must be greater than or equal to 1 and less than or equal to 22.5. In the risk tree, the wells being at this risk level are those without casing, abandoned before 2000 and with cementing to surface, and those without casing and abandoned after 2000

Considering the MIN and MAX cases, we obtain the distribution of risk levels shown in Table 3. In the MAX case where the lack of information maximizes the risk, we note that 56.02% of wells have a high risk of CO 2 surface leakage. This distribution is explained by the fact that 122 of the 216 wells (56.5% of wells) have a casing that may intercept the Dogger reservoir, 84% of wells have cement to surface and 95% have been abandoned before 2000 (because of the manner we have taken into account the lack of information) and these three criteria are very disadvantageous for wells that verify them.

In the MIN case, the manner to take into account the lack of information leads to minimize the risk (value of lacked criterion is 1), we notice only 26.39% of wells with a high risk of CO 2 surface leakage, while the number of well with a low risk rises from 6.48% to 28.70% of total number of wells. To have a spatial vision of the impact of lack of information, the results (MIN and MAX cases) are putted into a GIS. The resulting maps in Figure 2 emphasize that the classification is strongly influenced by the way that the lack of information is managed. Indeed, in the MIN cases, a CO 2 storage south-east of the Paris region appears to induce a risk of CO 2 leakage by wells from low to very low (only a few wells at high risk of CO 2 leakage should be taken into account), while in the MAX case, a medium to high risk would seem to prevail. These results clearly show the dependence of the results on how the lack of information has been taken into account. Thus, one must be careful in using the results of such a classification, rather than entering each well in one risk level, we recommend to attribute to each well the two risk levels associated with the "MIN" and "MAX" cases. So two maps would be generated for each CO 2 storage site, and it would be left to the decision makers, according to their knowledge of the studied Basin, to choose the most appropriate map.

Conclusion

This paper shows the application of a wellbore classification system within the Paris Basin, France. The input data necessary to establish the classification are mostly available. When data is missing, we have considered either: 1) a most prejudicial option assuming that the value of the criterion linked to this data is maximum; 2) or less prejudicial option where the criterion is not taken into account in the computation of the risk category (i.e., considered value is equal to 1). For the most prejudicial case, the results show that 55% of the wellbores that reached the selected deep saline aquifer (i.e. the Calcareous Dogger of the Paris Basin), are in the worst risk category. This high risk level can be explained by the fact that 56.5% of the selected wellbore have a casing that intercepts the aquifer, 84% are not cemented to the ground level and 95% were abandoned before 2000. On the contrary, if we consider the less prejudicial case, only 25% of the studied wellbore are in the worst risk category. This emphasizes that the management of the lack of data is of primary importance in this kind of study. Thus, without any additional information, the evaluation of the CO 2 leakage risk along the abandoned wells should at least exhibit the two extreme values that respectively correspond to the less and the more prejudicial options.

Unlike the Alberta basin, there are very few records of gas leakage from abandoned wells in the Paris basin and in France in general, making it difficult to derive a classification of abandoned wells according to their potential to leak. This is why we believe that using the adapted classification of Watson et al. [START_REF] Watson | Evaluation of the potential for gas and CO 2 leakage along wellbores[END_REF][START_REF] Watson | Identification of Wells with high CO 2 -Leackage Potential in Mature Oil Fields Developed for CO 2 -Enhanced Oil Recovery[END_REF] in the Paris Basin, might be very useful.

However the methodology used to classify the wells according to their risk level of leakage has its limitations since it does not take into account any changes in the integrity of the well in the long run. In future works, we therefore propose to investigate through numerical simulations the impact and thresholds of new criteria related to the environment around the well (chemical composition of fluid initially saturating the aquifer, the type of rocks intersected by the well, the in situ temperature and pressure, and stress regime).

Fig. 1 .

 1 Fig. 1. The risk tree used to derive the risk levels

Fig. 2 .

 2 Fig. 2. Mapping of results of the analysis of the risk of leakage through abandoned wells in the Seine-Normandy basin

Table 1 .

 1 Availability of data in the 216 abandoned wells likely to intercept the Dogger aquifer

	Criteria	Well type	Abandonment	Cement to	Well deviation Spud date
			date	surface

Table 2 .

 2 Weights assigned to criteria for the calculation of CO 2 surface leakage risk, taking into account the lack of information.

	According to Watson and Bachu [8]
	Criterion	Value
	Cased well	8
	Uncased well	1
	Lack of information	8 or 1
	Abandoned before 2000	5
	Abandoned after 2000	1
	Lack of information	5 or 1
	No cement to surface	5
	Cementing intervals unknown	4
	Other cases	1
	Deviated well	1.5
	Non deviated well	1
	Spud date before 1980	3
	Spud date after 1980	1

Table 3 .

 3 Comparison of distributions obtained on the risk levels in the two approaches

		Elevated risk	Mean risk	Low risk	Very low risk
	"MIN" case	26.39%	34.26%	10.65%	28.70%
	"MAX" case	56.02%	25.46%	12.04%	6.48%

Acknowledgements

This work is supported by the BRGM through the research program INTEGMECA10.

Availability

85.6% 57.87% 100% 100% 100%


Audun Bakk
email: audun.bakk@sintef.no

Jean-Francois Girard

Erik Lindeberg

Eyvind Aker

Frederic Wertz

Maike Buddensiek

Maria Barrio

David Jones

6th Trondheim CCS Conference (TCCS-6) CO 2 Field Lab at Svelvik ridge: Site suitability

Keywords: CO2, injection, monitoring, migration

A field laboratory for monitoring CO 2 migration has been established in the Holocene deposit of the Svelvik ridge, located at the sill of the Drammensfjord 50 km south of Oslo. Initial characterization of the site shows that the formation is very suitable for studying migration and associated leakage, particularly within the top 50 m where the sediments consist of relatively homogenous sand. In the deeper part of the deposit, the formation seems to have more structure possibly enabling monitoring of horizontal migration and more complex migration patterns.

Introduction

The CO 2 Field Lab project is a consortium of industry and research institutes that have had a strong involvement in CO 2 storage site management, monitoring and certification from the early stages. For both the planning and safety of any CO 2 storage operation, it is important to monitor the location and possible migration of CO 2 within and around the reservoir. The overall credibility of any storage project will ultimately rely on proper monitoring of CO 2 and understanding its behaviour.

The EU directives on geological storage of carbon dioxide [1] and on the emissions trading scheme [2] both require the ability to detect and quantify CO 2 within and outside the storage complex. However, specific instructions on monitoring are not yet developed. Therefore, the CO 2 Field Lab project plans controlled injection of CO 2 into the aquifer in the Svelvik ridge (Fig 1). The CO 2 displacement in the subsurface and at the surface will be monitored with an extensive set of techniques to test their sensitivity to the presence of CO 2 underground with respect to saturation, the size of the plume and the detection of leakage. Moreover, through assessing monitoring technologies for detection of CO 2 migration and ultimately leakage, the project aims to provide a validated monitoring system through a protocol and certification scheme. The Svelvik ridge is located about 50 km south-west of Oslo and forms the sill of the Drammensfjord. It is a glaciofluvial-glaciomarine terminal deposit formed during the Ski stage of the Holocene deglaciation [START_REF] Sørensen | Foreløpig beskrivelse til kvartaergeologisk kart SVELVIK -CL 083[END_REF][START_REF] Melø | Hydrogeology of the shallow aquifer at the Svelvik ridge[END_REF]. The depth to the bedrock is between 300 and 400 m. The central part of the ridge is subaerially exposed with the top about 70 m above sea level. It forms a phreatic aquifer. In this part of the ridge, sand has been excavated since 1915. Clay layers onlap both flanks of the ridge below sea level. To the south the thick clay/silt layer fills the bedrock basin up to a few meters below sea level, while to the north the thinner clay/silt layer is at water depths of 100 -120 m. This site was chosen as a field laboratory on the assumption that the sand ridge contains more or less homogeneous, unconsolidated, highly permeable sand, which offers well constrained conditions for controlled gas injection experiments. Naturally, there had to be an extensive site characterization to assess whether the conditions would enable to achieve the objectives.

The site characterization was conducted during 2009 and 2010. A geological model was built for use in flow simulations and geophysical modelling. A shallow (20 m deep) CO 2 injection experiment is planned for during the fall of 2011 and, we hope to subsequently carry out a deep CO 2 injection experiment (100 -300 m deep). This paper summarizes briefly the main results of the site characterization and discusses the potential of the site as a field laboratory for CO 2 monitoring.

Results and discussion

Site characterization

Beginning in November 2009 and throughout 2010, a series of surveys were conducted to characterize the site. This included drilling, sampling and logging of a 333 m deep exploration well, analysis of core and flow-line samples, geophysical surveys including resistivity, seismic reflection and ground penetrating radar along two 2D lines, and hydrodynamical, geochemical and soil gas surveys.

Due to the reverse circulation method used for drilling of the exploration well, the depth reference and quality of the flow-line samples were very satisfactory. Therefore, the visual inspection of these samples provides vital information about the sub-surface. The flow-line samples and cores show that the Svelvik ridge consists, as assumed, of permeable sand down to approximately 50 m. However, below 50 m results indicate layers or lenses of considerable vertical extent that consist of variable proportions of sand, silt and clay. For the shallow injection the feasibility and suitability of three different areas within the CO 2 Field Lab site experiment were surveyed with ground penetrating radar and surface gas measurements. Two areas were found to be problematical, due to either hard ground conditions which make gas sampling difficult, or due to the high water level (saturation at 25-65 cm depth). The third area did not have these problems and its suitability was confirmed by ground penetrating radar measurements. The surface gas measurements were used to assess baseline conditions. Soil gas concentrations and fluxes are low with CO 2 concentrations mostly below 1.0 % and fluxes less than 8 gm -2 d -1 . This is consistent with the scarcity of vegetation.

Modelling of the shallow CO 2 injection

The shallow injection will be conducted at 20 m depth where the formation is interpreted to be a relatively homogeneous sand body. The objective is to study migration and surface leakage under well controlled conditions.

Using the results from the core analyses, an Eclipse reservoir model was developed with an isotropic sand with a permeability of 2 Darcy equivalent to the average value of direct permeability measurements on cores taken in the 333 m deep exploration well [START_REF] Bøe | Petrographical and petrophysical analysis of a core and drill mud samples in a 330 m deep well in the Svelvik sand ridge[END_REF]. The modelling algorithm uses the equation of state of Span and Wagner [START_REF] Span | A new equation of state for carbon dioxide covering the fluid region from triple-point temperature to 1100 K at pressures up to 800 MPa[END_REF]. To model the pVT properties and the phase behaviour of the CO 2 /water system, the properties of CO 2 and the solubility model for CO 2 and water described in Spycher et al. [START_REF] Spycher | CO2-H2O Mixtures in the Geological Sequestration of CO2. II. Partitioning in Chloride Brines at 12-100°C and up to 600 bar[END_REF] were used. The relative permeability and capillary curves used in the Eclipse simulations are similar to experimental data from the Utsira sand [START_REF] Akervoll | Feasibility of reproduction of stored co2 from the Utsira formation at the Sleipner gas field[END_REF].

The modelling of various flow rates and durations, shows that an injection rate of 187.5 kg/d for 4 days would provide suitable conditions, in terms of plume growth, extension and migration, to perform a sensitivity study of the monitoring techniques planned for the shallow experiment.

The saturation profiles of CO 2 are shown in Fig. 3. Injection point is at 20 m depth. After approximately 55 hours injection, CO 2 arrived at the surface and a steady state occurred soon after. At this point, the escape rate can be adjusted simply by adjusting the injection rate. The maximum diameter of the exposed surface area is 7.5 m.

In this simulation 750 kg was injected. 45 % of the injected CO 2 (336 kg) escaped to the atmosphere, 46 % was dissolved and 9 % trapped as residual gas after 14 days. The escape profile is illustrated in Fig. 4. The peak escape rate of 7.5 kg/hour (lasting for about 1.4 days) corresponds to a flux of 4.1 kgm -2 d -1 . The escape rate falls to almost zero after the injection stops.

For the same shallow injection experiment, a different scenario including anisotropy was modelled using Tough2 simulations. In this simulation, the injection point was also 20 m below surface with a CO 2 injection rate of 173 kg/d. For this anisotropic model, the horizontal permeability was set to 20 Darcy while the vertical permeability was 2 Darcy. To obtain insight into the sensitivity of the modelling, the Sleipner [START_REF] Audigane | Two-dimensional reactive transport modelling of CO2 Injection in a Saline Aquifer at the Sleipner site, North Sea[END_REF] and Viking [START_REF] Bennion | Relative permeability Characteristics for supercritical CO2 displacing water in a variety of potential sequestration zones in the Western Canada Sedimentary Basin[END_REF] sandstone models were applied in addition to testing both models with and without capillary entry pressures (0 and 1500 Pascal, respectively). Fig. 3. The vertical saturation profile of the CO 2 plume in a shallow injection scenario simulated with Eclipse. The sand is assumed to be isotropic sand body (2 Darcy permeability).750 kg of CO 2 is injected totally The injection point is 20 m below surface and the total thickness of the sand body is 30 m. The colour bar indicates CO 2 saturation. The saturation varies form; dark blue = 0 % to red = 76 %. The total horizontal extension in the images is 30 m. The Tough2 simulation shows that injection of 200 -500 kg CO 2 is needed for it to reach the surface. The model with capillary entry pressure significantly increases the quantity of CO 2 needed for it to reach the surface and results in a more laterally extended plume (Fig. 5).

hours of injection 38 hours of injection 62 hours of injection

Injection stops 750 kg CO 2 injected 40 hours of injection 62 hours of injection 24 days after stop of injection

highly on the injection rate, intrinsic permeability, anisotropy factor, and relative permeability. This sphere contains the minimum required gaseous saturation that allows CO 2 mobility. As injection continues, the CO 2 sphere stops growing and starts moving upward, saturating the overlying porous medium with CO 2 , until it reaches the surface. Finally, almost immediately after breaking through the surface, the CO 2 flux attains a steady state, where the leakage rate is equal to the injection rate.

Conclusion

A field laboratory for monitoring CO 2 migration and leakage has been established at the Svelvik ridge (Norway). The characterization and associated modelling have shown that the site is very suitable for studying both migration and leakage at shallow depths. The shallow aquifer (down to approximately 50 m) consists of a relatively homogeneous sand body. Moreover, the site offers the opportunity to study migration signatures in more complex structures at semi-shallow depths (50-200 m) possibly with associated surface seepage. The deeper strata (below approximately 200 m) are most probably impermeable to CO 2 migration to the surface. However, these beds may be suitable for studying horizontal migration. By further confirmation of the site suitability, the established field laboratory may persist beyond the scope of the CO 2 Field Lab project.

Fig. 1 .

 1 Fig. 1. Aerial photo looking northwards to the Svelvik ridge at the outlet of the Drammensfjord. The sill is formed by deglaciation deposits. The photo shows the sand excavation on the ridge. The laboratory location is indicated by the white rectangle (extension is approximately 300 m × 150 m).

Fig. 2

 2 Fig. 2 shows the north-south vertical seismic profile. The depth migrated seismic section shows a number of reflectors, particularly below 100 m. Note, that reflectors in the deposits seem more continuous below 200 m.

Fig. 2 .

 2 Fig. 2. Seismic profile through the central part (oriented from north (left) to south (right)) of the Field Laboratory at the Svelvik ridge. The upper boundary of the bedrock is indicated by solid black lines.

Acknowledgements

The CO 2 Field Lab project is funded by Climit research programme (through Gassnova), DGCIS (French ministry), BGS, BRGM, Bureau Veritas, CNRS, imaGeau, Schlumberger Services Pétroliers and WesternGeco. Andy Barkwith, Bob Lister and Neil Breward from BGS are acknowledged for their contribution to surface gas measurements.