Memory in Non-Monotonic Stress Response of an Athermal Disordered Solid

Athermal systems across a large range of length scales, ranging from foams and granular bead packings to crumpled metallic sheets, exhibit slow stress relaxation when compressed. Experimentally they show a non-monotonic stress response when decompressed somewhat after an initial compression, i.e. under a two-step, Kovacs-like protocol. It turns out that from this response one can tell for how long the system was in a compressed state, suggesting an interpretation as a memory effect. In this work we use a model of an athermal jammed solid, specifically a binary mixture of soft harmonic particles, to explore this phenomenon through in-silico experiments. Using extensive simulations under conditions analogous to those in experiment, we observe identical phenomenology in the stress response under a two-step protocol. Our model system also recovers the behaviour under a more recently studied three-step protocol, which consists of a compression followed by a decompression and then a final compression. We show that the observed response in both two-step and three-step protocols can be understood using Linear Response Theory. In particular, a linear scaling with age for the two-step protocol arises generically for slow linear responses with power law or logarithmic decay and does not in itself point to any underlying aging dynamics.

Memory in the context of physical systems refers to the ability to encode, access, and erase signatures of past history in the state of a system [1].This ability is absent when the system is in thermal equilibrium and consequently it is tied to out-of-equilibrium behavior.The origin and process of memory formation have been an intense area of research in last decade.In the context of jamming, the origin of memory in hard sphere glasses has been linked to a Gardner transition [2] very recently.Researchers have also explored memory formation in disordered media through directed aging [3,4] or simply through periodic driving [5,6].In structural glasses, memory effects have also been studied [7] through temperature cycles, which have helped to unveil the complex hierarchical structure of the free energy landscape; see Ref. [8] for a discussion of recent advances in the study of memory and rejuvenation effects in structural glasses.
A remarkable experimental protocol introduced by Kovacs [9] for the analysis of polymer glasses revealed that the time-dependent evolution of a glassy system can depend sensitively on its thermal history.A prime example is the non-monotonic evolution of a macroscopic observable (e.g.volume) after rapidly cooling an initially equilibrated sample, allowing it to relax for some duration and then warming it up instantaneously to a higher temperature; this has since then been dubbed "Kovacs effect" and has been analyzed theoretically for different models of glasses [10][11][12], also within the framework of Linear Response Theory (LRT) [13].
A generalization of the protocol described above has received increasing attention over the last few years, whereby the control parameter is generally different from temperature and the dynamics can be athermal.The systems that have been considered in this context are quite diverse, ranging from foams to crumpled metallic sheets and even jammed glass beads [14][15][16][17][18][19][20].These systems typically show logarithmic stress relaxation when compressed.This slow relaxation can be understood as a collective effect coming from a broad distribution of relaxation rates of the system (see [21] for details).The common effect observed in all of those systems is the nonmonotonic relaxation of the relevant response variable associated with the control parameter.To be specific, in a two-step protocol, where a compression is held for a waiting time t w and then the system is decompressed, a non-monotonic response of the pressure is observed after the second perturbation, i.e. for t > t w .It turns out that the time t p where the peak in the response occurs grows with the waiting time t w .In Refs.[14,17] this scaling is reported as linear.In a recent study [19] a three-step protocol (which adds a final compression after the two previous step perturbations) was explored in athermal systems and found to produce a more complicated behaviour, with the response exhibiting two extrema (one maximum and one minimum) instead of one.Nonetheless, the position of the maximum again scaled with the first waiting time, therefore showing some similarity with the two-step protocol as regards memory effects.
In spite of a substantial research effort aimed at understanding such memory effects and non-monotonic responses, the field still lacks a simple particle-based simulation model.In this paper we show that a well-known athermal model system [22] efficiently and accurately recovers the phenomenology of most relevant experimental scenarios, i.e. one-step, two-step and three-step deformation protocols, using simple numerical simulations.To be specific we recover the logarithmic stress relaxation when the system is compressed and the non-monotonic stress response for the Kovacs-like two-step protocol.Our simulations also exhibit a linear scaling of t p with t w to leading order, as well as further aspects of the phe- nomenology recently observed under three-step protocols.We emphasise the usefulness of having such an insilico experimental system as the simulations are easily reproducible and one can in principle measure every dynamical and structural aspect, which could be very challenging or infeasible in an actual experiment with a real athermal solid like a foam or crumpled metallic sheet.
Remarkably, we are able to rationalize our results solely using Linear Response Theory, using an approach similar to [23] but of course with the relevant response functions, which for step strain perturbations show an initial elastic response followed by a slow relaxation.An important theoretical insight will be that these response functions are time translation-invariant, which means that the linear memory effects we see here are, in spite of their superficial similarities, quite distinct from memory caused by an underlying aging dynamics as in glasses [9][10][11][12] We emphasize finally that we use the term LRT here in its broad sense of a linear relation between response and the corresponding perturbation.Such a linear response is generally expected for small perturbations, but this is not a consequence of the Fluctuation-Dissipation Theorem here: the latter is inapplicable for our athermal systems as there are no spontaneous fluctuations whose correlations could be measured.This paper is structured as follows.In section I we discuss the model (subsection I A) and the different perturbation protocols (subsection I B).In section II we set up the Linear Response formalism.In section III we show the results from simulations for each of the protocols (one-step, two-step and three-step) and we compare those with the predictions from LRT.Finally, in section IV we conclude with a summary and discussion.

A. Model
To model an athermal solid we use a binary mixture of soft particles interacting through a pairwise harmonic interaction (as used in Ref. [24,25] and more recently in Ref. [22], see Fig. 1 for a typical snapshot).The particles interact only when they overlap; explicitly the interaction potential is where σ ij is an additive interaction diameter computed as Here R i and R j are the radii of particles i and j, which can be either R A = R or R B = 1.4R depending on whether the particle concerned is of A-type or B-type (cf.Fig. 1).We fix our units by choosing for simplicity the force constant k = 1 and the radius scale R = 1.The time-dependent control parameter h(t) determines how compressed the system is and Θ(•) is the Heaviside function.We simulate N = 10 4 particles in a two-dimensional periodic square box with side length L = 216, with equal numbers N A = N B = N/2 of A-and B-particles.For the dynamics we assume a standard athermal overdamped equation of motion with implicit friction against a stationary solvent, where r i is the position of particle i and ∇ i = ∂/∂r i .We fix the drag coefficient ζ = 1 from now on.This, combined with the unit choices for R and k, then also sets the time unit ζ/(kR) = 1.To solve the particle dynamics given by Eq. 2 numerically we use an Euler scheme with fixed time step ∆t = 0.01.The packing fraction in our system is defined as All dynamical quantities are averaged over 128 independent simulation runs unless explicitly stated otherwise, starting as in Ref. [22] with randomly distributed positions of the particles in the simulation box.We first prepare the system at area fraction φ = 1 (corresponding to h = 1) and run the dynamics for a long time (t = 10 6 ) to guarantee that the system reaches an energy minimum, i.e. a fully relaxed state.These minimised configurations are then used as the initial conditions (with t reset to zero) for all perturbation protocols considered below.

B. Protocols
To explore the memory-like response we deform the system using three different protocols.These are described below and shown schematically in Figs. 1 and 2. w ), the value h2 is maintained for a time t (2) w after which a final change of h(t) to h3 is performed.
The protocols are named according to the number of times in which the control parameter h(t) is changed instantaneously.

One-Step Protocol
In the one-step protocol the system starts in the fully relaxed state associated with h 0 = 1.0 and then at time t = 0 this parameter is suddenly changed to h 1 > h 0 (right panel of Fig. 1).This causes a compression of the system: the increase in h increases the effective particle sizes and so leads to a higher area fraction as the box size remains the same.
The conjugate variable to h is the virial pressure, which is the negative diagonal component of the virial stress tensor where F i the net force on the i-th particle.From symmetry Σ xx = Σ yy ; numerically we therefore use only Σ xx to evaluate P (t).In our simulations the relaxation of this quantity is tracked as a function of time.This protocol provides the relaxation curves (see Fig. 3) that will be used for the Linear Response prediction, as discussed below in section II.

Two-Step Protocol
In the two-step protocol one starts again with the fully relaxed system at h 0 = 1.0.The control parameter h(t) is switched at t = 0 from h 0 to a higher value h 1 , as in the one-step protocol.The system is then allowed to relax in that state for a waiting time t w .At this point one makes a further change by decompressing the system, reducing h(t) from h 1 to a new value h 2 at t = t w ; we consider h 0 < h 2 < h 1 (see Fig. 2(a) for a sketch).Similar to the one-step protocol we monitor the pressure P (t) as an observable.The slow decay of P (t) up to time t w , which is the same as in the one-step protocol, is "interrupted" by the sudden drop in area fraction at t = t w .We analyze the evolution of P (t) after this point through direct numerical simulation and also using LRT (see section III).

Three-Step Protocol
Finally in the three-step protocol one proceeds as in the two-step protocol, with the time during which h(t) is held at h 1 now denoted t (1) w .After the second step, h(t) is left constant at h 2 for a further time interval t (2) w and at time w the parameter h(t) is switched for a third time, from h 2 to h 3 .We assume again that this value lies between the two previous ones, h 0 < h 2 < h 3 < h 1 (see Fig. 2(b) for a schematic).Hence the system remains at w .The pressure P (t) is again monitored after the final perturbation through direct numerical simulation and compared with LRT (see section III).
In real experiments, particle sizes -or equivalently system volume at fixed particle size -cannot be changed instantaneously.Nonetheless these changes can often be performed rapidly enough so that internal relaxation processes can be neglected while they take place.The instantaneous changes of h(t) that we work with constitute the idealized limit of such a scenario.

II. LINEAR RESPONSE PREDICTION
LRT is a well grounded framework to analyze memory effects in thermal systems [13,[26][27][28].In particular, Ref. [23] illustrates the theory for Markovian athermal systems, as is the case here.For a general timedependent perturbation one has the general linear response form [27,29] where δP (t) = P (t) − P (0), ḣ denotes the time derivative and χ(t) is called the step response function (also known as time-dependent susceptibility).This function is obtained from the one-step protocol and is used to predict the response for the multi-step protocols as we show shortly.FIG. 3. Linear response regime.After complete relaxation at h0 = 1.0 (P (0) = 0.01563) the scale factor h is changed to a new value h1 at time t = 0 following the one-step protocol (see Fig. 1).The resulting response function δP (t)/δh (see Eq. ( 6)) is plotted as a function of time (markers: data, dashed line: double logarithmic fit (Eq.( 9))); results for different h1 overlap almost completely, demonstrating linear response.Fit parameters: a = 0.2125, b = −0.00533,c = 0.00462, t0 = 176.21.
One-step protocol.For this protocol we set h(t) = h 0 + Θ(t)δh where δh = h 1 − h 0 .Substitution of this in Eq. ( 8) retrieves In LRT this response function is independent of δh, the size of the perturbation.Two-step protocol.For small changes h 1 = h 0 + δh 1 and h 2 = h 1 − δh 2 with |δh i /h i | 1, it follows from Eq. ( 8) that the pressure for t > t w is given by the superposition of the relaxation function evaluated at different times, that is Three-step protocol.In addition to the previous steps we consider here a further perturbation Using again the linear superposition principle Eq. ( 8), the pressure for t > t w ) + δh 3 χ(t − t (1)  w − t (2)  w ) (8) in terms of the same relaxation function χ(t) as before.

III. RESULTS
In this section we show that our simulation results recover all relevant observations for multi-step protocols (two and three-step) as seen in recent experiments [14][15][16][17][18][19][20].We first demonstrate that our simulations lie in the FIG. 4. Evolution of the pressure P (t) for a two-step protocol with tw = 10 2 and h1 = 1.002, h2 = 1.001.Dashed lines correspond to the evolution for the one-step protocol with the associated values of h.Inset: Non-monotonic response after the final (second) perturbation, showing pressure normalized by its value at t + w := tw + 1, i.e. just after the second perturbation; the time of the maximum of the response tp and its corresponding normalized height hp are indicated.
Linear Response limit, and then we show that using LRT we can correctly predict the observed phenomenology via Eqns.( 7) and ( 8).
In addition we make an asymptotic prediction for the scaling with t w of the time of the maximum of the response t p and its corresponding normalized height h p for the two-step protocol (see Fig. 4) using LRT.For this purpose we first fit the step response from the one-step protocol to a double logarithm.This fit was already suggested in the analysis of experimental data (see Ref. [14]).We further improve the fit by including the effects of the eventual saturation of the response function.Finally, we show that according to our numerical results we do not need an extra prefactor 'C' in front of t w (in Eq. ( 8)) to fit our response function in the framework of LRT.
One-step protocol.After the perturbation at t = 0 we observe a sudden pressure increase that comes from the elastic response of the amorphous solid, followed by a slow monotonic decay of the pressure.In Fig. 3 we show that the observed response falls within the linear regime.This is established by the collapse of the ratio δP (t)/δh for different values of h 1 .This ratio corresponds to the step response (Eq.( 6)).We then fit it to a double logarithmic curve of the form This fit is correct in a restricted time window: Fig. 3 shows that it becomes inaccurate at small times and Fig. 7 demonstrates that it also fails for very long times; see also the discussion around Eq. (11).Note that as long as we can measure and evaluate the linear step response χ(t) in some way, such a fit is not required to be able to use Eq.(7) or Eq. ( 8).However, a closed form expression  10)), dash-dotted line: asymptotic prediction with first order correction (Eq.( 11)).(c) Dependence of hp − 1 on tw (markers: data, solid line: linear response prediction).
for the step response allows analytical predictions as we will discuss shortly.
Two-step protocol.After the sudden drop at t w , again due to the elastic response to the step change in h(t), the pressure P (t) rises for some time and then decreases, eventually merging into the asymptotic curve from the one-step protocol for the same final h 2 .This can be seen in Fig. 4, where the solid line shows the non-monotonic response for the two-step protocol while the dashed lines correspond to the relaxations one finds for the one-step protocol with final parameters h 1 and h 2 respectively.The location of the maximum of the two-step relaxation curve depends strongly on the age of the sample t w ; intuitively, then, the system seems to "remember" how long it has been kept at the compression set by h 1 (see Fig. 5 for a quantitative comparison).Such memory effect has been observed in experiments on various athermal systems [14][15][16]19].
We can use Eq. ( 7) and the data for χ(t) derived from the simulation of the one-step protocol via Eq.( 6) to predict this nonmonotonic response (dashed lines in Fig. 5(a)) as a function of t w .As the figure demonstrates, the prediction from LRT is very accurate, clearly supporting the validity of LRT for our simulation setup.We remark that this prediction is made directly from the measured one-step response and does not require e.g. a parametric fit of this function.Such a fit is necessary, however, to predict analytically the scaling of t p and h p with t w (solid lines in Fig. 5(b) and Fig. 5(c)).Here we define t p as the length of time between t w and the maximum in the response.To find this, one solves for the maximum of Eq. ( 7) and subtracts t w .For the double logarithmic fit in Eq. ( 9) this yields a complicated but closed form expression that can be inserted back into Eq.(7).
Following [14] we define the height of the maximum h p by normalizing the result by P (t + w ), the value of the pressure just after the time t w of the second perturbation.Explicitly we use P (t w + 1.0); the unit time difference is chosen here because it is the smallest time for which the logarithmic fit works well (see Fig. 3).The resulting expressions for t p and h p are rather long, but for large t w (compared to the fit parameter t 0 , which is of order 10 2 ) one recovers a single logarithmic step response function (from Eq. ( 9)) for which the scaling can easily be found.In this case the solution for the position of the maximum of Eq. ( 7) is given by This relation was derived in Ref. [14] from another perspective and is recovered here within the LRT.Two main points can be made from the result in Eq. ( 10): (i) the asymptotic prediction for the scaling of t p is linear in t w , (ii) the scaling factor is independent of any fitting parameters of the step response function and only depends on the relative strength of the two perturbations.For our data the relative ratio δh2 δh1−δh2 = 1 and it predicts that t p = t w (see dashed line in Fig. 5(b)).As can be seen from Fig. 5(b), however, corrections to this leading order asymptotic prediction remain visible within the time window considered in the simulations.
The asymptotic prediction for the scaling of t p with t w can be improved by considering a more realistic fit for the step response χ(t), which incorporates the physical expectation that the step response function should approach a constant (which can be zero) at long times.The simple fit in Eq. ( 9) does not have this property in the generic case b = −c.We choose then the modified fit function χ(t) = ã + b log(t) 1+c log(t+ t0) . In Appendix V A we show that this function provides an adequate fit of the data for the step response for intermediate and long times (Fig. 7).The asymptotic prediction based on this fit (worked out in Appendix V B) is This agrees with the leading order term in Eq. ( 10) but also contains a relative correction of order 1/ log(t w ).For our chosen set of parameters it quantitatively evaluates to the dash-dotted line in Figure 5(b) and gives a significantly improved prediction for t p that matches the numerical data in the long waiting time regime.Three-step protocol.In the three-step protocol we have two control parameters for the timing of the perturbations; t w .Depending on the combination of these waiting times a non-monotonic evolution of the pressure P (t) can again be observed (see Fig. 6(a),(b)).
In the first scenario we keep t w constant and vary t (2) w .The non-monotonicity of the response is evident for small values of t (2) w ; this then crosses over to a monotonic decay for long values of the second waiting time (see Fig. 6(a)).Still the LRT works perfectly well (dashed lines in Fig. 6(a)).The linear response prediction is obtained by again using the response measured for a single step, together with Eq. ( 8).
In a second scenario we keep t w constant and vary t (1) w instead.Here we observe two extrema in the response after the final perturbation, a minimum followed by a maximum.The maximum shifts to larger time differences as we increase the first waiting time t w , which is reminiscent of the behavior we saw for the two-step protocol.The position of the minimum, on the other hand, remains constant (see Fig. 6(b)).Again, the linear response prediction works well without any fitting (dashed lines in Fig. 6(b)).In this case analytical closed form expressions cannot be obtained for the dependence of t p and h p on t (1) w , not even for asymptotically large t Nonetheless we can of course solve numerically for the maximum of Eq. ( 8) to arrive at the predictions shown in Fig. 6(c).The qualitative behaviour of t p and h p is the same as in the two-step protocol: both increase with the waiting time t (1) w , thus the system also exhibits memory effects under this three-step protocol.

IV. DISCUSSION
In this work we demonstrate that a model athermal solid made of a binary mixture of soft (harmonic) particles can serve as a canonical model for understanding the non-monotonic response seen recently in a diverse type of athermal systems subjected to multi-step variations of the control parameter.Although the particular model we study has been used before to understand the rheology [24,25] and aging of athermal systems [22], we have deployed it here for the first time to explore and understand two-step and three-step non-equilibrium experiments similar in spirit to the Kovacs protocol.
Our results show that this athermal system exhibits memory-like effects that can be understood from Linear Response Theory.Indeed, as long as we are in the regime of small perturbations, the only ingredient needed for the LRT is the step response function (Eq.( 6)).If the decay of this function is slow enough (compared to e.g. an exponential) a non-monotonic response is expected for a two-step or three-step protocol.This follows from the structure of Eq. ( 7) and Eq. ( 8) (the same observation was also pointed out in Ref. [28]).We note that we have checked the generality of this statement by considering a particle model with a different, Hertzian interaction.We found that it also exhibits a logarithmic decay of the pressure response to compression, and would therefore expect qualitatively the same phenomenology under multi-step protocols as for the harmonic repulsion considered above.
For the Kovacs-like two-step protocol, the way in which memory is displayed is via a relation between t w and t p (or h p ).Using for the susceptibility the logarithmic fit function in Eq. ( 9), as proposed previously for experimental data in Ref. [14], an asymptotic relation for large t w can be derived (Eq.( 10)).This relation is linear and independent of the parameters used in the logarithmic fit of the step response function χ(t).This linear scaling has already been observed in athermal systems and its origin lies in the logarithmic behavior of the step response function (see for instance Refs.[14,17]).
Here we have derived this scaling on the basis of Linear Response Theory, which constitutes a new result in itself.
Two comments are in order as regards the above observations.Firstly, one may wonder how generic the linear scaling of t p with t w is.We show in appendix V B that this linear scaling applies surprisingly broadly, i.e. whenever the step response function decays as a power law or a power of a logarithm.Secondly, we expect that for long waiting times the system actually reaches a steady state so that the response will reach a limiting value; a double logarithmic fit as in Eq. ( 9) therefore has to be an approximation that applies only for a finite time window, e.g. the one that is observable in the simulations.On this basis we proposed a fit function that takes into account the saturation of the response; for this a subleading correction to the linear scaling can then be derived, see Eq. (11).We showed that this modified fit function does an appreciably better job in terms of prediction of our simulation data (see Fig. 5(b)).
In previous studies, the non-monotonic response observed in the two-step protocol was explained via the assumption of a broad distribution of relaxation times [14,19].This provides intuition for the initially unexpected non-monotonicity: consider a system (as for instance a glass [21]) that is composed of slow and fast elements (with respect to experimental timescales), and that is subjected to the two-step protocol.Directly after the second perturbation, the system remembers -via the slow elements -its state before the perturbation, while the fast elements adapt rapidly to the new value of the external parameter.This creates a situation in which the two types of elements relax in opposite directions.The fast elements relax first, causing the unexpected increase in the response (in our case: pressure) with time.Eventually the slow elements then turn the relaxation around, giving rise to non-monotonic dynamics.This picture, though physically appealing, is arguably not needed to w as shown in legend.For (a) and (b), solid lines: direct numerical simulations, dashed lines: linear response prediction.The pressure has been normalized by its value just after the last step perturbation, i.e. at (t w + 1.0.(c) tw-scaling of tp and hp for the data in (b) (markers: simulation data, solid line: linear response prediction).The correlations between the small statistical fluctuations appearing at very long times ∼ 10 3 − 10 4 in subfigures (a) and (b) arise from the fact that the same set of initial configurations was used for the simulations with different t deduce the explicit formula in Eq. ( 7): the only ingredient required here is linear response, and we would argue that the earlier approaches were effectively using LRT, at least implicitly.To be specific, Eq. (3) from Ref. [14] can be recast in terms of the Linear Response formula (Eq.( 7)) with the susceptibility a single logarithmic decay function.On the other hand, it is easy to check that LRT with an exponential relaxation function cannot produce a non-monotonic response in a two-step protocol; what is needed is a slow relaxation, e.g. of power law or logarithmic type.In this sense the LRT framework helps to sharpen the physical picture of slow and fast elements into quantitative conditions on the relaxation function.
Overall, the memory effect we observe, which has also been reported for a range of different of systems [14][15][16][17][18][19][20], can be explained using linear response as long as the perturbations are small.A natural question is what happens beyond the linear regime.As shown in Appendix V C, qualitatively the same phenomenology is present in the non-linear regime.This indicates that the memory effect is a generic response of the system under the chosen protocols and not merely an artifact of linear response.A quantitative analysis of the t p scaling in this regime is left for future work.
Finally, it is worth discussing to what extent the physics of non-monotonic response to two-step or threestep protocols is related to memory and aging.As mentioned above, a linear response picture with a slow response function can be interpreted in terms of the dynamics of slow and fast modes: the slow modes can then reasonably be said to provide 'memory' of the history of the system, even though this is in the form of a weak (linear) perturbation of the system.In the non-linear regime discussed above, the different perturbations can cause irreversible structural changes, which are even more intuitive carriers of memory (see Fig. 9 and Appendix V C for more details).
An interpretation in terms of aging may also seem tempting, given that the maximum position t p grows with the age t w of the system, in fact with a linear scaling that is often referred to as simple aging [30,31].However, given that our results are very well represented by a linear response theory with a response function that is time-translation invariant, i.e. depends only on time differences, the dynamics does not meet the conventional definition of aging, where correlations and response become genuine two-time functions [32].In this context it is interesting to note that in recent experimental studies [14,19] a constant factor C had to be introduced ad hoc in front of t w in Eq. (7) or in front of t (1) w and t (2) w in Eq. ( 8) to fit the results of two-step and three-step protocols.This could indicate genuine aging behaviour, or alternatively effects that go beyond the linear response regime.
We believe that in order to arrive at a full understanding of the type of memory analyzed in this paper for athermal systems, the gap between experiments and theory needs to be closed by simulation models that allow detailed inspection of the microscopic dynamics.We hope that our work paves the way for future work in this direction.An obvious avenue in that direction would be to extend our current study to initial conditions that are not yet arrested and lie in the athermal aging regime [22], to understand the extent and impact of true aging effects on memory and non-monotonic response.3, where now fit provided by Eq. 12 has been added (dotted line).In addition, the xaxis has been extended to show the fit prediction in the large time limit.In the derivation of Eq. ( 11) we introduced a modified fit function for the susceptibility, namely In Fig. 7 we show the fit of the simulation data to this functional form and its extrapolation to long times to compared with the double logarithm (Eq.( 9)).
For long times t t0 , one may approximate Eq. ( 12) as From this expression one easily reads off that the step response saturates at the value χ(t → ∞) = ã + b/c.This result is physically sensible: in contrast to a liquid (where the susceptibility decays to zero in the long time limit), an athermal solid can only relax part of the stress generated by a perturbation.For short times, both expressions (Eq.( 9) and Eq. ( 12)) are inconsistent with the data: as the material has a short time elastic response with a finite bulk modulus, χ(t = 0) should be finite.We ignore this point as it has no impact on the predictions for multi-step protocols that we discuss in the main text.
B. Generic linear relation between tp and tw for some non-exponential response functions In this section we consider different non-exponential relaxation (step response) functions and obtain the associated scaling of t p with t w in the large t w limit.
We start by considering a simple power law of the form χ(t) = a + b/t.Substitution of this in the response for the two-step protocol (Eq.( 7)) yields The position t * of the maximum of this function (P (t * ) = 0) is easy to determine; we subtract t w to arrive at with r := δh 1 /δh 2 .We can repeat these steps for a generic power law of the form χ(t) = a+bt −α with α > 0. Then one arrives at the same linear form as in Eq. ( 15) but with a modified prefactor For α = 1 we recover Eq. ( 16) as expected.
Another realistic non-exponential response functions is an inverse logarithm of the form χ(t) = A + B/log(t).This is the asymptotic limit of Eq. ( 12) as derived in Eq. ( 13) and yields the formula given in Eq. (11) in the manuscript as we are going to show.
In this case one cannot solve explicitly for the maximum of the two-step response P (t) (Eq.( 7)); instead the following implicit condition for t p is obtained: By extending the simple expression from Eq. ( 15) one can use as an ansatz for t p a linear term plus a correction term as the inverse of the logarithm of the age: and then one has to solve for the unknowns c 0 and c 1 .Substitution of Eq. ( 19) into expression given in Eq. ( 18) leads to It is then convenient to perform the change of variable t w = exp(l): Since we are interested in an asymptotic relation (for large t w ) we now perform a series expansion in 1/l, which yields up to first order 1 r Solving now for the unknowns c 0 and c 1 , we get first from the zeroth order term Notice that this is consistent with Eq. ( 17) for α = 0.
From the first order term (in 1/l) the coefficient of the leading correction term becomes Putting Eqs. ( 23) and ( 24) back into Eq.( 19) one then obtains Eq. ( 11) given in the main text.
As a further step we note that the previous result can be generalized to step response functions χ(t) = A + B(log(t)) −β with β > 0. Following the previous steps with the same ansatz (Eq.( 19)) one arrives at the following expressions for c 0 and c 1 : The linear term is the same as for the simple inverse logarithm (Eq.( 23)), whereas the subleading correction differs in the prefactor 1 + β; for β = 1 the results are of course consistent with Eq. ( 24) for β = 1.Finally one can generalize the above arguments to relaxation functions that decay as power laws with logarithmic corrections, χ(t) = a + bt −α (log t) −β .We omit the details and only note that in this case the ansatz of Eq. ( 19) still works, with c 0 given by Eq. ( 17) and c 1 a (rather complicated) function of α and β.

C. Non-linear Response
The non-linear regime is detected by varying the strength of the perturbation δh = h 1 −h 0 in the one-step protocol.In Fig. 8 we show the response δP/δh for different h 1 .The results suggest that the system enters the non-linear regime around δh = 0.002 where deviations from the linear response limit become clearly visible.
We first explored a special two-step protocol where h 2 = h 0 and h 1 = 1.1.In this case one finds that even though the perturbation parameter is brought back to the original value, the transient perturbation causes irreversible changes in pressure (see Fig. 9 (a)); the perturbation also causes nonzero particle displacements with a strongly non-affine pattern (see Fig. 9 (b)), a clear indication of irreversible changes in the configuration.
To explore the memory phenomenology in the nonlinear regime, we have carried out the two-step protocol for different waiting times (t w ).The results show a similar scenario as in the linear regime, both in terms of nonmonotonic response (stress hump) and memory (linear dependence of t p on t w ); see Fig. 10.This also points towards the fact that the memory effect we observe is quite generic and not restricted to the range of validity of linear response theory.

FIG. 1 .
FIG. 1. (Left) Typical snapshot of a portion of the athermal binary mixture of harmonic particles.Compression and decompression are implemented by scaling the diameter of all particles by a time dependent scaling factor, h(t) (sketch top right).(Right) Schematic of the control parameter h(t) for the one-step protocol.

FIG. 2 .
FIG. 2. (a) Two-step protocol.At time t = 0 the control parameter h is suddenly switched from h0 to h1 and kept at that value for a time duration tw.Then it is changed again to the final value h2.(b) Three-step protocol.In addition to the steps in (a) (with tw → t (1)

FIG. 9 .
FIG. 9. (a) Pressure P (t) (in blue solid line) of the system after a special two-step protocol with large deformation h0 = 1.0, h1 = 1.1, h2 = 1.0 and tw = 2000; with these parameters the protocol brings the perturbation parameter back to its initial state.Nonetheless the pressure does not return to its initial value (red dashed line) before the perturbation.(b) Displacement field of particles (scaled by a factor of 10 for better visualisation) between the first perturbation and the last one clearly shows the irreversible changes in the particle configuration.