# Subdiffusion and intermittent dynamic fluctuations in the aging regime of concentrated hard spheres

###### Abstract

We study the nonequilibrium aging dynamics in a system of quasi-hard spheres at large density by means of computer simulations. We find that, after a sudden quench to large density, the relaxation time initially increases exponentially with the age of the system. After a surprisingly large crossover time, the system enters the asymptotic aging regime characterized by a linear increase of the relaxation time with age. In this aging regime, single particle motion is strongly non-Fickian, with a mean-squared displacement increasing subdiffusively, associated to broad, non-Gaussian tails in the distribution of particle displacements. We find that the system ages through temporally intermittent relaxation events, and a detailed finite size analysis of these collective dynamic fluctuations reveals that these events are not spanning the entire system, but remain spatially localized.

## I Introduction

Aging refers to the slow evolution with time of the physical properties of a disordered material suddenly quenched into a glass phase young (). This might refer to the evolution of the density of a polymer glass, of the dielectric response in an organic liquid, or the height of a gently shaken sandpile. Aging is more easily detected by focusing on the dynamics of the system, for instance by measuring correlation or response functions. One can think of measuring the decay of density fluctuations using light scattering techniques in colloidal glasses, or the time-dependent magnetic response in a spin glass. After several decades of aging studies in glassy materials, these generic features are very well documented young (); leszouches (); lucareview ().

Much less is known, however, about the microscopic mechanisms involved during aging, and how these evolve with time, for two main reasons. First, the important theoretical developments from the last decade about aging dynamics mainly stemmed from ‘mean-field’ types of approaches, which provided detailed predictions about the behavior of averaged dynamic quantities, but very little information about microscopic motion young (). Second, it is only relatively recently that the microscopic dynamics, of, say, molecules in a supercooled liquid, has been characterized in detail in equilibrium conditions above the glass transition ediger (), and, by comparison, much less work has been done about the corresponding aging regime at lower temperatures.

In this paper, we use numerical simulations to study the aging dynamics of a model system designed to understand the behavior of dense suspensions of colloidal hard spheres. Hard spheres represent one of the simplest, and thus most studied, models to study the glass transition. It is easily studied in simulations, and the model finds its experimental realization using well-controlled colloidal suspensions naturepusey (); ericluca (). As opposed to molecular liquids, colloidal hard spheres can be studied at the particle scale using microscopy techniques weeks (); KegelScience2000 (), while dynamic light scattering provides a convenient way to characterize their dynamics in great detail martinez (); djamel (). It should be noted that while a relatively small set of experiments have been reported on the aging of colloidal hard spheres martinez (); djamel (); SimeonovaPRL2004 (); weeksaging (); weeksaging2 (), a much broader literature exists on the out-of-equilibrium dynamics of a variety of more complex colloidal systems. Throughout this paper, we will in particular make reference to findings for colloidal attractive gels lucaaging1 (); DuriEPL2006 (); DuriPRL2009 () and to systems with soft repulsive interactions, such as closely packed soft spheres yodh (); MazoyerPRL2006 (); MazoyerPRE2009 () or Laponite clay platelets interacting via Coulomb repulsion bonn (); kaloun (); bob1 (); joshi (); lequeux ().

In this work, we will be concerned with three main questions.

(1) How does the structural relaxation slow down with the aging time? While numerical simulations of model glasses generically find that the structural relaxation time increases roughly linearly (or sub-linearly) with the sample age kob-barrat (), some light scattering experiments on gels and Laponite report an unexpected exponential growth of the relaxation time with sample age lucaaging1 (); bonn (); kaloun (), at least at short times.

(2) How do particles move during the aging regime? Confocal microscopy experiments report that colloidal particles move very little in concentrated hard spheres, to the point that distinguishing thermal vibrations from genuine relaxation becomes an experimental challenge weeksaging (); weeksaging2 (). In simulations as well, particles move very little, leading to a very slow, typically algebraic, decay of time correlation functions kob-barrat (); kob-barrat2 (); puertas (). This is in stark contrast with several experimental reports of sharply decaying time correlation functions in aging molecular liquids seeded with colloidal particles or in colloidal gels and Laponite, where compressed exponential decay and ballistic particle motion over large distances were reported lucareview (); lucaaging1 (); bob1 (); CaronnaPRL2008 (); bob2 ().

(3) How heterogeneous is the dynamics? Due to the key role played by dynamic heterogeneity in equilibrium studies of the glass transition ediger (), similar signatures have been sought in aging materials. Optical and confocal microscopy studies revealed the existence of rare small-scale relaxation events involving extended clusters of particles weeks (), whose size is modest and grows very little yodh (), or not at all weeksaging (); weeksaging2 (), during aging. This agrees with numerical simulations on model glasses, where four-point dynamic susceptibility appear to grow at a very slow rate with the sample age parisi (); castillo (), suggesting that if collective displacements occur, they are correlated over relatively short length scales. This set of results is however in contrast with another series of observations. Sudden relaxation events termed (rather dramatically) ‘earthquakes’ kob-barrat2 () or ‘avalanches’ lee1 (); lee2 () were reported in numerical simulations of Lennard-Jones glasses and these were suggested to be spanning the entire system. Similarly, highly intermittent dynamic fluctuations were experimentally recorded in several aging systems from polymer glasses buisson () to colloidal gels lucaaging1 (); DuriEPL2006 (); lequeux () and soft spheres MazoyerPRL2006 (), the latter typically involving ballistic particle motions correlated over very large distances DuriEPL2006 (); DuriPRL2009 (); MazoyerPRL2006 (); MazoyerPRE2009 ().

The paper is organized as follows. In Sec. II, we present the numerical model and techniques used in the present study. In Sec. III, we describe the average aging dynamics for the evolution of the energy and relaxation timescale. In Sec. IV, we focus on the subdiffusive and heterogeneous dynamics at the single particle level. In Sec. V we present results concerning the intermittent collective dynamics of the system. In Sec. VI we discuss our results and conclude the paper with some perspective for future research.

## Ii Numerical model and techniques

We perform molecular dynamics simulation of dense systems of strongly repulsive particles interacting with a very steep pair potential designed to model the behavior of colloidal hard spheres Voigtmann ():

(1) |

where is an energy scale, , being the position of particle , and , where represents the diameter of particle . To prevent crystallization occurring in dense systems of hard spheres, we introduce a size polydispersity and draw the particle diameters from a flat size distribution, , so that the average diameter is , and the polydispersity is given by

(2) |

To study the dynamics of the system we solve Newton’s equations for a system composed of particles, using a velocity Verlet algorithm allen () in a cubic box of linear size . We use periodic boundary conditions in the three directions of space. We have studied two system sizes, and . When dealing with averages, we shall report results for the largest system size, while a comparison between the two systems will allow us to perform a detailed comparison of the dynamic fluctuations occurring during the aging dynamics as a function of system size. We perform simulations in the ensemble, and we control the temperature by rescaling velocities every 100 molecular dynamics timesteps in order to maintain the kinetic energy to the desired constant value.

For the inverse power law potential in Eq. (1), density and temperature cannot be controlled independently, as rescaling the density by a factor is equivalent to rescaling the temperature by a factor . Thus we fix the temperature and energy scales, , and simply vary the density, , of the system Voigtmann (). To ease the comparison with hard sphere experiments, we express our results using the volume fraction, , rather than density , where

(3) |

In the following, we express length scales in units of , and timescales in units of , where is the mass of the particles. We use a time discretization , which ensures a proper integration of the equations of motion.

We shall study the aging dynamics of samples in the range . To obtain reproducible results, we need to produce fully disordered initial states. To this end, we first equilibrate the system at very low volume fraction, . We then compress the system very rapidly to the desired final volume fraction. The compression is done in small successive steps in order to avoid large overlaps leading to very large repulsive forces. The age of the system is counted from the time when the system reaches the final volume fraction. To increase the statistical significance of our results, we have performed 5 independent runs at each volume fraction with particles, starting from independent configurations compressed from the fluid at . During the course of the simulations we found no sign of incipient crystallization in the system, while crystallization was a major obstacle in an initial set of studies using a smaller polydispersity, , for which we do not present results.

We have previously studied the equilibrium dynamics of this system djamel (). We found that the dynamics slows down considerably when increases above . Fitting the increase of the relaxation to a power-law divergence at some critical volume fraction , we obtained . This fit then locates the (apparent) mode-coupling singularity for this system brambilla (), which should serve as a reference volume fraction for the aging studies below, since it becomes very difficult to reach thermal equilibrium within the accessible numerical timescales when increases beyond .

## Iii Towards the asymptotic aging regime

We now present the results of our study, starting from the time evolution of the simplest quantity we can monitor in our study, the potential energy density, defined as

(4) |

where the brackets stand for an average over independent initial configurations. We present representative results for the time evolution of at large densities in Fig. 1.

The central observation from this figure is that, at these large densities, the potential energy keeps evolving over the 6 decades of the simulation and never reaches an asymptotic plateau. This simply confirms that thermal equilibrium cannot be reached at these volume fractions, because the equilibrium relaxation time is much larger than the maximum time accessible to our simulations. More in detail, we also observe that the energy evolves initially quite rapidly, decreasing by a factor of about 2 when time increases from to , while it evolves only by a few percent when time further increases by two additional decades. Thus, we clearly see the effect of ‘aging’, since dynamical evolution slows down considerably as the age of the system increases, as is well-known from decades of experimental aging studies in many different materials.

To describe the slow evolution of the potential energy, we fitted our data to both a power law decay, , and a logarithmically slow evolution, . We find that both fits give an equally good description of the final, slow evolution of the potential energy. The fit parameters, and confirm that the evolution of the energy is indeed extremely slow.

We also find that the fits only hold when is very large: the first few decades of the simulation, corresponding to a faster evolution of are not well described by this asymptotic slow decay. This implies that a relatively long time is needed for the system to enter the asymptotic aging regime. Moreover, we find that the time required to enter the asymptotic regime increases when density increases and the system is quenched deeper into the glassy state. This implies that considerable care must be taken when performing data analysis of the aging regime, since it already takes a large part of the simulation simply to enter the final regime, indicated by the arrows in Fig. 1. It is likely that the first regime corresponds to faster processes where the largest overlaps created during the compression are removed, producing heat which is then removed by our thermostatting procedure. Thus, no ‘universal’ characteristics are to be expected in this time regime, which might well depend quite strongly on the details of the preparation procedure, or the chosen thermostat. Note that for the data presented in Fig. 1, meaning that scaling behavior might only be observed in the demanding limit of

(5) |

It is crucial to recognize the existence of two distinct aging regimes in order to analyze properly the scaling properties of dynamic functions in the asymptotic aging regime, as we find that these two distinct regimes in fact affect most of the measurements we made in our simulations. In Fig. 2 we show the behavior of the self-intermediate scattering function following a quench at . This two-time quantity is defined as

(6) |

and we perform measurements at , close to the first peak in the static structure factor. As is well-known, two-time quantities reveal the non-stationary evolution of the system in a very direct manner. In agreement with previous work kob-barrat (); puertas (), we find that the decay with the delay time of occurs in a two-step manner, the slow decay being a strong growing function of the waiting time . The faster initial decay is much less sensitive to the waiting time, and corresponds physically to particle vibrations in a nearly frozen amorphous structure.

The existence of two distinct aging regimes is obvious from these data, since the slow decay of clearly shows a crossover for time delay such that , as shown in Fig. 2, while the data become simpler to describe when as the entire decay then takes place in the asymptotic long-time regime defined by Eq. (5). As found in other systems kob-barrat2 (), we find that the long-time decay of the self-intermediate scattering function is well described, in the aging regime, by a power law, , as indicated in Fig. 2 by the full lines. As noted before, such an algebraic decay is in contrast with equilibrium measurements at lower volume fraction which typically display stretched exponential relaxations. Stretched exponential decays are only seen for shallow quenches when a crossover towards thermal equilibrium is possible. The small value of the exponent, means that relaxation occurs over a very broad time window, and that it is not possible to define a ‘typical’ relaxation time for the system.

Nevertheless, to quantify the aging of the dynamics, we follow the common practice and define an -relaxation time from the time decay of as , as indicated by the horizontal dashed line in Fig. 2. In Fig. 3, we report the evolution of with waiting time for all volume fractions investigated in this work, from to . Again, we observe two distinct regimes for the evolution of . While for , seems to increase exponentially with waiting time,

(7) |

with a numerical constant, its growth is better described by an algebraic law at long times, , which is well compatible with a so-called ‘simple aging’ behavior young (),

(8) |

This simple aging behavior is shown with full lines in Fig. 3. In agreement with the observations in Fig. 1, we find that it takes an increasingly long time for the system to enter the asymptotic (simple) aging regime when increases. For the larger studied, , the simple aging regime is not reached during the course of our simulations, and the system appears to be in the crossover between (7) and (8), and the system seems to undergo an effective ‘super-aging’, i.e. its relaxation time increases faster than its age. Finally, the opposite ‘sub-aging’ behavior is found when volume fraction is not too large, for instance , because the system is crossing over towards thermal equilibrium, such that should saturate at long waiting times to the equilibrium relaxation time. These observations show that a complex behavior of , as often reported in numerical works kob-barrat (); lee3 () and experiments on colloidal gels and Laponite lucaaging1 (); bonn (); kaloun (), may be due to the occurrence of multiple crossovers which are highly sensitive to volume fraction.

These observations are also in agreement with the common observation of sub-aging behavior in aging molecular liquids, for which experiments are traditionally performed not very far below the glass temperature young (); struik (). Exponential growth of the relaxation time sometimes (but not always) followed by algebraic growth, has been reported in a few experiments on colloidal glasses or gels as well lucaaging1 (); bonn (); kaloun (). Our results thus suggest that this exponential growth might well be a transient behavior, which can persist, however, over a very large time window, in particular for very deep quenches. Our results also highlight the difficulty in analyzing the scaling properties of two-time dynamic quantities in numerical simulations, since the asymptotic aging regime is only accessed after a time which can be very large, thus decreasing the window where ‘universal’ aging properties can be probed. This might well explain some conflicting results and analysis reported in the literature regarding, e.g., the proper scaling behavior of the self intermediate scattering function in Lennard-Jones glasses kob-barrat3 (); rieger ().

## Iv Subdiffusive and heterogeneous dynamics

The above findings that the energy decays very slowly, while time correlation functions decay in an algebraic manner indicate that there is actually very little dynamics taking place in the system. To confirm this idea, we measured the averaged mean-squared displacement,

(9) |

In Fig. 4, we present the time evolution of the mean-squared displacements at volume fraction measured at different waiting times, using both a standard log-log representation (inset) and a lin-log representation (main plot). As found in Fig. 2 for the self-intermediate scattering function, the dynamics proceeds again in a two-step fashion with a rapid, ballistic regime at short-time, corresponding to vibrations of the particles in a frozen amorphous structure, followed by a much slower, waiting time dependent structural relaxation, which becomes slower when increases. This implies that the aging of the system corresponds, at the microscopic scale, to a dramatic slowing down of single particle displacements, as found in experiments weeksaging (); yodh (); SimeonovaPRL2004 ().

It is significant that in order to represent graphically the behavior of the mean-squared displacements, we had to use a linear scale for its amplitude, with a range which remains smaller than . This implies that over the entire duration of the simulation, particles, on average, move a distance which is smaller than the mean particle diameter: despite the non-trivial aging dynamics we discuss, it should be clear that the structure of the material remains in fact almost frozen as soon as the system is quenched in the glassy phase, with essentially no large-scale dynamics taking place.

The linear representation in Fig. 4 also confirms the existence of the two distinct aging regimes discussed above, thus directly revealing the influence on the dynamics at the particle scale of the timescale . Again, the data obtained for present a crossover for a delay given by , as seen in Fig. 4. Thus, in the following we again focus on waiting times that are larger than , which, we believe, reflect better the ‘universal’, quench-independent nature of the aging regime in concentrated hard spheres.

In the long time regime, the dashed line in Fig. 4, which corresponds to a diffusive growth, , is obviously an incorrect description of our data since the displacements seem to increase much more slowly than diffusively. As shown in the inset, a subdiffusive law describes our results very satisfactorily,

(10) |

where is an exponent characterizing the subdiffusion.

We have measured over a broad range of waiting times and volume fractions, and we present its evolution in Fig. 5. Apart from the low volume fractions where (diffusive) thermal equilibrium is reached rapidly, we systematically find that the system obeys subdiffusive rather than diffusive behavior. At a given volume fraction, the system gets closer to equilibrium for larger , and correspondingly we find that the exponent increases with the age of the system, although it always remains very far from its equilibrium, diffusive value . Note that an increasing does not imply that particles move faster at large waiting times, as the prefactor in the power law (10) is itself a decreasing function of , which implies that the total amplitude of the particle displacements indeed decreases with the age of the system.

As volume fraction increases, the system is quenched deeper and deeper into its glass phase, and deviations from diffusive behavior are correspondingly larger, which translates into a subdiffusive exponent which gets smaller when increases, see Fig. 5. This simply means that particles move less and less when becomes larger, which is physically expected.

While the mean-squared displacement quantifies the average dynamical behavior of the system, it does not convey much information on dynamic fluctuations, which are recognized as an important feature of the single particle dynamics in glassy materials. To quantify this ‘dynamical heterogeneity’ during aging, we measure the probability distribution function of the single particle displacement, also known as the self-part of the van-Hove correlation function:

(11) |

where represents the projection of along the -axis. To improve the statistics, we use the isotropy of the system and further average along the three directions of space. For simplicity, we keep the notation for the van-Hove function averaged over the , , and directions.

To gain deeper insight into the subdiffusive behavior described above, we present in Fig. 6 the typical shape and evolution of the van-Hove function measured for , for a large waiting time, in the asymptotic aging regime, and a delay within the long-time subdiffusive regime. The data are presented in a lin-log scale, where Fickian behavior leading to a Gaussian shape of the particle distribution would appear as an inverted parabola, as shown by the dashed lines. These data are highly reminiscent of the equilibrium findings that van-Hove functions are well-described at short displacements by a Gaussian form, with tails that are much broader than the Gaussian prediction weeks (); KegelScience2000 (). This implies that a small fraction of the particles move significantly farther than what is expected from a purely diffusive and homogeneous process. This is the most basic observation that the system is dynamically heterogeneous and can be described, as being composed of distinct families of ‘slow’ and ‘fast’ particles, a well-established distinctive feature of many disordered, glassy materials ediger ().

How should we describe the functional form of the tails of these distributions? At thermal equilibrium, it is now understood that the tails are well-described by an exponential (rather than Gaussian) decay, , for corresponding to the alpha-relaxation timescale pinaki (). The physical origin of the exponential, detailed in Refs. pinaki (); epl (); pinaki2 (), is due to the stochastic nature of the diffusion process in disordered materials, which is well described by the formalism of continuous time random walks (CTRW) MW (). In this description, particles in a dense supercooled liquid undergo a succession of long periods of vibrations within the cages formed by their neighbors, followed by rapid ‘jumps’. Both the size of the jumps and, more importantly, the timescale separating them are statistically distributed quantities. The existence of these statistical distributions directly accounts for the exponential form of the tails in the van-Hove functions pinaki ().

The natural extension of these considerations to the non-equilibrium aging regime studied in the present work is the aging continuous time random walk (ACTRW) formalism actrw (), whose main features are those of the CTRW recalled above. The only difference lies in the functional form used for the distribution of times separating the jumps, which acquires ‘fat’, non-normalizable tails in the aging regime, in direct analogy with the trap model introduced by Bouchaud to describe aging dynamics in glasses jptrap (). In this approach, the tails of the jump time distribution decay as , where is the subdiffusion exponent introduced in Eq. (10), while the tails of the self-part of the van-Hove function are not exponential anymore, but are instead asymptotically described by metzler ()

(12) |

where the exponent is also connected to via the relation

(13) |

so that (Gaussian) is recovered when (Fickian diffusion). For subdiffusive processes, , one expects instead .

To test in detail the ACTRW picture, one should measure the distribution of jump times , and use it to directly predict the mean-squared displacements and van-Hove functions rottler (). This is an interesting project, but is beyond the scope of the present study. Instead, we more simply use Eq. (12) as a theoretically motivated fitting formula for the tails of the distributions. As shown in Fig. 6, such a fit describes the tails very well. We find a similarly good description for a broad range of volume fractions and timescales. Unfortunately our data is not accurate enough to allow for a very precise determination of all the fitting parameters. Thus, the following results should be discussed at a qualitative, rather than quantitative, level.

In Fig. 7 we present the evolution of the exponent describing the tails of the van-Hove functions in the asymptotic subdiffusive aging regime. We remark that these data have considerably more scatter than the data for the subdiffusive exponent , confirming the difficulty to obtain a reliable determination of from our data. Regardless of the noise, we note that values are in the interval , as expected from Eq. (13), suggesting that our description is physically sound. Moreover, the data at large density present a systematic evolution of , which increases towards 2 with . Assuming that Eq. (13) is correct, this would imply that increases towards unity with , as indeed is observed in Fig. 5. The data at moderate volume fractions, , are more difficult to interpret as seems to decrease with in this regime, while was found to increase. This could be due to the fact that for , the subdiffusive aging regime is crossing over towards equilibrium during the simulation, and so the results could be a subtle combination of the exponential tail reported at intermediate times for equilibrium dynamics pinaki (), and the subdiffusive regime found deeper in the glass and described in the present study.

## V Intermittent dynamic fluctuations

Our study of single particle dynamics suggests that particles undergo very little dynamics, leading to a subdiffusive growth of mean-squared displacements, associated to a broad distribution of particle displacements. The picture of ACTRW which we used to describe our data is very similar in spirit to Bouchaud’s trap model, and both suggest that the aging dynamics of concentrated hard spheres is a temporally intermittent process monthus (). It is the aim of this section to establish whether this picture is indeed correct.

To this end, we must turn to collective dynamic observables, and resolve the aging dynamics in space and time, and focus on , the instantaneous (unaveraged) value of . We present in Fig. 8(a) the results of independent realizations of the dynamics at a given waiting time, , and volume fraction, . For this system, containing particles, we observe small run-to-run fluctuations, which show that resolving the temporal evolution of the dynamics in a system of linear size is not sufficient to reveal significant dynamic fluctuations.

Thus, we improve our spatial resolution and repeat this analysis for a smaller system with which contains particles, see Fig. 8(b). The data in Fig. 8 show that run-to-run fluctuations of become larger when is smaller, as expected. However, we emphasize that the nature of the dynamic fluctuations seems to change qualitatively when the size is reduced. For a single realization of a quench with a small enough number of particles, the decay of the self-intermediate scattering function in the long-time regime is highly intermittent, as shown by sudden ‘drops’ of kob-barrat2 (). From one run to another, both the temporal location and the amplitude of these drops vary considerably, indicating that in each run the dynamics proceeds via temporally intermittent relaxation events that mobilize a finite fraction of the whole system. However, the comparison with the fluctuations obtained with suggests that the spatial extension of these intermittent events does not increase with the system size, but remains instead spatially localized.

Note that this finding could be due to a spurious finite size effect. An alternative interpretation of the data in Fig. 8 could be that sudden decorrelations in the system only occur when the system is too small, and these events are not observed in a bigger system because they disappear altogether. To disprove this interpretation, we reanalyze a single quench with particles, and further decompose the computation of the correlation function in 8 distinct sub-systems, each comprising 500 particles and corresponding to a cubic sub-box of linear size . In Fig. 9, we present the result of this analysis. Remarkably, we find that in a single quench, a small part of the system might indeed undergo the type of sudden rearrangement seen in particles runs, but decorrelation events are not necessarily seen in other subsystems. This is visible in the top graph, where representative for various sub-boxes are shown, and in the bottom snapshot, where the most mobile particles between the two times shown by the arrows in the main plot are depicted as large spheres. Clearly, individual rearrangement events are confined to a fraction of the total volume (in that case a corner of the box), thus affecting significantly only the correlation function of the corresponding sub-box. We conclude that in a large system the aging dynamics occurs through an intermittent succession of spatially localized decorrelation events, with no obvious ‘confinement’ effect imposed by the use of too small a system size, at least for the two system sizes used in this study.

To confirm quantitatively both the absence of serious finite size effects, and the spatially localized nature of decorrelation events, we have measured the variance of the spontaneous fluctuations of the self-intermediate scattering function, also known as a ‘four-point dynamic susceptibility parisi (); castillo ():

(14) |

where represents the instantaneous value of . With this normalization, should become independent of the system size in the thermodynamic limit, with an amplitude that gives direct access to the typical number of particles relaxing in a correlated manner during aging TWBBB (); dalle (). The data for and shown in Fig. 10 are essentially the same, within our statistical accuracy, for both system sizes, and reaches a value near for timescales corresponding to the slow relaxation of the correlation function, confirming the spatially localized nature of relaxation events.

## Vi Discussion

We have studied numerically the aging dynamics of a concentrated system of nearly hard spheres over a broad range of volume fractions, and a large time window. We have addressed the three central questions presented in the introduction, concerning the evolution of the dynamics, single particle motion, and the heterogeneous nature of the dynamics.

We found that over the 6 decades of aging times covered by our simulations, the structural relaxation continuously slows down, as expected for any aging system, but we showed that an asymptotic aging regime is only accessed after a timescale , which can become quite large for dense systems. Interestingly, we found that the relaxation time first increases exponentially, as observed in some experiments on colloidal gels and repulsive Laponite platelets lucaaging1 (); bonn (); kaloun (), before crossing over to a simple aging form, , at large times, . In the initial regime, we also found that the energy density decays very quickly, suggesting that the system evolves rapidly from its highly disordered initial state to form a nearly frozen structure, which then ages more slowly. Thus, this initial regime is likely quite sensitive to the detailed initialization procedure of the system both in simulations and experiments. The ‘universal’ features of the aging dynamics of concentrated hard spheres should be discussed for the second regime only, and our simulations indicate that concentrated hard spheres follow a simple aging form, as indeed found for very many glassy materials young ().

In the aging regime, we found that particle motions are very restrained, with particles moving on average less than their own diameter over the entire duration of simulation. In the aging regime, , we found that the self-intermediate functions decay algebraically at long times, in agreement with several simulations of soft and hard particle systems kob-barrat (); puertas (). This is also in agreement with recent experiments performed on colloidal suspensions that are dense enough, so that the aging regime does not cross over at long times towards equilibrium behavior djamel (); martinez ().

We also found that single particle motion is sub-diffusive, , with . This behavior is actually quite consistent with the results obtained by optical and confocal microsopy on hard and soft colloids, where diffusive behavior is not reached in the experimental time window weeksaging (); weeksaging2 (); yodh (). This result is also in broad agreement with earlier numerical work kob-barrat (), although subdiffusion was never described in detail before. We have suggested that a natural theoretical framework to analyze our results should the Aging Continuous Time Random Walk (ACTRW), i.e. the nonequilibrium extension of the random walk picture used to describe single particle motion at thermal equilibrium near the glass transition pinaki (); epl (). This formalism makes a number of detailed predictions for the aging regime. Here, for lack of statistics, we have simply been able to show that our data are in qualitatively good agreement with this picture, which can link the shape of the distribution of particle displacements to the subdiffusion exponent. It would be extremely interesting to extend this analysis in future work and check in more detail how far the ACTRW picture can be pushed to describe the aging dynamics of concentrated hard spheres.

We finally showed, by resolving the measurement of time correlation functions in space and time, that the aging dynamics of concentrated hard spheres is highly intermittent, with very sudden relaxation events separating long periods of time where very little dynamics occurs. Similar dynamic events have been observed in numerical work before, and were coined ‘earthquakes’ kob-barrat2 () or ‘avalanches’ lee2 (). Although these names and additional numerical evidence suggested that these events could well be system-spanning lee2 (), we showed by using much larger system sizes that these events remained actually spatially localized and do not grow with system size. This is in fact quite consistent with the experimental report of dynamic clusters that grow rather modestly in aging colloidal samples weeksaging2 (); yodh (), and with numerical work reporting a slow growth of four-point dynamic susceptibilities in aging Lennard-Jones glasses castillo ().

Thus, we find that in concentrated hard spheres the dynamics are intermittent as it occurs in several more complex colloidal materials, but this leads neither to anomalous compressed exponential relaxation for time correlation functions lucaaging1 (), nor to ballistic motion over large distances lucaaging1 (); DuriPRL2009 (); MazoyerPRL2006 (); MazoyerPRE2009 (); maccarrone (), whose origin thus remains largely mysterious. This means that model systems such as hard spheres or Lennard-Jones glasses are not good starting points to gain insight into the nature of these intriguing aging dynamics. It is not clear what ingredient these models are missing, since similar anomalous aging dynamics was recently reported in simple ‘hard’ molecular glasses approaching the glass transition CaronnaPRL2008 (); bob2 (). It would be very interesting to discover a model system displaying the same type of anomalous aging dynamics, with a particle-scale dynamics that can be followed in the simulations in the way similar to what we did for hard spheres.

###### Acknowledgements.

This work was financed by Grants NWO-SRON, ANR Dynhet, Région Languedoc-Roussillon ‘Chercheurs d’avenir’, and ACI Jeunes Chercheurs.## References

- (1) Spin glasses and random fields, Ed.: A. P. Young (World Scientific, Singapore, 1998).
- (2) Slow relaxations and nonequilibrium dynamics in condensed matter, Eds: J.-L. Barrat, J. Dalibard, M. Feigelman, J. Kurchan (Springer, Berlin, 2003).
- (3) L. Cipelletti and L. Ramos, J. Phys.: Condens. Matter 17, R253 (2005).
- (4) M. D. Ediger, Annu. Rev. Phys. Chem. 51, 99 (2000).
- (5) P. N. Pusey and W. van Megen, Nature 320, 340 (1986).
- (6) E. Weeks and L. Cipelletti, to be published.
- (7) E. Weeks, J. C. Crocker, A. C. Levitt, A. Schofield, and D. A. Weitz, Science 287, 627 (2000).
- (8) W. K. Kegel and A. van Blaaderen, Science, 287, 290, 2000.
- (9) V. A. Martinez, G. Bryant, and W. van Megen, Phys. Rev. Lett. 101, 135702 (2008).
- (10) D. El Masri, G. Brambilla, M. Pierno, G. Petekidis, A. Schofield, L. Berthier, and L. Cipelletti, J. Stat. Mech. P07015 (2009).
- (11) N. B. Simeonova and W. K. Kegel, Phys. Rev. Lett. 93, 035701 (2004).
- (12) R. E. Courtland and E. R. Weeks, J. Phys.: Condens. Matter 15, S359 (2003).
- (13) J. M. Lynch, G. C. Cianci, and E. R. Weeks, Phys. Rev. E 78, 031410 (2008).
- (14) L. Cipelletti, S. Manley, R. C. Ball, and D. A. Weitz, Phys. Rev. Lett. 84, 2275 (2000).
- (15) A. Duri and L. Cipelletti, Europhys. Lett., 76, 972, 2006.
- (16) A. Duri, D. A. Sessoms, V. Trappe and L. Cipelletti, Phys. Rev. Lett., 102, 085702, 2009.
- (17) P. Yunker, Z. Zhang, K. B. Aptowicz, and A. G. Yodh, Phys. Rev. Lett. 103, 115701 (2009).
- (18) S. Mazoyer, L. Cipelletti and L. Ramos, Phys. Rev. Lett., 97, 238301, 2006.
- (19) S. Mazoyer, L. Cipelletti and L. Ramos, Phys. Rev. E, 79, 011501, 2009.
- (20) B. Abou, D. Bonn, and J. Meunier, Phys. Rev. E 64, 021510 (2001).
- (21) S. Kaloun, R. Skouri, M. Skouri, J. P. Munch, and F. Schosseler, Phys, Rev. E 72, 011403 (2005).
- (22) R. Bandyopadhyay, D. Liang, H. Yardimci, D. A. Sessoms, M. A. Borthwick, S. G. J. Mochrie, J. L. Harden and R. L. Leheny, Phys. Rev. Lett. 93, 228302 (2004).
- (23) Y. M. Joshi and G. R. K. Reddy, Phys. Rev. E 77, 021501 (2008).
- (24) A. Mamane, C. Fretigny, F. Lequeux, and L. Talini, EPL 5, 58002 (2009).
- (25) W. Kob and J.-L. Barrat, Phys. Rev. Lett. 78, 4581 (1997).
- (26) W. Kob and J.-L. Barrat, Eur. Phys. J. B 13, 319 (2000).
- (27) A. M. Puertas, J. Phys.: Condens. Matter 22, 104121 (2010).
- (28) C. Caronna, Y. Chushkin, A. Madsen and A. Cupane, Phys. Rev. Lett., 1, 055702, 2008.
- (29) H. Y. Guo, G. Bourret, M. K. Corbierre, S. Rucareanu, R. B. Lenox, K. Laaziri, L. Piche, M. Sutton, J. L. Harden, and R. L. Leheny, Phys. Rev. Lett. 102, 075702 (2009).
- (30) G. Parisi, J. Chem. Phys. B 103, 4128 (1999).
- (31) A. Parsaeian and H. E. Castillo, Phys. Rev. E 78, 060105(R) (2008).
- (32) K. Vollmayer-Lee, W. Kob, K. Binder and A. Zippelius, J. Chem. Phys. 116, 5158 (2002).
- (33) K. Vollmayer-Lee and E. A. Baker, Europhys. Lett. 76, 1130 (2006).
- (34) L. Buisson, L. Bellon, and S. Ciliberto, J. Phys.: Condens. Matter 15, S1163 (2003).
- (35) T. Voigtmann, A. M. Puertas, and M. Fuchs, Phys. Rev. E 70, 061506 (2004).
- (36) M. Allen and D. Tildesley, Computer Simulation of Liquids (Oxford University Press, Oxford, 1987).
- (37) G. Brambilla, D. El Masri, M. Pierno, G. Petekidis, A. B. Schofield, L. Berthier, and L. Cipelletti, Phys. Rev. Lett. 102, 085703 (2009).
- (38) K. Vollmayr-Lee, J. A. Roman, and J. Horbach, arXiv:1001.0715.
- (39) A. S. Negi and C. O. Osuji, Phys. Rev. E 80, 010404 (2009).
- (40) L. C. E. Struik, Physical aging in amorphous polymers and other materials (Elsevier, Amsterdam, 1978).
- (41) U. Müssel and H. Rieger, Phys. Rev. Lett. 81, 930 (1998).
- (42) W. Kob and J.-L. Barrat, Phys. Rev. Lett. 81, 931 (1998).
- (43) P. Chaudhuri, L. Berthier, and W. Kob, Phys. Rev. Lett. 99, 060604 (2007).
- (44) L. Berthier, D. Chandler, and J. P. Garrahan, Europhys. Lett. 69, 320 (2005).
- (45) P. Chaudhuri, Y. Gao, L. Berthier, M. Kilfoil, and W. Kob, J. Phys.: Condens. Matter 20, 244126 (2008).
- (46) E. W. Montroll and G. H. Weiss, J. Math. Phys. 6, 167 (1965).
- (47) E. Barkai and Y. Cheng, J. Chem. Phys. 118, 6167 (2003).
- (48) J. P. Bouchaud, J. Phys. I France 2, 1705 (1992).
- (49) R. Metzler and J. Klafter, Phys. Rep. 339, 1 (2000).
- (50) M. Warren and J. Rottler, EPL 88, 58005 (2009).
- (51) C. Monthus and J. P. Bouchaud, J. Phys. A 29, 3847 (1996).
- (52) C. Toninelli, M. Wyart, L. Berthier, G. Biroli, and J.-P. Bouchaud, Phys. Rev. E 71, 041505 (2005).
- (53) Dalle-Ferrier C. Dalle-Ferrier, C. Thibierge, C. Alba-Simionesco, L. Berthier, G. Biroli, J.-P. Bouchaud, F. Ladieu, D. L’Hôte, and G. Tarjus Phys. Rev. E 76, 041510 (2007).
- (54) S. Maccarrone et al., to appear in Soft Matter.