# A multi-body finite element model for hydrogel packings: Linear response to
shear

Ahmed Elgailani<EMAIL_ADDRESS>Craig E Maloney
<EMAIL_ADDRESS>Department of Mechanical and Industrial
Engineering, Northeastern University

###### Abstract

We study a multi-body finite element model of a packing of hydrogel particles
using the Flory-Rehner constitutive law to model the deformation of the
swollen polymer network. We show that while the dependence of the pressure,
$\Pi$, on the effective volume fraction, $\phi$, is virtually identical to a
monolithic Flory material, the shear modulus, $\mu$, behaves in a non-trivial
way. $\mu$ increases monotonically with $\Pi$ from zero and remains below
about $80\%$ of the monolithic Flory value at the largest $\Pi$ we study here.
The local shear strain in the particles has a large spatial variation. Local
strains near the centers of the particles are all roughly equal to the applied
shear strain, but the local strains near the contact facets are much smaller
and depend on the orientation of the facet. We show that the slip between
particles at the facets depends strongly on the orientation of the facet and
is, on average, proportional to the component of the applied shear strain
resolved onto the facet orientation. This slip screens the stress transmission
and results in a reduction of the shear modulus relative to what one would
obtain if the particles were welded together at the facet. Surprisingly, given
the reduction in the shear modulus arising from the facet slip, and the
spatial variations in the local shear strain inside the particles themselves,
the deformation of the particle centroids is rather homogeneous with the
strains of the Delaunay triangles having fluctuations of only order $\pm 5\%$.
These results should open the way to construction of quantitative estimates of
the shear modulus in highly compressed packings via mean-field, effective-
medium type approaches.

##  I Introduction

Particles made of hydrogels – chemically cross-linked polymeric materials
which swell when they absorb water – are important in industries such as
pharmaceuticals, bioengineering, agriculture, food science, and cosmetics [1,
2, 3, 4]. Owing to their softness and the in situ tunability of their size via
controlled swelling, they have also been used more recently as a test bed for
basic ideas in the physics of jamming and the glass transition [5, 6, 7, 8, 9,
10, 11]. The basic statistical mechanical description of the mechanical
properties of macroscopic, monolithic pieces of hydrogel material goes back to
Flory and Rehner [12, 13, 14, 15], however, the properties of suspensions or
packings of _particles_ made out of the hydrogel are much more subtle and have
been the subject of active research only over the past couple of decades.

Typically, in experiments, the particles in a suspension are initially allowed
to freely swell at low number density, then, once fully swollen, are subjected
to osmotic confinement with a membrane [11, 16, 17] or to centrifugation [18],
forcing the swollen particles into persistent contact and making a jammed
solid where the particles are no longer free to diffuse and pass each other
without incurring an energetic penalty associated with deformation of the
polymer network. At low degrees of confinement, the particles remain
essentially spherical, and contacts between them are only slightly deformed.
The mechanics and energetics of the interactions in this regime should be
reasonably well described by standard Hertzian contact mechanics [19] treating
one contact at a time resulting in a pairwise description of the energy of the
whole system. However, as the confinement increases, the inter-particle
contacts will develop into curved facets [20, 21], and, at large enough
confinement, the void space with pure solvent will completely disappear
altogether. In foams and emulsions at such large particle volume fractions, it
is known that not only does linear contact mechanics fail to describe the
interactions at a given contact, but even worse, pairwise descriptions of any
kind become qualitatively inaccurate in this regime [22, 23, 24], and there is
no reason to think that the micro gel packings of interest here will be any
different. One must consider in detail the deformation inside the particles.
It is this strongly confined, fully-faceted, regime which is of primary
interest to us here.

Several groups [11, 16, 17] have measured the elastic modulus, $G$, of highly
compressed microgel packings under osmotic confining pressure, $\Pi$, at
increasing nominal volume fraction, $\phi$. In all cases, one observes a rapid
increase of $\Pi$ and $G$ from, essentially, zero as the particles are first
forced into contact at the jamming transition when $\phi$ reaches the random
close packing point $\phi_{J}$. The observed behavior of $\Pi$ and $G$ is
qualitatively consistent with the jamming picture [25] where $\Pi$ and $G$ are
both zero below the jamming point, $\phi_{J}$, and, near jamming, scale as
non-trivial power of $\delta\phi=\phi-\phi_{J}$. Precise scaling laws for
$\Pi$ and $G$ are difficult to obtain in experiments, but in all experiments,
there is an obvious and dramatic onset of both $\Pi$ and $G$ near a critical
$\phi$ value. At higher $\phi$, away from jamming on approach to the fully
faceted regime, the situation becomes more complicated. Cloitre et. al. [11],
have shown that there is a transition from a more pronounced to a less
pronounced dependence of $G$ on $\phi$ and argued that it occurs at the point
at which the assembly becomes fully faceted with no void regions of pure
solvent remaining. Liétor-Santos et. al. [17] showed that $G$, $\Pi$, and, the
compression modulus, $K$, all scaled with $K_{p}$, the compression modulus of
an isolated particle osmotically confined to the same average particle size as
in the packing, so that $G/K_{p}$, $K/K_{p}$, and $\Pi/K_{p}$ were all
universal constants in the compressed packing independent of the degree of
compression of the packing or the degree of crosslinking of the particles.
Menut et. al. [16] have studied a variety of particles with different cross-
linking densities and sizes. They argued that sufficiently far above the
jamming point, the density dependence of the modulus followed the trend that
would be expected for a monolithic Flory material [12, 13]. Here, we find,
what is perhaps the simplest non-trivial outcome one could have expected: a
pressure vs. $\phi$ behavior almost precisely the same as the monolithic Flory
material and a universal shear modulus vs. pressure curve for systems with
different cross-linking density when the modulus and pressure are both scaled
by the Flory pressure, $NkT$, where $N$ is the density of cross-links in the
dry reference state.

##  II Model

Previous modeling work has proceeded along many fronts. Many groups assume
pair-wise additive interactions between particles even in a very high-density
regime far from jamming [26, 27, 20, 28, 29], sometimes modified to attempt to
account for $\phi$ dependent effective interactions [30, 31, 32, 33, 34, 24,
35], which are, however, still _pair-wise_ contact interactions. As we mention
above, this approach is known to fail in foams and emulsions due to strong
many-body effects in the particle-particle forces [24], and we would expect it
to fail as well here for the microgel packings in the high $\phi$, fully
faceted, regime. More realistic and appropriate models take into account the
nature of the deformation of the polymer network itself [36, 37, 27, 38, 39,
21]. Nikolova et. al. [39] studied a bead-spring coarse-grained model with
dissipative particle dynamics (DPD) to model both the elasticity of the gel
network and the dynamics of solvent expulsion out of the network as the
packing is confined. They allowed their particles to swell and subjected them
to isotropic osmotic confinement. They found that, above $\phi=1$, $K/K_{p}$
approached a constant value of about $0.8$, in agreement with the experiments
of Liétor-Santos et. al. [17], but they did not study the shear modulus.

Here, we take a different approach. Rather than using a bead-spring model to
explicitly represent the gel network, we represent the gel network as a
homogenized continuum using the Flory-Rehner constitutive law. This approach
is standard in the solid mechanic’s community where various non-linear elastic
properties of macroscopic swollen gels are of interest [40, 41, 42, 43],
however, perhaps surprisingly, it has not been applied to particle packings.
The model can easily incorporate spatial variations in local cross-linking
densities within a particle which may arise in various particle synthesis
procedures (e.g. a hard, moderately swollen, core with high cross-link density
enclosed by a soft, highly swollen, corona with low cross-link density), but
in this preliminary work, we assume a homogeneous cross-linking density across
the particle, and we assume the same cross-linking density for large and small
particles. The main disadvantage of our approach is that we cannot
realistically study the fluid dynamics of solvent uptake/expulsion from the
network and are restricted to compression and shear rates which are slow
enough that we can assume the hydrodynamic forces generated by solvent flows
are negligible. However, as we are primarily interested in the quasi-static
regime in this work where the packing is sheared slowly enough to allow the
fluid to be fully expelled/absorbed before further shearing, we are not
adversely affected by this limitation. Another potential disadvantage is that
we are not able to capture inter-digitation of polymer segments across inter-
particle facets [44, 39, 21, 45]. However, the effects of inter-digitation are
largely uncontrolled in bead-spring models; this will give rise to tangential
forces at the facets in a way that is not particularly well controlled or
calibrated. In our opinion, it is best to study and characterize a model which
at first excludes these effects explicitly and only later introduces them in a
controlled way. Furthermore, in [39], the authors tune the properties of their
bead-spring network to make sure that it behaves mechanically as a Flory
solid, whereas here, we simply start with a Flory solid by construction.

Flory’s constitutive law is governed by a total free energy density,
$W_{tot}$, which gives the free energy per unit (unswollen reference) volume
where $W_{tot}=W_{e}+W_{m}$ receives independent contributions from non-linear
elastic deformation of the polymer network, $W_{e}$, and from the free energy
of mixing of the polymer and solvent, $W_{m}$. We have, for $W_{e}$ and
$W_{m}$ [12, 13, 15, 14, 40],

$W_{e}=\frac{Nk_{b}T}{2}[F_{\alpha\beta}F_{\alpha\beta}-d-2\log(J)]$ (1)
$W_{m}=\frac{k_{b}T}{\Omega}(J-1)\left[\log\left(\frac{J-1}{J}\right)+\frac{\chi}{J}\right]$
(2)

where $N$ is the number of cross-links per unit volume _in the dry reference
state_ , $k_{b}$ is Boltzmann’s constant, $T$ is the temperature,
$F_{\alpha\beta}=\delta_{\alpha\beta}-\partial_{\alpha}u_{\beta}$ is the
deformation gradient where $u$ is the displacement of the polymer network away
from its unswollen reference configuration, $d$ is the dimension of the
system, $J$ is the determinant of $F_{\alpha\beta}$ (so that $J$ gives the
volumetric expansion at a given point relative to the dry reference and
$J^{1/d}=\lambda$ gives the linear stretch ratio), $\Omega$ is the volume of a
solvent particle (which is taken to be constant) and $\chi$ gives the
energetic contribution to the free energy of mixing. In this preliminary work,
we consider entropic swelling only and set $\chi=0$, so _there is one and only
one dimensionless parameter in the model_ : $N\Omega$, the number of cross-
links in the dry reference configuration per unit solvent particle volume. For
large $N\Omega$ the particles will be stiffer and undergo less free swelling,
while for small $N\Omega$, the particles will undergo more free swelling and
be softer. If we consider the isotropic swelling case with a linear stretch
factor, $\lambda$, we have, $F_{\alpha\beta}=\lambda\delta_{\alpha\beta}$ and
$J=\lambda^{d}$ so $F_{\alpha\beta}F_{\alpha\beta}=d\lambda^{2}=dJ^{2/d}$. The
osmotic pressure, $\Pi$, is simply the derivative of the free energy with
respect to the logarithm of the volume, $\Pi=J\frac{\partial W}{\partial J}$,
so

$\frac{\Pi}{Nk_{b}T}=(J^{2/d}-1)+\frac{1}{N\Omega}(1+J\log[(J-1)/J])$ (3)

The degree of free swelling is determined by setting $\Pi=0$ and solving for
the equilibrium, $J_{0}$. $J_{0}$ will depend on $N\Omega$, with larger
$N\Omega$ (more cross-linking) giving less unconstrained swelling [12, 13, 15,
14, 40]. We study three different $N\Omega$ values: $1/20,1/133,$ and $1/754$,
and solving for $\Pi=0$ gives: $J_{0}=4.00762,8.99301,20.2519$ respectively.
That is: our least densely cross-linked (softest) particles freely swell to
roughly $20$ times their dry area, while our most densely cross-linked
(hardest) particles freely swell to roughly $4$ times their dry area.

We simulate the system using standard finite element method (FEM) techniques
using constant strain triangle (CST) elements to mesh the particles [46] and
use a simple gradient descent to re-equilibrate the system. We use a 50:50
mixture of two species of circular particles with the ratio between the radii
equal to $1.4$ [25]. The disks are initialized on a square grid, allowed to
swell freely, and then slowly compressed.

To enforce impenetrability constraints, we simply introduce a strong power-law
repulsion between surface nodes of opposing particle FEM meshes:
$V_{ij}=V_{0}(r_{ij}/R)^{-a}$ where $r_{ij}$ is the distance between the
surface nodes on opposing particle meshes, $R$ is the length-scale of the
repulsion, $V_{0}$ is the energy scale of the repulsion, taken here to simply
be $NkT$, and $a=12$ is an exponent which is chosen to be large enough that it
prevents inter-penetration and does not otherwise affect the results. To avoid
artificial interlocking of surface nodes and the resulting tangential forces,
we choose a lengthscale for the power-law repulsion which is roughly $6$ times
the surface node spacing in the dry reference mesh, and we find that this
sufficiently suppresses tangential traction forces at the facets. Because of
the stiffness of the surface repulsion, the surface interactions do not
contribute significantly to the overall energy or stress, so we neglect them
when reporting stresses and moduli. However, the numerical scheme results in a
skin at the facets between the surface nodes which has a constant size which
does not change as the packing is subjected to contraction of the space. This
results in a larger and larger fraction of the space being taken up by the
spurious skin as the system is compressed. To compensate for this artifact, we
define the Cauchy stress, whose trace is the pressure, as the derivative with
respect to an infinitesimal strain increment of the total energy of the
deformed elastic network divided by the area occupied by the gel _excluding
the area of the skin region_. We have checked that the contribution of the
power-law repulsion to the total energy is negligible as long as the exponent
is chosen to be sufficiently large. We similarly define $\phi=J_{0}/\langle
J\rangle$ as the ratio of the total freely-swollen particle area to the
current total particle area occupied _excluding the skin_. This definition of
$\phi$ – which is, by construction, greater than or equal to $1$ for
compressed particles – is somewhat un-natural near the jamming onset where one
would want to divide the freely swollen particle volume by the area of the
simulation cell rather than just the area occupied by the polymer network. To
measure $\phi=J_{0}/\langle J\rangle$ in an experiment, one would need to have
an independent measurement of the pore volume in the sample to infer the
current total particle volume, but, nevertheless, it is the natural quantity
to use in the model to make a direct connection between the mechanical
response of the packing and the mechanical response of the equivalent
monolithic continuum.

##  III Results

Figure 1:  Local pressure, $\Pi$, in units of $NkT$ for a typical system of
$16$ particles with $N\Omega=1/133$ at a nominal $\phi=1.73$. Particles are
labeled with a unique identification number for discussion in the text.

In figure 1, we show an image of the local pressure in units of $NkT$ for a
typical system of $16$ particles with $N\Omega=1/133$ at a nominal $\phi=1.73$
which is in the fully faceted regime with no remaining voids of pure solvent.
Each triangular element in the FEM mesh is colored according to its pressure
(with the standard solid mechanics convention that compressive stresses are
negative). We observe some high-frequency oscillations from element to
element, but we have checked that these are suppressed using higher-order
finite element schemes and do not affect any of the main results of this work.
Some general trends can be seen here which are observed in other members of
the ensemble and at various other $N\Omega$ and various other $\phi$. i) Small
particles tend to have higher pressure than large particles. ii) Pressure
tends to be larger near the centers of facets than near the facet junctions.
iii) Facets between large and small particles tend to be convex on the small-
particle side and concave on the large-particle side: that is, the small
particles tend to protrude into the large ones. iv) Particles with fewer
facets tend to be at larger pressure and vice versa. For instance, particle
$14$, a small particle which has four neighbors, is at a larger pressure than
the other small particles which have $5$ or $6$ neighbors, while particle
$10$, a large particle which has eight neighbors, is at a smaller pressure
than the other large particles which have $7$ or $6$ neighbors.

Figure 2:  Main: Shear modulus vs. osmotic pressure, $\Pi$, for various cross-
link densities. The red line is the shear modulus of a two-dimensional
monolithic Flory material. Inset: $\Pi/NkT$ vs $\phi$. Dots are data from the
FEM simulation. Solid lines are for a monolithic Flory solid.

In the inset of figure 2, we show the pressure, $\Pi$, vs $\phi$ (recall the
definition of $\phi=J_{0}/J$ excludes regions of pure solvent). The results
are for an ensemble of $36$-particle systems. The symbols represent the
pressure computed in the actual packing, and the curves simply represent the
pressure for an isotropic monolith of Flory material as computed from equation
3 at the given $J=J_{0}/\phi$. The agreement is striking! The pressure field
inside the particles is inhomogeneous: it is large near the facets and smaller
near the vertices in the facet network with overall fluctuations of order tens
of percent. Yet despite this inhomogeneity, the _spatial average_ of the
pressure is very close to that of the isotropically confined Flory monolith at
the same density as the _average density_ of the packing. There is no
fundamental reason why this relation must be an exact identity, and we do see
deviations from the monolithic result of the order of a few percent. It is a
manifestation of the fact that the distribution of local $J$ values has a
width that is small relative to its average and that the $\Pi(J)$ function
from equation 3 varies relatively slowly on that scale so that
$\langle\Pi(J)\rangle\approx\Pi(\langle J\rangle)$.

The pressure diverges locally as $J\rightarrow 1$ from above, and this sets an
upper bound on $\phi$ at $\phi_{\text{max}}=J_{0}$ as the solvent is
completely expelled from the system. Accordingly, we see $\Pi$ begin to
diverge first for the system with the highest $N\Omega$ and last for the
system with the lowest $N\Omega$. Our numerical scheme becomes unstable at a
pressure of about $\Pi=10Nk_{b}T$ for all systems, so we are able to reach
slightly higher $\phi$ for the systems with lower $N\Omega$.

To measure the shear modulus, we make a small axial deformation to all FEM
nodes and the periodic boundaries with extension along the horizontal, $x$,
and contraction along the vertical, $y$: $x\rightarrow
e^{\epsilon}x\approx(1+\epsilon)x,y\rightarrow
e^{-\epsilon}y\approx(1-\epsilon)y$, re-equilibrate, and then measure the
tangent modulus as $\mu=\Delta\sigma/\epsilon$ where $\Delta\sigma$ is the
change in the shear component of the Cauchy stress. Since we hold the cell in
a square shape during the initial swelling of the circular particles, there
will be a random residual shear stress which is distributed normally about
zero in our ensemble. It is very small compared to $\mu$ and does not affect
the results.

In the main plot of figure 2, we show the shear modulus, $\mu$, vs the osmotic
pressure, $\Pi$. For a monolithic Flory material in 2D, $\mu=Nk_{b}T$
regardless of $N\Omega$ and completely independent of $\Pi$ (and/or $\phi$).
For our particle packing, the curves all start at zero near the jamming point
at $\Pi=0$ and increase monotonically remaining below the monolithic Flory
value. Strikingly, the curves collapse, indicating that the shear modulus of
the packing is a function of the pressure alone and is independent of the
cross-linking density (after scaling by $NkT$). This observation for the shear
modulus is in a similar spirit to the observation of Liétor-Santos et. al.
[17] that the pressure and shear modulus scale like the single-particle
compression modulus, but must differ in detail as we argue below.

Figure 3: (a) Non-affine displacement field normalized by the imposed shear
strain. The imposed shear strain is axial with vertical contraction
(horizontal extension) and area-preserving at $\phi$ = 2.1 for a system with
cross-linking density $N\Omega$ = 1/20. Arrows are drawn to scale to precisely
give the non-affine displacement per unit applied strain. (b) Local shear
strain, $(\partial_{x}u_{x}-\partial_{y}u_{y})/2$, scaled by the imposed shear
strain, $e1$.

In figure 3 (a), we show a typical $N\Omega=1/20$ system with $36$ particles
at $\phi=2.09$ subjected to a small strain step as described above. We show
only the non-affine component of the nodal displacements. In figure 3 (b), we
show the combination of the components of the gradient of the displacement
corresponding to the applied strain: $(\partial_{x}u_{x}-\partial_{y}u_{y})/2$
normalized by the applied strain so that a value of $1$ indicates that the
material is locally shearing precisely according to the globally imposed shear
strain. The strain is quite inhomogeneous across the particles. In the centers
of all the particles, the shear strain is approximately equal to the globally
imposed shear strain, while the deformation near facets and vertices depends
on the orientation. Facets that are roughly vertical or horizontal have
essentially no relative displacement, while facets that are roughly diagonal
across the cell have a significant amount of relative tangential motion
indicating sliding along those facets.

Figure 4: Closeup on the facet between particles $25$ and $31$ from figure 3.
The triple junction with particle $26$ on the top and $30$ on the bottom are
labeled $T$ and $B$ respectively for ease of discussion in the text.

In figure 4, to illustrate this point, we zoom in on the facet between
particles $25$ and $31$, since this facet is nearly vertical. For ease of
discussion, we label the triple junction with particle $26$ on top and with
particle $30$ on the bottom as $T$ and $B$ respectively. Since the facet is a
contact between a large and small particle, it has a slight curvature toward
the small particle, but this should not significantly affect the response. The
facet is nearly vertical, and there is essentially no non-affine displacement
on either side. However, the displacement on particle $26$ at the $T$ vertex
is downward, and the displacement on particle $30$ at the $B$ vertex is
upward, so the facet is _shortening_. As particles $26$ and $30$ are advancing
into the facet, the strain near their vertices at the triple junctions is
nearly zero: green in the color scheme of the image. On the other hand, the
material near the four vertices on particles $25$ and $31$ at the $T$ and $B$
triple junctions is very strongly sheared – with local shear strains greater
than the applied shear strain – to accommodate the shortening of the facet.
The two other facets connected to $T$ and the two other facets connected to
$B$ are all oriented approximately $30$ degrees from horizontal and have a
significant amount of slipping. Since the tangential traction forces at the
slipping interface are zero, and the shear stress must be zero along the
facet, the slip results in a screened region where the strains in the
particles are quite small and nearly zero at the facet.

The same scenario plays out in reverse for facets which are oriented
horizontally. The result is that all particle vertices in the packing whose
opposing facet is nearly horizontal or vertical should have a very low strain
in their neighborhood, and a quick scan of the packing shows this to be true.
For example, the facet between particles $34$ and $35$ is nearly horizontal,
so it is lengthening, and the corresponding vertices on particles $28$ and $4$
have a very low shear strain on them: the fact that particle $28$ has $5$
neighbors and is compressed into a pentagon is incidental and has little
impact on this result.

Figure 5:  a) Facet slip , $\omega$ (defined in text) scaled by applied strain
$\epsilon$ vs. facet angle $\theta$ for three different $N\Omega$ systems at
$\phi=1.5,1.9$ and $2.1$. b) Facet slip in real space with arrow length
indicating the amount of slip with red and blue indicating counter-clockwise
and clockwise vorticity respectively.

Figure 6: a) Facet strain, $\epsilon_{L}$ (defined in text) scaled by applied
strain $\epsilon$ vs. facet angle $\theta$ for three different $N\Omega$
systems at $\phi=1.5,1.9$ and $2.1$. b) Facet strain in real space with bar
length indicating the amount of relative lengthening (red) or shortening
(blue).

In figures 5 and 6, we show the slip and strain on the facets as a function of
the facet angle. We approximate the facet network as polygonal and identify
the vertices of the polygonal network with the surface nodes of the mesh
located at the triple junctions. Then, to define the facet slip, $\omega$, we
simply take the difference in the average tangential component of displacement
on either side of the facet. Similarly, to define the facet strain,
$\epsilon_{L}$, we take the difference in the component of the end-to-end
displacement along the end-to-end separation and divide by the current facet
length. What we see agrees remarkably well with the anecdotal description of
the displacement fields in figure 3 where facets along the box diagonals slip
significantly but neither shorten nor lengthen while facets along the axes do
not slip but shorten (the vertical ones) or lengthen (the horizontal ones). In
fact, we can make some simple assumptions about the deformation kinematics. If
we assume the slip on a facet is equal to the locally resolved transverse
component of the applied displacement gradient:
$\omega=\hat{n}_{\alpha}\hat{t}_{\beta}\partial_{\alpha}u_{\beta}$ where
$\partial_{\alpha}u_{\beta}$ is the imposed deformation and $\hat{n}$ and
$\hat{t}$ are the unit normal and tangent to the facet, then we would get:
$\omega=e_{1}\sin(2\theta)$ where $e_{1}$ is the amplitude of the imposed
strain and $\theta$ is the angle of the tangent. For the strain, if one
imagines that the triple junctions in the facet network deform affinely with
the imposed shear, one would expect that
$\epsilon_{L}=\hat{t}_{\alpha}\hat{t}_{\beta}\partial_{\alpha}u_{\beta}=e_{1}\cos(2\theta)$.
We see that the data follow that trend remarkably well regardless of $\phi$
with no discernible trend with $\phi$. Our argument seems to slightly
overestimate the slip and underestimate the strain, but all things considered,
this prediction with no adjustable parameters seems to work out well: the
deformation of the facet network is essentially affine and equal to the
homogeneously imposed deformation.

Figure 7: Shear strain of the particle centroids scaled by the applied strain.
The strain is computed for each Delaunay triangle in the contact network by
linear interpolation. Slip vectors are reproduced from the figure for
convenience 5.

In figure 7, we show the strains associated with the motion of the particle
centroids. We define the strain as piecewise constant on the Delaunay
triangulation [47] of the contact network using linear interpolation of the
three centroid displacements. What we see is that the resulting deformation of
the particle centroids is quite homogeneous despite the large inhomogeneities
_within_ any given particle. It is hard to draw any statistical trends or make
correlations between the Delaunay strain and any obvious geometrical or
topological properties of the packing.

##  IV Discussion and summary

Figure 8:  a) Compression modulus, $K/NkT$, vs. osmotic pressure, $\Pi/NkT$,
for a monolithic Flory material in the range of pressures of interest in this
work $0<\Pi/NkT<10$. b) Their ratio, $\Pi/K$, vs $\Pi/NkT$.

We have shown here that the microgel packings behave in one of the simplest
ways one could have imagined: the pressure has a density dependence as if it
were a monolithic Flory material, while the shear modulus, although it has a
non-trivial $\phi$ dependence and starts from zero at the jamming point, is a
universal function of the dimensionless pressure, $\Pi/NkT$, _independent of
cross-linking densities_. The $\mu/NkT$ vs $\Pi/NkT$ curve shows a transition
from a more strongly pressure-dependent regime at low pressure near jamming
onset to a much weaker pressure dependence in the fully-faceted regime and is
consistent with the conjecture in Cloitre et. al. [11] that the transition to
the weaker concentration dependence of $G$ occurred at the onset of full-
faceting where solvent-pure void space completely disappeared. This is also
completely consistent with the arguments of Menut et. al. [16] who claimed
that their systems showed monolithic Flory behavior at high density (although
we point out that the scaling behavior for a Flory monolith in $3D$ is
$\mu\sim\phi^{1/3}$ rather than the $\mu\sim\phi^{1}$ suggested in that work
[15, 14]

At face value, our results would seem to be inconsistent with Liétor-Santos
et. al [17] where $\Pi$ became independent of $\phi$ when scaled by the
particle compression modulus, $K_{p}$, rather than the Flory modulus, $NkT$.
However, in figure 8, we plot the compression modulus, $K$ vs. the pressure,
$\Pi$, for our three $N\Omega$ values for a homogeneously compressed Flory
solid where $K=J\frac{\partial\Pi}{\partial J}$. The compression modulus for
the zero pressure freely swollen state is around $K_{eq}\approx 2NkT$ and has
only a weak dependence on cross-linking density $N\Omega$. We see that $\Pi/K$
does not vary dramatically over the pressure range, and for the two more
weakly cross-linked systems, $\Pi/K\approx 0.4$ above about $\Pi=4NkT$. Since
our $\Pi$ vs. $\phi$ curves for the packing are virtually identical to the
monolith, this means that our results are roughly consistent with reference
[17] where $\Pi/K$ was roughly constant above $\phi=1$. Our results for the
shear modulus, in distinction with the pressure and compression modulus, are
not as easy to reconcile with reference [17]. If we were to scale the shear
modulus by $K$ rather than $NkT$, because of the $N\Omega$ dependence of $K$
illustrated in figure 8, we would break the good collapse obtained in figure
2. Furthermore, the $\phi$-independent values of $G/K_{p}$ obtained in [17]
are of the order of a few percent, whereas here, our values of $\mu/NkT$ are
around $0.8$ for most of the pressure range of interest; our values for
$\mu/K$ would vary much more strongly with pressure than our $\mu/NkT$ do.

What we’ve shown here is that, once the system is in the fully faceted regime,
the kinematics of shear deformation becomes largely insensitive to the
compression. The vertices in the facet network, on average, follow the imposed
homogeneous deformation. The resulting slip at the facets results in a
reduction of the modulus of the packing away from the monolithic Flory value
by some tens of percent but remains of the order of the monolithic value. In
the future, it would be important to use the kinematic information about the
affine motion of the facet-vertex network to try to construct quantitative
estimates for the modulus reduction away from the Flory value. Since each
particle is assumed to be homogeneous, the result of imposing a homogeneous
deformation on the particle’s boundary would be a homogeneous deformation of
the particle interior, and the result would be a recovery of the full Flory
shear modulus as if the particles were welded together at the facets. So it
seems, going beyond considering the average kinematic behavior of the facet
network would be necessary to obtain any corrections to the Flory shear
modulus.

In this work, we made many simplifying assumptions. We assumed that: i) the
particles were homogeneous disks, ii) the free energy of mixing was completely
entropic with $\chi=0$, iii) the facets were free of friction and neglected
any possible effects of inter-digitation of the polymer networks across the
facets, iv) the deformation was slow enough that kinetics and viscous flow
effects were negligible and that the solvent was able to freely flow into and
out of the polymer network with no resistance v) the particles were large
enough that Brownian effects were negligible, vi) the system was 2D. It will
be useful to go beyond these assumptions in future models. It would also be
important to know how the many-body interactions studied here at large volume
fraction affect the plastic yielding behavior at large strains and the glass
transition for the Brownian case, and, more generally, the overall rheological
response at arbitrary strains and strain rates. Despite these interesting
future directions for study, our explanation for the corrections to the
monolithic Flory behavior induced by the affine facet slip provides a starting
point for future quantitative models of the mechanical response of these
packings.

##  V Acknowledgements

We thank Alberto Fernandez-Nieves, Emanuela Zaccarelli, Lorenzo Rovigatti,
Michel Cloitre, Robin Selinger, and especially Massimo Pica Ciamarra, Max Bi,
and Jeffrey Sokoloff for useful discussions in the early stages of this work.
We thank Manuel Valera for their help with code development and optimization.
We thank MGHPCC for computing support. This material is based upon work
supported by the National Science Foundation under Grant No. CMMI-1822020.

## References

  * [1] Eugene Pashkovski.  Applications of Biopolymer Microgels, chapter 17, pages 423–450.  John Wiley and Sons, Ltd, 2011.
  * [2] L. Andrew Lyon, Grant R. Hendrickson, Zhiyong Meng, and Ashlee N. St. John Iyer.  Exploiting the Optical Properties of Microgels and Hydrogels as Microlenses and Photonic Crystals in Sensing Applications, chapter 14, pages 355–374.  John Wiley and Sons, Ltd, 2011.
  * [3] Yuxing Ben, Ian Robb, Peng Tonmukayakul, and Qiang Wang.  Microgels for Oil Recovery, chapter 16, pages 407–422.  John Wiley and Sons, Ltd, 2011.
  * [4] Martin Malmsten.  Microgels in Drug Delivery, chapter 15, pages 375–405.  John Wiley and Sons, Ltd, 2011.
  * [5] Praveen Agarwal and Lynden A. Archer.  Strain-accelerated dynamics of soft colloidal glasses.  PHYSICAL REVIEW E, 83(4, 1), APR 11 2011.
  * [6] Johan Mattsson, Hans M. Wyss, Alberto Fernandez-Nieves, Kunimasa Miyazaki, Zhibing Hu, David R. Reichman, and David A. Weitz.  Soft colloids make strong glasses.  Nature, 462(7269):83–86, 2009.
  * [7] Anindita Basu, Ye Xu, Tim Still, P. E. Arratia, Zexin Zhang, K. N. Nordstrom, Jennifer M. Rieser, J. P. Gollub, D. J. Durian, and A. G. Yodh.  Rheology of soft colloids across the onset of rigidity: scaling behavior, thermal, and non-thermal responses.  SOFT MATTER, 10(17):3027–3035, 2014.
  * [8] Ludovic Berthier, Elijah Flenner, Hugo Jacquin, and Grzegorz Szamel.  Scaling of the glassy dynamics of soft repulsive particles: A mode-coupling approach.  PHYSICAL REVIEW E, 81(3, 1), MAR 2010.
  * [9] M Cloitre, R Borrega, F Monti, and L Leibler.  Glassy dynamics and flow properties of soft colloidal pastes.  PHYSICAL REVIEW LETTERS, 90(6), FEB 14 2003.
  * [10] M Cloitre, R Borrega, and L Leibler.  Rheological aging and rejuvenation in microgel pastes.  PHYSICAL REVIEW LETTERS, 85(22):4819–4822, NOV 27 2000.
  * [11] M Cloitre, R Borrega, F Monti, and L Leibler.  Structure and flow of polyelectrolyte microgels: from suspensions to glasses.  COMPTES RENDUS PHYSIQUE, 4(2):221–230, MAR 2003.
  * [12] Paul J. Flory and John Rehner.  Statistical mechanics of cross‐linked polymer networks i. rubberlike elasticity.  The Journal of Chemical Physics, 11(11):512–520, 1943.
  * [13] Paul J. Flory and John Rehner.  Statistical mechanics of cross‐linked polymer networks ii. swelling.  The Journal of Chemical Physics, 11(11):521–526, 1943.
  * [14] Masao Doi.  Soft Matter Physics.  Oxford, 2013.
  * [15] M Rubinstein and Ralph H. Colby.  Polymer Physics.  Oxford, 2003.
  * [16] Paul Menut, Sebastian Seiffert, Joris Sprakel, and David A. Weitz.  Does size matter? Elasticity of compressed suspensions of colloidal- and granular-scale microgels.  SOFT MATTER, 8(1):156–164, 2012.
  * [17] Juan Jose Lietor-Santos, Benjamin Sierra-Martin, and Alberto Fernandez-Nieves.  Bulk and shear moduli of compressed microgel suspensions.  PHYSICAL REVIEW E, 84(6, 1), DEC 8 2011.
  * [18] K. N. Nordstrom, E. Verneuil, W. G. Ellenbroek, T. C. Lubensky, J. P. Gollub, and D. J. Durian.  Centrifugal compression of soft particle packings: Theory and experiment.  PHYSICAL REVIEW E, 82(4, 1), OCT 8 2010.
  * [19] K. L. Johnson.  Contact Mechanics.  Cambridge University Press, Cambridge, 1985.
  * [20] Frank Scheffold.  Pathways and challenges towards a complete characterization of microgels.  Nature Communications, 11(1):4315, 2020.
  * [21] A. Scotti, M. Brugnoni, A. A. Rudov, J. E. Houston, I. I. Potemkin, and W. Richtering.  Hollow microgels squeezed in overcrowded environments.  The Journal of Chemical Physics, 148(17):174903, 2018.
  * [22] Reinhard Hoehler and Sylvie Cohen-Addad.  Many-body interactions in soft jammed materials.  SOFT MATTER, 13(7):1371–1383, FEB 21 2017.
  * [23] D. Weaire, R. Höhler, and S. Hutzler.  Bubble-bubble interactions in a 2d foam, close to the wet limit.  Advances in Colloid and Interface Science, 247:491–495, 2017.  Dominique Langevin Festschrift: Four Decades Opening Gates in Colloid and Interface Science.
  * [24] Reinhard Höhler and Denis Weaire.  Can liquid foams and emulsions be modeled as packings of soft elastic particles?  Advances in Colloid and Interface Science, 263:19–37, 2019.
  * [25] CS O’Hern, LE Silbert, AJ Liu, and SR Nagel.  Jamming at zero temperature and zero applied stress: The epitome of disorder.  PHYSICAL REVIEW E, 68(1, 1), JUL 2003.
  * [26] Maxime J Bergman, Nicoletta Gnan, Marc Obiols-Rabasa, Janne-Mieke Meijer, Lorenzo Rovigatti, Emanuela Zaccarelli, and Peter Schurtenberger.  A new look at effective interactions between microgel particles.  Nat Commun, 9(1):5039, Nov 2018.
  * [27] I. Bouhid de Aguiar, T. van de laar, M. Meireles, A. Bouchoux, J. Sprakel, and K. Schroen.  Deswelling and deformation of microgels in concentrated packings.  SCIENTIFIC REPORTS, 7, AUG 31 2017.
  * [28] Gaurasundar M. Conley, Chi Zhang, Philippe Aebischer, James L. Harden, and Frank Scheffold.  Relationship between rheology and structure of interpenetrating, deforming and compressing microgels.  Nature Communications, 10(1):2436, 2019.
  * [29] Ludovic Berthier, Angel J. Moreno, and Grzegorz Szamel.  Increasing the density melts ultrasoft colloidal glasses.  PHYSICAL REVIEW E, 82(6, 1), DEC 14 2010.
  * [30] Matthew Urich and Alan R. Denton.  Swelling, structure, and phase stability of compressible microgels.  SOFT MATTER, 12(44):9086–9094, 2016.
  * [31] Tianfei Liu, Fardin Khabaz, Roger T. Bonnecaze, and Michel Cloitre.  On the universality of the flow properties of soft-particle glasses.  Soft Matter, 14:7064–7074, 2018.
  * [32] Ruben Higler and Joris Sprakel.  Apparent strength versus universality in glasses of soft compressible colloids.  Scientific Reports, 8(1):16817, 2018.
  * [33] Martin-D. Lacasse, Gary S. Grest, Dov Levine, T. G. Mason, and D. A. Weitz.  Model for the elasticity of compressed emulsions.  Phys. Rev. Lett., 76:3448–3451, Apr 1996.
  * [34] T. G. Mason, Martin-D. Lacasse, Gary S. Grest, Dov Levine, J. Bibette, and D. A. Weitz.  Osmotic pressure and viscoelastic shear moduli of concentrated emulsions.  Phys. Rev. E, 56:3150–3166, Sep 1997.
  * [35] Reinhard Höhler and Sylvie Cohen-Addad.  Many-body interactions in soft jammed materials.  Soft Matter, 13:1371–1383, 2017.
  * [36] Giovanni Romeo and Massimo Pica Ciamarra.  Elasticity of compressed microgel suspensions.  SOFT MATTER, 9(22):5401–5406, 2013.
  * [37] Lorenzo Rovigatti, Nicoletta Gnan, Letizia Tavagnacco, Angel J. Moreno, and Emanuela Zaccarelli.  Numerical modelling of non-ionic microgels: an overview.  SOFT MATTER, 15(6):1108–1119, FEB 14 2019.
  * [38] Lorenzo Rovigatti, Nicoletta Gnan, and Emanuela Zaccarelli.  Internal structure and swelling behaviour of in silico microgel particles.  JOURNAL OF PHYSICS-CONDENSED MATTER, 30(4), JAN 31 2018.
  * [39] Svetoslav V. Nikolov, Alberto Fernandez-Nieves, and Alexander Alexeev.  Behavior and mechanics of dense microgel suspensions.  Proceedings of the National Academy of Sciences, 117(44):27096–27103, 2020.
  * [40] Shengqiang Cai and Zhigang Suo.  Equations of state for ideal elastomeric gels.  EPL, 97(3), FEB 2012.
  * [41] Zishun Liu, Wei Hong, Zhigang Suo, Somsak Swaddiwudhipong, and Yongwei Zhang.  Modeling and simulation of buckling of polymeric membrane thin film gel.  COMPUTATIONAL MATERIALS SCIENCE, 49(1):S60–S64, JUL 2010.  Symposium Q: Computational Materials Design at All Scales: From Theory to Application International Conference on Materials for Advanced Technologies (ICMAT), undefined, SINGAPORE, JUN 28-JUL 03, 2009.
  * [42] Shengqiang Cai and Zhigang Suo.  Mechanics and chemical thermodynamics of phase transition in temperature-sensitive hydrogels.  JOURNAL OF THE MECHANICS AND PHYSICS OF SOLIDS, 59(11):2259–2278, NOV 2011.
  * [43] Wei Hong, Zishun Liu, and Zhigang Suo.  Inhomogeneous swelling of a gel in equilibrium with a solvent and mechanical load.  INTERNATIONAL JOURNAL OF SOLIDS AND STRUCTURES, 46(17):3282–3289, AUG 15 2009.
  * [44] Priti S. Mohanty, Sofi Nöjd, Kitty van Gruijthuijsen, Jérôme J. Crassous, Marc Obiols-Rabasa, Ralf Schweins, Anna Stradner, and Peter Schurtenberger.  Interpenetration of polymeric microgels at ultrahigh densities.  Scientific Reports, 7(1):1487, 2017.
  * [45] Andrea Scotti, M. Friederike Schulte, Carlos G. Lopez, Jérôme J. Crassous, Steffen Bochenek, and Walter Richtering.  How softness matters in soft nanogels and nanogel assemblies.  Chemical Reviews, 122(13):11675–11700, 07 2022.
  * [46] David V. Hutton.  Fundamentals of finite element analysis.  2003\.
  * [47] Franco P. Preparata and Michael I. Shamos.  Computational Geometry: An Introduction.  Springer, 1985.