# Optimizing properties on the critical rigidity manifold of underconstrained
central-force networks

Tyler Hain  Department of Physics and BioInspired Institute, Syracuse
University, Syracuse NY 13210  Chris Santangelo<EMAIL_ADDRESS>Department
of Physics and BioInspired Institute, Syracuse University, Syracuse NY 13210
M. Lisa Manning<EMAIL_ADDRESS>Department of Physics and BioInspired
Institute, Syracuse University, Syracuse NY 13210

###### Abstract

Our goal is to develop a design framework for multifunctional mechanical
metamaterials that can tune their rigidity while optimizing other desired
properties. Towards this goal, we first demonstrate that underconstrained
central force networks possess a critical rigidity manifold of codimension one
in the space of their physical constraints. We describe how the geometry of
this manifold generates a natural parameterization in terms of the states of
self-stress, and then use this parameterization to numerically generate
disordered network structures that are on the critical rigidity manifold and
also optimize various objective functions, such as maximizing the bulk
stiffness under dilation, or minimizing length variance to find networks that
can be self-assembled from equal-length parts. This framework can be used to
design mechanical metamaterials that can tune their rigidity and also exhibit
other desired properties.

##  I Introduction

Mechanical metamaterials are designed with specific meso-scale structures that
imbue them with macroscopic mechanical properties or behaviors distinct from
those of the material from which they are constructed goh_nanolatticed_2021 ;
wei_multi-bionic_2022 . The design space for such materials is vast, yet
designed materials cannot match the breadth and versatility of evolved
biological materials, which exhibit an incredible ability to adapt their
mechanical properties in response to external stimulus. For example, confluent
tissues in developing embryos undergo controlled fluidization to facilitate
the large-scale shape changes associated with development
mongera_fluid–solid_2018 ; bi_density-independent_2015 ; bi_motility-
driven_2016 . On the other hand, biopolymer networks of crosslinked fibers can
be very soft for small amounts of strain, but suddenly increase their
stiffness by orders of magnitude at a particular critical strain
mackintosh_elasticity_1995 ; storm_nonlinear_2005 ; sharma_strain-
controlled_2016 ; rens_micromechanical_2018 . While many metamaterials are
designed as regular lattices of modular unit cells, biological materials are
often disordered, and consequently have the potential for a wider space of
functionality as well as increased mechanical robustness in comparison to
ordered structures, for example by avoiding shear banding and fracture
zaiser_disordered_2023 . Similar adaptability and functionality in a
metamaterial could be achieved by designing structures poised near the
rigidity transition. Yet the question of how to design disordered mechanical
metamaterials that take advantage of rigidity transitions, so that small
changes to the geometry of the structure significantly change the macroscopic
mechanical properties like biological examples, remains completely open.

In this paper, we introduce a method to completely enumerate and explore the
entire space of _geometric_ rigidity transitions of large networks. Past work
has shown that underconstrained tissue and biopolymers can be driven to a
critical configuration by applying sufficient strain merkel_minimal-
length_2019 , but it has not been clear exactly how common these critical
configurations are or whether they exhibit any emergent order or symmetry,
especially in large disordered systems. In addition, the critical
configuration that results from straining to the critical point is an emergent
feature of the energy minimization, so that it has not been possible to design
a specific desired critical configuration. One of our main results is that the
space of critical configurations is a smooth manifold almost everywhere and
consequently we can explore disordered networks beyond those that are poised
at an unprogrammed strain-induced critical point. Specifically, we derive an
algorithm that finds critical networks with other optimized desired
properties, such as maximizing the jump in bulk modulus or being composed of
fibers of the same length. These results shed new light on an under-explored
problem in metamaterials: which material properties can be designed and to
what degree? And in contrast to much previous work focused on ordered
structures, such as miura-ori patterns in origami silverberg_applied_2014 ,
our method explicitly encompasses a much broader set of disordered networks.

Our result takes advantage of the distinction between “first-order” and
“second-order” rigidity. Indeed, models for both types of biological materials
in the previous paragraph – vertex models for confluent tissues and spring
networks for biopolymer networks – can undergo rigidity transitions without
any change to their network topology. This is in stark contrast to the jamming
transition in granular systems, where the contact network topology changes as
the density increases, and the material rigidifies when the number of contacts
equals the number of degrees of freedom. In this case, first-order
perturbations to the lengths of the contacts constrain all possible
infinitesimal motions of the system, and so the system is said to be “first
order rigid”. But in the examples of tissues and biopolymer networks, the
number of constraints is often less that the number of degrees of freedom –
i.e. the systems are _underconstrained_. In those cases, rigidity occurs at
the onset of geometric incompatibility merkel_minimal-length_2019 . These
systems are called “second order rigid” because they possess an extensive
number of modes that cost zero energy with respect to first-order
perturbations to constraints, but do cost energy at second order. Although
enumerating states near a first-order rigid jamming transition is notorious
difficult due to the larger number of nearly “kissing contacts”
donev_underconstrained_2007  and the related Gardner transition
berthier_gardner_2019 , second-order rigidity occurs even with a fixed network
topology and, consequently, depends only on its geometry.

As we show in more detail, the space of critical configurations is naturally
parameterized using the states of self-stress. This parameterization is
related to the force density method, which has been used for form-finding of
pre-tensioned nets schek_force_1974 ; malerba_extended_2012  and tensegrities
tran_advanced_2010 ; wang_form-finding_2021  in regular lattices. The
parameterization enables the ability to optimize a broad array of objective
functions that can be written in terms of the states of self stress and
therefore minimized via existing gradient descent algorithms. We demonstrate
this by numerically generating and investigating optimal configurations for
several objective functions. To find networks with enhanced elastic response
at the transition, we maximize either the bulk modulus or shear modulus, and
to obtain networks with tailored structure we minimize the fluctuations in the
edges lengths or the edge tensions. We also compare these optimized
configurations to networks taken from random samples of two different
probability densities over the critical manifold. This provides insight into
design rules for multifunctional mechanical metamaterials that could
eventually guide manufacturing or self-assembly.

##  II Mathematical Preliminaries

An important aspect of any mechanical system is its rigidity: are there
deformations that cost zero energy? Traditionally this question is described
in terms of the degrees of freedom of a system and a set of constraints
$f_{\alpha}$, each of which arise from a quadratic term in the energy

$\displaystyle E=\frac{1}{2}\sum_{\alpha}f_{\alpha}^{2}.$ (1)

For example, consider a network of $N_{b}$ Hookean springs connecting $N_{v}$
vertices with the energy functional
$E=\sum_{\alpha}(L_{\alpha}-\ell_{\alpha})^{2}/2$ where $L_{\alpha}$ and
$\ell_{\alpha}$ are the actual and prescribed lengths of edge
$\alpha=1,...,N_{b}$. The degrees of freedom are the positions of the vertices
$x_{i\mu}$, where $i=1,...,N_{v}$ indexes the vertices and $\mu=1,..,d$
indexes the spatial dimensions. Each edge contributes one constraint,
$f_{\alpha}=L_{\alpha}-\ell_{\alpha}$.

In some cases, one can simply linearize the system and consider whether any
infinitesimal motions are permitted by the constraints. This information is
contained in the rigidity matrix $R_{\alpha i\mu}=\partial f_{\alpha}/\partial
x_{i\mu}$, which relates a small perturbation to the degrees of freedom
$\delta x_{i\mu}$ with the resulting first-order change to the constraints:

$\displaystyle\delta f_{\alpha}=\sum_{i\mu}\frac{\partial f_{\alpha}}{\partial
x_{i\mu}}\delta x_{i\mu}=\sum_{i\mu}R_{\alpha i\mu}\delta x_{i\mu}=0.$ (2)

The rigidity matrix specialized to Hookean springs, $\tilde{R}$, can be
written as

$\displaystyle\tilde{R}_{\alpha i\mu}=\frac{\partial f_{\alpha}}{\partial
x_{i\mu}}=g_{\alpha i}\frac{L_{\alpha\mu}}{L_{\alpha}},$ (3)

where $g_{\alpha i}$ is the signed incidence matrix of an underlying directed
graph representing the network that describes which vertices are connected by
each edge, and $L_{\alpha\mu}$ is the $\mu$ component of the vector along edge
$\alpha$. These bond vectors can in turn be written in terms of the vertex
coordinates as

$\displaystyle L_{\alpha\mu}=\sum_{i}g_{\alpha i}x_{i\mu}+b_{\alpha\mu},$ (4)

where $b_{\alpha\mu}$ is a constant matrix that encodes the boundary
conditions of the network, which can come from being in a periodic box or from
having a subset of vertices pinned in place. Additional details on how to
derive $b_{\alpha\mu}$ can be found in the Supplemental Materials.

Vectors in the right nullspace of $R$ are called linear zero modes (LZM’s) as
they result in no first-order changes to the constraints. $R$ also relates the
tensions on the edges of the network to the resulting forces on the vertices
$F_{i\mu}$:

$\displaystyle-F_{i\mu}=\frac{\partial E}{\partial
x_{i\mu}}=\sum_{\alpha}\frac{\partial E}{\partial f_{\alpha}}\frac{\partial
f_{\alpha}}{\partial x_{i\mu}}=\sum_{\alpha}f_{\alpha}{R}_{\alpha i\mu}.$ (5)

If the constraints are not all satisfied $(f_{\alpha}\neq 0)$ we say the
network is prestressed, and Eq. (5) implies that $f_{\alpha}$ must be in the
left nullspace of $R$ in order for the network to be in equilibrium. For a
given configuration, any vector $\bm{\sigma}$ in the left nullspace of $R$

$\displaystyle\sum_{\alpha}\sigma_{\alpha}R_{\alpha i\mu},$ (6)

which is called a state of self-stress, corresponds with a possible set of
tensions on the edges of the network that result in force balance on the
vertices. Note that the self-stresses are purely geometric objects that merely
describe the possible ways that a given configuration can support internal
forces.

The rank-nullity theorem applied to the rigidity matrix results in Maxwell-
Calladine constraint counting, which relates the number of independent LZM’s
($N_{0}$) and self-stresses ($N_{s}$) to the the number of degrees of freedom
($dN_{v}$) and constraints ($N_{b}$):

$\displaystyle N_{0}-N_{s}=dN_{v}-N_{b}.$ (7)

If a system has no linear zero modes (excluding trivial rigid-body
translations and rotations), then it is called first-order or infinitesimally
rigid. Systems that undergo a jamming transition are often first-order rigid:
as the free volume decreases, the number of contact between particles (and
hence the number of constraints) increases until all LZM’s are eliminated.

However, underconstrained systems can also become rigid despite having an
extensize number of LZM’s, when nonlinear effects cause the LZM’s to change
the constraints at higher order. For example, sub-isostatic spring/fiber
networks and vertex/Voronoi models for confluent tissues undergo a
strain/shape-controlled rigidity transition despite being underconstrained
bi_density-independent_2015 ; merkel_geometrically_2018 ; merkel_minimal-
length_2019 , and packings of non-spherical particles become jammed below the
isostatic point donev_underconstrained_2007 ; hecke_jamming_2009 . In these
systems the LZM’s change the constraints at quadratic order, so they are said
to be second-order rigid.

Underconstrained systems in a generic configuration do not have any self-
stresses. In this case, any linear zero mode can always be extended into a
non-linear floppy mode connelly2022frameworks . However, in certain
configurations the system has at least one self-stress. We call these special
states critical configurations, and the set of all critical configurations we
call the critical manifold. The existence of a self-stress leads to an
additional necessary (but not sufficient) condition that a LZM
$\dot{x}_{i\mu}$ must satisfy in order to be extended into a non-linear mode
connelly_rigidity_1980 ; connelly_second-order_1996 ; damavandi_energetic_2022
:

$\displaystyle\sum_{\alpha
ij}\sigma_{\alpha}^{I}\frac{\partial^{2}f_{\alpha}}{\partial x_{i\mu}\partial
x_{j\nu}}\dot{x}_{i\mu}\dot{x}_{j\nu}=\dot{\bm{x}}^{T}P^{I}\dot{\bm{x}}=0,$
(8)

for all self-stresses $\sigma^{I}$. If no LZM’s satisfy these conditions, the
system is second-order rigid. In practice, many second-order rigid systems
satisfy a weaker condition called prestress stability, where there is a
particular self-stress $\sigma_{\alpha}$ for which the associated prestress
matrix $P_{i\mu
j\nu}=\sum_{\alpha}\sigma_{\alpha}(\partial^{2}f_{\alpha}/\partial
x_{i\mu}\partial x_{j\nu})$ is positive-definite on the subspace of linear
zero modes, which clearly implies second-order rigidity. In this work, we
study networks that have at most one independent self-stress, for which
second-order rigidity and prestress stability are equivalent connelly_second-
order_1996 .

Underconstrained networks are generically floppy because the lengths of edges
for which critical configurations can occur are a set of measure zero.
Consider the three-bar linkage in Fig. 1, a simple example of an
underconstrained system. In a generic configuration, the linkage has no self-
stress, so the LZM corresponding with a shearing motion can always be extended
into a non-linear zero mode. However, if all vertices are colinear, as shown
in Fig 1A-C, a state of self-stress appears which can rigidify the linkage or
induce a branch point. While these critical configurations are only a measure-
zero subset of configuration space, the system can be driven to these special
points due to geometric incompatibility. For example, if the distance between
the ends of the linkage is larger than the sum of the preferred lengths of the
three edges, the constraints cannot all be satisfied. This leads to inevitable
prestress, and energy minimization then drives the system to a critical
configuration. This mechanism of geometric incompatibility in the preferred
geometry has been shown to be the cause of rigidity in general
underconstrained spring networks as well as in 2- and 3-dimensional vertex and
Voronoi models merkel_minimal-length_2019 . For the three bar linkage, we can
visualize the rigid portion of the critical manifold as the boundary of
geometric incompatibility. In Fig. 1F we show the full critical manifold of
the three-bar linkage in the space of the squared lengths of the three edges
berry_controlling_2022 . The grey surfaces correspond to branch points, and
the orange and purple surfaces correspond with rigid configurations. The sub-
surface of rigid configurations divides the space into two separate regions,
which are shown in Fig. 1G. Sets of squared edge lengths that are longer than
those on the colored manifold, corresponding to points in the upper right
corner of Fig 1G, are compatible with the boundary conditions, so choosing the
rest lengths from this region will result in a floppy configuration like 1A.
However, choices of squared edge lengths interior to the colored manifold in
the lower left-hand corner of Fig 1G are impossible to achieve, and the
resulting prestress then drives the system to a rigid critical configuration.

Figure 1: The critical manifold of the three bar linkage. A,B) Examples of two
different rigid critical configurations of the three bar linkage. C) Example
of a critical configuration that is a non-rigid branch point. D) All points on
the critical manifold have a corresponding geometric stress. The sphere
represents the space of these geometric stresses for the 3-bar linkage: the
orange and purple regions contain the rigid configurations, which look like A
and B respectively; the grey region contains the branch points like C; the red
line contains choices of geometric stress for which the parameterization is
not defined. As rescaling the geometric stress by the same factor results in
the same real-space configuration (Eq. 17), which is why the space of
geometric stresses is an $N_{b}-1$ -dimensional space, here represented as a
2-sphere. E) A stereographic projection of the space of geometric stresses as
described in Eq. (22). F) The critical manifold embedded in the space of
squared lengths of the edges of the three bar linkage, where the different
colored regions are the images of the corresponding regions of geometric
stresses. G) Same as panel F, except just the rigid parts of the critical
manifold are highlighted for clarity.

##  III Analytical Results

###  III.1 Geometric Stress Parameterization

When reaching the second-order rigidity transition by straining a network to
the point of geometric incompatibility, the resulting self-stress is an
emergent feature of the energy minimization, and there is no clear way to
predict the properties of the resulting rigid network. In the rest of this
work, we describe a method that characterizes the entire critical manifold of
any network with a given graph structure and boundary conditions using a
natural parameterization that can be used to rationally explore the critical
manifold and obtain configurations with specific properties.

Let us consider a network with an arbitrary energy functional
$E(x_{i\mu},p_{n})$ that depends on the node coordinates and any other control
parameters $p_{n}$, such as rest lengths or stiffnesses. Any equilibrium
configuration of this network must satisfy force balance:

$\displaystyle\frac{\partial E}{\partial x_{i\mu}}=\sum_{\alpha}\frac{\partial
E}{\partial h_{\alpha}}\frac{\partial h_{\alpha}}{\partial x_{i\mu}}=0,$ (9)

where $h_{\alpha}(x_{i\mu})$ can be any set of functions of the coordinates
that can be used to completely determine the energy:
$E(h_{\alpha}(x_{i\mu}))$. Here $h_{\alpha}$ plays a similar role to the
constraint $f_{\alpha}$ in the previous section; however, while $f_{\alpha}$
is completely determined by the specific energy function being used,
$h_{\alpha}$ can be chosen to be a more general coarse-graining of the node
coordinates into higher-order geometric quantities that determine the energy.
This allows us to write Eq. (9) separately in terms of the conjugate variables
$\partial E/\partial h_{\alpha}$, which are the generalized tensions
associated with the $h_{\alpha}$’s and contain all information about $E$, and
$\partial h_{\alpha}/\partial x_{i\mu}$, which plays the role of the rigidity
matrix by connecting the higher-level geometry of the $h_{\alpha}$’s to the
lower-level geometry of the vertex coordinates.We term this generalized
rigidity matrix $R^{h}_{\alpha i\mu}=\partial h_{\alpha}/\partial x_{i\mu}$
the “geometric response matrix”.

We can then find critical configurations by looking for states that satisfy
the force balance of Eq. (9) despite having non-zero internal stresses
$\partial E/\partial h_{\alpha}\neq 0$. To do this, we will use the
generalized tensions as a different set degrees of freedom, defining
$\sigma^{h}_{\alpha}=\partial E/\partial h_{\alpha}$. We call
$\sigma^{h}_{\alpha}$ the “generalized states of self-stress” corresponding to
the geometric response matrix for a given $h$. Then the force balance
condition $\sum_{\alpha}\sigma^{h}_{\alpha}(\partial h_{\alpha}/\partial
x_{i\mu})=0$ can be thought of as implicitly defining the critical
configuration $x_{i\mu}(\sigma^{h}_{\alpha})$ corresponding to each choice of
generalized tensions. For brevity, in remainder of this paper we sometimes
refer to these as simply ”states of self-stress”, though they are distinct
from the standard $\sigma$ defined in Eq. 6 when $h_{\alpha}\neq f_{\alpha}$.

However, most choices of $h_{\alpha}$ do not result in conjugate variables
that generate a proper set of degrees of freedom. To see this, let us now
assume that we have a central force network whose energy functional depends
only on the edge lengths: $E(L_{\alpha},p_{n})$. One might think to use
$h_{\alpha}=L_{\alpha}$. In this case, the conjugate variables $\partial
E/\partial L_{\alpha}=\tau_{\alpha}$ are just the actual tensions on the
edges, and the force balance condition is

$\displaystyle\frac{\partial E}{\partial x_{i\mu}}=\sum_{\alpha}\frac{\partial
E}{\partial L_{\alpha}}\frac{\partial L_{\alpha}}{\partial
x_{i\mu}}=\sum_{\alpha}\tau_{\alpha}\left(g_{\alpha
i}\frac{L_{\alpha\mu}}{L_{\alpha}}\right)=0.$ (10)

For a fixed $\tau_{\alpha}$, Eq. (10) cannot be solved to get a corresponding
configuration $x_{i\mu}$. Consider the three-bar linkage in Fig. 1A; any such
critical configuration with all edges under tension must have equal tension on
all edges to maintain force balance: $\tau_{1}=\tau_{2}=\tau_{3}$. Therefore,
if we choose $\tau_{\alpha}\propto[1\,1\,1]^{T}$ Eq. (10) will have infinitely
many solutions, and for any other choice of $\tau_{\alpha}$ Eq. (10) will have
no solutions. Therefore we reject this choice of $h_{\alpha}$ and its
conjugate self-stress $\tau_{\alpha}$.

Instead, to ensure we can solve the force balance equation, we should choose
$h_{\alpha}$ to be the simplest function of $L_{\alpha}$ that is quadratic in
$x_{i\mu}$ so that taking a derivative results in an equation that is linear
in the coordinates. We can do this using $h_{\alpha}=L_{\alpha}^{2}/2$, which
results in the conjugate variable

$\displaystyle\sigma^{h}_{\alpha}=\frac{\partial
E}{\partial(L_{\alpha}^{2}/2)}=\frac{\partial E}{\partial
L_{\alpha}}\frac{1}{L_{\alpha}}=\frac{\tau_{\alpha}}{L_{\alpha}},$ (11)

which corresponds with the force density (tension per length) of each edge. We
use this choice of $h_{\alpha}$ and the corresponding generalized self-stress
for the remainder of this work. We will refer to this particular self-stress
as the geometric stress, as we demonstrate below that it naturally encodes the
geometry of the critical manifold. The force balance condition then becomes

$\displaystyle\frac{\partial E}{\partial x_{i\mu}}$
$\displaystyle=\sum_{\alpha}\frac{\partial
E}{\partial(L_{\alpha}^{2}/2)}\frac{\partial(L_{\alpha}^{2}/2)}{\partial
x_{i\mu}}$ $\displaystyle=\sum_{\alpha}\sigma^{h}_{\alpha}\left(g_{\alpha
i}L_{\alpha\mu}\right)=0,$ (12)

which eliminates the factor of $1/L_{\alpha}$ in Eq. (10). Now, for any choice
of $\sigma^{h}_{\alpha}$, we can write Eq. (12) as a linear equation in the
coordinates:

$\displaystyle\frac{\partial E}{\partial x_{i\mu}}$
$\displaystyle=\sum_{\alpha}\sigma^{h}_{\alpha}g_{\alpha
i}L_{\alpha\mu}=\sum_{\alpha j}\sigma^{h}_{\alpha}g_{\alpha i}(g_{\alpha
j}x_{j\mu}+b_{\alpha\mu})$ (13)
$\displaystyle=\sum_{j\nu}\left(\sum_{\alpha}\sigma^{h}_{\alpha}g_{\alpha
i}g_{\alpha
j}\delta_{\mu\nu}\right)x_{j\nu}+\sum_{\alpha}\sigma^{h}_{\alpha}g_{\alpha
i}b_{\alpha\mu}$ (14) $\displaystyle=\sum_{j\nu}P_{i\mu
j\nu}\,x_{j\nu}+B_{i\mu}=0,$ (15)

where $P$ plays the role of the prestress matrix, and is given by

$\displaystyle P_{i\mu
j\nu}=\sum_{\alpha}\sigma^{h}_{\alpha}\frac{\partial^{2}h_{\alpha}}{\partial
x_{i\mu}\partial x_{j\nu}}=\sum_{\alpha}\sigma^{h}_{\alpha}g_{\alpha
i}g_{\alpha j}\delta_{\mu\nu}.$ (16)

Notice that this geometric prestress matrix has no explicit dependence on the
coordinates $x_{i\mu}$, and depends only on the geometric stress
$\bm{\sigma}^{h}$ and the structure of the underlying graph encoded by $g$.
$P$ is also block-diagonal with the respect to the spatial directions
($P_{i\mu j\nu}=P_{ij}\delta_{\mu\nu}$), and the non-zero blocks are exactly
the weighted Laplacian matrix associated with with the underlying graph with
the geometric stress as the weights. We will sometimes refer to just this non-
zero block as the prestress matrix, in which case it will have only two
indices: $P_{ij}=\sum_{\alpha}\sigma^{h}_{\alpha}g_{\alpha i}g_{\alpha j}$.

This matrix always has the trivial zero mode $[11\dots 1]^{T}$ corresponding
to rigid-body translations. In this work we will choose the components
$\sigma^{h}_{\alpha}$ to be all positive, corresponding to networks with all
edges under tension. In this case, $P$ has no other zero modes, so Eq. (15)
can always be solved to obtain a unique critical configuration up to
translations:

$\displaystyle
x_{i\mu}(\bm{\sigma}^{h})=-\sum_{j}\left(P^{-1}\right)_{ij}\left(\sum_{\alpha}\sigma^{h}_{\alpha}g_{\alpha
j}b_{\alpha\mu}\right),$ (17)

where ${P}^{-1}$ should be understood as the pseudo-inverse of $P$. Note that
the translational zero mode $[11...1]^{T}$ is also in the nullspace of the
incidence matrix $g_{\alpha j}$, which guarantees that Eq. (15) has a
solution.

To describe each critical configuration uniquely, we can instead write the
parameterization directly for the vectors along the edges:

$\displaystyle
L_{\alpha\mu}(\bm{\sigma}^{h})=\sum_{\beta}\left\\{\delta_{\alpha\beta}-\left[\sum_{ij}g_{\alpha
j}\left({P}^{-1}\right)_{ij}g_{\beta
j}\right]\sigma^{h}_{\beta}\right\\}b_{\beta\mu}.$ (18)

The mapping in Eq. (17) provides a natural parameterization of the critical
manifold: while a generic choice of vertex coordinates will not result in a
critical configuration, almost any choice of $\sigma^{h}_{\alpha}$ will
produce a unique critical configuration (up to rigid-body motions). This is
essentially the force-density method, which has been used by engineers to find
equilibrium shapes for tensioned nets and tensegrity structures that are
stabilized by internal stresses schek_force_1974 ; wang_form-finding_2021 . Of
course, for the special case where the energy functional represents springs
with no rest length, $f_{\alpha}=L^{2}_{\alpha}/2$, the geometric response
matrix $\partial h_{\alpha}/\partial x_{i\mu}$ is the standard rigidity
matrix.

Note that this geometric stress parameterization can still be used for
overconstrained networks, although it is significantly less useful: since
constraint counting implies that every configuration possesses a self-stress,
the critical manifold fills all of configuration space, negating the need for
an alternate parameterization. In addition, for networks with multiple self-
stresses, the mapping is no longer one-to-one, as multiple independent
geometric stresses get mapped onto the same critical configuration. In future
work, it may be interesting to use this method to find any relationship that
may exist between different self-stresses of the same configuration. However,
for this work we will focus only on networks with exactly one self-stress.

Once the configuration corresponding to a given $\bm{\sigma}^{h}$ has been
obtained, we still need to choose any control parameters such that force
balance is actually achieved with the desired energy functional. For example,
if we want a network with a Hookean energy $E=\sum
k_{\alpha}/2(L_{\alpha}-\ell_{\alpha})^{2}$, we can ensure it is in mechanical
equilibrium by choosing the rest lengths and stiffnesses that satisfy
$\partial
E/\partial(L_{\alpha}^{2}/2)=k_{\alpha}(L_{\alpha}-\ell_{\alpha})/L_{\alpha}=c\sigma^{h}_{\alpha}$.
Choosing $c=0$ (setting $\ell_{\alpha}=L_{\alpha}$) results in the
configuration being at the critical point, as it possesses a self-stress but
is not under any actual stress.

Different choices of the energy functional result in different effective
stiffness of the resulting configuration, which (following
guest_stiffness_2006 ) we can see by writing the Hessian matrix of a generic
$E(L_{\alpha})$:

$\displaystyle H_{i\mu j\nu}$ $\displaystyle=\frac{\partial^{2}E}{\partial
x_{i\mu}\partial x_{j\nu}}{}$ $\displaystyle=\sum_{\alpha\beta}\frac{\partial
h_{\alpha}}{\partial x_{i\mu}}\frac{\partial^{2}E}{\partial h_{\alpha}\partial
h_{\beta}}\frac{\partial h_{\beta}}{\partial x_{j\nu}}+\frac{\partial
E}{\partial h_{\alpha}}\frac{\partial^{2}h_{\alpha}}{\partial x_{i\mu}\partial
x_{j\nu}}{}$ $\displaystyle={R^{h}}^{T}KR^{h}+P.$ (19)

$K$ is a matrix of “modified stiffnesses”, which contains all information
about the specific energy function: $R^{h}$ depends only on the geometry of
the current configuration, and $P$ depends only on the current geometric
stress. For a Hookean spring network, $K$ is a diagonal matrix with entries

$\displaystyle K_{\alpha\alpha}=\frac{\partial^{2}E}{\partial
h_{\alpha}^{2}}=\frac{1}{L_{\alpha}^{2}}\left(\frac{\partial^{2}E}{\partial
L_{\alpha}^{2}}-\frac{1}{L_{\alpha}}\frac{\partial E}{\partial
L_{\alpha}}\right)=\frac{k_{\alpha}-\sigma^{h}_{\alpha}}{L_{\alpha}^{2}},$
(20)

where $k_{\alpha}=\partial^{2}E/\partial L_{\alpha}^{2}$ is the actual
stiffness of edge $\alpha$, which is constant for Hookean springs but could be
a function of the edge length for a nonlinear spring potential. Other models
can have non-zero entries off the main diagonal of $K$. For example, a vertex
model without an area term would have $K_{\alpha\beta}\neq 0$ whenever edges
$\alpha$ and $\beta$ are part of the same cell.

From Eq. (19), we can see that the for small amounts of prestress, the low-
frequency spectrum of the Hessian is controlled solely by the prestress
matrix, and is not affected by the specific energy functional. The linear zero
modes, which are determined only by the geometry of the current configuration,
are zero modes of $R^{h}$, and so are also zero modes of the first term in the
Hessian. Therefore, their stiffnesses depend only on $P$, which in turn only
depends on the geometric stress and the graph structure of the network. In the
Supplemental Materials, we compare these expressions for the Gram and
prestress terms of the Hessian to those derived from expanding the derivatives
in terms of the constraints $f_{\alpha}=L_{\alpha}-\ell_{\alpha}$ rather than
$h_{\alpha}=L_{\alpha}^{2}/2$.

Figure 2: Geometry of the critical manifold of the three bar linkage.
Green/purple lines correspond with evenly spaced grid lines in the space of
the Cartesian coordinates $s_{1},s_{2}$ in Fig. 1E. For each point on the
manifold, the geometric stress $\sigma^{h}$ (in orange) is always normal to
the surface.

###  III.2 Geometry of the Critical Manifold

Given that we have an analytic expression for the geometric stress
parameterization, we can describe the geometry of the critical manifold, which
we consider as a subset of the space of squared edge lengths,

$\displaystyle
h_{\alpha}(\bm{\sigma}^{h})=\frac{1}{2}L^{2}_{\alpha}(\bm{\sigma}^{h})=\frac{1}{2}\sum_{\mu}L^{2}_{\alpha\mu}(\bm{\sigma}^{h}).$
(21)

This surface is a smooth manifold almost everywhere, except for choices of
$\bm{\sigma}^{h}$ for which the prestress matrix has non-trivial zero modes
and thus Eq. (15) has no unique solution (up to translations). As this space
is parameterized by the $N_{b}$ components of the geometric stress, one might
expect the critical manifold to be $N_{b}$-dimensional. However, the
parameterization in Eq. (18) is invariant under rescaling the geometric
stress: $L_{\alpha\mu}(c\bm{\sigma}^{h})=L_{\alpha\mu}(\bm{\sigma}^{h})$. This
implies that we should think of the space of geometric stresses as the
$(N_{b}-1)$-dimensional projective space $\mathcal{RP}^{N_{b}-1}$, for which
the $N_{b}$ components $\sigma^{h}_{\alpha}$ form a set of homogeneous
coordinates. In Fig 1D we represent this manifold for the 3-bar linkage as a
2-sphere. We can define a unique coordinate for each point on the manifold by
taking a stereographic projection of the geometric stresses. We can choose a
given unit-normalized geometric stress $\hat{\bm{\sigma}}^{0}$ as the “north
pole” of the projection and choose an orthonormal basis
$\\{\hat{\bm{\sigma}}^{0},\bm{v}^{1},\dots,\bm{v}^{N_{b}-1}\\}$ for
$\mathcal{R}^{N_{b}}$. Then we can specify each geometric stress by $N_{b}-1$
Cartesian coordinates $s_{\gamma}$:

$\displaystyle\sigma^{h}_{\alpha}(s_{\gamma})=\frac{\sum_{\gamma}v_{\alpha}^{\gamma}s_{\gamma}+(1-|\bm{s}|^{2}/4)\hat{\sigma}_{\alpha}^{0}}{1+|\bm{s}|^{2}/4}.$
(22)

In Figure 1E we show this projection for the three bar linkage, using
$\hat{\bm{\sigma}}^{0}=(1/\sqrt{3})[1\,1\,1]^{T}$. When writing geometric
quantities we will use the derivative of this projection

$\displaystyle V_{\alpha\gamma}:=\frac{\partial\sigma^{h}_{\alpha}}{\partial
s_{\gamma}}=\frac{v_{\alpha}^{\gamma}-s_{\gamma}[\hat{\sigma}_{\alpha}^{0}+\sigma^{h}_{\alpha}(\bm{s})]/2}{1+|\bm{s}|^{2}/4}.$
(23)

While it is sometimes easier to use the projected coordinates $s_{\gamma}$ to
write some geometric quantities (as we do below), we will typically just use
the redundant homogeneous coordinates $\sigma^{h}_{\alpha}$.

Interestingly, in addition to acting as the coordinates for the critical
manifold, the geometric stress $\bm{\sigma}^{h}$ is also normal to the surface
berry_controlling_2022 . Consider any smooth trajectory
$\sigma^{h}_{\alpha}(t)$ through the space of geometric stresses, and its
image in the space of lengths $h_{\alpha}(\sigma^{h}(t))$. The tangent vectors
along this curve are given by

$\displaystyle\frac{\partial h_{\alpha}}{\partial
t}=\sum_{i\mu\beta}\frac{\partial h_{\alpha}}{\partial x_{i\mu}}\frac{\partial
x_{i\mu}}{\partial\sigma^{h}_{\beta}}\frac{\partial\sigma^{h}_{\beta}}{\partial
t}.$ (24)

Then because the geometric stresses are defined as being in the nullspace of
the geometric response matrix,

$\displaystyle\sum_{\alpha}\sigma^{h}_{\alpha}\frac{\partial
h_{\alpha}}{\partial x_{i\mu}}\Bigr{|}_{\sigma^{h}_{\alpha}}=0,$ (25)

the geometric stresses are orthogonal to the tangent vectors at each point on
the critical manifold

$\displaystyle\sum_{\alpha}\sigma^{h}_{\alpha}\frac{\partial
h_{\alpha}}{\partial t}\Bigr{|}_{\sigma^{h}_{\alpha}}=0,$ (26)

so the geometric stresses are normal to the surface.

We want to be able to search the critical manifold for particular
configurations. To do this, we need to know how a small step in geometric
stress space changes the resulting critical configuration: $\delta
L_{\alpha\mu}\sim(\partial
L_{\alpha\mu}/\partial\sigma^{h}_{\beta})\delta\sigma^{h}_{\beta}$. We show in
the Supplemental Materials that this derivative of the parameterization with
respect to the geometric stress is given by

$\displaystyle\frac{\partial
L_{\alpha\mu}}{\partial\sigma^{h}_{\beta}}=-\left[\sum_{ij}g_{\alpha
i}\left(P^{-1}\right)_{ij}g_{\beta j}\right]L_{\beta\mu},$ (27)

which we can then use to write the derivatives of the squared lengths:

$\displaystyle\frac{\partial h_{\alpha}}{\partial\sigma^{h}_{\beta}}$
$\displaystyle=\sum_{\mu}L_{\alpha\mu}\frac{\partial
L_{\alpha\mu}}{\partial\sigma^{h}_{\beta}}=-\sum_{ij\mu}L_{\alpha\mu}g_{\alpha
i}\left(P^{-1}\right)_{ij}g_{\beta j}L_{\beta\mu}$
$\displaystyle=-\left(R^{h}P^{-1}{R^{h}}^{T}\right)_{\alpha\beta}:=B_{\alpha\beta},$
(28)

This is a key equation that allows us to search the geometric stress space in
order to optimize material properties that are naturally written as a function
of $L_{\alpha\mu}$, as shown in Eq. (34) and discussed in the next section.

Notice that this matrix is symmetric, because it is closely related to the
bond operator which has been used to quantify the elastic response of
amorphous particle packings giannini_bond-space_2021 ; lerner_breakdown_2014 :

$\displaystyle\Xi_{\alpha\beta}=\left(RH^{-1}R^{T}\right)_{\alpha\beta},$ (29)

which gives the change in length of bond $\alpha$ due to a force dipole on
bond $\beta$. Here $H$ is the full Hessian matrix and $R=\partial
L_{\alpha}/\partial x_{i}$ is the usual rigidity matrix. Since Eq. (28) is a
function of only the geometric prestress matrix rather than the full Hessian,
it does not account for the full linear response of the network, but just
describes how the structure must change to accommodate an external force
dipole _before_ relaxation.

The matrix $B$ in Eq. (28) also gives the tangent space at each point on the
critical manifold. The columns of Eq. (28) span the tangent space, but there
is one redundant column because we are using the homogeneous coordinates
$\sigma^{h}_{\alpha}$. We can get rid of this by using a projection as in Eq.
(22) to write the $N_{b}-1$ tangent vectors:

$\displaystyle T_{\alpha}^{\gamma}=\frac{\partial h_{\alpha}}{\partial
s_{\gamma}}=\sum_{\beta}\frac{\partial
h_{\alpha}}{\partial\sigma^{h}_{\beta}}\frac{\partial\sigma^{h}_{\beta}}{\partial
s_{\gamma}}=(BV)_{\alpha\gamma}.$ (30)

We can then write the metric as the inner products of the tangent vectors:

$\displaystyle
g_{\gamma\delta}=\sum_{\alpha}T_{\alpha}^{\gamma}T_{\alpha}^{\delta}=(V^{T}B^{2}V)_{\gamma\delta}.$
(31)

The fact that the geometric stress is also the normal vector to the smooth
part of the critical manifold results in an interesting coincidence. We show
in the Supplemental Materials that the second fundamental form of the critical
manifold, which describes the local curvature information, is given by

$\displaystyle\mathrm{I\\!I}_{\gamma\delta}=\sum_{\alpha}\frac{\partial^{2}h_{\alpha}}{\partial
s_{\gamma}\partial
s_{\delta}}\frac{\sigma^{h}_{\alpha}}{|\bm{\sigma}^{h}|}=\frac{1}{|\bm{\sigma}^{h}|}(V^{T}BV)_{\gamma\delta}.$
(32)

This means essentially all geometric information about the critical manifold
is contained in the matrix $B$ from Eq. (28), along with whatever projection
is chosen from the homogeneous coordinates $\sigma^{h}_{\alpha}$ to the
Cartesian coordinates $s_{\gamma}$.

##  IV Numerical Results

Figure 3: Representatives of all groups of optimized and random critical
configurations from a single network structure. A) The configuration that
minimizes the length fluctuations $V_{L}$. B,C) Two degenerate optimal
configurations of the tension fluctuations $V_{\tau}$ that are related by an
isogonal transformation. D) The configuration that maximizes the bulk modulus
at the transition $B_{0}$ obtained by setting the geometric stress equal to
the edge stiffnesses $(\sigma^{h}_{\alpha}=k_{\alpha}=1)$. E,F) Two other
configurations with optimum $B_{0}$ that are related to the configuration in
D) by two different isogonal modes. G,H,I) Configurations that maximize the
shear modulus at the transition with various lower thresholds on the geometric
stress. Note that lowering the threshold allows the network to become more
aligned, increasing the resulting shear modulus. J,K) Two configurations
generated by randomly assigning rest lengths and volumetrically straining to
the critical point. L,M) Two configurations generated by sampling a random
positive-definite geometric stress and using the geometric stress
parameterization to obtain the corresponding configuration. As the tensions in
the networks are defined only up to an overall scale, we fix the normalization
of the geometric stress by setting the pressure equal to unity
$(P=\sum_{\alpha}\sigma^{h}_{\alpha}L_{\alpha}^{2}/2V=1)$ and coloring the
edges by the resulting tensions $\tau_{\alpha}=\sigma^{h}_{\alpha}L_{\alpha}$.

Here, we specialize the results from the previous section to two-dimensional
periodic random regular networks generated from a Voronoi tessellation, which
ensures that the networks are random regular graphs with a uniform
connectivity $z=3$ and thus possess at most one self-stress. To keep the
lengths of the edges $\mathcal{O}(1)$, we scale the size of the periodic box
with the number of vertices: $V=N_{v}$. We also limit our networks to have
positive-definite geometric stresses, equivalent to all edges being under
tension rather than compression, although other possibilities would be
interesting for future work. We assume all networks have a Hookean energy
function

$\displaystyle
E=\frac{1}{2}\sum_{\alpha}k_{\alpha}(L_{\alpha}-\ell_{\alpha})^{2},$ (33)

where the edge stiffnesses $k_{\alpha}$ are set to unity unless stated
otherwise.

Using the derivative of the geometric stress parameterization in Eq. (27), we
can calculate the total derivative of any function
$O(\sigma^{h}_{\alpha},L_{\alpha\mu})$ with respect to the geometric stress:

$\displaystyle\frac{\text{d}O}{\text{d}\sigma^{h}_{\alpha}}=\frac{\partial
O}{\partial\sigma^{h}_{\alpha}}+\sum_{\beta\mu}\frac{\partial O}{\partial
L_{\beta\mu}}\frac{\partial L_{\beta\mu}}{\partial\sigma^{h}_{\alpha}}.$ (34)

We then use the FIRE gradient descent algorithm guenole_assessment_2020  to
optimize the objective function on the critical manifold (details in the
Supplemental Materials). We consider two different types of objective
functions: structure-based functions that describe the geometry of a network
or its dual force network, and response-based functions that quantify the
network’s mechanical behavior.

###  IV.1 Generating Critical Networks

We generate three classes of networks on the second order rigidity critical
manifold.

  1. 1.

Response objective functions: networks optimized to maximize a specific
mechanical response

  2. 2.

Structural objective functions: networks that optimize a specific structural
quantity, and

  3. 3.

Random samples: networks that are constructed from random samples of a
distribution over the critical manifold.

####  IV.1.1 Response Objective Functions

Underconstrained networks on the floppy side of the critical manifold
generically have all of their elastic moduli vanish, but upon reaching the
critical manifold the moduli jump discontinuously to a finite value. In the
Supplemental Materials, we use linear response (following
maloney_amorphous_2006 ) to calculate this discontinuity in the elastic
modulus associated with a bulk affine deformation parameterized by a strain
$\gamma$:

$\displaystyle\mu=\frac{1}{V}\frac{\text{d}^{2}E}{\text{d}\gamma^{2}}$
$\displaystyle=\frac{1}{2V}\frac{\left(\sum_{\alpha}\sigma^{h}_{\alpha}\dot{h}_{\alpha}\right)^{2}}{\sum_{\alpha}(\sigma^{h}_{\alpha})^{2}h_{\alpha}/k_{\alpha}},$
(35)

where $\dot{h}_{\alpha}=\partial h_{\alpha}/\partial\gamma$ quantifies how the
critical manifold changes due to the bulk deformation to the periodic box.

Past the critical manifold, the presence of prestress leads to additional
corrections to the modulus, but right at the transition, where the prestress
vanishes, we can write the modulus as a simple function of the geometric
stress, squared lengths, and edge stiffnesses. Note that the geometric stress
$\sigma^{h}_{\alpha}$ in Eq. (35) is explicitly acting as the coordinates for
the critical manifold and not describing any actual tension on the edges
(hence why $\mu$ is independent of the normalization of $\bm{\sigma}^{h}$).

Eq. (35) has a nice interpretation in terms of the geometry of the critical
manifold. The numerator of Eq. (35) is maximized at the particular point on
the critical manifold (that is, the particular critical configuration) where
the bulk deformation causes the critical manifold to move parallel to the
geometric stress, which is normal to the critical manifold. Then this
configuration would be pushed furthest away from the critical manifold by a
given strain. However, this argument does not account for the non-affine
relaxation that occurs during a bulk deformation. If it costs a lot of energy
to strain a configuration away from the critical manifold along the normal
direction, then there will be a more non-affine deformation that decreases the
elastic modulus. The denominator of Eq. (35) quantifies this energy cost.
Since moving normal to the critical manifold has the effect of simply scaling
the tensions, we can write

$\displaystyle\frac{\partial E}{\partial
h_{\alpha}}=k_{\alpha}\frac{L_{\alpha}-\ell_{\alpha}}{L_{\alpha}}=c\sigma^{h}_{\alpha},$
(36)

where $c>0$ quantifies how far we have moved away from the manifold. Then we
can write

$\displaystyle\sum_{\alpha}(\sigma^{h}_{\alpha})^{2}h_{\alpha}/k_{\alpha}$
$\displaystyle=\frac{1}{c^{2}}\sum_{\alpha}k^{2}_{\alpha}\frac{(L_{\alpha}-\ell_{\alpha})^{2}}{L^{2}_{\alpha}}\frac{L_{\alpha}^{2}}{2k_{\alpha}}$
$\displaystyle=\frac{1}{c^{2}}\sum_{\alpha}\frac{k_{\alpha}}{2}(L_{\alpha}-\ell_{\alpha})^{2}=\frac{E(c)}{c^{2}}.$
(37)

Since the energy grows quadratically with $c$, this quantity is constant and
finite even as $c$ goes to 0. Therefore, we can think of the process of
maximizing Eq. (35) as looking for a configuration that can be moved the
largest distance away from the critical manifold for the smallest energy cost.

First, we optimize the bulk modulus at the rigidity transition, which is given
by

$\displaystyle B_{0}$
$\displaystyle=V\frac{\text{d}^{2}E}{\text{d}V^{2}}=\frac{1}{4V}\frac{\left(\sum_{\alpha}\sigma^{h}_{\alpha}L_{\alpha}^{2}\right)^{2}}{\sum_{\alpha}(\sigma^{h}_{\alpha})^{2}L_{\alpha}^{2}},$
(38)

since $\dot{h}_{\alpha}=\partial h_{\alpha}/\partial V=h_{\alpha}/V$, and
where for we have taken the stiffnesses of the edges $k_{\alpha}=1$ as the
natural unit.

As with the structural objective functions, we consider this quantity as a
function only of the geometric stress, since the geometric stress
parameterization gives us $L_{\alpha}(\bm{\sigma}^{h})$. We then calculate the
total derivative $\text{d}B_{0}/\text{d}\sigma^{h}_{\alpha}$ and use gradient
ascent to find a critical configuration that maximizes $B_{0}$.

We find that there is a degenerate subset of critical configurations that have
the same optimum value of $B_{0}$. All degenerate optima of $B_{0}$ possess
the same dual force network, and thus are related by isogonal modes
noll_active_2017 ; brauns_geometric_2024 . To demonstrate this, we calculate
the Hessian matrix of $B_{0}$ with respect to the geometric stress
$B_{0}^{\prime\prime}=\partial^{2}B_{0}/\partial\sigma^{h}_{\alpha}\partial\sigma^{h}_{\beta}$
in order to quantify the shape of the cost landscape. We find that this matrix
is negative semi-definite, with the number of flat directions equal to the
number of polygons in the network.

$\displaystyle\text{Null}(B_{0}^{\prime\prime})=N_{v}/2-1=N_{f}.$ (39)

Since any isogonal mode can be written as a sum of modes that expand or
contract the individual polygons of the network, this counting argument is
consistent with the degenerate subspace of optima possessing the same force
network. $B_{0}$ cannot be written as a function only of the tensions, so this
is a somewhat surprising symmetry that should be studied in future work.

To parameterize the degenerate subset of configurations with the same bulk
modulus $B_{0}$, we write the modulus as a sum of affine and non-affine
contributions

$\displaystyle B_{0}$
$\displaystyle=B_{0}^{\text{aff}}+(B_{0}-B_{0}^{\text{aff}}),$ (40)

where $B_{0}^{\text{aff}}=V\frac{\partial^{2}E}{\partial
V^{2}}=\sum_{\alpha}k_{\alpha}L_{\alpha}^{2}/4V$. Configurations in the
degenerate subset differ in how $B_{0}$ is split between the affine and non-
affine terms. Interestingly, the most affine configuration, with non-affine
component equal to zero, is given by

$\displaystyle 0$
$\displaystyle=B_{0}-B_{0}^{\text{aff}}=\frac{1}{4V}\frac{\left(\sum_{\alpha}\sigma^{h}_{\alpha}L_{\alpha}^{2}\right)^{2}}{\sum_{\alpha}(\sigma^{h}_{\alpha})^{2}L_{\alpha}^{2}/k_{\alpha}}-\sum_{\alpha}k_{\alpha}L_{\alpha}^{2}/4V,$
(41)

which occurs precisely when the geometric stress is proportional to the the
edge stiffness, $\sigma^{h}_{\alpha}=k_{\alpha}$. The most affine example
shown in Fig. 3D, and several other isogonal modes are shown in Fig. 3E,F. One
can think of this as an algorithm to construct the degenerate space of maximal
$B_{0}$ configurations: they are the set of all configurations with the same
force network as the state where $\sigma^{h}_{\alpha}=k_{\alpha}$.

We also attempt to maximize the shear modulus of a network at the rigidity
transition. While the shear modulus of an isotropic second-order rigid network
vanishes at the transition point in the thermodynamic limit, anisotropic
networks can have a finite shear modulus. We choose the modulus associated
with the pure shear deformation

$A_{\mu\nu}=I+\gamma\left[\begin{array}[]{cc}0&1/2\\\
1/2&0\end{array}\right],$

which is given by

$\displaystyle G_{0}$
$\displaystyle=\frac{1}{V}\frac{\text{d}^{2}E}{\text{d}\gamma^{2}}=\frac{1}{V}\frac{\left(\sum_{\alpha}\sigma^{h}_{\alpha}L_{\alpha
x}L_{\alpha
y}\right)^{2}}{\sum_{\alpha}(\sigma^{h}_{\alpha})^{2}L_{\alpha}^{2}},$ (42)

where we have again set the edge stiffnesses to $k_{\alpha}=1$.

However, this process is not as well-behaved as optimizing the bulk modulus,
which results in a rigid configuration with a positive-definite geometric
stress and edge lengths clustered around the natural length scale $\bar{L}=1$.
In contrast, when walking through the space of geometric stresses along the
direction of increasing shear modulus, a subset of edges that are aligned with
the imposed shear direction begin to grow until they become on the order of
the linear box size $\sqrt{N_{v}}$. In Fig. 4A-B we plot the shear modulus and
two structural indicators during the optimization procedure for an ensemble of
100 networks. The colored regions corresponds with standard deviations over
the ensemble, and the fluctuations are due to an adaptive time step in the
gradient descent algorithm. In Fig. 4B we plot the length of the longest edge
in the network normalized by the linear box dimension
$\text{max}(L_{\alpha})/\sqrt{N_{v}}$ in blue, and the minimum tension in the
network normalized by the largest-magnitude component of tension
$\text{min}(\tau_{\alpha})/\text{max}(|\tau_{\alpha}|)$ in red (we normalize
the tensions because, like $\bm{\sigma}^{h}$, they are defined only up to an
overall scale factor). Then if $\text{max}(L_{\alpha})/\sqrt{N_{v}}=1$ the
longest edge spans the entire system, and if
$\text{min}(\tau_{\alpha})/\text{max}(|\tau_{\alpha}|)=-1$ the edge with the
largest magnitude of tension is under compression.

Fig. 4C-E show a network at various points along the optimization, where red
edges are under tension and blue edges are under compression. The final
configurations reached after the optimization are essentially tensegrity
structures, which in two dimensions require a significant number of edges to
wrap around the periodic box and cross over each other. For this reason we do
not show these final configurations here, see Movie 1 in the Supplemental
Materials for a full animation of a typical optimization. In Fig. 4G-I we
investigate the size scaling of the subset of compressive edges. The number of
edges with negative tension scales with the linear box size
$N(\tau<0)\propto\sqrt{N_{v}}$, while the sum of lengths of negative-tension
edges scales with the system area $L(\tau<0)\propto N_{v}$. This implies that
the average length of the negative-tension edges $L(\tau<0)/N(\tau<0)$ also
scales with the linear box size.

This makes sense in the context of our interpretation of the numerator and
denominator of the elastic moduli. We saw that any modulus becomes larger when
the associated deformation of the critical manifold $\dot{\bm{h}}$ is close to
being aligned with the normal vector $\bm{\sigma}^{h}$. For a volumetric
expansion, $\dot{h}_{\alpha}=L_{\alpha}^{2}/2V$ is always positive definite,
so it makes sense that the resulting optimal geometric stresses are also
positive. However, for pure shear along the direction $45\degree$ from the
$x$-axis of the periodic box we have $\dot{h}_{\alpha}=L_{\alpha x}L_{\alpha
y}$ which necessarily has many negative components, so gradient ascent will
bring the network to a configuration with many negative geometric stress
components.

We can see that the edges begin growing to an extent that would be difficult
to realize in an experiment, with system-spanning edges and edge overlaps, at
the same point in the optimization that the smallest tensions become negative.
This suggests that to study configurations that are possible to realize in an
experiment, we might add additional constraints on the geometric stress to
keep the edges under tension rather than compression. Therefore, we add an
additional term to the objective function that penalizes components of
$\sigma^{h}_{\alpha}$ that are less than a given threshold value $r>0$,

$\displaystyle\mathcal{O}(\sigma^{h},r)=G_{0}(\sigma^{h})-\frac{1}{N_{b}}\sum_{\alpha}\left(\frac{\sigma^{h}_{\alpha}}{r}\right)^{-10},$
(43)

where the normalization of the geometric stress is fixed as
$\sum_{\alpha}(\sigma^{h}_{\alpha})^{2}=N_{b}$. An illustration of this
constraint in geometric stress space is shown in the inset to Fig. 5B, with
the allowed region highlighted in teal. In Fig. 3 we show networks that result
from maximizing Eq. (43) with various values of the threshold $r$. The open
teal circles in Fig. 5A demonstrates that this constrained optimization does
generate systems with finite $G_{0}$, independent of system size $N_{v}$, in
contrast to other critical networks (closed circles in Fig. 5A) which scale as
$1/N_{v}$. Lower values of this threshold result in the network becoming
increasingly aligned along the shear direction.

Figure 4: Unconstrained optimization of the shear modulus $G_{0}$ A) Shear
modulus as a function of effective number of optimization steps $s$ during the
optimization procedure without any additional constraints on
$\sigma^{h}_{\alpha}$, averaged over an ensemble of 100 networks. B) For the
same ensemble, the length of the longest edge of the network in units of the
linear box dimension $\sqrt{V}$ (blue) and the minimum component of the edge
tensions $\tau_{\alpha}=L_{\alpha}\sigma^{h}_{\alpha}$ normalized by the
largest-magnitude component of tension. C-E) Images of a network at different
points along the optimization. Red edges are under tension, blue edges are
under compression. See Movie 1 in the Supplemental Materials for an animation
of the optimization. (F,G) Features of the final optimized networks as a
function of system size: F) System size scaling of the number of edges under
compression $N(\tau<0)$. G) System size scaling of the sum of lengths of edges
under compression $L(\tau<0)$. Figure 5: Constrained optimization of the shear
modulus $G_{0}$. A) Size scaling of the optimal $G_{0}$ for all groups of
critical configurations. Configurations that were not optimized for shear
modulus (filled circles) have the expected finite-size scaling for isotropic
systems $G_{0}\sim 1/N$, while configurations that were optimized for $G_{0}$
(open circles) do not scale with system size. B) Optimal values of $G_{0}$ for
various thresholds $r$ and for system sizes $N_{v}=100$ (dotted), $500$
(dashed) and $1000$ (solid). As an example, the inset illustrates the region
of geometric stresses of the three bar linkage that satisfy $\sigma^{h}>r$
(colored green).

####  IV.1.2 Structural Objective Functions

Suppose we wanted to make a second-order rigid network that can be assembled
by combining many identical components. In this case, we would like a critical
configuration where the edges lengths are as similar as possible. We can look
for such configurations by minimizing the relative fluctuations of the
lengths:

$\displaystyle V_{L}$ $\displaystyle=\frac{\bm{\langle}(L-\langle
L\rangle)^{2}\bm{\rangle}}{\langle L\rangle^{2}}=\frac{\left\langle
L^{2}\right\rangle}{\langle L\rangle^{2}}-1,$ (44)

where $\langle\rangle$ indicates an average over the edges of the network.
Similarly we could minimize the fluctuations of the putative edge tensions
$\tau_{\alpha}=\sigma^{h}_{\alpha}L_{\alpha}$, which are the edge lengths of
the dual force network.

$\displaystyle V_{\tau}$
$\displaystyle=\frac{\bm{\langle}(\tau-\langle\tau\rangle)^{2}\bm{\rangle}}{\langle\tau\rangle^{2}}=\frac{\left\langle\tau^{2}\right\rangle}{\left\langle\tau\right\rangle^{2}}-1=\frac{\langle(\sigma^{h}L)^{2}\rangle}{\left\langle\sigma^{h}L\right\rangle^{2}}-1.$
(45)

This is potentially useful for design, as if the edges of the network can only
support a finite amount of tension before breaking, a network with smaller
fluctuations in tensions would be more resistant to fracture.

We find that each network structure has a single unique critical configuration
that minimizes $V_{L}$, shown in Fig 3A. In contrast, there are many different
configurations that minimize $V_{\tau}$, and several examples with the same
$V_{\tau}$ are shown in Fig 3B,C. This is because $V_{\tau}$ depends only on
the dual force network, and two different configurations of the physical
network possess the same force network when they are related by an isogonal
mode which changes the lengths of the edges without changing the angles
between them noll_active_2017 ; brauns_geometric_2024 .

####  IV.1.3 Random Critical Configurations

We also study two populations of critical configurations from different random
samples of critical configurations. First, we randomly assign rest lengths to
the edges of the network and uniformly decrease them (equivalent to performing
a volumetric expansion) until the system hits the critical manifold and
becomes rigid. This is the standard method used to study strain stiffening in
spring and fiber network model. Examples are shown in Fig. 3J,K. We also
randomly assign the component of geometric stress associated with each edge
and use the geometric stress parameterization to obtain the corresponding
critical configuration. Details of these sampling procedures can be found in
the Supplemental Materials. Examples are shown in Fig. 3L,M.

Figure 6: Structural Metrics of Different Populations of Critical
Configurations: The horizontal axis labels $V_{L}$, $V_{\tau}$, $B_{0}$, and
$G_{0}$ refer to the objective function being optimized, while $\ell_{\alpha}$
and $\sigma^{h}_{\alpha}$ refer to the populations derived by randomly
assigning rest lengths and geometric stresses respectively. The images of
networks next to the violin plots show examples of network structures that
correspond with high and low values of each of the structural metrics.

###  IV.2 Analyzing Ensemble Network Properties

We seek to understand ensemble-averaged properties for these different types
of configurations. For each instance of 100 unique network structures
described by incidence matrix $g_{\alpha i}$ with $N_{v}=1000$ vertices, we
generate six populations of critical configurations: minimized $V_{L}$,
minimized $V_{\tau}$, maximized $B_{0}$, maximized $G_{0}$, random rest
lengths, and random geometric stresses. Note that for the purpose of comparing
the configurations that were optimized for $G_{0}$ with the other networks, we
use $r=0.1$ for the threshold on the geometric stress, an example of which is
shown in Fig. 3G.

We then compare these populations in order to investigate the full range of
behaviors and properties that these second-order rigid networks can achieve,
and look for emergent properties of the optimized networks that could
eventually be used to generate such structures using local rules instead of
global optimization.

####  IV.2.1 Local Structure

First, we compare the local structure of the configurations in both physical
space and in the dual force space. In Figure 6, we show for each population of
critical configurations violin plots of the distributions of several
structural parameters. First the objective function $V_{L}$, which quantifies
the fluctuations in the edge lengths. As this quantity does not measure the
orientations of the edges, we also plot the average shape anisotropy of the
faces of the networks $\langle\kappa_{\text{phys}}\rangle$. The shape
anisotropy of a face with vertices indexed by $i=1,..,n$ with coordinates
$r^{i}_{\mu}$ can be written in terms of the eigenvalues $\lambda_{1}$,
$\lambda_{2}$ of the gyration tensor
$S_{\mu\nu}=1/n\sum_{i}r^{i}_{\mu}r^{i}_{\nu}$:

$\displaystyle\kappa_{\text{phys}}=2\frac{\lambda_{1}^{4}+\lambda_{2}^{4}}{(\lambda_{1}^{2}+\lambda_{2}^{2})^{2}}-1$
(46)

This quantity is 0 when the face is a regular polygon, and close to 1 if the
face is elongated along a particular direction sarkar_shear-induced_2016 .
Note that the average $\langle\kappa_{\text{phys}}\rangle$ is taken over all
faces of a single member of a population, and the violin plots in 6 show the
distribution of this quantity over each ensemble. We then plot the equivalents
of the length fluctuations and shape anisotropy for the dual force networks
$V_{\tau}$ and $\langle\kappa_{\text{force}}\rangle$.

We see that for all metrics, networks with optimized $V_{L}$ are very similar
those with random rest lengths, which makes sense as the random rest lengths
are all chosen from the same Gaussian distribution. Our goal when optimizing
$V_{L}$ was to find a critical configuration with equal lengths, and thus a
variance $V_{L}=0$. With our chosen network structure, the optimization
procedure is able to find critical configurations with very small variances in
edge lengths on the order of $V_{L}\lesssim 10^{-7}$. In other words, for any
given incidence matrix in our ensemble, it is generically possible for us to
identify a single spring length $L^{*}$ so that if the network were assembled
of units of that length (with any type of nonlinear spring response), the
material would be precisely at the rigidity transition. Note that in this work
we study a limited ensemble of incidence matrices that are three-fold random
regular coordinated, with initial structure generated by a random sequential
addition protocol that further regularizes the network topology, and future
work could expand this ensemble.

Networks with optimized $B_{0}$ are also largely similar to those with
optimized $V_{\tau}$, as both have small variations in edge tensions, leading
to regular force networks and hence small values of
$\langle\kappa_{\text{force}}\rangle$. However, unlike the configurations that
were able to achieve negligible values of $V_{L}$, the networks with optimal
$V_{\tau}$ still have a non-zero variance in edge tensions $V_{\tau}\approx
10^{-2}$. This is because any network with equal edge tensions must have all
adjacent edges meet at $120\degree$ angles in order to maintain force balance,
and such an arrangement is only geometrically possible in a regular hexagonal
lattice and not in a generic 3-regular graph. The other populations of
configurations have edges that are more aligned, which causes a concentration
of tension. Configurations with optimized $G_{0}$ (with a geometric stress
tolerance of $r=0.1$) are highly aligned, and
$\langle\kappa_{\text{force}}\rangle$ is close to its maximum value of 1 as
the triangles that make up the dual force network are nearly one-dimensional.
This avoids the near-complete alignment of edges that occurs as this threshold
is decreased, as can be seen in 3H-I.

####  IV.2.2 Mechanical Response

We also characterize the mechanical behavior of the different populations of
critical configurations as they are strained away from the critical manifold.
Using the same ensembles as in Fig. 6, we start each configuration on the
critical manifold then subject the networks to a volumetric expansion by
increasing the volume as

$\displaystyle V(\gamma)=V(0)(1+\gamma)=N_{v}(1+\gamma),$ (47)

in strain steps of $\delta\gamma=0.001$. At each strain step we minimize the
elastic energy in Eq. (33) as a function of the vertex coordinates $x_{i\mu}$
using the FIRE algorithm until the largest component of the unbalanced forces
$\partial E/\partial x_{i\mu}$ has a magnitude less than $10^{-15}$. Using
linear response, we then calculate the bulk and shear modulus at each strain
step.

In Fig. 7A & C we plot the distributions of the bulk modulus $B_{0}$ and shear
modulus $G_{0}$ at the transition point. As expected, the networks that were
optimized for $B_{0}$ do in fact have the largest bulk moduli, followed
closely by the networks optimized for $V_{\tau}$. This maximum value of the
bulk modulus is about $10\%$ greater than that of networks that were assigned
random rest lengths or were optimized for $V_{L}$, and about $50\%$ higher
that of configurations optimized for $G_{0}$ or derived from a random
geometric stress. Networks optimized for shear modulus have significantly
larger $G_{0}$ than the other populations of approximately isotropic
configurations, which have negligible shear moduli due to finite size effects
that vanish in the thermodynamic limit as $G_{0}\sim N^{-1}$ as shown in Fig.
5.

We also investigate how the elastic moduli change as the networks are strained
past the rigidity transition. Close to the critical point, the moduli should
scale linearly with volumetric strain lee_stiffening_2022

$\displaystyle B(\gamma)$ $\displaystyle=B_{0}+C_{B}\gamma,$ $\displaystyle
G(\gamma)$ $\displaystyle=G_{0}+C_{G}\gamma.$ (48)

Fig. 7B & D show the distributions of the coefficients $C_{B}$ and $C_{G}$,
which were obtained from a least-squares fit of the moduli measured near the
transition. We observe that populations with larger moduli have smaller
scaling coefficients. Networks with optimized $B_{0}$ or $V_{\tau}$ have the
largest bulk moduli at the transition, but this decreases with additional
strain (i.e $C_{B}<0$), whereas the bulk modulus of the other types of
configurations increases with strain $(C_{B}>0)$. Similarly, networks
optimized for $G_{0}$ have their shear moduli decrease with strain, and all
other groups’ shear moduli increase past the transition point.

Figure 7: Mechanical Response Metrics of Different Populations of Critical
Configurations. Distributions of A) the bulk modulus $B_{0}$ at the transition
point, B) the bulk modulus strain scaling coefficient $C_{B}$, C) the shear
modulus $G_{0}$ at the transition point, and D) the shear modulus strain
scaling coefficient $C_{G}$ for the same six ensembles of critical
configurations as in Fig. 6. The red lines in B,D) are drawn at $C_{B,G}=0$ to
visually separate populations whose elastic moduli increase or decrease with
volumetric strain.

##  V Summary and Discussion

In this paper we have described a method to parameterize the space of all
critical configurations of any central force network. We showed that we can
use gradient descent to search the space of geometric stresses to find
critical configurations with particular structural properties (equal edge
lengths or tensions) or elastic responses (maximum bulk or shear modulus). For
the shear modulus, we found that global optimization pushes the system towards
extreme states that are not experimentally realizable, though it is
straightforward to constrain the minimization to reasonable states that
exhibit finite $G_{0}$ even in the large-$N$ limit, in contrast to random
disordered networks where $G_{0}$ vanishes.

The geometric stress parameterization allows us to enumerate all possible
rigid configurations of a given network structure. In this paper we have
focused on 3-coordinated random regular networks for simplicity, but future
work should investigate how changes to the topology of a network affect the
form of the critical manifold. This could be useful to study how
underconstrained networks fracture when edges with high tension are broken
driscoll2016role , or how biopolymer networks are remodelled by fiber-cleaving
enzymes which preferentially sever low-tension fibers saini2020tension .

One result that may be particularly useful for materials design is our
discovery that we can design materials perched precisely at the rigidity
transition that are composed of units of identical length (e.g. we can
optimize so that $V_{L}\sim 0$ for all network topologies in our limited
ensemble). This suggests that it may be possible to design self-assembled
materials where the network topology is programmed, e.g. by small DNA tags, to
link together ”bars” that are all a standard length. Because the network would
sit exactly at the edge of a rigidity transition, the resulting material could
be stiffened by orders of magnitude via very small changes to tuning
parameters, and this macroscopic, emergent property would be robust to small
differences in the mechanical properties of the individual bars. Moreover,
with small extensions to the work presented above, it should also be possible
to program materials that strain stiffen or become floppy after the material
experiences a precisely controlled amount of strain.

In overconstrained networks, the force transmission properties are relatively
divorced from the precise structural details, as any system will have many
different ways to distribute stress through the system. However, the geometric
stress parameterization shows that in underconstrained networks the structure
and force transmission are intimately connected, both arising from the
underlying geometry of the critical manifold. This has implications for the
possible states that underconstrained networks can access. For example,
several researchers have looked at emergent alignment of stress in
underconstrained networks that are strained to the critical point
abhilash2014remodeling ; taufalele2019fiber ; wang2014long . In future work we
plan to study how the geometry of the critical manifold makes aligned
structures highly probable.

Underconstrained networks also differ from overconstrained networks in that
their vibrational spectrum is separated into two distinct groups. There are
“longitudinal” modes that involve first-order deformations to the members of
the network and whose behavior are thus determined in part by material
properties, and there are “transverse” modes (the linear zero modes) which are
governed by geometry alone. It has been shown that in ordered networks with
sufficient symmetries, these populations of modes can be completely decoupled
rocks_integrating_2024 . The geometric stress parameterization gives a
framework for searching for such states in a disordered network, which would
then have interesting vibrational behaviour.

Mechanical networks have been trained to exhibit auxetic behavior
rocks_designing_2017 ; reid_auxetic_2018  as well as specific allostery-
inspired responses rocks_designing_2017 . However, as these learning methods
rely on the rigidity of the network, most studies of physical learning
algorithms in mechanical networks have focused on overconstrained networks.
The geometric stress parameterization should be very helpful in implementing
learning rules and global gradient descent that optimize desired objective
functions in underconstrained networks, and may therefore accelerate research
on physical learning in such systems.

##  VI Acknowledgements

MLM and TH were supported by Simons Foundation #454947, NSF-DMR-1951921 and
SU’s Orange Grid research computing cluster. CDS was supported by DMR-2217543.

## References

  * [1] Rubayn Goh, Scott P. O. Danielsen, Eric Schaible, Robert M. McMeeking, and J. Herbert Waite.  Nanolatticed Architecture Mitigates Damage in Shark Egg Cases.  Nano Letters, 21(19):8080–8085, October 2021.  Publisher: American Chemical Society.
  * [2] Yu-Ling Wei, Qing-Sheng Yang, Xia Liu, and Ran Tao.  Multi-bionic mechanical metamaterials: A composite of FCC lattice and bone structures.  International Journal of Mechanical Sciences, 213:106857, January 2022.
  * [3] Alessandro Mongera, Payam Rowghanian, Hannah J. Gustafson, Elijah Shelton, David A. Kealhofer, Emmet K. Carn, Friedhelm Serwane, Adam A. Lucio, James Giammona, and Otger Campàs.  A fluid-to-solid jamming transition underlies vertebrate body axis elongation.  Nature, 561(7723):401–405, September 2018.  Number: 7723 Publisher: Nature Publishing Group.
  * [4] Dapeng Bi, J. H. Lopez, J. M. Schwarz, and M. Lisa Manning.  A density-independent rigidity transition in biological tissues.  Nature Physics, 11(12):1074–1079, December 2015.  Number: 12 Publisher: Nature Publishing Group.
  * [5] Dapeng Bi, Xingbo Yang, M. Cristina Marchetti, and M. Lisa Manning.  Motility-Driven Glass and Jamming Transitions in Biological Tissues.  Physical Review X, 6(2):021011, April 2016.  Publisher: American Physical Society.
  * [6] F. C. MacKintosh, J. Käs, and P. A. Janmey.  Elasticity of Semiflexible Biopolymer Networks.  Physical Review Letters, 75(24):4425–4428, December 1995.  Publisher: American Physical Society.
  * [7] Cornelis Storm, Jennifer J. Pastore, F. C. MacKintosh, T. C. Lubensky, and Paul A. Janmey.  Nonlinear elasticity in biological gels.  Nature, 435(7039):191–194, May 2005.  Number: 7039 Publisher: Nature Publishing Group.
  * [8] A. Sharma, A. J. Licup, K. A. Jansen, R. Rens, M. Sheinman, G. H. Koenderink, and F. C. MacKintosh.  Strain-controlled criticality governs the nonlinear mechanics of fibre networks.  Nature Physics, 12(6):584–587, June 2016.  Number: 6 Publisher: Nature Publishing Group.
  * [9] Robbie Rens, Carlos Villarroel, Gustavo Düring, and Edan Lerner.  Micromechanical theory of strain stiffening of biopolymer networks.  Physical Review E, 98(6):062411, December 2018.  Publisher: American Physical Society.
  * [10] Michael Zaiser and Stefano Zapperi.  Disordered mechanical metamaterials.  Nature Reviews Physics, pages 1–10, September 2023.  Publisher: Nature Publishing Group.
  * [11] Matthias Merkel, Karsten Baumgarten, Brian P. Tighe, and M. Lisa Manning.  A minimal-length approach unifies rigidity in underconstrained materials.  Proceedings of the National Academy of Sciences, 116(14):6560–6568, April 2019.  Publisher: Proceedings of the National Academy of Sciences.
  * [12] Jesse Silverberg, Arthur Evans, Lauren McLeod, Ryan Hayward, Thomas Hull, Christian Santangelo, and Itai Cohen.  APPLIED ORIGAMI Using origami design principles to fold reprogrammable mechanical metamaterials.  Science (New York, N.Y.), 345:647–650, August 2014.
  * [13] Aleksandar Donev, Robert Connelly, Frank H. Stillinger, and Salvatore Torquato.  Underconstrained jammed packings of nonspherical hard particles: Ellipses and ellipsoids.  Physical Review E, 75(5):051304, May 2007.  Publisher: American Physical Society.
  * [14] Ludovic Berthier, Giulio Biroli, Patrick Charbonneau, Eric I. Corwin, Silvio Franz, and Francesco Zamponi.  Gardner physics in amorphous solids and beyond.  The Journal of Chemical Physics, 151(1):010901, July 2019.
  * [15] H. J. Schek.  The force density method for form finding and computation of general networks.  Computer Methods in Applied Mechanics and Engineering, 3(1):115–134, January 1974.
  * [16] Pier Malerba, Matteo Patelli, and Manuel Quagliaroli.  An Extended Force Density Method for the form finding of cable systems with new forms.  Structural engineering & mechanics, 42:191–210, March 2012.
  * [17] Hoang Chi Tran and Jaehong Lee.  Advanced form-finding of tensegrity structures.  Computers & Structures, 88(3):237–246, February 2010.
  * [18] Yafeng Wang, Xian Xu, and Yaozhi Luo.  Form-finding of tensegrity structures via rank minimization of force density matrix.  Engineering Structures, 227:111419, January 2021.
  * [19] Matthias Merkel and M. Lisa Manning.  A geometrically controlled rigidity transition in a model for confluent 3D tissues.  New Journal of Physics, 20(2):022002, February 2018.  Publisher: IOP Publishing.
  * [20] M. van Hecke.  Jamming of soft particles: geometry, mechanics, scaling and isostaticity.  Journal of Physics: Condensed Matter, 22(3):033101, December 2009.  Publisher: IOP Publishing.
  * [21] Robert Connelly and Simon D Guest.  Frameworks, tensegrities, and symmetry.  Cambridge University Press, 2022.
  * [22] Robert Connelly.  The rigidity of certain cabled frameworks and the second-order rigidity of arbitrarily triangulated convex surfaces.  Advances in Mathematics, 37(3):272–299, September 1980.
  * [23] R. Connelly and W. Whiteley.  Second-Order Rigidity and Prestress Stability for Tensegrity Frameworks.  SIAM J. Discret. Math., 1996.
  * [24] Ojan Khatib Damavandi, Varda F. Hagh, Christian D. Santangelo, and M. Lisa Manning.  Energetic rigidity. I. A unifying theory of mechanical stability.  Physical Review E, 105(2):025003, February 2022.  Publisher: American Physical Society.
  * [25] M. Berry, David Limberg, M. E. Lee-Trimble, Ryan Hayward, and C. D. Santangelo.  Controlling the configuration space topology of mechanisms.  Technical Report arXiv:2205.10180, arXiv, May 2022.  arXiv:2205.10180 [cond-mat] type: article.
  * [26] Simon Guest.  The stiffness of prestressed frameworks: A unifying approach.  International Journal of Solids and Structures, 43(3):842–854, February 2006.
  * [27] Julia A. Giannini, David Richard, M. Lisa Manning, and Edan Lerner.  Bond-space operator disentangles quasilocalized and phononic modes in structural glasses.  Physical Review E, 104(4):044905, October 2021.  Publisher: American Physical Society.
  * [28] Edan Lerner, Eric DeGiuli, Gustavo Düring, and Matthieu Wyart.  Breakdown of continuum elasticity in amorphous solids.  Soft Matter, 10(28):5085–5092, June 2014.  Publisher: The Royal Society of Chemistry.
  * [29] Julien Guénolé, Wolfram G. Nöhring, Aviral Vaid, Frédéric Houllé, Zhuocheng Xie, Aruna Prakash, and Erik Bitzek.  Assessment and optimization of the fast inertial relaxation engine (fire) for energy minimization in atomistic simulations and its implementation in lammps.  Computational Materials Science, 175:109584, April 2020.
  * [30] Craig Maloney and Anaël Lemaître.  Amorphous Systems in Athermal, Quasistatic Shear.  Physical review. E, Statistical, nonlinear, and soft matter physics, 74:016118, August 2006.
  * [31] Nicholas Noll, Madhav Mani, Idse Heemskerk, Sebastian J. Streichan, and Boris I. Shraiman.  Active tension network model suggests an exotic mechanical state realized in epithelial tissues.  Nature Physics, 13(12):1221–1226, December 2017.  Publisher: Nature Publishing Group.
  * [32] Fridtjof Brauns, Nikolas Claussen, Matthew F. Lefebvre, Eric Wieschaus, and Boris I. Shraiman.  The Geometric Basis of Epithelial Convergent Extension.  bioRxiv: The Preprint Server for Biology, page 2023.05.30.542935, May 2024.
  * [33] Sumantra Sarkar, Dapeng Bi, Jie Zhang, Jie Ren, R. P. Behringer, and Bulbul Chakraborty.  Shear-induced rigidity of frictional particles: Analysis of emergent order in stress space.  Physical Review E, 93(4):042901, April 2016.  Publisher: American Physical Society.
  * [34] C.-T. Lee and M. Merkel.  Stiffening of under-constrained spring networks under isotropic strain.  Soft Matter, 18(29):5410–5425, 2022.
  * [35] Michelle M Driscoll, Bryan Gin-ge Chen, Thomas H Beuman, Stephan Ulrich, Sidney R Nagel, and Vincenzo Vitelli.  The role of rigidity in controlling material failure.  Proceedings of the National Academy of Sciences, 113(39):10813–10817, 2016.
  * [36] Karanvir Saini, Sangkyun Cho, Lawrence J Dooling, and Dennis E Discher.  Tension in fibrils suppresses their enzymatic degradation–a molecular mechanism for ‘use it or lose it’.  Matrix Biology, 85:34–46, 2020.
  * [37] AS Abhilash, Brendon M Baker, Britta Trappmann, Christopher S Chen, and Vivek B Shenoy.  Remodeling of fibrous extracellular matrices by contractile cells: predictions from discrete fiber network simulations.  Biophysical journal, 107(8):1829–1840, 2014.
  * [38] Paul V Taufalele, Jacob A VanderBurgh, Adam Muñoz, Matthew R Zanotelli, and Cynthia A Reinhart-King.  Fiber alignment drives changes in architectural and mechanical features in collagen matrices.  Plos one, 14(5):e0216537, 2019.
  * [39] Hailong Wang, AS Abhilash, Christopher S Chen, Rebecca G Wells, and Vivek B Shenoy.  Long-range force transmission in fibrous matrices enabled by tension-driven alignment of fibers.  Biophysical journal, 107(11):2592–2603, 2014.
  * [40] Jason W. Rocks and Pankaj Mehta.  Integrating local energetics into Maxwell-Calladine constraint counting to design mechanical metamaterials, March 2024.  arXiv:2208.07419 [cond-mat].
  * [41] Jason W. Rocks, Nidhi Pashine, Irmgard Bischofberger, Carl P. Goodrich, Andrea J. Liu, and Sidney R. Nagel.  Designing allostery-inspired response in mechanical networks.  Proceedings of the National Academy of Sciences, 114(10):2520–2525, March 2017.  Publisher: Proceedings of the National Academy of Sciences.
  * [42] Daniel R. Reid, Nidhi Pashine, Justin M. Wozniak, Heinrich M. Jaeger, Andrea J. Liu, Sidney R. Nagel, and Juan J. de Pablo.  Auxetic metamaterials from disordered networks.  Proceedings of the National Academy of Sciences, 115(7):E1384–E1390, February 2018.  Publisher: Proceedings of the National Academy of Sciences.

##  VII Supplemental Materials

###  VII.1 Generating Random Underconstrained Networks

In the main text, we compared the emergent properties of four different
populations of critical configurations of 2 dimensional periodic central-force
networks. These networks were initially generated by making a Voronoi
tessellation using the algorithm in the Qhull library implemented in SciPy.
The sets of points used to produce the Voronoi diagrams were generated by a
random sequential addition procedure. The coordinates of the points are chosen
randomly from a unit square one at a time, and each successive point is
accepted only if it further than a distance $r$ away from all previously
accepted points. Varying the threshold distance $r$ changes the total number
of points that can be placed.

We use this procedure to produce networks that experience full rigidity
percolation under strain when all edges are given similar rest lengths, in
order to compare them with the fully rigid critical configurations obtained
from the geometric stress parameterization. If the seed points for the Voronoi
tessellation are chosen purely randomly, the resulting networks have large
fluctuations in the density of vertices, and regions with higher density then
become rigid at a larger critical strain than low-density regions when all
edges have similar rest lengths.

Once the Voronoi tessellation has been generated, giving us a network with
$N_{v}$ vertices and $N_{b}=3N_{v}/2$ edges, we describe the underlying graph
structure with the incidence matrix. After choosing an arbitrary orientation
for each edge, the incidence matrix $g_{\alpha i}$ is given by

$\displaystyle g_{\alpha i}=\begin{cases}+1\quad\text{if vertex }i\;\text{is
the head of edge }\alpha\\\ -1\quad\text{if vertex }i\;\text{is the tail of
edge }\alpha\\\ \;\;0\;\quad\text{else}\end{cases}$ (49)

If the $\mu=x,y$ coordinate of vertex $i=1,..,N_{v}$ is $x_{i\mu}$, we can
write the $\mu$ component of the vector along edge $\alpha$ as

$\displaystyle L_{\alpha\mu}=\sum_{i}g_{\alpha i}x_{i\mu}+b_{\alpha\mu},$ (50)

where $b_{\alpha\mu}$ encodes the periodic boundaries as shown in Figure 8.
The light blue interior edges have $b_{\alpha\mu}=0$, while the dark blue
edges that cross the periodic boundary have non-zero $b_{\alpha\mu}$. For
example, the edge connecting the principle domain to the periodic copy
directly to the right would have $b_{\alpha\mu}=[S,0]$ where $S$ is the linear
dimension of the periodic box, which we take to be $S=\sqrt{N_{b}}$ so that
the area of the box scales with the number of vertices and edges. If we want
to apply an affine deformation $A_{\mu\nu}$ to the network, we can simply
apply it to the boundary term: $b_{\alpha\mu}\rightarrow
A_{\mu\nu}b_{\alpha\nu}$.

Figure 8: A) A small periodic network. Light blue edges connect two vertices
within the principle domain, with their periodic copies in light green. Dark
blue edges connect vertices in different periodic copies, crossing the
periodic boundary (periodic copies shown in dark green). B) The same network
with all vertices placed at the same point, demonstrating the boundary term
$b_{\alpha\mu}$.

###  VII.2 Analytic Expressions

####  VII.2.1 Deriving Rigidity, Geometric Response, and Prestress Matrices

We can now derive the expressions for the geometric response matrix and the
prestress matrix by taking derivatives of the quantity
$h_{\alpha}=L_{\alpha}^{2}/2$.

$\displaystyle R^{h}_{\alpha i\mu}$ $\displaystyle=\frac{\partial
h_{\alpha}}{\partial x_{i\mu}}=\frac{1}{2}\sum_{\nu}\frac{\partial
L_{\alpha\nu}^{2}}{\partial x_{i\mu}}=\sum_{\nu}\frac{\partial
L_{\alpha\nu}}{\partial x_{i\mu}}L_{\alpha\nu}$ $\displaystyle=g_{\alpha
i}\delta_{\mu\nu}L_{\alpha\nu}=g_{\alpha i}L_{\alpha\mu},$ (51) $\displaystyle
P_{i\mu j\nu}$
$\displaystyle=\sum_{\alpha}\sigma^{h}_{\alpha}\frac{\partial^{2}h_{\alpha}}{\partial
x_{i\mu}\partial
x_{j\nu}}=\sum_{\alpha}\sigma^{h}_{\alpha}\frac{\partial}{\partial
x_{j\nu}}\left(g_{\alpha i}L_{\alpha\mu}\right)$
$\displaystyle=\sum_{\alpha}\sigma^{h}_{\alpha}g_{\alpha i}g_{\alpha
j}\delta_{\mu\nu}=P_{ij}\delta_{\mu\nu}.$ (52)

As an illustrative example, let us derive these matrices for the Hookean
constraint $f_{\alpha}=L_{\alpha}-\ell_{\alpha}$:

$\displaystyle\tilde{R}_{\alpha i\mu}$
$\displaystyle=\frac{\partial\tilde{f}_{\alpha}}{\partial
x_{i\mu}}=\frac{\partial L_{\alpha}}{\partial
x_{i\mu}}=\frac{1}{2L_{\alpha}}\frac{\partial L_{\alpha}^{2}}{\partial
x_{i\mu}}L_{\alpha\nu}=g_{\alpha i}\frac{L_{\alpha\mu}}{L_{\alpha}},$ (53)
$\displaystyle\tilde{P}_{i\mu j\nu}$
$\displaystyle=\sum_{\alpha}\tau_{\alpha}\frac{\partial^{2}\tilde{f}_{\alpha}}{\partial
x_{i\mu}\partial x_{j\nu}}=\sum_{\alpha}\tau_{\alpha}\frac{\partial}{\partial
x_{j\nu}}\left(g_{\alpha i}\frac{L_{\alpha\mu}}{L_{\alpha}}\right)$
$\displaystyle=\sum_{\alpha}\tau_{\alpha}g_{\alpha
i}\left(\frac{1}{L_{\alpha}}\frac{\partial L_{\alpha\mu}}{\partial
x_{j\nu}}-\frac{L_{\alpha\mu}}{L_{\alpha}^{2}}\frac{\partial
L_{\alpha}}{\partial x_{j\nu}}\right)$
$\displaystyle=\sum_{\alpha}\frac{\tau_{\alpha}}{L_{\alpha}}g_{\alpha
i}g_{\alpha
j}\delta_{\mu\nu}-\frac{\tau_{\alpha}}{L_{\alpha}^{3}}L_{\alpha\mu}g_{\alpha
i}g_{\alpha j}L_{\alpha\nu}$ $\displaystyle=P_{i\mu
j\nu}-\sum_{\alpha}\tilde{R}_{\alpha i\mu}\tau_{\alpha}\tilde{R}_{\alpha
j\nu}.$ (54)

Since $R$ and $\tilde{R}$ are identical up to a rescaling of the rows by
$L_{\alpha}$, they have the same right nullspace, so the linear zero modes are
exactly the same. The prestress matrix derived from the Hookean constraint
$\tilde{P}$ is the same as the prestress matrix of the squared lengths up to
an additional term that depends on the rigidity matrix. However, because the
rigidity of the network is determined by the projection of the prestress
matrix onto the subspace of linear zero modes (which are by definition in the
nullspace of $R$), this extra term adds no new information about the rigidity
of the network.

####  VII.2.2 Geometry of the Critical Manifold

In order to calculate total derivatives of an objective function, we need to
write the derivatives of the geometric stress parameterization. Let us write
the bond vector $L_{\beta\mu}$ as a function of $\sigma^{h}_{\alpha}$ using
the geometric stress parameterization:

$\displaystyle L_{\beta\mu}$ $\displaystyle=\sum_{i}g_{\beta
i}x_{i\mu}+b_{\beta\mu}$ (55) $\displaystyle=\sum_{i}g_{\beta
i}\left[\sum_{j\gamma}-\left(P^{-1}\right)_{ij}g_{\gamma
j}\sigma^{h}_{\gamma}b_{\gamma\mu}\right]+b_{\beta\mu}$ (56)
$\displaystyle=-\sum_{ij\gamma}g_{\beta i}\left(P^{-1}\right)_{ij}g_{\gamma
j}\sigma^{h}_{\gamma}b_{\gamma\mu}+b_{\beta\mu},$ (57)

where by $P$ we mean the non-zero $N_{v}\times N_{v}$ block
$P_{ij}=\sum_{\alpha}\sigma_{\alpha}^{h}g_{\alpha i}g_{\alpha j}$ of the full
$2N_{v}\times 2N_{v}$ prestress matrix $\bm{P}_{i\mu
j\nu}=P_{ij}\delta_{\mu\nu}$. Taking a derivative with respect to
$\sigma^{h}_{\alpha}$ gives

$\displaystyle\frac{\partial L_{\beta\mu}}{\partial\sigma^{h}_{\alpha}}$
$\displaystyle=\frac{\partial}{\partial\sigma^{h}_{\alpha}}\left[-\sum_{ij\gamma}g_{\beta
i}\left(P^{-1}\right)_{ij}g_{\gamma
j}\sigma^{h}_{\gamma}b_{\gamma\mu}+b_{\beta\mu}\right]$ (58)
$\displaystyle=\sum_{ijkl\gamma}g_{\beta
i}\left(P^{-1}\right)_{ij}\frac{\partial
P_{jk}}{\partial\sigma^{h}_{\alpha}}\left(P^{-1}\right)_{kl}g_{\gamma
l}\sigma^{h}_{\gamma}b_{\gamma\mu}-\sum_{ij}g_{\beta
i}\left(P^{-1}\right)_{ij}g_{\alpha j}b_{\alpha\mu}$ (59)
$\displaystyle=\left[\sum_{ij}g_{\beta i}\left(P^{-1}\right)_{ij}g_{\alpha
j}\right]\left[\sum_{kl\gamma}g_{\alpha k}\left(P^{-1}\right)_{kl}g_{\gamma
l}\sigma^{h}_{\gamma}b_{\gamma\mu}\right]-\left[\sum_{ij}g_{\beta
i}\left(P^{-1}\right)_{ij}g_{\alpha j}\right]b_{\alpha\mu}$ (60)
$\displaystyle=-\left[\sum_{ij}g_{\beta i}\left(P^{-1}\right)_{ij}g_{\alpha
j}\right]\left[-\sum_{kl\gamma}g_{\alpha k}\left(P^{-1}\right)_{kl}g_{\gamma
l}\sigma^{h}_{\gamma}b_{\gamma\mu}+b_{\alpha\mu}\right]=-\left[\sum_{ij}g_{\beta
i}\left(P^{-1}\right)_{ij}g_{\alpha j}\right]L_{\alpha\mu},$ (61)

where we have used the standard matrix identity

$\displaystyle\frac{\partial}{\partial x}P^{-1}=-P^{-1}\frac{\partial
P}{\partial x}P^{-1}.$ (62)

Now we can write the derivative of the squared lengths
$h_{\beta}=L_{\beta}^{2}/2$:

$\displaystyle\frac{\partial
h_{\beta}}{\partial\sigma^{h}_{\alpha}}=\sum_{\mu}L_{\beta\mu}\frac{\partial
L_{\beta\mu}}{\partial\sigma^{h}_{\alpha}}$ (63)
$\displaystyle=-\sum_{\mu}L_{\beta\mu}\left[\sum_{ij\mu}g_{\beta
i}\left(P^{-1}\right)_{ij}g_{\alpha j}\right]L_{\alpha\mu}.$ (64)

Since the geometric response matrix is $R^{h}_{\alpha i\mu}=g_{\alpha
i}L_{\alpha\mu}$, we can further write this as the matrix product

$\displaystyle\frac{\partial
h_{\beta}}{\partial\sigma^{h}_{\alpha}}=-\left(R^{h}\bm{P}^{-1}{R^{h}}^{T}\right)_{\alpha\beta}.$
(65)

Then $\sum_{\beta}\sigma^{h}_{\beta}(\partial
h_{\beta}/\partial\sigma^{h}_{\alpha})$ must vanish because
$\sigma^{h}_{\beta}$ is guaranteed to be in the left nullspace of $R^{h}$. As
discussed in the main text, the columns of this matrix span the tangent space
of the critical manifold $h_{\beta}(\bm{\sigma}^{h})$ at the current
configuration, so this implies that the geometric stress is normal to the
critical manifold.

We can also calculate the second fundamental form of the critical manifold,
which describes the local curvature and is given by

$\displaystyle\mathrm{I\\!I}_{\beta\gamma}=\sum_{\alpha}\frac{\partial^{2}h_{\alpha}}{\partial
s_{\beta}\partial
s_{\gamma}}\frac{\sigma^{h}_{\alpha}}{|\bm{\sigma}^{h}|}=\frac{1}{|\bm{\sigma}^{h}|}(V^{T}BV)_{\beta\gamma},$
(66)

where $s_{\beta}$ are the $N_{b}-1$ stereographic coordinates as described in
the main text and $\sigma^{h}_{\alpha}/|\bm{\sigma}^{h}|$ is the local unit
normal vector for the manifold. Defining
$V_{\alpha\beta}:=\partial\sigma^{h}_{\alpha}/\partial s_{\beta}$ we have

$\displaystyle\mathrm{I\\!I}_{\beta\gamma}$
$\displaystyle=\frac{1}{|\bm{\sigma}^{h}|}\sum_{\delta\epsilon}\frac{\partial\sigma^{h}_{\delta}}{\partial
s_{\beta}}\left(\sum_{\alpha}\sigma^{h}_{\alpha}\frac{\partial^{2}h_{\alpha}}{\partial\sigma^{h}_{\delta}\partial\sigma^{h}_{\epsilon}}\right)\frac{\partial\sigma^{h}_{\epsilon}}{\partial
s_{\beta}}$ (67)
$\displaystyle=\frac{1}{|\bm{\sigma}^{h}|}(V^{T}BV)_{\beta\gamma},$ (68)

where
$B_{\beta\gamma}=\sum_{\alpha}\sigma^{h}_{\alpha}(\partial^{2}h_{\alpha}/\partial\sigma^{h}_{\beta}\partial\sigma^{h}_{\gamma})$
satisfies the interesting equality

$\displaystyle B_{\beta\gamma}$
$\displaystyle=\sum_{\alpha}\sigma^{h}_{\alpha}\frac{\partial^{2}h_{\alpha}}{\partial\sigma^{h}_{\beta}\partial\sigma^{h}_{\gamma}}=\sum_{\alpha}\sigma^{h}_{\alpha}\frac{\partial}{\partial\sigma^{h}_{\beta}}\left(-R^{h}P^{-1}{R^{h}}^{T}\right)_{\alpha\gamma}$
(69)
$\displaystyle=-\sum_{\alpha}\sigma^{h}_{\alpha}\left[\left(\frac{\partial
R^{h}}{\partial\sigma^{h}_{\beta}}\right)P^{-1}{R^{h}}^{T}+R^{h}\left(\frac{\partial
P^{-1}}{\partial\sigma^{h}_{\beta}}\right){R^{h}}^{T}+R^{h}P^{-1}\left(\frac{\partial{R^{h}}^{T}}{\partial\sigma^{h}_{\beta}}\right)\right]_{\alpha\gamma}$
(70)
$\displaystyle=-\sum_{ij\mu}\sum_{\alpha}\sigma^{h}_{\alpha}\left[\frac{\partial}{\partial\sigma^{h}_{\beta}}\left(g_{\alpha
i}L_{\alpha\mu}\right)\right]\left(P^{-1}\right)_{ij}g_{\gamma
j}L_{\gamma\mu}=\sum_{ij\mu}\sum_{\alpha}\sigma^{h}_{\alpha}g_{\alpha
i}\left[\frac{\partial
L_{\alpha\mu}}{\partial\sigma^{h}_{\beta}}\right]\left(P^{-1}\right)_{ij}g_{\gamma
j}L_{\gamma\mu}$ (71)
$\displaystyle=\sum_{ij\mu}\sum_{\alpha}\sigma^{h}_{\alpha}g_{\alpha
i}\left[\sum_{kl}g_{\alpha k}\left(P^{-1}\right)_{kl}g_{\beta
l}L_{\beta\mu}\right]\left(P^{-1}\right)_{ij}g_{\gamma j}L_{\gamma\mu}$ (72)
$\displaystyle=\sum_{lj\mu}(g_{\beta
l}L_{\beta\mu})\left[\left(P^{-1}\right)_{kl}\left(\sum_{\alpha}\sigma^{h}_{\alpha}\;g_{\alpha
i}g_{\alpha k}\right)\left(P^{-1}\right)_{ij}\right](g_{\gamma
j}L_{\gamma\mu})$ (73) $\displaystyle=\sum_{lj\mu}R^{h}_{\beta
l\mu}\left(P^{-1}PP^{-1}\right)_{lj}R^{h}_{\gamma
j\mu}=\left(R^{h}\bm{P}^{-1}{R^{h}}^{T}\right)_{\beta\gamma}=\frac{\partial
h_{\beta}}{\partial\sigma^{h}_{\gamma}},$ (74)

using Eqs. 65 & 64 and the fact that by definition the geometric stress
$\sigma_{\alpha}^{h}$ is in the nullspace of the geometric response matrix
$R^{h}$. So for the critical manifold the matrix $B$ determines both the local
curvature through the second fundamental form as well as the gradients
$\partial h_{\beta}/\partial\sigma_{\gamma}^{h}$, which give the local tangent
space and hence the metric.

####  VII.2.3 Elastic Moduli at the Rigidity Transition

We want to write the elastic moduli of a network as a simple function of the
network configuration $L_{\alpha\mu}$ and the geometric stress
$\sigma^{h}_{\alpha}$ in order to take derivatives and perform gradient
descent. Consider a global affine deformation $I+\gamma A_{\mu\nu}$, which
changes the boundary term and the bond vectors as

$\displaystyle\frac{\partial b_{\alpha\mu}}{\partial\gamma}$
$\displaystyle=\sum_{\nu}A_{\mu\nu}b_{\alpha\nu}$ (75)
$\displaystyle\frac{\partial L_{\alpha\mu}}{\partial\gamma}$
$\displaystyle=\sum_{\nu}A_{\mu\nu}L_{\alpha\nu},$ (76)

which results in a change to the squared lengths
$h_{\alpha}=\sum_{\mu}L_{\alpha\mu}^{2}/2$

$\displaystyle\dot{h}:=\frac{\partial
h_{\alpha}}{\partial\gamma}=\sum_{\mu}L_{\alpha\mu}\frac{\partial
L_{\alpha\mu}}{\partial\gamma}=\sum_{\mu\nu}L_{\alpha\mu}A_{\mu\nu}L_{\alpha\nu}.$
(77)

The elastic modulus associated with this deformation can be calculated as

$\displaystyle\mu_{A}=\frac{1}{V}\frac{\text{d}^{2}E}{\text{d}\gamma^{2}}=\frac{1}{V}\left(\frac{\partial^{2}E}{\partial\gamma^{2}}-\bm{\Xi}^{T}H^{-1}\bm{\Xi}\right),$
(78)

where $V$ is the system volume (or area for 2D systems). The first term
contains the contributions from the affine motions of the system, while the
second term accounts for the non-affine relaxation that occurs in order to
achieve force balance. $H=\partial^{2}E/\partial x_{i}\partial x_{j}$ is the
Hessian matrix and $\Xi_{i}=\partial^{2}E/\partial x_{i}\partial\gamma$
contains the unbalanced forces that appear as a result of the affine
deformation. Note that $H^{-1}$ should be understood as the pseudo-inverse of
$H$, which projects out any zero modes. Let us calculate these quantities at
the second order rigidity transition, where a self-stress exists but the
prestress is zero $(\partial E/\partial h_{\alpha}=0)$. First the affine term:

$\displaystyle\frac{\partial^{2}E}{\partial^{2}\gamma}$
$\displaystyle=\frac{\partial}{\partial\gamma}\left(\sum_{\alpha}\frac{\partial
E}{\partial h_{\alpha}}\frac{\partial h_{\alpha}}{\partial\gamma}\right)$ (79)
$\displaystyle=\sum_{\alpha\beta}\dot{h}_{\alpha}\frac{\partial^{2}E}{\partial
h_{\alpha}\partial h_{\beta}}\dot{h}_{\beta}+\sum_{\alpha}\frac{\partial
E}{\partial h_{\alpha}}\frac{\partial^{2}h_{\alpha}}{\partial\gamma^{2}}$
$\displaystyle=\dot{\bm{h}}^{T}K\dot{\bm{h}},$

where $\dot{h}_{\alpha}=\partial h_{\alpha}/\partial\gamma$ gives the change
in squared lengths if the network due to the chosen affine deformation and
$K_{\alpha\beta}=\partial^{2}E/\partial h_{\alpha}\partial h_{\beta}$ is the
matrix of modified stiffnesses. For a network of Hookean springs with energy
$E=\sum_{\alpha}k_{\alpha}(L_{\alpha}-\ell_{\alpha})^{2}/2$ at zero prestress,
this is a positive-definite diagonal matrix:
$K_{\alpha\beta}=k_{\alpha}/L_{\alpha}^{2}\;\delta_{\alpha\beta}$. Next, at
the transition the Hessian is given by just the Gram term

$\displaystyle H={R^{h}}^{T}KR^{h},$ (80)

and the non-affine forces are

$\displaystyle\Xi_{i\mu}$
$\displaystyle=\frac{\partial^{2}E}{\partial\gamma\partial
x_{i\mu}}=\frac{\partial}{\partial\gamma}\left(\sum_{\alpha}\frac{\partial
E}{\partial h_{\alpha}}\frac{\partial h_{\alpha}}{\partial x_{i\mu}}\right)$
(81) $\displaystyle=\sum_{\alpha\beta}\frac{\partial h_{\alpha}}{\partial
x_{i\mu}}\frac{\partial^{2}E}{\partial h_{\alpha}\partial
h_{\beta}}\frac{\partial
h_{\beta}}{\partial\gamma}=({R^{h}}^{T}K\dot{\bm{h}})_{i\mu}.$

Then the non-affine contribution to the elastic modulus is

$\displaystyle\bm{\Xi}^{T}H^{-1}\bm{\Xi}=\dot{\bm{h}}^{T}KR^{h}({R^{h}}^{T}KR^{h})^{-1}{R^{h}}^{T}K\dot{\bm{h}}.$
(82)

Since $K$ is positive definite and diagonal we can scale it out by defining

$\displaystyle\bar{R}=\sqrt{K}R^{h}$ (83)

so that

$\displaystyle\bm{\Xi}^{T}H^{-1}\bm{\Xi}=\dot{\bm{h}}^{T}\sqrt{K}\bar{R}(\bar{R}^{T}\bar{R})^{-1}\bar{R}^{T}\sqrt{K}\dot{\bm{h}}.$
(84)

If the network is not in a critical configuration then $R$ (and hence
$\bar{R}$) are full rank, which implies that

$\displaystyle\bar{R}(\bar{R}^{T}\bar{R})^{-1}\bar{R}^{T}$ $\displaystyle=I$
(85)

so that

$\displaystyle\bm{\Xi}^{T}H^{-1}\bm{\Xi}$
$\displaystyle=\dot{\bm{h}}^{T}K\dot{\bm{h}}.$ (86)

Then if the network does not possess a self-stress, the affine part of the
elastic modulus is exactly cancelled by the non-affine part so that all
elastic moduli vanish. However, if a network has any self-stresses then $R$
and $\bar{R}$ are not full rank, so those directions get projected out of the
non-affine term. Assuming only a single geometric stress $\bm{\sigma}$ in the
left nullspace of $R$ for simplicity, we have

$\displaystyle\bar{R}(\bar{R}^{T}\bar{R})^{-1}\bar{R}^{T}$
$\displaystyle=I-\frac{\bar{\bm{\sigma}}{\bar{\bm{\sigma}}}^{T}}{|\bar{\bm{\sigma}}|^{2}},$
(87)

where $\bar{\bm{\sigma}}=\sqrt{K}^{-1}\bm{\sigma}$ is in the left nullspace of
$\bar{R}$, since

$\displaystyle\bar{R}^{T}\bar{\bm{\sigma}}=\left({R^{h}}^{T}\sqrt{K}\right)\left(\sqrt{K}^{-1}\bm{\sigma}\right)={R^{h}}^{T}\bm{\sigma}=0.$
(88)

Then

$\displaystyle\bm{\Xi}^{T}H^{-1}\bm{\Xi}$
$\displaystyle=\dot{\bm{h}}^{T}\sqrt{K}\left(I-\frac{\bar{\sigma}\bar{\sigma}^{T}}{|\bar{\sigma}|^{2}}\right)\sqrt{K}\dot{\bm{h}}$
(89)
$\displaystyle=\dot{\bm{h}}^{T}K\dot{\bm{h}}-\frac{\left(\dot{\bm{h}}^{T}\bm{\sigma}\right)^{2}}{\bm{\sigma}^{T}K^{-1}\bm{\sigma}}$
(90)
$\displaystyle=\dot{\bm{h}}^{T}K\dot{\bm{h}}-\frac{\left(\sum_{\alpha}\sigma_{\alpha}\dot{h}_{\alpha}\right)^{2}}{\sum_{\alpha}\sigma_{\alpha}^{2}/K_{\alpha\alpha}},$
(91)

so the full elastic modulus is

$\displaystyle\mu_{A}=\frac{1}{V}\frac{\left(\sum_{\alpha}\sigma_{\alpha}\dot{h}_{\alpha}\right)^{2}}{\sum_{\alpha}\sigma_{\alpha}^{2}/K_{\alpha\alpha}}=\frac{1}{V}\frac{\left(\sum_{\alpha}\sigma_{\alpha}\dot{h}_{\alpha}\right)^{2}}{\sum_{\alpha}\sigma_{\alpha}^{2}L_{\alpha}^{2}/k_{\alpha}}.$
(92)

###  VII.3 Numerical Methods for Optimization and Random Sampling on the
Critical Manifold

####  VII.3.1 Randomly Sampled Configurations

First, we sample the critical manifold by assigning random rest lengths to the
networks and subjecting them to a bulk volumetric strain until they reach the
critical manifold and become rigid. We sample the rest lengths from a Gaussian
distribution centered at $\ell_{\alpha}=1$ with standard deviation $s=0.3$,
which was chosen to give the largest amount of variance in rest lengths while
ensuring the networks experience full rigidity percolation at the critical
strain.

Once the rest lengths are chosen, we uniformly rescale them as
$\ell_{\alpha}(\gamma)=\ell_{\alpha}/(1+\gamma)$ (equivalent to a bulk
expansion), increasing the strain $\gamma$ from 0 in steps of 0.01. We do this
quasistatically, where after each strain step we minimize the energy of the
network $E=1/2\sum_{\alpha}(L_{\alpha}-\ell_{\alpha})^{2}$. For this we use
the FIRE gradient descent algorithm, stopping the minimization once all
components of the gradient $\partial E/\partial x_{i\mu}$ are less than a
tolerance of $10^{-16}$. We strain each network until the total energy is at
least $10^{-4}$, then record the configuration $x_{i\mu}$ and the geometric
stress $\sigma^{h}_{\alpha}=(L_{\alpha}-\ell_{\alpha})/L_{\alpha}$, which we
then normalize to $\sum_{\alpha}(\sigma_{\alpha}^{h})^{2}=N_{b}$ in order to
compare with the geometric stresses used in the parameterization.

Rather than randomly assigning rest lengths to the edges, we can uniformly
sample in the space of positive-definite geometric stresses and use the
geometric stress parameterization to get the corresponding critical
configurations. Since geometric stresses are defined only up to an overall
scaling, we take each component of the geometric stress from a standard normal
distribution: $\sigma^{h}_{\alpha}\sim\mathcal{N}(0,1)$, which results in a
uniform distribution over the hypersphere of unique geometric stresses. To
ensure the resulting critical configurations are rigid, and to compare
directly with the strained configurations, we then take the absolute value of
the geometric stress ($\sigma^{h}_{\alpha}\rightarrow|\sigma^{h}_{\alpha}|$),
which is a uniform sample over the positive definite orthant of geometric
stress space. Once the geometric stress is chosen, we form the prestress
matrix $P_{ij}=\sum_{\alpha}\sigma^{h}_{\alpha}g_{\alpha i}g_{\alpha j}$ then
remove the first row and column, which eliminates the trivial zero mode by
effectively pinning the first vertex at the origin. We then perform an LU
decomposition on $P$ and solve for the coordinates of the other vertices using
the geometric stress parameterization described in the main text.

####  VII.3.2 Optimized Configurations

We also generate four populations of critical configurations by optimizing
different objective functions ($V_{L}$, $V_{\tau}$, $B_{0}$, $G_{0}$) over the
space of geometric stresses. To optimize an objective function
$O(\bm{\sigma}^{h})$ we use the FIRE gradient descent algorithm, treating the
geometric stress as the degrees of freedom, for which we need to calculate the
gradient $\partial O/\partial\sigma^{h}_{\alpha}$. The normalization of the
geometric stress is not relevant for the parameterization, but we constrain
the normalization to $\sum_{\alpha}|\sigma^{h}_{\alpha}|^{2}=N_{b}$ so that
each $\sigma^{h}_{\alpha}$ is $\mathcal{O}(1)$ to prevent any numerical
issues. We stop the optimization when the magnitudes of all components of the
gradient are less than a threshold value $10^{-12}$.