Contact network changes in ordered and disordered disk packings

Philip J. Tuckman a, Kyle VanderWerf a, Ye Yuan bc, Shiyun Zhang cd, Jerry Zhang c, Mark D. Shattuck e and Corey S. O’Hern *acfg
aDepartment of Physics, Yale University, New Haven, Connecticut 06520, USA
bDepartment of Mechanics and Engineering Science, College of Engineering, Peking University, Beijing 100871, China. E-mail: yuanyepeking@pku.edu.cn
cDepartment of Mechanical Engineering and Materials Science, Yale University, New Haven, Connecticut 06520, USA
dDepartment of Physics, University of Science and Technology of China, Hefei, Anhui 230026, China. E-mail: zsy12@mail.ustc.edu.cn
eBenjamin Levich Institute and Physics Department, The City College of New York, New York, New York 10031, USA
fDepartment of Applied Physics, Yale University, New Haven, Connecticut 06520, USA. E-mail: corey.ohern@yale.edu
gGraduate Program in Computational Biology and Bioinformatics, Yale University, New Haven, Connecticut 06520, USA

Received 21st June 2020 , Accepted 1st September 2020

First published on 17th September 2020


We investigate the mechanical response of packings of purely repulsive, frictionless disks to quasistatic deformations. The deformations include simple shear strain at constant packing fraction and at constant pressure, “polydispersity” strain (in which we change the particle size distribution) at constant packing fraction and at constant pressure, and isotropic compression. For each deformation, we show that there are two classes of changes in the interparticle contact networks: jump changes and point changes. Jump changes occur when a contact network becomes mechanically unstable, particles “rearrange”, and the potential energy (when the strain is applied at constant packing fraction) or enthalpy (when the strain is applied at constant pressure) and all derivatives are discontinuous. During point changes, a single contact is either added to or removed from the contact network. For repulsive linear spring interactions, second- and higher-order derivatives of the potential energy/enthalpy are discontinuous at a point change, while for Hertzian interactions, third- and higher-order derivatives of the potential energy/enthalpy are discontinuous. We illustrate the importance of point changes by studying the transition from a hexagonal crystal to a disordered crystal induced by applying polydispersity strain. During this transition, the system only undergoes point changes, with no jump changes. We emphasize that one must understand point changes, as well as jump changes, to predict the mechanical properties of jammed packings.


1 Introduction

Granular materials, which are composed of macroscopic grains that interact via frictional contact forces, are ubiquitous in the natural world and industrial applications. Unless they are continuously driven, granular materials will come to rest and, when confined, they exist in jammed, solid-like states.1 The mechanical response of jammed granular materials is highly nonlinear, which gives rise to shear jamming,2 intermittency and avalanches,3,4 shear banding,5,6 and other collective behavior.7

Numerous theoretical and computational studies have focused on simplified descriptions of dry granular media, where they are modeled as packings of frictionless, purely repulsive spherical grains.8,9 These studies have provided significant insights into the jamming transition in packings of frictionless, spherical particles. Disordered packings of frictionless spherical particles are typically isostatic at jamming onset,10i.e. they possess the same number of interparticle contacts Nc as the number of non-trivial degrees of freedom: Nc = Nisoc, where Nisoc = dNd + 1 (for systems with periodic boundary conditions), N is the number of (non-rattler11) grains, and d = 2, 3 is the spatial dimension. Ordered or compressed jammed packings can be hyperstatic with NcNisoc.12 Each jammed packing exists in a local energy minimum in configuration space, and therefore possesses a percolating network of non-zero interparticle forces and nonzero bulk and shear moduli. In contrast, packings with fewer contacts than the isostatic value, Nc < Nisoc, are unjammed and all interparticle forces are zero.13 Several studies have shown that isostatic jammed packings possess unique structural and mechanical properties, such as an excess number of low-frequency vibrational modes above the Debye prediction for the density of states14,15 and the power-law scaling of the shear modulus with increasing pressure.16

In prior studies, we considered jammed packings of frictionless, spherical particles undergoing quasistatic deformation (i.e. steps of applied simple or pure shear strain with each step followed by energy minimization).17 During quasistatic deformation, grains in the packings undergo continuous motions along “geometric families”, in which the network of interparticle contacts does not change.18,19 The continuous geometric families are punctuated by particle rearrangements, which cause the contact networks to change. Such rearrangements determine the structural and mechanical properties of jammed packings. For example, particle rearrangements control the power-law scaling of the ensemble-averaged shear modulus as a function of pressure during isotropic compression.20 Prior studies of sheared particulate materials have shown that there are two types of changes in the contact networks.21 We refer to these contact network changes as (1) jump changes and (2) point changes. These previous studies also found that the relative frequency of jump and point changes is roughly constant with increasing system size.

In this work, we further investigate jump and point changes in the contact network and show that these two types of contact network changes occur during a wide range of quasistatic deformations in model granular materials. We carry out discrete element method simulations of purely repulsive, frictionless disks in 2D, focusing on several types of quasistatic deformations: simple shear strain, changes in the size polydispersity of the grains, and isotropic compression. For jump changes, jammed packings become mechanically unstable during quasistatic deformation,22 the particles rearrange, and as a result, the total energy, pressure, shear stress, and other thermodynamic quantities are discontinuous at the strain where the particle rearrangement occurs.23 At a point change, a contact is added or removed from the interparticle contact network at a given strain, but the particles do not move significantly. The positions of the particles are continuous with strain, but the derivatives of the particle positions with respect to strain are discontinuous. As a result, for point changes, the potential energy (in the case of strain applied at fixed packing fraction) or enthalpy (in the case of strain applied at fixed pressure) and their first derivatives are continuous as a function of strain.24 For repulsive linear spring interactions, second- and higher-order derivatives of the potential energy/enthalpy are discontinuous at a point change, while for Hertzian spring interactions, third- and higher-order derivatives of the potential energy/enthalpy are discontinuous. We illustrate the importance of point changes by starting with a perfectly ordered jammed disk packing, adding small increments of size polydisperity to the system, and minimizing the potential energy (at fixed packing fraction) or enthalpy (at fixed pressure). This system undergoes a series of point changes as it proceeds from a hyperstatic toward an isostatic state.25,26

The remainder of the article is organized as follows. In Section 2, we describe the numerical methods that we use to generate disk packings at jamming onset and that we use to deform the jammed packings. In Section 3, we show results for the coordination number (z = 2Nc/N), total potential energy, shear stress, pressure, and other thermodynamic properties of jammed packings as a function of strain for each type of deformation, which allows us to illustrate point and jump changes. These studies are performed for both ordered packings of monodisperse disks and disordered packings of polydisperse disks. In Section 4, we summarize the conclusions and provide several possible future research directions including determining how point and jump changes separately contribute to the power-law scaling of the shear modulus with pressure during isotropic compression and investigating the effects of point changes in disk packings that interact via repulsive Hertzian spring interactions27 and in jammed systems containing frictional and non-spherical particles.

2 Methods

We consider packings of N circular disks in rectangular cells with area A = LxLy and periodic boundary conditions in both the x- and y-directions. We study packings of monodisperse disks, for which there is significant positional order, as well as disordered packings of polydisperse disks. The monodisperse disk packings possess jammed packing fractions near the value for the hexagonal lattice, ϕx = 0.907, whereas the disordered polydisperse disk packings possess jammed packing fractions ϕJ ≈ 0.81–0.84.

The disks interact via the following purely repulsive pair potential:

 
image file: d0sm01137a-t1.tif(1)
where ε is the characteristic energy scale of the repulsive interaction potential, the exponent of the interaction potential α = 2 for repulsive linear springs and α = 5/2 for “Hertzian” springs, rij is the center-to-center distance between disks i and j, [r with combining circumflex]ij = [r with combining right harpoon above (vector)]ij/rij, σij = (σi + σj)/2 is the average diameter of disks i and j, and the Heaviside function Θ(·) ensures that the interaction is nonzero only when the disks overlap (rij < σij). The total potential energy is given by image file: d0sm01137a-t2.tif. The repulsive force on disk i, arising from an overlap with disk j, is image file: d0sm01137a-t3.tif. Studies have shown that disks interacting via the purely repulsive potential in eqn (1) recapitulate the structural and mechanical properties of hard-sphere systems near jamming onset.28

Note that the Hertzian theory for the force between two contacting elastic spherical particles depends on the spatial dimension. The theory gives an exponent of α = 5/2 for the interaction energy between two elastic spheres in 3D and an exponent of α = 2 for the interaction between two parallel cylinders,29 which can mimic interactions between elastic disks in 2D. Thus, formally, “Hertzian” interactions between elastic disks should consider α = 2 in 2D, not α = 5/2. However, our goal was to investigate the effect of variations of the power-law exponent in eqn (1) on contact changes. Thus, we study both α = 2 and 5/2 for disk packings in 2D, and refer to the 5/2 exponent as the “Hertzian” value since this is value of the exponent in 3D.8

To generate jammed packings, we first randomly place N disks in the simulation cell at small packing fraction ϕ0 ≈ 0.1. We set the particle diameters to be σi = 〈σ〉 + ηδi, where −0.5 ≤ δi/〈σ〉 ≤ 0.5 is uniformly distributed, 〈δi〉 = 0, image file: d0sm01137a-t4.tif is the standard deviation of the disk diameters, and image file: d0sm01137a-t5.tif defines the average diameter. For disordered packings, we employ a square box, whereas for crystalline packings, we employ a rectangular box with aspect ratio image file: d0sm01137a-t6.tif, which allows a hexagonal packing of contacting disks to fit in the simulation cell without any defects. We isotropically compress the system in small packing fraction steps, Δϕ, until the system develops a small nonzero pressure, image file: d0sm01137a-t7.tif. After each compression step, the total potential energy is minimized using the FIRE algorithm30 until the magnitude of the total net force on the disks, image file: d0sm01137a-t8.tif. We study the coordination number, total potential energy, pressure, shear stress, and elastic moduli in jammed packings as a function of the packing fraction and strain. We measure energy, stress, and force in units of ε, ε/〈σ2, and ε/〈σ〉, respectively.

To understand the effects of jump and point changes in the interparticle contact networks, we consider jammed disk packings undergoing several types of quasistatic deformations: (1) simple shear at constant packing fraction, (2) simple shear at constant pressure, (3) increments of increasing size polydispersity at constant packing fraction, (4) increments of increasing size polydispersity at constant pressure, and (5) isotropic compression.

2.1 Simple shear strain at fixed packing fraction

For simple shear deformations, the particle positions are transformed to (xi′,yi′) = (x0i + γLxy0i/Ly,y0i) consistent with Lees–Edwards boundary conditions, where (x0i,y0i) are the initial particle positions. After each simple shear strain step γ, we minimize the total potential energy at constant packing fraction until the system is in force balance, such that image file: d0sm01137a-t9.tif.

During the simple shear strain deformation, we calculate several quantities as a function of γ including the shear stress,

 
image file: d0sm01137a-t10.tif(2)
which becomes
 
image file: d0sm01137a-t11.tif(3)
for repulsive linear spring interactions (α = 2 in eqn (1)) (where yij = yiyj, xij = xixj, and dxij/dγ = yijLx/Ly),31 and the shear modulus,
 
image file: d0sm01137a-t12.tif(4)

The shear modulus can be decomposed into the affine and nonaffine contributions,32Gγ = Gaγ + Gnaγ, respectively. To calculate Gaγ, we assume that all particles move according to the affine deformation, (xi′,yi′) = (x0i + γLxy0i/Ly,y0i). Gnaγ includes the nonaffine particle motion in response to potential energy minimization at fixed packing fraction and boundary strain. For repulsive linear spring interactions (α = 2 in eqn (1)), the affine contribution to the shear modulus can be calculated analytically,

 
image file: d0sm01137a-t13.tif(5)
We monitor U, Σγ, Gγ, and Gaγ before and after jump and point changes during the applied simple shear strain.

2.2 Simple shear strain at fixed pressure

We also apply quasistatic simple shear strain as described in Section 2.1, except at constant pressure. At each strain increment, we set the target pressure pt and minimize the enthalpy, H = U + ptA. After each strain step, we terminate the minimization when image file: d0sm01137a-t14.tif. Minimizing the enthalpy ensures that we can maintain constant pressure pt as the system is strained. At each strain step, we measure the enthalpy and its derivative dH/dγ with respect to shear strain, and monitor jump and point changes in the interparticle contact network.

2.3 Polydispersity strain at fixed packing fraction

In this section, we describe simulations in which we start the system with monodisperse (η = 0) or polydisperse disks (ησ〉 = 0.08), and increase η in small steps Δη ∼ 10−5 to increase the polydispersity of the disks. After each increment, Δη, we reset the packing fraction to its desired value and minimize the total potential energy at constant packing fraction. We measure the “polydispersity stress” as a function of η,
 
image file: d0sm01137a-t15.tif(6)
which becomes
 
image file: d0sm01137a-t16.tif(7)
for the repulsive linear spring potential, and the associated elastic modulus,
 
image file: d0sm01137a-t17.tif(8)

As discussed for applied simple shear strain, Gη can also be decomposed into the affine and nonaffine contributions: Gη = Gaη + Gnaη. For repulsive linear spring interactions, the affine contribution can be calculated analytically, which becomes

 
image file: d0sm01137a-t18.tif(9)
We measure U, Ση, Gη, and Gaη as a function of η and identify jump and point changes in the interparticle contact network.

2.4 Polydispersity strain at fixed pressure

To increase the polydispersity at fixed pressure, we take small steps in η and minimize the enthalpy after each step until image file: d0sm01137a-t19.tif. During the applied strain, we measure the enthalpy, its derivative dH/dη, and changes in the interparticle contact network.

2.5 Isotropic compression

We also study the response of jammed packings to isotropic compression. We compress the system by decreasing the box size in both dimensions by −2ΔL/L = Δϕ/ϕ. At the same time, we transform the particle coordinates by xi′ = x0iΔL/L and yi′ = y0iΔL/L. After each compression step, we minimize the total potential energy until force balance is achieved. We measure the pressure,
 
image file: d0sm01137a-t20.tif(10)
which becomes
 
image file: d0sm01137a-t21.tif(11)
for repulsive linear spring interactions, and the bulk modulus,
 
image file: d0sm01137a-t22.tif(12)
B can be decomposed into the affine and nonaffine contributions: B = Ba + Bna, respectively. For repulsive linear spring interactions, the affine contribution can be calculated analytically,
 
image file: d0sm01137a-t23.tif(13)
We calculate B and Ba as a function of packing fraction and monitor changes in the contact network during isotropic compression. Results for changes in the total poential energy and bulk modulus across point and jump changes in the contact network for disk packings with repulsive linear spring interactions are provided in Appendix I.

3 Results

In this section, we present the results for the energy, stress, and elastic moduli for the five applied deformations described in Section 2. We first show that changes in the interparticle contact networks during applied strain are either point changes or jump changes. For a jump change, the positions of the particles are discontinuous at the particular strain where the system becomes mechanically unstable and a particle rearrangement occurs. In contrast, for a point change, an interparticle contact either breaks or a new contact forms as the particles move continuously during the applied strain. We show that at a point change the derivative of the particle motions with respect to strain are discontinuous as are the derivatives of the potential energy/enthalpy, but at an order that depends on the interparticle potential. At small, but nonzero pressure, point changes occur in pairs over a range in strain. The first point change involves the formation of a new contact and the second involves the breaking of an existing contact. The difference in strain between these point changes decreases with pressure, and thus the pair of point changes coincide in the zero-pressure limit. To illustrate their importance, we detect exclusively point changes as we add polydispersity to originally monodisperse, ordered disk packings. Lastly, we present the statistics for jump and point changes for polydispersity strain applied at fixed packing fraction.

3.1 Jump changes

We define a jump change as a change in the interparticle contact network for which the particle positions as a function of applied strain are discontinuous, i.e. the particles rearrange. The origin of the discontinuous particle motion stems from strain-induced changes in the energy or enthalpy landscape22,33 and is illustrated in Fig. 1 for a disk packing undergoing polydispersity strain at fixed pressure (Section 2.4). In Fig. 1(a)–(c), we show the disk configurations before, during, and after a jump change. The enthalpy H as a function of the polydispersity strain η and the distance λ along the path from the initial state before the jump change to the final state after the jump change is shown in Fig. 1(d). To calculate H(η,λ), we define a vector [small xi, Greek, vector] = (Lx,x1,…,xN,y1,…yN) that contains all of the degrees of freedom of the packing. If the path that the system takes from point b1 to b2 in Fig. 1(d) is given by [small xi, Greek, vector](η*,λ), where the jump change occurs at η*, 0 < λ < 1, and Δ[small xi, Greek, vector](η) = [small xi, Greek, vector](η,1) − [small xi, Greek, vector](η,0), then [small xi, Greek, vector](η,λ) = Δ[small xi, Greek, vector](η)(([small xi, Greek, vector](η*,λ) − [small xi, Greek, vector](η*,0))/Δ[small xi, Greek, vector](η*). λ parametrizes the path that the system takes in configuration space during enthalpy minimization from the initial state at b1 (λ = 0) to the final state at b2 (λ = 1). The system is strained by increasing η in small steps (moving from left to right), followed by enthalpy minimization (moving vertically). The system begins in the upper left region of the landscape (point a), moves to the right (increasing η), and is initially prevented from moving toward the deeper minimum at the bottom of the enthalpy landscape by a barrier. As the system is further strained, the enthalpy barrier shrinks until the system reaches point b1, where the barrier disappears, and the system evolves toward point b2 with lower enthalpy. The disappearance of the enthalpy barrier at a given strain gives rise to the discontinuous change in the particle positions. We then continue straining the system until it reaches point c. We find similar behavior for jump changes in the enthalpy landscape for systems undergoing simple shear strain at fixed pressure and in the energy landscape for systems undergoing simple shear strain or polydispersity strain at fixed packing fraction.
image file: d0sm01137a-f1.tif
Fig. 1 An example of a jump change in a disordered packing of N = 6 polydisperse disks during applied polydispersity strain at fixed pressure (Section 2.4). Panels (a)–(c) show the system, before, on both sides of, and after the change in the contact network. In (a) and (c), the red solid circles outline the particles, while the blue dashed lines represent the interparticle contact network. The arrows give the direction of the particle motion at the given value of strain. In (b), the blue solid circles (red dashed circles) represent the disk configuration and the blue dotted lines (red dashed lines) give the contact network after (before) the change. From the arrows and circles, we see that both the particle positions and directions of motion are discontinuous at the jump change. (d) The enthalpy H (increasing from blue to red) is plotted as a function of the polydispersity η and distance λ along the path from the initial to the final state. The system starts in the upper left of the enthalpy landscape in the configuration in (a). At every η, the system can move vertically as long as the enthalpy decreases. The system is strained (increasing η) until it reaches point b1, corresponding to the configuration in panel (b) with red dashed lines. After reaching b1, a path to b2 (the configuration in panel (b) with blue solid lines) opens and the system can reach a deeper local minimum without an increase in enthalpy during the trajectory. λ parametrizes the distance along this path from b1 to b2. The system is then strained until point c, which corresponds to the configuration in panel (c). The bold black lines with arrows indicate the path taken by the system.

3.2 Point changes

We define a point change as the addition or removal of an interparticle contact at a given strain without discontinuous motion of the particles. The origin of a point change is that the positions of all particles for two or more distinct interparticle contact networks are the same at a given strain. In Fig. 2, we illustrate two successive point changes for a disordered disk packing undergoing polydispersity strain at fixed pressure (Section 2.4). In panels (a) and (b), we show the disk configurations corresponding to a point change from an isostatic packing to a hyperstatic packing with one additional contact, and in panels (b) and (c), we show the disk configurations corresponding to a point change from the same hyperstatic packing to a different isostatic packing. The arrows indicate the direction of motion of the particles, which show that the directions of the particle motion are not continuous over a point change. In Fig. 2(d), we show the enthlapy of the configurations in panels (a)–(c) as a function of strain η for target pressure pt = 10−4. We assume that (in the absence of changes in the contact network) the direction of the particle motion is constant with strain to extrapolate H for the contact networks that are not enthalpy minima.
image file: d0sm01137a-f2.tif
Fig. 2 An example of an N = 8 polydisperse disk packing undergoing two successive point changes during applied polydispersity strain η at fixed pressure. Panels (a) and (b) illustrate the first point change from an isostatic packing to a hyperstatic packing (with one extra contact) and (b) and (c) illustrate the second point change from the same hyperstatic packing to a different isostatic packing. All three packings are at target pressure pt = 10−4. The arrows indicate the directions of particle motion at each strain. The number 1 (2) labels the interparticle contact that is removed (added) during the two point changes. (d) Enthalpy H plotted versus η for the isostatic (hyperstatic) contact networks indicated by solid (dashed) lines. η1* (η2*) labels the strain at which a contact is added (removed) from the contact network.

At small η, the isostatic network in Fig. 2(a) has the lowest enthalpy of the three contact networks. At 1.190 < η1* < 1.191, H of the configuration in (b) becomes less than that of the configuration in (a), and the system becomes hyperstatic with an additional interparticle contact. At a higher strain 1.191 < η2* < 1.192, H for the configuration in (c) becomes less than that of the configuration in (b), and the system transitions to a different isostatic contact network. Most importantly, the particle positions do not change discontinuously during each point change. In other words, the contact change happens between two energy minimized configurations. In contrast, for jump changes, as shown in Fig. 1(d), the contact change occurs between a non-minimized configuration (point b1) and a minimized configuration (point b2).

The changes of the particle trajectories in Fig. 2(a)–(c) demonstrate the importance of point changes. If contact 2 did not form in panel (b), the two particles that form that contact would continue to move towards each other as they do in panel (a). These particle trajectories would cause a dramatic increase in enthalpy, as shown by H(η) for the first isostatic contact network in panel (d). However, due to the formation of the new contact, the particle trajectories are altered following the point change as shown in panel (c). Despite the continuous particle motion that occurs during point changes, the particle trajectories are significantly altered with further strain.

Fig. 3 displays the values of the polydispersity strain η1* (η2*) at which several example polydisperse N = 8 packings transition from an isostatic packing to a hyperstatic packing (and from the same hyperstatic packing to an isostatic packing) as a function of the target pressure pt. For each packing, we find that both η1* and η2* are linear in pt with vertical intercept η0 = η1,2*(pt = 0). In Fig. 3, we show that the values of η1,2*, corresponding to when the packing either gains a contact or loses a contact, possess the same η0. Thus, the width of the strain region over which the system is hyperstatic between the two successive point changes (first from an isostatic packing to a hyperstatic packing and then from the same hyperstatic packing to another isostatic packing) tends to zero in the zero-pressure limit. We find similar behavior for disk packings undergoing simple shear strain, as well as for larger system sizes.


image file: d0sm01137a-f3.tif
Fig. 3 For three sample polydisperse N = 8 packings, we measure the polydispersity strain values at which the system transitions from an isostatic to a hyperstatic packing (η1*) and from the same hyperstatic packing to another isostatic packing (η2*) as shown in Fig. 2, at 10 target pressures pt. We plot η1,2*–η0, where η0 = η1,2*(pt = 0), versus pt for each contact change in each packing. The strain at which the packings transition from isostatic to hyperstatic, (i.e. between Fig. 2(a) and (b)), are represented by blue diamonds, red rightward triangles, and green downward triangles. The strains at which the packings transition from hyperstatic to isostatic, (i.e. between Fig. 2(b) and (c)) are represented by blue squares, red leftward triangles, and green upward triangles. Since all of the lines meet at η1,2* = η0, the width of the hyperstatic strain region tends to zero in the pt = 0 limit.

3.3 Generalization of jump and point changes to other strains

While the illustrations of jump and point changes in Sections 3.1 and 3.2 considered polydispersity strain at constant pressure, all interparticle contact changes that occur during the applied strains that we consider (i.e. simple shear strain at constant packing fraction and at constant pressure, polydispersity strain at constant packing fraction and at constant pressure, and isotropic compression) can be classified as jump or point changes. Further, we find that a point change at a given strain gives rise to continuous potential energy/enthalpy and its first derivatives, but causes discontinuities in the second derivatives of the potential energy/enthalpy at the given strain. The fact that the second derivatives of the potential energy/enthalpy are discontinuous (eqn (5)) is related to the repulsive linear spring interparticle potential that we employ; results for other finite-range repulsive potentials are discussed in Section 3.5. In contrast, all jump changes give rise to discontinuities in the potential energy/enthalpy, as well as all of its derivatives, independent of the interparticle potential.

As an example, in Fig. 4, we show the enthalpy H as a function of simple shear strain γγ* for an N = 8 packing (with repulsive linear spring interactions) undergoing simple shear at fixed pressure. For the jump change at γ*, H is discontinuous. For the point change at γ*, H and dH/dγ (in the inset) are both continuous, but d2H/dγ2 is discontinuous. The fact that the second derivative of the enthalpy, Gγ + pt(d2V/dγ2), is discontinuous at a point change can be illustrated by analyzing the affine contribution of the shear modulus, Gaγ in eqn (5), when contacts with zero overlap, rijσij, are added to or removed from the contact network. For the same reason, point changes give rise to discontinuities in the second derivatives with respect to strain of the potential energy/enthalpy for disk packings with repulsive linear spring interactions undergoing other applied strains.


image file: d0sm01137a-f4.tif
Fig. 4 The enthalpy HJ (left vertical axis) and HP (right vertical axis) for a N = 16 disk packing versus simple shear strain γγ* at constant pressure (p = 10−4) for a jump change (blue downward triangles) and a point change (red upward triangles) in the contact network at γ*. For the jump change, there is a discontinuity in H at γ*. For the point change, both the enthalpy and its first derivative dH/dγ (inset) are continuous at γ*. However, the slope of dH/dγ changes at γ*, which indicates that d2H/dγ2 is discontinuous.

3.4 Packing fraction-strain landscapes

We refer to jammed disk packings with the same contact network as geometrical families (continuous regions) in the packing fraction and applied strain plane.19 One can then consider contours of constant stress in the packing fraction and strain plane for each distinct contact network, and identify point and jump changes by calculating derivatives of the stress. In this section, we study the packing fraction and strain landscapes for both simple shear strain and polydispersity strain. To construct these landscapes, we first generate a series of unjammed packings (with ϕ ≈ 0.8) over a range of strains. We find similar results using other packing fractions ϕϕJ. We then isotropically compress these packings (quasistatically) at each strain to packing fractions above jamming onset. For the disk packings at each packing fraction and strain, we measure the contact network, coordination number, and stress. This protocol ensures that we can sample packings with both signs of the shear stress.17 For clarity, we show only a small portion of the strain-packing fraction landscape.

In Fig. 5, we visualize polydisperse N = 8 disk packings in the packing fraction ϕ and simple shear strain γ plane. The color of a region indicates the type of contact network: regions that are red indicate isostatic contact networks and regions that are green indicate hyperstatic contact networks. Regions with different hues of red and green correspond to different contact networks. The white regions represent unjammed states. The lines provide contours of constant shear stress Σγ. Σγ is discontinuous at jump changes, whereas it is continuous at point changes.


image file: d0sm01137a-f5.tif
Fig. 5 An example landscape in the packing fraction ϕ and simple shear strain γ plane for polydisperse N = 8 disk packings. We first apply simple shear strain (quasistatically) at ϕ = 0.79 (below jamming onset) to generate a series of unjammed packings over a range of γ with step size Δγ = 10−4. We then apply isotropic compression (quasistatically) with step size Δϕ = 10−4 to these packings at each strain to packing fractions above jamming onset. For each ϕ and γ, we show the contact network (color) and shear stress Σγ, where the lines are contours of constant Σγ and the difference between adjacent lines is ΔΣγ ≈ 7 × 10−5. Jump changes can be identified by discontinuities in Σγ, and point changes by discontinuities in the derivative of Σγ. Red regions indicate isostatic contact networks, green regions indicate hyperstatic contact networks, and white regions indicate unjammed systems. Each region with a distinct red or green hue indicates packings with the same contact networks.

The ϕγ landscape in Fig. 5 has two lines of point changes, which can be traversed by compressing or decompressing the packing at fixed γ, by applying simple shear strain at fixed ϕ, or by a combination of changes in ϕ and γ. The packing undergoes a point change when a contact is added (i.e. transitioning from an isostatic packing to a hyperstatic packing) or a contact is removed (i.e. transitioning from a hyperstatic packing to an isostatic packing). As discussed in Section 3.2, the two lines of point changes merge into a single point near (0.04, 0.81) in the zero-pressure limit. Traversing a point change in the forward direction leads to the same behavior as traversing it in the reverse direction.

Lines of jump changes in Fig. 5 occur when moving from an isostatic jammed region to an unjammed region. As we found for point changes, jump changes can be induced by compressing the packing at fixed γ, by applying simple shear strain at fixed ϕ, or by a combination of changes in ϕ and γ. When undergoing a jump change to an unjammed state, the total potential energy and shear stress drop discontinuously from a finite value to zero. In Fig. 5, there is also a line of jump changes between two different isostatic packings near (0.015, 0.81).

Note that in Fig. 5, the system can transition from a jammed packing to unjammed packing through isotropic compression. Indeed, in recent computational studies, we showed that “compression unjamming” occurs frequently near jamming onset. We also showed that the probability for compression unjamming (averaged over a finite range of strain) approaches a finite value in the large-system limit, and thus compression unjamming occurs in the large-system limit.20

In Fig. 6, we show a portion of the packing fraction and polydispersity strain landscape for N = 8 disk packings. The lines provide contours of constant polydispersity stress Ση. Ση is discontinuous at jump changes, whereas it is continuous at point changes. In Fig. 6, there are two lines of point changes, which can be traversed by compressing or decompressing the packing at fixed η, by applying polydispersity strain at fixed ϕ, or by a combination of changes in ϕ and η. Again, the two lines of point changes merge into a single point near (0.135, 0.815) in the zero-pressure limit. We find one line of jump changes in Fig. 6 that can cause a transition between two isostatic packings, between a hyperstatic and an isostatic packing, and between two hyperstatic packings.


image file: d0sm01137a-f6.tif
Fig. 6 An example landscape in the packing fraction ϕ and polydispersity strain η plane for N = 8 disk packings. We first apply polydispersity strain (quasistatically) at ϕ = 0.81 (below jamming onset) to generate unjammed packings over a range of η with step size Δη = 5 × 10−5. We then apply isotropic compression (with successive steps Δϕ10−5 followed by energy minimization) to these packings at each strain to packing fractions above jamming onset. For each ϕ and η, we show the contact network (color) and stress Ση, where the lines are contours of constant Ση and the difference between adjacent lines is ΔΣη ≈ 2.5 × 10−4. Jump changes can be identified by discontinuities in Ση, and point changes by discontinuities in the derivative of Ση. Red regions indicate isostatic contact networks, green regions indicate hyperstatic contact networks, and white regions indicate unjammed systems. Each region with a distinct red or green hue indicates packings with the same contact networks.

3.5 Hertzian spring interactions

In this section, we show preliminary results for frictionless disk packings that interact via repulsive Hertzian spring interactions (α = 5/2 in eqn (1)) undergoing simple shear strain at fixed packing fraction. In Fig. 7, we plot Gγversus γ for a N = 16 disk packing with repulsive Hertzian spring interactions across a point change. Gγ is continuous across the point change, but dGγ/dγ is discontinuous. This result can be anticipated by analyzing the affine contribution to the shear modulus for repulsive Hertzian spring interactions,
 
image file: d0sm01137a-t24.tif(14)
Gaγ for repulsive Hertzian spring interactions is similar to that for repulsive linear spring interactions (eqn (5)), but it has an additional factor of image file: d0sm01137a-t25.tif. Thus, when a new contact is added to or removed from the contact network (at rij = σij) during the applied strain, we expect that Gγ will be continuous. If we take an additional derivative of Gaγ with respect to γ, the factor of image file: d0sm01137a-t26.tif moves to the denominator, and thus we expect that dGγ/dγ will be discontinuous across point changes, as shown in Fig. 7.

image file: d0sm01137a-f7.tif
Fig. 7 Shear modulus Gγ as a function of simple shear strain γ (blue solid line) at fixed packing fraction ϕ = 0.88 for a N = 8 disk packing with repulsive Hertzian spring interactions. The point change in the contact network (vertical dotted black line at γ* ≈ 0.24) does not cause a discontinuity in Gγ, but does cause a discontinuity in dGγ/dγ.

3.6 Transition from a hexagonal crystal to a disordered crystal

To illustrate the importance of point changes, we investigate the transition from a hexagonal crystal to a disordered crystal25,26,34 as a function of applied polydispersity strain in disk packings with repulsive linear spring interactions. The disordered crystal has properties in common with a hexagonal crystal (such as the disk positions and packing fraction), whereas other properties, such as the coordination number, stress, and elastic moduli, are similar to disordered, isostatic packings. Here, we show that the transition from the hexagonal crystal to the disordered crystal can be understood as series of point changes as a function of polydispersity strain, with no jump changes. We note that the transition to the disordered crystal can also be induced by simple shear and other applied strains.

In Fig. 8, we plot the ensemble-averaged excess coordination number 〈zziso〉, where ziso = 2Nisoc/N, as a function of polydispersity strain η at fixed pt. 〈zziso〉 ≈ 2 at small η, and then begins to decrease toward zero at a characteristic ηc. As shown in the inset to Fig. 8, ηcpt since 〈zziso〉 collapses when plotted versus η/pt. Thus, in the zero-pressure limit, the hexagonal crystal at ϕ = ϕx becomes isostatic with z = ziso in the limit of zero applied strain.


image file: d0sm01137a-f8.tif
Fig. 8 The ensemble-averaged coordination number 〈zzisoversus the polydispersity strain η at constant pressure pt for N = 64 packings at pt = 10−4.5 (red solid line), 10−4.25 (green dashed line), 10−4 (blue dot-dashed line), and 10−3.75 (black dotted line). The system was initialized in a hexagonal crystal at η = 0. The inset shows that 〈zziso〉 can be collapsed by plotting it against η/pt. The data was averaged over 10 packings for each pressure.

We find similar behavior for the transition from a hexagonal crystal to a disordered crystal when we apply polydispersity strain at fixed packing fraction. In Fig. 9, we plot the total potential energy U and elastic modulus Gηversus η at fixed ϕ for an N = 16 packing initialized in a hexagonal crystal. We show that at each change in the contact network U is continuous, but Gη is discontinuous, which signals that the changes in the contact network are point changes. In Fig. 10, we show the ϕη landscape for an N = 16 packing initialized in a hexagonal crystal. There are many contact networks near the hexagaonal crystal, which are separated by point changes since there are no discontinuities in the polydispersity stress Ση. In the zero-pressure limit, all of the point changes coincide and the system transitions from a hexagonal network to an isostatic network at zero strain.


image file: d0sm01137a-f9.tif
Fig. 9 The total potential energy U (left vertical axis) and elastic modulus Gη (right vertical axis) versus polydispersity strain η at fixed packing fraction for a N = 16 packing initialized in a hexagonal crystal at pressure p = 10−4. At each change in the contact network (black dashed vertical lines), U (blue upward triangles) is continuous, while Gη (red downward triangles) is discontinuous.

image file: d0sm01137a-f10.tif
Fig. 10 The packing fraction ϕ and polydispersity strain η landscape for a N = 16 packing initialized in a hexagonal crystal. The color indicates the coordination number, ranging from isostatic with ziso ∼ 4 to crystalline with z = 6 (from blue to red). The white region corresponds to unjammed systems. The lines represent contours of constant polydispersity stress Ση and the difference between adjacent lines is approximately ΔΣη = 2 × 10−4. All of the changes in the contact networks are point changes, since there are no discontinuities in Ση.

3.7 Distinguishing point and jump contact changes

In this section, we discuss the changes in the total potential energy and elastic moduli that occur at point and jump changes for packings undergoing polydispersity strain at constant packing fraction. (Results for isotropic compression are provided in Appendix I.) In Fig. 11, we show a scatter plot of the absolute values of the changes in total potential energy |ΔU| and polydispersity modulus |ΔGη| at polydispersity strains that correspond to changes in the contact network. We also compare these values of |ΔU| and |ΔGη| to those obtained from successive polydispersity strains where there is no change in the contact network. (In this section, to illustrate the effects of changes in the contact network, we focus on a small system with N = 16 disks. In Appendix II, we show that the results are similar for much larger systems.) We find three distinct clusters of points: jump changes (with |ΔU| > 10−9 and large values of |ΔGη|), point changes (with |ΔGη| > 10−6 and small values of |ΔU|), and points with small values of |ΔU| and |ΔGη| where there are no changes in the contact network. This last set of points shifts to lower values of |ΔU| and |ΔGη| with decreasing Δη and improved force balance. All changes in the contact network during applied polydispersity strain can be classified as either point or jump changes. We find similar results for simple shear strain applied at fixed packing fraction and pressure, polydispersity strain applied at fixed pressure, and isotropic compression.
image file: d0sm01137a-f11.tif
Fig. 11 A scatter plot of the absolute values of changes in the potential energy |ΔU| and polydispersity modulus |ΔGη| between successive polydispersity strain steps Δη at constant packing fraction ϕ = 0.88 for 50 N = 16 packings. After every strain step, U and Gη were measured, and the difference between the values at the current step and the previous step was calculated to yield ΔU and ΔGη. The red triangles indicate a change in the contact network, whereas the black circles indicate strains where there was no change in the contact network.

In principle, one can also use particle displacements (i.e. nonaffine particle motion) to identify changes in the contact networks.35 For example, one could apply polydispersity strain from η1 to η2 yielding particle positions [r with combining right harpoon above (vector)](η1) and [r with combining right harpoon above (vector)](η2), and then reverse the strain from η2 to η1 to measure the new particle positions [r with combining right harpoon above (vector)]′(η1). The particle displacements Δr = |[r with combining right harpoon above (vector)](η1) − [r with combining right harpoon above (vector)]′(η1)| from this process will be large when there is a jump change between η1 and η2, whereas Δr → 0 (in the small strain limit) for strain intervals where there is no jump change. Thus, measuring non-affine particle motions cannot be used to identify point changes. For this reason, we recommend measurements of ΔG and ΔU to identify point and jump changes in particulate media.

4 Conclusions and future directions

In this article, we studied quasistatic deformations of jammed frictionless disk packings that interact via purely repulsive potentials as models of dense granular materials. The deformations included simple shear strain at fixed packing fraction and at fixed pressure, polydispersity strain at fixed packing fraction and at fixed pressure, and isotropic compression. We showed that there are two types of changes in the interparticle contact networks that occur during quasistatic deformation: point changes and jump changes. Jump changes involve changes in the contact network that are accompanied by discontinuous motion of the particles from one strain step to the next, whereas point changes involve small, continuous motion of the particles. It has been previously shown21 that the relative frequency of these two types of events is constant with increasing system size. Both types strongly affect the structural and mechanical properties of quasistatically deformed jammed granular systems. For jump changes, the total potential energy (when the deformation is applied at constant packing fraction), or the enthalpy (when the deformation is applied at fixed pressure), as well as their derivatives with respect to strain are discontinuous. In contrast, point changes give rise to discontinuities in higher-order derivatives with respect to strain of the potential energy/enthalpy. For example, for disk packings with repulsive linear spring interactions, point changes cause discontinuities in the elastic moduli, which are proportional to second-order derivatives with respect to strain of the potential energy (when the deformation is applied at constant packing fraction) or the enthalpy (when the deformation is applied at constant pressure). We then illustrated the important features of jump and point changes by showing contours of constant stress in the packing fraction and strain landscapes for the simple shear and polydispersity strain deformations. As a specific example of a physical phenomenon where point changes are dominant, we showed that the transition from a hexagonal crystal to a disordered crystal, which can possess an isostatic number of contacts, is caused by a series of point changes.

The fact that point changes cause discontinuities with respect to strain in the second derivative of the potential energy/enthalpy (for disk packings with repulsive linear spring interactions) stems from the shape of the interparticle potential energy (eqn (1)). The purely repulsive linear spring potential has a discontinuity in d2U/drij2 across a point change, and thus the elastic moduli, Gγ, Gη, and B, are discontinuous across a point change. For the purely repulsive Hertzian spring potential with α = 5/2 in eqn (1), d3U/drij3 is discontinuous across a point change, and thus the derivatives of the elastic moduli with respect to strain (not the moduli themselves) are discontinuous. The discontinuities caused by point changes will occur in higher-order derivatives of the potential energy (when the strain is applied at constant packing fraction) if higher-order derivatives of the interparticle potential are continuous. Similar results are found for the derivatives of the enthalpy when the strain is applied at fixed pressure.

These results raise several important questions for future research. First, how do jammed packings behave when the applied strain is reversed36–38 after point and jump changes occur in the interparticle contact networks? Point changes are completely reversible, since the particle motions are continuous during a point change. Jump changes, however, are not reversible in this way. As shown in Fig. 1, the packing immediately after the jump change has a lower potential energy (in the case of applied strain at constant packing fraction) than the packing immediately before the jump change. Thus, when the strain is reversed after the jump change, the system will follow a different path in the energy landscape (than the one followed during the forward strain). However, it is possible that the system can undergo a series of point changes or another jump change during the reversed strain and return to the path in the energy landscape that was traversed during the forward strain. This behavior was termed “loop reversibility” in ref. 39 and “limit cycle” behavior in ref. 40, both of which studied systems undergoing cyclic simple shear strain.

In recent studies,20 we found that changes in the contact network during isotropic compression of jammed packings give rise to the power-law scaling of the shear modulus with pressure, i.e. Gγp1/2 for repulsive linear spring interactions in d = 2 and 3. Since both point and jump changes cause jumps in the shear modulus, ΔGγ, an interesting question is to determine whether point changes, jump changes, or both contribute significantly to the increase in the shear modulus during isotropic compression. In addition, Gγp2/3 for Hertzian spring interactions undergoing isotropic compression in d = 2 and 3.8 In future studies, we will investigate how jump and point changes give rise to this behavior, given that point changes do not cause discontinuities in Gγ for Hertzian interactions.

To understand the mechanical response of jammed packings to applied strain, one must be able to predict the potential energy (and other physical quantities that depend on the particle positions) as the system evolves along geometrical families, as well as across point and jump changes. We emphasize that it is still important to study point changes in packings undergoing quasistatic deformation even if the interparticle potential does not possess discontinuities in its derivatives. Even if there are no discontinuities in the interparticle potential, the particle trajectories change directions when the system undergoes each point change, which influences the evolution of the potential energy, stress, and elastic moduli as a function of strain.

Another possible extension of the current studies is to investigate how point changes behave in packings of non-spherical particles. For example, in packings of circulo-lines in 2D, two particles with an “end–end” contact behave differently than two particles with an “end–middle” contact.41 It will be interesting to study packings of circulo-lines that transition between these two types of contacts and determine whether this process can be described as a generalized point change, even though the interparticle contact network does not change.

A similar effect can occur in packings of spherical particles with frictional interactions. Numerous studies have shown that in addition to the number of contacts per particle, the ratio of the tangential to the normal force, ζij, at each contact between particles i and j, plays an important role in determining the mechanical stability of frictional packings.42 Thus, it is possible that effective “point changes” can occur if ζij varies significantly during strain even though particles i and j remain in contact.

Conflicts of interest

There are no conflicts to declare.

Appendix I: isotropic compression

In this Appendix, we show that the results for isotropic compression are similar to the results for the other strains that we studied. In Fig. 12, we show a scatter plot of the absolute values of the changes in total potential energy |ΔU| and bulk modulus |ΔB| at compression values that correspond to changes in the contact network. We also compare these values of |ΔU| and |ΔB| to those obtained from successive compression steps where there is no change in the contact network. We find three distinct clusters of points: jump changes (with |ΔU| > 10−7 and large values of |ΔB|), point changes (with |ΔB| > 10−4 and small values of |ΔU|), and points with small values of |ΔU| and |ΔB| where there are no changes in the contact network. This last set of points shifts to lower values of |ΔU| and |ΔB| with decreasing compression step size and improved force balance. (See Appendix II.) All changes in the contact network during applied compression can be classified as either point or jump changes.
image file: d0sm01137a-f12.tif
Fig. 12 A scatter plot of the absolute values of changes in the potential energy |ΔU| and bulk modulus |ΔB| between successive compression steps Δϕ for 50 N = 16 packings. After every strain step, U and B were measured, and the difference between the values of the potential energy and bulk modulus at the current step and the previous step was calculated to yield ΔU and ΔB. The red triangles indicate a change in the contact network, whereas the black circles indicate strains where there was no change in the contact network.

Appendix II: system size dependence

In this Appendix, we show that the presence of point and jump changes and our method for distinguishing between them do not change with increasing system size. For most of the results in this article, we used small systems with N = 6 to 16 disks with periodic boundary conditions in the x- and y-directions. Since point and jump changes have not been described before in the literature, the main goal of this article is to illustrate the theoretical foundations of point and jump contact changes, not to provide statistics of point and jump changes in the large-system limit. In previous studies, it has been shown that the length of geometrical families decreases strongly with increasing system size,18 and thus it makes sense to illustrate point and jump changes using small systems, where one can clearly see the beginning and end of each family. Further, the threshold required on force balance on each particle necessary to identify point and jump changes decreases toward zero with increasing system size, and thus it is much less computationally costly to study point and jump changes in small systems.

Nevertheless, in Fig. 13, we show similar data as in Fig. 11, except for packings of N = 64, 128, and 256 disks undergoing simple shear (with step size Δγ = 7 × 10−13) at fixed packing fraction ϕ = 0.88. Again, we observe that there are three clusters of data points: one for jump changes (large |ΔU|/N and large |ΔGγ|), one for point changes (small |ΔU|/N and large |ΔGγ|), and one for the control group (small |ΔU|/N and small |ΔGγ|), for which point and jump changes do not occur. More importantly, we find that the location and spread of each of the three clusters remain the same for the three system sizes.


image file: d0sm01137a-f13.tif
Fig. 13 A scatter plot of the absolute values of changes in the potential energy per particle |ΔU|/N and shear modulus |ΔGγ| between successive shear steps Δγ = 7 × 10−13 for N = 64 (red upward triangles and black circles), N = 128 (green downward triangles and dark gray dots), and N = 256 (blue rightward triangles and light gray squares) packings. After every shear strain step, U and Gγ were measured, and the differences between the values at the current step and the previous step were calculated. The red, green, and blue triangles indicate a change in the contact network, whereas the black/gray points indicate strains where there was no change in the contact network.

In Fig. 14, we show the same plot as in Fig. 13 for the three system sizes N = 64, 128, and 256, except using a larger shear strain step size Δγ = 10−11. The data points for |ΔGγ| and |ΔU|/N corresponding to jump changes remain the same for the two shear strain step sizes. For the data points that correspond to point changes, the values of |ΔU|/N change with the shear strain step size, but the values of |ΔGγ| do not. In addition, for the points that do not correspond to changes in the contact network, both |ΔU|/N and |ΔGγ| shift to larger values with the larger shear strain step size. Thus, |ΔU|/N → 0 and |ΔGγ| → 0 in the limit Δγ → 0 for data points that do not correspond to changes in the contact network.


image file: d0sm01137a-f14.tif
Fig. 14 A scatter plot of the absolute values of changes in the potential energy per particle |ΔU|/N and shear modulus |ΔGγ| between successive shear steps Δγ = 10−11 for N = 64 (red upward triangles and black circles), N = 128 (green downward triangles and dark gray dots), and N = 256 (blue rightward triangles and light gray squares) packings. After every strain step, U and Gγ were measured, and the difference between the values at the current step and the previous step was calculated. The red, green, and blue triangles indicate a change in the contact network, whereas the black/gray points indicate strains where there was no change in the contact network.

Acknowledgements

We acknowledge support from NSF Grant No. CBET-1605178 (K. V. and C. O.), DBI-1755494 (P. T.), CBET-2002782 (C. O. and J. Z.), and CBET-2002797 (M. S.), National Natural Science Foundation of China Grants No. 11572004 (Y. Y.) and 11972047 (Y. Y.), and the China Scholarship Council Grant No. 201806010289 (Y. Y.) and 201906340202 (S. Z.). This work was also supported by the High Performance Computing facilities operated by Yale's Center for Research Computing and computing resources provided by the Army Research Laboratory Defense University Research Instrumentation Program Grant No. W911NF1810252.

References

  1. R. P. Behringer and B. Chakraborty, Rep. Prog. Phys., 2018, 82, 012601 CrossRef.
  2. D. Bi, J. Zhang, B. Chakraborty and R. P. Behringer, Nature, 2011, 480, 355 CrossRef CAS.
  3. J. Barés, D. Wang, D. Want, T. Bertrand, C. S. O'Hern and R. P. Behringer, Phys. Rev. E, 2017, 96, 052902 CrossRef.
  4. D. V. Denisov, K. A. Lörincz, J. T. Uhl, K. A. Dahmen and P. Schall, Nat. Commun., 2016, 7, 010641 CrossRef CAS.
  5. D. M. Mueth, G. F. Debregeas, G. S. Karczmar, P. J. Eng, S. R. Nagel and H. M. Jaeger, Nature, 2000, 406, 385 CrossRef CAS.
  6. K. Karimi and J.-L. Barrat, Sci. Rep., 2018, 8, 4021 CrossRef.
  7. I. S. Aranson and L. S. Tsimring, Rev. Mod. Phys., 2006, 78, 641 CrossRef CAS.
  8. 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 CrossRef.
  9. A. J. Liu and S. R. Nagel, Annu. Rev. Condens. Matter Phys., 2010, 1, 347 CrossRef.
  10. A. V. Tkachenko and T. A. Witten, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 60, 687 CrossRef CAS.
  11. F. Giacco, L. de Arcangelis, M. Pica Ciamarra and E. Lippiello, Soft Matter, 2017, 13, 9132 RSC.
  12. C. F. Schreck, C. S. O'Hern and L. Silbert, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 011305 CrossRef.
  13. T. Shen, C. S. O'Hern and M. D. Shattuck, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 011308 CrossRef.
  14. M. Wyart, S. R. Nagel and T. A. Witten, Europhys. Lett., 2005, 72, 486 CrossRef CAS.
  15. M. Wyart, L. E. Silbert, S. R. Nagel and T. A. Witten, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 051306 CrossRef.
  16. C. P. Goodrich, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2012, 109, 095704 CrossRef.
  17. S. Chen, T. Bertrand, W. Jin, M. D. Shattuck and C. S. O'Hern, Phys. Rev. E, 2018, 98, 042906 CrossRef CAS.
  18. G.-J. Gao, J. Blawzdziewicz and C. S. O'Hern, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 061303 CrossRef.
  19. T. Bertrand, R. P. Behringer, B. Chakraborty, C. S. O'Hern and M. D. Shattuck, Phys. Rev. E, 2016, 93, 012901 CrossRef.
  20. K. VanderWerf, A. Boromand, M. D. Shattuck and C. S. O'Hern, Phys. Rev. Lett., 2020, 124, 038004 CrossRef CAS.
  21. P. Morse, M. van Deen, S. Wijtmanns, M. van Heck and M. L. Manning, Phys. Rev. Res., 2020, 2, 023179 CrossRef CAS.
  22. D. L. Malandro and D. J. Lacks, J. Chem. Phys., 1999, 110, 4593 CrossRef CAS.
  23. Y. Cao, J. Li, B. Kou, C. Xia, Z. Li, R. Chen, H. Xie, T. Xiao, W. Kob, L. Hong, J. Zhang and Y. Wang, Nat. Commun., 2018, 9, 2911 CrossRef.
  24. H. Mizuno, L. Silbert, M. Sperl, S. Mossa and J.-L. Barrat, Phys. Rev. E, 2016, 93, 043314 CrossRef.
  25. C. P. Goodrich, A. J. Liu and S. R. Nagel, Nat. Phys., 2014, 10, 578 Search PubMed.
  26. H. Tong, P. Tan and N. Xu, Sci. Rep., 2015, 5, 15378 Search PubMed.
  27. C. F. Schreck, C. S. O'Hern and M. D. Shattuck, Granular Matter, 2014, 16, 209 Search PubMed.
  28. F. Arceri and E. I. Corwin, Phys. Rev. Lett., 2020, 124, 238002 CrossRef CAS.
  29. K. L. Johnson, Contact Mechanics, Cambridge University Press, 1985 Search PubMed.
  30. E. Bitzek, P. Koskinen, F. Gähler, M. Moseler and P. Gumbsch, Phys. Rev. Lett., 2006, 97, 170201 CrossRef.
  31. C. Maloney and A. Lemaître, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 016118 CrossRef.
  32. H. Mizuno, K. Saitoh and L. E. Silbert, Phys. Rev. E, 2016, 93, 062905 CrossRef.
  33. M. Blank-Burian and A. Heuer, Phys. Rev. E, 2018, 98, 033002 CrossRef CAS.
  34. P. Acharya, S. Sengupta, B. Chakraborty and K. Ramola, Phys. Rev. Lett., 2020, 124, 168004 CrossRef CAS.
  35. M. Fan, K. Zhang, J. Schroers, M. D. Shattuck and C. S. O'Hern, Phys. Rev. E, 2017, 96, 032602 CrossRef.
  36. M. Lundberg, K. Krishan, N. Xu, C. S. O'Hern and M. Dennin, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 041505 CrossRef.
  37. I. Regev, J. Weber, C. Reichhardt, K. A. Dahmen and T. Lookman, Nat. Commun., 2015, 6, 8805 CrossRef CAS.
  38. P. Das, H. A. Vinutha and S. Sastry, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 10203 CrossRef.
  39. 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 CrossRef.
  40. J. R. Royer and P. M. Chaikin, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 49 CrossRef CAS.
  41. K. VanderWerf, W. Jin, M. Shattuck and C. S. O'Hern, Phys. Rev. E, 2018, 97, 012909 CrossRef.
  42. L. E. Silbert, D. Ertas, G. S. Grest, T. C. Halsey and D. Levine, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 65, 031304 CrossRef.

This journal is © The Royal Society of Chemistry 2020