# Statistical Field Theory for Nonlinear Elasticity of Polymer Networks with
Excluded Volume Interactions

Pratik Khandagale<EMAIL_ADDRESS>Department of Mechanical
Engineering, Carnegie Mellon University  Timothy Breitzman  Air Force Research
Laboratory  Carmel Majidi  Department of Mechanical Engineering, Carnegie
Mellon University  Department of Civil and Environmental Engineering, Carnegie
Mellon University  Department of Materials Science and Engineering, Carnegie
Mellon University  Kaushik Dayal  Department of Civil and Environmental
Engineering, Carnegie Mellon University  Department of Mechanical Engineering,
Carnegie Mellon University  Center for Nonlinear Analysis, Department of
Mathematical Sciences, Carnegie Mellon University  Pittsburgh Quantum
Institute, University of Pittsburgh

###### Abstract

Polymer networks formed by cross-linking flexible polymer chains are
ubiquitous in many natural and synthetic soft matter systems. Current
micromechanics models generally do not account for excluded volume
interactions except, for instance, through imposing a phenomenological
incompressibility constraint at the continuum-scale. This work aims to examine
the role of excluded volume interactions on the mechanical response.

The approach is based on the framework of the self-consistent statistical
field theory of polymers, which provides an efficient mesoscale approach that
enables the accounting of excluded volume effects without the expense of
large-scale molecular modeling. A mesoscale representative volume element is
populated with multiple interacting chains, and the macroscale nonlinear
elastic deformation is imposed by mapping the end-to-end vectors of the chains
by this deformation. In the absence of excluded volume interactions, it
recovers the closed-form results of the classical theory of rubber elasticity.
With excluded volume interactions, the model is solved numerically in
3-dimensions using a finite element method to obtain the energy, stresses, and
linearized moduli under imposed macroscale deformation. Highlights of the
numerical study include: (1) the linearized Poisson’s ratio is very close to
the incompressible limit without a phenomenological imposition of
incompressibility; (2) despite the harmonic Gaussian chain as a starting
point, there is an emergent strain-softening and strain-stiffening response
that is characteristic of real polymer networks, driven by the interplay
between the entropy and the excluded volume interactions; and (3) the
emergence of a deformation-sensitive localization instability at large
excluded volumes.

††preprint: Phys. Rev. E 107, 064501 (DOI: 10.1103/PhysRevE.107.064501)

##  1 Introduction

A wide variety of soft-matter based systems are emerging as important for
engineering and scientific applications, and have been the focus of research
using both modeling and experiments, e.g. [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
31, 32]. Polymer-network based materials such as elastomers and hydrogels are
often at the heart of these soft-matter systems.

An important question for both fundamental understanding and application is to
predict the nonlinear elastic properties of polymer networks starting from a
micromechanical model of individual chains. The physics of polymer network
elasticity is governed by the conformational entropy of polymer chains and the
inter-segment excluded volume interactions. These contributions can be roughly
thought of as short-range and nonlocal interactions respectively. The short-
range interactions are associated with Gaussian polymer chain response and
depend on the relative configurations of adjacent segments in a chain. In
contrast, the nonlocal interactions are due to the interaction between polymer
segments that are nearby in space but nonlocal topologically (i.e., in terms
of their position along the chain (Fig. 1)).

While there are several useful phenomenological nonlinear elastic frame-
indifferent models, e.g. Mooney-Rivlin [33], Ogden [34] and Gent [35], they
lack a clear connection to the molecular structure of polymer network. An
important class of physics-based approaches to study the elasticity of polymer
networks are based on considering multiple Gaussian chains and then averaging
over the chains in different ways. These include the 3-chain model by James
and Guth [36], the 4-chain model by Flory and Rehner [37] and Treloar [38],
the affine full-network model by Treloar [39], the 8-chain model by Arruda and
Boyce [40], the non-homogeneous deformation based model by Wu and van Der
Giessen [41], and the non-affine microsphere model by Miehe [42]; the recent
work by Grasinger [43] provides a new perspective in which these myriad models
are shown to be special cases of a general approach. While these models have
provided important insights and prediction, they do not account for the
nonlocal excluded volume effects. Consequently, incompressibility of the
polymer network must be added as a phenomenological continuum-scale
approximation of the missing mesoscale physics. Another class of physics-based
molecular-statistical approaches are the constrained junction and constrained
segment theories, that aim to account for constraints arising due to chain
entanglements. Constrained junction theories, e.g. [44, 45, 46, 47, 48], apply
topological constraints on the fluctuations of chain cross-link junctions.
Constrained segment theories, e.g. [49, 50, 51, 52, 53], which are consistent
with the tube model of rubber elasticity, incorporate constraints on the
polymer segments along the chain contour. However, it is not easy to
incorporate the nonlocal excluded volume interactions in these approaches.

Figure 1: Excluded volume interactions are nonlocal in terms of the segment
coordinates.

###  1.A The Proposed Approach

Our approach is composed of 2 key elements: first, the statistical field
theory of polymers which provides an established and efficient approach to
account for the physics of polymer chain elasticity as well as excluded volume
interactions [54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69,
70]; and, second, the use of the 8-chain network averaging model that provides
a nonlinearly-elastic frame-indifferent approach to coarse-grain to the
continuum scale [40]. An important work in this direction is [71], wherein a
network with a simplified square lattice topology was studied using the field
theory approach to understand copolymers.

We begin by considering a representative volume element (RVE) of the polymer
network. A typical mesoscale RVE consists of several polymer chains that are
all interacting with configurations that are randomly distributed. While it is
a significant challenge to account for this randomness, we follow the 8-chain
RVE-averaging approach of Arruda and Boyce (Fig. 2, [40]) in approximating the
RVE in the undeformed state as composed of 8 polymer chains connecting the
center of a cube to each of the corners. The RVE then deforms under the action
of the macroscopic deformation tensor $\displaystyle{\mathbold F}$, i.e., the
chain end-to-end vectors are mapped by $\displaystyle{\mathbold F}$ from the
undeformed to the deformed state (Fig. 3, [72]). An important element of [40]
is that the RVE is oriented such that the cube is oriented along the principal
directions of the stretch tensor $\displaystyle{\mathbold U}$, where
$\displaystyle{\mathbold U}$ is the tensor square root of
$\displaystyle{\mathbold F}$, and can be obtained through the polar
decomposition $\displaystyle{\mathbold F}={\mathbold R}{\mathbold U}$ where
$\displaystyle{\mathbold R}\in SO(3)$.

Given the mapping of the end-to-end vectors of the chains, the polymer field
theory is then used to compute the partition function of the deformed state,
from which we can find the free energy and stress. Following [59], we use the
continuous Gaussian chain model for a single polymer chain. Next, we consider
chain segments that interact pairwise in real-space – and nonlocally in terms
of position along the chain contour (Fig. 1) – through a pairwise interaction
potential of mean force; these are given by Dirac potentials to model excluded
volume effects. Given the inter-segment interaction and end-to-end vectors,
the framework of polymer field theory enables us to compute, using the self-
consistent scheme, the partition function and consequently the free energy of
the RVE. We notice that because the ends of the polymer chain are constrained
by the macroscale deformation $\displaystyle{\mathbold F}$, this leads to a
restricted ensemble. Further, nonlinear elasticity provides the Piola-Kirchoff
stress tensor $\displaystyle{\mathbold P}$ as the energy-conjugate of
$\displaystyle{\mathbold F}$, enabling us to compute the stress-deformation
response of the polymer network.

Key results from the model are as follows. In the absence of excluded volume
interactions, we find that the closed-form orientationally-averaged elastic
response matches with classical rubber elasticity [72]. Considering excluded
volume interactions, closed-form solutions appear impossible, and we develop a
3-d finite element method (FEM) implementation to self-consistently solve the
equations of the polymer field theory. We find that the linearized Poisson’s
ratio $\displaystyle\nu\simeq 0.4943$, which is very close to the
incompressible limit $\displaystyle\nu\to 0.5$, without a phenomenological
imposition of incompressibility, and that the elastic moduli are in line with
typical polymer network gels. Further, despite the harmonic Gaussian chain as
a starting point, there is an emergent strain-softening and strain-stiffening
response that is characteristic of real polymer networks, driven by the
interplay between the excluded volume interactions and the entropy; it does
not require chains with limiting extensibility – such as the inverse Langevin
approximation – to model this behavior. Finally, we find the emergence of a
deformation-sensitive localization instability at large values of the excluded
volume parameter.

(a) Schematic of a polymer network.

(b) The 8-chain approximation.

Figure 2: The 8-chain approximation is obtained by averaging over a volume
element that aligns the polymer chains along the principal directions of the
deformation. Figure 3: The end-to-end vector in the undeformed state,
$\displaystyle\left({\mathbold X}_{\alpha}^{1}-{\mathbold
X}_{\alpha}^{0}\right)$, is mapped by the deformation $\displaystyle{\mathbold
F}$ to the end-to-end vector in the deformed state, $\displaystyle{\mathbold
x}_{\alpha}^{1}-{\mathbold x}_{\alpha}^{0}$, i.e., $\displaystyle{\mathbold
x}_{\alpha}^{1}-{\mathbold x}_{\alpha}^{0}={\mathbold F}\left({\mathbold
X}_{\alpha}^{1}-{\mathbold X}_{\alpha}^{0}\right)$.

##### Structure of the Paper.

Section 2 formulates the model. Section 3 summarizes the finite element
approach for the self-consistent solution. Section 4 presents numerical
results showing the predictions of the model.

##  2 Model Formulation

###  2.A Deformation of a Single Polymer Chain

Figure 4: Single polymer chain fixed at both ends.

We use the Continuous Gaussian Chain model for a single polymer chain [58]. In
the undeformed state, the coarse-grained trajectory of the
$\displaystyle\alpha$-th polymer chain is represented as a continuous 3-d
space curve $\displaystyle{\mathbold R}_{\alpha}(s)$, where $\displaystyle s$
is the chain contour coordinate and varies along the chain contour, and is
scaled such that $\displaystyle 0\leq s\leq 1$. The position vectors of the
beginning and end points of the chain in the undeformed state are
$\displaystyle{\mathbold X}_{\alpha}^{0}$ and $\displaystyle{\mathbold
X}_{\alpha}^{1}$.

The chain is deformed under the deformation gradient $\displaystyle{\mathbold
F}$. In the deformed state, $\displaystyle{\mathbold r}_{\alpha}(s)$ is a 3-d
curve that represents the coarse-grained trajectory of the
$\displaystyle\alpha$-th chain, as shown in Figure 4. The position vectors of
the beginning and end of the chain in the deformed state are
$\displaystyle{\mathbold x}_{\alpha}^{0}$ and $\displaystyle{\mathbold
x}_{\alpha}^{1}$.

Following [72], we use that the chain end-to-end vector is mapped under the
macroscale deformation $\displaystyle{\mathbold F}$:

${\mathbold x}_{\alpha}^{1}-{\mathbold x}_{\alpha}^{0}={\mathbold
F}\left({\mathbold X}_{\alpha}^{1}-{\mathbold X}_{\alpha}^{0}\right)$ (1)

We note that the affine deformation assumption depends strongly on the
assumption that there are no entanglements [73, 74, 75, 76].

####  2.A.1  Partition Function and Average Segment Density

Consider the $\displaystyle\alpha$-th chain that consists of $\displaystyle N$
coarse-grained polymer segments each of length $\displaystyle a$, and under
the influence of a field $\displaystyle w({\mathbold x})$ that will be used to
account for the excluded volume interactions [59].

From [58], the partition function, $\displaystyle Q_{\alpha}[w;{\mathbold
F}]$, and the average segment density,
$\displaystyle\langle\hat{\rho}_{\alpha}({\mathbold x};{\mathbold F})\rangle$,
are:

$\displaystyle\displaystyle Q_{\alpha}[w;{\mathbold F}]=\frac{1}{V}\int\
\mathrm{d}{\mathbold x}\ q_{\alpha}({\mathbold x},{\mathbold
x}_{\alpha}^{0},s)\ q_{\alpha}^{*}({\mathbold x},{\mathbold
x}_{\alpha}^{1},1-s),$ (2)
$\displaystyle\displaystyle\langle\hat{\rho}_{\alpha}({\mathbold x};{\mathbold
F})\rangle=\frac{1}{VQ_{\alpha}[w;{\mathbold F}]}\int\limits_{0}^{1}\
\mathrm{d}s\ q_{\alpha}({\mathbold x},{\mathbold x}_{\alpha}^{0},s)\
q_{\alpha}^{*}({\mathbold x},{\mathbold x}_{\alpha}^{1},1-s).$ (3)

Here, $\displaystyle q_{\alpha}({\mathbold x},{\mathbold x}_{\alpha}^{0},s)$
and $\displaystyle q_{\alpha}^{*}({\mathbold x},{\mathbold
x}_{\alpha}^{1},1-s)$ are the partial partition functions for the two chain
fragments, one from $\displaystyle 0$ to $\displaystyle s$ and the other from
$\displaystyle 1$ to $\displaystyle s$, respectively, as shown in Figure 4.

$\displaystyle q_{\alpha}({\mathbold x},{\mathbold x}_{\alpha}^{0},s)$ is
obtained by solving the following PDE with the initial condition:

$\frac{\partial q_{\alpha}({\mathbold x},{\mathbold
x}_{\alpha}^{0},s)}{\partial
s}=\frac{a^{2}N}{6}\nabla^{2}q_{\alpha}({\mathbold x},{\mathbold
x}_{\alpha}^{0},s)-w({\mathbold x})q_{\alpha}({\mathbold x},{\mathbold
x}_{\alpha}^{0},s),\quad q_{\alpha}({\mathbold x},{\mathbold
x}_{\alpha}^{0},s)\Big{|}_{s=0}=(aN^{1/2})^{3}\ \delta({\mathbold
x}-{\mathbold x}_{\alpha}^{0}).$ (4)

Similarly, $\displaystyle q_{\alpha}^{*}({\mathbold x},{\mathbold
x}_{\alpha}^{1},s^{\prime})$ is obtained by solving the same PDE as in (4),
but with the initial condition corresponding to keeping the other end fixed:

$\frac{\partial q_{\alpha}^{*}({\mathbold x},{\mathbold
x}_{\alpha}^{1},s^{\prime})}{\partial
s^{\prime}}=\frac{a^{2}N}{6}\nabla^{2}q_{\alpha}^{*}({\mathbold x},{\mathbold
x}_{\alpha}^{1},s^{\prime})-w({\mathbold x})q_{\alpha}^{*}({\mathbold
x},{\mathbold x}_{\alpha}^{1},s^{\prime}),\quad q_{\alpha}^{*}({\mathbold
x},{\mathbold
x}_{\alpha}^{1},s^{\prime})\Big{|}_{s^{\prime}=0}=(aN^{1/2})^{3}\
\delta({\mathbold x}-{\mathbold x}_{\alpha}^{1}).$ (5)

The initial conditions above correspond to the physical constraint that the
beginning and end points of the $\displaystyle\alpha$-th chain are fixed at
the given spatial positions $\displaystyle{\mathbold x}_{\alpha}^{0}$ and
$\displaystyle{\mathbold x}_{\alpha}^{1}$, respectively.

####  2.A.2 Reduction to Classical Rubber Elasticity

In the absence of excluded volume interactions, obtained by setting
$\displaystyle w({\mathbold x})\equiv 0$, we can find closed-form solutions
for $\displaystyle q_{\alpha}$ and $\displaystyle q_{\alpha}^{*}$:

$\displaystyle\displaystyle q_{\alpha}({\mathbold x},{\mathbold
x}_{\alpha}^{0},s)=\left(\frac{3}{2\pi
s}\right)^{3/2}\exp\left(-\frac{3|{\mathbold x}-{\mathbold
x}_{\alpha}^{0}|^{2}}{2a^{2}Ns}\right),$ (6) $\displaystyle\displaystyle
q_{\alpha}^{*}({\mathbold x},{\mathbold
x}_{\alpha}^{1},1-s)=\left(\frac{3}{2\pi(1-s)}\right)^{3/2}\exp\left(-\frac{3|{\mathbold
x}-{\mathbold x}_{\alpha}^{1}|^{2}}{2a^{2}N(1-s)}\right).$ (7)

The partition function $\displaystyle Q_{\alpha}[w;{\mathbold
F}]\Big{|}_{w=0}$ in (2) evaluates to the classical Gaussian distribution in
3-d:

$Q_{\alpha}[w;{\mathbold
F}]\Big{|}_{w=0}\propto\left(\frac{3}{\sqrt{\pi}}\right)^{3}\left(\frac{a^{2}N}{6}\right)^{3/2}\exp\left(-\frac{3}{2a^{2}N}|{\mathbold
x}_{\alpha}^{1}-{\mathbold x}_{\alpha}^{0}|^{2}\right).$ (8)

Because the chains do not interact, the free energy of the
$\displaystyle\alpha$-th polymer chain, $\displaystyle H_{\alpha}$, is
obtained from $\displaystyle Q_{\alpha}$ using $\displaystyle
H_{\alpha}=-k_{B}T\log Q_{\alpha}$ to be:

$H_{\alpha}[w;{\mathbold
F}]\Big{|}_{w=0}=\frac{1}{2}\left(\frac{3k_{B}T}{a^{2}N}\right)\left|{\mathbold
F}\left({\mathbold X}_{\alpha}^{1}-{\mathbold
X}_{\alpha}^{0}\right)\right|^{2}-\left(\frac{3k_{B}T}{2}\right).$ (9)

To account for the fact that chains are randomly oriented, we next average
$\displaystyle H_{\alpha}[w;{\mathbold F}]\Big{|}_{w=0}$ over all possible
orientations of the chain end-to-end vector by integrating (9) over all
orientations. That is, keeping $\displaystyle{\mathbold F}$ fixed, we
integrate $\displaystyle\left({\mathbold X}_{\alpha}^{1}-{\mathbold
X}_{\alpha}^{0}\right)$ over the sphere of appropriate radius. The resulting
expression for the orientationally-averaged free energy, $\displaystyle
H^{avg}_{\alpha}[w;{\mathbold F}]\Big{|}_{w=0}$, is:

$H^{avg}_{\alpha}[w;{\mathbold
F}]\Big{|}_{w=0}=\frac{k_{B}T}{2}\left(\operatorname{tr}({\mathbold
F}^{T}{\mathbold F})-3\right).$ (10)

This result recovers the classical rubber elasticity result [72], also known
as the incompressible Neo-Hookean elastic strain energy.

###  2.B Deformation of the Polymer Network

The pairwise excluded volume interactions are introduced through the field
$\displaystyle w({\mathbold x})$ following [59]. We introduce
$\displaystyle\bar{u}(|{\mathbold x}-{\mathbold x}^{\prime}|)$, which is the
pairwise interaction potential of mean force for two segments located at
spatial coordinates $\displaystyle{\mathbold x}$ and $\displaystyle{\mathbold
x}^{\prime}$. The corresponding partition function for the polymer network in
the deformed state, $\displaystyle Z({\mathbold F})$ in the field-theoretic
setting is:

$Z({\mathbold F})\propto\int\ \mathrm{D}\rho\int\ \mathrm{D}w\
\exp{\left(-\frac{H[w,\rho;{\mathbold F}]}{k_{B}T}\right)},$ (11)

where $\displaystyle H[w,\rho;{\mathbold F}]$ is the effective Hamiltonian of
the polymer network, and has the expression:

$\frac{H[w,\rho;{\mathbold F}]}{k_{B}T}=-\int\ \mathrm{d}{\mathbold x}\
w({\mathbold x})\rho({\mathbold x})+\ \frac{1}{2k_{B}T}\int\
\mathrm{d}{\mathbold x}\int\ \mathrm{d}{\mathbold x}^{\prime}\ \rho({\mathbold
x})\ \bar{u}(|{\mathbold x}-{\mathbold x}^{\prime}|)\ \rho({\mathbold
x}^{\prime})-\log\left(Q_{1}[w;{\mathbold F}]\ldots Q_{n}[w;{\mathbold
F}]\right).$ (12)

The auxiliary fields $\displaystyle w({\mathbold x})$ and
$\displaystyle\rho({\mathbold x})$ are interpreted as the fluctuating chemical
potential field generated internally because of the inter-segment interactions
and the fluctuating density of the polymer network, respectively [59].
$\displaystyle Q_{\alpha}[w;{\mathbold F}]$ is the partition function for the
$\displaystyle\alpha$-th chain in the polymer network under the influence of
$\displaystyle w({\mathbold x})$, and is calculated using (2).

The first term in (12) is the energy of interaction between the density and
the chemical potential. The second term is the inter-segment interaction
energy. The third term is the entropic contribution due to chain stretching.
The total Helmholtz free energy of polymer network in the deformed state,
$\displaystyle H({\mathbold F})$, is evaluated from the partition function
$\displaystyle Z({\mathbold F})$ in (11) using:

$H({\mathbold F})=-k_{B}T\log Z({\mathbold F}).$ (13)

In the deformed state, the average segment density of the polymer network,
$\displaystyle\langle\hat{\rho}({\mathbold x};{\mathbold F})\rangle$, is
obtained as:

$\begin{split}\langle\hat{\rho}({\mathbold x};{\mathbold
F})\rangle\propto&\frac{1}{Z({\mathbold F})}\int\ \mathrm{D}\rho\int\
\mathrm{D}w\ \Bigg{[}\exp\left(\int d{\mathbold x}\ w({\mathbold
x})\rho({\mathbold x})-\frac{1}{2k_{B}T}\int\ \mathrm{d}{\mathbold x}\int\
\mathrm{d}{\mathbold x}^{\prime}\rho({\mathbold x})\bar{u}(|{\mathbold
x}-{\mathbold x}^{\prime}|)\rho({\mathbold x}^{\prime})\right)\\\
&\left(\sum_{i=1}^{n}\left(\left(\int\limits_{0}^{1}\ \mathrm{d}s\
q_{i}({\mathbold x},{\mathbold x}_{i}^{0},s)\ q_{i}^{*}({\mathbold
x},{\mathbold x}_{i}^{1},1-s)\right)\ \prod\limits_{\begin{subarray}{c}k=1\\\
k\neq i\end{subarray}}^{n}Q_{k}[w;{\mathbold
F}]\right)\right)\Bigg{]}.\end{split}$ (14)

###  2.C Strain Energy Density of the Polymer Network

To obtain the continuum elastic response using nonlinear elasticity, we
introduce the elastic energy density (per undeformed unit volume)
$\displaystyle W({\mathbold F})$ and treat the RVE as a continuum material
point. This allows us to connect the total free energy $\displaystyle
H({\mathbold F})$ evaluated for the RVE to $\displaystyle W({\mathbold F})$ at
the corresponding spatial location:

$W({\mathbold F})=\frac{H({\mathbold F})}{V},$ (15)

where $\displaystyle V$ is the volume of the RVE in the undeformed state.

We can then find the Piola-Kirchhoff stress tensor $\displaystyle{\mathbold
P}$ and the fourth-order elasticity tensor
$\displaystyle\mathbb{C}=[C_{ijkl}]$ using:

$\displaystyle\displaystyle{\mathbold P}=\frac{\partial W}{\partial{\mathbold
F}},$ (16) $\displaystyle\displaystyle C_{ijkl}=\frac{\partial^{2}W}{\partial
F_{ij}\partial F_{kl}}\Bigg{|}_{{\mathbold F}={\mathbold I}},$ (17)

where $\displaystyle{\mathbold I}$ is the second-order identity tensor.

Applying (16), and substituting from (11), (12), (13), (15), we have:

$\begin{split}{\mathbold P}=&\frac{\partial W}{\partial{\mathbold
F}}=-\frac{k_{B}T}{V}\frac{1}{Z({\mathbold F})}\frac{\partial Z({\mathbold
F})}{\partial{\mathbold F}}=-\frac{k_{B}T}{V}\frac{1}{Z({\mathbold F})}\int\
\mathrm{D}\rho\int\ \mathrm{D}w\exp{\left(-\frac{H[w,\rho;{\mathbold
F}]}{k_{B}T}\right)}\frac{\partial H[w,\rho;{\mathbold F}]}{\partial{\mathbold
F}}\left(-\frac{1}{k_{B}T}\right)\\\ =&-\frac{k_{B}T}{V}\frac{1}{Z({\mathbold
F})}\int\ \mathrm{D}\rho\int\ \mathrm{D}w\exp{\left(-\frac{H[w,\rho;{\mathbold
F}]}{k_{B}T}\right)}\sum\limits_{\alpha=1}^{n}\frac{\partial}{\partial{\mathbold
F}}\log Q_{\alpha}[w;{\mathbold F}].\\\ \end{split}$ (18)

Then, defining the stress operator for the $\displaystyle\alpha$-th chain,
$\displaystyle\hat{{\mathbold P}}_{\alpha}$, as:

$\hat{{\mathbold P}}_{\alpha}:=-\frac{\partial\ }{\partial{\mathbold F}}\log
Q_{\alpha}[w;{\mathbold F}],$ (19)

we can write $\displaystyle{\mathbold P}$ as the statistical average [58,
Section 4.1.3] of $\displaystyle\hat{{\mathbold P}}_{\alpha}$:

${\mathbold P}=\frac{k_{B}T}{V}\left(\frac{1}{Z({\mathbold F})}\int\
\mathrm{D}\rho\int\ \mathrm{D}w\exp{\left(-\frac{H[w,\rho;{\mathbold
F}]}{k_{B}T}\right)}\sum\limits_{\alpha=1}^{n}\hat{{\mathbold
P}}_{\alpha}\right)=\frac{k_{B}T}{V}\sum\limits_{\alpha=1}^{n}\langle\hat{{\mathbold
P}}_{\alpha}\rangle.$ (20)

###  2.D Mean-field Assumption

The functional integration over the fields $\displaystyle w$ and
$\displaystyle\rho$ in (11) and (14) makes it expensive to evaluate
$\displaystyle H({\mathbold F})$ and
$\displaystyle\langle\hat{\rho}({\mathbold x};{\mathbold F})\rangle$.
Therefore, it is common to use a mean-field assumption to simplify the
functional integration in the expression for $\displaystyle Z({\mathbold F})$
in (11) [58]. This assumption implies that functional integration over the
fields $\displaystyle w$ and $\displaystyle\rho$ is dominated by the mean-
fields $\displaystyle\bar{w}$ and $\displaystyle\bar{\rho}$, respectively. The
mean-fields $\displaystyle\bar{w}$ and $\displaystyle\bar{\rho}$ are obtained
by requiring the effective Hamiltonian $\displaystyle H[w,\rho;{\mathbold F}]$
in (12) to be stationary with respect to variations in $\displaystyle
w({\mathbold x})$ and $\displaystyle\rho({\mathbold x})$. This gives the self-
consistent mean-field conditions:

$\displaystyle\displaystyle\frac{\delta H[w,\rho;{\mathbold F}]}{\delta
w}\Big{|}_{w=\bar{w}}=0,$ (21) $\displaystyle\displaystyle\frac{\delta
H[w,\rho;{\mathbold F}]}{\delta\rho}\Big{|}_{\rho=\bar{\rho}}=0.$ (22)

Using the mean-field assumption, $\displaystyle Z({\mathbold F})$ in (11)
simplifies to:

$Z({\mathbold F})\ \approx\exp\left(-\frac{H[\bar{w},\bar{\rho};{\mathbold
F}]}{k_{B}T}\right),$ (23)

where, $\displaystyle H[\bar{w},\bar{\rho};{\mathbold F}]$ is the effective
Hamiltonian in (12) evaluated using the mean-fields $\displaystyle\bar{w}$ and
$\displaystyle\bar{\rho}$. Using (13) and (23), the total free energy of
polymer network, $\displaystyle H({\mathbold F})$ under the mean-field
assumption is:

$\frac{H({\mathbold F})}{k_{B}T}=-\int\ \mathrm{d}{\mathbold x}\
\bar{w}({\mathbold x})\bar{\rho}({\mathbold x})\ +\ \frac{1}{2k_{B}T}\int\
\mathrm{d}{\mathbold x}\int\ \mathrm{d}{\mathbold x}^{\prime}\
\bar{\rho}({\mathbold x})\ \bar{u}(|{\mathbold x}-{\mathbold x}^{\prime}|)\
\bar{\rho}({\mathbold x}^{\prime})-\log\left(Q_{1}[\bar{w};{\mathbold
F}]\cdots Q_{n}[\bar{w};{\mathbold F}]\right).$ (24)

Further, the average segment density,
$\displaystyle\langle{\hat{\rho}}({\mathbold x};{\mathbold F})\rangle$, in
(14) simplifies to:

$\langle\hat{\rho}({\mathbold x};{\mathbold
F})\rangle\approx\sum\limits_{\alpha=1}^{n}\langle\hat{\rho}_{\alpha}({\mathbold
x};{\mathbold F})\rangle\Big{|}_{w=\bar{w}},$ (25)

where $\displaystyle\langle\hat{\rho}_{\alpha}({\mathbold x};{\mathbold
F})\rangle$ is the average segment density of the $\displaystyle\alpha$-th
chain in the polymer network, obtained using (3).

###  2.E Excluded Volume Interaction

Polymer segments in the network are considered to interact with each other
according to a pairwise interaction potential of mean force
$\displaystyle\bar{u}$ whose physical origin is due to excluded volume
effects. We account for the excluded volume effects by modeling a pairwise
inter-segment interaction using a simple Dirac delta potential of mean force
[77, 57]:

$\bar{u}(|{\mathbold x}-{\mathbold x}^{\prime}|)=k_{B}T\ u_{0}\
\delta(|{\mathbold x}-{\mathbold x}^{\prime}|).$ (26)

where $\displaystyle u_{0}$ is the excluded volume parameter. This form of
inter-segment interaction potential assumes the presence of a solvent in the
polymer network system with low density [49, 57]. The solvent mediates the
interactions among polymer segments. For $\displaystyle u_{0}>0$, implying
repulsion between the segments, the excluded volume potential
$\displaystyle\bar{u}$ in (26) is positive-definite and has an inverse;
following [58], this simplifies the field theory equations in (11) to:

$Z({\mathbold F})\propto\int\ \mathrm{D}w\exp\left(-\frac{H[w;{\mathbold
F}]}{k_{B}T}\right),$ (27)

where $\displaystyle H[w;{\mathbold F}]$ is the effective Hamiltonian of
polymer network in the simplified field theory:

$\frac{H[w;{\mathbold F}]}{k_{B}T}=\frac{1}{2u_{0}}\int\ \mathrm{d}{\mathbold
x}\ (w({\mathbold x}))^{2}-\log(Q_{1}[w;{\mathbold F}]\cdots
Q_{n}[w;{\mathbold F}]).$ (28)

Equations (27) and (28) present the simplified field theory for the
deformation of polymer network that is used in this work. The partition
function in (27) is evaluated using the mean-field assumption as:

$Z({\mathbold F})\approx\exp\left(-\frac{H[\bar{w};{\mathbold
F}]}{k_{B}T}\right),$ (29)

where $\displaystyle H[\bar{w};{\mathbold F}]$ is the effective Hamiltonian in
(28) evaluated using the mean-field $\displaystyle\bar{w}$. The mean-field
$\displaystyle\bar{w}$ is obtained by solving the stationarity condition for
the effective Hamiltonian $\displaystyle H[w;{\mathbold F}]$:

$\frac{\delta H[w;{\mathbold F}]}{\delta w}\Big{|}_{w=\bar{w}}=0.$ (30)

For the assumed form of the excluded volume interaction potential as in (26),
there is alternatively an expression for the average segment density [58]:

$\langle\hat{\rho}({\mathbold x};{\mathbold F})\rangle=\frac{1}{u_{0}}\langle
w({\mathbold x})\rangle=\frac{\bar{w}({\mathbold x})}{u_{0}}$ (31)

where $\displaystyle\langle w({\mathbold x})\rangle$ is the statistical
average of the fluctuating field $\displaystyle w({\mathbold x})$, and we use
that under the mean-field assumption $\displaystyle\langle w({\mathbold
x})\rangle=\bar{w}({\mathbold x})$. Finally, using (29) and (13), we obtain
the total free energy of polymer network, $\displaystyle H({\mathbold F})$ for
the simplified field theory as:

$\frac{H({\mathbold F})}{k_{B}T}=\frac{1}{2u_{0}}\int\ \mathrm{d}{\mathbold
x}(\bar{w}({\mathbold x}))^{2}-\log(Q_{1}[\bar{w};{\mathbold F}]\cdots
Q_{n}[\bar{w};{\mathbold F}]),$ (32)

where $\displaystyle\bar{w}({\mathbold x})$ is obtained by self-consistently
solving (31) and the mean-field condition in (30).

###  2.F Representative Volume Element Averaging: 8-Chain Model

A typical polymer network consists of a large number of cross-linked polymer
chains (Fig. 2) with random orientations at each continuum point, and is very
challenging to directly solve. To simplify this problem, we adopt the 8-chain
model for the RVE [40]. The 3-d RVE is assumed to be a cube in the undeformed
configuration, with 8 chains running between the center and each of the
corners. The cube is assumed to be oriented along the principal directions of
the macroscale right stretch tensor $\displaystyle{\mathbold U}$, where
$\displaystyle{\mathbold U}$ is the tensor square root of the deformation
$\displaystyle{\mathbold F}$ or alternatively is the positive-definite part of
the right polar decomposition of $\displaystyle{\mathbold F}$.

We assume that each chain begins ($\displaystyle s=0$) at the center of the
cube which is also taken to be the origin, and the chains terminate
($\displaystyle s=1$) at the corners. Denoting the terminating point of the
chains in the undeformed and deformed state, respectively, by
$\displaystyle{\mathbold X}_{1}^{1},\ldots{\mathbold X}_{8}^{1}$ and
$\displaystyle{\mathbold x}_{1}^{1},\ldots{\mathbold x}_{8}^{1}$, the relation
between the end-to-end vectors in the undeformed and deformed configurations
from (1) is:

${\mathbold x}_{\alpha}^{1}={\mathbold F}\ {\mathbold
X}_{\alpha}^{1},\qquad\alpha=1,\ldots 8.$ (33)

For a given value of $\displaystyle{\mathbold F}$, the right stretch tensor
$\displaystyle{\mathbold U}$ is used to orient the cube, and the equation
above provides the initial conditions for the partial partition functions
$\displaystyle q$ and $\displaystyle q^{*}$ of each chain in (4) and (5).

##  3 Numerical Method

Since the model with excluded volume interactions is not amenable to simple
closed-form solutions, we turn to numerical solutions. The goal is to evaluate
the total free energy $\displaystyle H({\mathbold F})$ and average segment
density $\displaystyle\langle\hat{\rho}({\mathbold x};{\mathbold F})\rangle$.
Our numerical method has the following steps:

  1. 1.

We first generate an initial field $\displaystyle w({\mathbold
x})=w^{0}({\mathbold x})$. The initial guess can be based on heuristics when
possible.

  2. 2.

We next solve for the average segment density
$\displaystyle\langle\hat{\rho}({\mathbold x};{\mathbold F})\rangle$ – using
$\displaystyle q$ and $\displaystyle q^{*}$ obtained by solving (4) and (5)
with the given $\displaystyle w({\mathbold x})$ – and the total free energy
$\displaystyle H({\mathbold F})$ using (30).

  3. 3.

The field $\displaystyle w$ is updated using (31) as $\displaystyle
w({\mathbold x})=u_{0}\langle\hat{\rho}({\mathbold x};{\mathbold F})\rangle$.

  4. 4.

In turn, we update $\displaystyle\langle\hat{\rho}({\mathbold x};{\mathbold
F})\rangle$ and $\displaystyle H({\mathbold F})$ as above.

This iteration continues until we reach convergence, which we define as a
relative change in the total free energy of less than $\displaystyle 0.1\%$.

For the solution of $\displaystyle q$ and $\displaystyle q^{*}$ in (4) and
(5), we use the finite element method (FEM) in the open-source FEniCS
framework [78]. The spatial domain is discretized using first-order Lagrange
family finite elements. The integration along the chain with respect to
$\displaystyle s$ in (4) and in (5) is performed using the implicit Crank-
Nicholson finite difference method with $\displaystyle 100$ steps for
$\displaystyle s\in(0,1)$. We test convergence of the FEM discretization as in
Figure 5(b).

While the RVE averaging nominally requires only 8 chains, these chains
interact not only with each other, but also with other chains that are not
contained in the RVE. To account for this, we use periodic images of the RVE;
we find that one image on each face of the cubic RVE is sufficient, giving us
27 cubes over which we must perform various integrations; Figure 5(a) shows a
schematic projection of this in 2-d, where the highlighted central RVE is used
for the energy calculations. When the deformation is applied, the image RVEs
are deformed following the central RVE. Similarly, $\displaystyle w({\mathbold
x})$ is defined over the larger cluster of RVEs for performing, but need only
be solved on the central RVE using periodicity.

(a)

(b)

Figure 5:  (a) Front view schematic of a bigger physical domain in 3-d that
consists of 27 RVEs, each with $\displaystyle 8$ chains. The central
$\displaystyle 8$ chain RVE (highlighted) is used for the free energy
calculation. (b) The convergence of the energy as the mesh is refined, shown
by plotting the free energy of a single chain without external field as a
function of mesh size. The converged mesh size of $\displaystyle 35^{3}$
elements is used for the numerical computations. Note that the saddle point
nature of the problem can lead to non-monotonic convergence.

##  4 Elastic Response with Excluded Volume Interactions

In the calculations reported here, we use the following model parameters:
total chain contour length $\displaystyle L=0.12\unit{\micro}$, number of
polymer segments in single chain $\displaystyle N=100$, excluded volume
parameter $\displaystyle u_{0}=0.005\ v_{seg}$ where $\displaystyle
v_{seg}=a^{3}$ is the volume of an individual monomer segment, and temperature
$\displaystyle T=303\unit{}$. This choice of $\displaystyle u_{0}$ gives an
excluded cube with side $\displaystyle 0.17a$, and corresponds to the Flory-
Huggins interaction parameter $\displaystyle\chi=0.4975$, using $\displaystyle
u_{0}=(1-2\chi)v_{seg}$. This characterizes a good solvent that is very close
to the $\displaystyle\Theta$-point [54, 58, 79, 80].

For the numerical calculations, larger values of $\displaystyle u_{0}$ in 3-d
require an excessively fine mesh to converge; however, the calculations
described below qualitatively agree with 2-d calculations – where much larger
values of $\displaystyle u_{0}$ can be used – and with a few representative
3-d calculations that were conducted with a larger value of $\displaystyle
u_{0}$ and a very fine mesh.

###  4.A Identifying the Stress-Free Free Energy Minimizing State

We first note that without excluded volume effects, all models based on the
Gaussian chain approximation would predict that the polymer network shrinks to
a point ($\displaystyle{\mathbold F}={\bf 0}$) because that would maximize the
configurational entropy of the chain. Therefore, as we add excluded volume
effects, we find that the equilibrium – stress-free or minimum energy – volume
of the polymer network RVE increases. Specifically, the equilibrium state is
achieved as a balance between the competing effects of entropic shrinkage and
excluded volume repulsion between monomers; [81, 82] examine related issues in
greater depth.

We denote the initial side of the cubic RVE by $\displaystyle
L_{uc}^{0}:=2aN^{1/2}/\sqrt{3}$, where $\displaystyle aN^{1/2}$ is the RMS
average diameter of an unconstrained chain, and the pre-factor accounts for
the geometry of the chain aligned along the diagonal of the RVE. However, we
emphasize that this is not the free energy minimizing state, i.e., it is not
stress-free, in the presence of excluded volume interactions. To this initial
state, we apply a deformation of the form $\displaystyle{\mathbold
F}=\lambda{\mathbold I}$, and compute the free energy for various values of
$\displaystyle\lambda$. Figure 6 shows the total free energy $\displaystyle H$
as a function of $\displaystyle\lambda=L_{uc}/L_{uc}^{0}$ for $\displaystyle
u_{0}=0.001\ v_{seg}$ and $\displaystyle u_{0}=0.005\ v_{seg}$. We find that
the stress-free stretches, i.e. the free-energy minimizing stretches, for
$\displaystyle u_{0}=0.001\ v_{seg}$ and $\displaystyle u_{0}=0.005\ v_{seg}$
are respectively at $\displaystyle\lambda=0.7$ and $\displaystyle\lambda=1$
respectively. The large increase in the equilibrium volume of the polymer
network system with an increase in the excluded volume parameter
$\displaystyle u_{0}$ is consistent with the phenomena of equilibrium swelling
for polymeric gels [83, 84, 85, 86].

In all of the subsequent calculations discussed below, we set the undeformed
state to correspond to the free energy minimizing stress-free state.

Figure 6: The total free energy as a function of the relative stretch for
$\displaystyle u_{0}=0.001\ v_{seg}$ and $\displaystyle u_{0}=0.005\ v_{seg}$.

###  4.B Volumetric and Shear Response, and Near-Incompressibility

We use the strain energy density $\displaystyle W({\mathbold F})$ to obtain
the mesoscale elastic response using nonlinear elasticity; specifically, (16)
and (17) are used to obtain the stress-stretch response and the elastic moduli
respectively. We assume below that the network can be treated as approximately
isotropic despite the 8-chain model.

To obtain the bulk and shear moduli, we impose deformations of the form:

${\mathbold F}_{v}=\begin{bmatrix}\lambda&0&0\\\ 0&\lambda&0\\\
0&0&\lambda\end{bmatrix}\quad,\quad{\mathbold
F}_{s}=\begin{bmatrix}1&\kappa_{s}&0\\\ 0&1&0\\\ 0&0&1\end{bmatrix}.$ (34)

The bulk and shear moduli, $\displaystyle K$ and $\displaystyle G$,
respectively can be computed using:

$K=\frac{1}{9}\frac{\partial^{2}W}{\partial\lambda^{2}}\Big{|}_{\lambda\approx
1}\quad,\quad
G=\frac{\partial^{2}W}{\partial\kappa_{s}^{2}}\Big{|}_{\kappa_{s}\approx 0}.$
(35)

We find $\displaystyle K=52.06\unit{\kilo}$ and $\displaystyle
G=0.60\unit{\kilo}$. Using isotropic linearized elasticity, this gives the
Poisson’s ratio $\displaystyle\nu=\frac{3K-2G}{6K+2G}=0.4943$ and the elastic
modulus $\displaystyle E=\frac{9KG}{3K+G}=1.79\unit{\kilo}$. These elastic
moduli are consistent with polymer network gels [87, 88, 22, 89, 90, 91, 92,
93, 94]. We highlight that $\displaystyle K$ is 2 orders of magnitude larger
than $\displaystyle G$, and $\displaystyle\nu$ is very close to the
incompressible limit of $\displaystyle 0.5$.

We next examine the shear stress vs. shear strain curve. In principle, the
shear stress can be computed using $\displaystyle\tau=\frac{\partial
W}{\partial\kappa_{s}}$. To avoid a lot of noise from numerical
differentiation, we fit $\displaystyle W$ by a polynomial and then
differentiate the polynomial to obtain the curve shown in Figure 7.

Figure 7: Shear stress $\displaystyle\tau$ vs. shear strain
$\displaystyle\kappa_{s}$ in simple shear. The segment density for the RVE at
various stretches are shown in the insets. We notice that the chains have
higher concentrations at the ends when the deformation is small, but are more
uniformly distributed as the deformation increases. Note that the RVE itself
does not appear to be sheared because the chain-averaging approach aligns the
averaging RVE along the principal directions that correspond to the maximum
and minimum elongation directions (Fig. 2).

###  4.C Extensional Response: Emergent Strain-Softening and Strain-Stiffening

We next examine extensional loading where the deformation has the form:

${\mathbold F}_{e}=\begin{bmatrix}\lambda_{1}&0&0\\\ 0&\lambda_{2}&0\\\
0&0&\lambda_{3}\end{bmatrix}$ (36)

Here, we consider $\displaystyle\lambda_{1}$ as the extensional stretch of
interest. We make 2 different choices for the transverse stretches
$\displaystyle\lambda_{2}$ and $\displaystyle\lambda_{3}$: the “constrained
case” where they are constrained such that
$\displaystyle\lambda_{2}=\lambda_{3}=1$; and the “volume-preserving case”
where they are set to be volume-preserving such that
$\displaystyle\lambda_{2}=\lambda_{3}=\lambda_{1}^{-\frac{1}{2}}$. Note that
the second case is approximately equivalent to having no transverse stress. To
obtain the extensional stress, we use $\displaystyle\sigma=\frac{\partial
W}{\partial\lambda_{1}}$. Figure 8 compares the stress-strain response of
these cases. We notice that the constrained case has significantly higher
stresses and tangent moduli. Figure 9 shows the evolution of the chain density
with stretch for the constrained case.

Figure 8: The stress-strain curves for the constrained and volume-preserving
extensional loadings. We notice that the stresses and the tangent moduli are
both significantly larger when the system is constrained to undergo
deformations that do not preserve the volume. Figure 9: Extensional stress
$\displaystyle\sigma$ vs. extensional stretch $\displaystyle\lambda_{t}$ for
the constrained case from Figure 8. The segment density for the RVE at various
stretches are shown in the insets. As with shearing, we notice that the chains
have higher concentrations at the ends when the deformation is small, but are
more uniformly distributed as the deformation increases.

From Figures 8 and 9, we notice that both cases show a pronounced strain-
softening and strain-stiffening behavior that is characteristic of many real
polymer networks such as polymeric gels. However, Gaussian chains do not show
such behavior, and it is typical to use chains with limiting extensibility –
such as the inverse Langevin approximation – to model this behavior. Here, we
find that it is a consequence of the competition between the excluded volume
parameter and the entropy. Figure 10 shows the decomposition of the free
energy $\displaystyle W$ into entropic $\displaystyle W_{entropic}$ and
excluded volume interaction $\displaystyle W_{interaction}$ contributions. We
observe that the excluded volume contribution is less than the entropic
contribution in both cases. Further, we notice that $\displaystyle
W_{entropic}$ monotonically increases with $\displaystyle\lambda_{1}>1$, and
is consistent with the stretching of the Gaussian polymer chains; however,
$\displaystyle W_{interaction}$ monotonically decreases with
$\displaystyle\lambda_{1}>1$, and is consistent with the chains being more
oriented, and hence having fewer excluded volume interactions. We notice that
in the approximate range $\displaystyle 1.5<\lambda_{1}<2.5$ where we see
strain softening, the decrease in the excluded volume interaction is faster
than the rise in the entropic contribution, causing softening. For
$\displaystyle\lambda_{1}>3$, we have the opposite trend in that the entropy
increases faster than the decrease in the excluded volume interaction, causing
stiffening. In summary, strain softening occurs because of the initial
decrease in excluded volume interactions, and subsequent strain stiffening
occurs because of the later increase in entropic effects.

(a) Free energy density $\displaystyle W$ vs stretch
$\displaystyle\lambda_{1}$, with $\displaystyle\lambda_{2}=\lambda_{3}=1$

(b) Free energy density $\displaystyle W$ vs stretch
$\displaystyle\lambda_{1}$ with
$\displaystyle\lambda_{2}=\lambda_{3}=\lambda_{1}^{-1/2}$

Figure 10: The free energy $\displaystyle W$ is decomposed into entropic
$\displaystyle W_{entropic}$ and excluded volume interaction $\displaystyle
W_{interaction}$ contributions for the constrained and volume-preserving
cases. The symbols show the simulation results, and the lines show best fits
that are differentiated to obtain stress-strain curves.

###  4.D Effect of Chain Length

Figure 11 shows the effect of chain contour length on the elastic moduli of
the polymer network. The chain contour length $\displaystyle L$ is varied from
$\displaystyle 0.01\unit{\micro}$ to $\displaystyle 0.3\unit{\micro}$ while
keeping $\displaystyle N$ fixed. We observe that both the elastic modulus
$\displaystyle E$ and the shear modulus $\displaystyle G$ decrease with
increased chain contour length. An increase in chain contour length
corresponds to an increase in the average molecular weight ($\displaystyle
M_{c}$) between the cross-links, and these results are consistent with
experiments that show that an increase in $\displaystyle M_{c}$ corresponds to
a decrease in the elastic moduli [95, 96]. The range of elastic and shear
moduli obtained using the model by varying chain contour length is consistent
with the experimental values for polymer network soft matter such as
elastomers and polymeric gels [97, 98, 99, 100, 101, 102, 87, 88, 22, 89, 90,
91, 92, 93, 94].

Figure 11: Elastic and shear moduli as a function of chain contour length; as
expected from polymer theory, these scale as $\displaystyle a^{-3}$ [103].

###  4.E Interactions between Deformation and an Excluded Volume Instability

We next examine an instability driven by the increase in the excluded volume
parameter. For computational feasability – because we aim to numerically
confirm that the instability is sharp by a large number of calculations near
the instability – we focus on 2-d systems; however, a few representative
calculations suggest that 3-d is qualitatively similar. Because it is in 2-d,
the undeformed RVE is a square with 4 chains.

Figure 12 presents the average segment density field for various equi-biaxial
stretches and various excluded volume parameter values. For a fixed RVE
stretch of $\displaystyle L_{uc}/L_{uc}^{0}=1$, where $\displaystyle
L^{0}_{uc}=\sqrt{2}\ aN^{1/2}$ in 2-d, we observe an instability at
$\displaystyle u_{0}\approx 0.7\ v_{seg}$ $\displaystyle(v_{seg}=a^{2}$ in
2-d), leading to the localization of chains. Physically, the chains strongly
repel each other and hence are highly restricted in the volume available. The
instability is symmetry breaking, in that the originally square-symmetric
chain configuration transitions to localize either vertically or horizontally
from the original square symmetry; in our numerical simulations, we find that
these occur essentially randomly due to numerical noise. As noted above, the
instability is a sharp transition.

We examine the effect of an imposed equi-biaxial stretch by setting
$\displaystyle L_{uc}/L_{uc}^{0}$ to $\displaystyle 1.5$ and $\displaystyle 2$
respectively. We notice that the critical values of $\displaystyle u_{0}$ for
the instability are, respectively, $\displaystyle u_{0}\approx 0.8\ v_{seg}$
and $\displaystyle u_{0}\approx 1.2\ v_{seg}$. This coupling between the
deformation and chain localization suggests new routes to obtain patterning in
polymer networks.

Figure 12:  Excluded volume-driven instability observed in 2-d, and the effect
of equi-biaxial deformation. Each subfigure shows the average segment density
plotted over 9 RVEs, with each subfigure corresponding to different values of
biaxial stretch $\displaystyle L_{uc}/L_{uc}^{0}$and $\displaystyle u_{0}$.
The instability corresponds to a sharp transition in the chain configuration:
it goes from being fairly uniform away from the crosslinking point to being
concentrated along the horizontal and vertical directions.

As the stretch increases, the critical value of the excluded volume parameter
at which the instability occurs also increases.

##  5 Discussion

We have used the statistical field theory of polymers in combination with the
8-chain network averaging approach to study the mechanical response of polymer
networks. The framework of polymer field theory provides a physics-based
approach to accounting for excluded volume interactions, which are imposed
phenomenologically in micromechanical models. In the absence of excluded
volume interactions, we find that that the closed-form orientationally-
averaged elastic response matches with the classical rubber elasticity [72].
With excluded volume effects, self-consistent numerical solutions using finite
elements find that the predicted elastic moduli are in line with typical
polymer network gels; particularly, the linearized Poisson’s ratio
$\displaystyle\nu\simeq 0.4943$, which is very close to the incompressible
limit $\displaystyle\nu\to 0.5$, without a phenomenological imposition of
incompressibility. Though the equilibrium state depends on the value of
$\displaystyle u_{0}$, the incompressible behavior is independent of the
specific value of $\displaystyle u_{0}$ for the values studied here. This can
be physically understood by considering that $\displaystyle\nu$ is computed
around the equilibrium stress-free state. Due to entropic effects, the chains
tend to reduce their end-to-end distance and would collapse to a point, while
the excluded volume effects prevent the chains from collapsing completely. The
equilibrium state is achieved when these opposing effects balance out. Around
this equilibrium state, we find $\displaystyle\nu\simeq 0.5$, which reflects
the role of the solvent in preventing further reduction in volume. Despite the
seeming presence of voids or open space in the density fields, the chains are
constrained due to the excluded volume interactions; the voids will be
occupied by the solvent which makes the system close to incompressible since
the solvent cannot leave the polymer network upon deformation. This is
consistent with the results of [104], wherein it was found that
$\displaystyle\nu\simeq 0.5$ at short times when the solvent has not had time
to diffuse out of the polymer network; at longer times, the solvent diffuses
out and the long-time equilibrium value of $\displaystyle\nu$ depends on the
shear modulus. While we have not considered this time-dependent behavior here,
it is related to similar effects in poromechanics, namely the Terzaghi and
Mandel model problems, wherein the stress response is closely tied to the
drainage of the pore fluid [105, 106, 107, 108, 109].

Another interesting finding is that, despite the harmonic Gaussian chain as a
starting point, there is an emergent strain-softening and strain-stiffening
response that is characteristic of real polymer network gels, driven by the
interplay between the excluded volume interactions and the entropy; it does
not require chains with limiting extensibility – such as the inverse Langevin
approximation – to model this behavior. We also we find the emergence of a
deformation-sensitive localization instability at large values of the excluded
volume parameter.

A natural question for future is to examine the interplay between chain-scale
instabilities such as microbuckling [110, 111, 112] and the network-scale
instabilities observed here. We highlight, however, that these examples of
instabilities are in fibrous networks, and it is possible that such
instabilities occur here because of the regular network structure that has
been assumed and will not appear in random networks. An important challenge,
however, is that the isolated harmonic Gaussian chain does not display
buckling or other instabilities; other nonlinear chain models are required to
capture this behavior. Further, as noted in [8], electrical field interactions
provide an effective compressive stiffness, and can induce new types of
instabilities [12]. Incorporating chain models that go beyond the Gaussian
approximation in polymer field theory is an interesting theoretical question.
Along similar lines, while we capture strain stiffening and strain softening
without the use of a chain model with limiting extensibility – such as the
inverse Langevin model – it would be interesting to incorporate such models in
polymer field theory to enable studying the interplay between entropic,
excluded volume, and limited extensibility effects.

The concept of chain topology or entanglement is an important aspect that is
not taken into consideration in this study. As highlighted in [73, 74], these
effects can play a significant role in the response of polymer networks. The
mean-field framework employed in this study is unable to account for such
effects directly. However, our inclusion of excluded volume effects provides
some insight into the effects of entanglement. Additionally, we have observed
that even without accounting for entanglements, excluded volume effects give
rise to many interesting physical characteristics that are relevant to the
response of real polymer networks. It should be noted that while entanglements
are crucial for large deformation response, the linearized properties such as
the Poisson’s ratio are expected to be relatively unaltered.

## Code Availability

A version of the code developed for this work is available at
https://github.com/pkhandag/polymer-network.git

###### Acknowledgements.

We thank Carlos Garcia Cervera for useful discussions; NSF XSEDE for computing
resources provided by the Pittsburgh Supercomputing Center; AFRL for hosting
visits by Kaushik Dayal; NSF (DMREF 1921857, DMS 2108784), ONR
(N00014-18-1-2528), BSF (2018183), and AFOSR (MURI FA9550-18-1-0095) for
financial support; and the anonymous reviewers for comments that improved the
paper significantly.

## References

  * [1] Qian Deng, Liping Liu, and Pradeep Sharma.  Electrets in soft materials: Nonlinearity, size effects, and giant electromechanical coupling.  Physical Review E, 90(1):012603, 2014.
  * [2] Lingling Chen, Xu Yang, Binglei Wang, Shengyou Yang, Kaushik Dayal, and Pradeep Sharma.  The interplay between symmetry-breaking and symmetry-preserving bifurcations in soft dielectric films and the emergence of giant electro-actuation.  Extreme Mechanics Letters, 43:101151, 2021.
  * [3] Fatemeh Ahmadpoor and Pradeep Sharma.  Flexoelectricity in two-dimensional crystalline and biological membranes.  Nanoscale, 7(40):16555–16570, 2015.
  * [4] Oleg V Kim, Xiaojun Liang, Rustem I Litvinov, John W Weisel, Mark S Alber, and Prashant K Purohit.  Foam-like compression behavior of fibrin networks.  Biomechanics and modeling in mechanobiology, 15(1):213–228, 2016\.
  * [5] André EX Brown, Rustem I Litvinov, Dennis E Discher, Prashant K Purohit, and John W Weisel.  Multiscale mechanics of fibrin polymer: gel stretching with protein unfolding and loss of water.  science, 325(5941):741–744, 2009.
  * [6] Tianxiang Su and Prashant K Purohit.  Semiflexible filament networks viewed as fluctuating beam-frames.  Soft Matter, 8(17):4664–4674, 2012.
  * [7] Matthew Grasinger and Kaushik Dayal.  Architected elastomer networks for optimal electromechanical response.  Journal of the Mechanics and Physics of Solids, 146:104171, 2021\.
  * [8] Matthew Grasinger, Carmel Majidi, and Kaushik Dayal.  Nonlinear statistical mechanics drives intrinsic electrostriction and volumetric torque in polymer networks.  Physical Review E, 103(4):042504, 2021.
  * [9] Matthew Grasinger and Kaushik Dayal.  Statistical mechanical analysis of the electromechanical coupling in an electrically-responsive polymer chain.  Soft Matter, 16(27):6265–6284, 2020.
  * [10] Noy Cohen, Kaushik Dayal, and Gal deBotton.  Electroelasticity of polymer networks.  Journal of the Mechanics and Physics of Solids, 92:105–126, 2016\.
  * [11] Faezeh Darbaniyan, Kaushik Dayal, Liping Liu, and Pradeep Sharma.  Designing soft pyroelectric and electrocaloric materials using electrets.  Soft matter, 15(2):262–277, 2019.
  * [12] Matthew Grasinger, Kosar Mozaffari, and Pradeep Sharma.  Flexoelectricity in soft elastomers and the molecular mechanisms underpinning the design and emergence of giant flexoelectricity.  Proceedings of the National Academy of Sciences, 118(21), 2021.
  * [13] Matthew Grasinger.  Group action markov chain monte carlo for accelerated sampling of energy landscapes with discrete symmetries and energy barriers.  arXiv preprint arXiv:2205.00028, 2022.
  * [14] Eric J Markvicka, Michael D Bartlett, Xiaonan Huang, and Carmel Majidi.  An autonomously electrically self-healing liquid metal–elastomer composite for robust soft-matter robotics and electronics.  Nature materials, 17(7):618–624, 2018.
  * [15] Michael D Bartlett, Navid Kazem, Matthew J Powell-Palm, Xiaonan Huang, Wenhuan Sun, Jonathan A Malen, and Carmel Majidi.  High thermal conductivity in soft elastomers with elongated liquid metal inclusions.  Proceedings of the National Academy of Sciences, 114(9):2143–2148, 2017.
  * [16] Navid Kazem, Tess Hellebrekers, and Carmel Majidi.  Soft multifunctional composites and emulsions with liquid metals.  Advanced Materials, 29(27):1605985, 2017.
  * [17] Carmel Majidi.  Soft-matter engineering for soft robotics.  Advanced Materials Technologies, 4(2):1800477, 2019.
  * [18] Yongyi Zhao, Pratik Khandagale, and Carmel Majidi.  Modeling electromechanical coupling of liquid metal embedded elastomers while accounting stochasticity in 3d percolation.  Extreme Mechanics Letters, 48:101443, 2021.
  * [19] Michael D Bartlett, Michael D Dickey, and Carmel Majidi.  Self-healing materials for soft-matter machines and electronics.  NPG Asia Materials, 11(1):1–4, 2019.
  * [20] Michael J Ford, Cedric P Ambulo, Teresa A Kent, Eric J Markvicka, Chengfeng Pan, Jonathan Malen, Taylor H Ware, and Carmel Majidi.  A multifunctional shape-morphing elastomer with liquid metal inclusions.  Proceedings of the National Academy of Sciences, 116(43):21438–21444, 2019.
  * [21] Navid Zolfaghari, Pratik Khandagale, Michael J Ford, Kaushik Dayal, and Carmel Majidi.  Network topologies dictate electromechanical coupling in liquid metal–elastomer composites.  Soft Matter, 16(38):8818–8825, 2020.
  * [22] Yunsik Ohm, Chengfeng Pan, Michael J Ford, Xiaonan Huang, Jiahe Liao, and Carmel Majidi.  An electrically conductive silver–polyacrylamide–alginate hydrogel composite for soft electronics.  Nature Electronics, 4(3):185–192, 2021.
  * [23] Mohammad H Malakooti, Michael R Bockstaller, Krzysztof Matyjaszewski, and Carmel Majidi.  Liquid metal nanocomposites.  Nanoscale Advances, 2(7):2668–2677, 2020.
  * [24] Taylor H Ware, John S Biggins, Andreas F Shick, Mark Warner, and Timothy J White.  Localized soft elasticity in liquid crystal elastomers.  Nature communications, 7(1):1–7, 2016.
  * [25] Timothy J White and Dirk J Broer.  Programmable and adaptive mechanics with liquid crystal polymer networks and elastomers.  Nature materials, 14(11):1087–1098, 2015.
  * [26] Taylor H Ware, Michael E McConney, Jeong Jae Wie, Vincent P Tondiglia, and Timothy J White.  Voxelated liquid crystal elastomers.  Science, 347(6225):982–984, 2015.
  * [27] Cedric P Ambulo, Julia J Burroughs, Jennifer M Boothby, Hyun Kim, M Ravi Shankar, and Taylor H Ware.  Four-dimensional printing of liquid crystal elastomers.  ACS applied materials & interfaces, 9(42):37332–37339, 2017.
  * [28] Jiuke Mu, Mônica Jung De Andrade, Shaoli Fang, Xuemin Wang, Enlai Gao, Na Li, Shi Hyeong Kim, Hongzhi Wang, Chengyi Hou, Qinghong Zhang, et al.  Sheath-run artificial muscles.  Science, 365(6449):150–155, 2019.
  * [29] Mohand O Saed, Cedric P Ambulo, Hyun Kim, Rohit De, Vyom Raval, Kyle Searles, Danyal A Siddiqui, John Michael O Cue, Mihaela C Stefan, M Ravi Shankar, et al.  Molecularly-engineered, 4d-printed liquid crystal elastomer actuators.  Advanced Functional Materials, 29(3):1806412, 2019.
  * [30] Jeong Jae Wie, M Ravi Shankar, and Timothy J White.  Photomotility of polymers.  Nature communications, 7(1):1–8, 2016.
  * [31] Mahnoush Babaei, J Arul Clement, Kaushik Dayal, and M Ravi Shankar.  Steering with light: indexable photomotility in liquid crystalline polymers.  RSC advances, 7(83):52510–52516, 2017.
  * [32] Mahnoush Babaei, Junfeng Gao, Arul Clement, Kaushik Dayal, and M Ravi Shankar.  Torque-dense photomechanical actuation.  Soft Matter, 17(5):1258–1266, 2021.
  * [33] Melvin Mooney.  A theory of large elastic deformation.  Journal of applied physics, 11(9):582–592, 1940.
  * [34] Raymond W Ogden.  Non-linear elastic deformations.  Courier Corporation, 1997.
  * [35] Alan N Gent.  A new constitutive relation for rubber.  Rubber chemistry and technology, 69(1):59–61, 1996.
  * [36] Hubert M James and Eugene Guth.  Theory of the elastic properties of rubber.  The Journal of Chemical Physics, 11(10):455–481, 1943.
  * [37] Paul J Flory and John Rehner Jr.  Statistical mechanics of cross-linked polymer networks i. rubberlike elasticity.  The journal of chemical physics, 11(11):512–520, 1943.
  * [38] LRG Treloar.  The elasticity of a network of long-chain molecules.—iii.  Transactions of the Faraday Society, 42:83–94, 1946.
  * [39] LRG Treloar.  The photoelastic properties of short-chain molecular networks.  Transactions of the Faraday Society, 50:881–896, 1954.
  * [40] Ellen M Arruda and Mary C Boyce.  A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials.  Journal of the Mechanics and Physics of Solids, 41(2):389–412, 1993\.
  * [41] PD Wu and Erik Van Der Giessen.  On improved network models for rubber elasticity and their applications to orientation hardening in glassy polymers.  Journal of the Mechanics and Physics of Solids, 41(3):427–456, 1993\.
  * [42] C Miehe, Serdar Göktepe, and F Lulei.  A micro-macro approach to rubber-like materials—part i: the non-affine micro-sphere model of rubber elasticity.  Journal of the Mechanics and Physics of Solids, 52(11):2617–2660, 2004.
  * [43] Matthew Grasinger.  Polymer networks which locally rotate to accommodate stresses, torques, and deformation.  Journal of the Mechanics and Physics of Solids, 175:105289, 2023\.
  * [44] G Ronca and G Allegra.  An approach to rubber elasticity with internal constraints.  The Journal of Chemical Physics, 63(11):4990–4997, 1975.
  * [45] Paul J Flory.  Statistical thermodynamics of random networks.  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 351(1666):351–380, 1976.
  * [46] PJ Flory.  Theory of elasticity of polymer networks. the effect of local constraints on junctions.  The Journal of Chemical Physics, 66(12):5720–5729, 1977.
  * [47] Burak Erman and Paul J Flory.  Theory of elasticity of polymer networks. ii. the effect of geometric constraints on junctions.  The Journal of Chemical Physics, 68(12):5363–5369, 1978.
  * [48] Paul J Flory and Burak Erman.  Theory of elasticity of polymer networks. 3.  Macromolecules, 15(3):800–806, 1982.
  * [49] RT Deam and Samuel Frederick Edwards.  The theory of rubber elasticity.  Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 280(1296):317–353, 1976.
  * [50] SF Edwards and Th A Vilgis.  The tube model theory of rubber elasticity.  Reports on Progress in Physics, 51(2):243, 1988.
  * [51] G Heinrich and E Straube.  On the strength and deformation dependence of the tube-like topological constraints of polymer networks, melts and concentrated solutions. i. the polymer network case.  Acta polymerica, 34(9):589–594, 1983.
  * [52] G Heinrich and E Straube.  On the strength and deformation dependence of the tube-like topological constraints of polymer networks, melts and concentrated solutions. ii. polymer melts and concentrated solutions.  Acta polymerica, 35(2):115–119, 1984.
  * [53] G Heinrich, E Straube, and G Helmis.  Rubber elasticity of polymer networks: Theories.  Polymer physics, pages 33–87, 1988.
  * [54] P-G De Gennes.  Some conformation problems for long macromolecules.  Reports on Progress in Physics, 32(1):187, 1969.
  * [55] Pierre-Gilles De Gennes.  Scaling concepts in polymer physics.  Cornell university press, 1979.
  * [56] Mark W Matsen.  Self-consistent field theory and its applications.  Soft Matter, 1, 2006.
  * [57] Masao Doi and Samuel Frederick Edwards.  The theory of polymer dynamics.  oxford university press, 1986.
  * [58] Glenn Fredrickson.  The equilibrium theory of inhomogeneous polymers.  Oxford University Press, 2006.
  * [59] Glenn H Fredrickson, Venkat Ganesan, and François Drolet.  Field-theoretic computer simulation methods for polymers and complex fluids.  Macromolecules, 35(1):16–39, 2002.
  * [60] Tanya L Chantawansri, August W Bosse, Alexander Hexemer, Hector D Ceniceros, Carlos J Garcia-Cervera, Edward J Kramer, and Glenn H Fredrickson.  Self-consistent field theory simulations of block copolymer assembly on a sphere.  Physical Review E, 75(3):031802, 2007.
  * [61] Erin M Lennon, Kirill Katsov, and Glenn H Fredrickson.  Free energy evaluation in field-theoretic polymer simulations.  Physical review letters, 101(13):138302, 2008.
  * [62] François Drolet and Glenn H Fredrickson.  Combinatorial screening of complex block copolymer assembly with self-consistent field theory.  Physical Review Letters, 83(21):4317, 1999.
  * [63] Scott W Sides, Bumjoon J Kim, Edward J Kramer, and Glenn H Fredrickson.  Hybrid particle-field simulations of polymer nanocomposites.  Physical review letters, 96(25):250601, 2006.
  * [64] David M Ackerman, Kris Delaney, Glenn H Fredrickson, and Baskar Ganapathysubramanian.  A finite element approach to self-consistent field theory calculations of multiblock polymers.  Journal of Computational Physics, 331:280–296, 2017.
  * [65] Glenn H Fredrickson.  Computational field theory of polymers: opportunities and challenges.  Soft Matter, 3(11):1329–1334, 2007.
  * [66] Kris T Delaney and Glenn H Fredrickson.  Recent developments in fully fluctuating field-theoretic simulations of polymer melts and solutions.  The Journal of Physical Chemistry B, 120(31):7615–7634, 2016.
  * [67] Eric W Cochran, Carlos J Garcia-Cervera, and Glenn H Fredrickson.  Stability of the gyroid phase in diblock copolymers at strong segregation.  Macromolecules, 39(7):2449–2451, 2006.
  * [68] Erin M Lennon, George O Mohler, Hector D Ceniceros, Carlos J García-Cervera, and Glenn H Fredrickson.  Numerical solutions of the complex langevin equations in polymer field theory.  Multiscale Modeling & Simulation, 6(4):1347–1370, 2008.
  * [69] Mark W Matsen and Michael Schick.  Stable and unstable phases of a diblock copolymer melt.  Physical Review Letters, 72(16):2660, 1994.
  * [70] Mark W Matsen.  The standard gaussian model for block copolymer melts.  Journal of Physics: Condensed Matter, 14(2):R21, 2001.
  * [71] Friederike Schmid.  Self-consistent field approach for cross-linked copolymer materials.  Physical Review Letters, 111(2):028303, 2013.
  * [72] Leslie Ronald George Treloar.  The physics of rubber elasticity.  Oxford University Press, USA, 1975.
  * [73] Jacob D Davidson and Nakhiah C Goulbourne.  A nonaffine network model for elastomers undergoing finite deformations.  Journal of the Mechanics and Physics of Solids, 61(8):1784–1797, 2013.
  * [74] Jacob D Davidson and NC Goulbourne.  Nonaffine chain and primitive path deformation in crosslinked polymers.  Modelling and Simulation in Materials Science and Engineering, 24(6):065002, 2016.
  * [75] Michael Rubinstein and Sergei Panyukov.  Nonaffine deformation and elasticity of polymer networks.  Macromolecules, 30(25):8036–8044, 1997.
  * [76] Michael Rubinstein and Sergei Panyukov.  Elasticity of polymer networks.  Macromolecules, 35(17):6670–6686, 2002.
  * [77] BeHe Zimm, WH Stockmayer, and M Fixman.  Excluded volume in polymer chains.  The Journal of Chemical Physics, 21(10):1716–1723, 1953.
  * [78] Hans Petter Langtangen and Anders Logg.  Solving PDEs in python: the FEniCS tutorial I.  Springer Nature, 2017.
  * [79] Scott T Milner, Martin-Daniel Lacasse, and William W Graessley.  Why $\displaystyle\chi$ is seldom zero for polymer- solvent mixtures.  Macromolecules, 42(3):876–886, 2009.
  * [80] Janhavi Nistane, Lihua Chen, Youngjoo Lee, Ryan Lively, and Rampi Ramprasad.  Estimation of the flory-huggins interaction parameter of polymer-solvent mixtures using machine learning.  MRS Communications, pages 1–7, 2022.
  * [81] JP Wittmer, Philippe Beckrich, Hendrik Meyer, Anna Cavallo, Albert Johner, and Joerg Baschnagel.  Intramolecular long-range correlations in polymer melts: The segmental size distribution and its moments.  Physical Review E, 76(1):011803, 2007.
  * [82] Michael Lang, Michael Rubinstein, and Jens-Uwe Sommer.  Conformations of a long polymer in a melt of shorter chains: Generalizations of the flory theorem.  ACS Macro Letters, 4(2):177–181, 2015.
  * [83] Wei Hong, Xuanhe Zhao, Jinxiong Zhou, and Zhigang Suo.  A theory of coupled diffusion and large deformation in polymeric gels.  Journal of the Mechanics and Physics of Solids, 56(5):1779–1793, 2008.
  * [84] Hiroyuki Kamata, Xiang Li, Ung-il Chung, and Takamasa Sakai.  Design of hydrogels for biomedical applications.  Advanced healthcare materials, 4(16):2360–2374, 2015.
  * [85] Toyoichi Tanaka and David J Fillmore.  Kinetics of swelling of gels.  The Journal of Chemical Physics, 70(3):1214–1218, 1979.
  * [86] Joseph Jagur-Grodzinski.  Polymeric gels and hydrogels for biomedical and pharmaceutical applications.  Polymers for Advanced Technologies, 21(1):27–47, 2010.
  * [87] Michael DA Norman, Silvia A Ferreira, Geraldine M Jowett, Laurent Bozec, and Eileen Gentleman.  Measuring the elastic modulus of soft culture surfaces and three-dimensional hydrogels using atomic force microscopy.  Nature Protocols, 16(5):2418–2449, 2021.
  * [88] Vivian R Feig, Helen Tran, Minah Lee, and Zhenan Bao.  Mechanically tunable conductive interpenetrating network hydrogels that mimic the elastic moduli of biological tissue.  Nature communications, 9(1):2740, 2018.
  * [89] Shuai Xu, Zidi Zhou, Zishun Liu, and Pradeep Sharma.  Concurrent stiffening and softening in hydrogels under dehydration.  Science Advances, 9(1):eade3240, 2023.
  * [90] Stephen Coyle, Carmel Majidi, Philip LeDuc, and K Jimmy Hsia.  Bio-inspired soft robotics: Material selection, actuation, and design.  Extreme Mechanics Letters, 22:51–59, 2018.
  * [91] Jianyu Li, Zhigang Suo, and Joost J Vlassak.  Stiff, strong, and tough hydrogels with good chemical stability.  Journal of Materials Chemistry B, 2(39):6708–6713, 2014.
  * [92] Chuanpeng Sun and Prashant K Purohit.  Rheology of fibrous gels under compression.  Extreme Mechanics Letters, 54:101757, 2022.
  * [93] Shengqiang Cai, Derek Breid, Alfred J Crosby, Zhigang Suo, and John W Hutchinson.  Periodic patterns and energy states of buckled films on compliant substrates.  Journal of the Mechanics and Physics of Solids, 59(5):1094–1114, 2011.
  * [94] Wanliang Shan, Stuart Diller, Abbas Tutcuoglu, and Carmel Majidi.  Rigidity-tuning conductive elastomer.  Smart Materials and Structures, 24(6):065001, 2015.
  * [95] Lale G Lovell and Christopher N Bowman.  The effect of kinetic chain length on the mechanical relaxation of crosslinked photopolymers.  Polymer, 44(1):39–47, 2003.
  * [96] Emilia I Wisotzki, Marcel Hennes, Carsten Schuldt, Florian Engert, Wolfgang Knolle, Ulrich Decker, Josef A Käs, Mareike Zink, and Stefan G Mayr.  Tailoring the material properties of gelatin hydrogels by high energy electron irradiation.  Journal of Materials Chemistry B, 2(27):4297–4309, 2014.
  * [97] Yongyi Zhao, Yunsik Ohm, Jiahe Liao, Yichi Luo, Huai-Yu Cheng, Phillip Won, Peter Roberts, Manuel Reis Carneiro, Mohammad F Islam, Jung Hyun Ahn, et al.  A self-healing electrically conductive organogel composite.  Nature Electronics, pages 1–10, 2023.
  * [98] Brad A Krajina, Audrey Zhu, Sarah C Heilshorn, and Andrew J Spakowitz.  Active dna olympic hydrogels driven by topoisomerase activity.  Physical review letters, 121(14):148001, 2018.
  * [99] Jianyu Li, Widusha RK Illeperuma, Zhigang Suo, and Joost J Vlassak.  Hybrid hydrogels with extremely high stiffness and toughness.  ACS Macro Letters, 3(6):520–523, 2014.
  * [100] Ziqian Li, Zishun Liu, Teng Yong Ng, and Pradeep Sharma.  The effect of water content on the elastic modulus and fracture energy of hydrogel.  Extreme Mechanics Letters, 35:100617, 2020.
  * [101] Max C Darnell, Jeong-Yun Sun, Manav Mehta, Christopher Johnson, Praveen R Arany, Zhigang Suo, and David J Mooney.  Performance and biocompatibility of extremely tough alginate/polyacrylamide hydrogels.  Biomaterials, 34(33):8042–8048, 2013.
  * [102] Dhriti Nepal, Saewon Kang, Katarina M Adstedt, Krishan Kanhaiya, Michael R Bockstaller, L Catherine Brinson, Markus J Buehler, Peter V Coveney, Kaushik Dayal, Jaafar A El-Awady, et al.  Hierarchically structured bioinspired nanocomposites.  Nature materials, 22(1):18–35, 2023.
  * [103] Michael Rubinstein and Ralph H Colby.  Polymer physics, volume 23.  Oxford university press, 2003.
  * [104] Robyn H Pritchard and Eugene M Terentjev.  Swelling and de-swelling of gels under external elastic deformation.  Polymer, 54(26):6954–6960, 2013.
  * [105] Maurice A Biot.  Nonlinear and semilinear rheology of porous solids.  Journal of Geophysical Research, 78(23):4924–4937, 1973.
  * [106] Olivier Coussy.  Poromechanics.  John Wiley & Sons, 2004.
  * [107] JW Rudnicki.  Coupled deformation-diffusion effects in the mechanics of faulting and failure of geomaterials.  Appl. Mech. Rev., 54(6):483–502, 2001.
  * [108] Mina Karimi, Mehrdad Massoudi, Noel Walkington, Matteo Pozzi, and Kaushik Dayal.  Energetic formulation of large-deformation poroelasticity.  International Journal for Numerical and Analytical Methods in Geomechanics, 46(5):910–932, 2022.
  * [109] Herbert Wang.  Theory of linear poroelasticity with applications to geomechanics and hydrogeology, volume 2.  Princeton university press, 2000.
  * [110] Georgios Grekas, Maria Proestaki, Phoebus Rosakis, Jacob Notbohm, Charalambos Makridakis, and Guruswami Ravichandran.  Cells exploit a phase transition to mechanically remodel the fibrous extracellular matrix.  Journal of the Royal Society Interface, 18(175):20200823, 2021.
  * [111] R Lakes, P Rosakis, and A Ruina.  Microbuckling instability in elastomeric cellular solids.  Journal of materials science, 28(17):4667–4672, 1993.
  * [112] Chuanpeng Sun, Irina N Chernysh, John W Weisel, and Prashant K Purohit.  Fibrous gels modelled as fluid-filled continua with double-well energy landscape.  Proceedings of the Royal Society A, 476(2244):20200643, 2020.