Mode-coupling analysis of residual stresses in colloidal glasses

We present results from computer simulation and mode-coupling theory of the glass transition for the nonequilibrium relaxation of stresses in a colloidal glass former after the cessation of shear flow. In the ideal glass, persistent residual stresses are found that depend on the flow history. The partial decay of stresses from the steady state to this residual stress is governed by the previous shear rate. We rationalize this observation in a schematic model of mode-coupling theory. The results from Brownian-dynamics simulations of a glassy two-dimensional hard-disk system are in qualitative agreement with the predictions of the theory.


Introduction
Amorphous so solids that are produced by owing them into shape (entailing a quench into a nonequilibrium glassy state) display residual stresses. 1,2 The internal stresses that build up during ow do not relax fully, so that some part of them persists in the solid that is formed by kinetic arrest in the uid; indefinitely, in the ideal-glass case. This was recently demonstrated qualitatively for a simple setup where glass-forming colloidal suspensions of nearly-hard-sphere particles were observed aer cessation of steady shear from their nonequilibrium stationary state (NESS), combining macroscopic and microscopic experiment, computer simulation, and mode-coupling theory of the glass transition (MCT). 1 In fact, the appearance of persistent residual stresses is known empirically since centuries, and is neither restricted to so solids, nor to amorphous ones. An early demonstration involves small drops of molten (ordinary window) glass that fall into cold water and thus solidify very rapidly. These drops, known as Prince Rupert's drops or Dutch tears since the 17th century, 3 have a surprisingly shock-resistant body (capable of withstanding the blow of a hammer) but explode dramatically upon the slightest damage done to their tail. The high mechanical stability arises from a shield of compressive residual stresses along the drop's surface that compensates high tensile stresses in the less rapidly cooled inner material and acts to stop crack formation. [3][4][5] Clipping the tail releases this frozen-in network of residual stresses, rendering the material less strong, and in this case even unstable. 6 Under normal conditions, the residual-stress network on the other hand is stable over decades in time. Safety glass and modern smartphone cover glasses (where chemical processes are used to impose residual stresses) exploit this effect to ne-tune the desired mechanical properties of the product. 7-9 Aside from the amorphous state, residual stresses are decisive for the fatigue crack growth and hence the long-term stability of railway rails under the strong external forces caused by trains. 10,11 One encounters residual-stress related material phenomena even in biophysics. For example the cytoskeleton of cells is prestressed (mainly due to the action of myosin motors). 12 The unsurpassed material properties of spider silk are attributed to residual stresses. 13 The control of residual stresses during the production stage is decisive in tuning the long-term stability of certain thin polymer lms. 14 Shear cessation experiments 15,16 arguably provide the cleanest setup to investigate residual stresses. Close to kinetic arrest, cessation experiments have addressed the relaxation of stresses 17 and the evolution of linear viscoelastic moduli aer the cessation of ow 18,19 in laponite gels, and aging effects in depletion gels. 20 We consider even simpler model systems with hard-sphere like interactions, at xed temperature and density, initially subject to homogeneous, stationary simple-shear ow. Residual stresses arise when the system does not relax back to equilibrium aer cessation of the ow, but gets trapped in a nonergodic state on the way. This is the case for (ideal) glasses, the quiescent state of which is nonergodic, but which uidize under shear, allowing the start of the cessation experiment from a well-dened, unique NESS.
We seek to improve the rst-principles theoretical understanding of residual stresses: they are denitely a nonequilibrium and nonlinear-response effect, and are found to depend strongly on the history of past material deformation. For colloidal suspensions as model so solids, the combination of an integration-through transients (ITT) approach to the nonlinear response with the mode-coupling theory of the glass transition (MCT) [21][22][23][24] allows us to address such issues starting from a particle-based microscopic description. The theory predicts shear-rate dependent residual stresses aer the cessation of steady shear, 1 in qualitative agreement with experiment and simulation.
It should be stressed that the theoretical modeling of residual stresses is far from trivial. A successful model used to describe the nonlinear time-dependent rheology of so materials is the so glassy rheology (SGR) model. 25 Although it explicitly addresses the aging dynamics of kinetically arrested solids, it implies the relaxation of stresses back to zero aer cessation of ow. 26 In this article, we discuss residual stresses in model hardsphere like glasses as explained by MCT. The theory predicts the persistence of a non-relaxing component of the stress aer ow cessation and captures its strong dependence on past ow history and on the mechanism of yielding. We compare with Brownian dynamics computer simulations on a two-dimensional hard-sphere glass to test the predictions of the theory, nding qualitative agreement.
This paper is organized as follows: we describe in Section 2 the schematic-MCT model and give details of the event-driven Brownian-dynamics (ED-BD) simulations we performed. In Section 3, the main results for the partial relaxation of stresses, and their analysis using the schematic-MCT model is given. Section 4 presents conclusions and outlook.

Methods
We discuss the evolution of the shear stress s(t) aer the cessation of steady simple-shear ow with shear rate _ g, instantaneously switched off at time t ¼ 0. The time is measured in units of the free-particle Brownian relaxation time s 0 ¼ d 2 /D 0 , where d is the diameter of the particles (the unit of length) and D 0 is their free diffusion coefficient. Stresses are reported in these natural units, k B T/d D , where k B T is the energy associated with thermal uctuations and D is the dimensionality (D ¼ 2 for our BD simulations, and D ¼ 3 for the theory, although this difference in dimensionality is not borne out in the schematic-MCT model).
Flow cessation is implemented in the theory by assuming that instantaneously, for all t > 0, the (assumed) homogeneous velocity-gradient eld vanishes everywhere. In our Browniandynamics simulations, this is reected, but in comparison with molecular-dynamics (MD) computer simulations in general, and with experiment, one has to keep in mind that this is a simpli-cation. If shear ow is imposed by moving boundaries, our discussion refers to the case where these boundaries are suddenly xed at their current position, and held there so that no further strain relaxation is allowed. This is to be distinguished from zero-stress boundary conditions, where stress relaxation by adjusting the total strain would occur in addition. 27,28 We also neglect the effect of transiently inhomogeneous ow elds as they might arise shortly aer stopping the ow. 29 The transients should be associated with a comparatively short timescale of transverse-momentum diffusion, as conrmed by thermostatted MD simulations. 30,31

Theory
The integration-through transient 32,33 approach to colloidal rheology assumes that the system is at rest and in Boltzmann equilibrium for the innitely distant past, t / ÀN. The ow is subsequently switched on and causes a nonequilibrium perturbation of the Smoluchwoski operator describing the temporal evolution of the N-particle distribution function. A formal manipulation allows us to derive the nonequilibrium, nonlinear generalization of the Green-Kubo relationship for the stress tensor. 23,24 MCT inspires the subsequent approximation in terms of density uctuations to wave vectork. For the shear stress in simple shear ow along the x-direction and with gradient in the y-direction, one gets Here,k(t, t 0 ) ¼k + k xẽy g tt 0 is the shear-advected wave vector, describing the affine deformation of a plane-wave uctuation through the accumulated strain g tt 0 ¼ _ gðsÞds between two points in time, t 0 # t. For the case of shear cessation, we set _ g(t) ¼ _ gQ(Àt) with the Heaviside step function Q. Eqn (1) consists of a history integral over transient density-correlation functions fk (t,t 0 ) (t, t 0 ) involving strain-dependent density-stress coupling vertices that are given in terms of the quiescent static structure factor S k . The transient correlation function is dened as the equilibrium average of uctuations evolved with the nonequilibrium dynamics. An important asset is that transient correlation functions only depend on the ow history between their two time arguments t 0 and t, and are oblivious of any changes in the ow before t 0 .
Eqn (1), together with the ITT-MCT equations for fk(t, t 0 ), allows us in principle to calculate the time-dependent shear stress in response to an arbitrary time-dependent shear ow _ g(t), given the static structural factors of the system. To reduce the computational effort, one oen employs an additional isotropic approximation; in ref. 1 such an isotropically sheared hard-sphere model (ISHSM) was solved for the case of shear cessation.
In schematic-MCT models, the full equations of motion of the theory are simplied ad hoc, dropping the dependence on wave vectors. This emphasizes the temporal correlations as the main origin for nonlinear rheology close to the mode-coupling glass transition. One possible schematic simplication of eqn (1), in the following referred to as model A, reads 34 with the schematic version of the generalized dynamic shear modulus Here, v s is a coupling constant that adjusts the energy scale of the schematic model. For a quantitative comparison with experimental data on hard-sphere suspensions in three dimensions, v s ¼ 100k B T/d 3 has been suggested. 34 Since we are mainly interested in a qualitative discussion here, we stick to this value even though our Brownian dynamics simulations are performed in 2D. The dynamic shear modulus G(t, t 0 ) entering eqn (1) is in the schematic model replaced by the square of the density correlation function. As a consequence, the integrand in eqn (2) is manifestly positive under steady shear. Yet, in many dense glass-forming systems, so-called stress overshoots are observed aer the startup of steady shear for accumulated strains g ¼ _ gt of a few percentage (assuming shear ow to be started at t 0 ¼ 0). The resulting stress-strain curves s(g) display an intermediate maximum before decreasing towards the steady-state value, at about g z 0.1 for hard-sphere systems. The position and strength of this overshoot depend on the details of the interaction and the sample preparation and age. 30,35,36 The overshoot corresponds to stress-over relaxation visible as a negative dip in G(t, 0) during the structural-relaxation process. It describes the capability of the viscoelastic system to transiently store more elastic energy than is maintained during plastic ow, and an elastic recoil process during yielding. Comparing the ITT Green-Kubo relationship, eqn (1), with its schematic simplication, eqn (2), one recognizes that reduction of the coupling coefficients driven by shear advection is missing from the schematic model. In the microscopic MCT model, a certain amount of negative coupling of density uctuations into the stress can occur. One can interpret this reversible suppression as the effect of anelasticity on the average structure, rather than that of plasticity, which is expressed through the strain-induced irreversible decay of the correlation function.
For the purpose of quantitative ts to such stress-strain curves, a schematic model (here referred to as model B) has been proposed, including an empirical term that captures the strain dependence of density-uctuation-stress couplings by setting, 37 We denote here the obvious generalization of ref. 37 to nonsteady shear. The parameters g * and g ** model the position and width of the overshoot, and are adjusted to t the experimental data. They show a weak dependence on the shear rate themselves in ts of the steady-shear rheology of a hard-sphere suspension. 37 To keep the discussion simple, we neglect this and set g * ¼ 0.105g c and g ** ¼ 0.14g c . This produces a pronounced stress overshoot in startup ow, allowing us to discuss the inuence of the elastic-recoil effect on the stress relaxation aer cessation in a qualitative way. The scale parameter v 0 s is adjusted to match the dynamic yield stress with that of model A with a xed v s .
The equation of motion for the transient density correlation functions in the schematic MCT (for both models) reads 34 v t fðt; t 0 Þ þ fðt; t 0 Þ þ where we have set the initial relaxation time to unity dening the unit of time. The memory kernel m(t, t 00 , t 0 ) describes structural relaxation and its modication through the applied ow. The mode-coupling theory of the glass transition closes the equation of motion by approximating the memory kernel as a quadratic functional of the density correlation functions themselves, assuming density uctuations to be the dominant slow dynamical variable. Mimicking the microscopic-MCT expression, we employ the common F 12 model for time-dependent shear, 34 Here, h(t, t 0 ) is a function describing the suppression due to shear of the memory effects that cause the slow structural relaxation. v 1 and v 2 are the coupling coefficients of the model. They describe the thermodynamic state of the system. For small coupling coefficients, the quiescent solution of the F 12 model is liquid like, where density uctuations fully relax: ) identied as the ideal glass transition of the model, where a nonergodic solution f(t, t 0 ) / f > 0 rst appears. 38 Choosing the point determined by v 2 c ¼ 2 ensures that the asymptotic dynamics of the F 12 model matches with that expected for hard-sphere-like suspensions. For simplicity, we set v 2 ¼ 2 throughout. We introduce a separation parameter 3 to indicate the distance to the glass transition: Accumulated strain decorrelates the MCT memory kernel; this is described by the empirical strain-reduction function h(t, t 0 ) ¼ 1/[1 + (g tt 0 /g c ) 2 ]. Here, g c ¼ 0.1 is a parameter that sets the scale for the elastic limit of the solid.
The ITT-MCT equations of motion are solved numerically. We need to keep the full dependence of correlation functions on their two time arguments. However, for the transient correlation functions the calculation is eased by noting that for cessation of steady ow at t ¼ 0, the time half-plane (t, t 0 < t) can be split into three regions. In particular, for the two regions t > t 0 > 0 and t 0 < t < 0, the equations for the transient correlation functions restore time-translational invariance since the shear rate is constant between any pair (t 0 , t) there. The nontrivial predictions of the model all stem from the behavior of the two-time correlation functions spanning the cessation point, t 0 < 0 < t. A numerical scheme to deal with these correlation functions in the case of step changes in shear rate has been devised. 39

Simulation
We perform event-driven Brownian dynamics simulations 40 for a 2D binary mixture of hard disks. The system consists of N A ¼ 500 particles with diameter d ¼ 1, and N B ¼ 500 particles with diameter d B ¼ 1.4, and has been widely studied as a model glass former, 41,42 including its nonlinear rheology in the steady state, large-amplitude oscillatory shear, and the transient dynamics aer shear startup. 37,[43][44][45][46] The only control parameter is the packing fraction 4 giving the ratio of the volume occupied by all the disks to the volume of the system. The mode-coupling glass transition point of this binary mixture has been determined from extensive simulations, 41 4 c z 0.795. For comparison, jamming in this system occurs around 4 rcp z 0.84, 47 signicantly above the density regime we are interested in.
In this contribution, we focus on the shear-cessation dynamics of a glassy state, and set 4 ¼ 0.81 throughout. Brownian dynamics is approximated by a sequence of eventdriven simulation steps of length s B , ensuring that no particle overlaps occur, aer which the particles are assigned to new Gaussian random displacement vectors. 40 For times t [ s B , the dynamics reects that expected from the N-particle diffusion equation for hard spheres with a free-diffusion coefficient The shear ow is imposed on the simulation by Lees-Edwards periodic boundary conditions, 48 and by imposing a corresponding deterministic bias for the random displacements. The resulting stationary state is homogeneous for all the shear rates considered below.
Stresses are measured in ED-BD by associating with each "collision" of two Brownian particles i and j a momentum transfer Dṽ ij , where each Brownian displacement Dx i is assigned to a velocityṽ i ¼ Dx i /s B . The shear stress is then given by where Dr ij (t c ) is the separation vector between two particles during their collision. The sum runs over all collision times t c in a small time interval Dt used to smoothen the result. We choose Dt in the range [25s B , 150s B ] h [0.00125s 0 , 0.75s 0 ] depending on the shear rate, and the time aer startup and cessation (using smaller time intervals close to the points where the shear rate changes). To improve statistics at the largest shear rates ( _ gs 0 $ 4 Â 10 À1 ), averaging intervals up to Dt ¼ 10s 0 have been used.
Aer startup of shear from a well-aged glassy conguration, the simulations were run up to an accumulated strain of g ¼ 1 to ensure that a steady state had been reached. The nal steadystate conguration is used to set the origin of time thereaer. Aer cessation, runs were performed with up to 5 Â 10 6 BD simulation steps, corresponding to a time t/s 0 ¼ 2.5 Â 10 4 . To obtain sufficient statistics for the stresses, a large number N r of independent simulation runs is required. We have used N r ¼ 393 (524, 801, 615, 392) runs for the calculation of s(t) aer cessation from shear rates _ gs 0 ¼ 4 Â 10 Àn with n ¼ 0 (1, 2, 3, 4), supplemented with 359 (200, 187, 193) shorter runs (with 1 Â 10 6 BD steps) for n ¼ 1 (2,3,4).
The results for s(t) from the ED-BD simulations are shown in Fig. 1. The region _ gt < 0 corresponds to simulations of startup shear, as discussed previously for this system 37 and for MD simulations of similar model glass formers. 30 This initial part of the s-versus-_ gt curves reects the stress-strain curves s(g) discussed in these references. They exhibit pronounced stress overshoot phenomena, a fact that will be related to the stress relaxation aer cessation below. At the highest shear rate, _ gs 0 ¼ 4, one notices a small undershoot aer the maximum in the stress-strain curve, and some smaller oscillations in s(t) until the point where the shear was switched off. Therefore, some care has to be taken to interpret the results for this shear rate, and we will only discuss generic features that are qualitatively identical also for the lower shear rates. While the present ED-BD simulation algorithm remains a well-dened, overlap-free, rejection-free Monte Carlo scheme, 40 its approximation of the Smoluchowski equation for the Brownian-particle motion becomes questionable for shear rates much larger than _ gs 0 z 4.

Results and discussion
We rst present results for the stress relaxation aer simpleshear cessation in our two-dimensional Brownian dynamics simulation. Fig. 2 shows the shear stress s(t) as a function of time t since cessation, for various shear rates. With increasing shear rate, the steady-state shear stress s ss ( _ g) increases; this is seen as the initial value of the s(t) curves in Fig. 2. Some of this shear stress relaxes aer cessation of ow, on a time scale set by the past shear rate. This is evidenced by the lower panel in the gure, where scaling of the normalized s(t)/s ss -versus-_ gt curves onto an apparent master curve is observed for an initial time window, _ gt ( 0.1. Curves for lower pre-shear rates deviate from the master curve at previous rescaled times. For the highest shear rates, the decay at intermediate times is approximately described by a 1/ _ gt form (dash-dotted line in the gure). For long times, a nonzero plateau value is indicated in the simulations. Note that large uctuations hamper the exploration of the long-time regime in the simulation; this is in particular true for the larger shear rates shown in Fig. 2. Here, the initial part of the relaxation covers around two orders of magnitude in stress, and hence reaches the noise limit of the statistical sampling. Nevertheless, the simulation data are compatible with the emergence of a nite long-time limit in the ideal glass, , a persistent residual stress. To demonstrate this, we have shown in Fig. 2 as dashed lines the time-averaged stress s over the interval t/s 0˛[ 9800, 25 000] for the simulations with shear rates _ gs 0 # 4 Â 10 À2 . For these cases, s is signicantly different from zero. Dotted lines in the upper panel of Fig. 2 show exemplary results from schematic MCT, discussed in more detail below, to highlight the expected behavior in the ideal glass. The theory result overestimates the amount of residual stress that persists in the glass, but describes the partial relaxation from the steady state. Our result corroborates previous MD simulation results for a binary Yukawa mixture with a dissipative-particle-dynamics thermostat and temperature as the control variable. 1 In MCT, a nite residual stress emerges since s(t) is determined from the past ow history, and innitely long-lasting memory effects in the ideal glass cause a nonvanishing contribution to the generalized dynamical modulus G(t, t 0 ) from times t 0 < 0 (i.e., before cessation) even as t / N. We will investigate the shape of these correlation functions in more detail below. Fig. 3 demonstrates the stress relaxation aer the cessation of ow in the schematic-MCT models introduced above. A state in the ideal glass (distance parameter 3 ¼ 0.01; solid lines) was chosen. For comparison, we also show stress relaxation curves for a state point in the liquid (3 ¼ À0.01; dashed lines). There, stresses ultimately relax to zero as the system reaches an equilibrium state again. This nal relaxation is determined by the structural-relaxation time scale s a of the quiescent equilibrium correlation function, which is independent of the pre-shear rate _ g. Since the simulation data shown in Fig. 2 show relaxation essentially only on a time scale connected with _ g, this can be taken as an indicator that they are not connected with equilibrium-liquid like relaxation.
Comparing with the simulation results, Fig. 2, one recognizes two main differences from schematic model A: in the simulation, the initial transient stress relaxation is stronger, so that the nal residual-stress plateau is hardly discernible in a semi-log representation (upper panel of Fig. 2); only in a doublelogarithmic representation (lower panel), the nite residual stress becomes apparent. The second difference concerns the dependence of the residual stress on the pre-shear rate: in the simulations, we observe s N ( _ g) to be a decreasing function of increasing _ g. In schematic model A with xed vertex v s , the residual stress s N instead increases with increasing pre-shear  rate; this is also the case in the isotropic approximation to the microscopic ITT-MCT. 1 The trend exhibited by the simulation results typically leads to a crossing of s(t)-versus-t curves belonging to different _ g. Stronger shear causes higher steady-state stresses (so that the initial value of the s(t) curve increases), but a larger amount of these stresses can relax aer cessation from strong pre-shear, than from weaker pre-shear. Intuitively, stronger shear uidizes the glassy structure more severely, allowing for more effective particle rearrangements in order to relax stresses aer the driving is switched off. Note however that the decrease of s N with increasing pre-shear rate is a general, but not a universal effect; it is observed in the MD simulations mentioned above, in experiments on hard-sphere-like suspensions of core-shell PS-PNIPAM particles, and also for hard-sphere suspensions of PMMA particles close to the glass transition. 1 However, these PMMA suspensions at higher packing fraction displayed an increase of the residual stress with increasing pre-shear rate. Coincidentally, the systems for which a decreasing s N ( _ g) was observed are also those with strong stress overshoot phenomena under startup ow.
The different trends of s N ( _ g) can indeed be connected to the amount of elastic recoil that causes the stress overshoots under startup ow. As explained above, schematic model A does not contain this mechanism. It only describes the relaxation of stresses due to irreversible, plastic rearrangements, in the theory identied as those that cause strain-induced decorrelation of advected density uctuations. The introduction of a time-dependent coupling between stresses and those density modes, introduced into the schematic model by a time-dependent vertex v s (t, t 0 ), accounts for a further relaxation caused by anelastic deformation due to the affine effects of shear. Including this mechanism of stress overshoots indeed has the effect of reducing the residual stress s N ( _ g) such that it becomes a decreasing function of increasing _ g. The lower panel of Fig. 3 demonstrates this for model B with the time-dependent vertex with parameters g * and g ** from eqn (4). With this choice, schematic model B reproduces the crossing of s-versus-t curves observed in the simulation. Curves from model B corresponding to two shear rates used in the simulation, _ gs 0 ¼ 4 Â 10 À4 and 4 Â 10 À3 , albeit at a smaller distance parameter 3 ¼ 0.0001 (corresponding to the simulated state which is rather close to the glass transition), have also been added to the simulation results for comparison (the dotted line in Fig. 2). Note that we did not attempt to t the simulation data by adjusting the model parameters.
The qualitative similarity between the MCT model and the simulation results becomes clearer in a double-logarithmic representation of s(t) as a function of t or rescaled time _ gt. This is shown for the schematic-MCT models in Fig. 4. As in the simulation results, one notes reasonably good data collapse for the normalized stress s(t)/s ss for intermediate rescaled times _ gt: all liquid-and glass-state curves for the shear rates shown almost coincide with their initial relaxation from the steadystate value. It should however be noted that in the schematic-MCT models, this scaling does not hold strictly. This can in particular be noted for model B with the strain-dependent vertex, shown in the lower panel of Fig. 4. The relative amount of stress relaxation aer cessation, Ds N ¼ (s ss À s N )/s ss , is described qualitatively and correctly by both variants of the schematic model: Ds N increases as a function of increasing preshear rate, in agreement with all known experimental and simulation data.
While the simulation results shown in Fig. 2 show indications of a decay $1/ _ gt, the schematic-MCT curves are signicantly less steep. In model A with a xed v s , the initial decay approximately follows s $ t Àa with a z 0.33 as expected from the asymptotic power law for relaxation onto the nonergodic plateau derived in MCT. 38 The s(t)/s ss -versus-_ gt curves from the experimental data of ref. 1 similarly show a power-law-like relaxation regime, s(t) $ ( _ gt) Àx . For the PS-PNIPAM system, one estimates x z 0.5, while the PMMA suspension is described by an exponent x T 0.67 close to the glass transition (4 ¼ 0.587), and x z 0.4 for the state that was investigated deeper in the glass (4 ¼ 0.614).
It is worth noting that the so glassy rheology model predicts such power law decay, with an exponent x that is identical to the SGR's effective temperature parameter, 26 and independent of the shear rate. In the SGR model, one thus expects the stress relaxation to be slower and deeper in the glass, with x z 1 close to the transition. In this respect, the initial decay of s(t) observed in our simulations and in the PMMA suspensions follow this expectation, while the value x z 0.5 observed for the PS-PNIPAM system is lower than expected (given that these experimental data correspond to a state not too far from the glass transition). The _ g-dependent slowing down of the stress relaxation (visible, e.g., in the lower panel of Fig. 2) interpreted as the approach to a nite residual stress is a feature that differs signicantly from SGR. Fig. 5 shows the dependence of the residual stress s N ¼ s( _ g;t ¼ N) as a function of the previous shear rate _ g. Symbols show experimental and simulation data. Values for s N for three different hard-sphere-like colloidal suspensions were taken from ref. 1, where they had been determined as the stress reached at a xed time towards the end of the measurement. We add our results from ED-BD computer simulations of the 2D hard-disk mixture (square symbols in the gure). Here, to determine s N , the simulated s(t) curves have been averaged over a certain time window, as explained in connection with Fig. 2 (dashed lines here). In line with the discussion above, the residual stress values s N ( _ g) decrease with increasing pre-shear rate for most experimental and simulation data, with the exception of the highest-density PMMA suspension where an increase is observed. Note again that this latter result is a direct consequence of the experimental protocol, i.e., of keeping the overall strain xed aer ow cessation. Otherwise, the system should start to ow since the frozen-in residual stress s N exceeds the dynamic yield stressand hence likely also the static yield stress that determines the maximum stress the solid can sustain. Even for the cases where s N < s ss , slow creep should be the result of free-strain boundary conditions, leading to additional relaxation phenomena. 36,49 Lines in Fig. 5 show the results for the schematic-MCT model for several state points in the ideal glass as labeled by their distance parameter 3. As observed above, model A with a constant v s predicts an increase of the residual stress with increasing pre-shear rate (dashed lines). Including the straindependent correction to the vertex according to eqn (4), model B predicts the observed decrease of s N as a function of increasing _ g (solid lines). Note that with our choice of model parameters g * and g ** , eqn (4) produces a marked stress overshoot in startup ow. As the density increases beyond the glass transition, the magnitude of this stress overshoot in hard-sphere colloidal suspensions diminishes; 36 a feature that is likely connected to the approach to random-close packing in these systems. Hence, g * and g ** should acquire a density dependence that we have not accounted for. Let us further note that with our choice of xed g * and g ** , the integral in eqn (2) together with eqn (4) is no longer guaranteed to be positive at large shear rates. The corresponding s N curves in Fig. 5 have thus been restricted to _ gs 0 < 0.05. The agreement between the schematic-MCT model and the available experimental and simulation data is surprisingly good, although this may be fortuitous.
Independent of quantitative details, there are two asymptotic regimes for s N ( _ g) that are predicted by schematic MCT. As _ g / 0, the residual stress approaches the yield stress s y ¼ lim _ g/0 s ss ( _ g) > 0. This is also indicated in the experiment (cf. open and lled circles in Fig. 5). For large pre-shear rates, on the other hand, s ss is predicted to reach a second plateau as _ g / N. Experimental data for the PS-PNIPAM suspension are compatible with these predictions. But note that the high-shear limit entails physics of particle contacts and solvent-induced interactions that are not captured in the present MCT. For the cases where s N grows with increasing _ g, this growth is much slower than that of the corresponding ow curve. The latter (dotted lines in Fig. 5) approaches a high-shear Newtonian regime (s ss $ _ g for _ g / N) in the schematic-MCT model. As a result, aer the cessation of stronger shear ow, a larger amount of stress is released than for weak ow.
The fact that the residual stress s N approaches a constant for both _ g / 0 and _ g / N can be understood on the basis of the ITT formula and an asymptotic consideration of the schematic-MCT correlation functions. For this, note that the ITT-MCT description is based on transient correlation functions f(t, t 0 ) that are modied solely by deformations occurring between their two time arguments t 0 # t. This is in distinct difference to the waiting-time dependent correlation functions F(t + t w , t w ) measured in nonequilibrium since there, the nonequilibrium distribution function at time t w determines the ensemble average, while for the transient correlators it is always the (deformation-independent) equilibrium distribution function.
In the integral determining s(t) aer cessation, eqn (2), only the regime t 0 < 0 t enters. Schematically, this is represented as the shaded area in Fig. 6. The MCT equation of motion for f(t, t 0 ) is dominated by the history-dependent contribution arising through the memory-kernel integral involving m(t, t 00 , t 0 ) and _ f(t 00 , t 0 ) for t 00˛[ t 0 , t], eqn (5). The paths for t 00 followed in this integral are indicated in Fig. 6 by dashed lines for the two functions. Fig. 7 shows the two-time "cessation correlator" f(t, t 0 ) for t > 0 and t 0 < 0 for an exemplary case.
The qualitative behavior of this correlation function can be understood from the structure of the memory-kernel equation and from the boundary conditions. Note that correlation functions are at least continuous, so that for t 0 ¼ 0 the cessation correlator equals the quiescent equilibrium one (cf. Fig. 6), the transient correlation function under steady shear rate _ g. These two boundary cases are shown in Fig. 7 as the red lines towards the back of the plot.
The MCT expression for s(t) suggests discussion of correlation functions at xed rst argument t, as a function of Àt 0 > 0. Exemplary results are shown in Fig. 8. Since the accumulated strain that decorrelates the MCT memory kernel increases with increasing |t 0 |, one expects the correlation functions to monotonically decay from their initial value for each t with increasing |t 0 |. (The situation may be different when, say, considering cessation of a large-amplitude oscillatory ow, where the transient correlation functions are nonmonotonic. 43 ) For liquid states (3 < 0 in the model), the equilibrium correlation function decays to zero on a structural-relaxation time scale s a . Hence, for t [ s a , the transient correlator f(t, t 0 ) for t 0 < 0 vanishes, and there is no contribution to s(t) in eqn (2). This explains the observation shown in Fig. 3 that in the liquid, stresses relax back to zero aer cessation on the time scale s a , independent of the previous shear rate.
In the glass, the quiescent correlator decays to a nite plateau, the nonergodicity parameter f ¼ lim t/N f eq (t). For a sufficiently large t, the cessation correlator f(t, t 0 ), viewed as a function of |t 0 |, hence has an initial value f. To understand its decay to zero, note rst that the steady-state transient correlation function in the shear-melted glass features a decay to zero on a time scale s _ g $ 1/ _ g (up to a prefactor set in the model by g c ), as long as _ gs 0 ( 1. For _ g / 0, a scaling limit is approached where f ss (t) $ f ss ( _ gt) for t / N. In this scaling limit, the short-time decay of the correlator on the time scale s 0 becomes irrelevant, and the scaling function obeysf ss (s) ¼ f for s / 0. Since the steady-state stress s ss is then determined essentially only by the integral over f ss ( _ gt), a dynamic yield stress arises that is qualitatively given by v s f 2 .
For times _ gt 0 ( 1, f ss (t 0 ) (or eqiuvalently f(0, Àt 0 )) remains close to the equilibrium curve. The decay of f(t, t 0 ) for t / N can now be understood considering the structure of the MCT memory integral shown in Fig. 6. At large t and small |t 0 |, the integral is dominated by equilibrium-correlator contributions, and extends the latter. The contributions only change signicantly once |t 0 | ¼ O(s _ g ), since then the contribution to the   memory kernel arising from _ f(t 00 , t 0 ) (the horizontal dashed line in Fig. 6) starts to decay. Therefore, also f(t, t 0 ) decays to zero on a time scale s 0 . If t [ s 0 , the two-time correlator hence resembles the scaling curvef ss ( _ gt) that determines the yield stress. This is visualized in Fig. 7 where we show f(t, t 0 ) as a function of the two time arguments t > 0 and t 0 < 0. At t 0 , t / 0, the microscopic relaxation from unity to f can be seen. For the small pre-shear rate chosen in the gure, this microscopic relaxation has a small weight in the integral determining s(t). For a large |t 0 |, the correlator becomes effectively independent of t, featuring decay on the time scale s _ g . This explains that for _ g / 0, both the residual stress s N and the steady-state stress s ss approach the dynamic yield stress.
For large _ gs 0 , the shear-induced relaxation of the steady-state correlator no longer obeys s _ g $ 1/ _ g. In the schematic model, it instead approaches s _ g / N $ s 0 . As a result, s ss f _ g, and the owcurve exhibits a high-shear Newtonian regime. The twotime cessation correlator f(t, t 0 ) is still governed by the same qualitative argument as given above: the initial value for t / N is still f, and the memory integral is still cut off by the decay of the steady-state correlator on a time scale O(1/ _ g) (the horizontal dashed line in Fig. 6). Thus, s N is still determined by the integral over a function whose decay scales as 1/ _ g, resulting in a constant expression in eqn (2). The numerical value of s N will differ somewhat from the yield stress s y , since the shape of the cessation correlator for t / N changes somewhat with changing _ g. The situation is demonstrated for two different shear rates, _ gs 0 ¼ 10 À4 and _ gs 0 ¼ 10, in Fig. 8 (solid and dashed lines). Here, the cessation correlator f(t, t 0 ) is shown for xed t > 0 as a function of rescaled second time argument, | _ gt 0 |. For the lower shear rate, the nal decay to zero is set by s _ g for all t, which is also the case demonstrated in Fig. 7. For _ gs 0 [ 1, the steadyshear correlation function decays on the microscopic time scale s 0 , and hence does not scale with _ g. As however t increases, the microscopic contributions disappear from f(t, t 0 ), and a decay on a shear-induced time scale, s 0 _ g $ 1/ _ g re-emerges. In this case one notices that s 0 _ g > s _ g , since the memory integral picks up microscopic transients that are, in the case of strong shear, no longer negligible. As a result, the irreversible plastic relaxation described by the transient correlation functions describe a residual stress that increases with increasing _ g, as noted above, but becomes independent of the pre-shear rate as _ g / N. So far, we have discussed constant-t cuts of the transient correlation functions, i.e., those that determine the shear stress within ITT-MCT. On discussing two-time-dependent correlation functions, it is oen more convenient to consider cuts for the xed second argument, i.e., constant-t 0 cuts, since this facilitates the interpretation of t 0 as a waiting time from which correlations into the future are measured. For the case of shear cessation, the transient correlation functions are of interest only for t 0 < 0. Fig. 9 shows an exemplary case for the schematic-MCT model. The qualitative behavior can again be understood from Fig. 7. For |t 0 | z 0, the quiescent correlator is recovered, while for |t 0 | / N, one approaches the steady-shear correlation function. Since this steady-state correlator decays on the time scale s _ g , the region of nontrivial t 0 -dependence is | _ gt 0 | ( 0.1. By continuity, the transient correlator then follows the steady-state correlator initially (for t < 0). Since for t > 0, no further relaxation mechanisms contribute to its decay, the transient two-time correlator then quickly crosses over to a t 0 -dependent plateau where f(t, t 0 ) is independent of t. In Fig. 7, this is exhibited by the fact that for a |t 0 | large enough such that the transient correlation function assumes a value less than the plateau f, it is independent of t in the glass. The curves shown in the main panel of Fig. 9 hence obey f(t, t 0 ) /f (t 0 ) < f for t / N and xed t 0 < 0. In the non-ideal glass, one expects this plateau to ultimately decay due to relaxation processes that are not captured in the MCT approximation. In the liquid, the curves are modulated by the nal decay on the structural-relaxation time scale s a . This is demonstrated in the inset of Fig. 9, where we show the same constant-t 0 cuts as in the main panel, but for the liquid state with separation parameter 3 ¼ À0.01. The transient cessation correlators are not straightforwardly accessible in molecular-dynamics simulations, since they require averaging over the equilibrium distribution function while the trajectories correspond to those inuenced by the past shear rate. In the simulation, the averaging is performed over a set of congurations that one assumes to be representative of the statistical ensemble as it evolves over time. This a priori only gives access to the waiting-time dependent correlation function F(t + t w , t w ). Evaluating this quantity within ITT-MCT involves additional approximations (see ESI of ref. 1) and is outside the scope of the present discussion.

Conclusions
We have analyzed the transient decay of the shear stress aer the cessation of steady shear ow in an ideal-glass model using schematic models of the mode-coupling theory combined with the integration-through transients framework (ITT-MCT). The theory predicts the partial relaxation of stresses to a long-time plateau value: a ow-history dependent residual stress s N emerges that is sustained by the glass produced from the shearmelted initial steady state. Brownian-dynamics computer simulations of a two-dimensional hard-disk mixture in the glass conrm this picture qualitatively. Our simulation results further corroborate similar results from previous experiments and 3D molecular-dynamics simulations. 1 The residual stress depends on the preparation history of the glass, in our case exemplied by its dependence on the preshear rate _ g. Two variants of the schematic-MCT model have been discussed that differ by whether or not they take account of anelastic strain-induced decorrelation of the overlap between stress and density uctuations. Taking into account only the relaxation of stresses through irreversible decorrelation of density uctuations caused by thermal noise (model A), residual stresses remain that are larger than the dynamical yield stress s y . Further stress relaxation is provided by reduction of stressdensity couplings that originate in the affine shear advection and is also the cause of pronounced stress overshoots in startup ow. In the schematic model, it needs to be captured empirically (model B), guided by a numerical discussion of the full, wave-vector dependent ITT-MCT. 37 The yield stress is the smallest stress the shear-melted glassy state needs to maintain in order to ow homogeneously. In the limit of innitesimally slow ow, none of this stress can relax. For faster ow, the steady-state stress rises above the dynamic yield stress, and for models with no or sufficiently weak startupstress-overshoot phenomena, these additional stresses only relax partially, due to the effect of irreversible thermal motion. A further stress-relaxation mechanism is provided by the same mechanism that is responsible for stress overshoots in ow startup: the capability of cages to transiently store a large amount of elastic energy that is released under plastic deformation. In the case of shear cessation, this mechanism allows for a further relaxation of past-ow-induced stresses, so that the remaining residual stresses decrease with increasing shear rate. Nevertheless, a certain amount of stresses never relaxes in the arrested glass, signalling the nonequilibrium nature of the amorphous solid that is being produced.
From the point of view of statistical physics, the appearance of persistent residual stresses is a clear deviation from Onsager's regression hypothesis, and hence a genuine nonlinearresponse effect. Consider Onsager's reasoning in the present context: 50 if an external perturbing eld h 0 coupling to the dynamical variable X, is switched on adiabatically in the innite past, h(t) ¼ h 0 exp[3t] (with 3 / 0+), the normalized deviation in X obeys: hDX(t)i ne /hDX(t ¼ 0)i ne ¼ F x (t) + O(h 0 ), where the relaxation function, F x (t), is the normalized correlation function of the uctuations of the variable X in the unperturbed system. It is thus independent of the external eld. In our case, the imposed shear ow with the rate _ g couples to the shear stress, yet we nd that the normalized deviatoric shear stress s(t)/s ss sensitively depends on small shear rates, even in the asymptotic long-time limit. The origin of the violation of Onsager's result lies in the time-dependence of the initial stationary state. Its decay is shear-induced and thus becomes arbitrarily slow in the limit of _ g / 0. In ITT-MCT this slow decay is the origin of a dynamic yield stress in the steady state, and of the persistent residual stress aer ow cessation. The ITT-MCT theory describes time-dependent single-point averages (such as the macroscopic stress s(t)) in nonequilibrium through Green-Kubo-like integrals, based on two-point transient correlation functions f(t, t 0 ). The latter are a convenient tool for theoretical analysis, but one has to keep in mind that they do not correspond to the correlation functions accessible in experiment or simulation. Since they are formed with the Boltzmann equilibrium distribution by denition, they describe the nonequilibrium evolution contingent on only those external elds that are present between their two time arguments t 0 and t. In ordinary correlation functions, averaging is performed with the time-dependent nonequilibrium distribution function at time t 0 , which is then denoted as a waiting time t w . The nonequilibrium distribution function contains a further ow-history dependence, so that the correlation functions F(t + t w , t w ) seen in experiment and simulation depend also on perturbations at times t 00 < t w . The ITT formalism can in principle be extended to describe these waiting-time dependent correlators, as has been demonstrated for the case of the MSD at t w ¼ 0 aer cessation. 1 We leave a more in-depth discussion of waiting-time dependences for a future publication. 51 Since glassy states contain internal stresses caused by past perturbations, they are no longer characterized solely in terms of the thermodynamic state variables. One expects them to show material properties that again depend on the past ow history; for example, the Maxwell shear modulus depends on the prestrain applied to the glass, 18,52,53 and also on the thermal history of the glass. 54 The time-dependent changes observed in the (frequency-dependent) dynamical shear moduli aer cessation provide a rheological probe of aging-like phenomena. 55 Understanding the dependence of material coefficients on the previous strain history is of great conceptual importance for applications. The further development of rst-principles microscopic theory will provide a basis for improving the predictability of materialsdesign processes.