{'fraction_non_alphanumeric': 0.043870785724361684, 'fraction_numerical': 0.012766257582169559, 'mean_word_length': 4.076987347815708, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 26, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spin Glasses on the Hypercube

We present a mean field model for spin glasses with a natural notion of distance built in, namely, the Edwards-Anderson model on the diluted D-dimensional unit hypercube in the limit of large D. We show that finite D effects are strongly dependent on the connectivity, being much smaller for a fixed coordination number. We solve the non trivial problem of generating these lattices. Afterwards, we numerically study the nonequilibrium dynamics of the mean field spin glass. Our three main findings are: (i) the dynamics is ruled by an infinite number of time-sectors, (ii) the aging dynamics consists on the growth of coherent domains with a non vanishing surface-volume ratio, and (iii) the propagator in Fourier space follows the p^4 law. We study as well finite D effects in the nonequilibrium dynamics, finding that a naive finite size scaling ansatz works surprisingly well.


I. INTRODUCTION
Spin Glasses (SG) are highly disordered magnetic systems [1]. Rather than by their practical usefulness, SG are often studied as a paradigmatic example of a complex system. Indeed, they display an extremely slow dynamics on a complex free-energy landscape, with many degenerate states. In addition, SG are a convenient experimental model for glassy behavior, due to the comparatively fast microscopic spin dynamics, as compared, for instance, with supercooled liquids. In fact, nowadays, SG are considered as a playground to learn about general glassy behavior, minimization problems in Computing Science, biology, financial markets, etc.
Maybe the most conspicuous feature of SG is Aging: they never reach thermal equilibrium in experimental times. Here we will only consider the simplest possible experimental protocol, the temperature quench (see [2] for very interesting, more sophisticated experimental procedures): the sample is cooled below the critical temperature, T c , and it is let to relax for a time t w at the working temperature T . Its properties are studied at a later time t + t w . It turns out that, if a magnetic field was applied from the temperature quenched until t w , when it is switched off, the thermoremanent magnetization M (t, t w ) decays with t/t w (at least for 10 −3 < t/t w < 10 and 50 s < t w < 10 4 s [3]). This lack of a characteristic time beyond t w , the glassy system age, is known as Full Aging. We now know that Full Aging is an effective description of the dynamics, no longer valid for t/t w ∼ 10 4 and t w ∼ 10 s [4].
Nowadays, we know that the slow dynamics in SG is originated by a thermodynamic phase transition at T c [5]. Below T c , the spins associate in coherent domains, whose size, ξ(t w ), grows with time. The lower the temperature, the slower the growth of ξ(t w ) is [experimentally, ξ(t w = 100 s; T = 0.9T c ) ∼ 100 atomic spacings].
There is a lively theoretical debate on the properties of the low temperature phase. Surprisingly enough, this controversy on the equilibrium properties of a non-accessible (in human time scales) spin glass phase is relevant to nonequilibrium experiments [6].
Mainly, there are two competing theories, the droplets [7], and the replica symmetry breaking (RSB) one [8]. According to the droplets picture, the SG phase would be ferromagnetic-like, with a complicated spin texture, but essentially with only two equilibrium states. On the other hand, the RSB theory predicts an infinite number of degenerated states with an ultrametric organization. For both theories, Aging would be a coarsening process, in the sense that coherent domains of low temperature phase grow with time. The two theories disagree in their predictions for these domains properties. According to droplets, the domains would be compact objects, with a surface-volume ratio that vanishes in the high ξ(t w ) limit [9]. The SG order parameter is non zero inside of each domain. On the contrary, the RSB theory expects noncompact domains with a surface-volume ratio constant for large ξ(t w ). Furthermore, in a RSB system, the SG order parameter vanishes inside those domains. In recent times, a somehow intermediate theory (TNT), has been proposed [10], but, to our knowledge, no detailed dynamic predictions have been provided.
The RSB theory is based on the Mean Field approximation (MF) which, unlike in the ferromagnetic case, is highly non trivial. Indeed, after 30 years of study, Aging is not yet fully quantitatively understood, not even within MF approximations. Therefore, non perturbative tools, such as Monte Carlo (MC) calculations, appear as an appealing alternative. Furthermore, MC calculations are called for even at the MF level. Hence, one is interested in mean field models, that is to say models where the MF approximation becomes exact in the thermodynamic limit (TL). The standard MF model, (the Sherrington-Kirkpatrick model, see [11] and Sect. II), has a number of disadvantages. It lacks a natural notion of distance (hence one cannot discuss a coherence length ξ(t w )), or coordination number. Furthermore, its numerical simulation is computationally heavier. In fact, recent advances on the analytical study of spin glasses on Bethe lattices [12] has shifted the attention to these far more numerically tractable models which share with experimental systems the notion of a coordination number.
Here we wish to present a new MF model for spinglasses: the spin-glass on a D-dimensional hypercube with fixed connectivity [13]. In the thermodynamic limit (which coincides with the large D limit for this model), Bethe approximation becomes exact. As a consequence, the statics is of Bethe-lattice type and can be computed. A nice feature of this new model, is that it has a natural definition of distance, which allow us to study spatial correlations within MF approximation. In other words: this MF model is more similar to a real D = 3 system than any other studied before or, at least, than those considered previously, since the space-time correlation functions can be studied.
The structure of this paper will be the following. In Sect. II we will describe the model and compare it with other MF models. In particular, in Sect. II A we will address the problem of fixing the connectivity in a diluted hypercube. In Sect. III we will briefly explain the numerical methods we have used and in Sect. IV we will introduce the observables measured during the simulations. In Sect. V and VI, the numerical results will be presented, both in equilibrium (as a test of the model) and nonequilibrium, respectively. In Sect. VII we will discuss finite size effects. The analysis will reveal the p 4 propagator [14], for the first time in a numerical investigation. Our conclusions will be presented in Sect. VIII. Finally, we include extended discussions in two appendices.

II. MODELS
The standard model in SG is the Edwards-Anderson (EA) model. We will work with two kinds of degrees of freedom: dynamical and quenched. The dynamical ones correspond to the spins, σ i , with i = 1, 2, . . . N . We will consider them as Ising variables ±1. The non dynamical (or quenched) represent the material impurities. We will consider here two types of them: the connectivity matrix, n ik = n ki = 1, 0 (n ik = 1 as long as spins i and k interact), and the coupling constants, J ik = J ki , which shall take only two opposite values (in general with certain exceptions that will be discussed in the next paragraph, we will consider J ik = ±1, which defines our energy scale). The interaction energy is Since the impurity diffusion time is huge compared to the spin-flip (picosecond), we will always work within the so-called quenched approximation: spins cannot have any kind of influence over the material impurities. Then, both the set of coupling constants in the Hamiltonian and its associated Gibbs free energy, will be considered random variables. Therefore, in order to rationalize the experiments, the useful free energy will be an average over the disorder. We will refer to each assignment of {n ik , J ik } as a sample. Its probability distribution defines the actual EA model. The average over samples will be represented as (. . .).
Only a few exact results are known for the Hamiltonian (1), and all them were obtained within the MF approximation [8]. This approximation becomes exact in a weak infinite-range interaction model (in a ferromagnet n ik = 1 and J ik = 1/N for each couple i, k). In SG it is usually represented as n ik = 1 for every couple i, k and, because of the random ferromagnetic and antiferromagnetic interaction character, {J ik } are independent gaussian random variables, with zero mean and variance 1/N [11].
Computer simulations of long-range models are extremely hard, because the energy evaluation for a system of N spins requires N 2 operations. The situation has improved since the discovery that EA models on Bethe lattices (not to be confused with Bethe trees) undergo Replica Symmetry Breaking at T c [12].
A very popular realization of a Bethe-lattice spin glass is the EA model on a Poisson graph. A simple way of drawing one of these graphs consists on connecting (n ik = 1) each possible pair of spins, i, k, (there are N (N −1)/2 possible couples), with probability z/(N −1). Thus, the number of neighbors of spin i, its coordination number n i follows in the large-N limit a Poisson distribution function with average z (the connectivity). We will consider z = 6 to mimic a three dimensional system. This kind of graphs are locally cycle-less: the mean shortest length among all the closed loops that passes through a given point is O(log N ), i.e. the system is still locally tree-like. Computationally convenient as they are, Poisson graphs still lack a natural notion of distance.
A simple alternative consists on formulating the model on a D-dimensional unit hypercube. Thus, the spins are located in each of the hypercube vertices (then, N = 2 D ) and the bonds lie on the edges. Therefore, each spin can be connected with, at most, D = log 2 N spins. By analogy with the Poissonian graph, we consider that a link is active (i.e. n ik = 1) over each edge with probability z/D. We call this model random connectivity hypercube. It is easy to prove that it is locally tree-like as well: the density of closed loops of length l decays, at least, with D −2 (i.e. with the squared logarithm of N , as it also happens in the Poisson graph). Incidentally, one could consider as well a non diluted hypercube, but this would have two shortcomings: the connection with three dimensional systems would get lost, and the numerical simulation would become computationally heavy for large D.
Note that, at variance with other infinite-dimensional graphs, the hypercube has at least two natural notions of distance: Euclidean metrics and the postman metrics [31]. The two distances are essentially equivalent, since the Euclidean distance between two sites in the hypercube is merely the square root of the postman distance. In the following we shall use the postman metrics, which has some amusing consequences. For instance, our correlation-length will be the square of the Euclidean one, thus yielding a critical exponent ν = 1, doubling the expected ν MF = 1/2. Of course, if we use the Euclidean metric we recover the usual exponent ν = 1/2.
However, it turns out that the random connectivity hypercube suffers a major disadvantage. The inverse of the critical temperature in a ferromagnet [15] or in a SG [16] can be computed within the Bethe approximation: In this expression n 1 is a conditional expectation value for n, the coordination number of a given site in the graph. This conditional expectation value is computed knowing for sure that our site is connected to another specific site (this is different from the average number of neighbors of a site that has at least one neighbor!). A simple calculation shows that n 1 = 1 + z − z D in the random connectivity model. Since D = log 2 N , we must expect huge finite size corrections (O(1/ log N )) at the critical point. Note that this problem is far less dramatic for a Poisson graph where n 1 = 1 + z − z N −1 . The cure seems rather obvious: place the occupied links in the hypercube in such a way that n = z (here, z = 6). Unfortunately, drawing these graphs poses a non trivial problem in Computer Science [17]. Our solution to this problem is discussed in the next paragraph.

A. The fixed connectivity hypercube
We have not found any systematic way of activating links in the hypercube that respects the fixed connectivity condition. Thus, we have adopted an operational approach: the distribution of bonds is obtained by means of a dynamic MC. We must define a MC procedure that generates a set of graphs that remains invariant under all symmetry transformations of the hypercube group.
Specifically, we start with an initial condition in which all bonds along the directions 1 to 6 are activated (of course, this procedure makes sense only for D ≥ 6). Clearly enough, the initial condition verifies the constraint n = 6. We shall modify the bond distribution by means of movements that do not change n. We perform what we called a "plaquette" transformation (a plaquette is the shortest possible loop in the hypercube, of length 4). We randomly pick, with uniform probability, one hypercube plaquette. In case this plaquette contains only two parallel active links (n ik = 1), these two links are deactivated at the same time that the other two are activated. On the opposite case, nothing is done [32]. This guarantees that the set of generated graphs is isotropic.
In order to this procedure to be useful, the dynamic MC correlations times must be short. In Fig. 1, we show the MC evolution of the system isotropy. We make kN plaquette transformations, and we control the density of occupied bonds in two directions: the first direction (initially occupied in every vertex) and the seventh direction (initially unoccupied). As we see, for two different system sizes, we get short isotropization exponential times (for D = 22 we get τ exp ≈ 4.7N ). For this reason, we assume that taking k = 100 is long enough to ensure that the configurations obtained are completely independent from the initial condition. At this point, the question arises of the completeness of the generated set of graphs. We first note that our set contains proper subsets that are also isotropic. However, finding them would require more involved algorithms which will not pay in a reduction of statistical errors (as we will show below, most of the sample dispersion is induced by the coupling matrix {J ik }). On the other hand, one could think that there are lacking graphs in our algorithm for a simple reason. The plaquette transformation cannot break loops: when we interchange neighboring links we can only either join two different loops or split up a loop into two loops. Due to the hypercube boundary conditions, in the initial configuration all sites belonged to closed loops. This situation cannot be changed by plaquette transformations. However, this objection does not resist a close inspection. In fact, a non-closed lattice path formed by occupied links should have an ending point with an odd coordination number, which violates the constraint n = z for any even z. Thus, all lattice paths compatible with our fixed connectivity constraint, do form closed loops. This argument, as well as the numerical checks reported below, make us confident that the set of generated graphs is general enough for our purposes. Actually, we conjecture that our algorithm generates all possible fixed connectivity graphs with z even.
One may worry as well about the applicability of the Bethe approximation to the fixed connectivity model, since all loops are closed. Actually, the crucial point to apply the Bethe approximation is that the probability of having a closed path of any fixed length should vanish in the large D limit. This is easy to prove for the random connectivity model. In the fixed connectivity case, one may argue as follows. Let us imagine a walk over the closed path. On the very first step, the probability that the chosen link is present is z/D, whereas in the following step the probability of finding the link is (z − 1)/(D − 1) in the limit of large D (since one of the z links available at the present site was already used to get there). This estimate implicitly assumes that the occupancy of different links is statistically independent. The independency approximately holds for large D and becomes exact in the D → ∞ limit, where occupied links form a diluted set. At this point, the estimate of the number of paths of any given fixed length in the large D limit can be performed as in the random-connectivity case. One finds as well, in the fixed connectivity case, that the number of closed loops per site of a given length decays at least as In addition to the above considerations, we may numerically compute in our graphs the length of the second shortest path that joins two connected nearest neighbors in the hypercube. In Fig. 2, we compare the probabilities for the length of such paths in the random (top) and fixed (bottom) connectivity models, for different system sizes [33]. In both cases, we note that the maximum of the probability shifts to larger length as D grows. We note as well that, for fixed connectivity, no tree-like graph arises [34].
A summary of our efforts is shown in Fig. 3, where we plot the evolution of the critical point with D for the ferromagnetic Ising model, defined on hypercubes with both random and fixed connectivity. As expected, the random connectivity model suffers very important finite volume corrections which make it essentially useless for numerical studies. The problem is solved using fixed con- nectivity hypercubes instead, where the finite volume effects are only caused by the residual presence of short closed loops.

III. NUMERICAL METHODS
We have simulated the Hamiltonian (1) using a Metropolis algorithm [18]. In addition, we use Multispin Coding: since spins are binary variables, we can simultaneously codify 64 systems in one single 64 bits word (all of them share the same connectivity matrix, n ik ). With this common matrix, we find errors which are ∼7 times smaller than those obtained with one single sample per matrix. This should be compared with the factor 8 that we would obtain in case no correlation were induced. Our program needs 0.29 ns/spin-flip in an Intel i7 at 2.93GHz (in Ref. [19] they report ∼ 1.2 ns/spin-flip on an Opteron at 2.0 GHz, for the simulation of the D = 3 EA model in the cubic lattice).
In a nonequilibrium dynamical study such as ours, one computes both one-time and two-times quantities, see Sect. IV. The calculation of two-times quantities implies the storage on disk of intermediate configurations. Disk capacity turned out to be the main limiting factor for the simulation. For this reason, we have worked in parallel with two program versions: one valid for measuring quantities at one and two times and another restricted to the computation of one-time quantities.
We have computed two-time quantities at temperature T = 0.7T c , on systems with D = 16, 18, 20 and 22. The number of simulated samples were 8×64 samples for each system size (hence, for self-averaging quantities the statistical quality of our data grow with D). Besides

IV. OBSERVABLES
The Hamiltonian (1) has a global symmetry Z 2 (σ i → −σ i for all i). Not as obvious is the gauge symmetry induced by the average over couplings. In fact, choosing randomly a sign for each position, ε i = ±1, the energy (1) is invariant under the transformation Now, since the transformed couplings ε i ε k J ik are just as probable as the original ones, we need to define observables that are invariant under the gauge transformation (3). With this aim we form gauge invariant fields from two systems at equal time, that evolve independently with the same couplings, {σ i } (real replicas) or, alternatively, from a single system considered at two different times: We can define three kinds of quantities with both fields. 1. One-time-quantities. The order parameter vanishes in the nonequilibrium regime (the system is much bigger than the coherence length, ξ(t w )). The non linear susceptibility is proportional to the SG susceptibility: that grows with a power of ξ(t w ). The Binder parameter provide us with information about the fluctuations In the Gaussian regime B = 3. In a ferromagnetic phase, B = 1. In the SG phase, in equilibrium (that for finite volume corresponds to t w → ∞), B grows with the temperature from B = 1 at T = 0. The equilibrium paramagnetic phase is in Gaussian regime.
2. Two-time-quantities. The correlation spin function tells us about the memory kept by the system, at time t + t w , about the configuration at t w [35]: The susceptibility is On the other hand, when t w is fixed, C(t, t w ) is just the thermoremanent magnetization [36]. The link correlation function carries the information of the density of the interfaces between coherent domains at t w , that at t + t w , have flipped. In case the surface-volume ratio decayed with a negative power of ξ(t w ) (droplets), C link would become t-independent [20]. On the contrary, in a RSB system, C link = a + b C 2 . Note that, for the Sherrington-Kirkpatrick model, one trivially finds C link = C 2 , but the linear relation is not straight-forward in fixed connectivity mean field models.
3. Spatial correlation functions. In the unit hypercube, the binary decomposition of the spin index i = 1, 2, . . . , 2 D can be equal to its Euclidean coordinates, x. The spatial correlation function is We consider r = |r| as the distance in the postman metrics. It would look rather natural to average all the c 4 (r, t w ) over the N r = D r displacements of length r = |r|: However, see Fig. 4, C 4 (r, t w ) does not present a limiting behavior with D for a given t w .
We can get a clue by looking at χ SG (t w ), Fig. 5, which does reach a thermodynamic limit. Since χ SG (t w ) is nothing but the integral of C 4 (r, t w ) with the Jacobian D r , it seems reasonable to define the spatial correlation function instead: We can see thatĈ 4 (r, t w ) does reach the high-D limit, Figure 6, at least for short t w . Besides, in the paramagnetic phase, it is possible to compute analytically   C 4 (r, t w ), see Appendix A, taking first the limit t w → ∞ and making afterwards D → ∞. The resulting correlation function, which is only valid in the paramagnetic phase, is a simple exponential. Hence, both the equilibrium and the nonequilibrium computations, suggest that one should focus onĈ 4 rather than on C 4 . We note in Fig. 6, that in the SG phase,Ĉ 4 is non monotonically decreasing with r, but rather presents a maximum. This maximum moves to bigger r with t w , then, the system has a characteristic length that increases with time. Thus, we can estimate the coherence length, by means of the integral estimator ξ 0,1 (t w ): A major advantage of ξ 0,1 over more heuristic definitions of the coherence length, is that it is computed from selfaveraging quantities (see details in [20,21], we note that, in this work, we have not tried to estimate the contribution to the integrals by the noise-induced long distance cutoff).
The existence of such a characteristic length is the main advantage of the hypercube model over other MF models.

V. EQUILIBRIUM RESULTS
Since the present work is the first study ever made of a EA model on a fixed connectivity hypercube it is necessary to make a few consistency checks. Equilibrium results are most convenient in this respect, since we have analytical computations (valid only for the large D limit) to compare with.
We will briefly study the spatial correlations in the paramagnetic phase. In addition, we will check, by approaching to T c from the SG phase, that the SG transition does lie on the predicted T c , Eq. (2).

A. Paramagnetic Phase
Our very first check will be the comparison between the Monte Carlo estimate of the SG susceptibility (at finite D) with the analytical computation for infinite D: see Appendix A. Our results are presented in Table I. We see that finite size effects increase while approaching T c . For our larger system, D = 22, the susceptibility significantly deviates from the asymptotic result only in the range T c < T < 1.1T c . After the fast convergence to the large D limit observed in the SG susceptibility, the results forĈ 4 are a little bit disappointing. In Fig. 7 we displayĈ 4 (r, D)−Ĉ 4 (r, ∞) as a function of r. We can see that finite size effects become more important once one approaches T c . Besides, finite D corrections as a function of r oscillate between positive and negative values. This is not surprising: the finite D corrections to the susceptibility, which are very small, are just the integral under these curves. More quantitatively, we see in Table II that

B. SG phase
In the SG phase, our test has been restricted to a check of Eq. (2), that predicts a SG phase transition for the

VI. NONEQUILIBRIUM RESULTS
In this section we will address the main features of the nonequilibrium dynamics obtained in our largest system, D = 22. The issue of finite D corrections will be postponed to Sect. VII.

A. The structure of isothermal aging
The picture of isothermal aging dynamics in MF models of SG behavior was largely drawn in [22] (see also [23]). The dynamics is ruled by an infinite number of time-sectors: The scaling functions f i are positive, monotonically decreasing and normalized, i.e. 1 = i f i (1). The unspecified functions h i are such that, in the large t w limit, In other words, the decay of C between values C i and C i+1 is ruled by the scaling function f i and takes place in the time-sector t ∼ t µi w . This picture is radically different to the Full Aging often found both in experiments and in 3D simulations. A full aging dynamics is ruled only by two sectors of time, µ i = 0, 1. Nevertheless, recent experimental studies [4] show that full aging is no longer fulfilled for t ≫ t w . Probably more time-sectors must be considered to rationalize these experiments.
However, Eq. (15) is probably an oversimplification, since the spectrum of exponents µ i might be continuous. An explicit realization of this idea was found in the critical dynamics of the trap model [24], where the correlation function behaves for large t w as C(t, t w ) = f (α(t, t w )) , α(t, t w ) = log t/ log t w . (16) Again, the scaling function f is positive and monotonically decreasing. Clearly enough, in the limit of large t w and for any positive exponent µ, if t = At µ w , the correlation function takes a value that depends only on µ, no matter the value of the amplitude A.
As expected, C(t, t w ) is clearly not a function of t/t w in our model, see Fig. 9. On the contrary, data seem to tend to a constant value when t w → ∞ in any finite range of the variable t/t w . This is precisely what one would expect in a time-sectors scheme. On the other hand, if we try (without any supporting argument) the Bertin-Bouchaud scaling, Eq. (16), see Fig. 10, the data collapse is surprisingly good. Therefore, the nonequilibrium dynamics in the SG phase seems ruled by a, not only infinite but continuous, spectrum of time-sectors.
We note en passant that the scaling (16) is ultrametric only if the scaling function vanishes for all α(t, t w ) > 1, for details see Appendix B. In fact, dynamic ultrametricity is a geometric property [22] that states that for all

C(t,t w )
t / t w t w = 2 8 t w = 2 10 t w = 2 12 t w = 2 14 t w = 2 16 triplet of times t 1 ≫ t 2 ≫ t 3 , one has in the limit t 3 → ∞: Finding dynamical ultrametricity in concrete models has been rather elusive up to now. An outstanding example is the critical trap model [24], where f (α > 1) = 0. It is amusing that the trap model is not ultrametric from the point of view of the equilibrium states [25]. Thus, the casual connections between static and dynamic ultrametricity are unclear to us. At any rate, since our scaling function in Fig. 10 does not show any tendency to vanish for α(t, t w ) > 1, we do not find compelling evidences for dynamic ultrametricity in this model. We have also looked directly to the plots of C(t 1 −t 2 , t 2 ) versus C(t 2 − t 3 , t 3 ) (see Appendix B) and we have not found convincing indications for the onset of dynamical ultrametricity. In this respect, it is worth to recall similarly inconclusive numerical investigations of the Sherrington-Kirkpatrick model [26]. There are two possible conclusions: 1. the model does not satisfy dynamical ultrametricity in spite of the fact that it satisfies (according to the standard wisdom) static ultrametricity. 2. Dynamical ultrametricity holds but its onset is terrible slow.
Both conclusions imply that it is rather difficult to use the dynamic experimental data (or any kind of data) to get conclusions on static ultrametricy. Of course it would be crucial to check if static ultrametricity is satisfied in this model, but this is beyond the scope of this paper.

B. Aging in C link
Just as in the 3D case [21], the aging dynamics in SG in the hypercube is a domain-growth process, see Fig 17. For any such process, the question of the ratio surfacevolume arises. When this ratio vanishes in the limit of large domain size, as it is the case for any RSB dynamics, one expects a linear relation between C link (t, t w ) and C 2 (t, t w ). This is precisely what we find in Fig. 11.

C. Thermoremanent magnetization
The experimental work indicates that for T < 0.9T c , the thermoremanent magnetization follows a power law with an exponent proportional to T c /T [27]. The data obtained in JANUS for a three dimensional SG (see Fig.  12 and [20]) agree with this statement. However, the data obtained in the hypercube model does not follow such power law, neither can them be rescaled with T log t.
This lack of an algebraic decay is surprising on the view of the exact results of Ref. [28]. Indeed, it was analytically shown there that, at T c , the thermoremanent magnetization of the SK model decays as t −5/4 . Universality strongly suggests that the same scaling should hold for our model. Although it seems not to be the case, at the first glance, Fig. 13-top, a closer inspection confirms our expectation. Indeed, when plotted as a function of t −5/4 , see inset in Fig. 13-top, the thermoremanent magnetization curve has a finite non-vanishing slope at the origin. As we show in bottom panel of Fig. 13, finite size effects do not contradict this claim. In summary, the magnetization decay for the hypercube suffers  from quite strong finite time effects, but asymptotically it scales with the proper exponent, at least at T c .

VII. NONEQUILIBRIUM CORRELATION FUNCTIONS AND FINITE SIZE EFFECTS
The importance of finite size effects in nonequilibrium dynamics has been emphasized recently [20,21]. In our case, we have encountered important size effects, both in C(t, t w ), Fig. 14, and in ξ(t w ), Fig. 17-top.
We compare in Fig. 15 the finite D effects in C(t, t w ) for two different MF models with fixed connectivity: the hypercube and a previously studied model (the random graph with connectivity z = 6, where each spin can interact with any other spin with uniform probability [29]). Clearly enough, the effects are much weaker in the hy-  percube model.
It is interesting to point out that, although the finite size effects seems to be important in C(t, t w ), they are largely absorbed when one eliminates the variable t in favor of C 2 (t, t w ), see Fig. 16. Hence, one of our main findings (the linear behavior of C link as function of C 2 ) seems not endangered by finite size effects.
A very clear finite size effect is in the coherence length, ξ(t w ). By definition, it cannot grow beyond D. Furthermore, what we find is that it hardly grows beyond D/2, Fig. 17-top. Nevertheless, at short times, we can identify a D-independent region, where it grows roughly as log t w . Hence, one is tempted to conclude that ξ D=∞ (t w ) ∝ log t w . At this point, finite size scaling suggests that both ξ D /D and log t w /D are dimensionless scaling variables. This is confirmed in Fig. 17-bottom, where a spectacular data collapse occurs. This is further confirmed by the Fourier transform ofĈ 4 (r), G(k). Note that, sinceĈ 4 (r) depends only on the length of the displacement vector r, also G(k) is rotationally invariant. Now, since k can range from 0 to D, it is clearly a dimensionless quantity (a dimesionful momentum would be p = k/D). It follows that G(k)/G(0) is a dimensionless quantity that may depend only on a dimensionless variable, such as log t w /D. Our data support this expectation, see Fig. 18. As for the k dependence of G(k), we expect a 1/p 4 behavior in the range of 1/ξ(t w ) ≪ p ≪ 1 [14] (note that we are in the q = 0 sector). Now, it is very important to recall that p 4 in Euclidean metrics translates into p 2 in the postman metrics. We have also seen that the dimensionful p (postman metrics) corresponds to k/D. Thus, since in our range of t w , ξ(t w ) ∼ log t w , the product G(k) p 2 + 1/ log 2 t w should be roughly constant as D grows. As we show in Fig. 19, the scaling is better for p of order 1 (k ∼ D), although it seems to improve for smaller p as D grows. As far as we know, this is the first observation of the p 4 propagator in a numerical work.

VIII. CONCLUSIONS
We have studied a spin glass model in the Ddimensional unit hypercube in the limit of large D, but with finite coordination number. We have shown that any short range model in such a lattice will behave as a mean field model in the thermodynamic limit (that coincides with the large D limit). An important advantage of this model is that it has a natural notion of spatial distance.
We have shown that any statistical mechanics model on the hypercube with random connectivity would be afflicted by huge finite size effects, for purely geometrical reasons. The obvious cure has consisted in restricting the connectivity graphs to those with a fixed number of neighbors. Unfortunately, constructing such graphs is far from trivial. We have generated a subset of them by means of a simple dynamic Monte Carlo. In this way, we obtain graphs that are isotropic. We have checked that the Edwards-Anderson model defined over these finite connectivity hypercubes verify some consistency checks, including comparison with the analytically computable correlation function in the paramagnetic phase.
We have numerically studied the nonequilibrium dynamics in the spin glass phase. The three main features found were: (i) aging dynamics consists in the growth of a coherence length, much as in 3D systems, (ii) the scal-ing of the two times correlation function implies infinitely many time-sectors, and (iii) the p 4 propagator has been observed. In addition, we have studied the finite size effects in our model, finding that a naive finite size scaling ansatz accounts for our data.

Acknowledgments
Computations have been carried out in PC clusters at BIFI and DFTI-UCM. We have been partly supported through Research Contracts No. FIS2006-08533 (MICINN, Spain) and UCM-BS, GR58/08. BS is an FPU fellow (Spain).

Appendix A: HIGH TEMPERATURE EXPANSION
For sake of clarity, we will firstly discuss the calculations for the random connectivity hypercube. Results for the fixed connectivity model will be then obtained by minor changes.
Using the identity (β = 1/T ) we can write the partition function and the spin propagator as: .
The high-temperature expansion (see, for instance [30]), expresses the propagator as a sum over lattice paths that join the points x and y, γ x→y : where l γ represents the length of the path γ x→y , J is the product of the couplings, J zw , along the path, and Z γ is a restricted partition function obtained by summing only over all closed paths that do not have any common link with the path γ x→y . However, when averaging over disorder, due to the randomness in the coupling signs, σ x σ y = 0. The spin glass propagator is obtained instead by averaging over disorder σ x σ y 2 . Clearly, the sum will be dominated by those diagrams where the go and return path are the same (thus, J 2 zw = 1): where K = tanh 2 β. In Bethe lattices, due to their cycleless nature, Z 2 γ /Z 2 = 1 in the thermodynamic limit. Hence, we are left with the problem of counting the average number of paths of length l γ that join x and y, p(l γ ). From it, we obtain The sum is restricted to l γ ≥ r because the length of the shortest path that joins x and y is given by their postman distance r.
In order to count the average number of paths, p(l γ ), let us distinguish two cases: l γ = r and l γ > r. The first will give the leading contribution in the large D limit.
The number of paths joining x and y in precisely r steps is r!, because the r steps are all taken along different directions and in a random order. For a given path, the probability of all the r links be active is (z/D) r . Hence Note that the D −r factor compensates exactly the divergence of the D r in Eq. (A6). In the case of l γ > r, one has l γ = r + 2k, with k > 0. Note that when l γ = r the path contains r different directions (namely, the Euclidean components in which x and y differ). Each of these directions appear only once. However, when l γ > r, other directions must be included, we call them unnecessary. Note that, if the path is to end at the desired point, any unnecessary step must be undone later on. Hence, l γ − r is always an even number 2k. Clearly, the number of such paths is bounded by Γ (r, k) D k , where Γ (r, k) is a D-independent amplitude. On the other hand, the probability of finding all the links active is (z/D) r+2k . Thus, we conclude that that results in a O D −k contribution toĈ 4 (r). Then, in the large D limit we obtain (A = zK): with finite size corrections of O D −1 . Thus, we encounter an exponential decay with an exponential correlation length given by Summing all up, we can compute the spin-glass susceptibility for the large D limit:  Figure  20: t3, t3)], t1 > t2 > t3 with t1 fixed by the condition C(t1 − t3, t3) = q and different t3. In the presence of dynamic ultrametricity, (B2), the parametric plot should tend for large t3 to the union of x = q and y = q. The panels correspond to q = 0.25 (top, nice BB scaling but not ultrametric), q = 0.35 (middle, nice BB scaling and ultrametric) and q = 0.5 (bottom, supposedly ultrametric but poor BB scaling).  Fig. 20, we take the intercept with x = y, i.e. C * = C(t1 − t2, t2) = C(t2 − t3, t3), and represent C * − q as a function of 1/ log t3. Lines: analogous plot for the toy model described in the text, where the BB scaling is exact.  where the scaling function f is smooth and monotonically decreasing. From now on, we shall refer to this scaling as BB scaling (after Bertin-Bouchaud).
Let us see under which conditions BB scaling implies the ultrametricity property where t 1 ≫ t 2 ≫ t 3 and t 3 tends to infinity.
The natural time dependency is a power law choice with µ 1 > µ 2 . In that case, the large t 3 limit for the argument of the scaling function are: α(t 1 − t 3 , t 3 ) = µ 1 , α(t 2 − t 3 , t 3 ) = µ 2 and α(t 1 − t 2 , t 2 ) = µ 1 if µ 2 < 1 and α(t 1 −t 2 , t 2 ) = µ 1 /µ 2 if µ 2 > 1. Then, the condition (B2) is only satisfied in case µ 2 < 1. If, as it is the case for the critical trap model [24], f (α > 1) = constant [37], the BB scaling would imply dynamic ultrametricity. This is not the case for a general scaling function f such as, for instance, the one we get in Fig. 10. Nevertheless, although this analysis implies that the dynamic ultrametricity is only present in our model in some range of parameters, let us try a more straight approach. We consider a fixed value for the correlation function, q. On the view of the previous considerations and of Fig. 10, we should expect ultrametricity only for q > f (α = 1) ≈ 0.35. Now, for each t 3 , we find t 1 such that C(t 1 − t 3 , t 3 ) = q. Then, we perform a parametric plot of C(t 1 − t 2 , t 2 ) vs. C(t 2 − t 3 , t 3 ), for t 3 < t 2 < t 1 . Ultrametricity predicts that, in the large t 3 limit, the curves should tend to a half square (e.g. the intersection of the straight lines x = q and y = q) and, in particular, when C(t 1 − t 2 , t 2 ) = C(t 2 − t 3 , t 3 ) = C * , C * should tend to q.
We present in Fig. 20 results for three different values of q: 0.5 (ultrametric region, but in our range of t w data do not scale according BB), 0.35 (ultrametric region and good BB scaling) and 0.25 (non ultrametric region but BB scaling works nicely). At the qualitative level, the parametric curves seem to tend to a corner, but the convergence is slow. Furthermore, there are no clear differences between the curves with q > f (α = 1) and those with q < f (α = 1). Hence, we may try a more quantitative analysis.
We obtain numerically C * , the point where C * = C(t 1 − t 2 , t 2 ) = C(t 2 − t 3 , t 3 ), and study C * − q as function of 1/ log t 3 in Fig. 21. This choice is due to the fact that in the ultrametric region BB scaling predicts Hence, we expect that C * − q will be of order 1/ log t 3 if ultrametricity holds. The numerical data confirms this picture only partly. For q = 0.35 the results are as expected, yet for q = 0.25 the difference is decreasing fast as t 3 grows and it is hard to tell whether the extrapolation will be zero or not. For q = 0.5 (where BB scaling is not working for our numerical data) the behavior is non monotonic.
To rationalize our finding, we consider a simplified model, where the BB scaling is supposed to hold exactly.
The master curve f (α) is taken from the numerical data for C(t, t 3 = 2 16 ) for D = 22. This toy model allows us consider ridiculously large values of t 3 . As we see in Fig. 21, the peculiarities of the master curve cause a non monotonic behavior in q for an ample range of t 3 .
The lack of monotonicity in q makes also on interest to focus on α, rather than on the correlation function. With this aim, we consider the time t 2 where C(t 1 − t 2 , t 2 ) = C(t 2 − t 3 , t 3 ) = C * , and compute 1 2 [α(t 1 − t 2 , t 2 ) + α(t 2 − t 3 , t 3 )] − α(t 1 − t 3 , t 3 ). BB scaling and ultrametricity combined, see Eq. (B5), imply that this quantity should be of order 1/ log t 3 (in the non ultrametric region, it should be of order one). Our results in Fig. 22 basically agree with these expectations.