An exact expression of three-body system for the complex shear modulus of frictional granular materials

Abstract

We propose a simple model comprising three particles to study the nonlinear mechanical response of jammed frictional granular materials under oscillatory shear. Owing to the introduction of the simple model, we obtain an exact analytical expression of the complex shear modulus for a system including many monodispersed disks, which satisfies a scaling law in the vicinity of the jamming point. These expressions perfectly reproduce the shear modulus of the many-body system with low strain amplitudes and friction coefficients. Even for disordered many-body systems, the model reproduces the results by introducing a single fitting parameter.

---

1 Introduction

The rheological property of densely dispersed grains, e.g., granular materials, colloidal suspensions, and emulsions, plays an important role in physics and engineering. This rheological property mainly depends on the packing fraction ϕ of the grains. The materials behave like fluids for ϕ < ϕ_J with jamming fraction ϕ_J and exhibit a solid-like elastic response above ϕ_J.^1,2 In the linear response regime (i.e., for small strains), the shear modulus is characterized by the density of states^3–5 and satisfies scaling laws.^6–9 However, the linear response region becomes narrower as ϕ approaches ϕ_J,^10,11 and the nonlinear response becomes relevant due to the plastic deformation associated with the yielding.^12–20

If we are interested in a nonlinear response to an applied oscillatory shear strain, it exhibits a complicated stress–strain curve. Although the storage and loss moduli G′ and G′′ were originally introduced to characterize the linear viscoelasticity of materials, they can be used to characterize nonlinear viscoelasticity or visco-elastoplastic responses to applied strains.²¹ In this case, G′ and G′′ are no longer constants but strongly depend on the strain amplitude γ₀. In particular, we have recognized that G′ decreases with γ₀^11,22–25 and G′′ remains non-zero in the low frequency limit^24,25 for densely dispersed grains.

The theoretical analysis of densely dispersed grains is challenging as a typical many-body problem in non-equilibrium systems. To date, a few theoretical approaches have been proposed for systems related to frictionless particles. The scaling laws for the linear elastic response were derived in terms of the vibrational density of states.^7,8 The Fourier analysis of particle trajectories helps to generate semi-analytical expressions for G′ and G′′.²⁵ Unfortunately, these theories cannot apply to frictional particles because of the history-dependent contact force.^9,24

It is helpful to analyze a simple model with small degrees of freedom to understand the behavior of many-body systems, including densely dispersed grains. This approach has been used in statistical mechanics. The mean-field approximation of the Ising model is a typical example in which the system contains only one Ising spin under the influence of a self-consistently determined mean field.²⁶ For atomic liquids, a cell model, in which a single atom exists in a cage, was used to derive the equation of state.^27,28 The coherent potential approximation for disordered solids has been used to understand electronic band structures.²⁹ The effective medium theory reveals the elastic response of random spring networks.³⁰ In addition, a simple model consisting of two particles was proposed to reproduce the liquid–solid phase transition.³¹ The advantage of such few-body models is that we can obtain exact solutions. The qualitative behavior of the corresponding many-body systems can be determined based on the solutions of the few-body models. Thus, we adopt this approach to determine the nonlinear responses of the frictional dispersed grains.

This study proposes a model consisting of three identical particles to describe the mechanical response of jammed frictional granular materials under oscillatory shear. In Section 2, we introduce the three-body system (TBS). This model can be analytically solved for low-strain amplitudes and friction coefficients near the jamming point in Section 3. In Section 4, we demonstrate that the analytical solution reproduces the storage and loss moduli of many-body systems (MBSs) without any fitting parameter if there is no disorder in the particle configuration. Even if disorder exists, a scaling law for the complex shear modulus for the TBS semi-quantitatively agrees with the numerical simulations of the MBS by introducing a fitting parameter. We discuss and conclude our results in Section 5. In Appendix A, we show the details of the MBS when the particles are initially placed on a triangular lattice. The effect of particle rotation is described in Appendix B. In Appendix C, we derive the analytical expressions for the shear stress and pressure in the TBS. In Appendix D, we relate the complex shear modulus with the hysteresis loop of the stress–strain curve. The details of the disordered MBS are presented in Appendix E. We present the numerical shear modulus for the TBS in Appendix F.

2 Three-body system

We consider two-dimensional granular materials consisting of many grains under oscillatory shear (Fig. 1). Here, the grains constituting granular materials are modeled as frictional spherical particles. Moreover, we introduce a system of three identical particles to simply describe the MBS (Fig. 2). The MBS can contain polydisperse particles, while we assume that the TBS is a monodisperse system. In the TBS, the position r_i(t) = (x_i(t), y_i(t)) of particle i with diameter d at time t is given bywhere [small script l] is the initial distance between particles. We also introduce ε:= 1 − [small script l]/d as the compressive strain. The compressive strain ε in the TBS corresponds to ϕ − ϕ_J in the MBS, as shown in Appendix A. We apply shear strain as

γ(θ) = γ₀ sin θ (4)

with strain amplitude γ₀, phase θ = ωt, and angular frequency ω. Note that we need at least three particles to realize a stable interlocking state.

Fig. 1 Schematics of the ordered MBS (a) and the disordered MBS (b).

Fig. 2 A schematic of the TBS.

We adopt the interaction force f_ij between particles i and j given by

f_ij = (f⁽ⁿ⁾_ijn_ij + f^(t)_ijt_ij)H(r_ij − d), (5)

where f⁽ⁿ⁾_ij and f^(t)_ij denote the normal and tangential forces between the particles i and j.³² The distance between the particles i and j is r_ij = |r_ij| with r_ij := r_i − r_j = (x_ij, y_ij). Here, H(x) is Heaviside's step function satisfying H(x) = 1 for x > 0 and H(x) = 0 otherwise. The normal and tangential unit vectors are denoted by n_ij := r_ij/r_ij = (n_ij,x, n_ij,y) and t_ij := (−n_ij,y, n_ij,x), respectively. For simplicity, we do not consider the torque balance and, thus, the rotation of the particles. See Appendix B for the effect of the rotation.

The normal force is assumed to be

f⁽ⁿ⁾_ij = −k_nu⁽ⁿ⁾_ij (6)

with the normal elastic constant k_n and normal relative displacement u⁽ⁿ⁾_ij := r_ij − d. Moreover, the tangential force is assumed to be

f^(t)_ij = min(|f̃^(t)_ij|,μf⁽ⁿ⁾_ij)sgn(f̃^(t)_ij), (7)

where f̃^(t)_ij = −k_tu^(t)_ij; k_t denotes the tangential elastic constant, and μ denotes the friction coefficient. Here, min(a,b) selects the smaller value between a and b, sgn(x) = 1 for x ≥ 0, and sgn(x) = −1 for x < 0. The tangential displacement u^(t)_ij satisfies for |f̃^(t)_ij| < μf⁽ⁿ⁾_ij with the tangential velocity whereas u^(t)_ij remains unchanged for |f̃^(t)_ij| ≥ μf⁽ⁿ⁾_ij. We refer to the contact with |f̃^(t)_ij| < μf⁽ⁿ⁾_ij as the stick contact and the contact with |f̃^(t)_ij| ≥ μf⁽ⁿ⁾_ij as the slip contact. The tangential displacement, u^(t)_ij, is initially set to zero.

The (symmetric contact) shear stress is given by

σ(θ; γ₀, μ) = σ⁽ⁿ⁾(θ; γ₀, μ) + σ^(t)(θ; γ₀, μ) (8)

with the normal component of σand tangential component of σHere, A corresponds to the area of the system, and we choose as shown in Appendix A. The pressure is given byIn the right-hand sides of eqn (9)–(11), we have omitted the arguments θ, γ₀, and μ. Similar abbreviations are used below. As we are interested in quasistatic processes, we do not consider the kinetic parts of σ and P and the dependence on ω. After several cycles of oscillatory shear, σ(θ) becomes periodic. The storage and loss moduli are given by³³

3 Theoretical analysis

Assuming γ₀ ≪ ε ≪ 1, we analytically obtain G′ and G′′ for the TBS. The derivation of the analytical results can be found in Appendix C.

First, the normal component of the shear stress is given by

The tangential component of the shear stress is given byfor γ₀ < γ_c(μ) with a critical amplitudewhich characterizes the transition from stick to slip states in the contact between the particles. For γ₀ ≥ γ_c(μ), the tangential component of the shear stress is given bywhere Θ = cos⁻¹(1 − 2γ_c(μ)/γ₀). Regions with and correspond to the stick state, and the other regions correspond to the slip state. Owing to this transition in the contact, the stress–strain curve given by eqn (14)–(17) exhibits a hysteresis loop. Eqn (17) does not exhibit a viscoelastic response but a typical elastoplastic response without viscous effect.

Fig. 3 shows the scaled shear stress σ/γ₀ against the scaled strain γ/γ₀ using eqn (4), (8) and (14)–(17) for various values of γ₀ with k_t/k_n = 1.0 and μ = 0.01. The shape of the scaled stress–strain curve is characterized by a parallelogram as a typical elastoplastic response. As γ₀ increases, the maximum value [small sigma, Greek, tilde]_max = (σ/γ₀)|_γ/γ0=1 decreases from a larger value to a smaller value . As shown in Appendix D, the storage modulus G′ is approximately given by [small sigma, Greek, tilde]_max. Hence, the decrease of [small sigma, Greek, tilde]_max in Fig. 3 indicates the decrease of G′. For γ₀ = 0.00003 and 0.0001, the hysteresis loop exists, but the area of the loop is negligible for γ₀ = 0.00001 and 0.001. The loss modulus G′′ is proportional to the area of the loop, as shown in Appendix D. Hence, the dependence of the area on γ₀ indicates that there is a peak in G′′ as γ₀ increases.

Fig. 3 Scaled shear stress σ/γ₀ against γ/γ₀ using eqn (4), (8) and (14)–(17) for various values of γ₀ with k_t/k_n = 1.0, ε = 0.001, and μ = 0.01.

Substituting eqn (8) and (14)–(17) into eqn (12), we obtain the storage modulus as

As γ₀ increases beyond γ_c(μ), G′ decreases from a higher value to a lower value. The corresponding behavior has been observed in the MBS in previous studies.^9,24

Substituting eqn (8) and (14)–(17) into eqn (13), the loss modulus is given by

The loss modulus G′′ is zero for γ₀ < γ_c(μ), whereas G′′ sharply increases with γ₀ when γ₀ exceeds γ_c(μ) and decreases to 0 after a peak. The behavior of G′′ for the TBS qualitatively reproduces that of the MBS in previous studies.²⁴

We adopt the abbreviation for the pressure at γ = 0 as

P₀ := P(θ = 0; γ₀, μ), (20)

which is also obtained asFrom eqn (16), (18), (19) and (21), we derive scaling laws for a given ε aswhere and denote scaling functions. The maximum values of G′ and G′′ are denoted as and , respectively. In the TBS, they are given aswith T(x) = cos⁻¹(1 − 2x_c/x), S(x) = sin(2T(x))/2, and .

4 Comparison with the MBS

We demonstrate the relevance of the TBS analysis based on the simulation of a two-dimensional MBS consisting of N frictional grains. First, we consider a system corresponding to the TBS, where all the particles are identical and initially placed on the triangular lattice with a unit length [small script l] (Fig. 1(a)). The details are shown in Appendix A. Next, we consider a bidisperse system where the number of particles with diameter d is equal to that of particles with diameter d/1.4, and the particles are randomly placed with packing fraction ϕ (Fig. 1(b)). The mass densities of the particles are identical. The details of the disordered MBS are shown in Appendix E. In both systems, the shear strain given by eqn (4) is applied for N_c cycles using the SLLOD equation under the Lees–Edwards boundary condition.³⁴ In the MBS, we replace the normal force as

f⁽ⁿ⁾_ij→ −(k_nu⁽ⁿ⁾_ij + η_nv⁽ⁿ⁾_ij) (27)

with the normal viscous constant η_n and the normal velocity to include the viscous force depending on the relative velocity. The tangential force is replaced by

f^(t)_ij → min(|f̃^(t)_ij|,μf^(n,el)_ij)sgn(f̃^(t)_ij), (28)

with

f̃^(t)_ij → −(k_tu^(t)_ij + η_tv^(t)_ij), (29)

where f^(n,el)_ij = −k_nu⁽ⁿ⁾_ij denotes the elastic part of the normal force with a tangential viscous constant η_t. We measure G′, G′′, and P₀ in the last cycle using eqn (11)–(13). For the ordered MBS, we use N = 64, k_t/k_n = 1.0, and ε = 0.001, whereas N = 1000, k_t/k_n = 0.2, and ϕ = 0.87 are used for the disordered MBS. In both systems, the other parameters are identical: N_c = 20, and with a mass m of larger particles.

As shown in Fig. 4, we plot G′ for the ordered MBS against γ₀ with k_t/k_n = 1.0 and ε = 0.001 for various values of μ as points. Moreover, we plot the analytical results of the TBS obtained using eqn (18) as thin solid lines. Surprisingly, the results of the TBS agree with those of the MBS for γ₀ < 0.003 without any fitting parameters. As γ₀ increases beyond γ_c(μ) shown by the vertical dashed lines, G′ for μ > 0 decreases and converges to a constant, which is equal to G′ for μ = 0. For larger γ₀, G′ for the MBS decreases again, whereas the theoretical G′ for the TBS is constant. This discrepancy results from the violation of condition γ₀ ≪ ε for the analytical calculation. If we numerically solve the TBS to obtain G′ without the assumption γ₀ ≪ ε, G′ decreases after a plateau again as in the case of MBS, although its value in the TBS for γ₀ → 0.1 slightly deviates from that of the MBS, as shown in Appendix F.

Fig. 4 Storage modulus G′ against γ₀ with k_t/k_n = 1.0 and ε = 0.001 for various values of μ. The points represent the results of the ordered MBS. The thin solid lines represent the analytical result given by eqn (18). The vertical dashed lines represent the critical amplitude γ_c(μ) given by eqn (16) for μ = 10⁻⁴, 10⁻³, 10⁻², 10⁻¹, and 1 from left to right.

As shown in Fig. 5, we plot G′′ for the MBS on the triangular lattice against γ₀ with k_t/k_n = 1.0 and ε = 0.001 for various values of μ as points. Moreover, we plot the analytical results of the TBS obtained using eqn (19) as thin solid lines. The analytical result agrees perfectly with the MBS for γ₀ < 0.003 without any fitting parameters. As γ₀ increases beyond γ_c(μ) shown by the vertical dashed lines, G′′ for μ > 0 increases from 0 and decreases after reaching a peak. The peak position of G′′ against γ₀ increases with μ. Thus, our analytical results fail to capture the behavior of G′′ for μ = 1.

Fig. 5 Loss modulus G′′ against γ₀ with k_t/k_n = 1.0 and ε = 0.001 for various values of μ. The points represent the results of the ordered MBS. The thin solid lines represent the analytical results obtained using eqn (19). The vertical dashed lines represent the critical amplitude γ_c(μ) given by eqn (16) for μ = 10⁻⁴, 10⁻³, 10⁻², 10⁻¹, and 1 from left to right.

Consider the disordered MBS shown in Fig. 1(b). Fig. 6 shows the scaled shear stress σ/γ₀ against the scaled strain γ/γ₀ in the disordered MBS with μ = 0.0001. The maximum value [small sigma, Greek, tilde]_max decreases as γ₀ increases. The area S of the curve is the largest for γ₀ = 0.00003. It is interesting that stress–strain curves are not characterized by parallelograms in this case in contrast to Fig. 3. This means that the disordered configuration of particles creates an effective viscosity, and thus, the response to an applied strain becomes visco-elastoplastic.

Fig. 6 Scaled shear stress σ/γ₀ against γ/γ₀ in the disordered MBS for various values of γ₀ with μ = 0.01 and ϕ = 0.870.

The behaviors of G′ and G′′ in the disordered MBS are similar to those of the TBS as shown in Appendix E. Therefore, it is expected that the scaling laws in eqn (22) and (23) for a given ε in the TBS can be used even in this system with corresponding ϕ. This expectation is verified by Fig. 7, in which we plot the scaled moduli and against the scaled strain k_tγ₀/(μP₀(γ₀, μ)) for various values of μ in the disordered MBS. Moreover, we plot the analytical results for the TBS obtained using eqn (25) and (26) as solid lines, which qualitatively reproduce the MBS results for small scaled strain, while the scaling is apparently violated for large scaled strain. Here, we choose k_t/k_n = 1.5 for the TBS to fit the second plateau to that of the MBS. At present, we do not know the relationship between ϕ and the fitting parameter.

Fig. 7 (a) Scaled storage modulus against the scaled strain k_tγ₀/(μP₀(γ₀, μ)) with ϕ = 0.87 and k_t/k_n = 0.2 for various values of μ in the disordered MBS. The solid line represents the analytical result of the TBS given by eqn (25) with k_t/k_n = 1.5. (b) Scaled loss modulus against the scaled strain k_tγ₀/(μP₀(γ₀, μ)) with ϕ = 0.87 and k_t/k_n = 0.2 for various values of μ in the disordered MBS. The solid line represents the analytical result of the TBS given by eqn (26) with k_t/k_n = 1.5.

5 Conclusions

We demonstrated the relevancy of a model of the TBS to describe the complex modulus of jammed frictional granular materials under oscillatory shear. We obtained the analytical expressions for the γ₀-dependence of G′ and G′′ as shown in eqn (16), (18), and (19), which predict the μ-dependence of the critical amplitude γ_c, the decrease of G′, and the peak of G′′ above γ_c for crystalline solids. The analytical expressions lead to the scaling laws given by eqn (22) and (23). Although we have ignored the non-affine motion for crystalline solids, these analytical results quantitatively agree with those of the ordered MBS. Surprisingly, some characteristic features of disordered solids for low strain (or high pressure) can be captured. These results indicate that the analysis of the toy model gives a basis for understanding the nonlinear rheology of frictional granular materials under small strain.

Although the values of the plateaus in G′ for disordered MBS depended on ϕ − ϕ_J,^6,9 the corresponding values of the TBS are independent of ϕ − ϕ_J, as expressed in eqn (18). In addition, our analytical expressions cannot reproduce the second decrease of G′ and increase of G′′ near γ₀ = 10⁻² in the MBS. The discrepancy should result from the disorder because it leads to the ϕ-dependence of G′.⁷ To include the disorder effect, we regarded k_t/k_n as a fitting parameter. In previous studies on models with small degrees of freedom, e.g., the coherent potential approximation,^26,29,30 the corresponding fitting parameters were self-consistently determined. In future studies, we will discuss the self-consistent determination of the parameter for the TBS.

Some researchers are interested in contributions from higher harmonics characterizing the nonlinear response to oscillatory shear,²¹ but the nonlinear viscoelastic moduli characterizing the higher harmonics are negligibly small for jammed frictionless particles.²⁵ However, the higher harmonics in the frictional granular materials require further careful investigation.

Conflicts of interest

There are no conflicts to declare.

Appendix A: details of Ordered MBS

This section explains the details of the ordered MBS consisting of monodispersed particles initially placed on a triangular lattice. We consider a two-dimensional assembly of N frictional particles in a periodic box with sizes along the x and y directions L_x and L_y, respectively. Here, we initially place N = 2N_xN_y particles of diameter d with integers n_x and n_y at r_i asfor 0 ≤ i < N_xN_y with integers n_x, n_y, and i = n_x + N_xn_y. For N_xN_y ≤ i < 2N_xN_y, r_i is defined aswith i = n_x + N_xn_y + N_xN_y. The initial configuration is illustrated in Fig. 8. We choose L_x = N_x[small script l] and with [small script l] = d(1 − ε).

Fig. 8 Initial configuration of mono-dispersed particles on a triangular lattice. The red rectangle, including interactions represented by the blue lines, corresponds to the TBS.

The position r_i and peculiar momentum p_i of particle i with mass m_i and diameter d_i are driven by the SLLOD equation under the Lees–Edwards boundary condition as³⁴

where [small gamma, Greek, dot above](t) = γ₀ω cos ωt and e_x = (1,0) is the unit vector along the x direction. The interaction force f_i is defined aswith d_ij = (d_i + d_j)/2, n_ij = r_ij/r_ij, t_ij = (−n_ij,y, n_ij,x), and r_ij = r_i − r_j = (x_ij, y_ij). The normal force is given by

f⁽ⁿ⁾_ij = −(k_nu⁽ⁿ⁾_ij + η_nv⁽ⁿ⁾_ij) (35)

with a normal viscous constant η_n and

v⁽ⁿ⁾_ij = (v_i − v_j)·n_ij, (36)

where the velocity of particle i is given by . The following model is adopted for the tangential force:

f^(t)_ij = min(|f̃^(t)_ij|,μf^(n,el)_ij)sgn(f̃^(t)_ij), (37)

where f^(n,el)_ij = −k_nu⁽ⁿ⁾_ij denotes the elastic part of the normal force. Here, f̃^(t)_ij is given by

f̃^(t)_ij = −(k_tu^(t)_ij + η_tv^(t)_ij) (38)

with a tangential viscous constant η_t. The tangential velocity v^(t)_ij is given by

v^(t)_ij = (v_i − v_j)·t_ij. (39)

The tangential displacement u^(t)_ij satisfies for |f̃^(t)_ij| < μf^(n,el)_ij, whereas u^(t)_ij remains unchanged for |f̃^(t)_ij| ≥ μf^(n,el)_ij. The tangential displacement u^(t)_ij is set to zero if i and j are detached.

If all the particles are separated, the packing fraction ϕ for the ordered MBS is defined as

Even if contacts exist between the particles, we use eqn (40) by assuming that the contact length d_ij − r_ij is sufficiently lower than d_ij. Using eqn (40), ϕ is defined asThe jamming point of this system iswith ε = 0. The distance from the jamming point is proportional to ε asfor ε ≪ 1.

The shear stress σ is defined by eqn (8) in the main article with the normal component

and tangential componentThe pressure is defined as

We use N_x = 8, N_y = 4, N_c = 20, k_t = k_n, and where m denotes the mass of a particle with diameter d. This model corresponds to a restitution coefficient e = 0.043. We adopt the leapfrog algorithm considering a time step of Δt = 0.05t₀. We choose as the quasistatic shear deformation because G′ and G′′ are almost independent of ω for .

As shown in Fig. 4 and 5, the behaviors of G′ and G′′ of the TBS agree with that of the MBS. We explain the theoretical background of the TBS. The initial configuration is shown in Fig. 8; it contains the unit cell represented by the red rectangle with length [small script l] and height . It contains interactions between the three particles represented by blue lines. Here, we assume that the particles move affinely as

r_i(t) = r_i(0) + γ(θ(t))y_i(0)e_x. (47)

In this case, the corresponding relative distances between the particles in any unit cell are identical.

In particular, in a unit cell containing particles i = i₁, i₂, and i₃ with i₁ = N_xN_y, i₂ = 0, and i₃ = 1, the positions of the particles are given by

The relative distances between these particles are identical to those of the TBS, given by eqn (1)–(3), which indicates that the TBS provides the interaction forces among the three particles. This system includes 2N_xN_y unit cells with identical interaction forces. Hence, the normal and tangential components of σ are given byThe pressure is also given by:Using the relation corresponding to σ⁽ⁿ⁾, σ^(t), and P coincide with eqn (8)–(11). Hence, if the assumptions of the affine motion, i.e., eqn (48)–(50), are satisfied, G′ and G′′ in the ordered MBS coincide with those in the TBS.

Appendix B: effect of particle rotation

In this section, we illustrate the effect of particle rotation, which was not described in Appendix A. In the model with rotation, the tangential velocity v^(t)_ij is given by

v^(t)_ij = (v_i − v_j)·t_ij − (d_iω_i + d_jω_j)/2 (54)

instead of eqn (39), where ω_i denotes the angular velocity of particle i. The time evolution of ω_i is given bywith the moment of inertia I_i = m_id_i²/8 and torque .

As shown in Fig. 9, we plot G′ in the ordered MBS with and without rotation with k_t/k_n = 1.0 and ε = 0.001 for various values of μ = 0.01. The values of other parameters are the same as those in Appendix A. The effect of particle rotation is negligible, except for the region near γ_c.

Fig. 9 Storage modulus G′ against γ₀ for the ordered MBS with μ = 0.01 and ε = 0.001. The open and filled symbols represent the results of the particles with and without rotation, respectively.

As shown in Fig. 10, we plot G′′ in the ordered MBS with and without rotation with k_t/k_n = 1.0 and ε = 0.001 for μ = 10⁻², 10⁻³, 10⁻⁴, 10⁻⁵ and 0.00001. The values of other parameters are the same as those in Appendix A. There are slight deviations in the peak position near γ_c between particles with and without rotation.

Fig. 10 Loss modulus G′′ against γ₀ for the ordered MBS with μ = 10⁻², 10⁻³, 10⁻⁴, 10⁻⁵ and ε = 0.001. The open and filled symbols represent the results of the particles with and without rotation, respectively.

In Fig. 11, we plot G′′ in the ordered MBS with and without rotation with k_t/k_n = 1.0 and ε = 0.001 for μ = 1.0 and 0.1. There are slight deviations in the peak position near γ_c between particles with and without rotation even for these higher μ. In addition, the second increase in G′′ around γ₀ = 0.1 for particles without rotation disappears for those with rotation.

Fig. 11 Loss modulus G′′ against γ₀ for the ordered MBS with μ = 1.0,0.1 and ε = 0.001. The open and filled symbols represent the results of the particles with and without rotation, respectively.

Appendix C: analytical calculation of shear stress and pressure

This section briefly explains the derivation of the normal and tangential components of shear stress and pressure for a small value of γ₀. From eqn (1)–(3), the relative distance r_ij is given by

r₂₃(θ(t)) = (−[small script l],0). (58)

Substituting these equations into u⁽ⁿ⁾_ij = r_ij − d, the normal displacements are given by

u⁽ⁿ⁾₂₃(t) = −εd. (61)

Substituting these equations into eqn (6), we obtain the normal force as

f⁽ⁿ⁾₂₃ = k_nεd (64)

up to O(γ₀).

By differentiating eqn (56)–(58) with time t, we obtain the relative velocity as

v₂₃(t) = (0,0) (67)

with the strain rate . The tangential unit vector is given by

t₂₃(t) = (0, −1). (70)

By considering the inner product of v_ij and t_ij, the tangential velocity is given by

v^(t)₂₃(t) = 0. (73)

If the transition from the stick state to the slip state does not occur under oscillatory shear, the tangential displacement is obtained by integrating v^(t)_ij(t) asu^(t)₂₃(t) = 0.(75)Substituting these equations into f^(t)_ij = −k_tu^(t)_ij yields

f^(t)₁₂ = f^(t)₁₃ = 3k_tγ(θ(t))[small script l]/4, (76)

f^(t)₂₃ = 0 (77)

up to O(γ₀). The condition that the transition does not occur is satisfied when f^(t)₁₂ < μf⁽ⁿ⁾₁₂ for γ = γ₀. Using eqn (64) and (77) with the assumption γ₀ ≪ ε, the condition is replaced by γ₀ < γ_c with γ_c given by eqn (16).

For γ₀ > γ_c, there exist regions where u^(t)_ij is unchanged in the slip state as

u^(t)₁₃ = u^(t)₁₂, (79)

u^(t)₂₃ = 0, (80)

where Θ satisfies

This equation provides Θ = cos⁻¹(1 − 2γ_c/γ₀). Substituting these equations into f^(t)_ij = −k_tu^(t)_ij yields

f^(t)₁₃ = f^(t)₁₂, (83)

f^(t)₂₃ = 0. (84)

The normal component of σ in eqn (9) is given by

σ⁽ⁿ⁾ = σ⁽ⁿ⁾₁₂ + σ⁽ⁿ⁾₁₃ (85)

withSubstituting eqn (58) and (57) with eqn (64) and (63) into eqn (86) and (87) and using eqn (85), we obtain σ⁽ⁿ⁾ as eqn (14).

The tangential component of σ in eqn (10) is given by

σ^(t) = σ^(t)₁₂ + σ^(t)₁₃ (88)

withSubstituting eqn (58) and (57) with eqn (77) into eqn (88), (89), and (90), we obtain σ^(t) as eqn (15) for γ₀ < γ_c. Using eqn (84) and (83) instead of eqn (77), we obtain σ^(t) as eqn (17) for γ₀ ≥ γ_c.

The pressure, i.e., P, in eqn (11) is defined as

P = P₁₂ + P₁₃ + P₂₃ (91)

withSubstituting eqn (56)–(58) with eqn (62)–(64) into eqn (91) and (92) with γ = 0, we obtain P₀(γ₀, μ) as eqn (21).

Appendix D: relation between shear modulus and stress–strain curve

In this section, we relate the shape of the stress–strain curve to the complex shear modulus. The shear stress σ(θ) is expanded using the Fourier series aswhere and with n > 1 denote the higher harmonics, , and . By neglecting and for n > 1,which is the maximum value of the scaled stress–strain curve illustrated in Fig. 3(b). This expression and the scaled stress–strain curve in Fig. 3(b) explain the decrease of G′ defined by eqn (18).

The area S of the curve for σ(θ)/γ₀ against γ(θ)/γ₀ is given by

Substituting eqn (4) into eqn (95) with eqn (12), we obtainwhich results in G′′ = S/π. As γ₀ increases, the area S of the scaled stress–strain curve in Fig. 3(b) increases first and decreases later, which explains the γ₀-dependence of G′′ provided by eqn (19).

Appendix E: details of disordered MBS

In this section, we present the details of the disordered MBS. This model is an extension of the monodisperse model used in Appendix A, including the dispersion of the particles and disordered initial configuration.

The system is bidisperse and includes an equal number of particles with diameters d and d/1.4. To simulate the disordered MBS, we randomly place the particles in a rectangular box with an initial packing fraction of ϕ_I = 0.75. The system is slowly compressed until the packing fraction reaches ϕ.²⁴ In each compression step, the packing fraction is increased by Δϕ = 1.0 × 10⁻⁴ with an affine transformation. Thereafter, the particles are relaxed to a mechanical equilibrium state with the kinetic temperature . Here, we choose T_th = 1.0 × 10⁻⁸k_nd². After compression, the oscillatory shear strain given by eqn (4) is applied for N_c cycles. In the last cycle, we measure G′ and G′′ using eqn (12) and (13) with eqn (8)–(10). The pressure, P₀(γ₀, μ) is obtained using eqn (11) after the last cycle. We use ϕ = 0.87, N = 1000, N_c = 20, L_y/L_x = 1, k_t = 0.2k_n, and .

Fig. 12 shows the storage modulus G′ against γ₀ in the disordered MBS for various values of μ. The storage modulus G′ is almost independent of γ₀ for a small γ₀ and decreases as γ₀ increases. The endpoint of the first plateau increases with μ except for μ = 0. A second plateau of G′ exists for μ = 10⁻⁴ and 10⁻⁵. The behavior of G′ for relatively small γ₀ is similar to that of crystalline solids as depicted in Fig. 4. On the other hand, the decrease of G′ for larger γ₀ cannot be captured by the analytical results of the TBS. Note that G′ for μ = 0.1 in the limit γ₀ → 0 is different from that for μ ≤ 0.01, which results from the μ dependence of the jamming point ϕ_J.⁹

Fig. 12 Storage modulus G′ in the disordered MBS against γ₀ with ϕ = 0.870 for various values of μ.

Fig. 13 shows the loss modulus G′′ in the disordered MBS against γ₀ for various values of μ. For sufficiently small γ₀, G′′ is zero, while G′′ becomes non-zero as γ₀ increases. The loss modulus G′′ starts to increase for smaller γ₀ as μ decreases. Similar to the case of G′, TBS captures only the behavior of relatively small γ₀ (see Fig. 5 and 13).

Fig. 13 Loss modulus G′′ in the disordered MBS against γ₀ with ϕ = 0.870 for various values of μ.

Appendix F: numerical shear modulus for TBS

In this section, we show the behaviors of G′ and G′′ in the TBS without the assumption used to obtain the analytical solution. Here, we numerically obtain G′ and G′′ under quasistatic oscillation using eqn (12) and (13) based on the left Riemann sum, where the integration of Ψ(θ), i.e.,is approximated aswith Δθ = 2π/M and θ_n = (n − 1)Δθ. We use ε = 0.001 and Δθ = 5.0 × 10⁻⁵ in our simulation.

As shown in Fig. 14, we plot the storage modulus G′ numerically obtained from the TBS against γ₀ with k_t/k_n = 1.0 for various values of μ as points. Moreover, we plot the analytical results derived from eqn (18) as thin solid lines. The numerical results agree with the analytical results for γ₀ < 0.003 and reproduce the second plateau of the MBS shown in Fig. 4.

Fig. 14 Storage modulus G′ against γ₀ with k_t/k_n = 1.0 and ε = 0.001 for various values of μ. The points represent the numerical results of the TBS, while the thin solid lines represent the analytical result given by eqn (18). The vertical dashed lines represent the critical amplitude γ_c(μ) given by eqn (16) for μ = 10⁻⁴, 10⁻³, 10⁻², 10⁻¹, 10⁰ from left to right.

Fig. 15 shows the loss modulus G′′ numerically obtained from the TBS against γ₀ with k_t/k_n = 1.0 for various values of μ as points. We also plot the analytical results given by eqn (19) as thin solid lines. The numerical results agree with the analytical results for γ₀ < 0.003.

Fig. 15 Loss modulus G′′ against γ₀ with k_t/k_n = 1.0 and ε = 0.001 for various values of μ. The points represent the numerical results of the TBS, while the thin solid lines represent the analytical results obtained from eqn (19). The vertical dashed lines represent the critical amplitude γ_c(μ) obtained from eqn (16) for μ = 10⁻⁴, 10⁻³, 10⁻², 10⁻¹, 10⁰ from left to right.

Acknowledgements

The authors thank K. Saitoh, D. Ishima, and S. Takada for fruitful discussions. This study was supported by JSPS KAKENHI under Grant No. JP19K03670 and JP21H01006.

References

1. M. van Hecke, J. Phys.: Condens. Matter, 2010, 22, 033101.
2. R. P. Behringer and B. Chakraborty, Rep. Prog. Phys., 2019, 82, 012601.
3. C. Maloney and A. Lemaître, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 93, 195501.
4. C. E. Maloney and A. Lemaître, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 016118.
5. D. Ishima, K. Saitoh, M. Otsuki and H. Hayakawa, 2022, arXiv:2207.06632.
6. C. S. O'Hern, S. A. Langer, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2002, 88, 075507.
7. M. Wyart, Ann. Phys., 2005, 30, 1.
8. B. P. Tighe, Phys. Rev. Lett., 2011, 107, 158303.
9. M. Otsuki and H. Hayakawa, Phys. Rev. E, 2017, 95, 062902.
10. C. Coulais, A. Seguin and O. Dauchot, Phys. Rev. Lett., 2014, 113, 198001.
11. M. Otsuki and H. Hayakawa, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 042202.
12. K. H. Nagamanasa, S. Gokhale, A. K. Sood and R. Ganapathy, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 062308.
13. E. D. Knowlton, D. J. Pine and L. Cipelletti, Soft Matter, 2014, 10, 6931–6940.
14. T. Kawasaki and L. Berthier, Phys. Rev. E, 2016, 94, 022615.
15. P. Leishangthem, A. D. S. Parmar and S. Sastry, Nat. Commun., 2017, 8, 14653.
16. A. H. Clark, J. D. Thompson, M. D. Shattuck, N. T. Ouellette and C. S. O'Hern, Phys. Rev. E, 2018, 97, 062901.
17. J. Boschan, S. Luding and B. P. Tighe, Granular Matter, 2019, 21, 58.
18. J. Boschan, D. VÅgberg, E. Somfai and B. P. Tighe, Soft Matter, 2016, 12, 5450–5460.
19. D. Nakayama, H. Yoshino and F. Zamponi, J. Stat. Mech.: Theory Exp., 2016, 2016, 104001.
20. T. Kawasaki and K. Miyazaki, arXiv, 2020, preprint, arXiv:2003.10716,  DOI:10.48550/arXiv.2003.10716.
21. K. Hyun, M. Wilhelm, C. O. Klein, K. S. Cho, J. G. Nam, K. H. Ahn, S. J. Lee, R. H. Ewoldt and G. H. McKinley, Prog. Polym. Sci., 2011, 36, 1697.
22. S. Dagois-Bohy, E. Somfai, B. P. Tighe and M. van Hecke, Soft Matter, 2017, 13, 9036–9045.
23. D. Ishima and H. Hayakawa, Phys. Rev. E, 2020, 101, 042902.
24. M. Otsuki and H. Hayakawa, Eur. Phys. J. E: Soft Matter Biol. Phys., 2021, 44, 70.
25. M. Otsuki and H. Hayakawa, Phys. Rev. Lett., 2022, 128, 208002.
26. N. Goldenfeld, Lectures on Phase Transitions and the Renormalization Group, CRC Press, 1992.
27. J. E. Lennard-Jones and A. F. Devonshire, Proc. R. Soc. A, 1937, 163, 53–70.
28. J. E. Lennard-Jones and A. F. Devonshire, Proc. R. Soc. A, 1938, 165, 1–11.
29. F. Yonezawa and K. Morigaki, Prog. Theor. Phys. Suppl., 1973, 53, 1–76.
30. S. Feng, M. F. Thorpe and E. Garboczi, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 31, 276–280.
31. A. Awazu, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 63, 032102.
32. P. A. Cundall and O. D. L. Strack, Géotechnique, 1979, 29, 47–65.
33. M. Doi and S. F. Edwards, The theory of polymer dynamics, Oxford University Press, 1986.
34. D. J. Evans and G. Morriss, Statistical Mechanics of Nonequilibrium Liquids, Cambridge University Press, 2008.

---