Pre-yielding mechanical response near the jamming transition

Abstract

The mechanical and rheological properties of jammed packings of frictionless particles under shear strain remain not fully understood, even when the strain amplitude is very small and well below the yielding threshold. Systems above the jamming transition point φ_J are known to display two anomalous mechanical behaviors with respect to the driving frequency ω (or time t) and the strain amplitude γ. In the linear-response regime (γ → 0), the complex modulus exhibits an algebraic scaling, G(ω) ∼ ω^1/2 (or G(t) ∼ t^−1/2 in the time representation). In contrast, in the quasi-static limit (ω → 0), the modulus shows the nonlinear behavior, G(γ) ∼ γ^−1/2, a phenomenon referred to as softening. The ranges of ω and γ over which these algebraic scalings hold broaden as φ_J is approached from above, whereas both G(ω) and G(γ) vanish for φ < φ_J. In this study, we investigate the mechanical response in the regime where these two anomalies coexist in the vicinity of φ_J. To this end, we perform numerical analyses using two rheological protocols: oscillatory shear and transient stress relaxation. Our results demonstrate that the mechanical responses are not simply described as a superposition of the two algebraic relaxations and instead exhibit rich nonlinear viscoelastic behavior both above and even below φ_J.

---

1. Introduction

Disordered packings of athermal particles, such as emulsions, non-Brownian suspensions, granular materials, and foams, become rigid when their density, or packing fraction φ, exceeds the jamming transition point φ_J.¹ For frictionless spherical particles, various scaling laws and universal behaviors near φ_J have been established.^2,3 In particular, the algebraic dependence of the contact number and pressure on δφ = φ−φ_J, the emergence of diverging length scales, and the universal spectrum of the vibrational density of states have been identified in both simulations^2,4,5 and experiments^6,7 and explained using theoretical approaches.^8–10 These results characterize the properties of packings in the absence of applied stress.

Jammed systems also exhibit rich mechanical and rheological behaviors when subjected to external stress or shear deformation. Thus far, studies of the mechanical response near the jamming transition have primarily followed two main directions: steady-shear rheology and the response to small shear strains. In the former case, the main focus is on the relationship between shear stress σ and the shear rate [small gamma, Greek, dot above]. For example, the shear viscosity, η ≡ σ/[small gamma, Greek, dot above], of non-Brownian suspensions diverges as η ∼ |δφ|^−β below φ_J, while, above φ_J, the shear stress is well described by the Herschel–Bulkley law, σ = σ_y + C[small gamma, Greek, dot above]ⁿ. The positive exponents β, n, together with the yield stress defined as σ_y ≡ σ([small gamma, Greek, dot above] → 0), serve as key indicators of the critical behavior associated with the jamming transition.^11–26

The latter case, the response to small shear strains, is the main focus of this study. A jammed packing behaves as an elastic solid at very small γ, but upon increasing γ, it yields and enters the plastic regime at a finite strain γ_y.^27–29 It is known that, even in the pre-yielding regime, where γ is well below γ_y, jammed packings exhibit highly nontrivial mechanical properties, owing to the marginal stability inherent to jamming criticality, in contrast to other amorphous solids. In particular, two distinct anomalous algebraic scalings for the mechanical response function are observed. One is the scale-free frequency dependence of the viscoelastic complex modulus. As the system approaches φ_J from above, the complex shear modulus, G*(ω) = G′(ω) + iG″(ω), develops a power-law behavior, G*(ω) ≈ Aω^1/2 (hereafter referred to as the “ω^1/2 scaling”).^30–32 The ω^1/2 scaling is universally observed, regardless of system details, over a wide range of frequencies. This scaling is observed for ω > ω_s, where ω_s is the onset frequency that scales as ω_s ∝ δφ for the harmonic potentials (and ∝δφ^α−1 for other soft-repulsive potentials defined in eqn (1)). The ω^1/2 scaling is understood as a consequence of the development of the plateau in the vibrational density of states near φ_J³ and the boson peak, universally observed in amorphous solids.^5,33–35 A real-time manifestation of the ω^1/2 scaling is observed in the transient stress relaxation after a step strain, where the modulus decays as G(t) ≡ σ(t)/γ ∼ t^−1/2.^36–38 The other anomalous algebraic scaling is the nonlinear stress response in the quasi-static limit (ω → 0), known as shear softening. The static storage modulus G_qs ≡ G′(ω → 0), which is independent of γ in the linear-response regime, starts to exhibit nonlinear softening behavior, G_qs ∼ γ^−1/2, as φ_J is approached from above.^39–43 We refer to the softening as the “γ^1/2 scaling”, hereafter. The softening exponent 1/2 is insensitive to system details. The onset strain γ_s that marks the crossover from the elastic to the softening regime is proportional to δφ. Recently, we have found a scaling law that unifies the elastic, softening, and yielding regimes.⁴³

In this study, we investigate the interplay between the ω^1/2 scaling and γ^1/2 scaling of jammed packings in the pre-yielding regime both above and below φ_J. We address several key questions in this work. First, how do the ω^1/2 scaling (or t^−1/2 scaling) and γ^1/2 scaling coexist? Previous studies of the ω^1/2 scaling assume small γ, whereas studies of the γ^1/2 scaling generally focus on the quasi-static limit. However, as φ_J is approached, the two scaling regimes inevitably broaden and overlap. Therefore, in the vicinity of φ_J, the two scalings become simultaneously relevant, and their interplay is important. If both scalings for ω and γ coexist, should we expect a simple superposition of the two scaling functions of ω and γ, analogous to the Cox–Merz rule⁴⁴ or the strain-rate frequency superposition?⁴⁵

Second, how does the mechanical response above φ_J cross over to that below φ_J? G*(ω) in the linear-response regime (γ → 0) and G_qs(γ) in the quasi-static limit (ω → 0) are finite only above φ_J. Both are strictly zero below φ_J, as a persistent contact network is absent. If ω and γ are finite, the stress may arise even below φ_J, because particles intermittently collide and form transient contact networks. The formation of the contact networks of frictionless spheres under oscillatory shear below—but close to—φ_J has been extensively studied in the context of the reversible–irreversible (or absorbing) transition.^46–50 However, the relationship between this nonequilibrium phase transition and the nonlinear mechanical response remains largely unexplored.

Third, how is the stress relaxation of jammed packings influenced by nonlinear perturbations? The relaxation dynamics of amorphous solids from perturbed configurations toward a ground state, i.e., toward the bottom of the potential-energy landscape, are important for understanding the hierarchical structure of rugged energy landscape in the context of the glass transition.^26,51–53 Studying the γ-dependence of G(t) for various γ values would thus provide insight into the connection between the nonlinear rheology and the topographical features of the energy landscape.

Our study is the first step toward addressing these questions. We perform numerical simulations and investigate the nonlinear mechanical responses under time-dependent shear strains γ(t) in frictionless, athermal jammed packings near φ_J using different shear protocols. The first protocol is the oscillatory-shear measurement, designed to address the first and second questions. By fixing the frequency within the ω^1/2 scaling regime and varying the strain amplitude over a wide range, we measure the complex shear modulus G*(ω) both above and below φ_J. While many studies have attempted to unify the mechanical responses above and below φ_J in the post-yielding regime, γ → ∞,^12–16,23 little attention has been paid to the pre-yielding regime γ < γ_y. The second is transient stress relaxation, in which the time dependence of the stress after a step shear strain is monitored. It is equivalent to the oscillatory-shear measurement in the linear-response regime, but deviates from it once the linear-response assumption breaks down. We analyze the transient relaxation of shear stress, energy, force, and displacement following instantaneous strain deformation over a broad range of initial strain amplitudes.

This study is organized as follows. Section 2 describes the simulation model. Sections 3 and 4 present our results for the two protocols. Finally, we summarize our findings and conclude in Section 5.

2. Numerical modeling and methods

We consider a two-dimensional equimolar binary mixture of frictionless particles with diameters D_S and D_L, where the size ratio is D_L/D_S = 1.4, a value commonly used in studies of the jamming transition.^2,4 The interaction potential between particles j and k is given bywhere Θ(x) is the Heaviside step function, r_jk = |r_j − r_k|, and D_jk = (D_j + D_k)/2, with D_j and D_k denoting the diameters of particles j and k, respectively.¹¹ In this study, we focus on the harmonic potential with α = 2. The number of particles is N = 1156. We choose D_S, ε, and ε/D_S² as the units of length, energy, and stress, respectively.

We generate initial configurations via mechanical training using quasi-static and cyclic volume changes, following the procedure of Kawasaki and Miyazaki,⁴³ in order to stabilize the jammed configurations against shear deformation. During the quasi-static deformation, we minimize the potential energy using the FIRE algorithm (see Appendix A for the detailed setup of the FIRE algorithm).^54,55 We regard the system as mechanically equilibrated when the magnitude of the average force acting on a particle is less than 10⁻¹⁴ε/D_S. To eliminate residual stress in the initial configurations, we apply shear stabilization,^43,56 under Lees–Edwards boundary conditions.⁵⁷ We first prepare a random configuration at φ = 0.8395 and increase φ incrementally by Δφ = 10⁻⁴ while maintaining mechanical stability, up to the maximum packing fraction φ = 0.9. Next, we decrease φ by Δφ = 10⁻⁴ when and by Δφ = 10⁻⁶ when E < 10⁻⁸. We identify the packing fraction at which E first satisfies E < 10⁻¹⁶ as the jamming transition point φ_J.

The dynamics of the jammed configuration under shear deformation are monitored via molecular dynamics simulations. Throughout the study, shear deformation is applied in the x-direction. We assume overdamped particle dynamics with Stokes drag, so that the equation of motion is given by

where ζ is the friction coefficient. The unit of time is given by t₀ = ζD_S²/ε. Under shear with a finite shear rate, the particle velocity ṙ_j(t) in eqn (2) is replaced by ṙ_j(t) − [small gamma, Greek, dot above](t)y_j(t)e_x, where e_x is the unit vector in the x-direction. We numerically integrate eqn (2) using the Euler method.⁵⁷ For the quasi-static process (ω → 0), we use the numerical protocol of Kawasaki and Miyazaki:⁴³ particle positions are relaxed to the mechanical equilibrium at every incremental shear strain Δγ. For γ < 10⁻³, Δγ is increased logarithmically from 10⁻⁹ to 10⁻³; for γ > 10⁻³, Δγ is fixed at 10⁻³. The shear stress is defined as⁵⁷where L is the system length, x_jk = x_j − x_k, and y_jk = y_j − y_k. All data presented in this paper are averages over at least 15 (typically 50) independent realizations.

3. Oscillatory shear

Consider the mechanical response of a jammed system under an oscillatory shear strain defined as γ(t) = γ sin(ωt), where γ is the strain amplitude and ω is the frequency. In the linear-response regime or when γ → 0, the complex modulus G*(ω) = G′(ω) + iG″(ω) becomes independent of γ and follows the scaling relation

ΔG′(ω) = Aω^1/2, (ω > ω_s) (4)

where ΔG′(ω) ≡ G′(ω) − G_qs and G_qs = G′(ω) (ω → 0) is the static storage modulus. For the harmonic potential, G_qs scales as δφ^1/2, the onset frequency ω_s scales as ω_s ∝ δφ, and the prefactor A is independent of δφ.^2,30 The loss modulus G″(ω) also follows the same scaling as eqn (4), except that the frequency window exhibiting the ω^1/2 scaling extends to lower frequencies than G′(ω).³²

As φ approaches φ_J from above, the linear-response regime narrows and softening sets in. In the quasi-static limit (ω → 0), we have found that the nonlinear static modulus G_qs(γ) follows a simple scaling law in the pre-yielding regime γ < γ_y,⁴³

where for x ≪ 1 and ∝x^−1/2 for x ≫ 1. The onset strain γ_s ∝ δφ is independent of the exponent α in eqn (1). Thus, one finds

G_qs(γ) ∼ δφγ^−1/2, (γ > γ_s) (6)

for the harmonic potential.⁴³ Note that, in both cases of (ω ≠ 0,γ → 0) and (γ ≠ 0,ω → 0), the moduli satisfy G*(ω) = G_qs(γ) = 0 below φ_J.

In this section, we investigate the case in which both ω and γ are finite and examine how the modulus G*(ω,γ) depends on these parameters. We systematically explore a broad range of ω values within the ω^1/2-scaling regime and consider systems both above and below φ_J.

The storage and loss moduli are evaluated as³⁹

where σ(t) is the stress in the stationary state. The loss modulus G″(ω,γ) is obtained by replacing sin(ωt) in eqn (7) with cos(ωt). Starting from a jammed packing generated using the protocol described in the previous section, we subject the system to multiple oscillatory shear cycles until the particle trajectories become periodic, ensuring that the system has reached the stationary state.

We first consider the case of φ > φ_J. We begin by confirming that G′(ω,γ → 0) = G′(ω) in the linear-response regime follows the ω^1/2 scaling given by eqn (4)³¹ and that G′(ω → 0,γ) = G_qs(γ) in the quasi-static limit follows the γ^1/2 scaling in eqn (6)⁴³ (see Fig. 8 in Appendix B). Next, we examine the dependence of the modulus on γ as ω increases from zero to a finite value larger than ω_s. As shown in eqn (4), the prefactor A is independent of δφ for the harmonic potential. The question is how A is modified with increasing γ in the softening regime. Fig. 1(a) shows ΔG′(ω,γ) as a function of γ for various δφ and ω values. The results show that ΔG′(ω,γ) decreases with increasing γ at fixed ω, a clear signature of softening, and increases with increasing ω at fixed γ. The softening saturates around γ ≈ 10⁻⁴. A similar trend is observed for G_qs(γ) in the quasi-static limit (see Fig. 8(b)). Kawasaki and Miyazaki⁴³ showed that this behavior arises from strain hardening when the initial jammed configurations are prepared by mechanical training. In Fig. 1(b), we plot the scaled quantity ΔG′(ω,γ) = G′(ω,γ) − G_qs(γ) normalized by ω^1/2 as a function of γ/δφ. The collapse of the data indicates that the prefactor A in eqn (4) remains independent of ω but becomes a nonlinear function of γ. The shear-softening exponent is found to be approximately −1/2. A qualitatively similar result for G″(ω,γ) is shown in Fig. 9 in Appendix B, although the softening exponent is slightly smaller. We attribute this to a small deviation from the ω^1/2 scaling of G″(ω) at high frequencies already seen in the linear-response regime as shown in Fig. 8(a).³¹

Fig. 1 (a) γ-Dependence of ΔG′(γ,ω) for several sets of δφ and ω. (b) Scaled ΔG′(ω)/ω^1/2 plotted as a function of γ/δφ. The error bars in all panels represent the standard error of the mean.

Fig. 1(b) demonstrates that ΔG′(ω,γ) follows both the ω^1/2 scaling of eqn (4) and the γ^1/2 scaling of eqn (5) simultaneously, implying the following scaling form:

We find that the scaling function is similar to, but not identical to, defined in eqn (5). In other words, ΔG′(ω,γ) is not simply a superposition of eqn (4) and (5). In Fig. 2, we plot the two moduli normalized by their values at γ = 0: (red filled circles) and (gray filled triangles). The packing fraction is δφ = 10⁻⁴. The figure shows that the two moduli follow similar scaling forms and deviate from them at a common strain value around γ ≳ 10⁻⁴. However, the onset strains γ_s differ. For G_qs(γ), we find γ_s ≈ 1.7 × 10⁻⁵, which is substantially smaller than γ_s ≈ 6 × 10⁻⁵ for ΔG′(ω,γ). This can be understood as follows. The fact that γ_s is proportional to δφ for both quantities (as shown in Fig. 1(b) and 8(b)) suggests a common physical mechanism underlying the softening. In the quasi-static limit, the softening occurs when the particle's non-affine displacement under strain becomes sufficiently large to release two particles' overlaps.⁴³ The same mechanism should also operate at finite ω. However, under oscillatory shear with finite ω, a substantial fraction of particle motion is affine-like. Configurations generated predominantly by affine displacements are less efficiently packed compared to those arising from non-affine displacements, which dominates in the quasi-static deformation. Therefore, the breaking of particle's bonds and release of the particle's overlaps occur at a smaller strain in the quasi-static case than in oscillatory shear. In the next section, we discuss the affine contribution to the softening in detail (see also other data in Fig. 2).

Fig. 2 γ-Dependence of the moduli normalized by their values at γ = 0 for δφ = 10⁻⁴. (red filled circles) at ω = 10⁻² taken from Fig. 1. G_qs(γ) (gray filled triangles) from Fig. 8(b). (open circles) evaluated by the Fourier transform of G(t) shown in Fig. 5(a) (see the next section). G(t = 0) (gray solid line) is taken from Fig. 5(a). The dotted vertical line indicates γ_s for G_qs(γ).

If the scaling relation in eqn (8) remains valid down to φ_J, ΔG′(ω,γ) should behave as

ΔG′(ω,γ) ∼ (ωδφ/γ)^1/2, (φ → φ_J+) (9)

for the harmonic potential. In other words, the modulus is expected to vanish at φ = φ_J. Indeed, this is the case when ω ≠ 0 with γ → 0 or γ ≠ 0 with ω → 0. However, as will be shown below, this no longer holds when both ω and γ are finite. The scaling relation in eqn (8) breaks down in the vicinity of φ_J and the modulus remains finite even below φ_J. Furthermore, the stress is no longer linear in strain at φ ≤ φ_J. This can be understood by noting that, for athermal, non-Brownian, and overdamped particles, no collisions occur during an oscillatory shear cycle if γ is sufficiently small. Only when γ exceeds a threshold value do particles begin to collide with their neighbors, producing a finite stress.

Fig. 3 shows ΔG′(ω,γ) as a function φ for several values of γ at a fixed frequency ω = 10⁻². For the φ -axis, we use 〈φ〉 ≡ 〈φ_J〉 + δφ, where 〈φ_J〉 is the sample-averaged jamming transition point, since δφ is the control parameter throughout this study. The range of 〈φ〉 values shown in the figure corresponds to −6 × 10⁻³ ≤ δφ ≤ 10⁻². The figure demonstrates that the scaling relation in eqn (4) holds in the limit γ → 0: at γ = 10⁻⁴, ΔG′(ω,γ) is independent of φ above φ_J and drops discontinuously to zero below φ_J. As γ increases, this discontinuous drop becomes increasingly rounded. As predicted in eqn (9), ΔG′(ω,γ) ∝ δφ^1/2 in the softening regime (γ > γ_s) at fixed ω, as indicated by the black dashed line in Fig. 3 at γ = 10⁻². However, as φ approaches φ_J for fixed γ, ΔG′(ω,γ) deviates from eqn (9), crosses smoothly through φ_J, and vanishes at a packing fraction below φ_J. The packing fraction φ₀, at which ΔG′(ω,γ) reaches zero, decreases with increasing γ. We find that the modulus obeys a simple algebraic form:

ΔG′(ω,γ) ∼ (φ − φ₀)^β (10)

Fig. 3 ΔG′(ω,γ) as a function of 〈φ〉 and γ at a fixed frequency of ω = 10⁻². The thin black line on the ΔG′(ω,γ) = 0-plane represents the jamming transition point φ_J ≈ 0.8455, while the thin red line indicates the onset strain of softening, γ_s. The symbols denote the modulus at γ = 10⁻⁴ (navy blue triangles), 10⁻³ (light blue squares), 10^−2.5 (green pentagons), and 10⁻² (orange hexagons). The dashed line at γ = 10⁻² represents a fit to δφ^1/2. The open circles indicate the packing fraction φ₀ at which ΔG′(ω,γ) vanishes. The thin gray lines on the ΔG′(ω,γ) = 0 plane connect φ_J and the data for the point-to-loop reversible transition line closest to φ_J reported in the study by Schreck et al.⁴⁷ (one between φ_J and φ₀) and Nagasawa et al.⁴⁸ (one below φ₀).

Our preliminary simulations yield a positive exponent β ≈ 1.2 for G′(ω,γ) at γ = 10⁻². β appears to vary with γ and ω, while it is independent of φ₀. For the loss modulus G″(ω,γ), the exponent is larger (Fig. 10 in Appendix B). A detailed analysis is needed, but it is challenging because the relaxation toward stationary states becomes prohibitively slow and is expected to depend sensitively on the system size near φ₀. We leave this important problem for future work.

The emergence of stress at φ ≥ φ₀ is caused by transient contact networks formed by colliding particles under oscillatory shear. The emergence of such networks with increasing γ below φ_J has been studied in the context of the reversible–irreversible (or absorbing) transition of particle trajectories.^46–50 These studies show that a system subjected to oscillatory shear exhibits distinct dynamical phases. For small γ, particles move affinely and return to their initial positions after each cycle. This small γ-regime is called the point-reversible phase, as the particles remain immobile in the co-moving frame. As γ increases, particles begin to interact with their neighbors and follow non-affine trajectories, yet still return to their original positions after each cycle. This regime is referred to as the loop-reversible phase. At even larger γ, the particles diffuse irreversibly and do not return to their original configurations, defining the irreversible phase.^46–50 Since the transition point from the point-reversible to loop-reversible phases, φ_PL, is the density at which the stress-bearing contact networks percolate, it is expected that φ_PL is the point at which the elastic response becomes finite, i.e., φ_PL = φ₀. The two thin gray lines on the ΔG′(ω,γ) = 0-plane in Fig. 3 connect φ_J with the available data for the point-to-loop reversible transition reported in the studies by Schreck et al.⁴⁷ and Nagasawa et al.⁴⁸ All values of φ₀ (open circles) lie between these two lines. The comparison remains semi-quantitative due to the limited resolution of earlier studies near φ_J.^46–50 A more precise determination of φ_PL (= φ₀), together with a detailed characterization of the critical behavior and its frequency dependence, is needed.

4. Transient stress relaxation

In this section, we investigate the transient stress relaxation of jammed packings, which represents the real-time analogue of the oscillatory-shear measurement and thus complements the results discussed in the previous section. A step shear strain γ is applied to the system at time t = 0, and the resulting shear stress σ(t) and other observables are monitored as functions of time t (>0). For sufficiently small γ, where the linear-response approximation holds, the constitutive equation for the shear stress can be written aswhere [small gamma, Greek, dot above](t) is the shear rate and G(t) is the relaxation modulus. The Fourier transform of G(t) is related to the complex modulus G*(ω) by:⁴⁴where ΔG(t) = G(t) − G_qs. For a step strain, where [small gamma, Greek, dot above](t) = γδ(t), eqn (11) reduces to σ(t) = G(t)γ. These equations imply that the ω^1/2 scaling, G*(ω) ∼ Aω^1/2 in eqn (4), corresponds to an algebraic decay of the stress:

G(t) = σ(t)/γ ∼ t^−1/2 (13)

The stress relaxation is expected to persist up to t_s = ω_s⁻¹, the inverse of the onset frequency of the ω^1/2-scaling regime. Such scale-free relaxation has been observed near the jamming transition φ_J.^36,38

The main question addressed in this section is how the algebraic decay in eqn (13) changes with increasing strain amplitude γ. Recall that the algebraic behaviors of G*(ω) ∼ ω^1/2 and G(t) ∼ t^−1/2 observed in the linear-response regime can be explained, within the harmonic approximation, as the superposition of the exponential relaxation modes of a jammed packing within a single harmonic basin of the potential energy landscape.³⁰ With increasing γ beyond the onset strain γ_s, nonlinear response or softening sets in. In this regime, the harmonic approximation is expected to break down, and the relaxation of the initial configuration toward the terminal state likely involves traversing saddles in the rugged energy landscape rather than a superposition of relaxations in a single basin. This naturally leads to deviations from the G(t) ∼ t^−1/2 behavior. If γ increases further, the configuration becomes so strongly distorted from the initial jammed packing that the resulting relaxation dynamics are expected to resemble those of completely random configurations. This is analogous to the steepest-descent dynamics of particle systems quenched from high temperature T = ∞ to T = 0, which have been recently studied through simulations.^26,51,53,58 These studies report that observables such as the energy E(t) decay algebraically with a nontrivial exponent β ranging between 0.7 and 1.3 depending on the spatial dimension and on the temperature at which initial configurations are prepared.^51,53

Below, we evaluate G(t) and E(t), together with several other observables, over a broad range of γ values and show that G(t) and E(t) decay as t^−1/2 at γ < γ_s and E(t) decays as t^−β at γ ≫ γ_y (where γ_y is the yielding threshold) with the same exponent β as reported for the steepest-descent dynamics from T = ∞.^51,53 Surprisingly, we find that G(t) ∼ t^−1/2 holds even in the softening regime, at γ_s < γ < γ_y. We show that this apparently robust exponent 1/2 in the softening region cannot be fully accounted for by the harmonic approximation. The dependence of the transient stress relaxation on γ has been investigated in pioneering work by Boschan et al.^37,40 Here, we extend their analysis and explore a much broader range of γ values, thereby providing continuous connection between the linear-response regime and the large-γ limit.

In our simulation, all particle positions are instantaneously sheared at t = 0, such that (x_j,y_j) → (x_j + γy_j,y_j) for j = 1, …, N. This sudden deformation introduces a force imbalance in the system. We then monitor the relaxation of the particles toward a new mechanical equilibrium for t > 0, while maintaining the applied strain. The particle dynamics follow the overdamped equation of motion given in eqn (2). We compute several relevant physical observables: the shear stress or the relaxation modulus, G(t) = σ(t)/γ, the potential energy change (with respect to the unsheared system) , where E_γ=0 is the energy before applying the step strain, the average magnitude of the force, , and the non-affine particle displacement, .

We first confirm the transient relaxation in the linear-response regime, consistent with previous studies.^36–38 Fig. 4 shows the observables for three values of δφ at a fixed value of γ = 10⁻⁶, which is sufficiently small to ensure the linear response. Fig. 4(a) and (b) show that both G(t) and E(t) decay as t^−1/2 at long times before eventually saturating at large t.^36–38 The insets present the normalized relaxation function, defined as

where . In Fig. 4(a), f_G(t) demonstrates that the prefactor of the t^−1/2 decay is independent of δφ, in agreement with the ω^1/2 scaling predicted in eqn (4). We also confirm that G(t → ∞) coincides with G_qs. The relaxation of E(t) in Fig. 4(b) is essentially identical to that of G(t), since E(t) ∼ G(t)γ² in the linear-response regime. These algebraic behaviors can be rationalized within the harmonic approximation, and the detailed derivations of these relations are summarized in Appendix C. Fig. 4(c) and (d) show that the average force relaxes as F(t) ∼ t^−3/4, while the non-affine particle displacement increases as δr(t) ∼ t^1/4. These exponents can be understood by noting that the interparticle force is balanced with the Stokes drag, F(t) ∼ ζṙ(t). Substituting this to the energy dissipation rate, dE/dt = −ζṙ²(t), one obtains F(t) ∼ t^−3/4. Similarly, integrating the velocity over time gives δr(t) ∼ t^1/4.

Fig. 4 Transient relaxation in the linear-response regime at a fixed strain γ = 10⁻⁶ for δφ = 10⁻⁵ (blue), 10⁻⁴ (green), and 10⁻³ (red). (a) Relaxation modulus G(t). The inset shows the normalized function f_G(t). (b) Potential energy E(t). The inset shows f_E(t). (c) Average force F(t). (d) Non-affine particle displacement δr(t). The dashed lines indicate the algebraic scaling.

Now, let us systematically investigate how the transient relaxations evolve with increasing γ. Fig. 5 shows the time evolution of the observables for a broad range of γ values from γ = 10⁻⁶, well within the linear-response regime, up to γ = 1, which extends far beyond the yielding threshold γ_y ≈ 2 × 10⁻². The packing fraction is fixed at δφ = 10⁻⁴. The quantities F(t) and δr(t) are scaled by γ, and ΔE(t) by γ², such that all observables become independent of γ in the linear-response regime. We confirm that G(t → ∞) in Fig. 5(a) agrees with the quasi-static modulus G_qs(γ) in Fig. 8(b). These results reveal several nontrivial features.

Fig. 5 The transient relaxation for various γ values (see the legend) at δφ = 10⁻⁴. (a) G(t) = σ/γ, (b) E(t)/γ², (c) F(t)/γ, and (d) δr(t)/γ. The insets in (a) and (b) show the time evolution of the normalized functions . The dashed lines represent the algebraic scaling.

First, G(t) and E(t)/γ² at short times decrease with increasing γ. In other words, not only the quasi-static values G(t → ∞) = G_qs(γ) and E(t → ∞)/γ² but also the instantaneous responses G(t = 0) and E(t = 0)/γ² exhibit shear softening. The short-time response to shear is dominated by the affine displacements of the constituent particles. Our results indicate that the deformation is affine, but that the response becomes progressively nonlinear. The γ-dependence of G(t = 0) is plotted as the solid gray line in Fig. 2. It shows that the softening of G(t = 0) is weaker than that of G_qs(γ) and ΔG′(ω,γ). Interestingly, the onset strain for softening γ_s for G(t = 0) is comparable to that of ΔG′(ω,γ). We also confirm that γ_s for G(t = 0) scales as δφ (results not shown). This observation further supports the argument in the previous section that the affine contribution to ΔG′(ω,γ) shifts the onset strain to a larger value than that of G_qs(γ).

Second, and more interestingly, all observables exhibit an algebraic time dependence over the entire range of γ values, except for G(t) in the large-γ regime beyond γ_y. We focus, in particular, on E(t) ∼t^−β_E and F(t) ∼t^−β_F and extract the positive exponents β_E and β_F by fitting the data shown in the inset of Fig. 5(b) and of Fig. 5(c), respectively. Fig. 6 shows the γ-dependence of β_E (red circles) and β_F (blue triangles). The onset strain for softening γ_s, obtained from G_qs(γ), and the yielding threshold γ_y are indicated by arrows.

Fig. 6 γ-Dependence of the exponents β_E and β_F for the algebraic relaxation of E(t) (red circles) and F(t) (blue triangles), with error bars representing the standard error, respectively. The horizontal solid lines represent the linear-response values, β_E = 1/2 and β_F = 3/4. The horizontal dotted lines denote the values reported for the steepest-descent dynamics in systems quenched from T = ∞ to T = 0:^51,53 β_F = 0.92 and β_E = 2β_F − 1 = 0.84. The two arrows at the top of the panel indicate the onset strain γ_s for G_qs(γ) and the yielding threshold γ_y.

As expected, we observe β_E ≈ 1/2 and β_F ≈ 3/4, in the small-γ regime where the linear-response approximation remains valid. On the other hand, as γ exceeds γ_y, both β_E and β_F increase sharply and then saturate at constant values around γ = 1. The saturated values at the largest γ are very close to β_F = 0.92 and β_E = 2β_F − 1 = 0.84 reported previously^51,53 for steepest-descent dynamics in two-dimensional systems obtained by quenching particle configurations prepared at T = ∞ down to T = 0. These observations confirm that configurations prepared under large shear strains are effectively as random as those prepared at T = ∞ and thus share the same relaxation dynamics. We further find that G(t) at large strain, e.g. γ = 1, decays rapidly and does not exhibit the algebraic relaxation observed in E(t). This can be explained by noting that, in contrast to the scalar nature of E(t), a tensorial observable such as the stress dephases quickly and vanishes when the initial configuration is random.

A particularly interesting aspect is the softening regime γ_s < γ < γ_y. In this range, the exponents β_E = 1/2 and β_F = 3/4 remain unchanged even though the linear-response approximation is no longer expected to hold in this regime. As shown in Fig. 5(a) and (d), the exponents for the relaxation modulus, β_G = 1/2, and for the non-affine displacement, β_δr = 1/4, also remain identical to their linear-response values, before they begin to deviate significantly at γ > γ_y. In particular, the collapse of f_G(t) for γ < γ_y shown in the inset of Fig. 5(a) implies a decoupling of the t-dependence and γ-dependence of G(t), such that ΔG(t) = G(t) −G(∞) can be written as

ΔG(t) = {G(0) −G(∞)}f_G(t) (15)

with f_G(t) being independent of γ. The γ dependence responsible for the softening is thus entirely encoded in the prefactor G(0) − G(∞). This superposition of t^−1/2 scaling and γ^1/2 scaling shown in eqn (15) stands in sharp contrast to the frequency-domain representation, ΔG′(ω,γ), given in eqn (8), for which the superposition picture breaks down at large γ near γ_y, as demonstrated in Fig. 1(b). To demonstrate how ΔG(t) in the time domain differs from that in the frequency domain, ΔG(ω,γ), we evaluate the Fourier transformation of eqn (15). The data for f_G(t) are taken from the inset of Fig. 5(a), G(0) from Fig. 2 and G(∞) from Fig. 8(b). The resulting modulus, denoted as (red open circles in Fig. 2), shows that the γ-dependence resembles that of the instantaneous modulus G(0). However, it does not reproduce obtained directly from oscillatory-shear simulations. The latter shows far more pronounced softening. This discrepancy demonstrates that the transient relaxation dynamics is not equivalent to the stationary dynamics under oscillatory shear in the nonlinear softening regime.

The robustness of the exponents could, in principle, be attributed to the robustness of the harmonic approximation even beyond the linear-response regime, because the characteristic t^−1/2 relaxation arises if the system lies within a single energy basin (even if that basin is distorted by the applied shear). However, we find that this is not the case. To show this, we compute the relaxation dynamics within the harmonic approximation using the configurations of packings obtained at t → ∞ after the step strain. The details of the method are provided in Appendix C. Fig. 7 shows the resulting ΔG(t) and ΔE(t), together with the corresponding simulation data from Fig. 5(a), for three representative values of γ. For the smallest values γ = 10⁻⁶, where the system is well within the linear-response regime, the simulation and harmonic-approximation data are essentially indistinguishable. For γ = 10⁻⁴, which lies in the softening regime, the harmonic approximation still reproduces the t^−1/2 algebraic decay of ΔG(t), but the terminal relaxation time, where the deviation from the power law begins, is substantially shorter than that in the simulations. For ΔE(t), as shown in Fig. 7(b), the discrepancy between the simulation and the harmonic approximation is even more pronounced: the approximation not only underestimates the terminal relaxation time but also significantly overestimates the amplitude of ΔE(t). At the largest strain, γ = 10⁻¹, which is above the yielding threshold, the harmonic approximation fails completely to reproduce the algebraic relaxation, whereas the simulation data still exhibit a power-law decay, albeit with a slightly larger exponent.

Fig. 7 (a) ΔG(t) for γ = 10⁻⁶ (blue), γ = 10⁻⁴ (green), and γ = 10⁻¹ (red) for δφ = 10⁻⁴. Dots represent the simulation data in Fig. 5(a) and the solid lines are the results obtained from the harmonic approximation (see the text). The dotted line represents t^−1/2. (b) The same data for ΔE(t)/γ².

These findings indicate that the algebraic relaxations ΔG(t) and ΔE(t) ∼t^−1/2 are robust, but they cannot be fully explained by the harmonic approximation. These behaviors call for a more careful analysis of how the dynamics depend on the initial degree of distortion.

5. Summary

In this study, we investigated the pre-yielding mechanical response of athermal jammed packings near the jamming transition point φ_J. Our primary goal was to provide a unified framework for the two algebraic scalings observed in the stress response to shear: one is the ω^1/2 scaling (and corresponding t^−1/2 scaling) in the linear-response limit (γ → 0) and the other is the γ^1/2 scaling in the quasi-static limit (ω → 0). Motivated by the work of the Tighe group^40,41 on the interplay between frequency and strain amplitude, we examined the two complementary rheological protocols: oscillatory shear and transient relaxation. The two protocols provide identical information in the linear-response regime, as they are related by a simple Fourier transformation. However, in the nonlinear regime, they are no longer equivalent and are expected to provide complementary information. The oscillatory shear protocol provides detailed information on scaling behavior in stationary states at long times, whereas the roles of higher harmonics, which are assumed to be small in the present study, remain elusive and difficult to interpret. In contrast, the transient relaxation protocol is advantageous for directly probing real-time dynamics, where all nonlinear effects are fully taken into account.

First, using oscillatory shear, we showed that the complex modulus, G*(ω,γ), exhibits the ω^1/2 scaling and γ^1/2 scaling simultaneously, as shown in eqn (8). However, this behavior does not correspond to an exact superposition of the two scalings, as evidenced by the quantitative differences in their onset strains. This finding suggests a nonlinear coupling between affine and non-affine dynamics when the frequency is finite. Moreover, we found that G*(ω,γ) remains finite even below φ_J and vanishes at a packing fraction slightly lower than φ_J. This critical density is close to the point-to-loop reversible transition, φ_PL, which is related to the reversible–irreversible (or absorbing) transition.^46–50 This correspondence implies that the underlying geometrical transition associated with the reversible–irreversible transition is closely linked to the emergence of nonlinear viscoelastic behavior near jamming.

Second, we measured the transient relaxation following a step strain. We observed power-law relaxation, such as t^−1/2, as a real-time manifestation of the ω^1/2 scaling for G(t) and other observables. This behavior is expected in the linear-response regime, where γ is small, and can be rationalized within the harmonic approximation. Remarkably, the t^−1/2 algebraic decay persists all the way up to the yielding threshold γ_y, beyond which qualitatively distinct algebraic behaviors emerge. Importantly, this persistence cannot be explained simply by the robustness of the harmonic approximation.

Our work serves as a starting point for understanding the nonlinear dynamics of jammed packings. However, we are left with more questions than answers. An important question concerns the connection between the mechanical response above and below the jamming transition point φ_J. From a geometrical standpoint, the reversible–irreversible transition predicts that the region where the stress remains finite corresponds to the loop-reversible phase, in which the particles return to their original positions after each oscillatory shear cycle, observed both above and below φ_J. However, the mechanical properties below φ_J differ markedly from those above. It would be particularly interesting to probe the distinction between the loop-reversible phases below and above φ_J (which is difficult to capture solely from the properties of particle trajectories) by analyzing the relaxation dynamics of configurations across φ_J.

The robustness of the t^−1/2 scaling at large γ is also interesting. The study of the relaxation under shear distortion shares a compelling analogy with steepest descent dynamics. The relaxation behavior in the steepest descent dynamics sensitively depends on the temperature T at which the initial configurations are prepared. If T = ∞, the relaxation is algebraic with the same exponents as we reported for the strongly distorted configurations. In the low-T limit, the single basin dynamics dominates. In the intermediate T limit, the behavior of the relaxation dynamics is neither of these limiting cases. We may encounter similarly hierarchical strain-dependent regimes, each characterized by distinct relaxation dynamics.

From this perspective, the nonlinear rheology subject to frequency or time dependent perturbation near φ_J provides a rich and fertile setting that connects several different disciplines: the reversible–irreversible transition, the steepest-descent dynamics, and nonlinear rheology. Pursuing this direction would be a worthwhile avenue for future study.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this article are available at ZENODO at https://zenodo.org/records/17746376.

Appendices

Appendix A. FIRE algorithm and shear-stabilization protocol

In this section, we describe the details of the FIRE algorithm and the shear-stabilization protocol used in this study.

The FIRE algorithm is employed for all energy-minimization procedures, including the preparation of initial configurations via compression–decompression cycles, as well as quasi-static simple shear. Within this algorithm, we perform molecular dynamics simulations using Newton's equations of motion,

We integrate this equation using a semi-implicit Euler scheme. The time-step width is chosen in the range Δt ∈ [10⁻⁶,0.05]. At each step, the total power is given by

and the particle velocities are updated as

ṙ_j = (1 − α)ṙ_j + αF̂_j|ṙ_j| (A3)

where F̂_j is the unit vector in the direction of − ∂U/∂r_j. If p > 0 for five consecutive steps, Δt and α are updated to min(1.1Δt,0.05) and 0.99α, respectively. If p ≤ 0, Δt and α are reset to 0.5Δt and 0.1, respectively. These parameters are identical to those used in ref. 43. The procedure is then repeated by returning to eqn (A1) and recalculating eqn (A2) and (A3).

We incorporate a shear-stabilization protocol into the FIRE algorithm during the compression–decompression process. In addition to the particle equations of motion, we solve the equation of motion for γ,

where W = (100L/σ_S)⁻¹, as proposed in ref. 43. Analogous to eqn (A3), the shear rate is updated aswhere σ/|σ| denotes the sign of the shear stress. We use the same definition of the power p as given in eqn (A2), while neglecting the contribution from σ[small gamma, Greek, dot above].

Appendix B. Supplemental data for oscillatory shear

We first reproduce the ω-dependence of the complex modulus, G*(ω) = G′(ω) + iG″(ω), in the linear-response regime. Fig. 8(a) shows the ω-dependence of the storage modulus, G′(ω), and the loss modulus, G″(ω), for several values of δφ at γ = 10⁻⁷, well below the onset of the softening regime γ_s. The qualitative features of the ω-dependence are identical to those reported in previous studies:^30–32 for ω > ω_s, both G′(ω) and G″(ω) exhibit ω^1/2 power-law behavior. Note that the ω^1/2-scaling regime for G″(ω) extends to smaller ω than ω_s, as reported by Hara et al.³² For ω ≪ ω_s, G′(ω) converges to its long-time limit G(t → ∞), whereas G″(ω) is proportional to ω. At high frequencies, ω ≳ 1, G′(ω) saturates to the short-time limit of the relaxation modulus, G(t = 0) (see Section 4), whereas G″(ω) starts to decrease as reported in ref. 31, which originates from the microscopic dynamics. In this study, we used Stokes drag in the equation of motion for the particles (eqn (2)), which is suitable for describing non-Brownian particles, instead of the Durian drag model used in other studies.^30,31,41 A noticeable feature of Fig. 8(a) is that the window over which the ω^1/2 scaling for G″(ω) holds is relatively narrow compared with that of G′(ω), owing to the crossover to the high-frequency regime. This limited scaling range makes the identification of the ω^1/2 scaling in our study more challenging, as discussed in Section 3.

Fig. 8 (a) ω dependence of G′(ω) (the filled symbols) and G″(ω) (the open symbols) for various δφ values at γ = 10⁻⁷. The dashed line represents ω^1/2. (inset) G*(ω)/δφ^1/2 vs. ω/δφ. The black solid and dashed lines represent G″(ω) ∝ ω and G′(ω), G″(ω) ∝ ω^1/2, respectively. (b) G_qs(γ)/δφ^1/2 as a function of γ/δφ for various δφ values. The dashed line indicates G_qs(γ) ∝ γ^−1/2. (inset) The raw data of G_qs(γ).

In Fig. 8(b), we show the γ-dependence of the modulus in the quasi-static limit, G_qs = G′(ω → 0, γ), which reproduces the results of the previous studies.^39,40,42,43 We observe that G_qs(γ) exhibits softening as G_qs(γ) ∼ γ^−1/2 at γ > γ_s and then saturates and eventually increases as the system enters the hardening regime at large γ.⁴³

Fig. 9 shows the γ-dependence of the loss modulus G″(ω,γ), obtained by using the same set of parameters as those for ΔG′(ω,γ) shown in Fig. 1. The γ-dependence of G″(ω,γ) is qualitatively similar to that of ΔG′(ω,γ), but slight quantitative differences are observed. The shear-softening exponent is smaller than 1/2, as shown in Fig. 9(b). This discrepancy originates from the narrow ω^1/2-scaling window already observed in the linear-response regime (Fig. 8(a)).

Fig. 9 (a) γ-Dependence of G″(ω,γ) for several sets of δφ and ω. (b) Scaled G″(ω,γ)/ω^1/2 plotted as a function of γ/δφ. The error bars in all panels represent the standard error of the mean.

Finally, we examine how the complex modulus G*(ω,γ) vanishes as φ is decreased. As shown in Fig. 3, G′(ω,γ) decreases to zero below φ_J, and the density at which it vanishes, φ₀, decreases with increasing γ. We plot G′(ω,γ) and G″(ω,γ) as a function of φ − φ₀ and fit them with a power law, G′(ω,γ), G″(ω,γ) ∝ (φ − φ₀)^β, where φ₀ and β are fitting parameters. Fig. 10 shows the results of the fitting at ω = 10⁻² and γ = 10⁻². Note that, similar to φ_J, φ₀ varies from sample to sample; therefore, we performed the fitting using the individual φ₀ for each configuration. We find that the exponent β for G′(ω,γ) (β ≈ 1.2) and for G″(ω,γ) (β ≈ 2) differs.

Fig. 10 G′(ω,γ) (filled symbols) and G″(ω,γ) (open symbols) as a function of φ − φ₀ at ω = 10⁻² and γ = 10⁻². The data in Fig. 3 are used. Different symbols correspond to different samples. Dashed lines represent (φ − φ₀)^β with β = 1.2 for G′(ω,γ) and β = 2 for G″(ω,γ).

Appendix C. Supplemental data for transient stress relaxation

We briefly review the derivation of the power-law relaxation for the shear stress and energy, σ(t), E(t) ∼t^−1/2, using the harmonic approximation based on the vibrational mode analysis.⁵³ We consider the dynamics of the shear stress in terms of the displacement vector u(t) and the Hessian matrix of the inherent structure of a packing. The potential energy can be expressed aswhere is the Hessian matrix and u(t) = r(t) − r^eq is the dN-dimensional displacement vector of the N particles from their mechanical equilibrium position r^eq = r(t → ∞), with d being the spatial dimension. The shear stress can be expanded in u(t) as

σ(r^eq + u(t)) ≈ σ(r^eq) + Ξ_xy(r^eq) × u(t) (C2)

where Ξ_xy is the shear stress gradient also known as the affine force field defined by

For the virial shear stress given in eqn (3), the gradient term with respect to particle j is

where κ_jk = d²U/drjk2, t_jk = dU/dr_jk, and n_jk = r_jk/r_jk. Using the eigenvalues {λ_k} and corresponding eigenvectors {e_k} of the Hessian matrix, the displacement vector can be expanded as , with the amplitude given by c_k(t) = u(t) × e_k due to the orthogonality of e_k. Under the harmonic approximation, the dynamics follow the linearized equation of motion which gives the time evolution of c_k(t) as c_k(t) = c_k(0)e^−λ_kt/ζ. Using the relationship, c_k(t) = u(t) × e_k, mentioned above, the relaxation of the shear stress, Δσ(t) = σ(r^eq + u(t)) − σ(r^eq), can be expressed aswhere is the density of eigenvalues and A(λ) = (Ξ_xy(r^eq) × e)(u(0) × e) represents the weight of each mode for the stress relaxation. The eigenvalue density is related to the vibrational density of states D(Ω) by ρ(λ) = D(Ω)(dΩ/dλ), where .

Similarly, the energy relaxation can be derived by expanding the energy using u(t) and (see eqn (C1)). With the linearized equation of motion described above, the energy relaxation follows

where B(λ) = (u(0)·e)² is the weight of each mode for the energy relaxation.⁵³

In the vicinity of φ_J, D(Ω) develops a plateau at low frequencies, so that ρ(λ) ∼ λ^−1/2. As shown below, A(λ) in eqn (C5) is approximately constant, whereas B ∼ λ⁻¹ in eqn (C6). Therefore, assuming that these integrals are dominated by the low-λ contribution, both Δσ(t) and ΔE(t) scale as t^−1/2, which reproduces the results of Fig. 4 in the linear-response regime.

Next, we perform the harmonic approximation analysis for sheared systems. We directly compute the eigenvalues of the Hessian matrix from numerically obtained jammed packings. Fig. 11 shows the density of eigenvalues ρ(λ) and the vibrational density of states D(Ω) (inset) for the packing at γ = 10⁻⁶, 10⁻⁴, and 10⁻¹, at δφ = 10⁻⁴. D(Ω) (and hence ρ(λ)) values for the first two values of γ are almost identical to those of the unstrained system, except for a small deviation at low-Ω. They exhibit a plateau, a feature universally observed near φ_J.^4,59 D(Ω) for γ = 10⁻¹ is also similar except for a larger onset frequency of the plateau, Ω_*, which is due to the decrease of φ_J at large γ.⁴³ Using the computed eigenvalues and eigenvectors shown in Fig. 11, we evaluate A(λ) and B(λ) from eqn (C5) and (C6). Fig. 12(a) and (b) show the λ -dependence of A(λ) and B(λ). For γ = 10⁻⁶ and 10⁻⁴, we find that A(λ) ∼ const. and B(λ) ∼λ⁻¹ over a broad range of λ values, corresponding to the plateau of D(Ω). The small difference in B(λ) for the two values of γ is reflected in the overestimation of ΔE(t) in the softening regime, as shown in Fig. 7. For γ = 10⁻¹, on the other hand, both A(λ) and B(λ) differ substantially from those for other γ.

Fig. 11 The density of eigenvalues ρ(λ) in mechanical equilibrium for γ = 10⁻⁶ (blue), 10⁻⁴ (green), and 10⁻¹ (red) at δφ = 10⁻⁴. The dashed line indicates λ^−1/2. (inset) The vibrational density of states D(Ω).

Fig. 12 λ-dependence of (a) A(λ) = (Ξ_xy(r^eq) × e)(u(0) × e) and (b) B(λ) = (u(0) × e)² scaled by γ for γ = 10⁻⁶ (blue), γ = 10⁻⁴ (green), and γ = 10⁻¹ (red) at δφ = 10⁻⁴. The dashed line in panel (b) indicates λ⁻¹ for B(λ).

Acknowledgements

We thank M. Otsuki, A. Ikeda, H. Mizuno, H. Ikeda, Y. Nishikawa, Y. Hara, K. Yoshii, H. Iwatsuki, and S. Tomioka for their valuable discussions. This work was supported by the JST FOREST Program (grant no. JPMJSP2125), JST SPRING (grant no. JPMJSP2125), and JSPS KAKENHI (grant no. JP22H04472, JP23H04503, JP23KJ1068, and JP24H00192). H. B. would like to take this opportunity to acknowledge support from the “THERS Make New Standards Program for the Next Generation Researchers”.

Notes and references

1. A. J. Liu and S. R. Nagel, Nature, 1998, 396, 21–22.
2. C. S. O'Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 011306.
3. M. van Hecke, J. Phys.: Condens. Matter, 2010, 22, 033101.
4. L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2005, 95, 098301.
5. H. Mizuno, H. Shiba and A. Ikeda, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, E9767–E9774.
6. T. S. Majmudar, M. Sperl, S. Luding and R. P. Behringer, Phys. Rev. Lett., 2007, 98, 058001.
7. G. Katgert and M. van Hecke, Europhys. Lett., 2010, 92, 34002.
8. M. Wyart, Ann. Phys. Fr., 2005, 30, 1.
9. E. DeGiuli, A. Laversanne-Finot, G. Düring, E. Lerner and M. Wyart, Soft Matter, 2014, 10, 5628–5644.
10. G. Parisi, P. Urbani and F. Zamponi, Theory of Simple Glasses: Exact Solutions in Infinite Dimensions, Cambridge University Press, 2020.
11. D. J. Durian, Phys. Rev. Lett., 1995, 75, 4780–4783.
12. P. Olsson and S. Teitel, Phys. Rev. Lett., 2007, 99, 178001.
13. M. Otsuki and H. Hayakawa, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 011308.
14. T. Hatano, J. Phys. Soc. Jpn., 2008, 77, 123002.
15. T. Hatano, Prog. Theor. Phys. Suppl., 2010, 184, 143–152.
16. B. P. Tighe, E. Woldhuis, J. J. C. Remmers, W. van Saarloos and M. van Hecke, Phys. Rev. Lett., 2010, 105, 088303.
17. F. Boyer, E. Guazzelli and O. Pouliquen, Phys. Rev. Lett., 2011, 107, 188301.
18. P. Olsson and S. Teitel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 030302.
19. E. Lerner, G. Düring and M. Wyart, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 4798–4803.
20. M. Dinkgreve, J. Paredes, M. A. J. Michels and D. Bonn, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 012305.
21. T. Kawasaki, D. Coslovich, A. Ikeda and L. Berthier, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 012203.
22. P. Olsson, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 062209.
23. D. Vågberg, P. Olsson and S. Teitel, Phys. Rev. E, 2016, 93, 052902.
24. E. DeGiuli, G. Düring, E. Lerner and M. Wyart, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 062206.
25. D. Bonn, M. M. Denn, L. Berthier, T. Divoux and S. Manneville, Rev. Mod. Phys., 2017, 89, 035005.
26. A. Ikeda, T. Kawasaki, L. Berthier, K. Saitoh and T. Hatano, Phys. Rev. Lett., 2020, 124, 058001.
27. D. Pan, Y. Wang, H. Yoshino, J. Zhang and Y. Jin, Phys. Rep., 2023, 1038, 1–18.
28. L. Berthier, G. Biroli, L. Manning and F. Zamponi, Nat. Rev. Phys., 2025, 7, 313–330.
29. T. Divoux, E. Agoritsas, S. Aime, C. Barentin, J.-L. Barrat, R. Benzi, L. Berthier, D. Bi, G. Biroli, D. Bonn, P. Bourrianne, M. Bouzid, E. Del Gado, H. Delanoë-Ayari, K. Farain, S. Fielding, M. Fuchs, J. van der Gucht, S. Henkes, M. Jalaal, Y. M. Joshi, A. Lemaiaetre, R. L. Leheny, S. Manneville, K. Martens, W. C. K. Poon, M. Popović, I. Procaccia, L. Ramos, J. A. Richards, S. Rogers, S. Rossi, M. Sbragaglia, G. Tarjus, F. Toschi, V. Trappe, J. Vermant, M. Wyart, F. Zamponi and D. Zare, Soft Matter, 2024, 20, 6868–6888.
30. B. P. Tighe, Phys. Rev. Lett., 2011, 107, 158303.
31. K. Baumgarten and B. P. Tighe, Soft Matter, 2017, 13, 8368–8378.
32. Y. Hara, R. Matsuoka, H. Ebata, D. Mizuno and A. Ikeda, Nat. Phys., 2025, 21, 262–268.
33. W. Phillips, Amorphous Solids: Low-temperature Properties, Springer-Verlag, 1981.
34. W. Schirmacher, G. Ruocco and T. Scopigno, Phys. Rev. Lett., 2007, 98, 025501.
35. P. Charbonneau, E. I. Corwin, G. Parisi, A. Poncet and F. Zamponi, Phys. Rev. Lett., 2016, 117, 045503.
36. T. Hatano, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 050301.
37. J. Boschan, S. A. Vasudevan, P. E. Boukany, E. Somfai and B. P. Tighe, Soft Matter, 2017, 13, 6870–6876.
38. K. Saitoh, T. Hatano, A. Ikeda and B. P. Tighe, Phys. Rev. Lett., 2020, 124, 118001.
39. M. Otsuki and H. Hayakawa, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 042202.
40. J. Boschan, D. Vågberg, E. Somfai and B. P. Tighe, Soft Matter, 2016, 12, 5450–5460.
41. S. Dagois-Bohy, E. Somfai, B. P. Tighe and M. van Hecke, Soft Matter, 2017, 13, 9036–9045.
42. M. Otsuki and H. Hayakawa, Phys. Rev. Lett., 2022, 128, 208002.
43. T. Kawasaki and K. Miyazaki, Phys. Rev. Lett., 2024, 132, 268201.
44. R. Larson, The Structure and Rheology of Complex Fluids, OUP USA, 1999.
45. H. M. Wyss, K. Miyazaki, J. Mattsson, Z. Hu, D. R. Reichman and D. A. Weitz, Phys. Rev. Lett., 2007, 98, 238303.
46. L. Milz and M. Schmiedeberg, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 062308.
47. C. F. Schreck, R. S. Hoy, M. D. Shattuck and C. S. O'Hern, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 052205.
48. K. Nagasawa, K. Miyazaki and T. Kawasaki, Soft Matter, 2019, 15, 7557–7566.
49. H. Matsuyama, M. Toyoda, T. Kurahashi, A. Ikeda, T. Kawasaki and K. Miyazaki, Eur. Phys. J. E, 2021, 44, 133.
50. P. Das, H. A. Vinutha and S. Sastry, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 10203–10209.
51. R. N. Chacko, P. Sollich and S. M. Fielding, Phys. Rev. Lett., 2019, 123, 108001.
52. G. Folena, S. Franz and F. Ricci-Tersenghi, Phys. Rev. X, 2020, 10, 031045.
53. Y. Nishikawa, M. Ozawa, A. Ikeda, P. Chaudhuri and L. Berthier, Phys. Rev. X, 2022, 12, 021001.
54. F. Shuang, P. Xiao, R. Shi, F. Ke and Y. Bai, Comput. Mater. Sci., 2019, 156, 135–141.
55. E. Bitzek, P. Koskinen, F. Gähler, M. Moseler and P. Gumbsch, Phys. Rev. Lett., 2006, 97, 170201.
56. S. Dagois-Bohy, B. P. Tighe, J. Simon, S. Henkes and M. van Hecke, Phys. Rev. Lett., 2012, 109, 095703.
57. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford, 1988.
58. Y. Nishikawa, A. Ikeda and L. Berthier, J. Stat. Phys., 2021, 182, 37.
59. L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 021308.

---