Stress relaxation in viscous soft spheres

We report the results of molecular dynamics simulations of stress relaxation tests in athermal viscous soft sphere packings close to their unjamming transition. By systematically and simultaneously varying both the amplitude of the applied strain step and the pressure of the initial condition, we access both linear and nonlinear response regimes and control the distance to jamming. Stress relaxation in viscoelastic solids is characterized by a relaxation time t * that separates short time scales, where viscous loss is substantial, from long time scales, where elastic storage dominates and the response is essentially quasistatic. We identify two distinct plateaus in the strain dependence of the relaxation time, one each in the linear and nonlinear regimes. The height of both plateaus scales as an inverse power law with the distance to jamming. By probing the time evolution of particle velocities during relaxation, we further identify a correlation between mechanical relaxation in the bulk and the degree of non-aﬃnity in the particle velocities on the micro scale.

Viscoelasticity is associated with one or more time scales that reflect the changing balance between viscous loss and elastic storage as a material's response to mechanical perturbations evolves in time. 1,2 Here we implement a standard rheometric test of viscoelasticity, namely stress relaxation in response to an instantaneous step strain, and apply it to a minimal numerical model for foams, emulsions, and soft colloidal suspensions. 3 Our focus is on athermal systems close to the nonequilibrium (un)jamming transition, where the material develops rigidity under compression. [4][5][6] Because the shear modulus vanishes continuously at the jamming point, weakly jammed states near the transition can be arbitrarily soft. 5 Intuition then suggests that their linear response window should also be narrow -small changes in strain amplitude should suffice to drive weakly jammed materials from linear to nonlinear response. Numerics confirm this expectation; under quasistatic shear, for example, the strain scales where the first contact change occurs and where bulk softening sets in both vanish as power laws with the pressure. [7][8][9] While there has recently been considerable interest in the nonlinear response near jamming, 7-21 the form of the relaxation time for large strain steps remains an important open question. Here we demonstrate for the first time that, as the system passes from linear to nonlinear response, relaxation times depend not only on the material constitution, but also on the amplitude of a shear perturbation.
A diverging relaxation time is an important mechanical property of soft amorphous matter near jamming. [22][23][24] In the jammed phase, the stress relaxation time t* describes the time needed to reach a new mechanical equilibrium after a sudden shear strain. 9 In the linear response, the divergence of t* as the confining pressure p is sent to zero signals the loss of rigidity. 23 The unjammed phase displays a similarly growing time scale, which marks a crossover from the power law to exponential stress relaxation. 25 Both linear and nonlinear stress relaxation can be characterized with the relaxation modulus G r (t,g 0 ), which describes the time evolution of the shear stress s(t,g 0 ) after a step strain with amplitude g 0 , For infinitesimal g 0 , the stress is directly proportional to the strain and G r is a function of time alone. In this limit the relaxation modulus is equivalent (i.e. related by standard mathematical transformations) to other common rheometric tests, including small amplitude oscillatory shear and flow start-up. 2 In the nonlinear regime this equivalence generally breaks down. Existing theoretical 23 and numerical 9 studies of G r near jamming are valid only in the linear response regime. Its typical form is illustrated in Fig. 1. After a brief plateau at short times, G r undergoes a power law decay before reaching a quasistatic plateau. The relaxation time is the time needed to reach the quasistatic plateau. In linear response it diverges as an inverse power law with p. 9,23,26,27 Numerical studies of nonlinear response near jamming typically neglect rate-dependent effects by focusing on quasistatic shear. [7][8][9]15,19 They have identified two important linearto-nonlinear crossover strain scales. The first corresponds to the breakdown of linear response on the scale of individual particle trajectories, which is driven by changes to the contact network. [7][8][9][10]12,14 The contact change strain scales as g cc B p 1/2 /N. [7][8][9] The second characteristic strain corresponds to softening, i.e. the loss of linearity in the average stress-strain curve. It scales as g † B p. 9,15,19 Note that these two strains scale differently with p; we will revisit this observation below.
In the present work we study the linear and nonlinear relaxation times of weakly jammed solids over a wide range in pressure and strain amplitude g 0 connecting linear and nonlinear response. Our central finding is that the relaxation time as a function of g 0 displays two plateaus: one in the linear regime, and a second, higher plateau at larger strains. The pressure dependence of these two plateaus is identical, i.e. they diverge as power laws with the same exponent. This is a surprising result, as there is no a priori reason for their exponents to be the same. We further relate the form of t* to the time evolution of floppy-like, non-affine particle motions during relaxation.

Methods and model
Foams are modeled with the Durian bubble model 3 in two dimensions. Bubbles are represented as disks that repel elastically when they overlap, with an additional dissipative force proportional to their relative velocity. The elastic force between particles i and j is proportional to their overlap d ij = (R i + R j ) À r ij , where R i and R j denote the radii and r ij is the length of the vector r ij , pointing from the center of particle i to the center of particle j, The viscous force depends on the relative velocity v ij of the touching particles evaluated at the contact, where t 0 is the microscopic relaxation time. All material properties are expressed in dimensionless units constructed from k, t 0 , and the mean bubble size. The stress tensor is where the Greek indices denote Cartesian coordinates. The contact stress term contains the total force at each contact, The inertial stress is dictated by the center of mass velocity v i . Each particle has unit density, so its mass m i is proportional to its area. V is the total area of the unit cell. The inertial stress term is negligible for times longer than the damping time t 0 .
Initial conditions are created by randomly populating the simulation box and then using an energy minimization protocol to quench instantaneously to a local minimum of the elastic potential energy at a fixed volume. The box is then allowed to undergo small changes in size and shape to achieve a target pressure p and zero shear stress -these are called ''shear-stabilized'' packings in the nomenclature of Dagois-Bohy et al. 28 The pressure provides a convenient measure of proximity to the (un)jamming point at p = 0. Packings are bidisperse to avoid crystallization; we use the standard 5,29 50-50 mixture of small and large particles and a radius ratio of 1 : 1.4. Once the initial state is prepared, we use molecular dynamics simulations to apply shear, which allows us to resolve the time evolution of the system. Newton's laws are integrated using a velocity Verlet algorithm.

Stress relaxation at finite strain
In order to describe the mechanical relaxation of soft sphere packings, we investigate the system's shear stress in response to an instantaneous step strain of amplitude g 0 applied at time t = 0. The strain is imposed using Lees-Edwards periodic boundary conditions while displacing the particles' coordinates (x i , y i ) affinely according to (x i , y i ) -(x i + g 0 y, y i ). In order to stay clear of any spurious periodic signatures in our results, we restrict applied strains to g 0 o 0.01; this is still large enough to observe the softening crossover for the highest pressure we simulate, as discussed below. For times t 4 0 after the instantaneous shear, the periodic boundaries are kept fixed in their strained position and the particles are allowed to relax to a new mechanical equilibrium. The resulting stress relaxation is illustrated in Fig. 1, which shows the relaxation modulus G r (t,g 0 ) as a function of time t for a single strain amplitude and pressure.
The relaxation modulus displays several noteworthy features. There is an initial plateau at times shorter than the damping time t 0 1, which occurs because viscous forces inhibit the Fig. 1 The time evolution of the shear relaxation modulus G r , calculated for a step strain with amplitude g 0 = 10 À6 at pressure p = 10 À4 and N = 1024. The characteristic relaxation time t* is identified as the point where G r reaches 1 + 1/e times its quasistatic plateau value.
system from relaxing at a rate faster than 1/t 0 . On longer time scales, the shear modulus decreases as a power law 1/t y with an exponent y = 1/2. 9,23 This relaxation continues until the stress reaches a second, long time plateau. The height of the plateau defines a quasistatic modulus G(g 0 ), which approaches the linear elastic shear modulus G 0 = G(0) in the limit of vanishing strain amplitude. The crossover between power law relaxation and the quasistatic plateau defines the relaxation time t*(g 0 ). Fig. 2a illustrates the evolution of the relaxation modulus with increasing strain amplitude at a pressure p = 10 À4 , which is representative of the entire range of pressures simulated here. All curves show qualitatively similar time evolution. However, there is a crossover with increasing g 0 . For the small values g 0 = 10 À6 and 10 À5 (solid and dashed curves), the relaxation modulus collapses, which is indicative of linear response. For higher strain amplitudes, beginning here around g 0 = 10 À4 , the entire curve shifts downward. This is strain softening. Softening is also evident in the quasistatic modulus G(g 0 ), estimated from G r (t = 10 6 ,g 0 ), which we plot in Fig. 2b for varying pressures (symbols). At low strains the modulus remains constant, consistent with linear response. Softening corresponds to a subsequent decrease in G(g 0 ) with increasing g 0 . This general trend is evident at all pressures.
Strain softening has been explored previously in ref. 9, where it was found that the onset of softening occurs at a strain scale proportional to p, after a finite fraction of the particles have undergone contact changes. There shear was built up incrementally using a quasistatic protocol, so that the final amplitude g 0 was reached via a large number of small steps Dg. Once the linear response has broken down and the system has begun to soften, however, there is no fundamental reason that the result of an incremental quasistatic protocol should correspond to the long time limit of viscoelastic relaxation after a single large step strain. It is therefore surprising that when we overplot the results for incremental strain from ref. 9 (lines), we find near perfect agreement between the two data sets. This suggests that, on average, the two protocols reach the same minimum in the energy landscape of the sheared system.

Relaxation time and strain dependence
We now investigate the time t*(g 0 ) needed to reach the quasistatic plateau after a strain of amplitude g 0 . While linear response can be accessed with careful numerical experiments, 9,23 one would prefer to have a complete characterization of the dependence of the relaxation time, not just on p, but also on the amplitude g 0 of the strain step. Our main result is the observation of a plateau in t* at large g 0 , with pressure dependence comparable to the relaxation time in linear response.
We identify t* as the time when the relaxation modulus reaches a value 1 + D times its value in the long time plateau. In the following we set D = 1/e; we have verified that our results are representative of a range of values for D. We simulate relaxation time measurements for stress relaxation over three decades in pressure, p = 10 À2 Á Á Á10 À5 , and four decades in strain amplitude, g 0 = 10 À6 Á Á Á10 À2 . Results are averaged over at least 500 realizations per condition. In MD simulations the total simulated time is limited by the available computational resources; especially for the lowest pressures and largest strain amplitudes, one can ask if the system might relax yet more at longer times. To exclude this possibility, we have also performed quasistatic simulations using the FIRE algorithm 30 to determine the long time limit of the shear modulus. We then recalculate the relaxation time using the quasistatic plateau value, in combination with the time evolution of the MD simulation. These results are in good agreement with the relaxation times calculated directly from MD. Hence we are confident that our results are representative of fully relaxed packings.
The evolution of the relaxation time, plotted in Fig. 3, can be separated into three stages. At low strains, the response is linear and the plot of t* versus g 0 plateaus, with the height of the plateau determined by the pressure. Next there is a second, intermediate regime, where linear response breaks down and Fig. 2 (a) The time evolution of the shear relaxation modulus G r for p = 10 À4 and N = 1024 at different strain amplitudes (see legend). (b) The quasistatic, long time shear modulus G as a function of strain. The data points show the long time response to instantaneous step strains. The lines are results from a separate set of simulations that reach the same total strain via a series of incremental steps applied using a quasistatic shear protocol. the relaxation time begins to grow with increasing strain amplitude. The crossover causes the relaxation time to increase by approximately one order of magnitude. Finally, there is a regime at comparably high strains where t* develops a second plateau. This trend continues throughout the studied pressure range, with the crossover shifting to higher strains with increasing pressure. As a result, the linear response window is at the edge of the sampled strain range for the lowest pressures, while the nonlinear plateau is only beginning to develop for the highest pressure.
In order to highlight pressure dependence, we seek to collapse the relaxation time data by plotting t*/p l versus g 0 /p n . We select l = 0.85, which is the relaxation time exponent identified numerically in our prior study of strictly linear response, 9 and close to the theoretically predicted value of 1. 23 For the strain axis rescaling we select n = 0.5, which is characteristic of the contact change strain scale g cc discussed above. [7][8][9] This choice is motivated by comparing Fig. 2b and 3, where one observes that the upturn in t* for increasing g 0 always occurs at a strain where the quasistatic shear modulus is still approximately flat, i.e. before the onset of softening. The rescaled data, plotted in Fig. 4, show good collapse over the entire range of strains and pressures. There is a small departure for the lowest pressure (i.e. closest to jamming) at the highest strain amplitudes, which may be associated with finite size effects.
The data collapse in Fig. 4 indicates one of our central results, namely that the relaxation time plateaus at low and high strains diverge as inverse power laws with p, with the same characteristic exponent l. We consider this result surprising, as there is no fundamental reason that the divergence of the relaxation times at finite strains should comply with the form for infinitesimal strain. The rescaling of the strain axis with p 0.5 , and the position of the crossover at a value g 0 /p 0.5 B O(1/N), strongly suggest that the increase in relaxation time is associated with the onset of contact changes, and therefore the breakdown of linearity in the particles' trajectories. We have verified that a plot with g 0 /p on the x-axis produces a significantly worse collapse, and also that reducing the system size shifts the crossover to higher strains.

Relaxation and non-affine particle motion
When jammed solids are sheared, particles primarily slide past their contacting neighbors, rather than interpenetrating. 11,23,26,31 This ''floppy-like'' motion is a precursor of true floppy modes, or zero frequency, non-rigid body eigenmodes, that appear below the jamming transition. Floppy-like motion is the physical origin of non-affine fluctuations. During floppy-like motion, relative displacements are predominantly perpendicular to the bond vector r ij pointing from the center of particle i to the center of particle j, not parallel to it. Floppy and nonaffine motion is well understood in linear elastic response. 22,31 However, little is known about how these displacements evolve in time, and/or in nonlinear response. Here we study the time evolution of the relative velocity of contacting particles during linear and nonlinear stress relaxation.
In order to analyze particle motions during relaxation, it is convenient to decompose each relative velocity v ij into longitudinal and transverse parts according to where the longitudinal velocity v 8 is parallel to the r ij direction, and the transverse velocity v > is along t ij = r ij Â ẑ, defined with respect to the unit vector ẑ pointing out of the plane. By construction the particles have zero velocity at t = 0, and they approach a new static state at long times. During the relaxation process we follow the full statistics of the longitudinal and transverse velocities. Fig. 5 shows the probability distribution functions (PDF's) of |v 8 | and |v > | for one pressure p = 10 À3 and several times, presented in units of the relaxation time t* (see legend). For both longitudinal and transverse velocities, the distribution grows as PDF B v for small v. The tails at large v are approximately exponential for short times t/t* { 1. At longer times the  distributions decay slower than an exponential and faster than a power law. Attempts to fit a stretched exponential do not yield a good fit. Rescaling velocities by their average value hv 8,> i at each time provides an approximate collapse for times t 4 t*, although some scatter remains. Due to this rough collapse, in the remainder we focus on average quantities, namely on the root mean squared (rms) velocities v rms  Fig. 6 for pressure p = 10 À3 and strain g 0 = 10 À4 , averaged over an ensemble of 100 packings. Note that v > is substantially larger than v 8 at all times, indicating that transverse motion is always dominant. After reaching their peak value at a time on the order of t 0 , the velocities steadily decrease as the packing relaxes, until eventually they drop sharply and simultaneously due to a fraction of the packings that fully arrest. This drop occurs long after the relaxation time, which is of the order t* B O(10 3 ) for this value of p and g 0 . Our interest here is primarily in the relaxation time t*, so in the remainder we focus on data at times prior to the drop.
To further assess the character of the particle motions at finite time, we introduce the ratio of rms velocities G measures the balance between motion that leads to sliding versus interpenetration. The value of G is of order unity for an affine velocity profile, while significantly larger values of G indicate strongly non-affine motion. In the following, we demonstrate that the relaxation of 1/G is strongly correlated with the relaxation of bulk stresses. Fig. 7 depicts 1/G for three values of the pressure and three values of the shear strain, for time intervals 10 À3 r t/t* r 10. In all cases 1/G decays, indicating that non-affinity increases with time. For further comparison, we overplot the corresponding G r in each panel (dashed lines). There is an evident similarity in their decay profiles; this strongly suggests a correlation between the mechanical relaxation time t* and the relaxation of floppylike, non-affine fluctuations.
In order to further probe the correlation between stress relaxation and non-affine fluctuations, we investigate the time evolution of 1/G for three pressures and two values of the strain amplitude, as shown in Fig. 8. The first strain amplitude, g 0 = 10 À6 , is in the linear regime for all values of p, while the second, g 0 = 4 Â 10 À3 , is in the second plateau of t* in Fig. 4. For both low and high strain amplitudes, we find good data collapse when time is rescaled by t* and G is rescaled with 1/p 0.4 . While the exponent 0.4 is determined empirically, it is in reasonable agreement with scaling arguments based on known relations in quasistatic linear response, 27 which give G B 1/p 1/2 .
The data collapse in Fig. 8 is further evidence that the same physics governs the relaxation of non-affine fluctuations and stress, in both the linear and nonlinear regimes. The data of Fig. 7 and 8 together indicate a strong correlation between non-affinity at the micro scale, and stress response on the macro scale. They establish a microscopic interpretation of the relaxation time: it is the time scale beyond which floppylike sliding motion (and hence non-affinity) fully dominates particle motion. Fig. 6 Longitudinal and transverse velocities versus time for p = 10 À3 and g 0 = 10 À4 . Fig. 7 A comparison of the shear relaxation modulus G r (dashed curves) with 1/G (solid curves) at three distinct pressures and strain amplitudes (see row and column labels).

Conclusion
We have used stress relaxation tests to determine the relaxation time of jammed solids as a function of strain and pressure. For sufficiently low strains, the linear response is valid and the relaxation time approaches a plateau determined solely by the pressure. Close to jamming, the strains needed to access linear response are extremely small, and many experimental protocols are likely to probe nonlinear response even if care is taken to apply a small strain. Beyond linear response, contact changes accumulate leading to softening, and the relaxation time grows. We find a second plateau in which t* is approximately independent of strain. To within the precision of our numerical measurement, this second plateau diverges at jamming with the same exponent that characterizes linear response. The crossover is associated with the onset of contact changes, and hence the post-crossover plateau should be accessible experimentally. Rheometry and simultaneous particle tracking in bubble rafts [32][33][34] could also access measures of non-affinity. In order to relate t* to the microscopic properties of the system, we have studied the statistics of floppy-like, non-affine motion, characterized by the time-dependent ratio G of the rms longitudinal and transverse velocities between particles in contact. We observe a strong correlation between G and the relaxation of shear stress in time. We infer that t* can be understood as the time needed to observe a fully-developed non-affine response; once non-affinity has reached its maximum, the system's subsequent response is quasistatic.
There are several likely directions for future work. A natural question is whether the observed behavior of the relaxation time persists in D = 3 spatial dimensions. D = 2 is the upper critical dimension for the jamming transition, 35,36 so we do not anticipate qualitative differences. One can also ask how the relaxation time develops for larger strains, up to and including the yielding crossover to steady plastic flow, which occurs for strains on the order of 10%. 9 We speculate that there exists some strain scale beyond which the instantaneously applied step strain is tantamount to thermalization of the system, hence it may be possible to make connections to the late stages of relaxation after a temperature quench. 37 Finally, it is also interesting to ask if there is any relationship between the relaxation time studied here and the duration of rearrangement events in steady plastic flow. 38,39

Conflicts of interest
There are no conflicts to declare.