††thanks: Present address: Google Cloud AI, Pittsburgh, PA, 15206, USA

# Prestressed elasticity of amorphous solids

Shang Zhang  Department of Physics, University of Michigan, Ann Arbor, MI
48109-1040, USA  Ethan Stanifer  Department of Physics, University of
Michigan, Ann Arbor, MI 48109-1040, USA  Vishwas Vasisht  Department of
Physics, Indian Institute of Technology, Palakkad, India  Leyou Zhang
Department of Physics, University of Michigan, Ann Arbor, MI 48109-1040, USA
Emanuela Del Gado  Department of Physics, Georgetown University, Washington
DC, USA  Xiaoming Mao<EMAIL_ADDRESS>Department of Physics, University of
Michigan, Ann Arbor, MI 48109-1040, USA

###### Abstract

Prestress in amorphous solids bears the memory of their formation, and plays a
profound role in their mechanical properties, from stiffening or softening
elastic moduli to shifting frequencies of vibrational modes, as well as
directing yielding and solidification in the nonlinear regime. Here we develop
a set of mathematical tools to investigate elasticity of prestressed discrete
networks, which disentangles the effects from disorder in configuration and
disorder in prestress. Applying these methods to prestressed triangular
lattices and a computational model of amorphous solids, we demonstrate the
importance of prestress on elasticity, and reveal a number of intriguing
effects caused by prestress, including strong spatial heterogeneity in stress-
response unique to prestressed solids, power-law distribution of minimal
dipole stiffness, and a new criterion to classify floppy modes.

##  I Introduction

Almost all solid materials are stressed. Amorphous solids exhibit quenched
residual stress from their preparation process, crystalline solids are
stressed by defects and grain-boundaries, and living matter experiences active
stress from biological processes. The ubiquity of stressed solids is
encapsulated in the fact that the stress tensor $\sigma$, as a $d\times d$
symmetrical matrix in $d$ dimensions, has $d(d+1)/2$ independent degrees of
freedom, but the force balance equation,

$\partial_{j}\sigma_{ij}=0,$ (1)

only poses $d$ constraints, leaving $d(d-1)/2$ unconstrained components in the
stress field. In absence of external load, all components have to vanish if
the material, and the stress field, have to be homogeneous, but they can be
nonzero at lengthscales over which heterogeneities or excess constraints are
present (“geometric frustration”), giving rise to prestress, (also known as
“residual stress”, “initial stress” or “eigenstress” in different contexts)
Alexander (1998); Ulm and Coussy (2003).

In structural engineering, prestress is proactively used to modify both
stability and load bearing capability of structures, from prestressed
reinforced concrete to tensegrity architectures [Fig. 1 (a), (h) and (j)]. In
materials, instead, prestresses can emerge spontaneously, as the direct
consequence of out-of-equilibrium processes through which they solidify, or of
the external load applied during processing. Prominent examples include
isotropic compressive prestress in jammed packings Wyart  _et al._ (2005) or
gels Bouzid  _et al._ (2017), shear prestress in shear jammed granular matter
Bi  _et al._ (2011) and shear thickened dense suspensions Sehgal  _et al._
(2019), isotropic prestresses in glasses Ballauff  _et al._ (2013), and rich
varieties of anisotropic prestress fields in prestressed/tensegrity
metamaterials Gei  _et al._ (2009); Chen  _et al._ (2012); Fraternali  _et
al._ (2014); Meeussen  _et al._ (2020); Merrigan  _et al._ (2021) and
biological systems Armon  _et al._ (2011); Wyczalkowski  _et al._ (2012); Feng
_et al._ (2018) [Fig.1 (b)-(g) and (i)]. Remarkably, very much in the same way
as for buildings and large scale structures, microscopic residual stresses in
amorphous solids may strongly affect their strength—stiffen or soften them,
and direct how they fail Alexander (1998); Zhang  _et al._ (2017); Ozawa  _et
al._ (2018); Vasisht  _et al._ (2020); Barlow  _et al._ (2020); Berthier  _et
al._ (2019). In addition, prestress stores elastic energy in materials, which
then feature a different energy landscape when compared to stress-free
materials Debenedetti and Stillinger (2001); Manning and Liu (2011); Ninarello
_et al._ (2017); Vasisht  _et al._ (2021); Bhaumik  _et al._ (2021); Richard
_et al._ (2020).

Figure 1: A series of examples of stressed materials including structures
under tension [a tensile net architecture in (a) and a model colloidal gel
Bouzid  _et al._ (2017)) in (b)], structures under compression [contact
network of a jammed packing of repulsive disks Bouchaud (2002) (c) and pile of
sand (d)], structures under shear [shear jamming of photoelastic discs (e) and
an acoustically tunable shear thickening fluid (image taken by Neil Lin) (f)],
structures under competing tensile and compressive forces [a Prince Rupert’s
drop Aben  _et al._ (2016) in (g) and a bridge made from prestressed concrete
in (h)], as well as systems with more complex modes of stress [a pea pod which
uses stress to open Armon  _et al._ (2011) in (i) and a tensegrity
metamaterial Wen  _et al._ (2020) in (j)].

A deep understanding of the consequences of prestress is therefore key to
predict elasticity and material properties in general, but microscopic
prestresses can be very elusive, since it is difficult to directly access them
in experiments and in most cases only their indirect consequences can be
detected. It has been shown that, in jamming of frictionless particles,
compressive prestress both stabilizes the packing by maintaining particle
contacts, and destabilizes it by shifting frequencies of vibrational modes
Wyart  _et al._ (2005); DeGiuli  _et al._ (2014), as well as control the
response to external forces Ellenbroek  _et al._ (2009). In biological systems
such as biopolymer gels and epithelial cell sheets, tensile stress has been
shown to provide stability Licup  _et al._ (2015); Feng  _et al._ (2015,
2016); Noll  _et al._ (2017); Hardin  _et al._ (2018); Yan and Bi (2019);
Merkel  _et al._ (2019); Liu  _et al._ (2021a), and active stress generated by
molecular motors has been shown to give rise to sophisticated effects on
rigidity Alvarado  _et al._ (2013) and mechanical signal transmission Ronceray
_et al._ (2016).

However, a general theoretical framework to disentangle effects of
microstructure and prestress, and map how external load is carried by
prestressed amorphous solids, is still missing. Many studies of elasticity of
amorphous solids have mainly focused on determining how the disorder in their
microstructure affects elasticity Tanguy  _et al._ (2002a); O’Hern  _et al._
(2003); Zaccone and Scossa-Romano (2011); Liu and Nagel (2010), whereas the
role of prestress is far less understood, so that there is a gap in the
theoretical approaches that can rationalize the mechanics of amorphous solids
when it comes to prestress. In a granular solid, the same configuration of
particles and the same value of the stress at the boundary are compatible with
multiple stress distributions, which correspond to the so-called “force
network ensemble” Snoeijer  _et al._ (2004); Tighe  _et al._ (2008, 2010); Bi
_et al._ (2015) and is related to the broader concept of the “Edwards
ensemble” Edwards and Oakeshott (1989) where both microscopic configuration
and stress are allowed to vary for a given boundary stress. Statistical
relations of prestress distributions in this type of ensembles have been
recently discussed, and force-balance constraints [Eq. (1)] have been shown to
produce intriguing long-range correlations Mao  _et al._ (2007); Henkes and
Chakraborty (2009); Mao  _et al._ (2009); Lemaître (2017); Maier  _et al._
(2017); DeGiuli (2018); Nampoothiri  _et al._ (2020). Nevertheless, unraveling
the role of disordered configurations and disordered prestress in the
mechanical properties of amorphous solids has remained difficult both
conceptually and computationally.

When considering the stress distributions in a mechanically stable structure,
the concept of states of self stress (SSSs) Calladine (1978); Pellegrino and
Calladine (1986) captures their eigenstates that leave all components of the
solid structure in force balance. In discrete mechanical networks, SSSs
originate from redundant constraints. SSSs were first introduced in mechanical
engineering to describe load bearing abilities of structures, and more
recently used by the condensed matter physics community for their role in
characterizing mechanical topological edge states Sun  _et al._ (2012); Kane
and Lubensky (2014); Lubensky  _et al._ (2015); Mao and Lubensky (2018).
Because SSSs span the linear space of all possible ways a system can carry
load, mapping out, and further programming, this SSSs linear space offers a
convenient handle to control the mechanical response of materials. Interesting
examples of this type include recent work on programming SSSs using
topological mechanics to direct buckling and fracturing of metamaterials
Paulose  _et al._ (2015); Zhang and Mao (2018). Intriguingly, SSSs in jammed
packings have been shown to exhibit rich spatial structures and a divergent
length scale near jamming Sussman  _et al._ (2016). Some of these findings
suggest that SSSs may be the natural tool to rigorously quantify the effect of
prestress distributions on material properties, but a direct connection
remains to be built.

Here we present a systematic method to investigate how prestress affects the
mechanical response of amorphous solids, based on the concept of SSSs, where
the configuration and the prestress are allowed to change _independently_ ,
while force-balance is always maintained. As we discuss below, prestress is a
linear combination of the SSSs of the stress-free version of the mechanical
network, and leads to a new set of SSSs, capturing the unique, prestress-
controlled, mechanical response of the system. By applying this method to a
prestressed triangular lattice model and a computational model of amorphous
solids, we show that prestress greatly enlarges the linear space of stress-
bearing states, and that a number of intriguing effects arise where prestress
strongly controls the mechanical response, both at the local and the global
level.

The general method we introduce here is applicable to a wide range of systems,
ordered and disordered, where prestress affects elasticity. It provides an
efficient computational algorithm for finding the stress distribution when the
system is under any load, as well as offering a platform to develop a field-
theoretic treatment of prestressed elasticity of amorphous solids. Because
this method allows the microstructure and the prestress field to vary
separately, it offers a pathway to investigate how amorphous solids evolve as
well as develop memory under stress without changing their configuration. We
demonstrate the method using a triangular lattice model with varying
prestress, and test it in amorphous configurations of compressed repulsive
particles, obtained through numerical simulations. In such model solids, we
show how prestress determines their response to both macroscopic shear strain
and local dipole forces, and show that they display qualitatively different
behaviors from unstressed spring networks with the same geometry. Finally, for
the model solids obtained through computer simulations, we use the new method
to study the dependence of the stress-bearing ability on the preparation
protocol, which changes the microscopic prestress distribution, as well as
various signatures of the spatial evolution of the stress under a shear
strain.

Overall, the new method proposed here is an ideal candidate to investigate the
nature of rigidity transitions in prestressed systems, and furthermore,
yielding and shear thickening/thinning systems such as granular matter or
jammed packings and dense suspensions, to potentially shed new light on the
dynamical interplay between stress and geometry in these complex materials.

The paper is organized as follows. In Sec. II we review the fundamental
theories of prestressed mechanical networks. In Sec. III we introduce a new
formulation of equilibrium and compatibility matrices in prestressed networks,
and discuss how they control both zero modes and stress response. In Sec. IV
we use a prestressed triangular lattice model to illustrate how prestress
controls the vibrational modes and the stress response. In Sec. V we apply
these methods to analyze the stress response of a computational model of
prestressed amorphous solids, and Sec. VI contains a summary and outlook.

##  II Prestressed mechanics

In this section we briefly review mechanics of prestressed networks, their
continuum elasticity limit, as well as how prestress affects rigidity of
discrete networks in general.

###  II.1 Prestressed mechanical networks

We consider a discrete network of point-like particles connected by pairwise
central-force potentials (bonds). When this network is deformed, particle
$\ell$, which was originally at $\vec{R}_{\ell,0}$ in the reference state,
undergoes a displacement $\vec{u}_{\ell}$ to a new position

$\vec{R}_{\ell}=\vec{R}_{\ell,0}+\vec{u}_{\ell}.$ (2)

The change of the elastic energy $V_{b}$ associated with this deformation, of
a bond $b$ connecting particles $\ell,\ell^{\prime}$, is $\delta
V_{b}=V_{b}(|\vec{R}_{b}|)-V_{b}(|\vec{R}_{b,0}|)$, where
$\vec{R}_{b,0}\equiv\vec{R}_{\ell,0}-\vec{R}_{\ell^{\prime},0}$,
$\vec{R}_{b}\equiv\vec{R}_{\ell}-\vec{R}_{\ell^{\prime}}$ are the vectors
along the bond in the reference state and the deformed state respectively.

This change in energy can be expanded to quadratic order in $u$ [using Eq.
(2)]

$\delta
V_{b}=\frac{V_{b}^{\prime\prime}}{2}|e_{b}^{\parallel}|^{2}+\frac{V_{b}^{\prime}}{2|\vec{R}_{b,0}|}|e_{b}^{\perp}|^{2},$
(3)

where the derivatives are taken at $|\vec{R}_{b,0}|$, and
$\vec{e}_{b}\equiv\vec{u}_{\ell}-\vec{u}_{\ell^{\prime}}$ is the difference of
displacement vectors of the two particles connected by bond $b$, and

$e_{b}^{\parallel}\equiv\hat{R}_{b,0}\hat{R}_{b,0}\cdot\vec{e}_{b},\quad
e_{b}^{\perp}\equiv(\mathbb{I}-\hat{R}_{b,0}\hat{R}_{b,0})\cdot\vec{e}_{b}$
(4)

are its components parallel and perpendicular to the original bond direction
$\hat{R}_{b,0}\equiv\vec{R}_{b,0}/|\vec{R}_{b,0}|$, respectively.

The derivatives of this potential,

$\displaystyle
V_{b}^{\prime\prime}(|\vec{R}_{b}|)\equiv\frac{d^{2}V_{b}(|\vec{R}_{b}|)}{d|\vec{R}_{b}|^{2}}=k_{b},$
$\displaystyle
V_{b}^{\prime}(|\vec{R}_{b}|)\equiv\frac{dV_{b}(|\vec{R}_{b}|)}{d|\vec{R}_{b}|}=t_{b,p},$
(5)

correspond to the spring constant $k_{b}$ and the pretension $t_{b,p}$
(tension of the bond in the reference state, where $p$ stands for “prestress”)
of bond $b$.

In a stress-free network, all bonds are at their rest length $R_{b,R}$ and
thus $V_{b}^{\prime}(|\vec{R}_{b,0}|)=0$, leaving only the
$V_{b}^{\prime\prime}$ term in Eq. (3). Vibrational modes of disordered
networks of this type have been extensively studied, yielding a rich set of
interesting phenomena including quasilocalized modes, anomalies of density of
states at low frequencies, etc. Tanguy  _et al._ (2002a); O’Hern  _et al._
(2003); Wyart  _et al._ (2005, 2008); Liu and Nagel (2010); Chen  _et al._
(2011); Lerner  _et al._ (2012); Mosayebi  _et al._ (2014); Stanifer  _et al._
(2018).

When a network is prestressed, $V_{b}^{\prime}(|\vec{R}_{b,0}|)\neq 0$, and
both terms contribute to the elastic energy, increasing the number of
microscopic constraints, as we discuss in detail in Sec. III. Interestingly,
in general the sign of $t_{b,p}$ can be either positive (tension) or negative
(compression). In the case of $t_{b,p}>0$, the prestress term contribute with
another complete square term to the elastic energy, clearly stabilizing the
system. In the case of $t_{b,p}<0$, naively, the prestress term appears to be
unstable. However, because $e_{b}^{\perp}$ are not variables independent of
$e_{b}^{\parallel}$ (of other bonds), the stability of the network needs to be
analyzed in terms of the collective modes, which we discuss more in Sec. II.3.

Here we assumed that the network in the reference state is in _force balance_
, i.e., the total force on each particle vanishes. As a result, there is no
$\mathcal{O}(e)$ term in the expansion of the elastic energy. Note that _the
force-balance condition and the stress-free condition are two distinct
conditions_ , where the later means $V_{b}^{\prime}(|\vec{R}_{b,0}|)=0$ on all
bonds, and is a much more stringent requirement than the force-balance
condition. We will revisit this distinction in continuum elasticity in Sec.
II.2.

###  II.2 Continuum elasticity with prestress

The discrete theory discussed above can be rigorously linked to continuum
elasticity using the relation between the discrete nonlinear strain $v_{b}$ of
bond $b$ defined as Born  _et al._ (1955); DiDonna and Lubensky (2005)

$v_{b}\equiv\frac{1}{2}\left(|\vec{R}_{b}|^{2}-|\vec{R}_{b,0}|^{2}\right)$ (6)

and the (continuum) nonlinear right Cauchy-Green strain tensor (repeated
indices are summed over)

$\epsilon_{ij}(\vec{x})\equiv\frac{1}{2}\left(\partial_{i}u_{j}+\partial_{j}u_{i}+\partial_{i}u_{l}\partial_{i}u_{l}\right)$
(7)

where the relation reads

$v_{b}=\vec{R}_{b,0}\cdot\epsilon_{ij}(\vec{x})\cdot\vec{R}_{b,0}.$ (8)

This follows from the definition of the nonlinear strain tensor, and we have
taken the continuum limit by assuming that the spatial variation of the strain
field is small over the scale of the particles and bonds, so that the
deformation of bond $b$ is determined by the strain at its location $\vec{x}$.

By expanding the bond length $|\vec{R}_{b}|$ in terms of $v_{b}$, we can
rewrite the expansion of the elastic energy of bond $b$ as

$\delta
V_{b}=V_{b}^{\prime}\frac{v_{b}}{|\vec{R}_{b,0}|}+\frac{1}{2}\left(V_{b}^{\prime\prime}-\frac{V_{b}^{\prime}}{|\vec{R}_{b,0}|}\right)\left(\frac{v_{b}}{|\vec{R}_{b,0}|}\right)^{2}.$
(9)

This is equivalent to the expansion in Eq. (3), by recognizing that

$v_{b}=\vec{R}_{b,0}\cdot\vec{e}_{b}+\frac{1}{2}\vec{e}_{b}\cdot\vec{e}_{b}=|\vec{R}_{b,0}|e_{b}^{\parallel}+\frac{(e_{b}^{\parallel})^{2}+(e_{b}^{\perp})^{2}}{2}.$
(10)

The combination of Eqs. (8) and (9) allows us to write the bond energy change
in terms of the strain tensor.

The elastic energy of the whole system can be taken to the continuum limit by
converting the sum over all bonds to an integral over space

$E=\frac{1}{2}\sum_{\ell}\sum_{\ell^{\prime}}V_{b=\langle\ell,\ell^{\prime}\rangle}=\frac{1}{2}\int
d^{d}\vec{x}\frac{1}{v(\vec{x})}\sum_{\ell^{\prime}}V_{b=\langle\ell,\ell^{\prime}\rangle}$
(11)

where the space is Voronoi tessellated according to the particles, and
$v(\vec{x})$ is the volume of the Voronoi cell at $\vec{x}$, and the overall
factor of $1/2$ comes from the fact that every bond is counted twice. The sum
$\sum_{\ell^{\prime}}$ in the continuum limit represent the sum over the
bonded neighbors of the particle at $\vec{x}$. Using this formulation, the
elastic energy of the system can be written in the conventional form

$E=\int
d^{d}\vec{x}\left[\frac{1}{2}K_{ijkl}(\vec{x})\epsilon_{ij}(\vec{x})\epsilon_{ij}(\vec{x})+\sigma_{p,ij}(\vec{x})\epsilon_{ij}(\vec{x})\right]$
(12)

where the local elastic-modulus tensor $K_{ijkl}(\vec{x})$ and prestress field
$\sigma_{p,ij}(\vec{x})$ are determined from the discrete network by

$\displaystyle K_{ijkl}(\vec{x})=$
$\displaystyle\frac{1}{2v(\vec{x})}\sum_{\ell^{\prime}}\left(V_{b}^{\prime\prime}|\vec{R}_{b,0}|^{2}-V_{b}^{\prime}|\vec{R}_{b,0}|\right)$
$\displaystyle\quad\cdot\hat{R}_{b,0,i}\hat{R}_{b,0,j}\hat{R}_{b,0,k}\hat{R}_{b,0,l},$
$\displaystyle\sigma_{p,ij}(\vec{x})=$
$\displaystyle\frac{1}{2v(\vec{x})}\sum_{\ell^{\prime}}(V_{b}^{\prime}|\vec{R}_{b,0}|)\hat{R}_{b,0,i}\hat{R}_{b,0,j},$
(13)

where the sum $\ell^{\prime}$ is over all connected neighbors of particle
$\ell$, which is at position $\vec{x}$, and
$b=\langle\ell,\ell^{\prime}\rangle$. From this relation, it is clear that the
prestress field $\sigma_{p}(\vec{x})$ in the continuum theory comes from
pretensions on bonds which are not at their rest length in the discrete
network, whereas the elastic-modulus tensor depends on both the spring
constants and the tensions.

The body-force in this continuum theory

$\displaystyle
f_{i}(\vec{x})=\partial_{j}\sigma_{p,ij}(\vec{x})=\frac{1}{2v(\vec{x})}\sum_{\ell^{\prime}}(V_{b}^{\prime}|\vec{R}_{b,0}|)\hat{R}_{b,0,i}$
(14)

corresponds to the total force on each particle (normalized by the volume of
the Voronoi cell) in the discrete network. Thus, the force-balance conditions
in the discrete network and the continuum theory are indeed the same
condition.

It is often useful to write this continuum theory in a quadratic expansion in
terms of the displacement field $\vec{u}(\vec{x})$,

$\displaystyle E=$ $\displaystyle\int
d^{d}\vec{x}\Big{[}\frac{1}{2}K_{ijkl}(\vec{x})\partial_{i}u_{j}(\vec{x})\partial_{j}u_{i}(\vec{x})$
$\displaystyle+\sigma_{p,ij}(\vec{x})\partial_{i}u_{l}(\vec{x})\partial_{j}u_{l}(\vec{x})\Big{]}$
(15)

where we have used the symmetry of $K_{ijkl}$ and the force balance condition
which eliminates $\mathcal{O}(u)$ terms in this elastic energy.

###  II.3 Prestressed rigidity

The concept of “rigidity” has been a central theme in the discussion of
mechanics of soft materials. In general, rigidity can be interpreted from two
viewpoints, (i) a microscopic one, where rigidity is attributed to the
vanishing of floppy modes (i.e. modes of deformation that cost no elastic
energy) Maxwell (1864); Tay and Whiteley (1984), and (ii) a macroscopic one,
where rigidity is defined as the emergence of a spanning rigid cluster that
can transmit stress Jacobs and Thorpe (1995). Remarkably, these two rigidity
criteria coincide for the jamming transition of frictionless repulsive
particles Liu and Nagel (2010).

In this section we focus on the first viewpoint and discuss how prestress
changes floppy modes. In Sec. III we discuss how prestress affects the ability
of the system to carry additional loads (second point of view) with the help
of the concept of SSSs.

Rigidity of stress-free discrete mechanical networks is often analyzed via the
comparison between the numbers of constraints and degrees of freedom Maxwell
(1864); Jacobs and Thorpe (1995). This is well summarized by the Maxwell-
Callandine index theorem Calladine (1978); Sun  _et al._ (2012). Rigorously
speaking, when this theorem determines that the number of floppy modes of a
system vanishes, it implies that the system is “first-order rigid”, where the
energy increases quadratically for all deformations. By contrast, in a
“second-order rigid” system the leading order expansion of the energy with
respect to some deformations is of higher order such as cubic or quartic
Connelly and Whiteley (1992, 1996).

Figure 2:  Prestressed rigidity. (a) A mechanical network of concentric
similar triangles Connelly and Whiteley (1992). When stress-free, this network
has one floppy mode (black dashed arrows). When prestressed, this floppy mode
is eliminated. Blue (red) bonds are under tension (compression). The thickness
of the bond is proportional to the magnitude of the prestress. (b) A honeycomb
lattice exhibits an extensive number of floppy modes when stress-free (one
example shown as dashed arrows). When prestressed with equal tension on all
bonds (blue), these floppy modes are eliminated. (c) Lowest phonon band of the
honeycomb lattice when stress-free (orange) and when stretched with tensile
prestress as shown in (b) (blue). (d) Eigenvalues $\lambda$ of the dynamical
matrix of a mechanical network (right). When prestress on the network increase
beyond a threshold, one normal mode (e) becomes unstable (mode number $6$ as
shown by the red arrow), demonstrating an example of type B ZMs discussed in
Sec. III.3 [black dashed arrows in (e)].

How does prestress change this paradigm of rigidity? It was pointed out in
Ref. Connelly and Whiteley (1992, 1996) that prestress can rigidify second-
order rigid systems and make all their modes first-order rigid. Interestingly,
this prestress stabilization effect can even take place in systems which are
underconstrained according to the (stress-free) Maxwell counting. One simple
example is the honeycomb lattice, which has coordination number $z=3$ and is
thus significantly below the Maxwell point of $z=2d$ ($d$ being the spatial
dimension) where degrees of freedom and constraints balance. When a honeycomb
lattice is under tension, all modes are lifted to finite frequency (except for
trivial translations). It is worth noting that this example does not violate
the theorems in Ref. Connelly and Whiteley (1996): although underconstrained,
floppy modes of the honeycomb lattice under periodic boundary conditions are
in fact second order rigid. In Fig. 2 we show examples of how networks become
rigid by prestress.

This effect has been studied in various disordered mechanical networks from
polymer gels to jamming of particles and biological tissues Wyart  _et al._
(2008); Licup  _et al._ (2015); Feng  _et al._ (2016); Vermeulen  _et al._
(2017); Merkel  _et al._ (2019); Damavandi  _et al._ (2021). In these models,
the underconstrained network typically exhibit an extensive number of floppy
modes before the application of stress. When certain types of strain (tensile
or shear) is applied to the network, a system spanning SSS emerges and becomes
stressed, providing (macroscopic) rigidity and thus elastic moduli. This
prestress also rigidifies floppy modes in the network, as the geometry of the
network evolves and these modes become second order rigid.

Interestingly, by increasing the magnitude of the prestress, a stable network
can also destabilize. This happens when the negative terms in Eq. (3), due to
compressed bonds, compete with the positive terms and creates negative
eigenvalues in the dynamical matrix. Some examples of this effect are also
shown in Fig. 2. As we introduce our new mathematical tools in Sec. III, we
discuss this effect in more detail.

##  III States of self-stress in stressed elasticity

Now we introduce a new $\mathbb{Q},\mathbb{C}$ (equilibrium and compatibility)
decomposition for the mechanics of prestressed systems, that can be used to
analyze their SSSs and zero modes (ZMs). We first briefly review this
decomposition for stress-free systems, and then discuss the more general case
of prestressed systems. We also develop a general formulation of these SSSs to
compute mechanical response to an external load, including both macroscopic
strain and local forces.

###  III.1 States of self-stress in stress-free systems

The elastic energy of a stress-free network of $N$ point-like particles
connected by $N_{b}$ central-force springs in $d$ dimensions can be written to
quadratic order in $u$ using the dynamical matrix $\mathbb{D}^{\parallel}$ as

$E=\sum_{b=1}^{N_{b}}\delta
V_{b}=\sum_{b=1}^{N_{b}}\frac{k_{b}}{2}(e_{b}^{\parallel})^{2}=\frac{1}{2}\langle
u|\mathbb{D}^{\parallel}|u\rangle$ (16)

where the inner product is taken in the $Nd$ dimensional space for particle
displacements $u$, and the superscript “$\parallel$” on the dynamical matrix
signifies that this dynamical matrix describes a stress-free system where only
$e^{\parallel}$ enters the elastic energy. The more general form for
prestressed systems will be discussed in Sec. III.2. This quadratic form can
be decomposed into two steps using the equilibrium and compatibility matrices,
defined as

$\displaystyle\mathbb{C}^{\parallel}|u\rangle$
$\displaystyle=|e^{\parallel}\rangle$ (17)
$\displaystyle\mathbb{Q}^{\parallel}|t^{\parallel}\rangle$
$\displaystyle=-|f\rangle.$ (18)

Here $|e^{\parallel}\rangle$ and $|t^{\parallel}\rangle$ represent the
extension and tension of every bond, which are both $N_{b}$ dimensional
vectors, and the superscript “$\parallel$” indicated that they are along the
bond direction. $|u\rangle$ and $|f\rangle$ represent the displacement and
total force on every site, which are both $Nd$ dimensional vectors. As a
result, the equilibrium matrix $\mathbb{Q}$ has the dimension $Nd\times
N_{b}$, and the compatibility matrix $\mathbb{C}$ has the dimension
$N_{b}\times Nd$. The two matrices are in fact transpose of one another,
$\mathbb{Q}=\mathbb{C}^{T}$. Note we take the convention that
$t^{\parallel}>0$ denotes tension and $t^{\parallel}<0$ denotes compression.
This is consistent with the minus sign of Eq. (18), as $f$ is the force
_exerted by the bonds in the network_ on the sites.

Using these relations in Eq. (16), it is straightforward to see that

$\mathbb{D}^{\parallel}=\mathbb{Q}^{\parallel}\cdot\mathbb{K}^{\parallel}\cdot\mathbb{C}^{\parallel}$
(19)

where $\mathbb{K}^{\parallel}$ is a diagonal matrix that contains all the
spring constants $k_{b}$. In the sign convention we use here, Hookian’s law on
the springs takes the form
$t^{\parallel}=\mathbb{K}^{\parallel}e^{\parallel}$, as extension causes
tension on the springs.

The null space of $\mathbb{Q}^{\parallel}$ is the set of tensions, called
SSSs, that produce no forces at any site,

$\displaystyle 0$
$\displaystyle=\mathbb{Q}^{\parallel}|t^{\parallel}_{\text{SSS}}\rangle.$ (20)

The null space of $\mathbb{C}^{\parallel}$ is the set of site displacements,
called ZMs, that produce no changes in bond lengths,

$\displaystyle 0$ $\displaystyle=\mathbb{C}^{\parallel}|u_{\text{ZM}}\rangle.$
(21)

Floppy modes are a subset of ZMs excluding trivial rigid body motions of the
whole system, and have been extensively studied in soft matter systems, due to
their obvious significance as capturing deformations with no cost of elastic
energy. SSSs have only recently been explored in soft matter, but also show
great potential in characterizing stress-bearing structures in mechanical
networks Lubensky  _et al._ (2015); Sussman  _et al._ (2016).

Applying rank-nullity theorem on $\mathbb{Q},\mathbb{C}$ matrices leads to the
Maxwell-Calladine index theorem Calladine (1978); Sun  _et al._ (2012)

$\displaystyle N_{0}-N_{S}=Nd-N_{b},$ (22)

where $N_{0},N_{S}$ are the numbers of ZMs and SSSs.

These concepts have found wide applications recently in the new field of
topological mechanics, as $t^{\parallel}_{\text{SSS}}$ and $u_{\text{ZM}}$ can
become topologically protected modes in Maxwell lattices and networks Kane and
Lubensky (2014); Lubensky  _et al._ (2015); Rocklin  _et al._ (2017); Zhang
and Mao (2018); Zhou  _et al._ (2018); Ma  _et al._ (2018); Mao and Lubensky
(2018); Zhou  _et al._ (2019, 2020); Sun and Mao (2020).

###  III.2 States of self-stress in prestressed systems

Prestress on a force-balanced mechanical network can always be viewed as
“exciting” an existing SSS in the “stress-free version” of the same network
(which is generated by turning off stress in the prestressed network but
keeping exactly the same geometry—requiring to redefine the rest-length of the
springs). Thus, the stress-free version of a prestressed network must exhibit
a SSS in the first place.

With prestress, the elastic energy includes both $e_{b}^{\parallel}$ and
$e_{b}^{\perp}$ terms as analyzed in Eq. (3),

$\displaystyle E=$ $\displaystyle\sum_{b=1}^{N_{b}}\delta
V_{b}=\sum_{b=1}^{N_{b}}\left[\frac{k_{b}}{2}(e_{b}^{\parallel})^{2}+\frac{t_{b,p}}{2|\vec{R}_{b,0}|}|e_{b}^{\perp}|^{2}\right]$
$\displaystyle=$ $\displaystyle\frac{1}{2}\langle
u|(\mathbb{D}^{\parallel}+\mathbb{D}^{\perp})|u\rangle.$ (23)

Similar to the stress-free case, using the fact that this elastic energy
consists only of complete square terms, this dynamical matrix can also be
decomposed into $\mathbb{Q},\mathbb{C}$ matrices,

$\mathbb{D}=\mathbb{Q}\cdot\mathbb{K}\cdot\mathbb{C},$ (24)

where

$\displaystyle\mathbb{C}|u\rangle$ $\displaystyle=\begin{pmatrix}e^{||}\\\
e^{\perp}\end{pmatrix}\equiv|e\rangle,$ (25)

defines the new $\mathbb{C}$ matrix and $\mathbb{Q}=\mathbb{C}^{T}$. Here we
have defined the new $N_{b}d$ dimensional $e$ vector which contains one
component from $\parallel$ and $d-1$ from $\perp$ for each of the $N_{b}$
bonds. It is worth noting here that the compatibility matrix $\mathbb{C}$ is
now $N_{b}d\times Nd$ dimensional instead of $N_{b}\times Nd$ dimensional,
because this $\mathbb{C}$ matrix maps the $Nd$ dimensional displacement vector
of the network into $e_{b}^{\parallel}$ and $e_{b}^{\perp}$ for each bond $b$.
At the same time, the equilibrium matrix $\mathbb{Q}$ is $Nd\times N_{b}d$
dimensional. The spring constant matrix is a $N_{b}d\times N_{b}d$ diagonal
matrix with spring constant $k_{b}$ for the $\parallel$ terms and pre-stress
$t_{b,p}$ for the $\perp$ terms,

$\displaystyle\mathbb{K}=\begin{pmatrix}k_{b}&0\\\
0&\frac{t_{b,p}}{|\vec{R}_{b,0}|}\end{pmatrix},$ (26)

where there are $d$ blocks of $N_{b}$ dimensional diagonal matrices (one
longitudinal and $d-1$ transverse).

These new $\mathbb{Q},\mathbb{C}$ matrices describe the mapping between the
($Nd$ dimensional) degrees of freedom space and the (now $N_{b}d$ dimensional)
constraint space. In parallel with Eq. (25) we have

$\displaystyle\mathbb{Q}\begin{pmatrix}t^{||}\\\ t^{\perp}\end{pmatrix}$
$\displaystyle\equiv\mathbb{Q}|t\rangle=-|f\rangle$ (27)

Correspondingly, $t^{\parallel}$ and $t^{\perp}$ are the parallel and
perpendicular components of $t$.

Figure 3:  Examples of SSSs in prestressed systems. (a) A stress-free
mechanical network with one SSS. When the network is prestressed with this
SSS, 5 new SSSs arise (b-f). Bonds with $t^{\parallel},t^{\perp}$ in these
SSSs are denoted by the color scheme shown on the top, where $t^{\parallel}$
is denoted by red-white-blue as it goes from negative (compression) to
positive (tension), and the relative strength of $t^{\perp}$ over the total
$t\equiv\sqrt{(t^{\parallel})^{2}+(t^{\perp})^{2}}$ is denoted by green. The
thickness of the bond is proportional to $t$, and the black (gray) arrows
denote the component of $t^{\parallel},(t^{\perp})$ on the nodes. (g) The
tensegrity T3-prism has one SSS when it is stress-free. When it is prestressed
with this SSS, 20 new SSSs arise, and one example of these SSSs is shown in
(h). This example is chosen such that it carries the load of torsion between
the top and bottom triangles (see Sec. III.4 for the projection of load onto
SSSs).

It might appear confusing how tension $t$ on a central-force spring can have a
component perpendicular to the spring. This can be understood by realizing
that the $\parallel$ and $\perp$ components are defined with respect to the
reference configuration, and $t$ denotes an increment of stress in addition to
the prestress $t_{p}$. In the prestressed reference state, $t_{b,p}$ is along
bond $b$ in the reference configuration $\vec{R}_{b,0}$. In the state after
the (infinitesimal) deformation, the total tension is along bond $b$ in the
deformed configuration $\vec{R}_{b}$. These two bond directions,
$\vec{R}_{b,0}$ and $\vec{R}_{b}$, are not parallel in general. The $t$ field
in Eq. (27) represents the increment from the pretension to the tension after
the deformation, and thus has components parallel and perpendicular to the
original bond direction. In other words, under an (infinitesimal) deformation,
both bond (parallel) extension $e^{\parallel}$ and bond rotation $e^{\perp}$
cause increment of stress [Eq. (III.2)], and the two components are
$t^{\parallel}$ and $t^{\perp}$ respectively.

SSSs in a prestressed system are thus defined as any vectors that satisfy

$\displaystyle\mathbb{Q}|t_{\textrm{SSS}}\rangle$
$\displaystyle=\mathbb{Q}\begin{pmatrix}t^{||}_{\textrm{SSS}}\\\
t^{\perp}_{\textrm{SSS}}\end{pmatrix}=-|f\rangle=0,$ (28)

which represent respectively parallel and perpendicular change of bond
tensions that leave all particles in force balance, for a prestressed network.
In Fig. 3 we show two examples of mechanical networks where prestress
introduces extra SSSs that involve $t^{\perp}$ components. This formulation
characterizes the tensegrity phenomena Calladine (1978); Connelly and Whiteley
(1992) where prestress generates new SSSs for the system to carry other types
of load.

###  III.3 Zero modes in prestressed systems

Prestressed systems have two types of ZMs, type $A$, which leave
$e^{||}=e^{\perp}=0$ on all bonds, and type $B$, which violate this relation
but are still zero energy modes. In this section we discuss both types of ZMs.

Type $A$ ZMs in prestressed networks can be defined as modes that live in the
null space of $\mathbb{C}$,

$\displaystyle\mathbb{C}|u_{\textrm{ZM}}^{(A)}\rangle$
$\displaystyle=\begin{pmatrix}e^{||}\\\ e^{\perp}\end{pmatrix}=|e\rangle=0.$
(29)

The number of this type of ZMs satisfy the Maxwell-Calladine index theorem for
a stressed network, which can be proven by applying the rank-nullity theorem
on the new $\mathbb{Q},\mathbb{C}$ matrices,

$\displaystyle N_{0}^{(A)}-N_{S}=Nd-N_{b}d.$ (30)

Note that the last term is now $N_{b}d$ for the prestressed network instead of
$N_{b}$ in the stress-free case. More rigorously, if we also consider the
possibility that not all bonds are stressed, this index theorem should be
written as

$\displaystyle N_{0}^{(A)}-N_{S}=Nd-
N_{b}^{\textrm{unstressed}}-N_{b}^{\textrm{stressed}}d,$ (31)

where each unstressed bonds only provide one constraint.

Type $B$ ZMs in prestressed systems are not included in Eq. (29). Instead,
they are ZMs due to the fact that some spring constants in $\mathbb{K}$ are
negative (from compressively prestressed bonds), and their contribution can
cancel out the positive terms. They are “fine-tuning” ZMs which are only zero
energy when the prestress satisfies certain conditions so that the positive
and negative terms exactly cancel. This is also the point when varying the
magnitude of the prestress starts to cause instability, as we discuss at the
end of Sec. II.3. One example of such mode is shown in Fig. 2. Interestingly,
because $\mathbb{C}|u_{\textrm{ZM}}^{(B)}\rangle\neq 0$ and
$\mathbb{Q}\mathbb{K}\mathbb{C}|u_{\textrm{ZM}}^{(B)}\rangle=0$, type $B$ ZMs
lead to a set of SSSs
$|t_{\textrm{SSS}}\rangle=\mathbb{K}\mathbb{C}|u_{\textrm{ZM}}^{(B)}\rangle$.
One trivial example of this type of ZMs is the rigid rotation of the network
if the system is under open boundary conditions: the rotation causes
$e^{\perp}$ and thus does not satisfy Eq. (29), but it is indeed a ZM of the
dynamical matrix (this ZM does not require fine tuning of stress but all other
type B ZMs do). A nontrivial example of such fine-tuning ZMs is shown in Fig.
2e where the sum of $e^{\parallel}$ and $e^{\perp}$ terms in the elastic
energy vanish. Type $B$ ZMs are not included in the Maxwell-Calladine index
theorem for prestressed systems [Eq. (30)].

Remarkably, although ZMs in prestressed systems are more complicated as we
discussed above, the definition of SSSs in Eq. (28) is robust as it relies
only on force balance, and is not affected by the positive definiteness of
$\mathbb{D}$. In our discussions below, we mainly focus on these SSSs and show
how they form a linear space that efficiently characterizes how load is
carried by a prestressed system.

###  III.4 Stress response to homogeneous load

One important property of SSSs is that they form a linear space containing
_all_ possible ways a network can carry stress while keeping all particles in
force balance. Thus, when the network is under load, actual stress
distributions must come from linear combinations of SSSs, and the linear space
of SSSs characterizes the capability of a system to carry any external load.
It is worth pointing out here that we refer to external loads that do not
“introduce new constraints” to the system, and the meaning of this condition
will become clear in the end of Sec. III.5.

Here we first illustrate this formulation using a simple shear on a stress-
free network with all spring constants being 1. The bond tension in response
to this shear can be decomposed in the SSSs linear space because this space
form a complete basis for all force-balanced stress distributions, as we
discussed above. The decomposition can be written as

$|t^{\parallel}\rangle=\sum_{i=1}^{N_{\textrm{SSS}}}|t_{\textrm{SSS},i}^{\parallel}\rangle\langle
t_{\textrm{SSS},i}^{\parallel}|e_{\textrm{affine}}^{\parallel}\rangle,$ (32)

where $e_{\textrm{affine}}^{\parallel}$ is the bond extension if the strain
were affine (i.e., homogeneous shear strain). The sum runs in the
$N_{\textrm{SSS}}$ dimensional linear space of all SSSs,
$|t_{\textrm{SSS},i}^{\parallel}\rangle$. This relation is straightforward to
prove as follows: when a strain is imposed on the system (e.g., via Lees-
Edwards boundary conditions in a computational model, or in the bulk of a
mechanical network strained from the boundary), in addition to affine bond
extensions $e_{\textrm{affine}}^{\parallel}$, the system responds by particle
displacements $u$ (e.g., nonaffine deformations) to minimize the elastic
energy while maintaining this macroscopic strain, so the resulting bond
extensions are

$|e^{\parallel}\rangle=|e_{\textrm{affine}}^{\parallel}\rangle+\mathbb{C}^{\parallel}|u\rangle.$
(33)

If all spring constants are 1, we have
$|t^{\parallel}\rangle=|e^{\parallel}\rangle$. As we mentioned above,
$|t^{\parallel}\rangle$ must belong to the SSSs linear space, and the linear
combination coefficients are

$c_{i}=\langle t_{\textrm{SSS},i}^{\parallel}|t^{\parallel}\rangle=\langle
t_{\textrm{SSS},i}^{\parallel}|e_{\textrm{affine}}^{\parallel}\rangle+\langle
t_{\textrm{SSS},i}^{\parallel}|\mathbb{C}^{\parallel}|u\rangle.$ (34)

The second term has to vanish because $\langle
t_{\textrm{SSS},i}^{\parallel}|\mathbb{C}^{\parallel}=(\mathbb{Q}^{\parallel}|t_{\textrm{SSS},i}^{\parallel}\rangle)^{T}=0$.
This proves the decomposition in Eq. (32).

Three comments can be made from this result. First, if the shear strain has no
overlap with any SSSs in the system, $\langle
t_{\textrm{SSS},i}^{\parallel}|e_{\textrm{affine}}^{\parallel}\rangle=0$ for
all $i$, the system can not carry the load. What would happen physically is
that the system yields until a new SSS, able to carry that load, emerges.
Second, the same approach can also be applied to other types of load, such as
a hydrostatic pressure. It is important to note here that any component in the
imposed strain that can be written in terms of
$\mathbb{C}^{\parallel}|u\rangle$ will not cause stress—it corresponds to
strain that will be relaxed by degrees of freedom available to the system.
Third, a similar formulation can also be developed with loads applied via
specifically controlled boundary displacements, where SSSs are defined as bond
tensions leaving internal sites in force balance, and the strain $e$ can be
applied from bonds connected to boundaries Zhang and Mao (2018).

The relation is slightly more complicated when the spring constants are
different for each spring,

$|t^{\parallel}\rangle=\sum_{i,j}^{N_{\textrm{SSS}}}|t_{\textrm{SSS},i}^{\parallel}\rangle[((\mathbb{K^{\parallel}})^{-1})_{ss}]^{-1}_{ij}\langle
t_{\textrm{SSS},j}^{\parallel}|e_{\textrm{affine}}^{\parallel}\rangle,$ (35)

where $((\mathbb{K^{\parallel}})^{-1})_{ss}$ is the inverse of the spring
constant matrix projected to the SSS linear space (see detailed definition in
App. A).

This relation can readily be generalized to the pre-stressed case, where

$|t\rangle=\sum_{i,j}^{N_{\textrm{SSS}}}|t_{\textrm{SSS},i}\rangle[(\mathbb{K}^{-1})_{ss}]^{-1}\langle
t_{\textrm{SSS},j}|e_{\textrm{affine}}\rangle.$ (36)

A detailed proof of this relation is included in App. A. It is worth noting
that the response $|t\rangle$ calculated here is the addition to the prestress
in response to the load, so the total stress in the system is
$|t_{p}\rangle+|t\rangle$ where $t_{p}$ only has longitudinal components and
$t$ has both longitudinal and transverse components (with respect to the bond
direction in the prestressed reference state).

###  III.5 Stress response to force dipoles

Similar projections can also be applied to compute the stress response of a
prestressed system to local forces, such as a force dipole acting on a bond.

Figure 4: Stress response to dipole forces. (a) Bond stress response of a
mechanical network to a force dipole $f_{\textrm{dipole}}$ (black arrows)
acting on a bond in the network. Blue and red marks the tension
$t_{\textrm{rsp}}>0$ and compression $t_{\textrm{rsp}}<0$ on the bonds in
response to the force dipole, with the thickness of the bond proportional to
the magnitude of the tension/compression. (b) The sum of the external dipole
and the network response, $t_{\textrm{dipole}}+t_{\textrm{rsp}}$, is a SSS of
the network. (c) Bond stress response of a mechanical network to a force
dipole acting on a pair of particles not connected by a bond. This force
dipole can be viewed as an additional bond in the network, such that
$t_{\textrm{dipole}}+t_{\textrm{rsp}}$ is a SSS of the network with the added
bond (d). The original structure in (c) does not have any SSS.

The main difference between this case of a force dipole and the external
strain case discussed in Sec. III.4 is that the force-balance condition here
needs to include the external forces,

$\displaystyle\mathbb{Q}|t_{\textrm{rsp}}\rangle=-|f_{\textrm{rsp}}\rangle=|f_{\textrm{external}}\rangle$
(37)

where $f_{\textrm{external}}$ is the externally imposed force and
$t_{\textrm{rsp}}$ is the tension response of the system, which causes
$f_{\textrm{rsp}}$. $f_{\textrm{rsp}}$ cancels with $f_{\textrm{external}}$,
leaving the system in force balance.

When the external forces take the form of a force dipole on a bond $b$, they
can be written as

$\displaystyle|f_{\textrm{external}}\rangle=|f_{\textrm{dipole}}\rangle=-\mathbb{Q}|t_{\textrm{dipole}}\rangle=-\mathbb{Q}\mathbb{K}|b\rangle$
(38)

where $|b\rangle$ is a vector in the bond space with 1 on the bond where the
dipole is imposed, and 0 on all other bonds. This sets the magnitude of the
dipole force to be $k_{b}\cdot 1$, the value of which is irrelevant for the
linear theory. The sign convention of this equation follows from Eq. (27)
where $t_{\textrm{dipole}}$ causes $f_{\textrm{dipole}}$. Combining Eqs. (37)
and (38) we have

$\displaystyle\mathbb{Q}(|t_{\textrm{dipole}}\rangle+|t_{\textrm{rsp}}\rangle)=0,$
(39)

indicating that the total tension field $t_{\textrm{dipole}}+t_{\textrm{rsp}}$
(sum of the external dipole and the network response) is a SSSs of the
network. The response can then be derived using a similar method as discussed
in Sec. III.4, and we have (as detailed in App. B)

$|t_{\textrm{rsp}}\rangle=-|t_{\textrm{dipole}}\rangle+\sum_{i,j}^{N_{\textrm{SSS}}}|t_{\textrm{SSS},i}\rangle[(\mathbb{K}^{-1})_{ss}]^{-1}_{ij}\langle
t_{\textrm{SSS},j}|b\rangle.$ (40)

In the simple case of $\mathbb{K}=k\mathbb{I}$ (all springs having the same
spring constant which only applies to stress-free networks), the second term
on the right hand side of the equation reduces to the “quasi-localized SSSs”
defined in the stress-free case, characterizing dipole responses Sussman  _et
al._ (2016); Lerner (2018). Here we extend this concept to prestressed
networks, and our formulation can be used to not only compute dipoles with
force along the bond, but also “transverse” or “tilted” dipoles where the
force pair is perpendicular or in an arbitrary angle to the bond (see example
in Fig. 4).

A “dipole stiffness” $\kappa$ can then be defined to characterize this
response

$\kappa\equiv\frac{\langle{f_{\textrm{dipole}}|f_{\textrm{dipole}}}\rangle}{\langle{f_{\textrm{dipole}}|u_{\textrm{rsp}}}\rangle},$
(41)

where $\langle{f_{\textrm{dipole}}}|$ serves the purpose of projecting the
force and displacements to the direction of the external force dipole. A
different version of (longitudinal only) dipole stiffness has been introduced
in Ref. Rainone  _et al._ (2020a) where $u_{\textrm{rsp}}$ of the whole system
is used instead of its projection to the dipole was used. In contrast, the
dipole stiffness defined here treats the whole system as a black box, and only
extracts the stiffness from the force-distance relation between the two
particles where the dipole acts on. The simplicity of this definition allows
it to be directly implemented in experiments via methods such as
microrheology.

Using the SSSs linear space, the dipole stiffness to a local force dipole on
bond $b$ is thus given by (details in App. B):

$\kappa_{b}=\frac{2k_{b}^{2}}{k_{b}-\displaystyle\sum_{i,j}^{N_{\textrm{SSS}}}\langle
b|t_{\textrm{SSS},i}\rangle[(\mathbb{K}^{-1})_{ss}]^{-1}\langle
t_{\textrm{SSS},j}|b\rangle},$ (42)

where $k_{b}$ is the spring stiffness for bond $b$.

Another interesting case of force dipole is when the dipole acts on a pair of
particles that are not connected in the network (which can be called “non-
local dipoles”). The analysis of force-balance in this case can be done by
considering the external force dipole as _introducing an additional constraint
into the network_. Thus, the stress distribution is a SSS of the network with
this additional “auxiliary bond”, instead of the original network. In App. B
we also derive the dipole stiffness of this case. Examples of both types of
dipoles are shown in Fig. 4.

##  IV Prestressed elasticity of triangular lattices

In this section we use prestressed triangular lattices to illustrate the
mathematical tools for prestressed elasticity we discussed in Secs. II and
III. Models of prestressed triangular lattices have been introduced in the
literature of force network ensembles, as an example system where multiple
internal stress distributions are compatible with given boundary stress
Snoeijer  _et al._ (2004); Tighe  _et al._ (2008, 2010). Here we discuss how
our SSSs formulation can be applied to these lattices to both systematically
generate all possible prestress distributions, and solve for vibrational modes
and stress response to external load, revealing interesting prestress effects.

In the stress-free case, triangular lattices have a coordination number
$z=6>2d$ and thus the number of SSSs per site is $\frac{z}{2}-2=1$. In fact, a
localized SSS $|t_{\ell}\rangle$ can be defined around each site $\ell$ as
shown in Fig. 5a. These site-localized SSSs, together with the homogeneously
state of self stress (all bonds carrying the same tension), form a complete
basis to decompose any SSSs on the triangular lattice.

Figure 5: Vibrational modes of prestressed triangular lattices. (a) One site-
localized SSS $|t_{\ell}\rangle$ on a triangular lattice, where bonds
surrounding the hexagon carry compression, and bonds in the hexagon carry
tension. The size of the blue disk on the center site represents the strength
of this SSS. (b) A prestressed triangular lattice with the prestress $t_{p}$
generated using a linear combination of site-localized SSSs with independent
coefficients on each site, drawing from a Gaussian distribution with mean 0
and standard deviation $\tilde{c}=0.12$. (c) Lowest 300 eigenvalues $\lambda$
of the dynamical matrix of prestressed triangular lattices of size $32\times
32$, at $\tilde{c}=0$ (blue), $\tilde{c}=0.1$ (orange), $\tilde{c}=0.2$
(green), from top to bottom. (d) The lowest eigenvalue
$\lambda_{\textrm{min}}$ as a function of $\tilde{c}$. (e) A triangular
lattice with positional disorder where the random displacements of the nodes
have standard deviation $\tilde{u}=0.12$. (f) Lowest 300 eigenvalues $\lambda$
of the dynamical matrix of a $32\times 32$ triangular lattices with positional
disorder, at $\tilde{u}=0$ (blue), $\tilde{u}=0.1$ (orange), $\tilde{u}=0.2$
(green), from top to bottom. (g) Phonon disperson relation of the lowest band
of triangular lattices with no prestress (yellow, upper), and critical
compressive prestress (blue, lower) where modes along $\Gamma M$ approach
instability.  Figure 6: Prestressed triangular lattices under load. (a-c) An
example of a prestressed triangular lattice (the one shown in Fig. 5b) under
shear load $\epsilon_{xy}$ (a), a longitudinal dipole (b), and a transverse
dipole (c), where the force dipoles are shown as black arrows. The color
scheme of bond tensions ($t^{\parallel}$ and $t^{\perp}$) is the same as in
Fig. 3. (d) Probability distributions of pretensions on bonds of the lattices
we tested for load bearing, with different mean $t_{0}$ and standard deviation
$\tilde{t}$ of pretension on bonds [legend shown as inset of (e)]. (e) Shear
modulus $G$ of these lattices as a function of mean pretension $t_{0}$. (f-g)
Probability distributions of the parallel (f) and transverse (g) local dipole
stiffness of these lattices. $20\times 20$ lattices are used in the
computations in (d-g), and we average over 100 realizations of disorder for
each set of parameters.

In particular, we can generate an ensemble of prestress $|t_{p}\rangle$ on
triangular lattices by taking an arbitrary coefficient $c_{\ell}$ for each
site-localized SSS $|t_{\ell}\rangle$ and sum them up. In Fig 5b we show an
example of such a prestressed state.

These prestressed lattices can be physically constructed by choosing the rest
length of each bond $b$ such that when they are at the length in the regular
triangular lattice, the tension/compression they carry is exactly the value
$t_{b,p}$ in the prescribed SSSs. Because the total force is balanced at each
site for any SSS, all bonds will stay at the length in the regular lattice,
and thus we obtain a triangular lattice with _regular geometry and a
prescribed SSS_.

What are the mechanical properties of such prestressed triangular lattices?
Here we study them from two perspectives: vibrational modes and load-bearing
capabilities. The vibrational modes can be calculated based on the quadratic
expansion in Eq. (3) for each bond, which leads to the dynamical matrix
defined in Eq. (24).

In Fig. 5c we show eigenvalues of the dynamical matrix of the triangular
lattice at three levels of prestress (by assigning $c_{\ell}$ from a Gaussian
distribution with mean at $0$ and standard deviation at $\tilde{c}$).
Prestress significantly affects the vibrational modes especially at low
frequencies. In particular, negative eigenvalues start to appear when the
fluctuations of the SSSs goes beyond a critical level,
$\tilde{c}>\tilde{c}^{*}$ (Fig. 5c,d). This indicates that by increasing the
fluctuation of prestress (where the mean remains 0), unstable modes appear,
although the system remains in force balance. This is in alignment with our
discussion on prestressed rigidity and type $B$ ZMs of prestressed systems in
Sec. II.3.

To contrast this result, we also consider triangular lattices with positional
disorder but no prestress (Fig. 5e). In this case, we move site $\ell$ by a
random displacement $\vec{u}_{\ell}$, the $x,y$ components of which being
generated from a Gaussian distribution with mean at $0$ and standard deviation
at $\tilde{u}$. We also plot the eigenvalues of the dynamical matrix of this
case in Fig. 5f. Although positional disorder also affects the eigenvalues, it
mostly smooths out the regular lattice eigenmodes, and does not lead to
qualitative change at low frequencies and does not generate unstable modes.
Note that to make a fair comparison we chose the same values for $\tilde{c}$
and $\tilde{u}$, which represents the same level of disorder, under the simple
convention we took where the bond length and spring constant on the triangular
lattice both being unity.

The load-bearing capabilities of prestressed triangular lattices can be
analyzed by the SSSs formalism discussed in Sec. III. At any given realization
of disordered prestress, we can compute the linear space of SSSs, and use Eq.
(36) and (40) to find the stress response of the system to external load (Fig.
6a-c). In particular, to characterize the homogeneous component of prestress
(hydrostatic pressure), we also introduce a constant pretension $t_{0}$ on all
bonds, in addition to the fluctuations from $|t_{\ell}\rangle$. As a result,
the tension on each bond obeys a normal distribution with mean at $t_{0}$ and
standard deviation at $\tilde{t}=2\tilde{c}$ (as each bond appears in 4 site-
localized SSSs around it), as shown in Fig. 6d. We apply two types of load,
simple shear and dipole forces, on these prestressed triangular lattices, and
the results are shown in Fig. 6. This analysis shows that for both shear
modulus and stiffness against dipole forces, their mean is mainly controlled
by $t_{0}$ (where compression decreases stiffness and tension increases
stiffness) whereas $\tilde{t}$ slightly decreases these stiffness, as well as
broadens them. In addition, the stress response demonstrates interesting
heterogeneities due to the disordered prestress, an effect we will discuss
more in Sec. V.

Effects of prestress on mechanics can also be seen from the perspective of
phonon structures of regular lattices (Figs. 2c and 5g). It is straightforward
to see that negative prestress (compression) destabilizes the originally
stable triangular lattice (Fig. 5g), and positive prestress (tension)
stabilizes the honeycomb lattice (Fig. 2c), which was originally unstable
against shear. Interestingly, the modes that first become unstable in the
triangular lattice as negative prestress increase are the modes along the
$\Gamma M$ direction in the first Brillouin zone, indicating that the type of
modes that first become unstable are the modes that zig-zag between the
straight lines of bonds, agreeing with recent studies of strain-localization
in lattice models Bordiga  _et al._ (2021).

##  V Prestressed elasticity of amorphous solids of soft repulsive particles

The new $\mathbb{Q},\mathbb{C}$ matrix approach described in Sec. III provides
a set of tools to analyze stress-bearing capabilities and mode structures of
prestressed systems, well beyond the triangular lattices just described. Here
we apply this method to a computational model for amorphous solids composed of
soft repulsive particles, where, akin to real materials, prestresses are
naturally introduced by the preparation (solidification) protocol. Given the
preparation protocol chosen for this study, the prestress is the result of
both compression and frozen-in structural disorder.

###  V.1 3D numerical simulations of amorphous solids of soft repulsive
particles

The numerical model describes a suspension of particles with soft repulsive
interactions given by a truncated and shifted Lennard-Jones potential Weeks
_et al._ (1971). The potential energy, for a pair of two particles $\ell$ and
$\ell^{\prime}$ separated by a center-to-center distance
$R_{b}\equiv|\vec{R}_{\ell}-\vec{R}_{\ell}^{\prime}|$, is
$V(R_{b})=4\epsilon[(a_{b}/R_{b})^{12}-(a_{b}/R_{b})^{6}]+\epsilon$ for
particle pairs closer than the cutoff distance $R_{b}\leq r_{c}\equiv
2^{1/6}a_{b}$, otherwise $V(R_{b})=0$, following the notation we defined in
Sec. II. Here $a_{b}=(a_{\ell}+a_{\ell^{\prime}})/2$ with $a_{\ell}$ and
$a_{\ell}$ being the diameters. The size poly-dispersity is introduced by
drawing the diameter of each particle $a_{i}$ from a Gaussian distribution
with mean $a$ and variance of $10\%$. The parameter $\epsilon$ is the unit
energy and $a$ is the unit length in the simulations. Particles closer than
the cutoff distance (beyond which the inter-particle force vanishes),
$R_{b}\leq r_{c}$, are considered in contact. All simulations used here have
volume fraction $\phi\approx 70\%$ and consist of $10^{4}(10976)$ particles in
a cubic box of linear size $L=20.36$, unless otherwise specified. The initial
samples are prepared by first melting face-centered-cubic (FCC) crystals of
particle at volume fraction $0.7$ at $T=5.0\epsilon/k_{B}$. We then cool the
particle assemblies down to a temperature $0.001\epsilon/k_{B}$ through an NVT
molecular dynamics (MD) protocol as described in Ref. Vasisht and Del Gado
(2020). The cooling rate $\Gamma$ varies from $5\times 10^{-2}$ to $5\times
10^{-6}~{}\epsilon/(k_{B}\tau_{0})$, where $\tau_{0}=a\sqrt{m/\epsilon}$ is
the MD time unit with $m$ the particle mass. Subsequently, each sample is
brought to the closest local energy minimum using a conjugate gradient (CG)
algorithm. The configurations prepared in this way are amorphous solids, as we
verify by measuring their viscoelastic response Vasisht and Del Gado (2020),
whose properties depend on the cooling rate utilized. These solids feature
prestress that comes both from the homogeneous compression and from the
structural disorder due to the size polydispersity.

From all contacting particle pairs, we compute the virial stress tensor Irving
and Kirkwood (1950) (neglecting the kinetic term):

$\displaystyle\sigma_{ij}=\frac{1}{V}\sum_{\ell}\sum_{\ell^{\prime}>\ell}R_{b,i}t_{b,j},$
(43)

where $V$ is the volume of the system, $R_{b}$ is the distance between
particle $\ell$ and $\ell^{\prime}$ (length of bond $b$), and $t_{b}$ is the
force on bond $b$.

In Fig. 7 we show a typical configuration where prestress on the bonds are
visualized. These $t_{b,p}$ terms are negative as the system is compressed.
For each cooling rate, we prepare 5 statistically independent samples. All
quantities investigated here have been averaged over this set of samples and
we obtain error bars from sample to sample fluctuations.

We also examine these samples as they are sheared using Lees-Edwards boundary
conditions and a shear rate $\dot{\gamma}$, by solving Newton equations of
motions with a drag force that guarantees minimal inertia effects as discussed
in Vasisht and Del Gado (2020). All simulations are performed with LAMMPS
Plimpton (1995), suitably modified to include the particle size polydispersity
and the interactions discussed above, and using a shear rate
$\dot{\gamma}=10^{-4}\tau_{0}^{-1}$.

With the bond network and particle coordinates, we build the compatibility and
equilibrium matrices $\mathbb{Q}$ and $\mathbb{C}$, to study the prestressed
elasticity using the general formalism described in the previous sections. In
this 3D network, each bond has one longitudinal $t^{||}$ and two transverse
directions $t^{\perp}$ for tension increments. We take the convention that the
first transverse direction
$\hat{R}_{b,0}^{\perp,1}\equiv\hat{R}_{b,0}\times\hat{z}/|\hat{R}_{b,0}\times\hat{z}|$
is the unit vector of the cross product between the bond unit vector
$\hat{R}_{b,0}$ and the $z$-axis $(0,0,1)$, and the second transverse
direction
$\hat{R}_{b,0}^{\perp,2}\equiv\hat{R}_{b,0}\times\hat{R}_{b,0}^{\perp,1}$ is
the unit vector of the cross product between the bond vector and the first
transverse direction [the subscript $0$ signifies the reference (undeformed
but prestressed) state].

Figure 7: Prestress in amorphous solids of repulsive particles. (a) Pretension
on bonds $t_{b,p}\equiv\frac{dU(r_{ij})}{dr_{ij}}$ where $b$ is the bond
connecting particles $i,j$. The system size is $N=10976$ with preparation
cooling rate $\Gamma=5\times 10^{-6}~{}\epsilon/(k_{B}\tau_{0})$. (b-c)
Distributions of $t_{b,p}$ (b), and spring constants $k_{\parallel}\equiv
k_{b}$ and $k_{\perp}\equiv\displaystyle\frac{t_{b,p}}{|\vec{R}_{b,0}|}$ (c),
for the configuration shown in (a). (d-g) $U$, $t$, $k_{\parallel}$, and
$k_{\perp}$ for the pairwise potential used in the model as functions of
distance $R_{b}$.

Once the $\mathbb{Q}$ and $\mathbb{C}$ matrices are constructed, we solve
their null space and obtain ZMs and SSSs of the computational model. The
samples we studied are in general deep in the solid phase and exhibit no ZMs
besides the 3 trivial translations. When solving for SSSs in large and dense
systems as the samples we generated, the null space for the equilibrium matrix
$\mathbb{Q}$ may have very high dimensions. We used the SPQR_RANK package
Foster and Davis (2013) from SuiteSparse sui , which is a high performance
sparse QR decomposition package and can provide a reliable determination of
null space basis vectors for large sparse matrices, to solve for SSSs
efficiently and reliably, especially when there are lots of degeneracies for
diagonalizing $\mathbb{Q}$. Similarly, SPQR_RANK is also suitable when solving
for the null space of $\mathbb{Q}^{T}$ to get ZMs in large systems.

###  V.2 Calculating stress response to shear using states of self-stress

We begin by demonstrating the power of the SSSs as “stress-eigenstates” that
characterize the stress-bearing capabilities of a system. To this end, we
compute all SSSs as described above, and use Eq. (36) to calculate the change
of tension, $t^{\parallel}$ and $t^{\perp}$, as the system is under
macroscopic shear, which we call “prestressed response” here, as this is the
difference between the stress under load and the prestress.

To contrast this result, we also calculated the tension increment treating the
system as a stress-free network, i.e., as characterized by equilibrium matrix
$\mathbb{Q}^{\parallel}$ instead of $\mathbb{Q}$ (where the geometry of the
contact network is kept the same). The resulting stress increment (which only
has the $t^{\parallel}$ component) is then calculated using Eq. (35), and we
will call this the “stress-free response”. At the same time, we also measure
the change of the stress in the system in a numerical experiment where the
system is under a small shear deformation imposed quasistatic or at finite
rate (as specified in captions of Figs. 8 and 9), and we call this the “actual
response”.

In the next two subsections we compare these different results on tension
increment in terms of both spatial heterogeneity and shear modulus.

###  V.3 Spatial heterogeneity of stress response

Figure 8: Spatial heterogeneity of tension increment of prestressed networks
under shear strain. We compare the actual tension response $t$ measured in
simulation to the tension increment calculated with and without prestress, for
a quasistatic shear deformation at strain 0.01% in (a) and a finite rate shear
deformation at strain 1% in (b) (system size $N=10976$ and cooling rate
$\Gamma=5\times 10^{-2}~{}\epsilon/(k_{B}\tau_{0})$). (c) Clustering tendency
index $H$ averaged over three different samples calculated using stress-free
and prestressed elasticity compared to the actual $H$ measured at quasistatic
shear strain $0.01\%$. The red dashed line indicates the threshold to
determine whether or not the tension change is clustered.

One significant effect of prestress on the mechanical response of amorphous
solids is the spatial heterogeneity of the stress increment, in contrast to
the stress-free response (where the geometry is kept the same). Fig. 8 shows a
visualization of the tension increment in the model amorphous solids we study
here.

Remarkably, this heterogeneity is accurately depicted from the SSSs
calculation when prestress is included. The tension increment is computed
following Eq. (36). The resulting $t$ field includes both $t^{||}$ and
$t^{\perp}$, although only $t^{||}$ is shown. To be precise, $t^{\perp}$
manifests as rotations of the bonds in the real response (which causes the
total tension to change direction). In contrast, the calculated stress
response of the system ignoring the prestress effect, i.e., using Eq. (35), is
almost completely homogeneous. Thus, in these dense systems, _the prestress
controls the spatial heterogeneity of the mechanical response of the system._
Treating the system as “stress-free” misses important signatures of
heterogeneities.

We characterize the spatial heterogeneity of the tension increment using a
clustering tendency index: the Hopkins statistic $H$ Hopkins and Skellam
(1954), with a value close to 1 indicating the data is highly clustered, and a
value around 0.5 from random data. Fig. 8(c) shows that the prestressed
elasticity captures the tension increment clustering tendency which exists in
real tension responses, while in stress-free elasticity this clustering is not
captured.

###  V.4 Shear modulus of prestressed amorphous solids

The shear modulus $G$ can be obtained from the shear component of the virial
stress. All the contributions to the virial stress in Eq. (43) can be obtained
from the SSSs projection calculation (which directly gives
$t^{\parallel}_{b},t^{\perp}_{b}$ of all bonds $b$) via

$\displaystyle\vec{t}_{b}$
$\displaystyle=(t_{b,p}+t^{\parallel}_{b})\hat{R}^{\parallel}_{b,0}+t^{\perp,1}_{b}\hat{R}^{\perp,1}_{b,0}+t^{\perp,2}_{b}\hat{R}^{\perp,2}_{b,0},$
(44) $\displaystyle\vec{R}_{b}$
$\displaystyle=\left(R_{b,R}+\frac{|\vec{t}_{b}|}{k_{b}}\right)\hat{t}_{b},$
(45)

where $\hat{t}_{b}\equiv\vec{t}_{b}/|\vec{t}_{b}|$ is the direction of the
bond after the deformation. Comparing with the continuum limit we discussed at
Eq. (II.2), this virial stress tensor is both a spatial average over the whole
system, and uses the bond directions after deformation $\vec{R}_{b}$ (which is
along the calculated $\vec{t}_{b}$) instead of original directions
$\vec{R}_{b,0}$. In the terminology of structural mechanics, this virial
stress [Eq. (43)] is the Cauchy stress (true stress, in the thermodynamic
limit), and Eq. (II.2) defines the Second Piola-Kirchhoff stress which uses
bond directions in the undeformed state. [In practice, using Eq. (II.2) would
just result in slightly higher errors.]

In Fig. 9 we show the agreement between the actual and the calculated tension
increment using the prestressed formulation on the bonds, as well as the
agreement of the shear modulus. In contrast, the tension increment calculated
from the stress-free formulation displays clear deviations.

This effect can also be viewed from the continuum elasticity formulation in
Sec. II.2: when the macroscopic shear deformation field $\partial_{y}u_{x}$ is
plugged into the continuum theory, the elastic energy density is

$E/V=\frac{1}{2}(K_{xyxy}+\sigma_{p,yy})(\partial_{y}u_{x})^{2},$ (46)

leading to a prestress-corrected shear modulus of $G=K_{xyxy}+\sigma_{p,yy}$.
It is straightforward from this relation that an isotropic compressional
prestress ($\sigma_{ij}=-p\delta_{ij}$) decreases the shear rigidity of a
network.

The results shown so far demonstrates that including prestress effects is
important to characterize the mechanical response of the model amorphous
solids obtained from the simulations. We next use the configurations of the
model amorphous solids to collect statistical information on the SSSs when
different preparation protocols, i.e. different cooling rates, are utilized.

Figure 9: Comparing tension increment calculated from SSSs projection and
numerical experiments (actual tension) of quasistatic $0.01\%$ shear strain
[system size N = 10976 and cooling rate $\Gamma=5\times
10^{-6}~{}\epsilon/(k_{B}\tau_{0})$]. (a-b) Comparison of calculated and
actual tension increment ($t^{\parallel},t^{\perp}$) of each bond, where good
agreement is found in the case of prestressed formulation. (c) Shear modulus
$G$ from stress-free response and prestressed response [Eq. (43)] compared to
the actual $G$. The shear modulus is averaged over three different
configurations at quasistatic $0.01\%$ shear strain, with the same system size
and preparation protocol to the system in (a-b).

###  V.5 General statistics of states of self-stress

By varying the cooling rate as described above, we prepare particle assemblies
that have different amount of disorder and compressive prestress, and have
shown elsewhere that the higher cooling rates lead to a higher degree of
disorder, characterized through a Voronoi analysis of the local particle
packing which shows a decrease in the local icosahedral order Vasisht  _et
al._ (2020); Vasisht and Del Gado (2020).

The total number of SSSs in a prestressed system $N_{S}$ is directly related
to the numbers of degrees of freedom and constraints via the Maxwell-Calladine
index theorem [Eq. (30)] as the system exhibits no ZMs except for the trivial
translations. In Fig. 10(a) we show the total number of SSSs in different
configurations of our model amorphous solid that were prepared by varying the
cooling rate. The number $N_{S}$ decreases with the increase of cooling rate.
Figs. 10(b)-(d) show that this effect is accompanied by a decrease in the
number of bonds and a decrease of the shear modulus, which is also associated
to a net increase of the average compressional prestress. These findings
support the idea that varying the preparation protocol by varying the cooling
rate, as we do, induces different amount of prestress into amorphous solid
configurations, which can be directly captured by the SSSs analysis. The data
also indicate that more aged samples, which tend to be stiffer, feature larger
amounts of SSSs. This is consistent with the finding in Ref. Vasisht  _et al._
(2020) that more aged samples (i.e., prepared with a lower cooling rate)
contained larger amounts of (overconstrained) icosahedrally packed domains,
identified through the Voronoi analysis of the particle packing. These
domains, which tend to be stiffer, were shown to favor the accumulation of
stress and promote dilation when the samples were driven towards yielding
under a shear deformation. The comparison of the results obtained here with
the analysis performed in Refs. Vasisht  _et al._ (2020); Vasisht and Del Gado
(2020) support the idea that those phenomena, i.e., the increase in stiffness
and in the tendency to accumulate stress and dilate under shear, should be due
to prestress, and that sizable changes in the number of SSSs, as a result of
sizable differences in prestress, eventually determine sizable changes in the
linear, and even nonlinear, response of amorphous solids.

Figure 10: Statistics of SSSs at different cooling rate. (a) Number of SSSs
($N_{\textrm{S}}$), (b) number of bonds ($N_{b}$), (c) measured shear modulus
($G$), (d) averaged normal stress ($P$) for different cooling rates $\Gamma$.
Those quantities are averaged over 5 different configurations for each cooling
rate, and the system size is $N=10976$.)

###  V.6 Dipole stiffness in prestressed glasses

SSSs have been also discussed in the context of the mechanical response of
amorphous solids to local perturbations (e.g., force dipoles), which could be
connected to the localized or quasi-localized plastic processes that are key
ingredients of the mechanical response of this class of materials Sussman  _et
al._ (2016); Lerner (2018). The analysis of nonaffine and plastic microscopic
displacements under shear and of their spatial correlations has led to the
notion of “soft spot”, i.e., the presence of localized regions that are
especially susceptible to plastic rearrangements under an external forceTanguy
_et al._ (2002b); Leonforte  _et al._ (2004); Maloney and Lemaître (2006);
Chen  _et al._ (2011); Schoenholz  _et al._ (2014).

Figure 11: Dipole responses in amorphous solids. (a-d) Tension responses to a
transverse (a) or longitudinal (c) force dipole applied on the a bond (black
arrows). The color scheme showing $t^{\parallel},t^{\perp}$ is the same as the
one used in Fig. 3. (b) and (d) are zoomed-in views near the bond in (a) and
(c) respectively. (e-f) Spatial decay of tension responses for the two force
dipoles in (a) and (c) respectively. (g) Distribution of minimum dipole
stiffness $\kappa_{b}^{\textrm{min}}$ for systems with different preparation
cooling rates. System size is $N=2916$. The distribution is averaged over 10
samples for each cooling rate. (h) Spatial distribution of
$\kappa_{b}^{\textrm{min}}$ on bonds for a configuration of system size
$N=2916$ and preparation cooling rate $\Gamma=5\times
10^{-6}~{}\epsilon/(k_{B}\tau_{0})$. Bonds with soft
$\kappa_{b}^{\textrm{min}}$ are thicker and more opaque. (i) Pretension
$t_{b,p}$ on bonds for the same configuration in (h).

The local response to force dipoles that we have discussed in Sec. III.5 has
been directly related to these soft spots in amorphous solids Sussman  _et
al._ (2016); Lerner (2018); Rainone  _et al._ (2020b) and such response
analysis has pointed to quasilocalized excitations (QLEs). In general, one
needs to follow the time evolution of a system under an external deformation
to identify the QLEs, however similar characteristics have been extracted, in
the case of model glasses and amorphous solids, even from the low-frequency
limit of the vibrational spectrum, i.e., from the linear response regime
Tanguy  _et al._ (2002b); Chen  _et al._ (2011); Schoenholz  _et al._ (2014);
Mosayebi  _et al._ (2014).

Here we use the newly developed methods to examine the response of the model
amorphous solids obtained in simulations to dipole forces, and collect the
statistics of the local stiffness these systems display. In particular, we
apply the approach for force dipoles discussed in Sec. III.5, and calculate
the stiffness as defined in Eq. (41) for each bond, to identify regions which
are soft. Measures of this kind lend themselves, in fact, to identify soft
regions with respect to various external perturbations.

Interestingly, the relative displacement that a pair of particles undergoes
during a plastic rearrangement is typically not along the direction of the
bond that initially connects them. Instead, this local deformation involves
several sliding motions. Therefore, a particularly interesting measure of
local softness, potentially complementing recently proposed ones Lerner
(2018); Cubuk  _et al._ (2015); Richard  _et al._ (2020), can be obtained by
quantifying the stiffness against a dipole, _whose forces are along a
different direction from the one of the bond on which they are applied._ Thus,
for each bond, in addition to computing the dipole stiffness to longitudinal
(along $t^{\parallel}$) and transverse (along $t^{\perp}$) forces, we compute
the softest dipole stiffness, $\kappa_{b}^{\textrm{min}}$, by minimizing it
with respect to the direction of the dipole forces (details in App. B).

Fig. 11(a) shows an example of applying longitudinal and transverse force
dipoles on a bond in our prestressed particle configurations, showing both
$\parallel$ and $\perp$ tension responses of all bonds. Remarkably, although
the near field response is mainly in the $\parallel$ direction, the far field
response contains significant $t^{\perp}$ components. This is a genuine effect
of prestress—these transverse stress response can only be captured when
prestress is included.

We collect the statistics of the minimum dipole stiffness
$\kappa^{\textrm{min}}_{b}$ on bonds in the networks obtained with two
different preparation cooling rates as shown in Fig. 11(g). These
distributions feature a power-law tail for small $\kappa^{\textrm{min}}_{b}$,
where $p(\kappa^{\textrm{min}}_{b})\sim(\kappa^{\textrm{min}}_{b})^{\beta}$
and $\beta$ appears to depend on the cooling rate. The results demonstrate the
impact of preparation history to the mechanical stiffness in prestressed
systems. The case with faster cooling rate leads to a smaller exponent and
thus the distribution extends more to soft $\kappa^{\textrm{min}}_{b}$,
agreeing with the general trend discussed in Sec. V.5, where faster cooling
rate leaves more heterogeneities regarding stress response in the sample. Our
results also provide an interesting comparison to (longitudinal) dipole
stiffness distributions studied in Ref. Rainone  _et al._ (2020a) for
different models of computational glasses.

Furthermore, we also characterize the spatial distribution of
$\kappa^{\textrm{min}}_{b}$ and compare it to the spatial distribution of
prestress $t_{b,p}$ (Fig. 11h,i). Remarkably, $\kappa^{\textrm{min}}_{b}$
shows interesting spatial patterns whose possible correlations with the amount
of prestress in these regions are however hard to tease out. Relating regions
with clustered soft $\kappa^{\textrm{min}}_{b}$ with soft spots will be an
interesting question for future studies.

##  VI Discussion

Through this paper, we have developed a set of mathematical tools to
investigate elasticity of prestressed discrete networks, which allows to
disentangle the prestress from the structural disorder, offering new insight
into the mechanical response of amorphous solids.

Our theory is based on the equilibrium and compatibility matrices, which
rigorously defines SSSs and ZMs in mechanical networks. SSSs provide a linear
space for the ensemble of force-balanced prestress distributions of the
network, therefore capturing the possibility that the prestress may vary
without altering the network configuration. Any mechanical network with SSSs
can become prestressed, e.g., by introducing mismatching bond lengths
(geometric frustration), which naturally arises in the solidification of
amorphous solids in general. These We have generalized the equilibrium and
compatibility matrices to these prestressed systems, which amounts to include
additional constraints on bonds from prestress. We have then used the new
mathematical tools to analyze the mechanical response of prestressed
triangular lattices and model amorphous solids obtained through 3D numerical
simulations, and found good agreement between the numerically measured
response and results obtained from our methods, whereas the time-complexity of
the computation is greatly reduced.

In particular, we not only characterize how prestress affects vibrational
modes, both stabilizing soft modes and introducing a new class of ZMs (type
$B$, which appear due to competing positive and negative terms in the energy),
but also reveal its profound role on how the network carries additional
stress. By analyzing prestressed SSSs of the network, and projecting external
load to this linear space, our methods conveniently map how external stress
transmits in prestressed lattices and bond networks obtained from model
amorphous solids. We show a number of intriguing effects which are unique to
prestressed networks, including the spatial heterogeneity in the stress-
response to homogeneous loads, the power-law distribution of the local
stiffness and the isotropic nature of the long-range stress response to local
perturbations such as force dipoles.

These results provide new tools and insight to unravel the complexity of
mechanical properties of amorphous solids. They also open new paths for
further exploration. For example, the new measure of minimum dipole stiffness
we propose is potentially a new direct probe to analyze soft spots in
amorphous solids. Relating this measure with existing characterizations of
soft spots, as well as yielding of amorphous solids, will be interesting
topics for future studies.

Moreover, the ability to separately vary the prestress and the microscopic
bond configuration in our theory offers a new perspective to analyze the
dynamics of amorphous solids under strain. In real amorphous solids, from
granular matter to colloids, particle-particle interactions are not simple
harmonic springs. Instead, they evolve under stress in complicated and often
unpredictable ways—micro-fracturing, crushing, sintering. This evolution is
further complicated by the role of the solvent, where contacts between
particles can be lubricated or frictional. Therefore, _stress distributions
can significantly change with little change in the configuration_. Our methods
provide a convenient set of tools to examine how the mechanical response of
these amorphous materials may evolve as the system explores an ensemble of
states where only internal stress is changed.

Prestress records the preparation history in amorphous solids and has an
important role in directing the mechanical response beyond the linear regime.
The implications of prestress and preparation history, in fact, have been
highlighted with respect to the development of flow inhomogeneities upon
yielding Vasisht  _et al._ (2020) and in connection with the fundamental
physics mechanisms controlling brittle or ductile yielding phenomena Ozawa
_et al._ (2018); Singh  _et al._ (2020); Barlow  _et al._ (2020). Our methods,
therefore, have wider relevance as they can be applied beyond the linear
response to investigate, through the SSSs characteristics, the role of
prestress in plasticity and yielding. The new tools discussed here can
potentially close the feedback loop between structure and stress, opening a
new pathway to study how amorphous solids yield and solidify under external
strain.

Finally, friction plays a central role in the complex behaviors of exotic
rigid states emerging in granular matter and dense suspensions, leading to
fascinating phenomena from shear jamming to discontinuous shear thickening Bi
_et al._ (2011); Seto  _et al._ (2013); Wyart and Cates (2014); Lin  _et al._
(2016); Guy  _et al._ (2018); Liu  _et al._ (2021b); Ramaswamy  _et al._
(2021). Interestingly, friction also lives in the $t^{\perp}$ channel of
stress, and a linear relation with $e^{\perp}$ can also be assumed, hence it
would be a natural next step to include friction into the SSSs description.
However, the Coulomb threshold imposes an upper limit to $t^{\perp}$ that is
not part of the methods developed here. How the frictional and the prestress
contributions to $t^{\perp}$ interplay and affect the macroscopic dynamics of
the material, will be of great interest for future studies.

###### Acknowledgements.

This work is supported in part by the National Science Foundation under grant
number NSF-DMR-2026825 (E.D.G and X.M.), NSF-EFRI-1741618 (S.Z., L.Z., and
X.M.), and the Office of Naval Research under grant number MURI
N00014-20-1-2479 (E.S. and X.M.)

##  Appendix A SSSs formulation for shear response of prestressed networks

###  A.1 Projection of a shear load to the SSSs linear space

When a mechanical network is subject to a shear load, the resulting bond
extensions and rotations can be written as the sum of contributions from the
affine shear field $e_{\textrm{affine}}$ and the nonaffine displacements
$\mathbb{C}|u_{\text{rsp}}\rangle$,

$\displaystyle|e\rangle$
$\displaystyle=|e_{\textrm{affine}}\rangle+\mathbb{C}|u_{\text{rsp}}\rangle,$
(47)

which is similar to Eq. (33) in the main text, but here we include both the
parallel and the perpendicular components of $e$.

When force-balance is reached, net force on each particle vanishes,

$\displaystyle-|f\rangle$
$\displaystyle=\mathbb{Q}|t\rangle=\mathbb{Q}\mathbb{K}(|e_{\text{affine}}\rangle+\mathbb{C}|u_{\text{rsp}}\rangle)=0.$
(48)

This means that
$|t\rangle=\mathbb{K}(|e_{\text{affine}}\rangle+\mathbb{C}|u_{\text{rsp}}\rangle)$
must be a vector that belongs to the null space of $\mathbb{Q}$. Thus it can
be written as a linear combinations of the SSSs of the system.

To facilitate the discussion of this SSSs linear combination, we define the
following notations. Let
$\\{\displaystyle\vec{t}_{s}^{~{}(1)},\ldots,\vec{t}_{s}^{~{}({N_{\textrm{S}}})}\\}$
be an orthonormal basis of the null space of $\mathbb{Q}$, and let
${\displaystyle P^{Q}_{s}}$ denote the $N_{b}d\times N_{\textrm{S}}$ matrix
whose columns are
$\displaystyle\vec{t}_{s}^{~{}(1)},\ldots,\vec{t}_{s}^{~{}({N_{\textrm{S}}})}$,
i.e., ${\displaystyle
P^{Q}_{s}={\begin{bmatrix}\vec{t}_{s}^{~{}(1)},\ldots,\vec{t}_{s}^{~{}({N_{\textrm{S}}})}\end{bmatrix}}}$.
One can also define the $N_{b}d\times(N_{b}d-N_{\textrm{S}})$ dimensional
matrix ${\displaystyle P^{Q}_{r}}$ whose columns are an orthonormal basis of
the orthogonal compliment of the null space of $\mathbb{Q}$. Similarly one can
define the $Nd\times N_{\textrm{0}}$ matrix ${\displaystyle P^{C}_{s}}$ whose
columns are an orthonormal basis of the null space of $\mathbb{C}$ (ZMs), and
the $Nd\times(Nd-N_{\textrm{0}})$ matrix ${\displaystyle P^{C}_{r}}$ whose
columns are an orthonormal basis of the orthogonal compliment of the null
space of $\mathbb{C}$. These matrices are represented as,

$\displaystyle P^{Q}_{s}$
$\displaystyle={\begin{bmatrix}\vec{t}_{s}^{~{}(1)},\ldots,\vec{t}_{s}^{~{}({N_{\textrm{S}}})}\end{bmatrix}},$
(49) $\displaystyle P^{Q}_{r}$
$\displaystyle={\begin{bmatrix}\vec{t}_{r}^{~{}(1)},\ldots,\vec{t}_{r}^{~{}({N_{b}d-N_{\textrm{S}}})}\end{bmatrix}},$
(50) $\displaystyle P^{C}_{s}$
$\displaystyle={\begin{bmatrix}\vec{u}_{s}^{~{}(1)},\ldots,\vec{u}_{s}^{~{}({N_{\textrm{0}}})}\end{bmatrix}},$
(51) $\displaystyle P^{C}_{r}$
$\displaystyle={\begin{bmatrix}\vec{u}_{r}^{~{}(1)},\ldots,\vec{u}_{r}^{~{}({Nd-
N_{\textrm{0}}})}\end{bmatrix}},$ (52)

and they satisfy the following identities,

$\displaystyle(P^{Q}_{s})^{T}\cdot P^{Q}_{s}$
$\displaystyle=\mathbb{I}_{(N_{\textrm{S}})},$ (53)
$\displaystyle(P^{Q}_{r})^{T}\cdot P^{Q}_{r}$
$\displaystyle=\mathbb{I}_{(N_{b}d-N_{\textrm{S}})},$ (54)
$\displaystyle(P^{C}_{s})^{T}\cdot P^{C}_{s}$
$\displaystyle=\mathbb{I}_{(N_{0})},$ (55) $\displaystyle(P^{C}_{r})^{T}\cdot
P^{C}_{r}$ $\displaystyle=\mathbb{I}_{(Nd-N_{0})},$ (56) $\displaystyle
P^{Q}_{s}\cdot(P^{Q}_{s})^{T}+P^{Q}_{r}\cdot(P^{Q}_{r})^{T}$
$\displaystyle=\mathbb{I}_{(N_{b}d)},$ (57) $\displaystyle
P^{C}_{s}\cdot(P^{C}_{s})^{T}+P^{C}_{r}\cdot(P^{C}_{r})^{T}$
$\displaystyle=\mathbb{I}_{(Nd)},$ (58)

where $\mathbb{I}_{(N_{\textrm{S}})}$ is an identity matrix of dimension
$N_{\textrm{S}}$, and other definitions follows similarly.

As we discussed above, $|t\rangle$ is a linear combination of the SSSs,

$\displaystyle|t\rangle$
$\displaystyle=\sum_{i}^{N_{\textrm{S}}}\alpha_{i}|t_{\textrm{SSS,
i}}\rangle=P_{s}^{Q}\cdot\vec{\alpha},$ (59)

where $\vec{\alpha}$ are coefficients of the linear combination of $|t\rangle$
as the SSSs.

Because the basis we use are orthonormal,

$\displaystyle\vec{\alpha}$
$\displaystyle=(P^{Q}_{s})^{T}|t\rangle=(P^{Q}_{s})^{T}\mathbb{K}(|e_{\text{affine}}\rangle+\mathbb{C}|u_{\text{rsp}}\rangle),$
(60) $\displaystyle\vec{0}$
$\displaystyle=(P^{Q}_{r})^{T}|t\rangle=(P^{Q}_{r})^{T}\mathbb{K}(|e_{\text{affine}}\rangle+\mathbb{C}|u_{\text{rsp}}\rangle),$
(61)

inserting identity matrix [Eq. (57)] into Eq. (61) and using the fact that,

$\displaystyle\mathbb{Q}\cdot P_{s}^{Q}=\vec{0},$

we have,

$\displaystyle\vec{0}$
$\displaystyle=(P^{Q}_{r})^{T}\mathbb{K}\cdot[P^{Q}_{s}\cdot(P^{Q}_{s})^{T}+P^{Q}_{r}\cdot(P^{Q}_{r})^{T}]\cdot(|e_{\text{affine}}\rangle+\mathbb{C}|u_{\text{rsp}}\rangle)$
$\displaystyle=(P^{Q}_{r})^{T}\mathbb{K}\cdot[P^{Q}_{s}\cdot(P^{Q}_{s})^{T}\cdot\mathbb{Q}^{T}+P^{Q}_{r}\cdot(P^{Q}_{r})^{T}\cdot\mathbb{Q}^{T}]\cdot|u_{\text{rsp}}\rangle$
$\displaystyle~{}~{}~{}+(P^{Q}_{r})^{T}\mathbb{K}\cdot[P^{Q}_{s}\cdot(P^{Q}_{s})^{T}+P^{Q}_{r}\cdot(P^{Q}_{r})^{T}]\cdot|e_{\text{affine}}\rangle$
$\displaystyle=(P^{Q}_{r})^{T}\mathbb{K}\cdot
P^{Q}_{r}\cdot(P^{Q}_{r})^{T}\cdot\mathbb{Q}^{T}\cdot|u_{\text{rsp}}\rangle+(P^{Q}_{r})^{T}\mathbb{K}\cdot[P^{Q}_{s}\cdot(P^{Q}_{s})^{T}+P^{Q}_{r}\cdot(P^{Q}_{r})^{T}]\cdot|e_{\text{affine}}\rangle$
$\displaystyle=\mathbb{K}_{rr}(P^{Q}_{r})^{T}\mathbb{Q}^{T}\cdot|u_{\text{rsp}}\rangle+[\mathbb{K}_{rs}(P^{Q}_{s})^{T}+\mathbb{K}_{rr}(P^{Q}_{r})^{T}]\cdot|e_{\text{affine}}\rangle$
$\displaystyle\implies(P^{Q}_{r})^{T}\mathbb{Q}^{T}\cdot|u_{\text{rsp}}\rangle=-(\mathbb{K}_{rr})^{-1}[\mathbb{K}_{rs}(P^{Q}_{s})^{T}+\mathbb{K}_{rr}(P^{Q}_{r})^{T}]\cdot|e_{\text{affine}}\rangle,$
(62)

where we defined the decomposition of $\mathbb{K}$ into the null and
orthogonal compliment space as

$\displaystyle\mathbb{K}\rightarrow\begin{pmatrix}(P^{Q}_{s})^{T}\cdot\mathbb{K}\cdot
P^{Q}_{s}&(P^{Q}_{s})^{T}\cdot\mathbb{K}\cdot P^{Q}_{r}\\\
(P^{Q}_{r})^{T}\cdot\mathbb{K}\cdot
P^{Q}_{s}&(P^{Q}_{r})^{T}\cdot\mathbb{K}\cdot P^{Q}_{r}\\\
\end{pmatrix}=\begin{pmatrix}\mathbb{K}_{ss}&\mathbb{K}_{sr}\\\
\mathbb{K}_{rs}&\mathbb{K}_{rr}\\\ \end{pmatrix}$ (63)

and also used the fact that $\mathbb{K}_{rr}$ is invertible.

The coefficients $\vec{\alpha}$ in Eq. (60) can then be solved as

$\displaystyle\vec{\alpha}$
$\displaystyle=(P^{Q}_{s})^{T}\mathbb{K}\cdot[P^{Q}_{s}\cdot(P^{Q}_{s})^{T}+P^{Q}_{r}\cdot(P^{Q}_{r})^{T}]\cdot(\mathbb{Q}^{T}|u_{\text{rsp}}\rangle+|e_{\text{affine}}\rangle)$
$\displaystyle=(P^{Q}_{s})^{T}\mathbb{K}\cdot
P^{Q}_{r}\cdot(P^{Q}_{r})^{T}\cdot\mathbb{Q}^{T}\cdot|u_{\text{rsp}}\rangle$
$\displaystyle~{}~{}~{}+(P^{Q}_{s})^{T}\mathbb{K}\cdot[P^{Q}_{s}\cdot(P^{Q}_{s})^{T}+P^{Q}_{r}\cdot(P^{Q}_{r})^{T}]\cdot|e_{\text{affine}}\rangle$
$\displaystyle=\mathbb{K}_{sr}(P^{Q}_{r})^{T}\mathbb{Q}^{T}\cdot|u_{\text{rsp}}\rangle+[\mathbb{K}_{ss}(P^{Q}_{s})^{T}+\mathbb{K}_{sr}(P^{Q}_{r})^{T}]\cdot|e_{\text{affine}}\rangle.$

Plug in Equation (62) to eliminate $u_{\text{rsp}}$,

$\displaystyle\vec{\alpha}$
$\displaystyle=\mathbb{K}_{sr}\\{-(\mathbb{K}_{rr})^{-1}[\mathbb{K}_{rs}(P^{Q}_{s})^{T}+\mathbb{K}_{rr}(P^{Q}_{r})^{T}]\cdot|e_{\text{affine}}\rangle\\}+[\mathbb{K}_{ss}(P^{Q}_{s})^{T}+\mathbb{K}_{sr}(P^{Q}_{r})^{T}]\cdot|e_{\text{affine}}\rangle$
$\displaystyle=\\{\mathbb{K}_{sr}(P^{Q}_{r})^{T}+\mathbb{K}_{ss}(P^{Q}_{s})^{T}-\mathbb{K}_{sr}(\mathbb{K}_{rr})^{-1}[\mathbb{K}_{rs}(P^{Q}_{s})^{T}+\mathbb{K}_{rr}(P^{Q}_{r})^{T}]\\}\cdot|e_{\text{affine}}\rangle$
$\displaystyle=[\mathbb{K}_{ss}(P^{Q}_{s})^{T}-\mathbb{K}_{sr}(\mathbb{K}_{rr})^{-1}\mathbb{K}_{rs}(P^{Q}_{s})^{T}]\cdot|e_{\text{affine}}\rangle$
$\displaystyle=[\mathbb{K}_{ss}-\mathbb{K}_{sr}(\mathbb{K}_{rr})^{-1}\mathbb{K}_{rs}]\cdot(P^{Q}_{s})^{T}\cdot|e_{\text{affine}}\rangle.$
(64)

This can be further simplified by letting $\mathbb{A}=\mathbb{K}^{-1}$, and
decompose $\mathbb{A}$ into the column-space and null-space of $\mathbb{Q}$
as,

$\displaystyle\mathbb{A}\rightarrow\begin{pmatrix}(P^{Q}_{s})^{T}\cdot\mathbb{A}\cdot
P^{Q}_{s}&(P^{Q}_{s})^{T}\cdot\mathbb{A}\cdot P^{Q}_{r}\\\
(P^{Q}_{r})^{T}\cdot\mathbb{A}\cdot
P^{Q}_{s}&(P^{Q}_{r})^{T}\cdot\mathbb{A}\cdot P^{Q}_{r}\\\
\end{pmatrix}=\begin{pmatrix}\mathbb{A}_{ss}&\mathbb{A}_{sr}\\\
\mathbb{A}_{rs}&\mathbb{A}_{rr}\\\ \end{pmatrix}.$ (65)

One can see that,

$\displaystyle\mathbb{K}_{ss}\cdot\mathbb{A}_{ss}+\mathbb{K}_{sr}\cdot\mathbb{A}_{rs}$
$\displaystyle=\mathbb{I}_{ss},$ (66)
$\displaystyle\mathbb{K}_{rs}\cdot\mathbb{A}_{ss}+\mathbb{K}_{rr}\cdot\mathbb{A}_{rs}$
$\displaystyle=\vec{0}_{rs}.$ (67)

Right multiply by $(\mathbb{A}_{ss})^{-1}$ on both sides of two equations,

$\displaystyle\mathbb{K}_{ss}+\mathbb{K}_{sr}\cdot\mathbb{A}_{rs}\cdot(\mathbb{A}_{ss})^{-1}=(\mathbb{A}_{ss})^{-1},$
(68)
$\displaystyle\mathbb{K}_{rs}=-\mathbb{K}_{rr}\cdot\mathbb{A}_{rs}\cdot(\mathbb{A}_{ss})^{-1}.$
(69)

Then,

$\displaystyle\mathbb{K}_{ss}+\mathbb{K}_{sr}\cdot\mathbb{A}_{rs}\cdot(\mathbb{A}_{ss})^{-1}=(\mathbb{A}_{ss})^{-1},$
(70)
$\displaystyle-(\mathbb{K}_{rr})^{-1}\cdot\mathbb{K}_{rs}=\mathbb{A}_{rs}\cdot(\mathbb{A}_{ss})^{-1}.$
(71)

Combining these two equations we have,

$\displaystyle((\mathbb{K}^{-1})_{ss})^{-1}$
$\displaystyle=(\mathbb{A}_{ss})^{-1}$ (72)
$\displaystyle=\mathbb{K}_{ss}-\mathbb{K}_{sr}\cdot(\mathbb{K}_{rr})^{-1}\cdot\mathbb{K}_{rs}.$
(73)

As a result, $\vec{\alpha}$ is simplified to

$\displaystyle\vec{\alpha}=((\mathbb{K}^{-1})_{ss})^{-1}\cdot(P^{Q}_{s})^{T}\cdot|e_{\text{affine}}\rangle$
(74)

and the tension response to this external shear is

$\displaystyle|t\rangle=P^{Q}_{s}\cdot((\mathbb{K}^{-1})_{ss})^{-1}\cdot(P^{Q}_{s})^{T}\cdot|e_{\text{affine}}\rangle.$
(75)

Note that this $t$ includes both $t^{\parallel}$ and $t^{\perp}$, and this
formulation applies to other types of homogeneous strain, such as hydrostatic
compression, as well.

###  A.2 Affine bond deformation in prestressed systems

In this subsection we derive the $e_{\textrm{affine}}$ field for any external
load represented by a strain tensor $\epsilon$.

We start from the (affine) deformation gradient

$\displaystyle\Lambda_{ij}\equiv\frac{\partial R_{i}}{\partial R_{0,j}}$ (76)

where the strain tensor $\epsilon=(\Lambda^{T}\Lambda-I)/2$. The affinely
deformed positions of each particle are then

$\displaystyle\vec{R}_{\ell,\textrm{affine}}=\Lambda\cdot\vec{R}_{\ell,0}.$
(77)

We can then use the formulation discussed in Sec. II [Eq.(4)] to calculate the
affine extensions/rotations
$e_{b,\text{affine}}^{\parallel},e_{b,\text{affine}}^{\perp}$ for each bond,
which consist $|e_{\text{affine}}\rangle$.

##  Appendix B Dipole stiffness $\kappa$ in prestressed systems

When a pair of dipole forces is applied on a mechanical network between two
particles that belong to the same rigid cluster, the network will show a
linear response with tension distributed on the bonds. In this Appendix we
derive the stress field of a prestressed network in response to the force
dipole, and obtain a computationally efficient formula that gives the
stiffness the system has against this force dipole.

###  B.1 Local dipole stiffness

As discussed in Sec. III.5, the sum of the external force dipole and the
tension response, $|t_{\textrm{dipole}}\rangle+|t_{\textrm{rsp}}\rangle$, must
be a SSS of the prestressed network,

$\displaystyle|t_{\textrm{dipole}}\rangle+|t_{\textrm{rsp}}\rangle=\sum_{i}^{N_{\textrm{SSS}}}\alpha_{i}|t_{\textrm{SSS,
i}}\rangle.$

The coefficients $\alpha_{i}$ are determined in the same way as discussed in
App. A, just by replacing $|e_{\textrm{aff}}\rangle$ with $|b\rangle$. As a
result

$\displaystyle|t_{\textrm{dipole}}\rangle+|t_{\textrm{rsp}}\rangle=P^{Q}_{s}\cdot((\mathbb{K}^{-1})_{ss})^{-1}\cdot(P^{Q}_{s})^{T}\cdot|b\rangle.$

We can thus use this in the expression for the dipole stiffness [Eq. (41)],
where the denominator is now

$\displaystyle\langle{f_{\textrm{dipole}}|u_{\textrm{rsp}}}\rangle$
$\displaystyle=-\langle
b|\mathbb{K}\mathbb{C}|u_{\textrm{rsp}}\rangle=-\langle
b|t_{\textrm{rsp}}\rangle$
$\displaystyle=k_{b}-\displaystyle\sum_{i,j}^{N_{\textrm{SSS}}}\langle
b|t_{\textrm{SSS},i}\rangle[(\mathbb{K}^{-1})_{ss}]^{-1}\langle
t_{\textrm{SSS},j}|b\rangle,$ (78)

where we used Eq. (38) and the fact that
$|t_{\textrm{rsp}}\rangle=\mathbb{K}\mathbb{C}|u_{\textrm{rsp}}\rangle$ (bond
extension causes tension).

Therefore the dipole stiffness is

$\displaystyle\kappa_{b}$
$\displaystyle=\frac{\langle{f_{\textrm{dipole}}|f_{\textrm{dipole}}}\rangle}{\langle{f_{\textrm{dipole}}|u_{\textrm{rsp}}}\rangle}$
$\displaystyle=\frac{\langle
b|\mathbb{K}\mathbb{C}\mathbb{Q}\mathbb{K}|b\rangle}{\langle
b|\mathbb{K}|b\rangle-\langle
b|\displaystyle\sum_{i,j}^{N_{\textrm{SSS}}}|t_{\textrm{SSS},i}\rangle[(\mathbb{K}^{-1})_{ss}]^{-1}\langle
t_{\textrm{SSS},j}|b\rangle}$ (79)
$\displaystyle=\frac{2k_{b}^{2}}{k_{b}-\displaystyle\sum_{i,j}^{N_{\textrm{SSS}}}\langle
b|t_{\textrm{SSS},i}\rangle[(\mathbb{K}^{-1})_{ss}]^{-1}\langle
t_{\textrm{SSS},j}|b\rangle},$

where the factor of 2 in the last line comes from the fact that
$f_{\textrm{dipole}}$ involves 2 particles.

###  B.2 Non-local dipole stiffness

Besides applying the force dipole on an arbitrary existing bond $b$ in the
system, one could also apply a force dipole between two particles which are
not connected. As we discuss below, this pair of dipole introduces a new
constraint into the system.

Force balance with imposed force dipole $|f_{\textrm{dipole}}\rangle$ can be
written as:

$\displaystyle 0$
$\displaystyle=|f_{\textrm{dipole}}\rangle+|f_{\textrm{rsp}}\rangle=|f_{\textrm{dipole}}\rangle-\mathbb{Q}\mathbb{K}\mathbb{C}|u_{\textrm{rsp}}\rangle,$
(80)

where $|u_{\textrm{rsp}}\rangle$ indicates the particle displacements in
response to the force dipole as the system reaches force balance. Unlike the
local dipole case discussed above, here $f_{\textrm{dipole}}$ can not be
written as $\mathbb{Q}|t_{\textrm{dipole}}\rangle$ on an existing bond, as the
two sites are not connected. Instead, we can introduce an auxiliary bond
between the two sites that carry $|f_{\textrm{dipole}}\rangle$. In this sense,
a new constraint and thus a new SSS is added to the system by the auxiliary
bond.

Here we introduce the new $\mathbb{Q},\mathbb{C}$ matrices after introducing
the auxiliary bond (which has zero spring constant so that it will not induce
tension responses) as

$\displaystyle\tilde{\mathbb{C}}=\begin{pmatrix}\mathbb{C}\\\
C_{a}\end{pmatrix},$
$\displaystyle\tilde{\mathbb{Q}}=\begin{pmatrix}\mathbb{Q}&Q_{a}\end{pmatrix},$
$\displaystyle\tilde{\mathbb{K}}=\begin{pmatrix}\mathbb{K}&\\\
&0\end{pmatrix},$
$\displaystyle\tilde{\mathbb{C}}\cdot|u\rangle=\begin{pmatrix}|e\rangle\\\
e_{a}\end{pmatrix},$
$\displaystyle\tilde{\mathbb{Q}}\cdot\begin{pmatrix}|t\rangle\\\
t_{a}\end{pmatrix}=-|f\rangle.$

The dimension of the bond space is extended with one additional component from
the auxiliary bond indexed as $a$.

Now the force-balanced total tension can be written as

$\displaystyle|\tilde{t}\rangle=\begin{pmatrix}|t_{\textrm{rsp}}\rangle\\\
t_{\textrm{dipole}}\end{pmatrix}=\begin{pmatrix}\mathbb{K}\mathbb{C}|u_{\textrm{rsp}}\rangle\\\
t_{\textrm{dipole}}\end{pmatrix}=\tilde{\mathbb{K}}\tilde{\mathbb{C}}|u_{\textrm{rsp}}\rangle+|\tilde{t}_{\textrm{dipole}}\rangle,$
(81)

where $t_{\textrm{dipole}}$ lives on the auxillary bond. This
$|\tilde{t}\rangle$ satisfies force balance and must be a SSS of the network
with the auxillary bond.

Similar to App. A, one can define
$P_{s}^{\tilde{Q}},P_{r}^{\tilde{Q}},P_{s}^{\tilde{C}},P_{r}^{\tilde{Q}}$ to
the new matrices. We can then decompose the total tension onto SSSs,

$\displaystyle|\tilde{t}\rangle$
$\displaystyle=\sum_{i}^{N_{\textrm{SSS}}}\alpha_{i}|t_{\textrm{SSS,
i}}\rangle=P_{s}^{\tilde{Q}}\cdot\vec{\alpha},$

where $\vec{\alpha}$’s are coefficients of the linear combination of the SSSs.

We can compute the tension in a similar way,

$\displaystyle(P^{\tilde{Q}}_{s})^{T}|\tilde{t}\rangle$
$\displaystyle=\vec{\alpha},$
$\displaystyle(P^{\tilde{Q}}_{r})^{T}|\tilde{t}\rangle$
$\displaystyle=\vec{0}.$

Thus

$\displaystyle\vec{0}$
$\displaystyle=(P^{\tilde{Q}}_{r})^{T}\left(\tilde{\mathbb{K}}\tilde{\mathbb{C}}|u_{\textrm{rsp}}\rangle+|\tilde{t}_{\textrm{dipole}}\rangle\right)$
(82)
$\displaystyle=(P^{\tilde{Q}}_{r})^{T}\tilde{\mathbb{K}}\left(P^{\tilde{Q}}_{r}(P^{\tilde{Q}}_{r})^{T}+P^{\tilde{Q}}_{s}(P^{\tilde{Q}}_{s})^{T}\right)\tilde{\mathbb{Q}}^{T}|u_{\textrm{rsp}}\rangle$
$\displaystyle\quad+(P^{\tilde{Q}}_{r})^{T}|\tilde{t}_{\textrm{dipole}}\rangle$
(83)
$\displaystyle=(P^{\tilde{Q}}_{r})^{T}\tilde{\mathbb{K}}P^{\tilde{Q}}_{r}(P^{\tilde{Q}}_{r})^{T}\tilde{\mathbb{Q}}^{T}|u_{\textrm{rsp}}\rangle+(P^{\tilde{Q}}_{r})^{T}|\tilde{t}_{\textrm{dipole}}\rangle$
(84)
$\displaystyle=\tilde{\mathbb{K}}_{rr}(P^{\tilde{Q}}_{r})^{T}\tilde{\mathbb{Q}}^{T}|u_{\textrm{rsp}}\rangle+(P^{\tilde{Q}}_{r})^{T}|\tilde{t}_{\textrm{dipole}}\rangle.$
(85)

As a result,

$(P^{\tilde{Q}}_{r})^{T}\tilde{\mathbb{Q}}^{T}|u_{\textrm{rsp}}\rangle=-\left(\tilde{\mathbb{K}}_{rr}\right)^{-1}(P^{\tilde{Q}}_{r})^{T}|\tilde{t}_{\textrm{dipole}}\rangle.$
(86)

Note that although the spring constant of the auxiliary bond is 0, the matrix
$\tilde{\mathbb{K}}_{rr}$ is still invertable. This is because the auxiliary
bond is a redundant bond (otherwise the network would yield, resulting in no
stress in linear response), which increases the dimension of the SSSs space
and does not introduce new vectors in the orthogonal compliment space. As a
result, $\tilde{\mathbb{K}}_{rr}=\mathbb{K}_{rr}$ and is invertable.

This can be used to find the coefficients for the SSSs

$\displaystyle\vec{\alpha}$
$\displaystyle=(P^{\tilde{Q}}_{s})^{T}|\tilde{t}\rangle$
$\displaystyle=(P^{\tilde{Q}}_{s})^{T}\left(\tilde{\mathbb{K}}\tilde{\mathbb{C}}|u_{\textrm{rsp}}\rangle+|\tilde{t}_{\textrm{dipole}}\rangle\right)$
$\displaystyle=\tilde{\mathbb{K}}_{sr}(P^{\tilde{Q}}_{r})^{T}\tilde{\mathbb{Q}}^{T}|u_{\textrm{rsp}}\rangle+(P^{\tilde{Q}}_{s})^{T}|\tilde{t}_{\textrm{dipole}}\rangle$
$\displaystyle=\left((P^{\tilde{Q}}_{s})^{T}-\tilde{\mathbb{K}}_{sr}\left(\tilde{\mathbb{K}}_{rr}\right)^{-1}(P^{\tilde{Q}}_{r})^{T}\right)|\tilde{t}_{\textrm{dipole}}\rangle,$
(87)

where in the last line we used Eq. (86). These coefficients gives the tension
field in response to a nonlocal dipole.

We can then proceed to calculate the dipole stiffness. Applying
$P^{\tilde{Q}}_{r}$ on both sides of Eq. (86), we have the left hand side,

$\displaystyle
P^{\tilde{Q}}_{r}(P^{\tilde{Q}}_{r})^{T}\tilde{\mathbb{Q}}^{T}|u\rangle$
$\displaystyle=\left(\mathbb{I}-P^{\tilde{Q}}_{s}(P^{\tilde{Q}}_{s})^{T}\right)\tilde{\mathbb{Q}}^{T}|u\rangle$
(88)
$\displaystyle=\tilde{\mathbb{Q}}^{T}|u\rangle-P^{\tilde{Q}}_{s}\left(\tilde{\mathbb{Q}}\cdot
P^{\tilde{Q}}_{s}\right)^{T}|u\rangle$ (89)
$\displaystyle=\tilde{\mathbb{Q}}^{T}|u\rangle.$ (90)

Equaling it to the right hand side we have

$\displaystyle\tilde{\mathbb{Q}}^{T}|u\rangle=-P^{\tilde{Q}}_{r}\left(\tilde{\mathbb{K}}_{rr}\right)^{-1}(P^{\tilde{Q}}_{r})^{T}|\tilde{t}_{\textrm{dipole}}\rangle,$
(91)

expressing the displacement field as a function of the imposed dipole.

One can then compute the non-local dipole stiffness $\kappa_{a}$ as (the force
dipole
$|f_{\textrm{dipole}}\rangle=-\mathbb{\tilde{Q}}|\tilde{t}_{\textrm{dipole}}\rangle\equiv-\mathbb{\tilde{Q}}|a\rangle$
where $|a\rangle$ is the vector in the labeling space of bonds which has zeros
in all bonds and unity on the $a$-th component (the auxiliary bond), and a
spring constant of unity is added (which is only used to represent the
external force dipole, where the actual spring constant of the auxillary bond
regarding its contribution to the network response is still 0 as discussed
above),

$\displaystyle\kappa_{a}$
$\displaystyle=\frac{\langle{f_{\textrm{dipole}}|f_{\textrm{dipole}}}\rangle}{\langle{f_{\textrm{dipole}}|u}\rangle}$
(92)
$\displaystyle=\frac{\langle{\tilde{t}_{\textrm{dipole}}|\mathbb{\tilde{Q}}^{T}\mathbb{\tilde{Q}}|\tilde{t}_{\textrm{dipole}}}\rangle}{-\langle{\tilde{t}_{\textrm{dipole}}|\mathbb{\tilde{Q}}^{T}|u}\rangle}$
(93)
$\displaystyle=\frac{\langle{\tilde{t}_{\textrm{dipole}}|\mathbb{\tilde{Q}}^{T}\mathbb{\tilde{Q}}|\tilde{t}_{\textrm{dipole}}}\rangle}{\langle{\tilde{t}_{\textrm{dipole}}|P^{\tilde{Q}}_{r}\left(\tilde{\mathbb{K}}_{rr}\right)^{-1}(P^{\tilde{Q}}_{r})^{T}|\tilde{t}_{\textrm{dipole}}}\rangle}$
(94)
$\displaystyle=\frac{\langle{a|\mathbb{\tilde{Q}}^{T}\mathbb{\tilde{Q}}|a}\rangle}{\langle{a|P^{\tilde{Q}}_{r}\left(\tilde{\mathbb{K}}_{rr}\right)^{-1}(P^{\tilde{Q}}_{r})^{T}|a}\rangle}$
(95)
$\displaystyle=\frac{2}{\langle{a|P^{\tilde{Q}}_{r}\left(\tilde{\mathbb{K}}_{rr}\right)^{-1}(P^{\tilde{Q}}_{r})^{T}|a}\rangle}$
(96)

The local force dipole response is a special case of non-local force dipole
response. When considering local force dipoles, the auxiliary bond $a$
overlaps with bond $b$. One can show that in this case, the non-local force
dipole stiffness reduces to local force dipole stiffness.

###  B.3 Minimum local dipole stiffness $\kappa_{b,\textrm{min}}$ in
prestressed systems

In this section we derive the minimum dipole stiffness on any given bond $b$,
which is the lowest with respect to the combination of $\parallel$ and $\perp$
directions. To do this we use $(\theta,\phi)$ to denote the direction of the
dipole forces, and the resulting $t_{\textrm{dipole}}$ can be written as

$\displaystyle|t_{\text{dipole}}\rangle$ $\displaystyle=\begin{pmatrix}0\\\
\vdots\\\ \cos{\theta}\\\ \sin{\theta}\cos{\phi}\\\ \sin{\theta}\sin{\phi}\\\
\vdots\\\ 0\\\ \end{pmatrix}\equiv|b\rangle$ (97)

where we defined the new $|b\rangle$ vector for this case allowing arbitrary
force directions. We can plug this $|b\rangle$ in Eq. (79) to find the dipole
stiffness for any given $(\theta,\phi)$, and represent the nontrivial term in
$\kappa_{b}$ as

$\displaystyle\eta_{b}$ $\displaystyle=\langle
b|\cdot\displaystyle\sum_{i,j}^{N_{\textrm{SSS}}}|t_{\textrm{SSS},i}\rangle[(\mathbb{K}^{-1})_{ss}]^{-1}\langle
t_{\textrm{SSS},j}|\cdot|b\rangle\equiv\langle b|\mathbb{P}|b\rangle$ (98)
$\displaystyle=\begin{pmatrix}0&\cdots&\cos{\theta}&\sin{\theta}\cos{\phi}&\sin{\theta}\sin{\phi}&\cdots&0\end{pmatrix}\cdot\mathbb{P}\cdot\begin{pmatrix}0\\\
\vdots\\\ \cos{\theta}\\\ \sin{\theta}\cos{\phi}\\\ \sin{\theta}\sin{\phi}\\\
\vdots\\\ 0\\\ \end{pmatrix}$ (99)
$\displaystyle=\begin{pmatrix}\cos{\theta}&\sin{\theta}\cos{\phi}&\sin{\theta}\sin{\phi}\\\
\end{pmatrix}\cdot\mathbb{P}_{b}\cdot\begin{pmatrix}\cos{\theta}\\\
\sin{\theta}\cos{\phi}\\\ \sin{\theta}\sin{\phi}\\\ \end{pmatrix}$ (100)

where $\mathbb{P}$ represents the matrix in between $|b\rangle$ vectors in Eq.
(98) and $\mathbb{P}_{b}$ is the $3\times 3$ part of $\mathbb{P}$ associated
with bond $b$.

Similarly,

$\displaystyle\langle b|\mathbb{K}\mathbb{C}\mathbb{Q}\mathbb{K}|b\rangle$
$\displaystyle=\begin{pmatrix}\cos{\theta}&\sin{\theta}\cos{\phi}&\sin{\theta}\sin{\phi}\\\
\end{pmatrix}\cdot{[\mathbb{K}\mathbb{C}\mathbb{Q}\mathbb{K}]}_{b}\cdot\begin{pmatrix}\cos{\theta}\\\
\sin{\theta}\cos{\phi}\\\ \sin{\theta}\sin{\phi}\\\ \end{pmatrix}$
$\displaystyle\langle b|\mathbb{K}|b\rangle$
$\displaystyle=\begin{pmatrix}\cos{\theta}&\sin{\theta}\cos{\phi}&\sin{\theta}\sin{\phi}\\\
\end{pmatrix}\cdot{[\mathbb{K}]}_{b}\cdot\begin{pmatrix}\cos{\theta}\\\
\sin{\theta}\cos{\phi}\\\ \sin{\theta}\sin{\phi}\\\ \end{pmatrix}$

where ${[\mathbb{K}\mathbb{C}\mathbb{Q}\mathbb{K}]}_{b}$ and
${[\mathbb{K}]}_{b}$ are the $3\times 3$ parts of
$\mathbb{K}\mathbb{C}\mathbb{Q}\mathbb{K}$ and $\mathbb{K}$ associated with
bond $b$, respectively.

As a result,

$\displaystyle\kappa_{b}$ $\displaystyle=\frac{\langle
b|\mathbb{K}\mathbb{C}\mathbb{Q}\mathbb{K}|b\rangle}{\langle
b|\mathbb{K}|b\rangle-\langle
b|\displaystyle\sum_{i,j}^{N_{\textrm{SSS}}}|t_{\textrm{SSS},i}\rangle[(\mathbb{K}^{-1})_{ss}]^{-1}\langle
t_{\textrm{SSS},j}|b\rangle}$
$\displaystyle=\frac{\begin{pmatrix}\cos{\theta}&\sin{\theta}\cos{\phi}&\sin{\theta}\sin{\phi}\\\
\end{pmatrix}\cdot{[\mathbb{K}\mathbb{C}\mathbb{Q}\mathbb{K}]}_{b}\cdot\begin{pmatrix}\cos{\theta}\\\
\sin{\theta}\cos{\phi}\\\ \sin{\theta}\sin{\phi}\\\
\end{pmatrix}}{\begin{pmatrix}\cos{\theta}&\sin{\theta}\cos{\phi}&\sin{\theta}\sin{\phi}\\\
\end{pmatrix}\cdot\\{{[\mathbb{K}]}_{b}-\mathbb{P}_{b}\\}\cdot\begin{pmatrix}\cos{\theta}\\\
\sin{\theta}\cos{\phi}\\\ \sin{\theta}\sin{\phi}\\\ \end{pmatrix}}.$

Minimizing $\kappa_{b}$ for a single bond $b$,

$\displaystyle\kappa_{b,\textrm{min}}$
$\displaystyle=\min_{\theta,\phi}\left[\frac{\begin{pmatrix}\cos{\theta}&\sin{\theta}\cos{\phi}&\sin{\theta}\sin{\phi}\\\
\end{pmatrix}\cdot{[\mathbb{K}\mathbb{C}\mathbb{Q}\mathbb{K}]}_{b}\cdot\begin{pmatrix}\cos{\theta}\\\
\sin{\theta}\cos{\phi}\\\ \sin{\theta}\sin{\phi}\\\
\end{pmatrix}}{\begin{pmatrix}\cos{\theta}&\sin{\theta}\cos{\phi}&\sin{\theta}\sin{\phi}\\\
\end{pmatrix}\cdot\\{{[\mathbb{K}]}_{b}-\mathbb{P}_{b}\\}\cdot\begin{pmatrix}\cos{\theta}\\\
\sin{\theta}\cos{\phi}\\\ \sin{\theta}\sin{\phi}\\\ \end{pmatrix}}\right]$
(101)

is an optimization problem with respect to the two variables $\theta$ and
$\phi$. To solve for such optimization problems in our system, we used the
Nelder-Mead simplex algorithm as described in Ref. Lagarias  _et al._ (1998).

By doing this minimization, we obtain $(\theta,\phi)$ as the softest dipole
stiffness direction for any bond $b$, which is typically transverse.

## References

  * Alexander (1998) S. Alexander, Physics reports 296, 65 (1998).
  * Ulm and Coussy (2003) F.-J. Ulm and O. Coussy,  _Mechanics and Durability of Solids. Solid mechanics_ (Pearson Education, 2003).
  * Wyart  _et al._ (2005) M. Wyart, L. E. Silbert, S. R. Nagel, and T. A. Witten, Phys. Rev. E 72, 051306 (2005).
  * Bouzid  _et al._ (2017) M. Bouzid, J. Colombo, L. V. Barbosa, and E. Del Gado, Nature Communications 8, 15846 (2017).
  * Bi  _et al._ (2011) D. Bi, J. Zhang, B. Chakraborty, and R. P. Behringer, Nature 480, 355–358 (2011).
  * Sehgal  _et al._ (2019) P. Sehgal, M. Ramaswamy, I. Cohen, and B. J. Kirby, Phys. Rev. Lett. 123, 128001 (2019).
  * Ballauff  _et al._ (2013) M. Ballauff, J. M. Brader, S. U. Egelhaaf, M. Fuchs, J. Horbach, N. Koumakis, M. Krüger, M. Laurati, K. J. Mutch, G. Petekidis, M. Siebenbürger, T. Voigtmann, and J. Zausch, Phys. Rev. Lett. 110, 215701 (2013).
  * Gei  _et al._ (2009) M. Gei, A. B. Movchan, and D. Bigoni, Journal of Applied Physics 105, 063507 (2009), https://doi.org/10.1063/1.3093694 .
  * Chen  _et al._ (2012) Z. Chen, Q. Guo, C. Majidi, W. Chen, D. J. Srolovitz, and M. P. Haataja, Phys. Rev. Lett. 109, 114302 (2012).
  * Fraternali  _et al._ (2014) F. Fraternali, G. Carpentieri, A. Amendola, R. E. Skelton, and V. F. Nesterenko, Applied Physics Letters 105, 201903 (2014), https://doi.org/10.1063/1.4902071 .
  * Meeussen  _et al._ (2020) A. S. Meeussen, E. C. Oğuz, M. van Hecke, and Y. Shokef, New Journal of Physics 22, 023030 (2020).
  * Merrigan  _et al._ (2021) C. Merrigan, C. Nisoli, and Y. Shokef, Phys. Rev. Research 3, 023174 (2021).
  * Armon  _et al._ (2011) S. Armon, E. Efrati, R. Kupferman, and E. Sharon, Science 333, 1726 (2011), https://science.sciencemag.org/content/333/6050/1726.full.pdf .
  * Wyczalkowski  _et al._ (2012) M. A. Wyczalkowski, Z. Chen, B. A. Filas, V. D. Varner, and L. A. Taber, Birth Defects Research Part C: Embryo Today: Reviews 96, 132 (2012), https://onlinelibrary.wiley.com/doi/pdf/10.1002/bdrc.21013 .
  * Feng  _et al._ (2018) D. Feng, J. Notbohm, A. Benjamin, S. He, M. Wang, L.-H. Ang, M. Bantawa, M. Bouzid, E. Del Gado, R. Krishnan, and M. R. Pollak, Proceedings of National Academy of Sciences 115, 1517 (2018).
  * Zhang  _et al._ (2017) L. Zhang, D. Z. Rocklin, L. M. Sander, and X. Mao, Physical Review Materials 1, 052602 (2017).
  * Ozawa  _et al._ (2018) M. Ozawa, L. Berthier, G. Biroli, A. Rosso, and G. Tarjus, Proceedings of the National Academy of Sciences 115, 6656 (2018), https://www.pnas.org/content/115/26/6656.full.pdf .
  * Vasisht  _et al._ (2020) V. V. Vasisht, G. Roberts, and E. Del Gado, Physical Review E 102, 010604 (2020).
  * Barlow  _et al._ (2020) H. J. Barlow, J. O. Cochran, and S. M. Fielding, Phys. Rev. Lett. 125, 168003 (2020).
  * Berthier  _et al._ (2019) E. Berthier, J. E. Kollmer, S. E. Henkes, K. Liu, J. M. Schwarz, and K. E. Daniels, Phys. Rev. Materials 3, 075602 (2019).
  * Debenedetti and Stillinger (2001) P. G. Debenedetti and F. H. Stillinger, Nature 410, 259 (2001).
  * Manning and Liu (2011) M. L. Manning and A. J. Liu, Phys. Rev. Lett. 107, 108302 (2011).
  * Ninarello  _et al._ (2017) A. Ninarello, L. Berthier, and D. Coslovich, Phys. Rev. X 7, 021039 (2017).
  * Vasisht  _et al._ (2021) V. V. Vasisht, P. Chaudhuri, and K. Martens, “Residual stress in athermal soft disordered solids: insights from microscopic and mesoscale models,” (2021), arXiv:2108.12782 [cond-mat.soft] .
  * Bhaumik  _et al._ (2021) H. Bhaumik, G. Foffi, and S. Sastry, Proceedings of the National Academy of Sciences 118 (2021), 10.1073/pnas.2100227118.
  * Richard  _et al._ (2020) D. Richard, M. Ozawa, S. Patinet, E. Stanifer, B. Shang, S. A. Ridout, B. Xu, G. Zhang, P. K. Morse, J.-L. Barrat, L. Berthier, M. L. Falk, P. Guan, A. J. Liu, K. Martens, S. Sastry, D. Vandembroucq, E. Lerner, and M. L. Manning, Phys. Rev. Materials 4, 113609 (2020).
  * Bouchaud (2002) J.-P. Bouchaud, arXiv preprint cond-mat/0211196 (2002).
  * Aben  _et al._ (2016) H. Aben, J. Anton, M. Õis, K. Viswanathan, S. Chandrasekar, and M. M. Chaudhri, Applied Physics Letters 109, 231903 (2016), https://doi.org/10.1063/1.4971339 .
  * Wen  _et al._ (2020) L. Wen, F. Pan, and X. Ding, Science Robotics 5, eabd9158 (2020).
  * DeGiuli  _et al._ (2014) E. DeGiuli, A. Laversanne-Finot, G. Düring, E. Lerner, and M. Wyart, Soft Matter 10, 5628 (2014).
  * Ellenbroek  _et al._ (2009) W. G. Ellenbroek, M. van Hecke, and W. van Saarloos, Phys. Rev. E 80, 061307 (2009).
  * Licup  _et al._ (2015) A. J. Licup, S. Münster, A. Sharma, M. Sheinman, L. M. Jawerth, B. Fabry, D. A. Weitz, and F. C. MacKintosh, Proceedings of the National Academy of Sciences 112, 9573 (2015), https://www.pnas.org/content/112/31/9573.full.pdf .
  * Feng  _et al._ (2015) J. Feng, H. Levine, X. Mao, and L. M. Sander, Physical Review E 91, 042710 (2015).
  * Feng  _et al._ (2016) J. Feng, H. Levine, X. Mao, and L. M. Sander, Soft matter 12, 1419 (2016).
  * Noll  _et al._ (2017) N. Noll, M. Mani, I. Heemskerk, S. J. Streichan, and B. I. Shraiman, Nature physics 13, 1221 (2017).
  * Hardin  _et al._ (2018) C. Hardin, J. Chattoraj, G. Manomohan, J. Colombo, T. Nguyen, D. Tambe, J. J. Fredberg, J. P. Birukov, K. Butler, E. Del Gado, and R. Krishnan, Biochem Biophys Res Commun 495, 749 (2018).
  * Yan and Bi (2019) L. Yan and D. Bi, Phys. Rev. X 9, 011029 (2019).
  * Merkel  _et al._ (2019) M. Merkel, K. Baumgarten, B. P. Tighe, and M. L. Manning, Proceedings of the National Academy of Sciences 116, 6560 (2019).
  * Liu  _et al._ (2021a) H. Liu, D. Zhou, L. Zhang, D. K. Lubensky, and X. Mao, arXiv preprint arXiv:2104.14743 (2021a).
  * Alvarado  _et al._ (2013) J. Alvarado, M. Sheinman, A. Sharma, F. C. MacKintosh, and G. H. Koenderink, Nature Physics 9, 591 (2013).
  * Ronceray  _et al._ (2016) P. Ronceray, C. P. Broedersz, and M. Lenz, Proceedings of the national academy of sciences 113, 2827 (2016).
  * Tanguy  _et al._ (2002a) A. Tanguy, J. P. Wittmer, F. Leonforte, and J.-L. Barrat, Phys. Rev. B 66, 174205 (2002a).
  * O’Hern  _et al._ (2003) C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E 68, 011306 (2003).
  * Zaccone and Scossa-Romano (2011) A. Zaccone and E. Scossa-Romano, Phys. Rev. B 83, 184205 (2011).
  * Liu and Nagel (2010) A. J. Liu and S. R. Nagel, Annu. Rev. Condens. Matter Phys. 1, 347 (2010).
  * Snoeijer  _et al._ (2004) J. H. Snoeijer, T. J. H. Vlugt, M. van Hecke, and W. van Saarloos, Phys. Rev. Lett. 92, 054302 (2004).
  * Tighe  _et al._ (2008) B. P. Tighe, A. R. T. van Eerd, and T. J. H. Vlugt, Phys. Rev. Lett. 100, 238001 (2008).
  * Tighe  _et al._ (2010) B. P. Tighe, J. H. Snoeijer, T. J. Vlugt, and M. van Hecke, Soft Matter 6, 2908 (2010).
  * Bi  _et al._ (2015) D. Bi, S. Henkes, K. E. Daniels, and B. Chakraborty, Annual Review of Condensed Matter Physics 6, 63 (2015), https://doi.org/10.1146/annurev-conmatphys-031214-014336 .
  * Edwards and Oakeshott (1989) S. F. Edwards and R. Oakeshott, Physica A: Statistical Mechanics and its Applications 157, 1080 (1989).
  * Mao  _et al._ (2007) X. Mao, P. M. Goldbart, X. Xing, and A. Zippelius, Europhysics Letters (EPL) 80, 26004 (2007).
  * Henkes and Chakraborty (2009) S. Henkes and B. Chakraborty, Phys. Rev. E 79, 061301 (2009).
  * Mao  _et al._ (2009) X. Mao, P. M. Goldbart, X. Xing, and A. Zippelius, Phys. Rev. E 80, 031140 (2009).
  * Lemaître (2017) A. Lemaître, Phys. Rev. E 96, 052101 (2017).
  * Maier  _et al._ (2017) M. Maier, A. Zippelius, and M. Fuchs, Phys. Rev. Lett. 119, 265701 (2017).
  * DeGiuli (2018) E. DeGiuli, Phys. Rev. Lett. 121, 118001 (2018).
  * Nampoothiri  _et al._ (2020) J. N. Nampoothiri, Y. Wang, K. Ramola, J. Zhang, S. Bhattacharjee, and B. Chakraborty, Phys. Rev. Lett. 125, 118002 (2020).
  * Calladine (1978) C. R. Calladine, International journal of solids and structures 14, 161 (1978).
  * Pellegrino and Calladine (1986) S. Pellegrino and C. R. Calladine, International Journal of Solids and Structures 22, 409 (1986).
  * Sun  _et al._ (2012) K. Sun, A. Souslov, X. Mao, and T. Lubensky, Proceedings of the National Academy of Sciences 109, 12369 (2012).
  * Kane and Lubensky (2014) C. Kane and T. Lubensky, Nature Physics 10, 39 (2014).
  * Lubensky  _et al._ (2015) T. Lubensky, C. Kane, X. Mao, A. Souslov, and K. Sun, Reports on Progress in Physics 78, 073901 (2015).
  * Mao and Lubensky (2018) X. Mao and T. C. Lubensky, Annual Review of Condensed Matter Physics 9, 413 (2018).
  * Paulose  _et al._ (2015) J. Paulose, A. S. Meeussen, and V. Vitelli, Proceedings of the National Academy of Sciences 112, 7639 (2015).
  * Zhang and Mao (2018) L. Zhang and X. Mao, New Journal of Physics 20, 063034 (2018).
  * Sussman  _et al._ (2016) D. M. Sussman, C. P. Goodrich, and A. J. Liu, Soft matter 12, 3982 (2016).
  * Wyart  _et al._ (2008) M. Wyart, H. Liang, A. Kabla, and L. Mahadevan, Physical review letters 101, 215501 (2008).
  * Chen  _et al._ (2011) K. Chen, M. L. Manning, P. J. Yunker, W. G. Ellenbroek, Z. Zhang, A. J. Liu, and A. G. Yodh, Phys. Rev. Lett. 107, 108301 (2011).
  * Lerner  _et al._ (2012) E. Lerner, G. Düring, and M. Wyart, Proceedings of the National Academy of Sciences 109, 4798 (2012).
  * Mosayebi  _et al._ (2014) M. Mosayebi, P. Ilg, A. Widmer-Cooper, and E. Del Gado, Phys. Rev. Lett. 112, 105503 (2014).
  * Stanifer  _et al._ (2018) E. Stanifer, P. K. Morse, A. A. Middleton, and M. L. Manning, Phys. Rev. E 98, 042908 (2018).
  * Born  _et al._ (1955) M. Born, K. Huang, and M. Lax, American Journal of Physics 23, 474 (1955).
  * DiDonna and Lubensky (2005) B. DiDonna and T. Lubensky, Physical Review E 72, 066619 (2005).
  * Maxwell (1864) J. C. Maxwell, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 27, 294 (1864).
  * Tay and Whiteley (1984) T.-S. Tay and W. Whiteley, Structural Topology, 1984, núm. 9 (1984).
  * Jacobs and Thorpe (1995) D. J. Jacobs and M. F. Thorpe, Physical review letters 75, 4051 (1995).
  * Connelly and Whiteley (1992) R. Connelly and W. Whiteley, International Journal of Space Structures 7, 153 (1992).
  * Connelly and Whiteley (1996) R. Connelly and W. Whiteley, SIAM Journal on Discrete Mathematics 9, 453 (1996).
  * Vermeulen  _et al._ (2017) M. F. Vermeulen, A. Bose, C. Storm, and W. G. Ellenbroek, Physical Review E 96, 053003 (2017).
  * Damavandi  _et al._ (2021) O. K. Damavandi, V. F. Hagh, C. D. Santangelo, and M. L. Manning, arXiv preprint arXiv:2102.11310 (2021).
  * Rocklin  _et al._ (2017) D. Z. Rocklin, S. Zhou, K. Sun, and X. Mao, Nature communications 8, 1 (2017).
  * Zhou  _et al._ (2018) D. Zhou, L. Zhang, and X. Mao, Phys. Rev. Lett. 120, 068003 (2018).
  * Ma  _et al._ (2018) J. Ma, D. Zhou, K. Sun, X. Mao, and S. Gonella, Phys. Rev. Lett. 121, 094301 (2018).
  * Zhou  _et al._ (2019) D. Zhou, L. Zhang, and X. Mao, Phys. Rev. X 9, 021054 (2019).
  * Zhou  _et al._ (2020) D. Zhou, J. Ma, K. Sun, S. Gonella, and X. Mao, Phys. Rev. B 101, 104106 (2020).
  * Sun and Mao (2020) K. Sun and X. Mao, Phys. Rev. Lett. 124, 207601 (2020).
  * Lerner (2018) E. Lerner, The European Physical Journal E 41, 1 (2018).
  * Rainone  _et al._ (2020a) C. Rainone, E. Bouchbinder, and E. Lerner, The Journal of Chemical Physics 152, 194503 (2020a).
  * Bordiga  _et al._ (2021) G. Bordiga, L. Cabras, A. Piccolroaz, and D. Bigoni, Journal of the Mechanics and Physics of Solids 146, 104198 (2021).
  * Weeks  _et al._ (1971) J. D. Weeks, D. Chandler, and H. C. Andersen, The Journal of chemical physics 54, 5237 (1971).
  * Vasisht and Del Gado (2020) V. V. Vasisht and E. Del Gado, Physical Review E 102, 012603 (2020).
  * Irving and Kirkwood (1950) J. H. Irving and J. G. Kirkwood, The Journal of Chemical Physics 18, 817 (1950).
  * Plimpton (1995) S. Plimpton, Journal of computational physics 117, 1 (1995).
  * Foster and Davis (2013) L. V. Foster and T. A. Davis, ACM Transactions on Mathematical Software (TOMS) 40, 1 (2013).
  * (95) https://people.engr.tamu.edu/davis/suitesparse.html.
  * Hopkins and Skellam (1954) B. Hopkins and J. G. Skellam, Annals of Botany 18, 213 (1954).
  * Tanguy  _et al._ (2002b) A. Tanguy, J. P. Wittmer, F. Leonforte, and J.-L. Barrat, Phys. Rev. B 66, 174205 (2002b).
  * Leonforte  _et al._ (2004) F. Leonforte, A. Tanguy, J. P. Wittmer, and J.-L. Barrat, Phys. Rev. B 70, 014203 (2004).
  * Maloney and Lemaître (2006) C. E. Maloney and A. Lemaître, Phys. Rev. E 74, 016118 (2006).
  * Schoenholz  _et al._ (2014) S. S. Schoenholz, A. J. Liu, R. A. Riggleman, and J. Rottler, Phys. Rev. X 4, 031014 (2014).
  * Rainone  _et al._ (2020b) C. Rainone, E. Bouchbinder, and E. Lerner, Proceedings of the National Academy of Sciences 117, 5228 (2020b).
  * Cubuk  _et al._ (2015) E. D. Cubuk, S. S. Schoenholz, J. M. Rieser, B. D. Malone, J. Rottler, D. J. Durian, E. Kaxiras, and A. J. Liu, Phys. Rev. Lett. 114, 108001 (2015).
  * Singh  _et al._ (2020) M. Singh, M. Ozawa, and L. Berthier, Phys. Rev. Materials 4, 025603 (2020).
  * Seto  _et al._ (2013) R. Seto, R. Mari, J. F. Morris, and M. M. Denn, Physical review letters 111, 218301 (2013).
  * Wyart and Cates (2014) M. Wyart and M. E. Cates, Physical review letters 112, 098302 (2014).
  * Lin  _et al._ (2016) N. Y. Lin, C. Ness, M. E. Cates, J. Sun, and I. Cohen, Proceedings of the National Academy of Sciences 113, 10774 (2016), https://www.pnas.org/content/113/39/10774.full.pdf .
  * Guy  _et al._ (2018) B. M. Guy, J. A. Richards, D. J. M. Hodgson, E. Blanco, and W. C. K. Poon, Phys. Rev. Lett. 121, 128001 (2018).
  * Liu  _et al._ (2021b) K. Liu, J. E. Kollmer, K. E. Daniels, J. M. Schwarz, and S. Henkes, Phys. Rev. Lett. 126, 088002 (2021b).
  * Ramaswamy  _et al._ (2021) M. Ramaswamy, I. Griniasty, D. B. Liarte, A. Shetty, E. Katifori, E. Del Gado, J. P. Sethna, B. Chakraborty, and I. Cohen, arXiv preprint arXiv:2107.13338 (2021).
  * Lagarias  _et al._ (1998) J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright, SIAM Journal on optimization 9, 112 (1998).