PHYSICAL REVIEW E 108, 034901 (2023)

Designing the pressure-dependent shear modulus using tessellated granular metamaterials
Jerry Zhang ,1,* Dong Wang ,1,* Weiwei Jin ,1 Annie Xia,1 Nidhi Pashine ,1 Rebecca Kramer-Bottiglio,1Mark D. Shattuck,2 and Corey S. O’Hern 1,3,4,5,†
1

Department of Mechanical Engineering & Materials Science, Yale University, New Haven, Connecticut 06520, USA
Benjamin Levich Institute and Physics Department, The City College of New York, New York, New York 10031, USA
3
Department of Physics, Yale University, New Haven, Connecticut 06520, USA
4
Department of Applied Physics, Yale University, New Haven, Connecticut 06520, USA
5
Integrated Graduate Program in Physical and Engineering Biology, Yale University, New Haven, Connecticut 06520, USA
2(Received 17 March 2023; accepted 22 August 2023; published 6 September 2023)
Jammed packings of granular materials display complex mechanical response. For example, the ensembleaveraged shear modulus G increases as a power law in pressure p for static packings of soft spherical particles
that can rearrange during compression. We seek to design granular materials with shear moduli that can either
increase or decrease with pressure without particle rearrangements even in the large-system limit. To do this,we construct tessellated granular metamaterials by joining multiple particle-filled cells together. We focus oncells that contain a small number of bidisperse disks in two dimensions. We first study the mechanical properties
of individual disk-filled cells with three types of boundaries: periodic boundary conditions (PBC), fixed-length
walls (FXW), and flexible walls (FLW). Hypostatic jammed packings are found for cells with FLW, but not incells with PBC and FXW, and they are stabilized by quartic modes of the dynamical matrix. The shear modulus
of a single cell depends linearly on p. We find that the slope of the shear modulus with pressure λc < 0 for all
packings in single cells with PBC where the number of particles per cell N  6. In contrast, single cells withFXW and FLW can possess λc > 0, as well as λc < 0, for N  16. We show that we can force the mechanicalproperties of multicell granular metamaterials to possess those of single cells by constraining the end pointsof the outer walls and enforcing an affine shear response. These studies demonstrate that tessellated granularmetamaterials provide a platform for the design of soft materials with specified mechanical properties.
DOI: 10.1103/PhysRevE.108.034901

I. INTRODUCTION

Granular materials represent an interesting class of physical systems that are composed of individual macroscopic
particles that interact via dissipative, contact forces [1]. As
a result of the dissipative particle interactions, granular materials come to rest in the absence of external driving, suchas applied shear or vibration. Because of this, they frequently
occur in amorphous states lacking long-range positional order.
Further, granular systems can undergo a jamming transition,
where they develop nonzero bulk and shear moduli when they
are compressed to large packing fractions [2–4].
There have been numerous computational [2,5–20] andexperimental [21–31] studies of the structural and mechanicalproperties of jammed granular materials. In particular, it has
been shown that the shear modulus depends sensitively on
structural disorder, nonaffine particle motion, the number ofcontacts, and anisotropy of the interparticle contact network[4,15,16,32–44]. For example, in jammed packings of frictionless spherical particles with purely repulsive linear springinteractions, we have shown that the shear modulus G decreases linearly with pressure p along “geometrical families,”*†

These authors contributed equally to this work.
corey.ohern@yale.edu

2470-0045/2023/108(3)/034901(14)where the network of interparticle contacts does not change
during isotropic compression [20]. If a particle rearrangement
occurs during the compression at p = p∗ , e.g., through theaddition of an interparticle contact, G jumps discontinuouslyat p∗ and the linear relation between G and p no longer
holds. Also, when a particle rearrangement occurs, it is difficult to predict the new interparticle contact network and themechanical properties of the jammed packing are no longer
reversible. The range of pressure p over which the contact
network does not change decreases with increasing systemsize p ∼ N −1 , where N is the number of particles in thesystem. Thus, in the large-N limit, granular packings undergo
frequent irreversible particle rearrangements to new jammed
packings after each p increment. During compression, eachnew contact network typically possesses an increased numberof contacts, and thus the shear modulus increases with pressure. In fact, studies have shown that the ensemble-averaged
shear modulus scales as G ∼ p0.5 in the large-pN 2 limit for
jammed packings of spherical particles with purely repulsive
linear spring interactions [15].
In this article, we design granular metamaterials for which
the shear modulus can either decrease or increase with increasing pressure with no particle rearrangements. In linear
elastic solids, the shear modulus does not depend on thepressure. In conventional atomic and molecular solids, boththe bulk and shear moduli increase with pressure at large

034901-1©2023 American Physical Society

JERRY ZHANG et al.

PHYSICAL REVIEW E 108, 034901 (2023)

FIG. 1. Illustration of a tessellated granular metamaterial, madeup of 36 individual cells. Each cell contains the same jammed
bidisperse packing of N = 4 disks that are confined by four freely
jointed, flexible walls. The interior cells share all four walls and theedge cells share three walls. To generate the collection of disk-filledcells, we first create a disk packing within a single cell, connect
multiple copies of this disk-filled cell, fix the outer blue vertices, andthen allow the disks and interior red vertices to relax during energy
minimization. The variation in the disk shading between different
cells indicates the types of cells based on their adjacent cells. For
example, cells on the right edge of the tessellation only have threeadjacent cells. The cell type is determined by its distance from thefour outer walls and four corner vertices.

pressures [45–47]. Similarly, in conventional granular materials, the shear modulus increases with pressure due to theformation of new contacts during compression, but it is history
dependent. Our design of granular metamaterials will leverage
the recent findings that particle rearrangements in granular
packings with small N are rare and the shear modulus depends linearly on pressure in the absence of rearrangements[20]. Preventing particle rearrangements ensures reversibility of the packing’s mechanical properties and improves ourability to predict them. We will first consider systems in two
dimensions (2D), but these concepts can easily be extended
to three dimensions (3D). For example, we have shown thatthe pressure-dependent shear modulus for jammed packingsof spherical particles is qualitatively the same in 2D and 3D[20].
We envision tessellated granular metamaterials that are
made up of many individual cells that each contain a small
number of grains, i.e., N < 16, and are bounded by four freely
jointed elastic walls. The disks within each cell are jammedwith typically an isostatic number of contacts (see Fig. 1). The
mechanical response of each cell is highly anisotropic, i.e., its
shear modulus depends on the angle θ of the applied shear
relative to the orientation of the confining walls. We find thatthe shear modulus of each cell obeys Gc = Gc0 + λc p, where
Gc = Gc0 at p = 0, and we determine the sign and magnitudeof λ as a function θ , N, and the ratio of the particle andwall stiffnesses. We vary the size of the tessellated granularmetamaterials by adding multiple copies of individual cells
together, e.g., by generating an n × n array of cells that sharethe confining walls. We identify the regimes where the shear
modulus of the full system is similar to that for the individual
cells. In particular, we find that large tessellated granularmetamaterials can possess shear moduli that decrease withincreasing pressure, i.e., the opposite behavior compared toconventional granular materials, and that these materials retain the anisotropy of the individual cells.
The remainder of the article is organized as follows. In
Sec. II, we describe the computational methods, including
the particle-particle, particle-wall, and wall-wall potential energies, the protocols for generating disk-filled single cells(henceforth referred to as “cells”) and collections of multiple cells, and the methods for calculating the pressure, shear
stress, and shear modulus of these structures. In Sec. III, we
present the results on how the boundary conditions, individual
disk packing configuration, and the ratio of the particle towall stiffness affect the relation between the shear modulus
and pressure in single cells, as well as coupled systems composed of Nc = n2 cells. In Sec. IV, we provide conclusions
and discuss promising directions of future research, such as
the mechanical response of tessellated granular metamaterials
in three dimensions. We also include three Appendixes. In
Appendix A, we show that Maxwell-type counting arguments
can be used to determine the minimum number of particleparticle and particle-wall contacts in jammed disk packings
within single cells with fixed length and flexible walls, andexplain the occurrence of “quartic modes” in cells with flexible walls. In Appendix B, we determine analytical expressions
for the dependence of the components of the stiffness matrix
on the angle of the applied simple shear strain for jammed
disk packings in single cells. In Appendix C, we verify thatthe pressure dependence of the single-cell shear modulus isrelated to the second derivative of the packing fraction at jamming onset φJ with respect to shear strain γ for an example
cell with fixed-length walls.
II. METHODS

We study individual cells containing jammed packings of
N bidisperse soft, frictionless disks: N/2 small and N/2 large
disks with diameter ratio σl /σs = 1.4. A diameter ratio ofσl /σs = 1.4 gives rise to disordered jammed disk packings[3,48]. We consider three types of boundary conditions for thecells as illustrated in Fig. 2: (a) periodic boundary conditions(PBC) in square cells with side length L0 , (b) cells with four
straight walls of fixed length L0 (FXW), and (c) cells with four
flexible walls (FLW) such that adjacent vertices are connectedby linear springs with preferred length L0 . For boundary condition (c), the connected walls are freely jointed such that theangle between them can change without energy cost.
We model the tessellated granular metamaterials using linear spring interaction potentials (either purely repulsive or
double sided), which are commonly used in discrete element
method simulations of granular materials [3]. Frictionlessdisks that interact via pairwise, purely repulsive linear springforces are placed within each cell. The corresponding interparticle potential energy is given bypp 2 
pp 
r jk
r jk
   pp
pp pp
1−
U r jk =
,
(1)
1−
2
σ jk
σ jk

034901-2

DESIGNING THE PRESSURE-DEPENDENT SHEAR …
(a)

PHYSICAL REVIEW E 108, 034901 (2023)

spring potential energy

(b)


bb 2
 bb  bb
ri,i+1
1−,
U ri,i+1 =
2
L0
bb

L0

(c)

L0

where bb is the characteristic energy scale of the linear springbb
is the distance between end points i and i + 1,
potential, ri,i+1
and L0 is equilibrium length for the ith wall. The force on
end point i from end point i + 1 in the ith wall is fibb =bbbbbb)r̂i,i+1, where r̂i,i+1
is the unit vector pointing−(dUibb /dri,i+1
from end point i + 1 to i.
We calculate the stress tensor αβ (with α, β = x, y) of
the tessellated granular metamaterials using the virial expression [50]. For cells with PBC, the total potential energypp) and stresses are generated only fromis U = Nj<k U pp (r jkpp,interparticle forces. The total stress tensor is thus αβ = αβ
where

(d)

L0

y
L0

L0

(3)

x
θ

N

1
pp
αβ =

L0

FIG. 2. Illustration of cells that contain N = 6 bidisperse diskswith three different boundary conditions: (a) periodic boundary conditions (PBC) in a square cell with side length L0 , (b) a cell withfour straight walls of fixed length L0 (FXW), and (c) a cell with four
flexible walls (FLW) such that adjacent vertices are connected bylinear springs with preferred length L0 . To generate jammed disk
packings within each cell, we successively compress the system,
fixing the blue vertices after the compression, and then allowing thedisks and red vertices to relax. (d) Illustration of the application of
simple shear strain γ = 0.2 to an originally square cell (solid line)
at angle θ relative to the x axis, which generates the parallelogramshaped cell indicated by the dashed-dotted line.
pp
where  pp gives the strength of the repulsive interactions, r jk
is the distance between the centers of disks j and k, σ jk is thesum of the radii of disks j and k, and (·) is the Heaviside
step function. The repulsive force on disk j from k is fjkpp =pp pppp)r̂ jk , where r̂ jk
is a unit vector pointing from−(dU pp /dr jk
the center of disk k to the center of disk j. Previous studies
have shown that the soft particle model in Eq. (1) generates
the same disk packings at jamming onset at those for rigid
disks [49].
For PBC boundary condition (a), there are only interparticle interactions. For boundary conditions (b) (FXW) and(c) (FLW), we also consider repulsive interactions between
the disks and walls using the purely repulsive linear spring
potential energy
2 r jipbr jipb   pbpb pb
U r ji =
1−
,
(2)
1−
2
Rj
Rj

where  pb is the strength of the repulsive interactions between
the disks and walls, r jipb is the shortest distance between thecenter of disk j and the ith wall, and R j is the radius of disk
j. The repulsive force on disk j from the ith wall is fjipb =−(dU pb /dr jipb )r̂ jipb , where r̂ jipb is the unit vector pointing to thecenter of disk j and perpendicular to the ith wall.
For FLW boundary conditions, we consider interactions
between the wall end points using the double-sided linear

pp pp
f jkα
r jkβ ,

A j>k

(4)

pp
is the α component of the force
A is the area of the cell, f jkαpp
on disk j from k, and r jkβ is the β component of the separation
vector from the center of disk k to the center of disk j.
For cells with physical walls as shown in Figs. 2(b) and2(c), the forces between the walls and particles also contributeto the stress tensor. For cells with FXW boundary conpp)+ditions, the total potential energy is U = Nj<k U pp (r jkN  4pb pb
j=1
i=1 U (r ji ). In this case, the total stress tensor is
αβ =

pp
αβ +

pb
αβ , where

1
pb
αβ =

4

N

A i=1 j=1

pb pb
f jiα
r jiβ ,

(5)

pb
f jiα
is the α component of the force on disk j from the ith
pbwall of the cell, and r jiβ
is the β component of the separation
vector from the contact point between wall i and disk j to thecenter of disk j.
For cells with FLW boundary conditions, in addition to the interparticle and particle-wall interactions,
the walls store potential energy. Thus, the total poten pp) + Nj=1 4i=1 U pb (r jipb ) +tial energy is U = Nj<k U pp (r jk
4
pp
bb bb
αβ = αβ +i=1 U (ri,i+1 ). The total stress tensor is
pb
αβ +

bb
αβ for cells with flexible walls, where

1
bb
αβ =

4

A i=1

bb
fiαbb riβ
,

(6)

bb
fiαbb is the α component of the spring force from wall i, and riβ
is the β component of the vector with the same length as wall ipointing in the same direction as fibb . The pressure of the cellis p = ( xx + yy )/2 and the shear stress is = − xy . We
use  pp /σs2 for the units of stress and shear modulus and  pp
for units of energy.
To generate jammed disk packings within a single cell, we
first place N disks randomly in the cell at a dilute packing
fraction φ < 10−3 . We then apply an affine isotropic compressive strain to the disk positions and decrease the length of the034901-3

JERRY ZHANG et al.

PHYSICAL REVIEW E 108, 034901 (2023)

walls by L to achieve a small packing fraction increment
φ/φ = 2L/L0 = 2 × 10−3 , followed by potential energy
minimization using the fast inertia relaxation engine (FIRE)
algorithm [51]. During energy minimization, the disk positions change, while the end points of the fixed-length walls
for the boundary conditions depicted in Figs. 2(a) and 2(b)are held fixed. However, the end points of the flexible walls
in Fig. 2(c) are allowed to move during energy minimization.
After energy minimization, we calculate the pressure p of
the cell. If p is less than the target pressure pt = ptj = 10−7 ,we again compress the system by φ/φ and perform energy minimization. If p > ptj , we return to the previous disk
and wall configuration, compress the system by φ/2, andperform energy minimization. We repeat this process until
the cell pressure satisfies |p − ptj |/ptj < 10−4 . Our results do
not change if we choose smaller values of ptj . For all threeboundary conditions in Figs. 2(a)–2(c), we generate 104 disk
packings at jamming onset.
To investigate the mechanical response of the cells as afunction of pressure, we apply isotropic compression to thecells at jamming onset to achieve a range of pressure values
pct that are logarithmically spaced between 10−7 to 10−2 .
To ensure that the shape of the cells does not significantly
deviate from the cell shape at jamming onset, we fix the end
points of the walls for all three types of boundary conditions
when generating the cells with pressures above jamming onset
ptj . Since we are using energy minimization to generate thepackings, nearly all packings that we obtain at various pct aremechanically stable.
To calculate the shear modulus of a single cell Gc at anangle θ relative to the x axis, we first rotate the cell clockwiseby θ , as shown in Fig. 2(d). Determining Gc (θ ) allows us toassess the anisotropy of the mechanical response of single
cells. We then apply successive small steps of simple shear
strain γ = 5 × 10−9 (where x is the shear direction and yis the shear gradient direction) to the disks and walls witheach strain step followed by potential energy minimization.
Note that after the applied simple shear strain, the walls remain fixed during energy minimization for all three boundary
conditions. We obtain the shear modulus for a single cell by
calculating Gc = d c /dγ , where c is the shear stress of asingle cell. Over the range of shear strain used to measure
Gc , xy (γ ) is linear. In practice, linear fits of xy (γ ) yieldR-squared values that satisfy 1 − R2 < 10−11 .
We build large-scale tessellated granular metamaterials by
joining multiple copies of a given cell with flexible walls
at jamming onset, e.g., the collection of 36 coupled cells in
Fig. 1. After joining the cells, we perform potential energy
minimization with the outermost (blue) wall end points held
fixed, while the internal (red) end points, as well as the disk
positions, are allowed to relax. Disks within a given cell onlyinteract with other disks and the walls of that cell. Interiorwall end points have four connections to other walls, while
exterior wall end points have either two or three connections
to other walls. We generate tessellated granular metamaterials
at jamming onset with p = ptj = 10−7 , compress the systems
to achieve pressures that are logarithmically spaced between
ptj and 10−2 , and measure the shear modulus at each pressure.
The shear modulus G of the collection of cells is calculated
in the same way as that for a single cell. In particular, wefirst rotate the aggregate by θ clockwise, and then we apply
small successive steps of simple shear strain, γ = 5 × 10−9 ,with each step followed by energy minimization, where theouter vertices are held fixed and the inner vertices, as wellas the disks, are allowed to relax. The total shear stressof the tessellated granular metamaterial is the sum of pp
and pb for all cells and the unique contributions to bb
for all of the cell walls. The shear modulus of the tessellated granular metamaterial is given by G = d /dγ . Whencompressing the tessellated metamaterials to increase p, orapplying simple shear to measure G, we impose a specified
displacement field on the particle and wall positions prior toenergy minimization. After imposing the displacement field,
we allow the particles and only the internal wall vertices torelax their positions under potential energy minimization.
After we apply each simple shear strain step followed byenergy minimization to tessellated granular metamaterials,
which can in principle give rise to nonaffine displacements,
we calculate the displacement field F pq of all cell wall end
points to characterize the nonaffine displacements of eachcell. We find the strain field that minimizes the total nonaffine
displacement of all end points for a given cell and simple shear
strain step [52]:F pq =Xps (Y −1 )sq ,(7)0
rci,p rci,s

(8)

0
0
rci,s
rci,q
.

(9)

s=x,y

where
4

Xps =
i=1

and
4

Ysq =
i=1

0
and rci,s are the sth component of the separation
Here, rci,svector from the center of mass of a given cell to its ith end
point before and after the applied simple shear strain, respectively. We subtract the applied simple shear strain γ from Fxy
to determine the nonaffine displacement field.

III. RESULTS

In this section, we describe the results for the mechanicalresponse of single cells, as well as large collections of cells. In
Sec. III A, we enumerate all of the distinct N = 4 bidisperse
disk packings in single cells at jamming onset for all threeboundary conditions. We determine whether the shear modulus for single cells Gc increases or decreases with pressure
over the full range of θ in Sec. III B. We find that Gc forcells with PBC nearly always decreases with pressure (for all
shear angles), while Gc can either decrease or increase withpressure for single cells with (both FXW and FLW) physicalwalls. We further show that the slope of the shear modulus
versus pressure λc = dGc /d p for single cells can be tuned
by varying the particle-wall interaction energy  pb and wall
stiffness bb . Finally, in Sec. III C, we emphasize that the signand magnitude of λc for a single cell can be maintained even ina large collection of cells since the assembly prevents particle
rearrangements. We then show that the mechanical response

034901-4

DESIGNING THE PRESSURE-DEPENDENT SHEAR …
(a)

1

(b)

1

2

2

3

4

3

5

6

PHYSICAL REVIEW E 108, 034901 (2023)
(c)
1

2

4

5

3

6

7

FIG. 3. Distinct N = 4 bidisperse cells at jamming onset for (a) periodic boundary conditions (PBC, dashed lines), (b) fixed-length physicalwalls (FXW), and (c) flexible physical walls (FLW). The solid red lines indicate interparticle contacts. The nongray background colors indicate
shared interparticle contact networks in (b) and (c). The disks shaded in dark gray are “rattler” particles with no interparticle contacts.

of large collections of cells can deviate from the single-cellbehavior when we allow the outer cell walls to relax andchange their positions during energy minimization.
A. Single cells with N = 4
1. Jammed configurations

We first illustrate the different types of jammed bidisperse
disk packings that occur in single cells with PBC and physicalwalls. In Fig. 3, we show all possible jammed cells with N =4. We find three distinct jammed packings for single cells with
PBC [53], six distinct packings for cells with FXW [54], andseven distinct packings for cells with FLW. For N = 2, there
is only one distinct jammed cell with the disks arranged along
the diagonal of the cell.
For jammed cells with FLW, the shape of the boundary
is not typically a square, as shown in Fig. 3(c), since theenergy function for the walls does not include a bending
energy term. Despite this, we show that several of the jammed
configurations in the cells with FXW and FLW share thesame interparticle contact networks, e.g., configuration 1 inFigs. 3(b) and 3(c).
For N = 4, we find that rattler particles occur in jammed
cells with FXW and FLW. See configurations 5 and 6 inFig. 3(b) and configuration 7 in Fig. 3(c). Rattler disks also occur for cells with PBC and physical walls for N > 4. Since our
focus is on jammed packings that do not undergo particle rearrangements during simple shear and isotropic compression,
we will not include calculations of the mechanical response
for cells with rattler disks.
2. Contact number

The boundary conditions of the cells affect the number ofinterparticle contacts that are required at the onset of jamming. The numbers of degrees of freedom for the cells are
the following: PBC, Nd = 2N  − 1, FXW, Nd = 2N  + 1, andFLW, Nd = 2N  + 2, where N  = N − Nr and Nr is the number of rattler disks. (see Appendix A.) For mechanically stabledisk packings [55,56], the number of contacts must satisfy

Nc  Nd . Packings with Nc = Nd are isostatic; packings with
Nc < Nd are hypostatic; and packings with Nc > Nd are hyperstatic. For N = 4, in PBC, we find that the jammed bidisperse
disk packings are either isostatic (Nc = Nd , N  = 4, Nc = 7
for configuration 1) or hyperstatic (Nc > Nd , N  = 4, Nc = 8
for configuration 2, and N  = 4, Nc = 9 for configuration 3).
In cells with FXW, all of the packings are isostatic (N  = 4,
Nc = 9 for configurations 1–4, N  = 3, Nc = 7 for configuration 5, and N  = 2, Nc = 5 for configuration 6). In thecells with FLW, most of the jammed bidisperse disk packings
are hypostatic (Nc < Nd , N  = 4, Nc = 9 for configurations
1–5 and N  = 3, Nc = 7 for configuration 7). As discussed inSec. III A 3, quartic modes are required to stabilize hypostatic
jammed packings. In contrast, configuration 6 in Fig. 3(c) is
hyperstatic (N  = 4, Nc = 13).
3. Quartic modes

Hypostatic jammed packings have only been reported for
packings of nonspherical particles [12,19] and particles withshape and size degrees of freedom [57–59]. Our results indicate that jammed packings of spherical particles can also
be hypostatic in cells with FLW. We have shown in previous studies that jammed hypostatic packings are stabilized
by quartic modes [12,19,58], which do not occur in isostaticand hyperstatic packings. Indeed, we find that hypostatic cells
at jamming onset possess Nd − Nc quartic modes (see Appendix A). For N > 4, we also find that jammed disk packings
are isostatic in cells with PBC, either isostatic or hyperstatic
for cells with FXW, and either isostatic, hyperstatic, or hypostatic for cells with FLW. At large N (N > 16), we findjammed disk packings are typically isostatic in all types ofboundary conditions studied (see Appendix A).
B. Shear modulus versus pressure for a single cell
1. Slope of shear modulus with pressure and anisotropy

In Figs. 4(a)–4(c), we show the shear modulus Gc (θ ) of
single cells as a function of pressure p over the full range

034901-5

JERRY ZHANG et al.

PHYSICAL REVIEW E 108, 034901 (2023)
π/2 (b)(a) 0.01

0.5

Gc (θ)θ

10-3
Gc (θ)+

1
G'c (θ)Gc (θ)

-0.02

0.12

0.02

0
-0.01

(c)

0.01

0.08

0.04

10-7

0.5
p

1

10-5
p

10-4

10-6

10-3
10-3

0

0
(e) 1

0

(f)
Conﬁg 1Conﬁg 2Conﬁg 3Conﬁg 40
λ c (θ)Conﬁg 1Conﬁg 2Conﬁg 3

-10

10-6

10-4
p

p

-5
λ c (θ)

0

10-4

10-6

Conﬁg 1Conﬁg 2Conﬁg 30.5

10-2

Conﬁg 4Conﬁg 5Conﬁg 60
λ c (θ)

(d)

0
0

+

-0.03

-0.5

-1

-15

-1
-2

-20

-1.5
0

0.1

0.2

0.3

0.4

0.5

0

0.1

0.2

0.3
θ/π

θ/π

0.4

0.5

0

0.1

0.2

0.3

0.4

0.5

θ/π

FIG. 4. (a)–(c) Shear modulus Gc (θ ) versus pressure p for a single cell in configuration 1 (PBC) in Fig. 3(a), configuration 4 (FXW)in Fig. 3(b), and configuration 5 (FLW) in Fig. 3(c). The color of each curve indicates the angle θ of the applied shear strain, which varies
from 0 (blue) to π /2 (red). The inset gives Gc (θ ) = [Gc (θ ) − Gc0 (θ )]/λc (θ ) for the data in the main panel of (a). (d)–(f) The slope of Gc (θ )versus p, λc (θ ), plotted as a function of θ /π for the configurations in Fig. 3 (except those with rattler disks). The different colors and symbols
correspond to different configurations of the cells with (d) PBC, (e) FXW, and (f) FLW (blue circles: configuration 1; red crosses: 2; purple
upward triangles: 3; magenta squares: 4; yellow diamonds: 5; and green left triangles: 6). The horizontal dashed lines indicate λc = 0.

of shear angles θ for cells with PBC, FXW, and FLW, respectively. In contrast to the behavior for large-N systems,
we find that the disks do not rearrange and Gc (θ ) varies
continuously with p over more than four orders of magnitude.
For cells with PBC, Gc (θ ) typically decreases with p as shownin Fig. 4(a). In contrast, for cells with FXW [Fig. 4(b)] andFLW [Fig. 4(c)], Gc (θ ) can either decrease or increase with p,depending on the value of θ .
As we showed previously for jammed packings of spherical particles with PBC, we find quite generally that Gc (θ )varies linearly with p [20,38], Gc (θ ) = Gc0 (θ ) + λc (θ )p, forcells with PBC and physical walls in the absence of particle rearrangements [see the inset to Fig. 4(a)]. Gc0 (θ ) gives
the single-cell shear modulus in the zero-pressure limit andλc (θ ) = dGc (θ )/d p gives the slope [20]. In Figs. 4(d)–4(f),we plot λc (θ ) as a function of θ for all N = 4 cells without rattlers. We show that λc (θ ) = λc,a sin[4(θ − θ0 )] + λc,dc varies
sinusoidally with period π /2, where λc,a is the amplitude,
θ0 is the phase shift, and λc,dc is the mean value of λc (θ )[40]. (Previous studies have shown that the shear modulus
of jammed packings of spherical particles is sinusoidal withperiod π /2 [15,44].) λc (θ ) < 0 for nearly all θ values andcells for PBC, except for configuration 2 [Fig. 3(a)] in therange 0.2  θ /π  0.3 [Fig. 4(d)]. For cells with FXW and
FLW, we observe similar sinusoidal behavior for λc (θ ), butthere are large θ ranges where λc (θ ) > 0. Our results showingthat λc,a ∼ λc,dc emphasize that cells at small N are highly
anisotropic. For cells with FLW, we do not find a correlation
between λc (θ ) > 0 and the occurrence of quartic modes as
discussed in Sec. III A 3. Isotropic, linearly elastic materials

in 2D possess only two elastic moduli, i.e., the bulk and shear
moduli. In this study, we focus on the response of the system
to simple shear strain. However, small, jammed disk packings
at low pressures are anisotropic [41,44], and thus they possess
more than two elastic constants that are inter-related and depend on θ . In Appendix B, we derive the θ dependence forall six stiffness matrix elements Ĉ and show that the elementsof Ĉ are related to each other under rotations. In contrast, forisotropic, linearly elastic materials, Ĉ33 ≡ G is independent of
θ . Therefore, in future studies we will consider all elements ofĈ to fully understand the pressure- and angle-dependent elastic moduli for anisotropic materials. In addition, our previous
studies have shown that the sign of λc (θ = 0) is determined by
the second derivative of the packing fraction at jamming onsetwith respect to γ in jammed packings of spherical particles in
PBC [20,38]. In Appendix C, we show that this relation is stilltrue at any θ in cells with fixed-length walls. In particular,
this relationship expresses Ĉ33 in terms of derivatives of thepacking fraction U and p with respect to shear strain for any
given θ and determines which terms in this expression depend
most strongly on p. Similar analyses will be carried out forthe other elements of Ĉ in future studies to gain a complete
understanding of the pressure-dependent elastic moduli ofsmall jammed disk packings.
2. Probability of λc (θ) > 0

As N increases, the probability P+ to obtain a cell withλc (θ ) > 0 decreases rapidly for PBC. As shown in Fig. 5,
we do not find λc (θ ) > 0 for cells with N  6 for PBC. For

034901-6

DESIGNING THE PRESSURE-DEPENDENT SHEAR …

PHYSICAL REVIEW E 108, 034901 (2023)
3. Effects of wall and particle stiffnesses on λc (θ)

1
PBC
FXW
FLW

0.8

0.6

0.4

0.2

0
4

2

6

8
N

32

16

FIG. 5. Probability P+ that cells possess λc (θ ) > 0 for any
nonzero range of θ, as a function of system size N for PBC (bluecircles), FXW (red crosses), and FLW (pink triangles).

We next investigate the dependence of λc (θ ) on theparticle-wall stiffness  pb / pp and wall stiffness bb / pp relative to the strength of the repulsive interparticle interactions
in cells with FXW and FLW. In Fig. 6(a), we show thatfor configuration 4 depicted in Fig. 3(b), λc,a , λc,dc , and θ0
undergo only small variations when  pb / pp changes by nearly
two orders of magnitude. λc,a and λc,dc converge for  pb / pp 
10 for all N = 4 cells with FXW [Figs. 6(b) and 6(c)]. For
configuration 4, we find that λc (θ ) > 0 for a finite range of θin the large  pb /bb limit. We can also fix  pb / pp and showthat λc,a , λc,dc , and θ0 converge in the large bb / pp limit [seeFigs. 6(d)–6(f)]. We find that λc (θ ) < 0 (for all θ ) at large
bb / pp for all N = 4 configurations with FLW since λc,dc < 0and λc,a < |λc,dc |. These results emphasize that cells with
FLW become similar to cells with PBC [with λc (θ ) < 0] in thelarge bb limit. Thus, particle-wall interactions are essential forλc (θ ) > 0.
C. Shear modulus versus pressure for tessellated
granular metamaterials
1. Lock-in of single-cell behavior

cells with physical walls, P+ also decreases with increasingN, but not as rapidly as that for cells with PBC. These results
emphasize that if one wants to tune the pressure dependenceof Gc (θ ) [i.e., between λc (θ ) > 0 and λc (θ ) < 0], one should
employ cells with small N.
(a) 1

λ c,a

λ c,dc0.9-10.8
Conﬁg 1Conﬁg 2Conﬁg 3Conﬁg 4

-1.5

-1

0.1

0.2

0.3

0.4

-2
10-2

0.5

θ/π
(e)

(d) 1
ε bb/ε pp = 0.01
ε bb/ε pp = 0.1
ε bb/ε pp = 1
ε bb/ε pp = 10
ε bb/ε pp = 50

0

0.1

0.2

0.3
θ/π

Conﬁg 1Conﬁg 2Conﬁg 3Conﬁg 4Conﬁg 5Conﬁg 6

0.4

0.5

100
ε bb / ε pp

102

0.6

-0.6

0.4

-0.8

0.2

-1
10-2

102

0.8

-0.4

-1

100
ε pb / ε pp

(f) 1

0
-0.2

0.6
10-2

102

λ c,dc0100
ε pb / ε pp

0.7

λ c,a

0

λ c (θ)1-0.5

λ c (θ)0(c) 1.10(b)
ε pb/ε pp = 0.01
ε pb/ε pp = 0.1
ε pb/ε pp = 1
ε pb/ε pp = 10
ε pb/ε pp = 50

We now study the pressure dependence of the shear modulus G(θ ) for tessellated granular metamaterials (Fig. 1)
constructed from multiple cells with flexible walls  pb / pp =1 and bb / pp = 0.1. In Fig. 7(a), we show G(θ ) versus p for
Nc = 36 cells that each contain configuration 5 from Fig. 3(c).

100
ε bb / ε pp

102

0
10-2

FIG. 6. (a) The slope of the shear modulus λc (θ ) versus pressure for a single cell with FXW [configuration 4 in Fig. 3(b)] versus the angle
θ of the applied shear strain for several values of the particle-wall interaction energy normalized by the strength of the interparticle repulsive
energy  pb / pp indicated by the different colors and symbols. (b) Mean value λc,dc and (c) amplitude λc,a of the slope of the shear modulus
as a function  pb / pp for the four single cells without rattlers and FXW shown in Fig. 3(b). The symbols and colors in (b) and (c) indicate
the specific configurations in Fig. 3(b). (d)–(f) Similar data to that in (a)–(c) except for single cells with FLW. Data for configuration 5 inFig. 3(c) is shown in (d). The symbols and colors in (e) and (f) indicate the specific configurations in Fig. 3(c). The dashed horizontal lines in(a) and (d) indicate λc = 0.
034901-7

PHYSICAL REVIEW E 108, 034901 (2023)()()

JERRY ZHANG et al.

FIG. 7. (a) Shear modulus G(θ ) versus pressure p for a tessellated granular metamaterial composed of Nc = 36 cells containing the sameN = 4 disk configuration [i.e., configuration 5 in Fig. 3(c)] with  pb / pp = 1 and bb / pp = 0.1. Different colors correspond to the angle θ ofthe applied simple shear strain, which varies from 0 (blue) to π /2 (red). (b) λ(θ ) for tessellated granular metamaterials with  pb / pp = 1 andbb / pp = 0.1 and different values of Nc (blue circles: 4, red crosses: 16, purple upward triangles: 36, and magenta squares: 64) and a single
cell with bb/ pp = 0.05 (yellow diamonds) and 0.1 (green left triangles). All cells contain the same disk configuration 5 in Fig. 3(c) with FLW.
(c) The off-diagonal term of the fitted displacement matrix Fxy − γ in Eq. (7) as a function of the applied simple shear strain γ at p = 10−7 foreach cell (indicated by different colors) for the tessellated granular metamaterial in (a).

Similar to the results for Gc (θ ) for single cells, the mechanicalresponse of tessellated granular metamaterials shows strong
shear angle dependence. In particular, for some values of θ ,the slope of G(θ ) versus p, λ(θ ) > 0, and for other values,λ(θ ) < 0. In Fig. 7(b), we show that λ(θ ) possesses weaksystem-size dependence as Nc is increased. λ(θ ) for the multicell system in the large-Nc limit converges to λc (θ ) for asingle cell with FLW with bb / pp = 0.05, which is half of thevalue for the multicell system. This result can be explained
because each wall in the tessellated granular metamaterial isshared by its neighboring cell except for those on the exterior.
Since λ(θ ) for the tessellated granular metamaterial mimics
that for single cells, these results indicate that we can lockin the behavior of the shear modulus versus pressure for asingle cell in tessellated granular metamaterials in large-Nc
limit. In particular, we find finite regions of θ where λ(θ ) > 0
in the large-Nc limit without particle rearrangements. The
similarity of the mechanical response between the multicell
and single-cell systems is caused by the fact that all of the single cells display similar displacements during applied simple
shear strains. In Fig. 7(c), we show that the nonaffine displacements of each cell, caused by energy minimization after
each applied simple shear strain, is negligible. The nonaffine
motion is obtained by measuring the difference in the offdiagonal element of the strain tensor of the interior vertices
and the expected value from simple shear, Fxy − γ , for eachcell [52]. Previous studies have shown that nonaffine particle
motion strongly affects the magnitude of the shear modulus[34,37,60]. Thus, near zero nonaffine displacements of thecells indicate that the strain of the tessellated metamatarials
has locked in the strain of each cell. In general, the results
for the tessellated granular metamaterials are qualitatively thesame for any jammed packings that are used to construct thetessellated granular metmaterial.
2. Constraints on outer wall vertices

We also investigate how many constraints on the outer wall
end points are necessary to enforce the mechanical responseof the single cells in tessellated granular metamaterials. To

address this question, we allow different wall end points on theouter boundary to relax during energy minimization following
the applied isotropic compression and simple shear strain [seeFig. 8(a)]. We find that G(θ ) versus p for tessellated granularmetamaterials with even a single mobile outer wall end point
deviates from that of the tessellated granular metamaterial
where all outer wall end points are fixed. We show G(0) versus
p as an example in Fig. 8(b). Allowing the outer wall end
points to move gives rise to buckling of tessellated granularmetamaterials caused by compression, as shown in Figs. 8(c)and 8(d). This buckling induces changes in the shape of single
cells compared to those of single cells with fixed outer wall
end points during compression and shear, which causes thedeviations in G(θ ) versus p. Therefore, to ensure that tessellated granular metamaterials lock in single-cell behavior, it isnecessary to constrain all of the outer wall end points.
IV. CONCLUSIONS AND FUTURE DIRECTIONSIn the large-system limit, the shear modulus G of static
packings of spherical particles increases with pressure due tofrequent particle rearrangements and nonaffine particle motion that enable the packings to increase their contact number
with increasing pressure. In this work, we investigate a classof granular materials, tessellated granular metamaterials, thatallow us to control the slope of the shear modulus versus
pressure by preventing particle rearrangements even in thelarge-system limit. We focus on tessellated granular metamaterials in two dimensions, which are collections of Nc coupledcells that each contain N bidisperse disks enclosed by fourphysical walls. In particular, we can design tessellated granular metamaterials with negative slope of the shear modulus
with pressure even in the large-Nc limit.
We first studied the mechanical properties of single cellswith three sets of boundary conditions: PBC, FXW, andFLW. Packings with small N do not undergo frequent particle
rearrangements, and thus we enumerated all possible mechanically stable cells with all three boundary conditions for N 
8. We find that the mechanically stable cells with PBC andFXW are either isostatic or hyperstatic, while those with FLW

034901-8

DESIGNING THE PRESSURE-DEPENDENT SHEAR …

PHYSICAL REVIEW E 108, 034901 (2023)

3

(a)

(b)

+

+
G(θ=0)

8

2

+

+

4

10-3
Fixed
1
2
3
4
1&3
2&4
1, 2, 3, & 4

6

4

2

0

0.002

0.004

0.006

0.008

0.01

p

+
1

(c)

(d)

FIG. 8. (a) Tessellated granular metamaterial with Nc = 36 described in Fig. 7, except now some of the wall end points marked by red
crosses are no longer fixed after the applied isotropic compression and simple shear strain. (b) Shear modulus G(0) measured at θ = 0 as afunction of pressure p when different wall end points on the outer boundary are switched from fixed to mobile (blue circles: all outer wall end
points are fixed, red crosses: end point 1 is mobile, magenta upper triangles: end point 2 is mobile, black squares: end point 3 is mobile, yellow
diamonds: end point 4 is mobile, green left triangles: end points 1 and 3 are mobile, cyan crosses: end points 2 and 4 are mobile, and purple
asterisks: end points 1–4 are mobile). (c), (d) Tessellated granular metamaterials at p = 0.01 when (c) end point 1 is mobile, and (e) end points
1–4 are mobile.

can be hypostatic, as well as isostatic and hyperstatic. The
hypostatic cells with FLW are stabilized by quartic modes, as
found for hypostatic packings of nonspherical and deformableparticles. Second, we showed that the angle-dependent shear
modulus of single cells depends linearly on pressure, Gc (θ ) =λc (θ )p + Gc0 (θ ). Further, the slope of the shear modulus
versus pressure for single cells is strongly anisotropic, i.e.,λc (θ ) = λc,dc + λc,a sin[4(θ − θ0 )] and λc,a ∼ λc,dc . We findthat λc (θ ) < 0 for single cells in PBC with N > 4. In contrast, cells with FXW and FLW and small N can possess
either λc (θ ) > 0 or λc (θ ) < 0. However, the probability ofobtaining cells with λc (θ ) > 0 vanishes in the large-N limit.
These results are summarized in Table I. After studying themechanical response of single cells, we investigated the shear
modulus of tessellated granular metamaterials formed by connecting many single cells with flexible walls. We showed that
we can lock in the mechanical response of single cells intessellated granular metamaterials. The ability to lock in themechanical response of single cells in multicell systems is
reduced if the outer wall end points are free to move during energy minimization after applied deformations. These
results demonstrate that we can build large-scale granular

TABLE I. Summary of the results for the number of contacts Nc and derivative of the shear modulus with respect to pressure λc of single
cells for all three types of boundary conditions: PBC, FXW, and FLW.

Number of contacts at jamming onset Nc
Sign of λc

PBC

FXW

FLW

Isotropic or hyperstaticλc < 0 for Nc > 4

Isotropic or hyperstaitc
λc < 0 or λc > 0

Isotropic, hyperstatic, or hypostaticλc < 0 or λc > 0034901-9

JERRY ZHANG et al.

PHYSICAL REVIEW E 108, 034901 (2023)

ACKNOWLEDGMENTS

The authors acknowledge support from ONR Grant No.
N00014-22-1-2652 and NSF Grant No. DMREF-2118988.
This work was also supported by the High Performance
Computing facilities operated by Yale’s Center for Research
Computing.

APPENDIX A: ISOSTATICITY AND QUARTIC
MODES IN A SINGLE CELL

In this Appendix, we discuss the number of degrees of
freedom Nd in cells for all three types of boundary conditions
and quartic modes in hypostatic disk packings with FLW.

1Hyperstatic, PBC
Hypostatic, PBC
Hyperstatic, FXW
Hypostatic, FXWHyperstatic, FLW
Hypostatic, FLW

0.8

(N)

metamaterials whose mechanical properties do not change
after repeated cycles of compression and decompression, as
well as positive and negative simple shear strain, since particle
rearrangements are eliminated.
These findings raise many interesting directions for future research. First, we found that the mechanical responseof both single- and multiple-cell granular systems is highly
anisotropic. To fully understand the pressure-dependent
mechanical properties of anisotropic materials in two dimensions, we must characterize all six stiffness matrix
components as a function of pressure, which requires additional mechanical tests in addition to simple shear as afunction of θ . Doing so will also aid in the developmentof a continuum elasticity theory for tessellated granular materials. Second, for the current studies, we fixed all of theouter wall end points during energy minimization to enforce
nearly affine simple shear of tessellated granular metamaterials. However, when the outer wall end points are not fixed,
the individual cells can change their shape during energy
minimization that follows the applied compression and simple
shear strain. Thus, it will be interesting to study and predictthe pressure dependence of Gc (θ ) of single cells when theouter wall end points are free to move or bending energy is
included between adjacent end points to generate cells with
arbitrary shapes. Third, we have focused on tessellated granular metamaterials composed of identical single cells. In futurestudies, we will consider tessellated granular metamaterials
composed of single cells with different disk configurations
and boundaries with varied  pb and bb to understand how
the mechanical properties of single cells determine the mechanical properties of the multicell system. Fourth, we have
only considered square tessellated metamaterials. We can also
investigate tessellated metamaterials with different numbersof rows and columns of cells, which can provide more tunability of the pressure-dependent mechanical response. Finally,
we will extend our studies of tessellated granular metamaterials to three dimensions. In three dimensions, there are threeprincipal simple shear directions, instead of one in two dimensions, which provides additional ways to design tessellated
granular metamaterials. For example, we can create strongly
anisotropic tessellated granular metamaterials by having somecells possess λc (θ ) > 0 in one shear direction, others possess λc (θ ) < 0 in another shear direction, and others possessλc (θ ) > 0 in the remaining shear direction.

0.6

0.4

0.2

0
2

4

6

8
N

16

32

FIG. 9. Probability P (N ) (over 104 single cells) of obtaining agiven number of contacts, hyperstatic (Nc > Nd ) or hypostatic (Nc <Nd ), for single cells with PBC, FXW, or FLW as a function of systemsize N.

In PBC, Nd = 2N  − 2 + 1 = 2N  − 1, where N  = N −Nr and Nr is the number of rattler disks. In this expression,the −2 comes from the two global translational degrees of
freedom in periodic boundary conditions; the +1 comes fromthe size degree of freedom of the disks while maintaining
fixed diameter ratio. The +1 can also be interpreted as astate of self-stress [56]. For the degrees of freedom in cellswith FXW and FLW, we must include the degrees of freedomof the wall end points, as well as the disks: Nd = 2N  +2Nv − NB − 3 + 1. Here, Nv = 4 is the number of wall end
points, NB is the number of constraints associated with thewalls, the −3 comes from the two rigid-body translational andone rotational degree of freedom, and the +1 again comesfrom the size degree of freedom of the disks. For FLW, four
springs connect the wall end points, and hence NB = 4. For
FXW, in addition to the length constraint for each wall, theangle between any two neighboring walls is also fixed, andthus NB = 5. Hence, Nd = 2N  + 1 for cells with FXW and
Nd = 2N  + 2 for FLW. For isostatic packings, the total number of contacts satisfies Nc = Nd . A packing is hyperstaticwhen Nc > Nd and hypostatic when Nc < Nd . We show thatthe probability of obtaining a hyperstatic or hypostatic celldecreases with increasing N for all three boundary conditions
as shown in Fig. 9. We note that there are a finite number ofhypostatic packings in cells with FLW even at N  16, whichhighlights the effect of soft physical walls on the structuraland mechanical properties of jammed granular materials.
“Quartic” modes stabilize hypostatic jammed packings.
They can be illustrated through the expansion of the total
potential energy of a jammed disk packing U in terms of small
displacements around its local energy minimum U (R0 ), whereR0 denotes the equilibrium positions of all particles and wallvertices. If we perturb the packing to a new set of positions
R = R0 + δ û, where δ is the amplitude of the perturbation andû is a unit vector characterizing the perturbation direction, U

034901-10

DESIGNING THE PRESSURE-DEPENDENT SHEAR …

PHYSICAL REVIEW E 108, 034901 (2023)

can be approximated by a Taylor expansion:

(a)
100

∂U
∂ U
1
ui + δ 2
ui u j∂Ri
2 i, j ∂Ri ∂R j
2

 = U (R0 ) + δ
U (R)
i

10-5
λM
i

∂ 2U
1
+ δ3
ui u j uk
6 i, j,k ∂Ri ∂R j ∂Rk

10-10

∂ 2U
1
+ δ4
ui u j uk ul + · · · , (A1)
24 i, j,k,l ∂Ri ∂R j ∂Rk ∂Rl

APPENDIX B: STIFFNESS MATRIX AFTER ROTATIONIn this Appendix, we show the angular dependence of
all elements of the stiffness matrix Ĉ, which relates stressand strain. At a predefined orientation with θ = 0, we cancalculate the stiffness matrix Ĉ(0). After rotating the configuration by an angle θ clockwise, the stiffness matrix becomesĈ(θ ) = R̂T (θ )Ĉ(0)R̂(θ ), where⎛⎞cos2 θ
sin2 θ
− 21 sin 2θ
⎜⎟1
2
R̂(θ ) = ⎜(B1)
cos2 θ
sin 2θ ⎟⎝ sin θ⎠.
2
sin 2θ − sin 2θ
cos 2θ

5

10-8

10-6

10
Mode Index i

15

(b)
10-5
10-10
10-15
10-20
10-4

10-2

δ
FIG. 10. (a) Sorted eigenvalues (from smallest to largest) of
the dynamical matrix λM for configuration 5 with FLW boundary
conditions in Fig. 3(c). The first three eigenvalues (corresponding
to two rigid-body translations and one rigid-body rotation) are lessthan 10−15 and decrease in magnitude as we improve the force
balance condition. The fourth eigenvalue (i = 4) corresponds to a“quartic mode” and those with i > 4 correspond to quadratic modes.
(b) Change in potential energy U in response to perturbations along
eigenvectors λ̂M plotted as a function of the perturbation amplitude δ.
The dashed and dotted-dashed lines have slope 2 and 4, respectively.
The colors of the symbols in (b) correspond to perturbations along
the eigenvectors associated with eigenvalues with the same coloredsymbols in (a).

Ĉ12 (θ ) = 41 (Ĉ11 (0) + Ĉ22 (0)) − Ĉ33 (0) sin2 (2θ )+ Ĉ12 (0)(sin4 θ + cos4 θ )+ 21 (Ĉ23 (0) − Ĉ13 (0)) sin(4θ ),+ 21 Ĉ33 (0) sin(4θ ) + 41 Ĉ12 (0) sin(4θ )+ Ĉ13 (0) cos2 θ (2 cos(2θ ) − 1)
+ Ĉ23 (0) sin2 θ (2 cos(2θ ) + 1),

Ĉ11 (θ ) = Ĉ11 (0) cos4 θ + Ĉ22 (0) sin4 θ+ Ĉ33 (0) sin2 (2θ ) + 21 Ĉ12 (0) sin2 (2θ )(B3)Ĉ13 (θ ) = − 21 Ĉ11 (0) sin(2θ ) cos2 θ + 21 Ĉ22 (0) sin(2θ ) sin2 θUsing Eq. (B1), we find the following angle-dependent stiffness matrix elements :(B4)

Ĉ22 (θ ) = Ĉ11 (0) sin4 θ + Ĉ22 (0) cos4 θ + Ĉ33 (0) sin2 (2θ )+ 21 Ĉ12 (0) sin2 (2θ ) − 2Ĉ13 (0) sin(2θ ) sin2 θ+ 2Ĉ13 (0) sin(2θ ) cos2 θ+ 2Ĉ23 (0) sin(2θ ) sin2 θ ,0

ΔU

where the derivatives are evaluated at R0 . The second term on
the right-hand side of this equation is zero since the jammed
packing is in force balance. The third term on the right2
U, whosehand side includes the dynamical matrix M = ∂R∂ i ∂R
jM
eigenvalues λ give the energy associated with perturbations
along the associated eigenvectors λ̂M . In isostatic and hyperstatic packings with Nc  Nd , all nontrivial eigenvalues λM
 − U (R0 ) increases quadraticallyare nonzero and U = U (R)with amplitude δ for perturbations along λ̂M : U = 21 λM δ 2 .
In hypostatic packings, some of the eigenvalues of M aremuch smaller than the nontrivial quadratic eigenvalues, yetthey are larger than the trivial eigenvalues corresponding torigid-body translations and rotations. [For example, see thefourth eigenvalue of M in Fig. 10(a).] The number of thesemodes matches the number of missing contacts, Nd − Nc .
When perturbing along eigenvectors λ̂M associated with oneof these small nontrivial eigenvalues, U first scales quadratically with δ and then scales quartically with δ at larger valuesof δ as shown in Fig. 10(b). In previous studies [61], we haveshown that the characteristic δc at which the scaling of Uchanges from quadratic to quartic decreases with pressure, andthus at zero pressure, U scales quartically with δ. (Theseprevious results emphasize that the fourth term in the Taylor
expansion containing third derivatives of U is small.) Since
these modes scale quartically with δ at zero pressure, we label
them as “quartic modes.”

10-15

(B2)
034901-11− 2Ĉ23 (0) sin(2θ ) cos2 θ ,(B5)JERRY ZHANG et al.

PHYSICAL REVIEW E 108, 034901 (2023)

Ĉ23 (θ ) = − 21 Ĉ11 (0) sin(2θ ) sin2 θ + 21 Ĉ22 (0) sin(2θ ) cos2 θ− 21 Ĉ33 (0) sin(4θ ) − 41 Ĉ12 (0) sin(4θ )+ Ĉ23 (0) cos2 θ (2 cos(2θ ) − 1),(B6)Ĉ33 (θ ) = 14 Ĉ11 (0) + Ĉ22 (0) − 2Ĉ12 (0) sin2 (2θ )+ 21 Ĉ23 (0) − Ĉ13 (0) sin (4θ )+ Ĉ33 (0) cos2 (2θ ).

( )+ Ĉ13 (0) sin2 θ (2 cos(2θ ) + 1)

(B7)

Equations (B2)–(B7) show that generally all six elements of
the reference stiffness matrix contribute to each Ĉ element
at a given angle θ . Therefore, in anisotropic materials, it isimportant to track all Ĉ elements to fully characterize theirmechanical properties.
APPENDIX C: RELATION BETWEEN SHEAR MODULUS
AND MIXED SHEAR STRAIN DERIVATIVESIn this Appendix, we verify that the pressure dependenceof the single-cell shear modulus is related to the variation of
the packing fraction at jamming onset φJ with simple shear
strain γ as shown in previous studies [38]. We illustrate thisrelationship using a single cell containing the N = 2 monodisperse disk packing with FXW in the inset to Fig. 11(a) sinceφJ (γ ) can be calculated analytically for this case. The shear
modulus can be written in terms of three mixed derivatives ofthe simple shear strain:  dU  
1 d dγ p
d
=
G=
dγ φ
A
dγ
φ


 dφ   dφ
p d dγ p
1 dp
−
−
.
(C1)
φ
dγ
φ dγ φ dγ p
φ

In Fig. 11(a), we demonstrate that Eq. (C1) still holds for
single cells with FXW.
The second term in Eq. (C1) is proportional to p, withdφ[d ( dγ) p /dγ ]φ > 0 for single cells with N  6 and PBC [20].
In contrast, the first and third terms in Eq. (C1) do not possess strong p dependence. Therefore, the second derivativeof φJ (γ ) typically determines whether G will increase or
decrease with p. In particular, if φJ (γ ) is concave upward,
G decreases with increasing p, and vice versa.
For single cells containing N = 2 monodisperse disks
at jamming onset with FXW boundary conditions, we cananalytically determine the packing fraction at jamming
onset φJ as a function of the simple shear strain γ andshear angle θ . Consider a square box (with initial side
lengths L0 = 1) whose vertices are located at (0, 0), (1, 0),(1, 1), and (0, 1). After rotation by θ clockwise about theorigin, these four points transform into (0, 0), (cos θ , sin θ ),(cos θ + sin θ , − sin θ + cos θ ), and (sin θ , cos θ ). After
we apply a simple shear strain γ to the rotated box,
the four points become (0, 0), (cos θ − γ sin θ , − sin θ ),( cos θ + sin θ + γ (− sin θ + cos θ ), − sin θ + cos θ ),and(sin θ + γ cos θ , cos θ ). The resulting box is a parallelogram

FIG. 11. (a) Shear modulus Gc (θ ) for a single cell containing amonodisperse N = 2 disk packing with FXW at shear angle θ = 0(red squares and red solid line) and π /4 (blue circles and blue solidline) as a function of pressure p. The squares and circles show Gc (θ )calculated using Gc = d c /dγ , while the solid lines show Gc (θ )calculated using Eq. (C1). The inset shows the N = 2 cell at simple
shear strain γ = 0 and θ = 0. (b) Packing fraction φJ of the single
cell in (a) at jamming onset as a function of γ at several values ofθ as indicated by the different colors and symbols. The symbols andsolid lines correspond to the results from the numerical simulations
and analytical calculations using Eqs. (C2)–(C9), respectively. The
inset provides an illustration of a monodisperse N = 2 packing ina cell with FXW boundary conditions at jamming onset. The disk
centers are specified by (x1 , y1 ) and (x2 , y2 ) and the lengths of theparallelogram side walls are a and b with angle 2α between them.

with unit area and side lengths,
a = 1 − 2γ sin θ cos θ + γ 2 sin2 θ ,(C2)b = 1 + 2γ sin θ cos θ + γ 2 cos2 θ .

(C3)To obtain φJ (θ , γ ), we need to find the radius of the circle r
such that the two disks fit within the parallelogram and aremechanically stable. [See Fig. 11(b) inset, where we haverotated the parallelogram by θ counterclockwise to simplify
the calculation.] The centers of the two disks are located at

034901-12

DESIGNING THE PRESSURE-DEPENDENT SHEAR …

PHYSICAL REVIEW E 108, 034901 (2023)(x1 , y1 ) and (x2 , y2 ), where x1 = r cot α, y1 = r, x2 = a +b cos 2α − r cot α, and y2 = b sin 2α − r cot α, with

The two solutions for r are(C7)√
B2 − 4AC.
(C8)
r1,2 =
2A
The first solution results in circles that are outside of theparallelogram. Thus, r = r2 and√(−B − B2 − 4AC)2
φJ (θ , γ ) = π
.
(C9)
2A2
We verify in Fig. 11(b) that φJ (θ , γ ) determined by the numerical simulations matches that predicted by Eq. (C9). At
γ = 0, which is where Gc (θ ) is measured throughout the main
text, we find that the φJ (γ ) is concave downward at θ = 0
and concave upward at θ = π /4 [Fig. 11(b)]. Thus, sinceλc (θ ) switches sign, we expect a saddle point to occur in theφJ (γ , θ ) plane between θ = 0 and π /4. In future studies, wewill apply a similar approach in Eq. (C1) to obtain the pressure
dependence of all elements of the stiffness matrix.

[1] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, Granular solids,
liquids, and gases, Rev. Mod. Phys. 68, 1259 (1996).
[2] C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Jammingat zero temperature and zero applied stress: The epitome ofdisorder, Phys. Rev. E 68, 011306 (2003).
[3] A. J. Liu and S. R. Nagel, The jamming transition and themarginally jammed solid, Annu. Rev. Condens. Matter Phys.
1, 347 (2010).
[4] M. van Hecke, Jamming of soft particles: Geometry, mechanics,
scaling and isostaticity, J. Phys.: Condens. Matter 22, 033101(2010).
[5] S. Edwards and R. Oakeshott, Theory of powders, Physica A(Amsterdam) 157, 1080 (1989).
[6] L. E. Silbert, A. J. Liu, and S. R. Nagel, Vibrations and Diverging Length Scales Near the Unjamming Transition, Phys. Rev.
Lett. 95, 098301 (2005).
[7] P. Olsson and S. Teitel, Critical Scaling of Shear Viscosity at the Jamming Transition, Phys. Rev. Lett. 99, 178001(2007).
[8] S. Henkes, C. S. O’Hern, and B. Chakraborty, Entropy andTemperature of a Static Granular Assembly: An Ab Initio Approach, Phys. Rev. Lett. 99, 038002 (2007).
[9] C. Song, P. Wang, and H. A. Makse, A phase diagram for
jammed matter, Nature (London) 453, 629 (2008).
[10] T. Aste, Volume Fluctuations and Geometrical Constraints in
Granular Packs, Phys. Rev. Lett. 96, 018002 (2006).
[11] S. Henkes and B. Chakraborty, Statistical mechanics frameworkfor static granular matter, Phys. Rev. E 79, 061301 (2009).
[12] M. Mailman, C. F. Schreck, C. S. O’Hern, and B. Chakraborty,
Jamming in Systems Composed of Frictionless Ellipse-Shaped
Particles, Phys. Rev. Lett. 102, 255501 (2009).
[13] C. F. Schreck, N. Xu, and C. S. O’Hern, A comparison of
jamming behavior in systems composed of dimer- and ellipseshaped particles, Soft Matter 6, 2960 (2010).
[14] M. L. Manning and A. J. Liu, Vibrational Modes Identify Soft
Spots in a Sheared Disordered Packing, Phys. Rev. Lett. 107,108302 (2011).
[15] C. P. Goodrich, S. Dagois-Bohy, B. P. Tighe, M. van Hecke,
A. J. Liu, and S. R. Nagel, Jamming in finite systems: Stability,

anisotropy, fluctuations, and scaling, Phys. Rev. E 90, 022138(2014).
[16] C. P. Goodrich, A. J. Liu, and J. P. Sethna, Scaling ansatz forthe jamming transition, Proc. Natl. Acad. Sci. USA 113, 9745
(2016).
[17] S. Patinet, D. Vandembroucq, and M. L. Falk, Connecting Local
Yield Stresses with Plastic Activity in Amorphous Solids, Phys.
Rev. Lett. 117, 045501 (2016).
[18] S. Martiniani, K. J. Schrenk, K. Ramola, B. Chakraborty, and D.
Frenkel, Numerical test of the edwards conjecture shows that all
packings are equally probable at jamming, Nat. Phys. 13, 848(2017).
[19] K. VanderWerf, W. Jin, M. D. Shattuck, and C. S. O’Hern, Hypostatic jammed packings of frictionless nonspherical particles,
Phys. Rev. E 97, 012909 (2018).
[20] K. VanderWerf, A. Boromand, M. D. Shattuck, and C. S.
O’Hern, Pressure Dependent Shear Response of Jammed Packings of Frictionless Spherical Particles, Phys. Rev. Lett. 124,
038004 (2020).
[21] A. Donev, I. Cisse, D. Sachs, E. A. Variano, F. H. Stillinger, R.
Connelly, S. Torquato, and P. M. Chaikin, Improving the density of jammed disordered packings using ellipsoids, Science
303, 990 (2004).
[22] E. I. Corwin, H. M. Jaeger, and S. R. Nagel, Structural signatureof jamming in granular media, Nature (London) 435, 1075
(2005).
[23] O. Dauchot, G. Marty, and G. Biroli, Dynamical Heterogeneity
Close to the Jamming Transition in a Sheared Granular Material, Phys. Rev. Lett. 95, 265701 (2005).
[24] T. S. Majmudar, M. Sperl, S. Luding, and R. P. Behringer,
Jamming Transition in Granular Systems, Phys. Rev. Lett. 98,
058001 (2007).
[25] A. S. Keys, A. R. Abate, S. C. Glotzer, and D. J. Durian, Measurement of growing dynamical length scales and prediction of
the jamming transition in a granular material, Nat. Phys. 3, 260(2007).
[26] M. Clusel, E. I. Corwin, A. O. N. Siemens, and J. Brujic, A “granocentric” model for random packing of jammed emulsions,
Nature (London) 460, 611 (2009).

α=

1
cos−1
2




γ (cos2 θ − sin2 θ ) − γ 2 sin θ cos θ. (C4)ab

The distance between the two disks must equal 2r, andthus (x1 − x2 )2 + (y1 − y2 )2 = (2r)2 . This expression yields aquadratic equation for r: Ar 2 + Br + C = 0, where A, B, andC are given by
A = 4 cot 2 α,
B = −4{a cot α + b[cot α cos(2α) + sin(2α)]},

(C5)
(C6)

and
C = a2 + 2ab cos(2α) + b2 .

−B ±034901-13

JERRY ZHANG et al.

PHYSICAL REVIEW E 108, 034901 (2023)[27] D. Bi, J. Zhang, B. Chakraborty, and R. P. Behringer, Jammingby shear, Nature (London) 480, 355 (2011).
[28] J. A. Dijksman, G. H. Wortel, L. T. H. van Dellen, O. Dauchot,and M. van Hecke, Jamming, Yielding, and Rheology of
Weakly Vibrated Granular Media, Phys. Rev. Lett. 107, 108303(2011).
[29] J. G. Puckett and K. E. Daniels, Equilibrating Temperaturelike
Variables in Jammed Granular Subsystems, Phys. Rev. Lett.
110, 058001 (2013).
[30] M. Z. Miskin and H. M. Jaeger, Adapting granular materials
through artificial evolution, Nat. Mater. 12, 326 (2013).
[31] J. Lin, I. Jorjadze, L.-L. Pontani, M. Wyart, and J. Brujic,
Evidence for Marginal Stability in Emulsions, Phys. Rev. Lett.
117, 208001 (2016).
[32] C. P. Goodrich, A. J. Liu, and S. R. Nagel, Solids between themechanical extremes of order and disorder, Nat. Phys. 10, 578(2014).
[33] C. E. Maloney and A. Lemaître, Amorphous systems in athermal, quasistatic shear, Phys. Rev. E 74, 016118 (2006).
[34] W. G. Ellenbroek, Z. Zeravcic, W. van Saarloos, and M. vanHecke, Non-affine response: Jammed packings vs. spring networks, Europhys. Lett. 87, 34004 (2009).
[35] W. G. Ellenbroek, E. Somfai, M. van Hecke, and W. van
Saarloos, Critical Scaling in Linear Response of Frictionless
Granular Packings near Jamming, Phys. Rev. Lett. 97, 258001(2006).
[36] A. Zaccone and E. Scossa-Romano, Approximate analytical description of the nonaffine response of amorphous solids, Phys.
Rev. B 83, 184205 (2011).
[37] P. Wang, S. Zhang, P. Tuckman, N. T. Ouellette, M. D. Shattuck,
and C. S. O’Hern, Shear response of granular packings compressed above jamming onset, Phys. Rev. E 103, 022902(2021).
[38] J. Zhang, K. VanderWerf, C. Li, S. Zhang, M. D. Shattuck, andC. S. O’Hern, Mechanical response of packings of nonsphericalparticles: A case study of two-dimensional packings of circulolines, Phys. Rev. E 104, 014901 (2021).
[39] S. Dagois-Bohy, B. P. Tighe, J. Simon, S. Henkes, and M. vanHecke, Soft-Sphere Packings at Finite Pressure but Unstable to
Shear, Phys. Rev. Lett. 109, 095703 (2012).
[40] C. P. Goodrich, A. J. Liu, and S. R. Nagel, Finite-Size Scaling at
the Jamming Transition, Phys. Rev. Lett. 109, 095704 (2012).
[41] H. Mizuno, L. E. Silbert, and M. Sperl, Spatial Distributions ofLocal Elastic Moduli Near the Jamming Transition, Phys. Rev.
Lett. 116, 068302 (2016).
[42] A. Zaccone and E. M. Terentjev, Short-range correlations control the g/k and poisson ratios of amorphous solids and metallic
glasses, J. Appl. Phys. 115, 033510 (2014).
[43] M. Schlegel, J. Brujic, E. M. Terentjev, and A. Zaccone, Local
structure controls the nonaffine shear and bulk moduli of disordered solids, Sci. Rep. 6, 18724 (2016).
[44] S. Zhang, W. Jin, D. Wang, D. Xu, J. Zhang, M. D. Shattuck,
and C. S. O’Hern, Local and global measures of the shear
moduli of jammed disk packings, Phys. Rev. E 107, 054903(2023).

[45] J. Thompson and K. Bailey, A survey of the elastic properties
of some semiconducting glasses under pressure, J. Non-Cryst.
Solids 27, 161 (1978).
[46] S. M. Antao, C. J. Benmore, B. Li, L. Wang, E. Bychkov, andJ. B. Parise, Network Rigidity in GeSe2 Glass at High Pressure,
Phys. Rev. Lett. 100, 115501 (2008).
[47] X. Zeng, R. Peng, Y. Yu, Z. Hu, Y. Wen, and L. Song, Pressure
effect on elastic constants and related properties of Ti3 Al intermetallic compound: A first-principles study, Materials 11, 10(2018).
[48] K. Zhang, B. Dice, Y. Liu, J. Schroers, M. D. Shattuck, andC. S. O’Hern, On the origin of multi-component bulk metallic
glasses: Atomic size mismatches and de-mixing, J. Chem. Phys.
143, 054501 (2015).
[49] G.-J. Gao, J. Blawzdziewicz, and C. S. O’Hern, Geometricalfamilies of mechanically stable granular packings, Phys. Rev. E
80, 061303 (2009).
[50] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids(Oxford University Press, Oxford, 2017).
[51] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch,
Structural Relaxation Made Simple, Phys. Rev. Lett. 97, 170201(2006).
[52] M. L. Falk and J. S. Langer, Dynamics of viscoplastic deformation in amorphous solids, Phys. Rev. E 57, 7192 (1998).
[53] G.-J. Gao, J. Bławzdziewicz, and C. S. O’Hern, Frequency
distribution of mechanically stable disk packings, Phys. Rev.
E 74, 061304 (2006).
[54] S. S. Ashwin, J. Blawzdziewicz, C. S. O’Hern, and M. D.
Shattuck, Calculations of the structure of basin volumes formechanically stable packings, Phys. Rev. E 85, 061307 (2012).
[55] J. C. Maxwell, On the calculation of the equilibrium and stiffness of frames, London, Edinburgh, Dublin Philos. Mag. J. Sci.
27, 294 (1864).
[56] S. Pellegrino and C. Calladine, Matrix analysis of statically andkinematically indeterminate frameworks, Int. J. Solids Struct.
22, 409 (1986).
[57] C. Brito, H. Ikeda, P. Urbani, M. Wyart, and F. Zamponi,
Universality of jamming of nonspherical particles, Proc. Natl.
Acad. Sci. USA 115, 11736 (2018).
[58] J. D. Treado, D. Wang, A. Boromand, M. P. Murrell, M. D.
Shattuck, and C. S. O’Hern, Bridging particle deformability andcollective response in soft solids, Phys. Rev. Mater. 5, 055605(2021).
[59] D. Wang, J. D. Treado, A. Boromand, B. Norwick, M. P.
Murrell, M. D. Shattuck, and C. S. O’Hern, The structural,
vibrational, and mechanical properties of jammed packings ofdeformable particles in three dimensions, Soft Matter 17, 9901
(2021).
[60] M. F. Staddon, A. Hernandez, M. J. Bowick, M. Moshe, and
M. C. Marchetti, The role of non-affine deformations in theelastic behavior of the cellular vertex model, Soft Matter 19,
3080 (2023).
[61] C. F. Schreck, M. Mailman, B. Chakraborty, and C. S. O’Hern,
Constraints and vibrations in static packings of ellipsoidal
particles, Phys. Rev. E 85, 061305 (2012).

034901-14