Recoverable strain in amorphous materials: The role of ongoing plastic events following initial elastic recoil Open Access

The concept of recoverable strain has a long history in the rheology literature. Reiner noted that, when stresses are removed from a body, part of the deformation will in general be recovered, and that this part of the deformation is elastic [1]. Weissenberg suggested that the stress in a flowing material can be expressed at any time in terms of a strain that relates the current state of the material to some reference state and that this reference state, which depends on the flow history, defines the recoverable state to which the material would ultimately deform if the stress were then removed [2]. In essence, then, recoverable strain is the amount of strain recovered once a stress is removed from a body, in the direction opposite to that in which the stress had acted.

In physical terms, recoverable strain has to date been understood as elastic in origin: polymer chains elastically stretched in the direction of an imposed stress show viscoelastic recoil after the stress is removed, for example. (In a material with a single Maxwell relaxation time scale, this recovery is predicted to occur instantaneously. A superposition of Maxwell modes with distinct relaxation time scales confers a finite recovery time [3].) Indeed, the term “elastic recovery” has often been used synonymously with “strain recovery” [3–6]. Any unrecoverable strain is instead attributed to irreversible plastic deformations.

The same distinction is also evident in more recent literature concerning the use of recoverable versus unrecoverable strain to characterize yielding [7–10]. Donley et al. [7] note that zero-stress recovery tests allow materials to recover their instantaneous “ground state,” enabling calculation of the relative amounts of viscosity and elasticity. Shim and Rogers [8] state that “the recoverable and unrecoverable [strain] components can be interpreted as the contributions from the viscoelastic solid and plastic properties respectively.” Kamani et al. [9] note that “the recoverable [strain] component is related to elastic processes, while the unrecoverable component is related to the plastic behaviour” and that “yield stress fluids change from being viscoelastic solids, where deformations are recoverable, to deforming plastically, where deformation is unrecoverable.” Lee et al. [10] state that “recoverable strain is elastic, while viscous properties are dictated by the rate at which strain is acquired unrecoverably.”

In experimental studies, strain recovery following an initial forward creep under an imposed stress, after the stress is subsequently switched off, has been observed in cold set gels [11], protein gels [12], fractal colloidal gels [13], hectorite clay [14], peptide gels [15], carbopol [16], and polystyrene aggregate concrete [17]. A high degree of recovery after a large deformation has been discussed as being an important property of tough hydrogels [18–21]. Shape memory effects in self-healing polymers are also attracting increasing attention [22], although here the recovery is typically induced by a change in temperature, or other control parameter besides stress or strain, which we do not consider here.

In this work, we study theoretically recoverable and unrecoverable strain within the class of materials that are variously (in different parts of the literature) referred to as soft glassy materials [23,24], amorphous materials [25], yield stress fluids [26], or elastoplastic materials [27]. These include dense foams and emulsions, glassy colloidal suspensions, jammed athermal suspensions and granular matter, as well as harder metallic and molecular glasses. We do so by performing simulations within the widely used SGR model [23,24], focusing, in particular, on the simple recovery protocol in which a material of age $t w$ is subject to the switch-on at some time zero of a shear stress of amplitude $σ 0$⁠, which is held constant until it is later set back to zero at some time $t stop$⁠, after which the strain recovery is observed as a function of the time $t− t stop$⁠.

After an initially fast recoil that is elastic in nature, we find, counterintuitively, that significant further slow strain recovery then occurs via ongoing plastic events in the reverse direction, opposite to that of the original shear. We provide a mechanistic understanding of this phenomenon in terms of the evolution of the SGR model’s population of elastoplastic elements. In particular, we show that the initial recoil, which elastically advects each element backwards in its trap, leaves in a state of negative stress those elements that had yielded plastically while the material was under load. The subsequent delayed plastic yielding of these elements with negative stress then leads to slow ongoing strain recovery. Because this yielding is of negative stresses, counter to the original shear, we call these “reverse plastic events.” Counterintuitively, therefore, those elements that had yielded plastically while the forward load was applied still make an important contribution to the material’s overall recovery after unloading. This has important consequences for constitutive modeling, because such behavior can only arise in a constitutive model that evolves a full distribution of local stresses (or multiple moments of such a distribution). A model that instead evolves only a single average stress will be unable to capture it. To illustrate this, we also simulate a simple fluidity model, which shares many features of SGR’s behavior, including a yield stress, rheological aging, and stress overshoot in shear startup. In evolving only an average stress rather than a stress distribution, however, it lacks the slow ongoing strain recovery predicted by the SGR model.

We show that the observation of slow ongoing strain recovery from reverse plastic events holds robustly across broad parameter regimes. Unexpectedly—although in relatively rare parameter regimes—we also find that the slow ongoing strain post switch-off does not always take place in the negative direction, counter to that of the originally imposed stress, but can sometimes continue to accumulate in the forward direction. The recovery process is then non-monotonic overall, with the initially backward elastic recoil followed by a (lesser) slow forward straining. This is reminiscent of recent observations of non-monotonic stress relaxation after straining [28–31], which have been attributed to complex material memory effects [32].

The notion of a distribution of local stresses in an amorphous material, rather than a single macroscopic average stress, has a well developed history in the statistical physics literature. Modeling approaches based on shear transformation zones [33] and lattice elastoplastic descriptions are based on it [27,34–36], as is the SGR model itself [23,24]. These approaches in turn draw on concepts originally developed by Spaepan [37] and Argon and Kuo [38] in the context of metallic glasses and bubbles rafts. More recently, frustrated local stresses have been recognized to play an important role in the elastically driven aging of soft solids [39], the memory of soft jammed solids to shear flow [40], the elasticity properties of amorphous solids [41], the directional memory of soft glasses [42], and the non-monotonic relaxation of shear stress after flow cessation [43]. A key contribution of this work is to recognize and explore the significance of a stress distribution in strain recovery after unloading and to emphasize that important effects such as non-monotonic recovery can be captured only in constitutive descriptions that indeed consider a full distribution of local stresses (or at least multiple moments of such a distribution), rather than a single average stress.

The paper is structured as follows. In Sec. II, we outline the SGR model [23,24], within which we shall study the phenomenon of strain recovery. The step stress protocol that we shall simulate is specified in Sec. III. In Sec. IV, we define our units, summarize the parameters of the model and protocol, and discuss the rheological quantities to be measured. Our results are presented in Sec. V. We set out our conclusions in Sec. VIII.

The SGR model [23,24] considers the dynamics of an ensemble of elastoplastic elements that explore a landscape of energy traps. Each element is notionally taken to represent a local mesoscopic region of a soft glassy material: a few tens of droplets in a dense emulsion, for example. Each element is assumed small relative to any macroscopic variations in the flow field, on the one hand, yet large enough to allow the definition of local continuum variables of shear strain $l$ and shear stress $kl$⁠, on the other hand. Here, $k$ is an elastic constant, assumed the same for all elements. The local strain variable $l$ is taken to describe the element’s state of elastic deformation relative to a locally undeformed equilibrium. Between local plastic yielding events, defined below, the strain of each element is assumed elastically to follow the macroscopically imposed strain rate $γ ˙$⁠, such that $l ˙= γ ˙$⁠.

Conversely, an imposed shear of rate $γ ˙$ will halt this aging process and rejuvenate the sample to a flowing state with an effective age set by the inverse shear rate, $1/ γ ˙$⁠. In its glass phase, $x< x g$⁠, the model predicts steady state flow curves with a yield stress $σ Y(x)$ in the limit of slow shear $γ ˙→0$⁠. For $x g<x<2 x g$⁠, the model predicts power law flow curves with a shear-thinning (sublinear) exponent, giving infinite viscosity in the limit of slow shear. For $x>2 x g$⁠, the flow curves are Newtonian.

In many soft glassy materials, the typical scale of energy barriers for local rearrangement events—the penalty for stretching soap films in a foam, for example—is far greater than the typical scale of thermal energy, $k BT$⁠. In such materials, the parameter $x$ is taken to be an effective noise temperature that models an assumed coupling with other yielding events elsewhere in the sample. This is intended to capture in a mean-field way the basic idea that a local yielding event in one part of the sample may trigger follow-on yielding events in other regions. To date, this remains a mean-field assumption that is yet to be made self consistent. In materials for which the energy barriers are instead on the scale of $k BT$⁠, the parameter $x$ can be taken as the true thermal temperature.

In this section, we define the recovery protocol that we shall simulate within the SGR model just described.

Sample preparation prior to shear is modeled by assuming a sudden deep quench at some time $t=− t w$⁠, from infinite effective temperature to the final working value $x$⁠, which is held constant thereafter. For times $− t w<t<0$⁠, the sample is allowed to age undisturbed, without any strain or stress being applied. Throughout most of the paper, we assume that all elemental strains are equal to zero during this initial aging stage, $l=0$⁠. An exception can be found in Figs. 3 and 4, where we assume an initial Gaussian distribution of $l$ values, of narrow standard deviation $l 0=0.05$⁠. This does not affect the basic physics of interest in this work, but merely allows us to show as a narrow Gaussian in Fig. 4 what would otherwise be an unplottable initial delta function distribution. We return at the end of Sec. V to consider in more detail the effect of non-zero $l 0$ values on the physics that we report.

At time $t=0$⁠, a step stress of amplitude $σ 0$ is imposed. At the time of this step, all elements instantaneously strain elastically forward within their traps by the same amount $l= σ 0/k$⁠, giving an instantaneous global elastic strain response $γ 0= σ 0/k$⁠. With the imposed stress held constant and equal to $σ 0$⁠, we then allow the sample to strain further forward by an amount $Δ γ f$⁠, beyond the elastic straining that occurred at the instant the stress was imposed. Note that this further straining arises from elements hopping between traps and accordingly is plastic in nature. Once the total strain has attained a value $γ(t)= γ 0+Δ γ f$⁠, we set the stress to zero, $σ=0$⁠. We define the time at which this occurs as $t= t stop$⁠. In principle, we could use either $t stop$ or $Δ γ f$ as a control parameter in our study. In practice, we have chosen to use $Δ γ f$ because it better reflects the physical state of the material at switch-off.

For all subsequent times, $t> t stop$⁠, we hold the stress equal to zero and measure the strain $γ$ as a function of the time interval $t− t stop$ since switch-off. We shall be interested, in particular, in the degree to which the strain recovers in the negative direction, whether it does so monotonically or non-monotonically, and in the basic physics that underlies this recovery process.

The parameters of the SGR model are the effective temperature $x$⁠, the attempt time for local yielding events $τ 0$⁠, the elemental modulus $k$⁠, and the glass transition temperature $x g$⁠. Recall Sec. II. The parameters of the rheological protocol are the time $t w$ for which the sample is aged before the stress is imposed, the amplitude of the imposed stress $σ 0$⁠, and the degree to which the sample is allowed to strain further forwards due to plastic yielding while the stress is held fixed, $Δ γ f$⁠. Recall Sec. III.

To solve the SGR model numerically, (at least) two different methods can be used. In the first, one derives integral equations from the governing master equation, Eq. (3), then uses numerical mathematics to solve these [24]. In the second method, which we adopt here, one instead directly simulates the stochastic hopping dynamics of a population of $M$ SGR elements. In any time step $dt$⁠, each element has its strain shifted elastically by an amount $l→l+ γ ˙(t)dt$⁠, where $γ ˙(t)$ is the strain rate at that time. Any element also plastically hops into a new trap with probability $rdt$⁠, with $r$ that element’s probability per unit time of hopping as given by Eq. (1). In the limit of many elements $M→∞$⁠, one recovers the solution of the governing equations. During any intervals in which the stress is held constant, the strain rate is calculated as $γ ˙=⟨lr ⟩ P$⁠. This ensures that the rate of stress loss via plastic hopping between traps is exactly counterbalanced by stress gain via elements shifting their local strains within their existing traps. At the instant of any step stress of size $Δσ$⁠, each element is simply shifted elastically in its trap by an amount $l→l+Δσ/k$⁠.

We work in units of stress in which the modulus $k=1$⁠, and of time in which $τ 0=1$⁠. We restrict ourselves to homogeneous shear, assuming that (for example) no shear bands form either during forward creep or backward recovery. With this simplification, there is no lengthscale implied by the flow geometry or model. We also set the glass transition temperature $x g=1$⁠. This amounts to setting the average trap depth in the landscape’s prior distribution to unity, thereby setting the typical scale of strain for local yielding events. To compare with any experimental data, therefore, all strains reported below would need to be rescaled by the typical strain for local yielding events in the material.

We set the noise temperature $x=0.3$⁠, deep within the glass phase. Results are presented for a number of SGR elements $M= 10 5$⁠, convergence checked against $M= 10 6$⁠. The numerical time step $dt=α/⟨|l|r(E,l)⟩$ with $α= 10 − 5$⁠, convergence checked against $α=2× 10 − 6$⁠.

Important physical parameters to be explored are then the sample age $t w$ prior to stress switch-on, the amplitude of the imposed stress $σ 0$⁠, and the degree of forward plastic strain that accumulates before the stress is switched off, $Δ γ f$⁠. Recall from Sec. IV that this further defines a switch-off time $t stop(x, t w, σ 0,Δ γ f)$⁠.

At the instant of stress switch-off, all elements instantaneously strain elastically backwards, giving a global strain change $Δγ=− σ 0$ (in our units). We then measure the subsequent strain evolution as a function of the time $t− t stop +$ since switch-off. Our interest, in particular, will be in the total reverse strain $Δ γ rec$ accumulated in the long time limit, $t− t stop→∞$⁠, beyond that which arose in the fast elastic recoil during the reverse step stress. As noted above, in rare cases, we shall find $Δ γ rec<0$⁠: the further straining that arises after the initial reverse elastic recoil, in fact, takes place in the forward direction.

We now present our results. As a preamble to understanding strain recovery post switch-off, we shall start in Sec. V A by analyzing the forward creep and/or yielding and flow that arise while the stress is held imposed. We then present our results for strain recovery in Sec. V B.

Following the application of a step stress of magnitude $σ 0$ at time $t=0$⁠, a material’s creep response is characterized by the strain curve $γ(t)$⁠. As noted above, in the SGR model, this shows an initially elastic step strain of magnitude $γ 0= σ 0/k(= σ 0$ in our units) at $t=0$⁠, followed by subsequent plastic straining for $t>0$⁠. Figure 1(a) shows this plastic part of the creep, beyond the initial elastic step, normalized by the stress amplitude $σ 0$⁠. Results are shown for a fixed sample age $t w$⁠, and several imposed stresses in curves upwards, spanning values from below the SGR model’s dynamic yield stress $σ Y$ to above it. (At the noise temperature $x=0.3$ considered in this work, the yield stress $σ Y=0.758$⁠.) The corresponding strain rate curves $γ ˙(t)$ are shown in panel (b).

For imposed stresses $σ 0< σ Y$⁠, the strain rate decreases perpetually as a power law $γ ˙∼ t − α$⁠. The exponent $0<α<1$ depends on the noise temperature $x$⁠. The strain correspondingly increases sublinearly as $γ∼ t 1 − α$⁠, corresponding to sustained power-law creep. In contrast, for imposed stresses $σ 0> σ Y$⁠, an initial interval of creep (for $σ 0$ only just above $σ Y$ at least) is interrupted by yielding, in which the strain rate curves upwards before finally settling to a time-independent value prescribed by the steady state flow curve $σ ss( γ ˙)$⁠, with $σ ss− σ Y(x)∼ γ ˙ 1 − x$⁠. The time delay before yielding scales as the sample age $t w$⁠, with a prefactor that diverges as $σ 0− σ Y→ 0 +$ [45].

Figure 2(b) shows the same data as in Fig. 1, but with the scaled strain rate $γ ˙(t)/ σ 0$ now plotted parametrically as a function of the forward plastic strain $(γ(t)− γ 0)/ σ 0$⁠. Accordingly, time $t$ is not explicitly represented in Fig. 2. The scaled strain on the abscissa, however, increases monotonically with increasing $t$⁠. The counterpart data for sample ages $t w= 10 1$ and $10 6$ are shown in Figs. 2(a) and 2(c), respectively. As can be seen, for imposed stresses above the yield stress, $σ 0> σ Y$⁠, yielding typically occurs at a scaled strain $(γ− γ 0)/ σ 0=O(1)$⁠, although the value increases slowly with $t w$ and diverges as $σ 0→ σ Y +$⁠.

A dot at each of several logarithmically spaced values of the abscissa in Fig. 2 shows the values of scaled forward plastic strain $(γ− γ 0)/ σ 0=Δ γ f/ σ 0$ at which we shall later, in Fig. 5, explore the amount of strain $Δ γ rec$ that is ultimately recovered plastically after the stress is switched off.

Having mapped out the forward creep and/or yielding and flow that arise while the stress is held imposed, we now explore the strain recovery after the stress is switched off. Our focus will be in particular on the part of this recovery that stems from slow ongoing plastic yielding of elements that were left with negative local stresses as a result of the fast initial elastic recoil. In Sec. VI, we shall consider the interplay of a non-zero time scale for the initial (visco)elastic recoil with the later slow plastic strain recovery.

The basic physics that we aim to elucidate is shown in Figs. 3 and 4. The imposed stress $σ$ as a function of time $t$ is shown in the top panel of Fig. 3. The corresponding strain response $γ(t)$ is shown in the bottom panel. As described above, a stress of amplitude $σ 0$ is switched on at time $t=0$⁠. (The preparation phase during which the sample is allowed to age undisturbed for a time $t w$ prior to $t=0$ is not shown.) At the instant of switch-on, the system shows an instantaneous forward elastic step strain of amplitude $γ 0= σ 0$⁠, in our units in which $k=1$⁠. (In experimental practice, this forward strain would be rapid but not instantaneous, with a rate limited by inertia or dissipation; we return to this point in Sec. VI.) The stress is then held constant for some time interval, over which the strain increases further due to plastic yielding events, this forward straining having been explored in Sec. V A. Once it has thus increased by an amount $Δ γ f$ beyond $γ 0$⁠, the stress is switched off. This defines the switch-off time $t stop$⁠. At the instant of switch-off, $t= t stop$⁠, the system shows an instantaneous backward elastic step strain, $Δγ=− σ 0=− γ 0$⁠. [We consider in Sec. VI the effect of a finite time scale for this initial (visco)elastic recoil.] The stress is then held at zero for all subsequent times $t> t stop$⁠. As a function of the time interval $t− t stop +$⁠, the strain further recovers in the negative direction by an amplitude $Δ γ rec$⁠, beyond the negative elastic step $Δγ=− γ 0$ that occurred at the instant of switch-off. This additional recovery $Δ γ rec$⁠, which has a significant amplitude compared with both $Δ γ f$ and $γ 0$⁠, arises from slowly ongoing plastic events, in which elements left with negative local stresses by the initial elastic recoil slowly reset these negative stresses to zero.

Consider the predictions of this equation for the imposed-stress protocol of interest here. Integrating across the positive (resp. negative) step in stress at $t=0$ (resp. $t= t stop$⁠) gives an elastic step strain of the same amplitude and sign as the step stress (in our units in which $k=1$⁠), as indeed discussed above and seen numerically in Fig. 3.

In order to fully understand the forward plastic straining during stress imposition and the strain recovery afterward; therefore, we must consider in more detail the evolution of the probability distribution of local strains, and the way this determines the strain rate $γ ˙(t)$ via Eq. (9). Accordingly, we plot this distribution $P(l,t)$ in Fig. 4 (having first integrated out the trap energy depth variable $E$⁠). We do so at several times $t$ while the stress is held applied (top panel in Fig. 4), and at several times after the stress is switched off (bottom panel). The time at which any distribution is plotted in Fig. 4 is indicated by a cross of corresponding color in Fig. 3.

The form of these distributions in Fig. 4 can be understood as follows. At the end of the sample preparation phase, immediately before the stress is applied, the distribution of local strains is (by assumption) a Gaussian of narrow width $l 0$⁠, centered at $l=0$⁠. (As noted above, our results in all figures except 3, 4, and 7(b) have $l 0=0$⁠, corresponding to a delta function.) At the instant the stress is switched on, this distribution shifts via the elastic advection term, without changing shape, to be now centered at $l=+ σ 0=+ γ 0$⁠. This distribution is shown by the black curve in Fig. 4. Over subsequent times while the stress is held imposed, plastic yielding events arise: recall that an element in a trap of depth $E$ with local strain $l$ has a hopping rate $r(E,l)$ given by Eq. (1). When any such plastic event occurs, the corresponding element resets its local strain $l→0$⁠. This slightly depletes the original spike in the distribution centered at $l= σ 0$ and leads to the development of a secondary lobe of probability around $l=0$⁠.

Were the strain held constant, with zero strain rate $γ ˙=0$⁠, the macroscopic consequence of these plastic yielding events would be a relaxation of the ensemble average stress at rate $σ ˙=⟨lr(E,l) ⟩ P ( t )$⁠. Conversely, to maintain the stress constant, $σ=+ σ 0$ with $σ ˙=0$⁠, requires a forward strain of corresponding rate $γ ˙=⟨lr(E,l) ⟩ P ( t )$⁠, ensuring sufficient forward advection of elements within their traps exactly to counterbalance this plastic relaxation. This is the origin of the forward plastic straining seen for times $0<t< t stop$ in Fig. 3. Its consequence is to advect forward at rate $l ˙= γ ˙$ both the spike in probability centered (originally) at $l= σ 0= γ 0$ (as it simultaneously depletes due to the plastic relaxation events just described), as well as the newly developing lobe of probability near the origin $l≈0$⁠. In this way, the distribution evolves over time from the black to cyan curve in the top panel of Fig. 4.

At the instant the stress is switched off, a fast elastic recoil occurs, in which the distribution of local strains instantaneously shifts by $l→l− γ 0$⁠, with $γ 0= σ 0$⁠, without changing shape, as each element elastically reduces its strain within its trap. Accordingly, the rightmost (cyan) curve in the top panel of Fig. 4, which obtained at $t= t stop −$⁠, shifts to become the rightmost (magenta) curve in the bottom panel at $t= t stop +$⁠.

Importantly, this distribution shown in magenta just after switch-off now has some weight for strains $− γ 0<l<− γ 0+Δ γ f$⁠, corresponding to those elements that had plastically yielded and reset their strains to zero while the stress was held applied, and then shifted elastically by $l→l− γ 0$ at switch-off. These elements then themselves plastically yield over the subsequent time $t− t stop$ since switch-off, resetting their strains from that negative strain interval $− γ 0<l<− γ 0+Δ γ f$ to $l=0$⁠. Accordingly, weight is progressively lost from the leftmost part of this probability lobe at negative $l$⁠, with probability instead accumulating nearer $l=0$⁠, in moving from the magenta to green curves in the bottom panel of Fig. 4.

If the strain were held constant, the macroscopic effect of these plastic yielding events from negative to zero $l$ would be a positive growth in stress, via Eq. (9), as elements hop from negative to zero local stress. Conversely, under the conditions of constant stress that pertain after stress switch off, $σ=0$ and $σ ˙=0$⁠, a corresponding negative strain rate is required to advect elements in the negative direction in their traps, in order to counterbalance this positive plastic stress change. Accordingly, as weight is lost from the negative $l$ part of the distribution in moving from the magenta to green curves, these curves also simultaneously advect leftwards. We, therefore, now understand the slow ongoing strain recovery that arises after the initial elastic recoil post switch-off to stem from the reverse local plastic yielding of elements left with negative stress values by that initial recoil. We further recognize these elements to be those that had previously yielded plastically in the forward direction while the material was under load. In this way, counterintuitively, stress relaxation that occurs plastically while a material is under load can still in fact contribute to strain recovery after the load is removed.

With this basic physics in mind, we shall now explore the dependence of the phenomenon on each of the three control parameters in our study: the amplitude of imposed stress $σ 0$⁠, the degree to which the sample plastically strains forward while the stress is held imposed $Δ γ f$⁠, and the age of the sample $t w$ before stress switch-on. We shall consider first in Fig. 5 the total strain $Δ γ rec$ recovered plastically at long times $t− t stop→∞$⁠, before examining the time-dependence of the recovery process in Fig. 6. Recall that $Δ γ f$ and $Δ γ rec$ are indicated on the left and hand sides, respectively, of Fig. 3(b): a value $Δ γ rec=0$ corresponds to zero plastic strain recovery (beyond the initial elastic recoil $Δγ=− γ 0$ at the instant of switch-off), whereas a value $Δ γ rec=Δ γ f$ would indicate full recovery.

Accordingly, we plot in Fig. 5 the scaled plastic strain recovery $Δ γ rec/Δ γ f$ as a function of the scaled forward plastic strain $Δ γ f/ σ 0$ that accumulated while the stress was imposed. Data are shown for imposed stress values $σ 0=0.1,0.2,…,2.0$ in curves upwards in each panel, for samples ages $t w= 10 1, 10 3$ and $10 6$ in panels downwards, corresponding to the counterpart forward plastic creep curves in Fig. 2. Indeed, dots marked left to right on each creep curve in Fig. 2 show the various scaled forward strain values at switch-off, $(γ(t)− γ 0)/ σ 0=Δ γ f/ σ 0$⁠, shown as corresponding data point dots on each curve in Fig. 5.

To understand the trends observed in Fig. 5 we recall that, immediately prior to stress switch-off, any elements that remain unyielded since the stress was imposed have local strain $l= σ 0+Δ γ f$⁠, as in the rightmost lobe of the cyan distribution in Fig. 4. (This actually has a slightly smeared distribution of strain values of small width $l 0$⁠. As noted above, this is simply for the purposes of being able to plot the distribution graphically and does not change the basic physics.) We shall call these “group I” elements. In contrast, those elements that yielded at least once during the initial forward creep (“group II”) have local strain values in the range $0<l<Δ γ f$⁠, giving the leftmost lobe of that distribution. At the instant of stress switch-off, the whole distribution of shifts to the left, with each $l→l− σ 0$⁠, corresponding to a fast elastic recoil. Immediately after switch-off, therefore, the unyielded group I elements have positive local strain $l=Δ γ f$⁠, as in the (slightly smeared) rightmost lobe of the magenta distribution in Fig. 4. They furthermore reside in traps of a depth set by the age $t w$⁠, and with a rate of yielding $r∼1/ t w$⁠. In contrast, the yielded group II elements have local strain values in the range $− σ 0<l<Δ γ f− σ 0$ immediately post switch-off, as in the leftmost lobe of the magenta distribution in Fig. 4. In being freshly yielded, they reside in shallow traps with a rate of yielding $r=O(1)$⁠. The relative proportion of elements in group II compared to those in group I increases with increasing forward strain $Δ γ f$⁠.

Strain recovery post switch-off then arises from a basic competition between the plastic yielding of elements in group I versus those in group II. In most parameter regimes, the yielding of group II outweighs that of group I. When these negatively strained group II elements yield, they reset their local strains to zero. As discussed above, to maintain the stress constant and equal to zero post-switch off, an exactly counterbalancing negative global strain—strain recovery —must arise.

Indeed, for small forward plastic strain values $Δ γ f< σ 0$⁠, all group II elements have negative local strains post switch-off. In this regime, the larger the imposed stress $σ 0$ prior to switch-off, the more negative overall will be the local strains of group II elements post switch-off, requiring more counterbalancing recovery. This explains the increasing levels of recovery with increasing imposed stress prior to switch-off at the left of Fig. 5. The larger the sample age $t w$⁠, the greater the dominance of group II over group I elements, which are stuck in deep traps with yielding rate $1/ t w$⁠. Accordingly, the degree of strain recovery increases with increasing $t w$ in moving from the top to the bottom panels of Fig. 5, for these small $Δ γ f/ σ 0$⁠.

A different trend is, however, observed for larger values of the scaled forward strain prior to switch-off, $Δ γ f/ σ 0$ just above 1. In this case, for large imposed stresses, the ultimate plastic strain recovery is actually observed to be negative: the system strains plastically further forwards as a function of time $t>0$ since switch-off. The overall strain recovery is therefore non-monotonic, with the initially backward elastic recoil at the instant of switch-off $t= t stop$ followed by this forward plastic straining for $t> t stop$⁠, reminiscent of observations of non-monotonic stress relaxation after straining [28–31]. The total recovery is still however negative, with the initial elastic recoil always exceeding any subsequent forward plastic straining, for all parameter regimes we have explored.

This counterintuitive observation of forward plastic straining post switch-off can be understood as follows. For these larger values of $Δ γ f σ 0$ just above 1, the local strain distribution $− σ 0<l<Δ γ f− σ 0$ of group II elements immediately post switch-off contains some weight at positive strain $l>0$⁠. Furthermore, group I elements are also strained further forwards in their traps, acquiring larger yielding rates. In consequence, as the time since switch-off elapses, the yielding of elements with positive local strains now outweighs that of elements with negative local strains. To maintain zero stress post switch-off, this must be exactly counter-balanced by an increase in strain, giving negative strain recovery. This effect is more pronounced in younger samples (small $t w$⁠) for which group I elements are in relatively shallow traps and yield more rapidly.

For large $Δ γ f/ σ 0$⁠, such significant forward straining takes place before switch off, with significant flow for imposed stresses $σ 0> σ Y$⁠, that the proportional degree of recovery post switch off is small. Accordingly, $Δ γ rec/Δ γ f→0$ as $Δ γ f/ σ 0→∞$⁠, as indeed seen in Fig. 5.

We show finally in Fig. 6 the strain recovery as a function of the time $t− t stop$ since switch-off. The six panels shown correspond to the smallest stress (left panel column) and largest stress (right panel column) and three different waiting times (increasing in panel rows downwards) for which the ultimate recovery in the limit $t− t stop→∞$ is plotted versus the scaled forward strain $Δ γ f/ σ 0$ in Fig. 5. Each panel in Fig. 6 shows the results for the time dependent strain recovery for several different values of forward strain $Δ γ f$⁠, increasing logarithmically in curves upwards, with the lowest curve in each panel of Fig. 6 corresponding to the leftmost datapoint in the counterpart curve of Fig. 5.

So far, we have assumed the distribution of local strains present in the sample at time $t=0$⁠, before shearing commences, to be a delta function of zero mean. (An exception can be found in Fig. 3, where we assumed a narrow Gaussian simply to facilitate plotting.) We have likewise assumed the distribution of strains selected by any element immediately post-hop to be a delta function of zero mean. We now consider the more general case in which the distribution before shearing commences is a Gaussian of zero mean and standard deviation $l 0$⁠, and the post-hop distribution is a Gaussian of zero mean and standard deviation $l p$⁠. (Our results so far therefore have $l 0= l p=0$⁠.) The parameter $l 0$ thus now represents the degree of local strain frustration present in the sample before shearing starts, arising from whatever complicated aging or annealing dynamics took place during sample preparation. Smaller values of $l 0$ correspond to better ordered, more strongly annealed samples. Likewise, $l p$ represents the frustration experienced by elements in attempting to relax their local strains during local rearrangement events.

In Fig. 7, we explore the effect of non-zero values of $l 0$ and $l p$ on the basic phenomenon of recoverable strain arising from reverse plastic events. Shown are results counterpart to those in Fig. 3, but now for a range of $l p$ values in panel (a) and $l 0$ in 3(b). As can be seen, the basic phenomenon is quite robust to varying these quantities. Indeed, the degree of strain recovery from reverse plastic events increases with increasing post-hop frustration $l p$⁠, as might be expected intuitively: the left hand lobe of the turquoise distribution (for example) in Fig. 4(a) will now be more smeared, including some weight to the left of the origin. When elastically shifted leftward to form the magenta distribution in (b), more weight will be present at more negative $l$ values, resulting in faster yielding of elements with more strongly negative $l$⁠, and more negative global straining (i.e., more strain recovery) to compensate. Although the effect of plastic strain recovery decreases with increasing initial frustration $l 0$ in panel Fig. 7(b), it nonetheless persists to values of this parameter $O(1)$⁠.

So far, we have used the SGR model in its original form, considering only the elastoplastic stress arising from the population of SGR elements. In this version of the model, the initial response to an applied step stress is instantaneous: elements shift elastically within their traps at the instant of the step. Accordingly, in Fig. 3 (say) we saw an instantaneous elastic recoil at the time the stress is switched off, followed by the later slow recovery arising from reverse plastic yielding of elements left with a negative local strain by the initial recoil. (The instantaneous elastic upward step during switch-on at $t=0$ is not seen on the logarithmic time axis in Fig. 3.)

In practice, however, drag against a background solvent will render the time scale for this first part of the strain response non-zero, such that the initial recoil is now viscoelastic rather than purely elastic. An obvious question is then to what extent the initial (visco)elastic recoil and later reverse plastic yielding of elements with negative local stresses during strain recovery after stress switch-off overlap in time, or whether they occur truly sequentially, with the fast (visco)elastic recoil simply acting to leave elements with negative local stresses that later relax to give reverse plastic events. To answer this question, we now consider the total stress $Σ=σ+η γ ˙$ to comprise the usual elastoplastic stress $σ$ as predicted by the SGR model in Eq. (5), plus a Newtonian solvent contribution of viscosity $η$⁠.

Figure 8 shows results for forward creep and reverse strain recovery with model parameter values otherwise as in Fig. 3, but now with a non-zero solvent viscosity $η$⁠. As can be seen, this indeed imposes a finite time scale for the initial rise in strain after stress switch-on [panel (a)] and for the initial (now visco)elastic strain recoil after stress switch-off [panel (b)]. The time scale for each of these processes scales as $η/k$⁠, tending to zero in the limit $η→0$⁠. It is in this limit that the initial elastic recoil and later reverse plastic events become truly distinct: for finite $η$⁠, the two processes merge somewhat.

Panels (c) and (d) show the number of plastic events $n$ accumulated per element on average during forward creep and reverse recovery. Consistent with plastic events being thermally activated (as well as strain-induced) in the SGR model (rather than purely strain-induced, as would be the case in an athermal system), these signals show only modest dependence on the strain (and so on $η$⁠). As can be seen, the plastic part of the forward creep and—importantly to the key message of this work—of the reverse recovery, arises in each case only once $n$ approaches $O(1)$⁠, denoting significant plastic yielding.

We have seen, then, that the initial (visco)elastic recoil post-switch off and the later slow plastic strain recovery will only be truly time-separated in the limit of zero solvent viscosity $η→0$⁠. In either case (for $η$ finite or zero), however, it is clear that these two stages interact, in the sense that the initial (visco)elastic recoil creates the state in which some elements have negative local strain values. It thereby creates the conditions out of which the slower plastic part of the strain recovery can later emerge.

The strain response of this simplified model to the stress protocol considered in this work is shown in Fig. 9. When the stress is imposed, we see an initial (visco)elastic loading on a time scale $O(η/G)$⁠, followed by a slower creep once the integrated rate of plasticity $n=∫dt1/τ$ approaches $O(1)$⁠, as in SGR. After stress switch-off, we see an initial (visco)elastic recoil on a time scale $O(η/G)$⁠, as in SGR. Crucially, however, this simplified model lacks a distribution of local stresses, having just one global averaged stress $σ$⁠. Accordingly, it lacks reverse plastic events and correspondingly lacks the slow ongoing strain recovery predicted by the SGR model: compare panel (b) of Fig. 9 for the fluidity model with the counterpart panel (b) of Fig. 8 for SGR.

This has important implications for constitutive modeling, in showing that the phenomenon studied in this work can be captured only by models that evolve a full distribution of stress values (or multiple moments of such a distribution), rather than a single sample-averaged stress.

In this work, we have studied theoretically recoverable and unrecoverable strain in amorphous and soft glassy materials such as dense foams and emulsions, glassy colloidal suspensions, jammed athermal suspensions, granular matter, and metallic and molecular glasses. Within the soft glassy rheology model [23,24], we have simulated the rheological protocol in which a material of age $t w$ is subject to the switch-on at time zero of a shear stress of amplitude $σ 0$⁠, which is later switched off at some time $t stop$⁠, and the strain recovery observed as a function of the subsequent time $t− t stop +$⁠.

We have reported the counterintuitive phenomena that a significant component of the strain recovery post switch-off arises, after an initially fast elastic recoil, via the slow reverse plastic yielding of elements left with negative local stresses by that recoil. We have also showed that the time evolution of strain post switch-off is not always monotonic but can be non-monotonic (albeit in relatively rare parameter regimes): the plastic straining that takes place for times after switch-off $t> t stop$ does not always occur in the negative direction, as in the initial elastic recoil at $t= t stop$⁠, but can sometimes in fact continue to accumulate in the same forward direction as the originally imposed stress. This provides a potential counterpart, in this imposed-stress protocol, to recent observations of non-monotonic stress relaxation after the switch off of an imposed strain [28–31].

We have discussed these observations in terms of the evolution of the SGR model’s population of elastoplastic elements. In particular, we have shown that plastic strain recovery (after the initial fast elastic recoil) arises via elements that locally yielded during the initial forward plastic creep, resetting their local strains from positive values to zero as they did so. These elements then acquire negative local stresses during the fast elastic recoil, and later yield again post switch-off, resetting these local stresses from negative values to zero. It is this plastic yielding from negative to zero stress values that leads to global plastic strain recovery, in order to maintain the global stress of the sample as a whole equal to zero post switch-off. In this way, elements that yield plastically in the forward direction while a material is under load can still contribute to strain recovery after load removal.

The phenomenon reported here has important implications for constitutive modeling, in showing that a full distribution of stress values (or multiple moments of such a distribution) is needed to capture it. Any model that evolves just a single average stress will fail to describe it.

The physics described in this work is reminiscent of but distinct from the concept of reversible plasticity, which has gained increasing attention in the recent physics literature. Indeed, a key concept in the rheology of amorphous materials is that of intermittent local plastic events, sometimes called T1 events, arising at soft spots or shear transformation zones [25], coupled by elastic stress propagation [46]. Within this framework, oscillatory shear experiments on jammed particle rafts identified three different regimes of imposed strain amplitude [47–49]. At low amplitudes, the material responds elastically, with no plastic events. At high amplitudes, irreversible plastic events occur within each cycle, with significant cycle-to-cycle stroboscopic changes in particle positions. More surprisingly, at intermediate strain amplitudes, these experiments report reversible local plastic events within each cycle. Such events dissipate energy, consistent with the notion of plasticity, but exactly reverse in the opposite part of the cycle, with zero stroboscopic cycle-to-cycle change in particle positions. Reversible plastic events have also been reported in experiments on foams [50], in particle simulations [50–54], and in elastoplastic models [55]. Reversible plastic limit cycles have likewise been studied in models of hysterons [56] and iterated maps [57]. For a review, see Ref. [58].

It should be noted, however, that the mean field SGR model as studied here does not provide any spatial resolution of elemental positions. Accordingly, we are not able to make any statement as to whether a plastic event at a specific location later exactly reverses at the same location. Furthermore, because plastic events in SGR are thermally activated as well as being strain-induced, exact reversibility under straining is unlikely, except in the athermal limit $x→0$—in which limit, however, predictions of sustained slow creep are difficult to obtain.

Therefore, we do not claim a close correspondence between the plasticity that arises in this work in the reverse direction during strain recovery post-switch off and the truly reversible plasticity in which individual elements exactly retrace their trajectories, as seen under oscillatory shear in athermal systems. Instead, we make only the weaker statement that some elements that yield during the initial forward creep later do so in reverse during recovery. For this reason, we have been careful to use the weaker term reverse plastic events rather than the term reversible plasticity, as coined in the literature concerning oscillatory shear in athermal systems. Indeed, we use the term “reverse” in this work to mean “in the opposite direction,” as distinct from “reversible” meaning “an event capable of being reversed.” Future work should study creep and recovery in spatially resolved athermal models of amorphous materials, to establish whether there is a stronger link to reversible plasticity, with its precise connotations of exact localized reversibility.

Related to this, we have assumed throughout that the shear remains macroscopically homogeneous across the sample, disallowing any possibility of shear band formation. However, it has been predicted that shear bands will tend to form in an imposed stress protocol at the end of any initial creep regime, as the strain rate dramatically increases and the sample yields [5]. This is likely to be an important effect for large values of $Δ γ f$ in particular. Future studies should consider extending this work to allow for the formation of shear bands.

The authors thank Mike Cates, Andrew Clarke, and Gareth McKinley for comments on an early draft of the manuscript and for helpful discussions. We also thank SLB (Schlumberger Cambridge Research Ltd.) for funding. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. 885146). We thank the anonymous reviewers for their insightful comments that improved the clarity of our descriptions.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.