Elastically driven intermittent microscopic dynamics in soft solids

Soft solids with tunable mechanical response are at the core of new material technologies, but a crucial limit for applications is their progressive aging over time, which dramatically affects their functionalities. The generally accepted paradigm is that such aging is gradual and its origin is in slower than exponential microscopic dynamics, akin to the ones in supercooled liquids or glasses. Nevertheless, time- and space-resolved measurements have provided contrasting evidence: dynamics faster than exponential, intermittency and abrupt structural changes. Here we use 3D computer simulations of a microscopic model to reveal that the timescales governing stress relaxation, respectively, through thermal fluctuations and elastic recovery are key for the aging dynamics. When thermal fluctuations are too weak, stress heterogeneities frozen-in upon solidification can still partially relax through elastically driven fluctuations. Such fluctuations are intermittent, because of strong correlations that persist over the timescale of experiments or simulations, leading to faster than exponential dynamics.

S oft condensed matter (proteins, colloids or polymers) easily self-assembles into gels or amorphous soft solids with diverse structure and mechanics [1][2][3][4][5] . The nanoscale size of particles or aggregates makes these solids sensitive to thermal fluctuations, with spontaneous and thermally activated processes leading to rich microscopic dynamics. Besides affecting the time evolution, or aging, of the material properties at rest, such dynamics interplay with an imposed deformation and are hence crucial also for the mechanical response [6][7][8][9][10] . In the material at rest, thermal fluctuations can trigger ruptures of parts of the microscopic structure that are under tension or can promote coarsening and compaction of initially loosely packed domains, depending on the conditions under which the material initially solidified. The diffusion of particles or aggregates following such local reorganizations (often referred to as 'micro-collapses') is governed by a wide distribution of relaxation times, due to the disordered tortuous environment in which it takes place, that is, the microstructure of a soft solid, be it a thin fractal gel or a densely packed emulsion 11,12 . Therefore, localization, caging and larger-scale cooperative rearrangements in such complex environment will lead to slow, subdiffusive dynamics, similar to the microscopic dynamics close to a glass transition 13,14 .
Several quasi-elastic scattering experiments or numerical simulations that can access micro-and nanoscale rearrangements in soft matter have indeed confirmed that the dynamics in aging soft solids are slower than exponential. Measuring the time decay of the correlations of the density fluctuations or of the local displacements of particles or aggregates, the experiments find that time correlations follow a stretched exponential decay e À ðt=tÞ b with bo1, akin to the slow cooperative dynamics of supercooled liquids and glasses, in a wide range of gels and other soft solids [15][16][17][18][19] . Nevertheless, in the past few years there has been emerging evidence, through a growing body of similar experiments, that microscopic dynamics in the aging of soft materials can be instead faster than exponential. In a wide range of soft amorphous solids including colloidal gels, biopolymer networks, foams and densely packed microgels, time correlations measured via quasi-elastic scattering or other time-resolved spectroscopy techniques appear to decay as e À ðt=tÞ b with b41 (and in most cases 1.3rbr1.5) [20][21][22][23][24][25][26][27][28][29][30] . When analysing the dependence of the relaxation time on the lengthscales probed, the data highlight the presence of super-diffusive rather than subdiffusive microscopic processes, with the relaxation time t increasing with decreasing the scattering wave vector q (and hence increasing the lengthscale being probed) as tp1/q rather than tp1/q 2 . Moreover, experiments specifically designed to quantify the fluctuations in time and space of the microscopic dynamics have revealed strong intermittency and pointed to the possibility of an abrupt and coherent reorganization of relatively large domains. Hence it has been suggested that the elastic rebound of parts of the material, after micro-collapses occur in its structure, should be, instead, the stress relaxation mechanism underlying the aging of soft solids [31][32][33] . As a matter of fact, viewing the soft glassy material as an elastic continuum and considering that the micro-collapses occur randomly in space (both strong assumptions in most of the materials of interest), the related disruption of the elastic strain field does lead to a compressed exponential relaxation regime, in relatively good agreement with experiments [33][34][35][36] . Understanding (and controlling) the nature of dynamical processes in aging soft solids is crucial to a wide range of technologies, from those involved in material processing to those relying on stable material properties. Moreover, the contrasting experimental findings also raise the very fundamental question of which microscopic mechanism should be prevalent in which material and under which conditions.
Here we use numerical simulations and a model particle gel which, during aging, undergoes structural changes by means of rupture of particle connections. Such ruptures are prototypical of the micro-collapses in the aging of such materials, occurring in regions where higher than average tensile stresses have been frozen-in during the solidification. By quantifying the consequences of the microscopic ruptures, we unravel that the origin of the different relaxation dynamics is in the dramatically different strain transmission if the stress heterogeneities are significantly larger than Brownian stresses. In those conditions, the characteristic timescale for stress relaxation through thermal fluctuations grows well beyond the duration of typical experiments or simulations, but the stress released in microscopic ruptures can still partially relax through elastically induced fluctuations, that are strongly correlated over time and through the structure. The emergent displacement distribution becomes scale-free, featuring a power law tail that points to the long-range elastic strain field in the gel as the main source of structural rearrangements. Thermal fluctuations and Brownian motion disrupt instead the spatial distributions of the local stresses, as well as the persistence of their correlations in time, which leads to the screening of the strain transmission over the structure. Our results provide a unifying framework to rationalize microscopic dynamical processes in a wide variety of soft solids, crucial to their mechanical performance.

Results
Model. We consider a minimal 3D model for a gel composed of N ¼ 62,500 colloidal particles of diameter s and mass m, with their motion described by a Langevin dynamics in the overdamped limit, at low volume fractions (C7%) (see 'Methods' section). The particles spontaneously self-assemble into a disordered interconnected network thanks to a short range attractive well U 2 (combined with a repulsive core) and an additional short range repulsion U 3 , which depends on the angle between bonds departing from the same particle and imparts an angular rigidity to the gel branches, as expected in such open structures [37][38][39][40] . Starting from the same gel configuration, we analyse the dynamics relevant to highly cohesive gels (where significant frozen-in stresses can be developed during solidification), by using different temperatures at fixed interaction strength E, to vary the ratio k B T/E of the Brownian stresses to stress heterogeneities frozen-in during solidification.
In spite of its simplicity, our model captures important physical features of real particle gels and can be used as a prototypical soft amorphous solid. The structural heterogeneities developed during the gel formation entail mechanical inhomogeneities and the simulations reveal the coexistence of stiffer regions (where tensile stresses tend to accumulate) with softer domains, where most of the relaxation occurs 41 . Hence aging is likely to occur via sudden ruptures of particle connections in the regions where the local tensile stresses are higher, as indeed suggested for the ultraslow aging observed in experiments 21 . In the simulations, they are quite rare even at the highest temperatures of interest here (roughly only a few over C10 6 MD steps at k B T/E ¼ 0.05 (ref. 39)) and become hardly observable on a reasonable simulation time window as we consider even smaller k B T/E. To investigate the consequences of the ruptures on a timescale computationally affordable, we periodically scan the local stresses in the gel structure and remove particle bonds where the local stress is the highest, with a fixed rate G (see 'Methods' section). Recombinations of the gel branches are possible but just not observed at the volume fractions and for the time window explored here.
In a rupture event, illustrated in Fig. 1a,b using two subsequent simulation snapshots, the particles depart from each other and move following the part of the network to which they are attached. The topological and elastic constraints exerted by the network, coupled to the local dissipation, result in a progressive slowing down of such motion. The average displacement u ¼ 1=N P N i¼1 jj r i ðt þ dtÞ À r i ðtÞ jj (dt is the integration timestep in the simulations) decreases exponentially between two rupture events, as shown in Fig. 1c for k B T/E ¼ 0. After each event, interparticle forces and/or thermal fluctuations drive the redistribution of the stresses over the structure, with the stress relaxation shown in Fig. 1d   stress tensor (that is, the stress computed from the microscopic particle configurations 42 , see 'Methods' section) plotted as a function of the time. All data shown in the following refer to a removal (or cutting) rate G¼0: is the MD time unit, and correspond to breaking C5% of the initial bonds in total over the simulation time window of 10 5 t 0 considered here. We have carefully investigated how the choice of G affects the results discussed below and determined that they do not significantly change as long as the removal rate is high enough to collect enough statistics (of the ruptures and their consequences) but slow enough to allow for a partial stress relaxation between two rupture events ( Fig. 1d) (see 'Methods' section).
Internal stresses. The analysis of the stress heterogeneities in the network helps disentangle how the stress reorganization and the presence of thermal fluctuations contribute to the relaxation dynamics. Because of the topological constraints produced in the network on its formation, in the initial configuration local stresses have significant spatial correlations. Particles located in the nodes of the gel network, in fact, contribute to the local stresses mostly in terms of tensions, as shown in Fig. 2a, where the gel is coloured according to the value of the local normal stresss n (that is, the normal virial stress computed over a small volume around each particle and containing B10 particles, see 'Methods' section). Analyzing the evolution over time of the largest among the contributions o i n of each particle i to the total normal stress s n , we find that it progressively decreases in all cases, but its fluctuations are clearly different in nature depending on the relative strength of thermal fluctuations (see Fig. 2b). The decrease rate of o i n clearly depends on k B T/E and is much slower than the cutting rate, indicating that the aging of the material corresponds to a slower, large-scale stress redistribution, which emerges over much longer timescales. In the athermal case, this slow evolution of the structure is found to be logarithmic in time. The final distribution of local stressess n (Fig. 2c) still features a few of the pronounced peaks present in the initial one, pointing to the persistence, during the aging, of the spatial correlations of the local stresses. Thermal fluctuations, instead, clearly disrupt such correlations and help redistribute the local stresses much more homogeneously over the gel structure.
Microscopic dynamics. The mechanical heterogeneities and their time evolution during the aging have an impact on the relaxation dynamics. To quantify it, we use the same observables measured in quasi-elastic scattering experiments. From the particle coordinates, we compute the coherent scattering function: where S(q) is the structure factor of the gel and t w the time at which the measurement starts (here t w ¼ 0 for all data sets). From the long time decay of F(q,t) ¼ F t w ¼0 q; t ð Þ ð Þ ; we extract the relaxation time t(q) for different wave vectors (see 'Methods' section) and in Fig. 3a we plot F(q,t) as a function of t/t(q), for q varying between 0.1 and 10 (in units of s À 1 ), at three different temperatures: k B T/E ¼ 10 À 3 (blue), k B T/E ¼ 10 À 4 (green) and for the athermal network (red). The data indicate that the form of the time decay is sensitive not only to q (as a consequence of the structure heterogeneities) 31,43,44 , but also to k B T/E. In all cases, the final part of F(q,t) is well fitted by a decay / exp À t=tðqÞ ð Þ b h i . The exponent b, extracted from the best fit of the past decade of F(q,t), indicate a transition from a stretched to a compressed exponential decay upon decreasing the ratio k B T/E (see Fig. 3b). The dependence of b on k B T/E found here is consistent with experimental observations in refs 25,26. The dependence on q is the same as typically observed in several experiments 22,28,31 , with the b decreasing rapidly with increasing q over lengthscales typically ranging between a few and several particle diameters, up to roughly the order of the mesh size of the gel network (the mean distance between two network nodes along a gel branch is 2p/qC5s at this volume fraction 39 ). For much larger distances, b does not change much with q. The particle mean squared displacement (MSD) is shown in Fig. 3c as a function of the time for the different values of k B T/E. After a localization process at intermediate times, the MSD displays a change in the time dependence at long times, from subdiffusive to super-diffusive, on decreasing k B T/E. The wave vector dependence of the relaxation time t(q) is reported, for different k B T/E, in Fig. 3d. For q corresponding to distances up to C5s, t(q) is strongly dependent on q and on k B T/E and, as long as k B T/E C10 À 3 , it follows the scaling p1/q 2 typically expected when diffusive processes are at play. Reducing k B T/E leads to a different scaling t(q)p1/q, which is consistent with the super-diffusive microscopic motion detected in Fig. 3c. The results obtained for k B T/Eo10 À 3 are consistent with several experiments, finding a compressed exponential dynamics and a 1/q scaling of the relaxation time. Such dependence was ascribed to a superdiffusive, quasi-ballistic particle motion in refs 21,31,45, a hypothesis not directly testable in experiments but proven in Fig. 3c. The ballistic motion reported here and in experiments persists over timescales corresponding to a large number of events. This feature points to the fact that the origin of the ballistic motion is not the single breaking event but the elastic relaxation accumulated in the network over many events.
The agreement between our scenario for small enough k B T/E and the experiments support the idea that frozen-in stresses control the microscopic dynamics also in the experimental samples 23,46 . The relaxation time t(q) from the simulations becomes less sensitive to the wave vector over distances beyond the mesh size of the network, as indeed observed in experiments 20 (see Fig. 3d). Whereas such finding could be interpreted as a loss of dynamical correlations over large distances, we have verified through mechanical tests that the gels are prevalently elastic for all values of k B T/E considered here, supporting the presence of long-range spatial correlations of the dynamics. Moreover, the relaxation time extracted from the incoherent scattering keeps instead increasing on decreasing q and also indicate long-range spatial correlations in the dynamics (see Supplementary Fig. 1). Hence the data in Fig. 3d can be better understood by considering that, although local rearrangements are possible, over larger distances the material is solid and elastically connected in spite of being structurally heterogeneous and sparse, and the time correlations of the density fluctuations over those lengthscales are therefore locked-in.
Beyond the capability of experiments, the simulations give us direct access to the distribution of the displacements over a lag time in between two rupture events, which, for k B T/E410 À 3 , is a bell-shaped curve with exponential tails (Fig. 4a), similar to those typically seen in supercooled liquids 38,47,48 . However, on decreasing k B T/E-0, the probability distribution function becomes instead progressively scale-free. Such findings support our interpretation of Fig. 3d and the idea that the dynamics are spatially correlated and cooperative over extended domains for all k B T/E. We visualize the displacement field in two representative snapshots of the gel network (Fig. 4b,c) where the colour code indicates the amplitude of the displacements measured along the network after a rupture event, for k B T/ E ¼ 10 À 3 and k B T/E ¼ 0, respectively. The visualizations highlight how the long-range elastic response of the network dominates the displacement pattern in the athermal limit, where it extends indeed over several network meshes, and it is instead screened by the small displacements induced by thermal fluctuations for larger k B T/E (see Supplementary Movies).
For a homogeneous elastic material, the amplitude of the strain u induced at a distance r from the rupture event decreases as u(r)p1/r 2 (ref. 50). Given the probability P(r) of having a rupture event at a distance r and if, in a mean field (MF) approximation, rupture events occur anywhere with the same probability, the probability distribution function for a displacement (or strain) of amplitude u can be estimated as PðuÞ $ PðrÞ j dr=du j $ r 2 r 3 $ u À 5=2 . In the athermal limit   (Fig. 4a), which is in fact compatible with a compressed exponential decay of the scattering functions with b ¼ 3/2 (refs 33,35), as shown in Fig. 3b. The agreement with the MF prediction may suggest that structural disorder and correlations between rupture events are in the end negligible. Nevertheless, here (as well as in many experiments) the structure is heterogeneous, the correlations of events are crucial and the MF assumptions are unlikely to truly hold 35 . We propose therefore a different explanation, that is, structural disorder and correlations between rupture events, in spite of being present, do not disrupt the elastic connectivity of the material because, when k B T/EB0, the relaxation of the stresses through elastic fluctuations does not weaken their spatio-temporal correlations. As a matter of fact, in addition to the spatial correlations discussed above, the total stress also features large fluctuations in time, that depend on k B T/E and become strongly correlated as k B T/E-0 (Fig. 5).
Stress fluctuations and stress relaxation. The fluctuations of one component of the total stress ds¼s xx À hs xx i ð Þ , measured over a time window significantly larger than the time interval between two rupture events, are plotted in Fig. 5a for k B T/E ¼ 10 À 3 and k B T/E ¼ 0. When Brownian stresses become negligible with respect to the stress heterogeneities, such fluctuations clearly deviates from randomly fluctuating values, even after many ruptures have occurred. Their probability distribution function is Gaussian for larger k B T/E, whereas it develops non-Gaussian tails as k B T/E-0 (Fig. 5b). The non-Gaussian processes and the intermittent fluctuations measured over the same time window correspond indeed to persistent time correlations (see Supplementary Fig.2).
When monitoring the stress fluctuations over the whole duration of the simulations, we find evidence of the progressive aging of the structure that eventually emerges from the correlated stress redistribution (see Fig. 5c, inset for k B T/E ¼ 0 data). Accordingly, the time correlations of the fluctuations of the total stress over these longer timescales develop a pronounced plateau on decreasing k B T/E and eventually extend to the whole duration of the simulation in the athermal limit (Fig. 5c, main frame).
Time-and space-resolved scattering experiments have indeed found evidence of highly intermittent dynamical processes 43 , hence Fig. 5a,b and the intermittency in the stress fluctuations detected here provide a first coherent physical picture also for these observations, beyond the aspects that can be captured by MF approaches. Our analysis reveals that indeed the spatio-temporal correlation patterns detected in the microscopic displacements and in the density fluctuations (Figs 3a-c and 4a) stem from the way thermal fluctuations, on increasing k B T/E, help redistribute elastic stresses (and strains) in the material. When thermal fluctuations are weak, the microscopic dynamics emerge instead from elastically driven stress fluctuations that remain strongly correlated in space and time.

Discussion
Our analysis provides a vivid microscopic picture for the hypothesis underlying recent theories of aging in soft solids 34,35 and has implications for a potentially wider range of materials.
The control parameter k B T/E simply reflects the ratio between the timescales governing stress relaxation, respectively, through thermal fluctuations (Z s 3 /k B T) and elastic recovery (Z s 3 /E) in the material (Z being the viscous damping). Hence it helps identify the conditions for which the elastically driven intermittent dynamics emerge in different jammed solids. First, we note that while we have considered here only one type of micro-collapses, the distruption of the elastic strain field due to the rupture of a branch of the gel is basically the same, in a first approximation, as the one induced by a recombination of the gel branches (or by other types of possible microscopic events that help relax frozen-in stresses) 33 . Therefore, we expect our general picture to have a much wider relevance. Second, the physical mechanisms we propose help rationalize several experimental observations in very different materials, ranging from biologically relevant soft solids to metallic glasses 31,45,46,50 . In a quasi-equilibrium scenario, enthalpic and thermal degrees of freedom may still couple and stress correlations decay relatively fast. When the material is deeply quenched and jammed, instead, recovering the coupling between the distinct degrees of freedom and restoring equilibrium will require timescales well beyond the ones accessible in typical experiments or simulations. The result will be intermittent dynamics and compressed exponential relaxations. The competition between Brownian motion and elastic effects through the relaxation of internal stresses illustrated in this work suggests different scenarios for the energy landscape underlying the aging of soft jammed materials. When thermal fluctuations screen the long-range elastic strain transmission, microscopic rearrangements may open paths to deeper and deeper local minima in a rugged energy landscape. Compressed exponential dynamics, instead, evoke the presence of flat regions and huge barriers, with the possibility of intermittent dynamics, abrupt rearrangements and avalanches 34,35 . Investigating how such different dynamical processes couple with imposed deformations will provide a new rationale, and have important implications, for designing mechanics, rheology and material performances.

Methods
Numerical model and viscoelastic parameters. The particles in the model gel interact through a potential composed of two terms. The first force contribution derives from a Lennard-Jones-like potential of the form: The second contribution confers an angular rigidity to the inter particles bonds and takes the form: The strength of the interaction is controlled by L(r) and vanishes over two particles diameters: where Y denotes the Heaviside function. The evolution of the gel over time is obtained by solving the following Langevin equation for each particle: U r i ; ::: where s is the particle diameter, x(t) is a random white noise that models the thermal fluctuations and is related to the friction coefficient Z f by means of its variance hx i ðtÞx j ðt 0 Þi ¼ 2Z f k B Td ij dðt À t 0 Þ. To be in the overdamped limit of the dynamics Z f is set to 10, and the timestep dt used for the numerical integration is dt ¼ 0.005. The parameters of the potential are chosen such that the disordered thin percolating network starts to self assemble at k B T/E ¼ 0.05. One convenient choice to achieve this configuration is given by this set of parameters: A ¼ 6.27, a ¼ 0.85, B ¼ 67.27, y ¼ 65°and w ¼ 0.3. The system is composed of N ¼ 62,500 particles in a cubic simulation box of a size L ¼ 76.43s with periodic boundary conditions, the number density N/L 3 is fixed at 0.14, which corresponds approximatively to a volume fraction of 7%. All initial gel configurations are the same, prepared with the protocol described in ref. 40, which consists in starting from a gas configuration (k B T/E ¼ 0.5) and letting the gel self-assemble upon slow cooling down to k B T/ E ¼ 0.05. We then quench this gel configuration by running a simulation with the dissipative dynamics m d 2 ri dt 2 ¼ À r ri U À Z f dri dt until the kinetic energy drops to zero(10 À 24 ). All simulations have been performed using a version of LAMMPS suitably modified by us 52 .
Stress calculation and cutting strategy. We let the initial gel configuration evolve with equation (4) for each of the different values of k B T/E considered here, while using the following procedure to cut network connections. At each timestep, we characterize the state of stress of a gel configuration by computing the virial stresses as s ab ¼ À 1 V P i w i ab , where the Greek subscripts stand for the Cartesian components x, y, z and w i ab represents the contribution to the stress tensor of all the interactions involving the particle i and V is the total volume of the simulation box 52 . w i ab contains, for each particle, the contributions of the two-body and the three-body forces evenly distributed among the particles that participate in them: The first term on the RHS denotes the contribution of the two-body interaction, where the sum runs over all the N 2 pair of interactions that involve the particle i. (r i , F i ) and (r 0 , F 0 ) denote, respectively, the position and the forces on the two interacting particles. The second term indicates the three-body interactions involving the particle i. We consider a coarse-graining volume O cg centered around the point of interest r and containing around 9-10 particles on average, and define the local coarse-grained stress based on the per-particle virial contribution ass ab r ð Þ ¼ À P i2Ocg w i ab =O cg . For a typical starting configuration of the gel, the local normal stresss n ¼ ðs xx þs yy þs zz Þ=3 reflect the heterogeneity of the structure and tend to be higher around the nodes, due to the topological frustration of the network. We consider that breaking of network connections underlying the aging of the gel is more prone to happen in the regions where local stresses tend to be higher, as found also in refs 40,41. Hence, to mimic the aging in the molecular dynamics simulations we scan the whole structure of the gel and remove one of the bonds (by turning off the well in U 2 ) whose contribution to the local normal stress w i n ¼ w i xx þ w i yy þ w i zz =3 is the largest (prevalently bonds between particles belonging to the network nodes). As the simulation proceeds, local internal stresses redistribute in the aging structure of the gel and the locations of more probable connection rupture (as well as their number) change over time. All simulations discussed here have been performed with a rate G ¼ 0:04t À 1 0 , corresponding to removing only B5% of the total network connections over the whole simulation time window. We have run several tests varying the removal rate in order to verify that changing the rate (having kept all other parameters constant) does not modify our outcomes and the emerging physical picture. Overall, varying G over nearly two orders of magnitudes, we recover the same results, as long as the t r ¼ 1/G between two rupture events allows for at least partial stress relaxation (see Fig. 1).
Stress autocorrelation function. The autocorrelation function of the stress fluctuations ds ¼ s xx À hs xx i is computed as follows: C ds ¼ f ds t ð Þ=f ds 0 ð Þ, where f ds is the auto-covariance and takes the form f ds t ð Þ¼ P n À t t¼0 s xx t þ t ð ÞÀhs xx i ð Þ s xx t ð Þ ð À hs xx iÞ and hs xx i¼ 1 n þ 1 P n t¼0 s xx t ð Þ. We have computed the autocorrelation function of the stress fluctuations over partial time series and over the whole duration of the simulations (see Fig. 5 and Supplementary Fig. 2).
Fitting procedure for the intermediate scattering functions. To extract the b exponent and the structural relaxation time t q , we fit the last decay of the coherent scattering function F(q,t) t w ¼ 0 simultaneously with the incoherent scattering function F s ðq; tÞ F s ðq; tÞ tw¼0 , defined as F s ðq; tÞ ¼ P N k¼1 exp½ À iq Á r k ðtÞ À r k ð0Þ ð Þ =N. This observable is less noisy and has the same exponent as the coherent scattering function shown in Fig. 3 (see also Supplementary Fig. 3). Having fitted F s (q,t) with a stretched exponential type of the form A þ B expð À t=t s Þ b and having extracted the exponent b, such exponent is used to initiate the fit of F(q,t) by locking this parameter and letting all the others free. This procedure helps us to improve the quality of the fit used to extract the relaxation time t q from the coherent scattering.
Data availability. The authors declare that all data supporting the findings of this study are available within the article and its Supplementary Information Files. All raw data can be accessed on request to the authors.