Damage spreading transition in glasses: a probe for the ruggedness of the configurational landscape

We consider damage spreading transitions in the framework of mode-coupling theory. This theory describes relaxation processes in glasses in the mean-field approximation which are known to be characterized by the presence of an exponentially large number of meta-stable states. For systems evolving under identical but arbitrarily correlated noises we demonstrate that there exists a critical temperature $T_0$ which separates two different dynamical regimes depending on whether damage spreads or not in the asymptotic long-time limit. This transition exists for generic noise correlations such that the zero damage solution is stable at high-temperatures being minimal for maximal noise correlations. Although this dynamical transition depends on the type of noise correlations we show that the asymptotic damage has the good properties of an dynamical order parameter such as: 1) Independence on the initial damage; 2) Independence on the class of initial condition and 3) Stability of the transition in the presence of asymmetric interactions which violate detailed balance. For maximally correlated noises we suggest that damage spreading occurs due to the presence of a divergent number of saddle points (as well as meta-stable states) in the thermodynamic limit consequence of the ruggedness of the free energy landscape which characterizes the glassy state. These results are then compared to extensive numerical simulations of a mean-field glass model (the Bernasconi model) with Monte Carlo heat-bath dynamics. The freedom of choosing arbitrary noise correlations for Langevin dynamics makes damage spreading a interesting tool to probe the ruggedness of the configurational landscape.

We consider damage spreading transitions in the framework of mode-coupling theory. This theory describes relaxation processes in glasses in the mean-field approximation which are known to be characterized by the presence of an exponentially large number of meta-stable states. For systems evolving under identical but arbitrarily correlated noises we demonstrate that there exists a critical temperature T0 which separates two different dynamical regimes depending on whether damage spreads or not in the asymptotic long-time limit. This transition exists for generic noise correlations such that the zero damage solution is stable at high-temperatures being minimal for maximal noise correlations. Although this dynamical transition depends on the type of noise correlations we show that the asymptotic damage has the good properties of an dynamical order parameter such as: 1) Independence on the initial damage; 2) Independence on the class of initial condition and 3) Stability of the transition in the presence of asymmetric interactions which violate detailed balance. For maximally correlated noises we suggest that damage spreading occurs due to the presence of a divergent number of saddle points (as well as meta-stable states) in the thermodynamic limit consequence of the ruggedness of the free energy landscape which characterizes the glassy state. These results are then compared to extensive numerical simulations of a mean-field glass model (the Bernasconi model) with Monte Carlo heat-bath dynamics. The freedom of choosing arbitrary noise correlations for Langevin dynamics makes damage spreading a interesting tool to probe the ruggedness of the configurational landscape. 05. 40.+j, 64.70.Pf, 75.50. Lk

I. INTRODUCTION
The theoretical understanding of the dynamical behavior of glasses is a long outstanding problem in statistical physics which has recently revealed new aspects on the underlying mechanism responsible of the glass transition 1,2,4 . The dynamical behavior of glasses is characterized by the fast grow of the characteristic time of relaxation processes in the vicinity of the glass temperature T g . This increase of the relaxation time, up to fifteen orders of magnitude in a relatively small range of temperatures, is usually termed as the viscosity anomaly. The first consideration of this anomaly, the Vogel-Tamman-Fulcher law, goes back to the twenties. However, there is still no satisfying and generally accepted theoretical explanation for this singular behavior.
Currently, there are two main approaches to understand the glass transition problem. One approach (the Adam-Gibbs-DiMarzio theory 6 ) focuses on thermodynamic considerations and proposes the existence of the ideal glass transition T g . This is a singularity where the configurational entropy of the under-cooled liquid vanishes at T g and a second order phase transition characterized by a finite jump in the specific heat occurs. This scenario has been rediscovered in the framework of mean-field spin-glasses with one step of replica symmetry breaking 7 . In the meanfield approach, T g is the temperature where configurational entropy vanishes (the so called Kauzmann temperature) and also replica symmetry breaks. In this paper we will denote the glass transition both by T g and T s (in this last case the subindex s stands for statics). The other approach relies on mode-coupling theory and describes the glass transition as a strongly non-linear dynamical effect which induces long-term memory properties in the correlation and response functions 3 . Consequence of these effects is the existence of a dynamical singularity T d where ergodicity breaks and correlation functions do not decay to zero. This dynamical transition T d is a consequence of the mean-field character of the mode-coupling approximation. Although these two approaches are apparently different they have in common their mean-field character.
One of the most distinct features of glasses is the presence of a complex free energy landscape. The viscosity anomaly is a signature of activated dynamics due to the existence of a rugged free energy landscape with several maxima and minima separated by energy barriers and saddle points which connect them. One could think that the existence of this type of landscape is a necessary ingredient to find the previous scenario. Well known results on the number of meta-stable states in spin-glasses reveal that the interesting spin-glass behavior emerges in systems with an exponential number of states 5 . For instance, models such as the spherical p-spin interactions spin glass (with p > 2), the Ising p-spin interactions spin glass (with p > 3), the Edwards-Anderson model and Potts glasses are characterized by an exponentially large number of meta-stable states. All these models are characterized by the presence of quenched disorder which facilitates analytical treatments always at the level of mean-field theory. In the absence of disorder, similar results are found [43][44][45] although exact calculations for the number of meta-stable states turn out to be more difficult. It is largely believed that an exponentially number of meta-stable states is a necessary condition for the existence of replica symmetry breaking.
There have been also recent studies of exactly solvable models which, in the absence of quenched disorder, also exhibit glassy behavior 8,9 (and, in particular, activated behavior of the relaxation time). These models are characterized by a small number of meta-stable states. What makes them to display glassy behavior is the presence of entropy barriers which leads to slow dynamics even in the absence of meta-stability. Consequently, one is tempted to conclude that a rugged free energy landscape with a large number of meta-stable states is not essential to find glassy behavior but the presence of an enormous number of flat directions in phase space.
A similar conclusion was reached in the study of the Sherrington-Kirkpatrick spherical spin-glass 10 . It was found that an enormous number of zero modes are responsible for the slow dynamics found in this model. Although that model does not have an exponential large number of meta-stable states and does not show activated dynamics it displays glassy behavior at low temperatures due to the existence of flat directions around the meta-stable states 11 .
A possible way to investigate the existence of flat directions in a rugged free energy landscape is the study of the damage spreading. Damage spreading consists in the study of the dynamical evolution of the distance D(t) (to be defined later) between two systems configurations evolving under the same dynamical rules and differing only in their initial condition (the so called initial damage).
The study of the damage spreading problem (hereafter referred as DS) was proposed by Kaufmann in the 60's for the study of the propagation of mutations in the genotype in the biological growth of individuals 12 . That is, how a small perturbation in the genotype (microscopic level) manifests itself in the long-time term in the fenotype (macroscopic level). Since then, such a problem has received considerable attention in the framework of statistical physics, particularly in the middle 80's 13 . Almost fifteen years ago it was realized that DS could be a powerful tool to distinguish different dynamical regimes in disordered systems, such as spin glasses 39,15 . Variants of damage spreading phenomena have also been proposed to numerically investigate equilibrium correlation functions in generic statistical systems 16 and lattice gauge theories 17,18 . However, the initial enthusiasm and exciting perspectives in the research of this problem decreased in subsequent years after realizing that this transition was dependent on the type of dynamics used. So the existence of the DS transition could have nothing to do with the presence of a thermodynamic phase transition. Physicists then started to systematically investigate the DS in well known ordered systems such as the Ising model. In particular, much work has been devoted to the study of the one-dimensional Domany-Kinzel automaton 19 and the one-dimensional Ising model 20,21 . The question of the non-universality of damage spreading has been also emphasized in the context of non-equilibrium phenomena such as domain growth by Graham, Hernandez-Garcia and Grant 33 .
So the question remains whether this transition has a truly physical meaning or not. In this direction, Hinrichsen and Domany tried to give a precise dynamic-independent definition for DS. To define a damage spreading phase one must consider all possible dynamical procedures which lead the system to thermal equilibrium. For discrete systems with small number of nearest neighbors this definition can be implemented but not in the general case (for instance, continuous systems) where an infinity of dynamical rules can be always implemented.
The purpose of the present work is to present a detailed study of the DS in a model with a rugged free energy landscape with an exponentially large number of meta-stable states. In particular, we will study the DS in the p-spin spherical spin glass, an exactly solvable model for the glass transition which is described by the Adam-Gibbs-DiMarzio scenario and whose dynamics is described by the mode-coupling equations. To be more specific, we will study damage spreading for Langevin dynamics, the simplest dynamics which is continuous in time and satisfies ergodicity and detailed balance. It must be stressed that although there are very few works on the DS problem using this dynamics (Stariolo 32 and Graham et al. 33 ), the majority of theoretical works in DS have considered discrete dynamics (in discrete systems).
We suggest that damage spreading can be used as a dynamical method to evidenciate the existence of a large number of flat directions or saddle points in phase space. Also the existence of a dynamical transition will be shown. Although we will check that damage spreading transitions are strictly non-universal still it is possible to use the asymptotic distance to define an order parameter for this dynamical transition. We also anticipate that by considering correlations between the noises of the two evolving systems an infinity of dynamical transitions can be obtained. For Langevin dynamics, the case of maximally correlated noises has a particularly interesting physical meaning.
Because our study considers the DS in the mode-coupling theory for glasses it is expected to be generally valid for Langevin dynamics in systems with a rugged free energy landscape such as realistic glasses. Although the DS transition is non-universal and depends on different dynamical rules (or cross-correlations between the stochastic noises), we believe that this transition gives an interesting information on the free energy landscape and could be investigated in structural glasses. Being a signature of the existence of saddle-points in phase space (i.e. points which separate stable and unstable phase space directions) it is natural to expect that real glasses are good systems to manifest these effects. This consideration makes our results more attractive from the viewpoint of numerical simulations of realistic glasses 23 . Nevertheless we call the attention to the reader that some of our claims in this paper are not generally proven (such as the connections between damage spreading transitions and saddle points in phase space) and the present reasarch should be seen as a first step for a better understanding of some of these questions.
The contents of the paper are divided as follows. The second section is devoted to general considerations and definitions about the DS problem. Section III demonstrates the existence of a dynamical transition T 0 for DS. Section III is divided in three subsections. Subsection III.A describes the mode-coupling equations for the p-spin spherical spin glass, starting from a random initial configuration for different correlations between the noises. This subsection also describes the different numerical methods we have used to analyze the mode-coupling equations for different cases. Subsection III.B analyzes the DS problem starting from an equilibrium initial condition. Subsection III.C analyzes the DS in the presence of asymmetric interactions. Section IV presents an analysis of the damage spreading for glassy models with discrete dynamics. Concretely, we study the DS transition in the Bernasconi model with heat-bath dynamics. Finally section V presents the conclusions. Two appendices are devoted to some technical issues.

II. GENERAL CONSIDERATIONS ABOUT DS
In the most general framework the DS problem can be stated as follows. Consider a dynamical system described by a generic variable C which denotes a given configuration. Suppose that the system evolves under a deterministic dynamical rule F . For sake of simplicity we take a continuous time dynamics. The equation of motion reads, On top of the configuration C and the dynamical rule F we also need to define a distance in the phase space of configurations (for instance, a Hamming distance for spin systems). This distance D needs to satisfy the usual good properties, in particular D(C t , C t ) = 0 at all times. Suppose that we take two initial configurations C 0 , C ′ 0 with initial distance D 0 = D(C 0 , C ′ 0 ) and consider the generic equal times distance, where C(t) and C ′ (t) start from configurations C 0 and C ′ 0 at time 0 and evolve under the same dynamical rule F eq.(1). Our main interest is to investigate the value of the asymptotic long-time distance D ∞ , Note that D(t) = 0 if D 0 = 0. Quite generally the asymptotic distance D ∞ will be a function of the type of initial configurations C 0 , C ′ 0 as well as their initial distance D 0 . The dependence of D ∞ on those parameters is governed by the dynamical properties of the deterministic rule, such as chaotic properties and Lyapunov exponents.
One could extend this general problem to stochastic systems, i.e. dynamical systems which evolve in the presence of a stochastic noise. Let us consider two systems described by the configuration variables C t , C ′ t at time t which evolve following a Langevin dynamics,Ċ following eq.(4) with the same realization of the stochastic noise and starting from two different initial conditions. We are thinking of noises statistically identical which coincide at equal times (i. where K is a generic function which satisfies the properties K(C, C) = 1 and −1 ≤ K(C, C ′ ) ≤ 1, ∀C, C ′ . In the presence of stochastic noise, again D(t) = 0 if D(0) = 0. Note that the role of the correlations is irrelevant for the evolution of the independent systems C, C ′ but crucial for their correlations and the equal times distance (2). Different choices of the cross-correlation K for Langevin dynamics are the analogous of different dynamical rules in discrete dynamics such as Monte Carlo (these rules could be Glauber, Metropolis or heat-bath among others). Now we are interested in the asymptotic long-time distance eq.(3). Quite generally, D ∞ will be a function of the type of initial condition (for instance, random or stationary), the initial distance D 0 , the intensity of the noise T , and the cross-correlator K. For the case in which η t = η ′ t (K = 1), we will find that there is a dynamical phase transition at a finite temperature T 0 below which the asymptotic distance is different from zero. The origin of this dynamical transition can be explained in plain physical words. In eq.(4) there is competition between two different terms. On the one hand, the force term F (C t ) propagates the error (or damage) in the initial configuration. Instead, the noise η t acts in the same way in both systems smearing out possible differences in the initial condition. In other words, the stochastic noise is the synchronizing force which tries to make both evolving configurations to merge in time while the force term amplifies the initial damage playing the role of a noise. This argument only applies if K = 1. In the general case −1 ≤ K ≤ 1 the noise does not necessarily synchronize both systems and its effect is similar to that of the force. In this case, the asymptotic distance will be also a function of the cross-correlation K (i.e. the dynamical rule).
To understand better the role of the cross correlation K let us consider as a starting point the simple problem of a particle which moves inside an harmonic potential V (x) = 1 2 x 2 following a Langevin dynamics, where η(t) is a stochastic white noise of variance 2T . The configuration C corresponds to the position x of the particle and we define a distance between two configurations x, y as D(x, y) = (x − y) 2 . Take now two identical particles x, y and make them follow eq.(5) both with the same stochastic noises η, η ′ and cross-correlation K(x, y). For simplicity we will take a x, y symmetric cross-correlation K = K(D). If < .. > stands for average over dynamical histories then the distance D(t) =< (x(t) − y(t)) 2 > satisfies the following equation, This equation has several stationary solutions depending on K. Obviously D = 0 is a stationary solution (remember K(0) = 1) which implies K ′ (0) < 0. Only if −1/2T < K ′ (0) < 0 the solution D = 0 is stable. On the other hand, for K ′ (0) < −1/2T the solution D = 0 is unstable and D(t) converges to another stationary solution (which can be shown it is always stable). So, there is a dynamical transition at T c = 1/(2|K ′ (0)|) where the asymptotic distance D ∞ changes from zero (T < T c ) to D ∞ = D * where D * satisfies the identity D * = 2T (1 − K(D * )). The asymptotic distance is then given by, D ∞ = lim t→∞ D(t) = D * Θ(T * − T ) independent of the value of the initial distance D 0 between the two particles. So when K ′ (0) does not vanish, already for the simple harmonic oscillator there is more than one stationary solution. In other words, the effect of the cross-correlation term K(D) manifests through the appearance of more than one stationary solution. As the reader can imagine the discussion turns out to be more difficult for other more complex potentials.
Although this dependence of the asymptotic distance D ∞ on the cross-correlation K is an intrinsic property of damage spreading it does not necessarily imply that this kind of dynamical phase transitions do not give any relevant information on the physical properties of the system. What this really means is that the results concerning DS may depend on particular forms of cross-correlation function between the noises (similarly as happens for discrete dynamics where different dynamical rules yield different results). Nevertheless, generic results for DS may be obtained for correlators which satisfy quite general conditions (as we will see later).
The major part of work in this paper will be devoted to the study of maximal cross-correlations, i.e K = 1 although the results are extrapolable to more generic cross-correlators. For Langevin dynamics, the case K = 1 is particularly appealing for two reasons.
On the one hand, it followed from the simple example of the harmonic oscillator that K = 1 implies K ′ = 0 everywhere. Below we will argue that this observation holds also for more complex potentials. More specific, we argue that in the case K = 1 there is a single stationary solution D ∞ = 0 for any finite system and for any confining potential which diverges in the boundaries (i.e V (x) → ∞ when x → ±∞). Clearly, the harmonic oscillator is a trivial case where the asymptotic distance always goes to zero independently of the temperature T and how far is x from y at t = 0. In order to justify our assertion let us take a more complex potential of two wells separated by a finite barrier, for instance a particle moving inside an harmonic plus a quartic term potential V (x) = − 1 2 x 2 + λ 4 x 4 . In this case the potential has two wells located at x = ± 1 √ λ . If two systems described by the variables x, y start to evolve within the same well (i.e. x, y > 0 or x, y < 0) they will always tend to finish in the same final configuration because the synchronizing effect of the noise is, at very long times, the dominant effect. If they start in different wells the conclusion is also the same because there is always a finite probability that a strong fluctuation in the stochastic noise drives both particles in the same well. This conclusion, which holds for maximal cross-correlations K = 1, can be generalized for any potential with a finite number of wells separated by finite energy barriers. There is always a finite probability that a fluctuation of the noise can take both particles into the same well and hereafter the distance D(t) between both configurations would tend exponentially fast to zero. Obviously this argument applies only for finite barriers, finite wells as well as finite temperature. At zero temperature the synchronizing effect of the stochastic noise is absent and the asymptotic distance may not vanish and show a non-trivial dependence on the initial distance.
From the discussion above it follows that the role of the thermodynamic limit N → ∞ in the DS is also crucial. In this limit the height of the barriers or the number of wells (i.e. meta-stable states) may diverge. The first case happens, for instance, in the ferromagnetic Ising model where the time-reversal symmetry of the Hamiltonian is broken below T c . The second case is realized in spin-glass models where the number of meta-stable states is exponentially large with N . In both cases, a finite fluctuation of the stochastic noise (even with K = 1) may not synchronize the system and D ∞ can be a non trivial function of both the temperature and the initial distance. To be more precise D ∞ is defined as follows, Note that for N finite we expect lim t→∞ D(t) = 0 so it is crucial that the thermodynamic limit N → ∞ is taken before the infinite-time limit. Taking the limits in the reverse order will result in that D ∞ will always vanish at finite temperature. Note that this discussion applies only when the D = 0 stationary solution is stable. This is indeed satisfied for K = 1 but may be also fulfilled in more general situations with K < 1 and D = 0 still being an stationary stable solution.
There is another property of the case K = 1 which makes it particularly interesting. Up to now our discussion was limited to different mathematical properties of damage spreading transitions. But, what about their physical significance? Suppose we take two generic statistical systems described by the set of variables {x i , y i }, i = 1, .., N which evolve under the Langevin dynamics, where, as before, η, η ′ are white noises with cross-correlation K which we will suppose is a generic function of the Hamming distance D. Let us suppose that the force derives from a potential F i ({x}) = − ∂V ({x}) ∂xi . If we define the new variables z i = x i − y i we may obtain, subtracting both equations in (8), where ν is a stochastic white noise of zero mean and variance 4T (1 − K(x, y)). The solution z i = 0 is a stationary solution of (9). A linear stability analysis around that solution yields the equation, where ∂yi∂yj is the Hessian matrix evaluated at the point (y 1 , .., y N ). The solution z i = 0, corresponding to vanishing Hamming distance D = 0, is stable if the Hessian is negative definite. The presence of the stochastic noise ν in (10) decreases the stability of the D = 0 solution. Because K ≤ 1 we conclude that the region where the D = 0 stationary solution is maximally stable corresponds to the case when K = 1 because ν vanishes. So K = 1 is the cross-correlation for which the stationary solution D = 0 is maximally stable. For K = 1 equation (10) is quite appealing and evidenciates the physical origin of the DS transition. An instability of the D = 0 solution may appear when an eigenvalue of the Hessian matrix H ij vanishes. This corresponds to a saddle point of the potential landscape V ({y}). Due to ergodicity, the x, y systems sample all the possible configurations. So the asymptotic distance D ∞ is a direct measure of the stability of the D = 0 solution along all possible configurations (weighted with their corresponding statistical Boltzmann weight). In other words, a instability in the solution D ∞ = 0 and the existence of a DS transition is indication of the presence of saddle points in the potential landscape of systems x, y. If K < 1 the temperature of the DS transition will depend on the particular form of K (actually we will see later, in the study of mode-coupling equations, that it depends through the value K ′ (D = 0)). Furthermore, due to the destabilizing effect of the noise the damage spreading transition will increase when K decreases so K = 1 yields the lowest damage spreading transition temperature among all possible cross-correlations K(D) for which the D = 0 solution is stable. In the presence of fixed points for the asymptotic distance other than D = 0 the stationary solution D ∞ = 0 may become unstable because the noise ν is too strong (similarly to what happens in the harmonic oscillator example) and damage spreading does not evidenciate anymore the existence of saddle points. In other words, saddle points may be observed only by studying the T dependence of the basin of attraction of the Q = 0 stationary solution supposing it stable at very high-temperatures (free case). Remark that a similar argument has already been presented by Loreto, Serva and Vulpiani 24 for systems described by a single variable x(t) in a potential field V (x).
We have argued above that for Langevin dynamics the maximal cross-correlator K = 1 is a special case, resulting in a simplification of the problem. We emphasize that in other dynamical systems it is unclear whether or not the maximal cross-correlator K = 1 plays the same role in the context of damage spreading. This is due to the complexity of cross-correlations.
In the next section, we will analyze in detail the DS dynamics in the mode-coupling theory of glasses. As has been previously said these equations describe the relaxation processes and dynamics in glasses (in the under-cooled liquid regime) in the mean-field approximation and represent the dynamical behavior of systems with an exponentially large number of meta-stable states.

III. DAMAGE SPREADING IN MODE-COUPLING THEORY
Mode-coupling theory describes relaxational processes in glasses. In the nutshell, mode-coupling theory corresponds to an exact re-summation of an infinite series of diagrams in the hydrodynamic theories. The kind of diagrams which the mode-coupling approximation selects are those which precisely survive in the mean-field limit of some realistic models. So, a way to get mode-coupling equations is by considering exact dynamical theories for mean-field disordered spin-glass models 7 . Spherical spins allow for an exact closure of the dynamical equations in terms of correlation and response functions (as was shown by Crisanti, Horner and Sommers 27 in the p-interactions spherical spin glass), leading to many analytical results. Although spherical spins are unrealistic (compared to Ising spins) they capture the essentials of the dynamics which is universally found in a large variety of models. Whereas for p = 2 the physical description of the model is quite simple 10 the behavior turns out to be much more interesting for p > 2 where an exponentially large number of meta-stable states are present 28 . In this type of models meta-stability plays a very important role so, according to the arguments of the previous section, we expect to get interesting results for the DS transition. Forthcoming subsections analyze this transition in detail. A good review on the main results obtained in this model have been collected and reported by Barrat 30 .

A. Random initial configuration
This section is devoted to the study of the DS problem in the mode-coupling equations. This section describes some preliminary work already presented in 25 but here we present a more extended research of the problem to include asymmetry, different class of initial conditions as well as general cross-correlations of the noises. The simplest solvable model whose dynamics is described by the off-equilibrium mode-coupling equations is the spherical p-spin glass model introduced by Crisanti, Horner and Sommers 27 . In this case, the configurations are described by N continuous spin variables The Langevin dynamics of the model is given by, where F i is the force acting on the spin σ i due to the interaction with the rest of the spins, The term µ in eq.(12) is a Lagrange multiplier which ensures that the spherical constraint is satisfied at all times and the noise η satisfies the fluctuation-dissipation relation η i (t)η j (s) = 2T δ(t − s)δ ij where ... denotes the noise average. The J i2,i3,..,ip i are quenched random variables with zero mean and variance p!/(2N p−1 ). The interactions J i2,i3,..,ip i are symmetric under the interchange of the super-indices i 2 , i 3 , .., i p but in the most general case may not be symmetric under the exchange of the subindex i with a generic super-index. So, for instance J i2,i3,..,ip i Most of the studies undertaken in this model concentrate on the symmetric case where J i2,i3,..,ip i is symmetric under the permutation of all possible indices. This case is particularly interesting because there exists an energy function such that the force F i derives from a Hamiltonian or potential function F i = − ∂H ∂σi so there exists an stationary state described by a Boltzmann-Gibbs distribution. Due to the mean-field character of the model, the dynamical equations depend on the statistical properties of the force only through its correlations. On the other hand, the statistical properties of the force F i depend on the correlations of the J ′ s. The simplest case 29 corresponds to correlations of the type, for every k. So if α = 1 we recover the symmetric case while for α = 0 we obtain the asymmetric case. Eq.(13) implies the following statistical properties for the force F i 29 , where f (q) = q p /2. In the asymmetric case α = 0 the forces are completely uncorrelated at different sites. Hence the set of equations (11) become uncorrelated and the problem can be partially solved. This particular case will be analyzed later. For α < 1 there does not exist an energy function H which drives the system to thermal equilibrium and the fluctuation-dissipation theorem is not fulfilled.
We define the overlap between two configurations of the spins σ, τ by the relation Q = 1 N N i=1 σ i τ i so the Hamming distance between these two configurations is, in such a way that identical configurations have zero distance and opposite configurations have maximal distance D = 1. Then we consider two copies of the system {σ i , τ i } which evolve under the same statistical noise (11) with cross-correlation K but with different initial conditions. We assume the cross-correlator to be a function of the Hamming distance D or the overlap Q. As a consequence, other choices for the cross-correlator (in general this could depend on both configurations C, C ′ ) may change our results obtained below for K(Q) = 1. The major part, however, is concerned with K = 1 (we will explain why) and will not suffer from this. In this section we restrict our attention to random initial configurations (i.e equilibrium configurations at infinite temperature) with initial overlap Q(0). The case of initial equilibrium configurations will be analyzed in the next subsection. The different set of correlation functions which describe the dynamics of the system are given by, where < .. > denotes the average over dynamical histories and h σ i , h τ i are fields coupled to the spins σ i , τ i respectively. In what follows we take the convention t > s. The previous correlation functions satisfy the boundary conditions, C(t, t) = 1, R(s, t) = 0, lim t→(s) + R(t, s) = 1 while the two replica overlap Q(t, s) defines the equal time overlap Q d (t) = Q(t, t) which yields the Hamming distance at equal times or damage D(t) through the relation (15). Following standard functional methods 31,30 it is possible to write a closed set of equations for the previous correlation functions. Some details of the computation are shown in Appendix A. The final result is: while the Lagrange multiplier µ(t) and the diagonal correlation function Q d (t) obey the equations, Note that the cross-correlation K(Q d ) only enters explicitly through the equation (23) so it does not affect the evolution of a single replica. The whole set of equations is quite involved. For the correlation C and response functions R eqs. (19,20,22) several results are known, in particular their behavior in the equilibrium regime (where time-translational invariance is satisfied and the fluctuation-dissipation theorem is obeyed) as well as in the nonstationary aging regime 31 .
In what follows we analyze different dynamical fixed points of eq.(23) and show the existence of a dynamical instability in the DS equations.
B. Fixed-point analysis for a generic cross-correlation K Different type of dynamical regimes may be distinguished depending on the cross-correlator K. Our analysis is similar to that performed in section II for the simple harmonic oscillator. Different fixed points for the dynamical equations can be analyzed from equation (23). If the temperature T is very large then eq.(23) becomes, where we have used µ = T using eq. (22). Equation (24) can be exactly solved. The stationary solutions are given by K(Q) = Q. In figure (1) we analyze the different solutions for a generic K. We find that there are different stationary solutions corresponding to all possible intersections of the two curves (Q and K(Q)). A linear stability analysis of (24) reveals that stationary solutions Q * are stable if K ′ (Q * ) < 1. A dynamical flow diagram can be constructed where the region of stability is indicated by different arrows. Regions where K(Q) > Q satisfy ∂Q d ∂t > 0 and regions where K(Q) < Q satisfy ∂Q d ∂t < 0. So in this case one may depict a diagram of all possible high-temperature dynamical phases which separate regions with different fixed-point attractors. Damage spreading transitions will strongly depend on the type of cross-correlator. The case K = 1 is shown in figure (2) where there is a unique attractor at Q = 1 at very high-temperatures. This analysis of the different dynamical phases is valid only at very high temperatures. As soon as the temperature is finite and starts to decrease some of the stable fixed points may become unstable and other unstable points may become stable. The damage spreading transition correspond to the appearance of an instability in one of these high-temperature fixed points. As we will see below, the damage spreading transition temperature may be different for different fixed-points since it depends on the value of K ′ (Q * ) which may vary for different fixed-points Q * . In what follows most of our discussion will concentrate on the particularly interesting case K = 1 which has a unique fixed point at Q * = 1. Although the analysis may be extended to other fixed-points this case is also the most interesting according to our preceding discussion in section II. As we will check below, this case also defines the lowest damage spreading temperature T 0 among all the possible cross-correlators K for which the fixed-point Q * = 1 is stable. A first glance to equations (21), (23) reveals that the overlap Q(t, s) and its diagonal part Q d (t) are coupled each other through the correlation C(t, s) and response function R(t, s). The trivial solution Q(t, s) = C(t, s) and Q d (t) = 1 corresponds to the case where the initial conditions are the same, Q d (0) = 1 and the distance D(t) = 0 for all times. This high-temperature fixed-point (hereafter we will denote it by HT) corresponds to D ∞ = 0 and is asymptotically reached by the dynamics for high enough temperatures under certain conditions of the cross-correlator K (see the preceding discussion). In what follows we concentrate our attention in the case α = 1 where there is a stationary solution for C, R and µ corresponding to the equilibrium measure. Numerical integration of those equations (see later) reveals that the typical time needed to reach that solution grows if temperature decreases. At a given temperature (which we identify with T 0 ) there is an instability in the dynamical equations (21),(23) and the asymptotic solution differs from the HT one. We did not succeed in finding an explicit expression for T 0 but we have been able to show its existence and find a lower and upper bound for its value (see Appendix B): Note that for the particular case p = 2 and K ′ (1) = 0 the inequality (25) yields T 0 ≤ 1. Taking into account that (25) was derived under the assumption T 0 ≥ T s = 1 (i.e. we supposed we were in the high-temperature regime) this is not inconsistent with the result T 0 = 1 derived by Stariolo 32 .
Note that both the lower and upper bound for T 0 diverge when K ′ (1) = 1. This limiting value sets a condition on the possible cross-correlations K(Q) where the Q = 1 solution is stable. Only those functions K(Q) for which K ′ (1) < 1 are those for which Q = 1 is linearly stable at very high-temperatures. According to our discussion in section II the appearance of DS in this case is related to the presence of a divergent number of saddle points which mark the onset of a dynamical instability. Because K ′ (1) < 0 then (according to (25)) the limiting case K ′ (1) = 0 (for instance, if K(Q) = 1) sets the lowest value of the temperature T 0 where DS appears. This is important because it means, that whatever correlator K we consider (such that the solution Q = 1 is stable for high-enough temperatures) damage spreads below T 0 = p−2 2 . Remark that in the general case p ≥ 3 the dynamical instability temperature T 0 stays well above any relevant critical temperature (T s or T d ).
In the next section we discuss the behavior of the asymptotic distance as a function of temperature. For simplicity, our analysis is restricted to the case K = 1 for which most of the numerical work has been done. We will see that D ∞ , for a given specification of the correlator K, seems indeed to play the role of a dynamical order parameter in DS transitions.

Numerical analysis
In general it is too complicated to obtain an analytical solution of the set of equations (19)- (23). We shall devote this section to a numerical study of the equations (19)-(23) for the DS problem. Although in some particular cases an exact solution can be found (see below) this is not the general situation.
First, one could investigate the long-time limit of Q d via a numerical integration of the set of equations (19)- (23). However, the CPU time and the memory needed to do this grows very fast with time because of the integrals occurring in the equations. Thus the spreading of damage at large times can only be obtained from the dynamical equations doing some extrapolations. This enlarges the error in the estimate of D ∞ , especially in cases where Q d (t) is a non-monotonic function of time. In figure (3) we show how the overlap Q d (t) depends on the initial condition. Although direct extrapolations from numerical data of the value of the asymptotic damage are difficult, the figure is not incompatible with an independence of D ∞ on the initial condition. Another more powerful technique is necessary to corroborate this result. An alternative approach to obtain the long-time behavior of time dependent variables with high accuracy was introduced by Franz, Marinari and Parisi 34 to study the long-time behavior of the energy. Here we extend their method to analyze the asymptotic behavior of D(t). In their method they first decompose the time dependent variables in a series expansion before extrapolating for large times with the help of Pade approximants. For the DS problem, it leads to a Taylor expansion of the correlation function C, the response function R and the overlap Q: where c 00 = r 00 = 1 since C(t, t) = 1 and lim t→(s) + R(t, s) = 1 35 . Moreover, µ(t) and the diagonal correlation function Q d (t) can be written as: with µ 0 = T and q 00 is a parameter identical to the value of Q d at t = 0. Assuming always t > s, the dynamical equations (19)-(23) can be transformed into recurrence relations for the coefficients of the expansion. To this end one first substitutes (26) and (27) into (19)- (23) under the constraint t > s, then calculate the integrals and finally rearrange terms.
Numerically, the coefficients of the expansion are now readily obtained. In case p = 3, the first 80 coefficients of the expansions can be computed on a RISC workstation in a few hours. However, for larger values of p the computational effort is larger.
To ensure high accuracy of the asymptotic extrapolation, one needs a large radius of convergence of the series expansion. A good method to enlarge the radius of convergence of a series expansion is using Pade approximants. In this method one introduces two polynomials U m (t) and V k (t) of degree at most m and k, respectively. The goal is to choose U m (t) and V k (t) for given m and k such that Q d (t) and U m (t)/V k (t) are equal at t = 0 and have as many equal derivatives as possible at t = 0.
The computations have been performed for the symmetric case α = 1 and K = 1. Moreover, three different initial conditions have been considered: a) anti-correlated random initial conditions with Q d (0) = −1, b) uncorrelated random initial conditions with Q d (0) = 0 and c) partially correlated random initial conditions with Q d (0) = 0.5.
To check that the extrapolations D ∞ using the Pade approximants are correct, the Pade series have been compared with numerical integrations of the dynamical equations. Indeed, the Pade series and the numerical integration fit closely 25 .
The estimate for D ∞ is obtained by division of the highest order coefficients of P m (t) and Q k (t), i.e., by a m /b k . Moreover, an asymptotic estimate can be obtained assuming a power law decay of the equal time overlap: Q d (t) = Q d (∞) + At −γ . The analysis of D ∞ suffers in some cases from a small radius of convergence (even with Pade) as well as from the presence of poles in the Pade expansion. The results are displayed for p = 3 in figure 4 and for p = 4 in figure 5 for cases a), b) and c) as a function of the temperature. Let us remark that a lower number of coefficients in the Taylor expansion in the case p = 4 with respect to p = 3 leads to a less accurate estimate of the asymptotic distance. Inspection of figures 4 and 5 reveals that the dynamical transition T 0 is in the predicted regime eq. (25). It can be estimated more accurately from the relaxation time τ relax associated to the decay of the distance D(t) to zero. Starting from high temperatures, we assume the relaxation time to diverge at T 0 according to a power law: τ relax ≃ (T − T 0 ) −γ . We thus have found T 0 (p = 3) = 1.04 ± 0.02 with γ ≃ 1.1 ± 0.1 and T 0 (p = 4) = 1.13 ± 0.02 with γ ≃ 1.1 ± 0.1.
We conclude that for all temperatures, both p = 3 and p = 4, the asymptotic distance is independent of the initial distance. This is in contrast with the case p = 2, where dependence of the initial distance is found for the low temperature region 32 . We must say that we have obtained the same results, as in cases p = 3, 4, for a model which is a combination of the p = 2 and p = 4 spherical spin-glass model 40 . For a certain range of parameters, this model is known to have a continuous phase transition with continuous replica symmetry breaking and without collapse of the configurational entropy. So the first-order character of the spin-glass transition found in our model for p > 3 is not essential for the appearance of the DS transition. Still, that model 40 , is also characterized by the presence of an exponentially large number of meta-stable states. From the point of view of the form of the dynamical equations, the fact that T 0 is present for p > 2 as well as in a model which combines p = 2 and p = 4 is a consequence of the non-linearity in the coupling between the damage Q d (t) and the two-time correlation function Q(t, s) which occurs for all p > 2. From the physical point of view, this independence of D ∞ on the initial distance is quite appealing. In general, one would expect that p = 4 is quite similar to p = 2 due to the presence of the time-reversal symmetry. The fact that the damage does not feel this symmetry for p > 2 means that the separation of dynamic trajectories does not occur in the borders or maxima which separate equilibrium states but within saddle points of the phase space which divide configurations separated by finite energy barriers. This is supported by the fact that the transition occurs at a temperature much higher than T s and, as we will see in the following section, from the fact that it happens starting already from an equilibrium configuration. The asymptotic value D ∞ can, on the basis of our computations, be regarded as an order parameter for the transition at T 0 . Although D ∞ and the transition T 0 itself do depend on the specific choice of the correlator K it is interesting that the asymptotic state does not depend on the initial distance. For a better understanding of the physical origin of this transition we shall consider the case of equilibrium initial conditions in the next section.

C. Equilibrium initial condition
We have seen in the previous section that the asymptotic distance D ∞ is a non-trivial function of the temperature which is finite below T 0 and vanishes above T 0 . The relevance of the existence of the meta-stable states has been already pointed out in previous sections specially when the cross-correlator is maximal K = 1. The fact that the DS transition exists suggests that the nature of this phase transition is related to the corrugated properties of the free energy landscape. To check this result it is convenient to investigate the DS transition starting from an equilibrium configuration. In this case the system starts from a stationary state and remains there forever. At high temperatures this state is the paramagnetic so in this case the DS is a direct check of the ruggedness of the paramagnetic state. In fact, we will find that if we start from an initial equilibrium condition 26 , then the DS transition persists and actually coincides with the previous T 0 found for a random initial configuration. This reinforces the idea of D ∞ as a dynamical order parameter for the DS transition for a given choice of the correlator K. Again, for sake of simplicity, we restrict our analysis here to the case K = 1.
The analysis of the dynamical equations for an equilibrium initial condition follows the same steps as for the random case but now we must impose a Gibbs distribution for the configurations σ and τ at time 0. Nevertheless, there is a point which must now be considered. Let us take T > T s (i.e we will suppose equilibrium configurations in the paramagnetic phase). Suppose that p is odd and we take an equilibrium configuration at temperature T . To impose Q 0 = −1 or D 0 = 1 we must take σ i = −τ i for all i. Because the energy is an odd function of the spin variables we have E({σ}) = −E({τ }). If the equilibrium energy is not zero (this happens everywhere except at β = 0) we cannot put both configurations at equilibrium at the same temperature (because the temperature uniquely determines the value of the equilibrium energy). Then, if both initial conditions σ and τ are equilibrium initial configurations their overlap Q(0) must vanish. Actually, for T > T s two equilibrium configurations do have overlap zero with probability 1 and overlap different from zero with probability exponentially small with N . So if we take the thermodynamic limit before the infinite time limit it is clear that we must start with zero initial overlap. To be more precise, the probability that two equilibrium configurations {σ}, {τ } in the paramagnetic phase do have overlap q = 1 N N i=1 σ i τ i is given by, where f (q) is the free energy cost to find a correlation q between the configurations. Clearly, because q = 0 corresponds to the equilibrium value in the paramagnetic phase, f (q) has a minimum at q = 0 so we can write f (q) = constant + q 2 /(2βχ SG ), with χ SG = N < q 2 > is the spin-glass susceptibility. Above T s the χ SG is finite and the probability to have q = 0 is exponentially small with N . Now the cost in free energy f (q) has two parts, a cost in energy u(q) and a cost in entropy s(q) = β(u(q) − f (q)). The cost in energy vanishes at infinite temperature and the full cost of f (q) is due to the entropy. So only at infinite temperature (i.e random initial configurations, the case considered in preceding sections) we can impose an initial condition with initial non-zero overlap.
The equivalent of the dynamical equations eqs. (19,20,21) can be easily obtained for T > T s in the replica symmetric approximation 30 . The correlation function C(t, s) and the response function R(t, s) are time-translational invariant and satisfy the fluctuation-dissipation theorem T R(t) = − ∂C(t) ∂t Θ(t). The C(t) satisfies the equation, with C(0) = 1. The two-times overlap satisfies the equation, with the initial condition Q(t, 0) = Q(0, t) = 0. The diagonal part, Q d (t) = Q(t, t) is given by, with the initial condition Q d (0) = 0. Now we are in equilibrium so µ(t) = T + pβ 2 30 . We have looked for a timetranslational invariant solution for Q(t, s) (i.e a solution of the type Q(t, s) = Q d (s)Q(t − s) for t > s) but we have not found it (even for p = 2). Our numerical results suggest that such an asymptotic solution does not exist. Using as before a series expansion in the time-dependent variables and Pade approximants, we have estimated the asymptotic distance for equilibrium initial conditions, i.e., for Q(0) = 0. The results are displayed in figure 6 for different temperatures. The divergence of the relaxation time leads to T 0 (p = 3) = 1.01 ± 0.04 with γ = 1.4 ± 0.3 which indicates that T 0 coincides with the result obtained for random initial conditions. This supports the idea that the transition at T 0 is of a dynamical nature and unrelated to the existence of a thermodynamic phase transition. D. The non-symmetric α = 1 case As we saw in the previous section, one of the most interesting results concerning the DS transition is the fact that D(∞) is a non-trivial quantity which does not depend on the initial condition (and depends only on the cross correlation K). In equilibrium thermodynamics, this is one of the features of order parameters which separate different equilibrium phases. In the stationary state, when fluctuation-dissipation theorem is obeyed, the order parameter is a a quantity which characterizes the equilibrium state and (in the absence of ergodicity breaking) does not depend on the initial condition.
In order to present a convincing proof on this result we have investigated the general non-symmetric case α = 1. A difficulty inherent to the extrapolations made from figures 1 and 2 is the fact that, below T 0 , the convergence of the distance D(t) towards its asymptotic value D ∞ is very slow (a power law in time). Consequently, both numerically or using the Pade method it is very difficult to extrapolate to the asymptotic value. As the asymmetry of the interactions eq.(13) is turned on (i.e if α < 1) the relaxation of the system to the stationary state turns out to be faster. Actually for the asymmetric case (α = 0) or the antisymmetric case (α = − 1 p−1 , see later) the relaxation of the distance D(t) is nearly exponential. We have no reasons to suppose that the independence of the asymptotic value on the initial distance is α dependent. Our analysis for α < 1 suggests that the independence of D(t) on the initial distance D 0 holds for all generic values of α.
It is important to note that, for α < 1, there is no equilibrium stationary state and the fluctuation-dissipation theorem is not satisfied. Still we expect, for p > 2, the DS transition to survive for any quantity of finite asymmetry. The reason is that the DS transition (even for α = 1) is inherently a non-equilibrium transition so the effect of the asymmetry may not change the character of that transition.
The non-symmetric case for p = 2 has been already considered by Crisanti and Sompolinsky ten years ago 36 . By assuming that, in the stationary state, the correlation and response functions are time-translational invariant they succeeded in showing that the thermodynamic transition T s = 1 for α = 1 turned out to be unstable against the asymmetry for any value of α < 1. They also derived the result µ(∞) = √ 1 + T 2 for α = 0. In the next paragraphs we study the cases α = 0, − 1 p−1 starting from a random initial configuration in some detail. Unless stated we will consider the case K = 1.

The asymmetric case α = 0
The case α = 0 is quite interesting. The equation for the response function (20) simplifies considerably, which can be readily written, using R(t, s) = z(s) z(t) with z(t) = exp t 0 µ(t ′ )dt ′ . The equation for the correlation function becomes, Define the new function A(t, s) = z(t)C(t, s)z(s). In terms of this new function eq.(34) is, From eq. (22) it is easy to derive a equation for z(t), Equations eq.(35),(36) form a closed set of equations which can be solved with the initial conditions A(t, t) = z 2 (t), z(0) = 1. Once this set of equations is solved one can also find a solution for the overlap Q(t, s) in eq. (21). Again, we define B(t, s) = z(t)Q(t, s)z(s) which satisfies the equation, (37) and the equal-times overlap eq.(23) b(t) = B(t, t) satisfies the equation, (38) with b(0) = Q d (0). Note that this set of equations is quite involved for p > 2. Only for p = 2 they dramatically simplify (the case considered by Crisanti and Sompolinsky) and become linear. For general p the previous equations are non-linear. We have not succeeded in finding the asymptotic solution of these equations although we have guessed the results from the numerical results. We find that the DS transition is still present at finite temperature for p > 2. The analytical expression for T c is given by Note that T c coincides precisely with the lower bound previously derived in eq. (25). The asymptotic value of µ(t) is given by µ(∞) = √ 1 + T 2 and is p-independent. The asymptotic distance for p = 3, 4 is given by 37 , A full theoretical derivation of this results remains an interesting open problem. In what follows we compare our results obtained by numerical integrations with time step ∆t = 0.01 with these theoretical guesses. In figure 7 we show the overlap Q d (t) as a function of time for different temperatures below T c for p = 3. Note that the asymptotic value clearly does not depend on the initial condition. The horizontal dotted lines correspond to the asymptotic value eq.(40). This figure unambiguously demonstrates that the asymptotic distance does not depend on the value of the initial overlap Q d (0). In figure 8 we show µ(t) for p = 3, 4 compared with the theoretical prediction µ(∞) = √ 1 + T 2 .   This is an extremal case where J i2,i3,...,ip i1 J i1,..,i k−1 ,i k+1 ,..,ip i k is maximally negative. Physically this means that the force felt by an spin i due to a multiplet M of p − 1 spins is as much opposite as possible as the force which feels another spin contained in that multiplet due to the action of another multiplet M ′ of p − 1 spins constructed from the rest of the p − 2 spins in the previous multiplet M plus the spin i. In the particular case p = 2 this can be easily achieved making J j i = −J j i which corresponds, according to eq. (13), to α = −1. But in general, α can never be equal to −1 for p > 2. Take p = 3 and three couplings J jk It can be easily shown that the minimum value for α is given by α = − 1 p−1 . Interestingly, this case can be exactly solved for the correlation and response function. Although it turns out to be quite difficult to solve for the overlap function we will analyze here a general correlator K. For α = − 1 (p−1) eqs. (19)(20)(21)(22)(23) considerably simplify because µ(t) in eq.(22) is time independent. Because the initial configuration was taken random at time 0 this means that the stationary state follows completely random configurations.
To solve the equations define the following correlators, In this case the dynamical equations for C and R simplify. A particular solution for c(t, s) and r(t, s) can be found which simplifies dramatically the dynamical equations. This solution is given by c(t, s) = r(t, s) = f (t − s). This time-translational invariant solution is consistent with all dynamical equations for all times. The final closed equation for f (t) is given by, with the initial condition f (0) = 1. But the equation for the overlap Q(t, s) is more complicated and cannot be reduced to a time-translational invariant solution. Writing Q(t, s) = q(t, s) exp (−T |t − s|) (with q(t, t) = Q d (t) we obtain the following equations, A time-translational invariant solution for q(t, s) does not exist for all times (contrarily to what happens for C and R) because the first integral in the right hand side of eq. (44) does not vanish.
Previous equations are solvable for p = 2 and K = 1. In this case we get, and there is no DS transition for p = 2 as expected (i.e. T 0 = 0). For a generic cross-correlation K let us note that the equation (45) reduces to equation (24) so there will be different asymptotic values depending on the value of K ′ (Q * ) at the different set of fixed points Q * = K(Q * ). It is notorious that the case p = 2, α = −1 at finite temperature reduces to the infinite temperature case for any α. This is closely related to the fact that the stationary solution in this case coincides with the random initial configuration although this is not true anymore for p > 2. For p > 2 the lowest T 0 temperature for maximal cross-correlation K = 1 becomes finite and this is due to the rugged structure of the force landscape.
Equations for the damage for p > 3 are difficult to solve. We have not succeeded to obtain an analytical expression for the asymptotic values as well as for T c and D ∞ . Numerical integrations of the equations reveal that the transition persists at finite temperatures for p > 2. Figure 9 shows the overlap Q d (t) as a function of time for different temperatures for p = 3, K = 1. In this case the transition is located at T 0 ≃ 0.595 ± 0.005. Note that again relaxation to the stationary state is faster than in the case α = 1 and D ∞ is again independent of D(0). These results are quite appealing. Here we find a DS transition in the presence of a time-translational invariant solution for C and R, i.e when the system starts already in the stationary state. This is in agreement with the results of previous subsection III.B for α = 1 where a DS transition was found (at the same temperature than starting from random initial conditions) when the system already started in the stationary state. Figure 10 summarizes our results. We show the T − α phase diagram of the DS transition for p = 3. Let us remark our final conclusion for this section. A DS transition is present for all models with p > 2 above the TAP temperature where an exponentially large number of states appear. For a given choice of K such that Q = 1 is stable at infinite temperature, this transition has the following interesting properties: a) The asymptotic distance is independent of the initial distance; b) Is also independent of the type of initial configuration; c) Is stable in the presence of asymmetry. (but is unstable for p = 2 in agreement with results derived by Crisanti and Sompolinsky). This suggests that D ∞ has some of the crucial properties to be a good dynamical order parameter. The correlator Q(t, s) is not time-translational invariant in the time scale in which the correlation and the response are.
As we said previously we expect the properties of this transition to strongly depend on the type of dynamics through the cross-correlator K. In the next section, we will discuss discrete (Monte Carlo) dynamics, the case in which different algorithms correspond to different cross-correlators.

IV. DAMAGE SPREADING IN DISCRETE GLASSY MODELS
Up to now we have considered the DS problem in case of a dynamics continuous in time such as Langevin dynamics.
Here we want to investigate damage spreading and in particular the existence of T 0 for discrete dynamics such as Monte Carlo algorithms.
The analogue of the functions K for continuous dynamics are the different algorithms used in Monte Carlo dynamics for discrete dynamics. So different Monte Carlo algorithms determine different types of correlations between the noises. As for the cross-correlator K this implies that the algorithms determine the structure of the high-temperature fixed points. Their instabilities determine the subsequent low-temperature behavior. As we have already commented in the introduction we therefore do not expect the DS transition to be universal and the results of this section aim to be compared with the results already obtained for the Langevin dynamics in a continuous system.
One of the essential ingredients for the DS transition is the presence of two competing effects: a synchronizing force (the stochastic noise η in the Langevin dynamics) and a landscape-dependent force which pulls configurations apart one from each other into different directions. In the case of a discrete (Monte Carlo) dynamics the equivalent role of the stochastic noise is played by the set of random numbers generated during the Monte Carlo updates. Now, the random number in the Monte Carlo algorithm (uniformly chosen between 0 and 1) determines the probability of a move depending also on the configuration of the system. This last dependence corresponds to the role played by the cross-correlator K in the Langevin case where the two noises η, η ′ may be different depending on the value of K(C, C ′ ). Roughly speaking the Metropolis algorithm for Monte Carlo dynamics corresponds to the case K(Q) = Q for Langevin. It is easy to check that, at infinite temperature, the fixed-points in both dynamics are the same.
For continuous (Langevin) dynamics we had the freedom to choose the maximal cross-correlation K = 1. For discrete dynamics, however, this is not the case. There are several well known algorithms in the Monte Carlo approach according to which updating rule they use. For instance, Metropolis, Glauber, or heat-bath. Among these, the last one is the only one which has a unique fixed point Q * = 1 at infinite temperature. So, heat-bath dynamics is the closest case (but different) to the K = 1 of Langevin dynamics. Here, our numerical investigation will focus on this type of discrete dynamics. Let us note that the other algorithms may show different behavior (due to the presence of other infinite-temperature fixed points) and consequently also different DS transitions. This non-universality of the DS transition (as in our previous analysis of the Langevin case) has received some attention in the literature 22 .

A. Damage spreading in the Bernasconi model
Here we will analyze the Monte Carlo dynamics with the heat-bath algorithm for the Bernasconi model 42 . This is a long-range interaction model without disorder which is known to have a glassy behavior being in the universality class of spin-glass models with one step of replica symmetry breaking 43 . Consequently, its dynamical behavior is the same as predicted by the mode-coupling theory.
The Bernasconi model (for simplicity we will consider the closed model, see 43 for more details) consists of N Ising spins σ i = ±1 in a one dimensional chain interacting through a long-range four-spin interaction. It is defined by the following Hamiltonian, with C k = N j=1 σ j σ j+k and where we take periodic boundary conditions σ i = σ i+N . In this model there are particular values of N for which the ground state is exactly known 43 . The interest of this model is that it behaves like a disordered spin glass in the absence of explicit quenched disorder in the Hamiltonian. Apparently, disorder is self-induced by the dynamics 41,[43][44][45] . This means that dynamics itself generates slow evolving variables which effectively act as quenched disordered fields. This model is characterized by three temperatures, a melting crystal-liquid first-order transformation temperature T M , a dynamical transition temperature T d ≃ 0.5 46 where the relaxation time diverges and ergodicity breaks and, finally, a static (or glass) transition temperature T s ≃ 0.25 where replica symmetry breaks and the configurational entropy collapses (this is the ideal glass transition predicted in the AGM theory).
In the heat-bath algorithm, to go from a configuration {σ i (t)}, at a given time t, to a configuration {σ i (t + ∆t)} at the next time-step t + ∆t with ∆t = 1/N , a spin σ k is chosen at random among the N spins to be updated. The probability to put the spin up or down is decided according to the intensity of the local field acting on that spin. More precisely, if we write the Hamiltonian (48) in terms of the local field H = − k h k σ k then the probability of putting the spin σ k up (σ k = 1) or down (σ k = −1) at time t + ∆t is given by, with σ = ±1 and h k (t) is the local field acting on the spin k at time t. Note that the probability (49) does only depend on the local field acting on the spin k and not on the actual value of that spin at time t. Then, a random number z(t) with a uniform distribution between 0 and 1 can be introduced and spins are sequentially updated according to the dynamical rule σ k (t + ∆t) = sign 1 2 With this rule (50) we have studied numerically the damage spreading of three different initial conditions, as in the p-spin model: a) anti-correlated random initial conditions with D(0) = 1, b) uncorrelated random initial conditions with D(0) = 0.5 and c) partially correlated random initial conditions with D(0) = 0.1. For each of these cases, the distance D(t) is computed up to 100000 and 10000 Monte-Carlo time-steps for N = 1000 and N = 5000 respectively. To analyze the data, the logarithmic time with base a = 1.1 is considered. Moreover, the data is averaged in intervals of the form (a k , a k+1 − 1) with k a positive integer. For T = 0.3 the evolution of D(t) is plotted in figure 11. One observes that the distance does not vanish for any finite temperature. Moreover, figure 12 indicates the existence of a temperature T 1 above which the asymptotic distance is independent of the initial distance. Below this temperature T 1 , D ∞ does seem to be dependent of the initial distance. This dependence is supported by a numerical extrapolation which could well fail when going to enormously large time scales. Still, what we certainly find is the appearance of a dynamical transition temperature T 1 ≃ 0.5 in very good agreement with the transition T d where ergodicity breaks 43,46 .
The behavior we find here, when compared to the previous Langevin analysis for the p-spin model, may appear quite different. But a careful analysis reveals that this is not the case. If we consider that K, for the Langevin case, is a generic function which may depend on the overlap as well as on the temperature T = 1/β, we may then imagine a situation such as depicted in figure (13) where the infinite-temperature fixed-point Q = 1 becomes unstable as soon as β is finite. In this case the asymptotic distance would be a non trivial function of β and the damage-spreading transition could well happen at the usual dynamical transition T d where ergodicity is broken. The dependence of the asymptotic overlap on the initial value could be a consequence of the presence of different fixed points at low temperatures.
In the most general case one could imagine a scenario with three possible different regimes. A high-temperature regime T > T 0 where D ∞ = 0 independently of the initial distance D 0 . An intermediate regime is not zero but independent of the initial distance (this regime would correspond to the appearance of temperature dependent fixed point for β finite as depicted in figure (13)). And finally a low-temperature regime T < T 1 where D ∞ = D ∞ (T, D 0 ) depends on both temperature and initial distance. The results we find for the Bernasconi model are the same as found by Derrida and Weisbuch 39 for the Sherrington-Kirkpatrick model. The fact that in this model T 0 → ∞ is related to the infinite-range character of the interaction. Actually, Derrida 15 has found numerical evidence that for finite-dimensional spin glasses there exists a range of temperatures where the asymptotic distance vanishes and T 0 is finite. The dependence of the asymptotic overlap on the initial condition found here and in 15 below T 1 could well be an artifact of the large-time extrapolation where the time window simulated and the size as well are not sufficiently large. Unfortunately, it is not easy to simulate very large times and sizes in infinite-ranged models like the present one. The results of this section show that, the DS transition is very close (and probably coincides) with the dynamical transition temperature T d (below which the system never attains equilibrium and ergodicity is broken). Nevertheless, in this case the asymptotic damage below T 1 apparently depends on the initial condition (and probably on the type of initial condition as well) although such a firm conclusion needs more understanding of damage-spreading transitions for generic updating rules.

V. CONCLUSIONS
In this paper we have studied the problem of damage spreading in the mode-coupling theory of glasses. Modecoupling theory is well known to describe relaxation processes in glasses in the mean-field approximation. A simple way to obtain the mode-coupling equations is analytically solving the dynamics of multipsin interaction spherical spinglass models. These models are characterized by the presence of a huge number of meta-stable states (exponential with the system size) which appear at a temperature T T AP higher than T d (where ergodicity breaks) and T s (where replica symmetry breaks). Because the phase space in this class of models is characterized by a extremely rugged and complex free energy landscape they are good candidates to study the landscape properties using techniques taken from dynamical systems.
A very interesting technique which is able to probe the topological features of the phase space landscape is damage spreading. This consists in the study of the distance between the configurations of two stochastic systems submitted to the same realization of the stochastic noise but differing in the initial conditions. By the same realization of the stochastic noise we mean noises that are statistically identical although generally correlated through a function K(Q) which satisfies the condition K(Q) ≤ K(Q = 1) = 1. In general, any choice for the correlator K alters the results. For Langevin dynamics, we have shown that interesting results appear for the case K = 1. In that case, both noises are identical for the two copies independently of their configurations. This yields the lowest damage spreading transition for which the Q * = 1 fixed-point is stable at high-temperatures. Depending on the value of K ′ (Q = 1) one finds a different damage spreading transition temperature up to the limiting case K ′ (Q = 1) = 1 (see eq.(25)) where T 0 = ∞ and the fixed-point Q * = 1 becomes unstable. Whether this holds for other type of dynamics is not studied and remains unclear.
An exhaustive study has been done for the case K = 1 although similar results are obtained for any K for which the solution Q = 1 is asymptotically stable. In this case, through functional methods and using the Pade series expansion method (to make safe extrapolations in the asymptotic long-time limit), we have shown the existence of a damage spreading transition T 0 in general mode-coupling equations with any degree of asymmetry in the interactions. In particular, in case of symmetric interactions (where detailed balance holds) we have found evidence for a damage spreading transition at a finite temperature. This transition occurs at temperatures T 0 higher than T T AP , this last one being the temperature where an exponentially large number of meta-stable states start to appear. The transition is characterized by a dynamical order parameter D ∞ which is the asymptotic distance between the two evolving replicas. Interestingly, D ∞ has the good properties of order parameters being able to distinguish different dynamical phases (in our case, there are two possible phases depending whether D ∞ vanishes or not). These properties are: a) D ∞ is independent on the initial distance D 0 for a given class of initial conditions, b) D ∞ does not depend also on the class of initial condition (whether they are random or thermalized) c) The DS transition is stable against the inclusion of asymmetry in the interactions (i.e. against the violation of detailed balance in the dynamics). Furthermore, regarding the mode-coupling equations with asymmetry we have obtained some exact results for the asymmetric case α = 0 and exactly solved the correlation and the response function for the antisymmetric case α = − 1 p−1 which interestingly turns out to be time-translational invariant. The existence of DS transition in this case reveals that this transition already appears when the system is time-translational invariant.
We insist on the fact that the precise value of the damage spreading temperature T 0 as well as the asymptotic distance D ∞ both depend on the correlator K considered. This fact expresses the non-universal character of this transition where the cross-correlator K plays the equivalent role of a stochastic noise for the dynamical order parameter Q d (t). Different functions K imply different dynamical phase transitions so their physical significance must be appropriately interpreted. In this direction we have tried to interpret our results in terms of saddle points in phase space for those cases where Q * = 1 is an asymptotically stable solution at high temperatures. The essentials of the argument where given in section II where it was shown that for close enough initial conditions an instability in the Q = 1 hightemperature solution is due to the presence of vanishing modes in the Hessian of the potential function. This signals the presence of saddle points in the free energy landscape, of which an infinite amount yields a DS transition. This certainly happens at T T AP but we have found that the transition occurs already at a temperature T 0 much above that temperature. How can we reconcile our results with that? At T T AP the number of meta-stable states is exponentially large with N . Probably for the DS transition to occur it is only necessary that this number be a big (for instance , a power of N ) but not exponentially large. Consequently, a divergent number of meta-stable states and saddle points (but not exponentially large with N ) should be enough to make the DS transition appear. Unfortunately, the analytical computation of the temperature T 0 by counting the number of meta-stable states is not so direct because such a calculation involves the estimate of finite-size corrections to the dominant saddle point calculation 47 . Actually, the evaluation of finite-size corrections in spin glasses (even in mean-field) is known to be quite difficult. Concerning the role of the cross-correlations between the noises we remind that the appearance of vanishing modes tends to be suppressed by cross-correlations K < 1 which play the role (according to eq.(10)) of a thermal noise on the system described by the effective distance variables z i = x i − y i . This explains why K = 1 yields the lowest damage spreading temperature. We must stress that this transition should persist beyond mean-field. Actually the presence of a huge number of saddle points in phase space is not an exclusive mean-field feature but should persist in the presence of short-range corrections. Short-range corrections dramatically modify the height of the barriers but not their number which still could remain as large as in the mean-field approximation.
To remark the advantage of considering Langevin dynamics for DS we have also investigated the DS in discrete glassy models without disorder such as the Bernasconi model with heat-bath dynamics. The difference between any discrete dynamics and Langevin dynamics is that the type of correlator K in the last case may be chosen at will. So the case K = 1 which has been studied for Langevin dynamics cannot be implemented for Monte Carlo dynamics. Each algorithm for Monte Carlo dynamics defines a given K so the study of damage spreading in those cases remains more speculative because K is essentially unknown. For the particular case of the heat-bath in the Bernasconi model we find a DS transition which separates two regions, a region T > T 1 where D ∞ is independent of the initial distance but finite and a region T < T 1 where D ∞ depends on the initial distance. T 1 coincides (within numerical precision) with the mode-coupling transition T d . Although for heat-bath dynamics D ∞ depends on the initial distance we may not exclude that convergence to the asymptotic limit is extremely slow and convergence poor.
We end up our discussion describing some opening problems. Much work on the DS has been devoted to the study of discrete dynamics in discrete systems the situation being the opposite for continuous-time dynamics. In this work we have shown that for Langevin dynamics we may choose specific noise correlations so it is easier to interpret what physical properties of the system we are looking at. It is difficult to adscribe any physical significance to properties measured for arbitrary discrete algorithms and this has been one of the major problems to interpret the large amount of numerical results obtained in the study of damage spreading in Monte Carlo simulations. Here we have seen that for Langevin dynamics such a task turns out to be easier. Still one of the major tasks which remain open is to understand better under which conditions there is a unique absorbing state for the damage. In another direction, one would like to get analytical proofs and equilibrium based analytical methods to investigate the connection between the DS transition and the topological properties of the phase space (such as the presence of saddle points in phase space). Finally, we would like to extend our research to real structural glasses where it could be interesting to study DS and investigate in which conditions such a transition is a precursor of the glass transition.

APPENDIX A:
In this appendix we describe some of the main steps necessary to derive equations (19)- (23). Consider two replicas {σ}, {τ } submitted to the dynamics (11), where F i is the force acting on the spin σ i due to the interaction with the rest of the spins eq. (12). The noises η,η ′ are generically correlated < η i (t)η ′ j (s) >= 2T K({σ}, {τ })δ(t − s)δ ij . Following Domany and Hinrichsen 22 we will suppose that the correlator K is a generic function of the overlap due to the infinite-ranged nature of the model. So we will take K({σ}, {τ }) a generic function of the equal-times overlap K(Q d ) where Q d = 1 N N i=1 σ i τ i . Following the same steps as is usually done in the study of the dynamics of a single replica (see 30 for details) we may write the generating functional for the dynamics in the Ito prescription, Introducing a new set of fieldsσ i ,τ i and averaging over the noise we get, Because Z dyn = 1 we may average the dynamical partition function over the disorder. We use the cumulant expansion and retain only the first and second order terms exp(V ) ≃ exp V + 1 2 (V 2 − (V ) 2 ) where (.) stands for disorder average. Using eq. (14) we obtain the final result, so, in the long-time limit t → ∞ where the C(t) vanishes we get,