# Partition sum of thermal, under-constrained systems

Cheng-Tai Lee  Aix Marseille Univ, Université de Toulon, CNRS, CPT (UMR 7332),
Turing Center for Living Systems, Marseille, France  Matthias Merkel
<EMAIL_ADDRESS>Aix Marseille Univ, Université de Toulon, CNRS, CPT
(UMR 7332), Turing Center for Living Systems, Marseille, France

###### Abstract

Athermal (i.e. zero-temperature) under-constrained systems are typically
floppy, but they can be rigidified by the application of external strain.
Following our recently developed analytical theory for the athermal limit,
here and in the companion paper, we extend this theory to under-constrained
systems at finite temperatures. Close to the athermal transition point, we
derive from first principles the partition sum for a broad class of under-
constrained systems, from which we obtain analytic expressions for elastic
material properties such as isotropic tension $t$ and shear modulus $G$ in
terms of isotropic strain $\varepsilon$, shear strain $\gamma$, and
temperature $T$. Our work unifies the physics of systems as diverse as polymer
fibers & networks, membranes, and vertex models for biological tissues.

##  I Introduction

Under-constrained systems possess more degrees of freedom (dofs),
$N_{\mathrm{dof}}$, than constraints (or springs), $N_{s}$. In the athermal
limit, i.e. at zero temperature, under-constrained systems are generally
floppy. Conversely, over-constrained systems, which have
$N_{\mathrm{dof}}<N_{s}$, are generally rigid [1, 2]. We recently developed a
generic analytic framework to predict the elastic properties of athermal,
under-constrained systems [3, 4]. However, such a framework does not yet exist
for finite temperatures.

In the athermal limit, under-constrained systems are typically floppy due to
the existence of zero modes. For instance, there are _infinitesimal_ zero
modes, i.e. collective displacements of the dofs that do not affect any of the
constraints _to linear order_. The number of linearly independent
infinitesimal zero modes, $N_{0}$, is given by [1, 2]:

$N_{0}=N_{\mathrm{dof}}-N_{s}+N_{SSS},$ (1)

where $N_{SSS}$ is the number of linearly independent states of self-stress
(SSS). A SSS is a combination of virtual tensions put on the springs of a
network that does not result in any net force on the dofs. According to Eq.
(1), under-constrained systems always have a positive number of infinitesimal
zero modes.

There are three possibilities regarding the nature of infinitesimal zero
modes. First, they could be _trivial_ zero modes, i.e. global translations or
rotations of the system, which leave the “shape” of the system unchanged.
Second, they could correspond to _mechanisms_ , i.e.  _finite_ displacements
that change the “shape” of the system, but do not affect any of the
constraints [5, 6]. Third, while inducing no changes to the constraints to
linear order, collective displacements could induce changes to higher order,
corresponding to higher-order rigidity [7, 8]. Here, we refer to this third
kind as _$K$ th-order zero modes_ with some integer $K>0$. Such modes do not
affect the constraints to $K$th order, but do so to $(K+1)$th order (precise
definitions in subsection II.4). Because many of the infinitesimal zero modes
often correspond to mechanisms, at zero temperature, under-constrained systems
are typically floppy.

Still, under-constrained systems can be rigidified by externally applied
strain, such as isotropic expansion or shear [2, 9, 10, 3, 7, 8, 11]. This is
because at some amount of externally applied strain, the constraints cannot be
fulfilled any more. For a spring network, this happens when the amount of
external strain becomes geometrically incompatible with the spring rest
lengths. Right at the rigidity transition, a SSS is created [5]. This allows
prestresses to form once the network is strained any further. Beyond the
transition point both SSS and prestresses rigidify the system [2, 9, 12, 5,
13, 3]. In previous work, we used these insights to analytically predict the
elastic material properties close to the rigidity transition in the athermal
limit [3, 4]. These ideas apply not only to spring networks, but also to other
systems such as polymer networks [14, 15, 16, 17] and vertex models for
biological tissues [18, 13, 19, 3, 20].

Despite these advances for the athermal limit, there is still relatively
little known about _thermal_ , i.e. finite-temperature, under-constrained
systems [21, 22, 23]. Recent work used effective-medium theory (EMT) to obtain
analytical expressions for the shear modulus of spring networks [24, 25], and
Mannattil et al. discuss the effect of singularities in the configuration
space on the mechanics of the system [6]. However, we still lack a generic
theory for the elastic properties of thermal, under-constrained systems.

Here, we analytically derive from first principles the partition sum and
elastic properties of thermal, under-constrained systems close to the athermal
transition for varying isotropic strain $\varepsilon$, shear strain $\gamma$,
and small temperatures $T$. After defining the class of under-constrained
systems that we discuss here (section II), we first revisit the athermal limit
(section III). This will provide the foundations for discussing at finite
temperature both the limit of infinitely stiff springs (section IV) and the
general case of finite spring stiffness (section V). In the companion paper
111Reference to be inserted later., we numerically verify the analytical
results obtained here.

##  II Model

###  II.1 Network of generalized springs

We consider a very general class of Hookean spring networks that display a
strain-induced rigidity transition in the athermal limit.

The network is embedded in $D$ spatial dimensions with periodic boundary
conditions with system volume $\mathcal{V}$. The periodic box can be sheared
with shear strain variable $\gamma$, where our framework applies both to
simple and pure shear. In appendix A, we show examples for the implementation
of both kinds of shear.

A network consists of $N_{\mathrm{node}}$ nodes whose position components
combine into a $N_{\mathrm{dof}}=DN_{\mathrm{node}}$-dimensional vector
$\vec{R}$. The nodes are connected by $N_{s}$ generalized springs with lengths
$L_{i}$, labeled by $i=1,\dots,N_{s}$. We consider an under-constrained spring
network, where $N_{\mathrm{dof}}>N_{s}$.

Each spring length depends on node positions, system size, and shear through
an arbitrary function $L_{i}(\vec{R},\mathcal{V},\gamma)$ that only needs to
fulfill the following two criteria:

  1. 1.

In the vicinity of the athermal transition, the function
$L_{i}(\vec{R},\mathcal{V},\gamma)$ is _analytic with respect to all three
parameters_.

  2. 2.

The spring length $L_{i}$ _has dimensionality $d_{i}$_ in the following sense:
For some non-negative integer $d_{i}$, the spring length function
$L_{i}(\vec{R},\mathcal{V},\gamma)$ is homogeneous of degree $d_{i}$ with
respect to an isotropic rescaling. This means that for any real rescaling
factor $c>0$:

$c^{d_{i}}L_{i}(\vec{R},\mathcal{V},\gamma)=L_{i}(c\vec{R},c^{D}\mathcal{V},\gamma).$
(2)

For example, angular springs have $d_{i}=0$, ordinary linear springs or
perimeter springs in the vertex model have $d_{i}=1$ (appendix A), and area
and volume springs in the vertex model have $d_{i}=2$ and $d_{i}=3$,
respectively.

Note that the functions $L_{i}$ are allowed to substantially differ among the
generalized springs of the network. For example, also models that combine 2D
or 3D vertex models with spring networks are accounted for [27].

The system energy is:

$E=\frac{1}{2}\sum_{i=1}^{N_{s}}{K_{i}\Big{[}L_{i}(\vec{R},\mathcal{V},\gamma)-L_{0i}\Big{]}^{2}},$
(3)

where we introduced spring constants $K_{i}$ and rest lengths $L_{0i}$.
Throughout this article, we treat all $K_{i}$ and $L_{0i}$ as constant and
study the behavior of the $L_{i}$ depending on system size $\mathcal{V}$,
shear $\gamma$, and temperature $T$.

###  II.2 Athermal transition and isotropic strain

We consider a network that in the athermal limit transitions from floppy (zero
shear modulus) to rigid (finite shear modulus) at some critical volume. For
spring networks or vertex models, this is typically the case when the network
percolates the periodic box (in the sense of connectivity percolation). For
zero shear strain, $\gamma=0$, we denote this critical volume by
$\mathcal{V}^{\ast}$, where, consistent with our earlier findings [3, 4], we
assume here that the system is floppy for $\mathcal{V}<\mathcal{V}^{\ast}$ and
rigid for $\mathcal{V}>\mathcal{V}^{\ast}$.

We define the isotropic strain $\varepsilon$ with respect to the critical
volume $\mathcal{V}^{\ast}$:

$\varepsilon:=\frac{1}{D}\log{\left(\frac{\mathcal{V}}{\mathcal{V}^{\ast}}\right)}.$
(4)

Thus, for $\gamma=0$, the transition occurs at $\varepsilon=0$.

For simplicity, we assume here and in the following that there are no states
of self-stress (SSS) in the floppy regime. This implies that in the floppy
regime and at the transition, each spring attains its rest length,
$L_{i}=L_{0i}$. Also, we assume that only a single SSS forms at the
transition, which is generally true for disordered networks.

We denote the values of the dofs $\vec{R}$ at the transition by
$\vec{R}^{\ast}$. For simplicity, we assume here that $\vec{R}^{\ast}$ is
uniquely defined. A degeneracy in the athermal transition configuration
$\vec{R}^{\ast}$ could be created either by global translations or by
mechanisms. Global translations do not change any of our arguments in this
paper, and we deliberate more on the possibility of a $\vec{R}^{\ast}$
degeneracy by mechanisms in the discussion section.

###  II.3 Dimensionless quantities

We introduce dimensionless dofs $\vec{r}$ and dimensionless spring lengths
$\ell_{i}$ as:

$\displaystyle\vec{r}:=\frac{1}{\sqrt[D]{\mathcal{V}}}\vec{R}=\frac{1}{e^{\varepsilon}\sqrt[D]{\mathcal{V}^{\ast}}}\vec{R},$
(5)
$\displaystyle\ell_{i}:=e^{-d_{i}\varepsilon}\frac{L_{i}(\vec{R},\mathcal{V},\gamma)}{L_{0i}}.$
(6)

Using Eqs. (5) and (2) with
$c=1/\sqrt[D]{\mathcal{V}}=1/(e^{\varepsilon}\sqrt[D]{\mathcal{V}^{\ast}})$,
we obtain from Eq. (6):

$\ell_{i}(\vec{r},\gamma)=\frac{(\mathcal{V}^{\ast})^{d_{i}/D}}{L_{0i}}L_{i}(\vec{r},1,\gamma).$
(7)

In other words, when expressed using the dimensionless node positions
$\vec{r}$, the dimensionless spring lengths do not depend on isotropic strain
any more.

As discussed in subsection II.2, for $\gamma=0$, the transition occurs at
$\varepsilon=0$ and $\vec{r}=\vec{r}^{\ast}$, and each spring attains its rest
length, $L_{i}=L_{0i}$. From Eq. (6) thus follows that at the transition point
for $\gamma=0$:

$\ell_{i}\big{(}\vec{r}=\vec{r}\,^{\ast},\gamma=0\big{)}=1.$ (8)

Thus, the homogeneity property in Eq. (2) allows us to drop the need for a
common rest lengths and equal dimensions of all springs, which was required in
our earlier work [3, 4]. Instead, the common rest length of all rescaled
spring lengths is just one.

For convenience, we define the dimensionless weights

$w_{i}:=\sqrt{\frac{K_{i}L_{0i}^{2}}{E_{0}}},$ (9)

where we have introduced the energy scale
$E_{0}:=(\sum_{i}{K_{i}L_{0i}^{2}})/N_{s}$ such that
$\sum_{i}w_{i}^{2}=N_{s}$. This allows us to rewrite the system energy as

$E=\frac{E_{0}}{2}\sum_{i=1}^{N_{s}}{w_{i}^{2}\big{(}e^{d_{i}\varepsilon}\ell_{i}-1\big{)}^{2}}.$
(10)

Our approach of exploiting spring length homogeneity with respect to spatial
rescaling, Eq. (2), thus allows us to include the overall isotropic strain
explicitly into the energy functional, Eq. (17). This will help us to obtain
our analytical results more directly than in our previous work [3, 4].

###  II.4 Zero modes

We close this section with some definitions of different kinds of zero modes
in the system, as they will turn out to be important for its elastic behavior.

As usual [2, 5], for some configuration $(\vec{r},\varepsilon,\gamma)$, we
call a vector $\vec{z}\in\mathbb{R}^{N_{\mathrm{dof}}}$ _infinitesimal zero
mode_ whenever for all $i=1,\dots,N_{s}$:

$\frac{\partial\ell_{i}(\vec{r},\gamma)}{\partial r_{n}}z_{n}=0.$ (11)

Here, $r_{n}$ and $z_{n}$ with $n=1,\dots,N_{\mathrm{dof}}$ are the components
of the vectors $\vec{r}$ and $\vec{z}$, respectively, and we imply summation
over repeated indices. For a given configuration, we denote the dimension of
the vector space of infinitesimal zero modes by $N_{0}$.

Any given infinitesimal zero mode is either trivial, it belongs to a
mechanism, or it belongs to a $K$th-order zero mode. For our periodic boundary
conditions, we call an infinitesimal zero mode $\vec{z}$ _trivial_ whenever it
corresponds to one of the $D$ global translations.

A _mechanism_ is a finite displacement of the dofs during which all $N_{s}$
springs keep their rest lengths. More formally, we say an infinitesimal zero
mode $\vec{z}$ belongs to a mechanism whenever it is not trivial and there is
a continuous function
$\vec{x}:u\mapsto\vec{x}(u),\mathbb{R}\rightarrow\mathbb{R}^{N_{\mathrm{dof}}}$
with $\vec{x}(0)=\vec{r}$ and some range $\eta>0$ such that for all
$i=1,\dots,N_{s}$:

$\displaystyle\vec{x}^{\prime}(0)$ $\displaystyle=\vec{z}$ and (12)
$\displaystyle e^{d_{i}\varepsilon}\ell_{i}\big{(}\vec{x}(u),\gamma\big{)}$
$\displaystyle=1$ for $|u|<\eta$, (13)

where the prime denotes a derivative with respect to $u$. The second
condition, Eq. (13), is equivalent to $L_{i}=L_{0i}$ (compare Eq. (6)). For
given $\varepsilon$ and $\gamma$, we denote the number of mechanisms of a
system by $N_{\mathrm{mech}}$ 222 More precisely, we define the number of
mechanisms $N_{\mathrm{mech}}$ of a system as the maximal integer $N_{m}$ for
which there is a continuous function
$\vec{x}:\vec{u}\mapsto\vec{x}(\vec{u}),\mathbb{R}^{N_{m}}\rightarrow\mathbb{R}^{N_{\mathrm{dof}}}$
and some $\eta>0$ such that Eq. (13) holds within some open ball of radius
$\eta$, i.e. for $|\vec{u}|<\eta$. . As discussed in subsection II.2, we focus
for simplicity on a system where _at the transition point_ ,
$\vec{r}=\vec{r}\,^{\ast}$ and $\varepsilon=\gamma=0$, there are no
mechanisms, $N_{\mathrm{mech}}^{\ast}=0$. Note however, that this does not
preclude the existence of mechanisms away from the transition point.

A _$K$ th-order zero mode_ does not affect the spring length to $K$th order,
but to $(K+1)$th order. More precisely, we say the infinitesimal zero mode
$\vec{z}$ belongs to a $K$th-order zero mode with $K>0$ whenever there is a
$K$ times differentiable function
$\vec{x}:u\mapsto\vec{x}(u),\mathbb{R}\rightarrow\mathbb{R}^{N_{\mathrm{dof}}}$
with $\vec{x}(0)=\vec{r}$ such that for all $i=1,\dots,N_{s}$:

$\vec{x}^{\prime}(0)=\vec{z}\qquad\text{and}\qquad\left.\frac{\mathrm{d}^{K}\ell_{i}\big{(}\vec{x}(u),\gamma\big{)}}{\mathrm{d}u^{K}}\right|_{u=0}=0,$
(14)

_while these conditions are not fulfilled for $K+1$._ I.e. there is _no_ $K+1$
times differentiable function
$\vec{y}:u\mapsto\vec{y}(u),\mathbb{R}\rightarrow\mathbb{R}^{N_{\mathrm{dof}}}$
with $\vec{y}(0)=\vec{r}$ such that for all $i=1,\dots,N_{s}$:

$\vec{y}^{\prime}(0)=\vec{z}\qquad\text{and}\qquad\left.\frac{\mathrm{d}^{K+1}\ell_{i}\big{(}\vec{y}(u),\gamma\big{)}}{\mathrm{d}u^{K+1}}\right|_{u=0}=0.$
(15)

We stress that due to the last condition, a mechanism is never a $K$th-order
zero mode. For the same reason, while a first-order zero mode always
corresponds to an infinitesimal zero mode, the converse is not always true.
For some configuration $(\vec{r},\varepsilon,\gamma)$, we denote the number of
$K$th-order zero modes by $N_{\mathrm{1st}},N_{\mathrm{2nd}},\dots$.
Specifically, at the transition point, $\vec{r}=\vec{r}\,^{\ast}$ and
$\varepsilon=\gamma=0$, we denote those numbers by
$N_{\mathrm{1st}}^{\ast},N_{\mathrm{2nd}}^{\ast},\dots$.

While there is generally a finite number of first-order zero modes at the
transition, $N_{\mathrm{1st}}^{\ast}$ (see Fig. 1 for examples), we expect the
occurrence of $K$th-order zero modes with $K>1$ to be an exception in
disordered systems. For simplicity we will focus in this article on systems
without such modes at the transition point, $N_{K\mathrm{th}}^{\ast}=0$ for
$K>1$. In appendix F, we discuss how our results can be affected by the
presence of such modes.

##  III Athermal limit ($T\rightarrow 0$)

While in our previous work we derived the elastic system properties in the
athermal limit based on a minimal-length function [3, 4], we present here an
alternative derivation, which will prepare our discussions in the subsequent
sections.

We will minimize the energy varying the node positions, $\vec{r}$, but keeping
the strain variables $\varepsilon,\gamma$ fixed. To this end, we first derive
an explicit expression of the energy in terms of the dofs $\vec{r}$ and the
strain variables $\varepsilon$ and $\gamma$ in the vicinity of the transition.

###  III.1 Explicit expression of the system energy in terms of dofs and
strain

To obtain an explicit expression for the energy in Eq. (10) in terms of the
dofs $\vec{r}$ and strain $\varepsilon,\gamma$ in the vicinity of the
transition, we first define the deviation of the spring lengths from their
transition point values (Eq. (8)) as

$\Delta\ell_{i}(\vec{r},\gamma):=\ell_{i}-1.$ (16)

This allows us to rewrite the energy, Eq. (10), to lowest order in
$\varepsilon$ as:

$E=\frac{E_{0}}{2}\sum_{i=1}^{N_{s}}{w_{i}^{2}\big{(}\Delta\ell_{i}+d_{i}\varepsilon\big{)}^{2}}.$
(17)

In deriving this expression from Eq. (10), we have factored out the prefactor
$e^{d_{i}\varepsilon}$ from the parentheses. Because the parentheses scale to
lowest order as $\sim\varepsilon$, the resulting prefactor of
$e^{2d_{i}\varepsilon}$ only adds higher-order terms and is thus neglected.

To obtain an explicit expression for the $\Delta\ell_{i}$ in Eq. (17) in terms
of the dofs $\vec{r}$, we expand the dimensionless spring lengths $\ell_{i}$
around the transition point $\vec{r}\,^{\ast}$ and for small $\gamma$.
Denoting the components of $\vec{r}$ by $r_{n}$ with
$n=1,\dots,N_{\mathrm{dof}}$, we can write to second order in $\Delta
r_{n}:=r_{n}-r_{n}^{\ast}$ and $\gamma$:

$w_{\underaccent{\dot}{i}}\Delta\ell_{\underaccent{\dot}{i}}=C_{in}\Delta
r_{n}+\frac{1}{2}M_{imn}\Delta r_{m}\Delta
r_{n}+B_{i}^{(1)}\gamma+\frac{1}{2}B_{i}^{(2)}\gamma^{2}.$ (18)

Here and in the following, we imply summation over repeated indices, except
for indices with underdots. In the absence of explicit sum symbols, implicit
sums over the index $i$ always run over all springs ($i=1,\dots,N_{s}$), and
implicit sums over the indices $m,n$ always run over all dofs
($m,n=1,\dots,N_{\mathrm{dof}}$). On the right-hand side,
$C_{in}:=w_{\underaccent{\dot}{i}}(\partial\ell_{\underaccent{\dot}{i}}/\partial
r_{n})$ is a generalized compatibility matrix and the matrix
$M_{imn}:=w_{\underaccent{\dot}{i}}(\partial^{2}\ell_{\underaccent{\dot}{i}}/\partial
r_{m}\partial r_{n})$ is symmetric in the last two indices. We introduce the
shear strain dependence through the coefficients
$B_{i}^{(1)}:=w_{\underaccent{\dot}{i}}(\partial\ell_{\underaccent{\dot}{i}}/\gamma)$
and
$B_{i}^{(2)}:=w_{\underaccent{\dot}{i}}(\partial^{2}\ell_{\underaccent{\dot}{i}}/\gamma^{2})$.
All derivatives are evaluated at the transition point,
$\vec{r}=\vec{r}\,^{\ast}$, and $\gamma=0$. We found that mixed terms in
$\Delta\vec{r}$ and $\gamma$ turn out not to be important to lowest order in
our final results, both here and in the following sections (appendix C). There
are no $\varepsilon$-dependent terms in Eq. (18), because the $\ell_{i}$
functions are independent of $\varepsilon$ (subsection II.3).

Equations (17) and (18) provide an explicit expression of the system energy in
terms of dofs $\vec{r}$, and strain variables $\varepsilon,\gamma$. To
simplify minimizing this energy with respect to the dofs, we now rewrite
$\Delta\ell_{i}$ and the energy $E$ in terms of eigenmodes of the
compatibility matrix $C_{in}$. To this end, we first perform a singular-value
decomposition (SVD) on it:

$C_{in}=\sum_{p=1}^{N_{s}}{U_{ip}s_{p}V_{np}},$ (19)

where $U_{ip}$ and $V_{np}$ are orthogonal square matrices of dimensions
$N_{s}\times N_{s}$ and $N_{\mathrm{dof}}\times N_{\mathrm{dof}}$,
respectively. The $N_{s}$ singular values $s_{p}$ are sorted in decreasing
order: $s_{1}\geq\dots\geq s_{N_{s}}$. Because at the transition a single SSS
is created, we have $s_{N_{s}-1}>0$ and $s_{N_{s}}=0$, where $U_{iN_{s}}$ is
the corresponding SSS. Furthermore, all $V_{np}$ with fixed
$p=N_{s},\dots,N_{\mathrm{dof}}$ are infinitesimal zero modes.

Based on the SVD, Eq. (19), we define changes of spring lengths and of degrees
of freedom in the left and right eigenbases of $C_{in}$, respectively:

$\displaystyle\Delta\tilde{\ell}_{p}$
$\displaystyle:=\sum_{i}{U_{ip}w_{i}\Delta\ell_{i}}$ for $p=1,\dots,N_{s}$,
(20) $\displaystyle\Delta\tilde{r}_{p}$ $\displaystyle:=V_{np}\Delta r_{n}$
for $p=1,\dots,N_{\mathrm{dof}}$. (21)

As a consequence of Eq. (19), $s_{N_{s}}=0$, and Eq. (21), all
$\Delta\tilde{r}_{p}$ with $p=N_{s},\dots,N_{\mathrm{dof}}$ are infinitesimal
zero modes.

We can now express the Taylor expansion, Eq. (18), in terms of the eigenmodes
of $C$. For $p=1,\dots,N_{s}$:

$\Delta\tilde{\ell}_{p}=s_{\underaccent{\dot}{p}}\Delta\tilde{r}_{\underaccent{\dot}{p}}+\frac{1}{2}\tilde{M}_{pmn}\Delta\tilde{r}_{m}\Delta\tilde{r}_{n}+\tilde{B}_{p}^{(1)}\gamma+\frac{1}{2}\tilde{B}_{p}^{(2)}\gamma^{2},$
(22)

where we have also transformed $M$ and $B^{(1/2)}$ into the eigenbases of $C$:
$\tilde{M}_{pmn}:=U_{ip}M_{im^{\prime}n^{\prime}}V_{m^{\prime}m}V_{n^{\prime}n}$
and $\tilde{B}_{p}^{(1/2)}:=U_{ip}B_{i}^{(1/2)}$.

We can now express the system energy, Eq. (17), in terms of the eigenmodes of
$C$:

$\begin{split}E&=\frac{E_{0}}{2}\sum_{i=1}^{N_{s}}{\Big{(}w_{i}\Delta\ell_{i}+d_{i}w_{i}\varepsilon\Big{)}^{2}}\\\
&=\frac{E_{0}}{2}\sum_{i,j=1}^{N_{s}}{\Big{(}w_{i}\Delta\ell_{i}+d_{i}w_{i}\varepsilon\Big{)}\delta_{ij}\Big{(}w_{j}\Delta\ell_{j}+d_{j}w_{j}\varepsilon\Big{)}}\\\
&=\frac{E_{0}}{2}\sum_{i,j,p=1}^{N_{s}}{\Big{(}w_{i}\Delta\ell_{i}+d_{i}w_{i}\varepsilon\Big{)}U_{ip}U_{jp}\Big{(}w_{j}\Delta\ell_{j}+d_{j}w_{j}\varepsilon\Big{)}}\\\
&=\frac{E_{0}}{2}\sum_{p=1}^{N_{s}}{\Big{(}\Delta\tilde{\ell}_{p}+\tilde{w}_{p}\varepsilon\Big{)}^{2}}.\end{split}$
(23)

In the second step of Eq. (23), we introduce the Kronecker delta
$\delta_{ij}$, and in the third step, we have used the orthogonality of $U$:
$U_{ip}U_{jp}=\delta_{ij}$. In the fourth step we have used Eq. (20) and
introduced $\tilde{w}_{p}:=\sum_{i=1}^{N_{s}}{U_{ip}d_{i}w_{i}}$. Note that
Eq. (19) is invariant with respect to a sign flip in the SSS, whose components
are given by $U_{iN_{s}}$ for each spring $i$. For convenience, we choose this
sign such that $\tilde{w}_{N_{s}}=\sum_{i=1}^{N_{s}}{U_{iN_{s}}d_{i}w_{i}}\geq
0$.

To more clearly express the dependency of the energy on the dofs
$\Delta\tilde{r}_{n}$ and strain variables $\varepsilon,\gamma$, we define for
$p=1,\dots,N_{s}$:

$\begin{split}\Lambda_{p}&:=\Delta\tilde{\ell}_{p}+\tilde{w}_{p}\varepsilon\\\
&=s_{\underaccent{\dot}{p}}\Delta\tilde{r}_{\underaccent{\dot}{p}}+\frac{1}{2}\tilde{M}_{pmn}\Delta\tilde{r}_{m}\Delta\tilde{r}_{n}\\\
&\qquad\qquad\qquad\qquad+\tilde{w}_{p}\varepsilon+\tilde{B}_{p}^{(1)}\gamma+\frac{1}{2}\tilde{B}_{p}^{(2)}\gamma^{2}.\end{split}$
(24)

In the second step, we used Eq. (22). Inserting the definition of
$\Lambda_{p}$ into the energy, Eq. (23), we obtain:

$E=\frac{E_{0}}{2}\sum_{p=1}^{N_{s}}\Lambda_{p}^{2}.$ (25)

Thus, each $\Lambda_{p}$ indicates how much the corresponding left eigenmode
$p$ of $C$ contributes to the system’s energy.

###  III.2 Energy minimum

We now discuss the minimum of the energy given by Eqs. (24) and (25) for fixed
$\varepsilon,\gamma$ while varying the dofs $\Delta\tilde{r}_{n}$. For
clarity, we focus here on the case where $\tilde{B}_{p}^{(1)}=0$ for all $p$.
The general case of arbitrary $\tilde{B}_{p}^{(1)}$ is discussed in appendix
C.

At $\varepsilon=\gamma=0$, the system is at its transition point
$\vec{r}=\vec{r}\,^{\ast}$ (subsection II.2), which implies
$\Delta\tilde{r}_{n}=0$ and thus, using Eq. (24), $\Lambda_{p}=0$. However, as
$\varepsilon$ and $\gamma$ become nonzero, also the $\Delta\tilde{r}_{n}$ and
$\Lambda_{p}$ generally become nonzero. We provide the detailed arguments for
the scaling of $\Lambda_{p}$ and $\Delta\tilde{r}_{n}$ with $\varepsilon$ and
$\gamma$ in appendix B.1, and provide here a brief sketch of these arguments.
Minimizing the energy in Eq. (25) can be achieved by minimizing the absolute
values of all the $\Lambda_{p}$. Indeed, to lowest relevant order in
$\varepsilon$ and $\gamma$, the values of the $\Lambda_{p}$ with $p<N_{s}$ can
be set to zero by adjusting the value of the corresponding
$\Delta\tilde{r}_{p}$ to (compare Eq. (102) in appendix B.1):

$\begin{split}\Delta\tilde{r}_{p}=-\frac{1}{s_{\underaccent{\dot}{p}}}\Bigg{(}&\frac{1}{2}\sum_{m,n=N_{s}}^{N_{\mathrm{dof}}}{\tilde{M}_{{\underaccent{\dot}{p}}mn}\Delta\tilde{r}_{m}\Delta\tilde{r}_{n}}\\\
&\qquad\qquad\qquad+\tilde{w}_{\underaccent{\dot}{p}}\varepsilon+\frac{1}{2}\tilde{B}_{\underaccent{\dot}{p}}^{(2)}\gamma^{2}\Bigg{)}.\end{split}$
(26)

However, $\Lambda_{N_{s}}$ cannot be set to zero in this way, because
$s_{N_{s}}=0$. As a consequence, the energy is dominated by the
$\Lambda_{N_{s}}$ contribution, $E=E_{0}\Lambda_{N_{s}}^{2}/2$, where,
according to Eq. (24) we have:

$\Lambda_{N_{s}}=\frac{1}{2}\sum_{m,n=N_{s}}^{N_{\mathrm{dof}}}{\tilde{M}_{N_{s}mn}\Delta\tilde{r}_{m}\Delta\tilde{r}_{n}}+\tilde{w}_{N_{s}}\varepsilon+\frac{1}{2}\tilde{B}_{N_{s}}^{(2)}\gamma^{2}.$
(27)

In the sum, the contributions by the terms with $m<N_{s}$ or $n<N_{s}$ can be
neglected to lowest order, since they only generate higher order terms in
$\varepsilon$ and $\gamma$ (Eq. (26) and appendix B.1).

We know that for $\gamma=0$ and $\varepsilon>0$, where the system is in the
rigid regime, the energy minimum is strictly positive,
$E=E_{0}\Lambda_{N_{s}}^{2}/2>0$. This is only possible if the matrix
$\tilde{M}_{N_{s}mn}$ (restricted to $m,n=N_{s},\dots,N_{\mathrm{dof}}$) is
positive semi-definite; otherwise $\Lambda_{N_{s}}$ could become zero. As a
consequence, for the athermal rigid regime,
$\varepsilon+b_{\varepsilon}\gamma^{2}>0$, the energy is minimal if
$\Delta\tilde{r}_{n}=0$ for all $N_{s}\leq n\leq N_{\mathrm{dof}}$, and it has
the value

$E=\frac{E_{0}\tilde{w}_{N_{s}}^{2}}{2}\Big{[}\varepsilon+b_{\varepsilon}\gamma^{2}\Big{]}^{2},$
(28)

where we introduced

$b_{\varepsilon}:=\frac{B_{N_{s}}^{(2)}}{2\tilde{w}_{N_{s}}}\qquad\text{if
{\hbox{\tilde{B}^{(1)}_p=0}} for all {\hbox{p}}}.$ (29)

In the more general case where the $\tilde{B}_{p}^{(1)}$ do not vanish, the
definition of $b_{\varepsilon}$ includes more terms, while Eq. (28) remains
the same up to a shift in $\gamma$ (appendix C).

Note that the matrix $\tilde{M}_{N_{s}mn}$ is closely linked to the Hessian of
the system and the emergence of rigidity in the athermal limit (appendix D).

###  III.3 Comparison to earlier results

The expression we published earlier for the energy of a network with only
$d$-dimensional springs in the rigid regime corresponds to [3, 4]:

$E=\frac{E_{0}N_{s}}{2(L_{0}^{\ast})^{2}(1+a_{\ell}^{2})}\Big{[}-\Delta
L_{0}+b\gamma^{2}\Big{]}^{2},$ (30)

where $\Delta L_{0}=-dL_{0}^{\ast}\varepsilon$, and $L_{0}^{\ast}$,
$a_{\ell}$, and $b$ are constants that depend on the network structure. We
included here the prefactor of 1/2 in our energy definition as compared to [3,
4]. Comparing Eqs. (28) and (30), we find:

$\displaystyle\tilde{w}_{N_{s}}$
$\displaystyle=d\left[\frac{N_{s}}{1+a_{\ell}^{2}}\right]^{1/2}$ (31)
$\displaystyle b_{\varepsilon}$ $\displaystyle=\frac{b}{dL_{0}^{\ast}}.$ (32)

In addition to our earlier work [3, 4], our approach here also clarifies why
the result is analytic in $\gamma$; it is inherited from the analytic nature
of the spring length functions $L_{i}(\vec{R},\mathcal{V},\gamma)$.

###  III.4 Elastic properties

In the rigid regime of the athermal limit, i.e. for
$\varepsilon+b_{\varepsilon}\gamma^{2}>0$, tension $t_{E}=\partial
E/\partial\mathcal{V}=(\partial E/\partial\varepsilon)/D\mathcal{V}$ and shear
modulus $G_{E}=(\partial^{2}E/\partial\gamma^{2})/\mathcal{V}$ are to lowest
order in $\varepsilon$ and $\gamma$:

$\displaystyle t_{E}$
$\displaystyle=\kappa_{E}\Big{[}\varepsilon+b_{\varepsilon}\gamma^{2}\Big{]}$
(33) $\displaystyle G_{E}$
$\displaystyle=2Db_{\varepsilon}\kappa_{E}\Big{[}\varepsilon+3b_{\varepsilon}\gamma^{2}\Big{]}$
(34)
$\displaystyle=2Db_{\varepsilon}t\left(1+\frac{2b_{\varepsilon}\gamma^{2}}{|\varepsilon+b_{\varepsilon}\gamma^{2}|}\right),$
(35)

where we defined

$\displaystyle\kappa_{E}$
$\displaystyle:=\frac{E_{0}\tilde{w}_{N_{s}}^{2}}{D\mathcal{V}^{\ast}}.$ (36)

The index $E$ indicates purely energetic rigidity.

##  IV Stiff-spring limit ($K_{i}\rightarrow\infty$)

In the limit of infinitely stiff springs but at finite temperature $T$,
elasticity is created purely by entropic effects. For
$\varepsilon+b_{\varepsilon}\gamma^{2}>0$ (rigid regime in the athermal
limit), there are no configurations where all springs can attain their
respective rest lengths. Thus, this regime is inaccessible in the stiff-spring
limit. Conversely, for $\varepsilon+b_{\varepsilon}\gamma^{2}\leq 0$ (floppy
regime in the athermal limit) there _are_ configurations where the springs can
attain their rest lengths. In this section, we derive an expression for the
accessible configurational phase space volume $\Omega$, from which we can
directly derive the elastic network properties. We focus again on the limit of
small $\varepsilon$ and $\gamma$.

###  IV.1 Configurational phase space volume

In appendix E we explicitly take the limit of the partition sum $Z$ for
$K_{i}\rightarrow\infty$ to obtain an expression for the configurational phase
space volume $\Omega$ 333In the stiff-spring limit, one might naively replace
each spring by an exact constraint and compute the partition sum using
Lagrangian mechanics. However, it has been shown that such an approach
generally leads to a different result than taking the partition sum $Z$ of a
system with finite spring stiffnesses $K_{i}$ and taking the limit
$K_{i}\rightarrow\infty$. [37]. The latter is of course the physically more
relevant approach, since there are no exact geometric constraints in nature..
Up to a constant factor, it is given by:

$\Omega\sim\int{\left(\prod_{n=1}^{N_{\mathrm{dof}}}{\mathrm{d}R_{n}}\right)\prod_{i=1}^{N_{\mathrm{s}}}{\delta(L_{i}-L_{0i})}},$
(37)

where $R_{n}$ for $n=1,\dots,N_{\mathrm{dof}}$ are the components of
$\vec{R}$. Using the dimensionless node positions and spring lengths, Eqs. (5)
and (6), this transforms into

$\Omega\sim\int{\left(\prod_{n=1}^{N_{\mathrm{dof}}}{\mathrm{d}r_{n}}\right)\prod_{i=1}^{N_{\mathrm{s}}}{\delta(\ell_{i}-e^{-d_{i}\varepsilon})}}.$
(38)

Here and in the following, we ignore any prefactors in $\Omega$ that are
powers of $e^{\varepsilon}$. This is because these factors only create terms
$\sim\varepsilon$ in the free energy $F_{S}\sim\log{\Omega}$. For small
$\varepsilon$, such terms are negligible compared to the leading-order term in
the free energy, which, as we will show, scales as $\sim\log{(-\varepsilon)}$
(Eq. (54)).

We aim to obtain an analytical expression for $\Omega$ in the vicinity of the
transition point, i.e. for $|\varepsilon|\ll 1$ and $|\gamma|\ll 1$. We thus
rewrite $\Omega$ as:

$\Omega\sim\int{\left(\prod_{n=1}^{N_{\mathrm{dof}}}{\mathrm{d}\Delta
r_{n}}\right)\prod_{i=1}^{N_{\mathrm{s}}}{\delta\Big{(}w_{i}\Delta\ell_{i}+d_{i}w_{i}\varepsilon\Big{)}}}.$
(39)

Here, we have included a constant factor of $w_{i}$ in each of the Dirac
deltas, and we again ignored prefactors of $e^{\varepsilon}$ in $\Omega$.

To better understand which set of points $\vec{r}$ is selected by the Dirac
deltas in Eq. (39), we again use the eigenmodes of the compatibility matrix.
To this end, we apply a rotation by $V_{np}$ to the integration variables
$\Delta r_{n}$, and a rotation by $U_{ip}$ on the argument of the multi-
dimensional Dirac delta, neither of which yields an additional factor in
$\Omega$:

$\begin{split}\Omega\sim&\int{\left(\prod_{n=1}^{N_{\mathrm{dof}}}{\mathrm{d}\Delta\tilde{r}_{n}}\right)\prod_{p=1}^{N_{s}}{\delta\Big{(}\Delta\tilde{\ell}_{p}+\tilde{w}_{p}\varepsilon\Big{)}}}\\\
&=\int{\left(\prod_{n=1}^{N_{\mathrm{dof}}}{\mathrm{d}\Delta\tilde{r}_{n}}\right)\prod_{p=1}^{N_{s}}{\delta\big{(}\Lambda_{p}\big{)}}}.\end{split}$
(40)

In the second step, we have used the definition of the $\Lambda_{p}$, Eq.
(24). The constraints by the Dirac deltas in this equations thus imply
$\Lambda_{p}=0$, where close to the transition, the $\Lambda_{p}$ are given by
Eq. (24).

We now discuss which values of the dofs, $\Delta\tilde{r}_{n}$, satisfy the
Dirac delta conditions in Eq. (40). In the case where all $B_{p}^{(1)}=0$,
using similar arguments as for the athermal limit, one can show that in the
$\Lambda_{p}$ with $p=1,\dots,N_{s}$, all terms of order
$\sim\Delta\tilde{r}_{m}\Delta\tilde{r}_{n}$ with $m<N_{s}$ or $n<N_{s}$ can
be neglected (appendix B.2). Thus, we can split the integrations in Eq. (40)
as follows:

$\Omega\sim\int{\left(\prod_{n=N_{s}}^{N_{\mathrm{dof}}}{\mathrm{d}\Delta\tilde{r}_{n}}\right)\delta(\Lambda_{N_{s}})\prod_{p=1}^{N_{s}-1}{\int{\mathrm{d}\Delta\tilde{r}_{p}\;\delta(\Lambda_{p})}}}.$
(41)

This is because to lowest order, each of the $\Delta\tilde{r}_{p}$ with
$p<N_{s}$ appears only once, namely as linear term in the corresponding
$\Lambda_{p}$. Thus, each of the $N_{s}-1$ inner integrals evaluates to a
constant factor of $1/s_{p}$. We are thus left with:

$\Omega\sim\Omega_{N_{s}}(\varepsilon,\gamma),$ (42)

where

$\Omega_{N_{s}}(\varepsilon,\gamma):=\int{\left(\prod_{n=N_{s}}^{N_{\mathrm{dof}}}{\mathrm{d}\Delta\tilde{r}_{n}}\right)\delta\big{(}\Lambda_{N_{s}}\big{)}}$
(43)

with $\Lambda_{N_{s}}$ given by Eq. (27).

To evaluate $\Omega_{N_{s}}(\varepsilon,\gamma)$, we first consider the
symmetric matrix $\tilde{M}_{N_{s}mn}$ restricted to
$m,n=N_{s},\dots,N_{\mathrm{dof}}$, which appears in $\Lambda_{N_{s}}$ (Eq.
(27)). We diagonalize this
$(N_{\mathrm{dof}}-N_{s}+1)\times(N_{\mathrm{dof}}-N_{s}+1)$ sub matrix:

$\tilde{M}_{N_{s}mn}=\sum_{k=1}^{N_{\mathrm{dof}}-N_{s}+1}{\mu_{k}v^{k}_{m}v^{k}_{n}},$
(44)

where $\mu_{k}$ are the eigenvalues and $v^{k}_{n}$ are the associated
orthonormal eigenvectors. We sort the eigenvalues in decreasing order,
$\mu_{1}\geq\dots\geq\mu_{N_{\mathrm{dof}}-N_{s}+1}\geq 0$. All eigenvalues
are non-negative, because the matrix $\tilde{M}_{N_{s}mn}$ is positive semi-
definite.

With

$\Delta
q_{k}:=\sum_{n=N_{s}}^{N_{\mathrm{dof}}}{v^{k}_{n}\Delta\tilde{r}_{n}},$ (45)

the equation for $\Lambda_{N_{s}}$, Eq. (27), transforms into:

$\Lambda_{N_{s}}=\frac{1}{2}\sum_{k=1}^{N_{\mathrm{dof}}-N_{s}+1}{\mu_{k}\Delta
q_{k}^{2}}+\tilde{w}_{N_{s}}\Big{[}\varepsilon+b_{\varepsilon}\gamma^{2}\Big{]},$
(46)

where we also used Eq. (29). We now substitute $\Lambda_{N_{s}}$ into
$\Omega_{N_{s}}$, Eq. (43):

$\begin{split}\Omega_{N_{s}}&=\int{\left(\prod_{k=1}^{N_{\mathrm{dof}}-N_{s}+1}{\mathrm{d}\Delta
q_{k}}\right)}\\\
&\times\delta\Bigg{(}\frac{1}{2}\sum_{k=1}^{N_{\mathrm{dof}}-N_{s}+1}{\mu_{k}\Delta
q_{k}^{2}}+\tilde{w}_{N_{s}}\Big{[}\varepsilon+b_{\varepsilon}\gamma^{2}\Big{]}\Bigg{)}.\end{split}$
(47)

Because the $v^{k}_{n}$ represent a rotation, the differential transforms
without prefactor.

Figure 1: Zero mode counting for four example systems, each both at the
transition point ($\varepsilon=0$) and under compression ($\varepsilon<0$). In
each case, we respectively indicate the number of spatial dimensions $D$, the
number of degrees of freedom $N_{\mathrm{dof}}$, the number of springs
$N_{s}$, the number of states of self-stress $N_{SSS}$, the number of
infinitesimal zero modes $N_{0}$, the number of 1st-order zero modes
$N_{\mathrm{1st}}$ (definition in subsection II.4), and the number of
mechanisms $N_{\mathrm{mech}}$. At the transition, $\varepsilon=0$, a finite
number $N^{\ast}_{\mathrm{1st}}$ of 1st-order zero modes exists. Upon
compression, $\varepsilon<0$, exactly $N^{\ast}_{\mathrm{1st}}-1$ of these
modes turn into mechanisms. The remaining mode ceases to be an infinitesimal
zero mode for $\varepsilon<0$, together with the disappearing SSS. While for
simplicity, we focus in this article on the case where there is no mechanism
at the transition point, we do not expect the existence of such mechanisms to
substantially change our results. An example is shown by the fourth system,
where a dangling spring creates additional mechanisms at the transition, which
however do not change the physics of the system as compared to the second
system.

###  IV.2 Eigenmodes of $\tilde{M}_{N_{s}mn}$

To evaluate the integral in Eq. (47), we first discuss how the matrix
$\tilde{M}_{N_{s}mn}$ relates to the different kinds of zero modes in the
system (see also subsection II.4).

At the transition point, there are $N_{\mathrm{dof}}-N_{s}+1$ linearly
independent infinitesimal zero modes. We consider here the zero modes
$\vec{z}_{k}$ with $k=1,\dots,N_{\mathrm{dof}}-N_{s}+1$ that have the
components $z_{k,m}:=\sum_{n=N_{s}}^{N_{\mathrm{dof}}}{V_{mn}v_{n}^{k}}$.
Insertion in Eq. (11) with Eq. (19) confirms that they are infinitesimal zero
modes. According to Eqs. (21) and (45) these modes contribute to
$\Delta\vec{r}=\vec{r}-\vec{r}^{\ast}$ with the respective amplitudes $\Delta
q_{k}$.

The $\Delta q_{k}$ modes that correspond to zero eigenvalues, $\mu_{k}=0$, of
the matrix $\tilde{M}_{N_{s}mn}$ do not affect the Dirac delta in Eq. (47).
There are three possibilities for the corresponding infinitesimal zero modes
$\vec{z}_{k}$:

  * •

$\vec{z}_{k}$ could be one of the $D$ _global translations_. In this case, the
integral in Eq. (47) leads to a factor of $e^{\varepsilon}$, which can again
be neglected in the end result.

  * •

$\vec{z}_{k}$ could belong to a _mechanism_ existing at the transition point,
but here we focus on systems without such mechanisms (subsection II.2).
However, even if we took such mechanisms into account, the integral over the
corresponding manifold would yield a constant to lowest order in $\varepsilon$
and $\gamma$. Thus, any possible variation of the contribution by a mechanism
with $\varepsilon$ and $\gamma$ would anyway only contribute higher-order
terms.

  * •

$\vec{z}_{k}$ belongs to a _$K$ th-order zero mode with $K>1$_. We believe
this is an exceptional case and do not consider it further in the main text,
but we partly discuss it in appendix F.

However, none of the infinitesimal zero modes $\vec{z}_{k}$ with $\mu_{k}=0$
can be a first-order zero mode 444 One only needs to show that for any given
$k$ with $\mu_{k}=0$, there exists a two-times differentiable function
$u\mapsto\vec{y}(u)$ with $\vec{y}(0)=\vec{r}\,^{\ast}$,
$\vec{y}^{\prime}(0)=\vec{z}_{k}$, and
$\mathrm{d}^{2}\ell_{i}(\vec{y}(u),\gamma)/\mathrm{d}u^{2}=0$. Indeed, such a
function is given by:
$\begin{split}y_{m}(u)&:=r_{m}^{\ast}+uV_{mn}v_{n}^{k}\\\
&\qquad\qquad-\frac{u^{2}}{2}\sum_{p=1}^{N_{s}-1}{\frac{V_{mp}}{s_{p}}\tilde{M}_{pn^{\prime}n^{\prime\prime}}v_{n^{\prime}}^{k}v_{n^{\prime\prime}}^{k}}\end{split}$
with $m=1,\dots,N_{\mathrm{dof}}$. .

Conversely, all of the infinitesimal zero modes $\vec{z}_{k}$ with $\mu_{k}>0$
_do_ affect the Dirac delta in Eq. (47). All of them are _first-order zero
modes_. This is because these modes affect $\Delta\tilde{\ell}_{N_{s}}$ to
second order 555For any two-times differentiable function $u\mapsto\vec{y}(u)$
with $\vec{y}(0)=\vec{r}\,^{\ast}$ and $\vec{y}^{\prime}(0)=\vec{z}_{k}$, it
is impossible to have
$\mathrm{d}^{2}\ell_{i}(\vec{y}(u),\gamma)/\mathrm{d}u^{2}=0$ for all
$i=1,\dots,N_{s}$, because
$\mathrm{d}^{2}\Delta\tilde{\ell}_{N_{s}}(\vec{y}(u),\gamma)/\mathrm{d}u^{2}=\mu_{k}>0$..
These first-order zero modes dominate the scaling behavior of $\Omega_{N_{s}}$
close to the transition, as we will see in the following subsection. We recall
from subsection II.4 that we denote the number of first-order zero modes at
the transition point by $N_{\mathrm{1st}}^{\ast}$.

In this subsection, we so far discussed the infinitesimal zero modes at the
transition point, $\varepsilon=0$. However, the situation changes as soon as
we move away from the transition point. We illustrate this this in Fig. 1,
where we show a few examples for $N_{\mathrm{1st}}^{\ast}=1,\dots,3$. As soon
as we apply isotropic compression, $\varepsilon<0$ (while formally keeping
$\gamma=0$), exactly $N_{\mathrm{1st}}^{\ast}-1$ of the first-order zero modes
turn into mechanisms. The remaining first-order zero mode ceases to be an
infinitesimal zero mode for $\varepsilon<0$. Indeed, the number of
infinitesimal zero modes needs to decrease by one, since the SSS that existed
at $\varepsilon=0$ disappears for $\varepsilon<0$ (compare Eq. (1)).

###  IV.3 Scaling of the the phase space volume

Based on our discussion in the previous section, integrating out the modes
$\Delta q_{k}$ for $\mu_{k}=0$ in Eq. (47) only leads to factors that are to
lowest order finite constants, and we obtain:

$\begin{split}\Omega_{N_{s}}&\sim\int{\left(\prod_{k=1}^{N_{\mathrm{1st}}^{\ast}}{\mathrm{d}\Delta
q_{k}}\right)}\\\
&\qquad\quad\times\delta\left(\frac{1}{2}\sum_{k=1}^{N_{\mathrm{1st}}^{\ast}}{\mu_{k}\Delta
q_{k}^{2}}+\tilde{w}_{N_{s}}\Big{[}\varepsilon+b_{\varepsilon}\gamma^{2}\Big{]}\right).\end{split}$
(48)

The Dirac delta ensures that the integral is over the surface
$\partial\mathcal{E}$ of a hyper-ellispoid $\mathcal{E}$ in a
$N_{\mathrm{1st}}^{\ast}$-dimensional space. The surface is defined by the
equation:

$\sum_{k=1}^{N_{\mathrm{1st}}^{\ast}}{\frac{\Delta
q_{k}^{2}}{\rho_{k}^{2}}}=1$ (49)

with half axes
$\rho_{k}:=(-2\tilde{w}_{N_{s}}[\varepsilon+b_{\varepsilon}\gamma^{2}]/\mu_{k})^{1/2}\sim[-(\varepsilon+b_{\varepsilon}\gamma^{2})]^{1/2}$.
In particular, for $\gamma=0$, the hyper-ellipsoid axes scale as
$\rho_{k}\sim\sqrt{-\varepsilon}$. Thus, consistent with the previous section
(subsection IV.2), as soon as compression is applied, $\varepsilon<0$, there
are $N_{\mathrm{1st}}^{\ast}-1$ mechanisms, which correspond to the surface
manifold of the hyper-ellipsoid. Meanwhile, movement perpendicular to the
ellipsoid surface is not allowed, corresponding to the first-order zero mode
that ceases to be an infinitesimal zero mode for $\varepsilon<0$.

Using standard transformation rules of the Dirac delta, we get:

$\Omega_{N_{s}}\sim\int_{\partial\mathcal{E}}{\mathrm{d}S\;\left|\frac{\partial\Lambda_{N_{s}}}{\partial\Delta\vec{q}}\right|^{-1}},$
(50)

where

$\begin{split}\left|\frac{\partial\Lambda_{N_{s}}}{\partial\Delta\vec{q}}\right|^{2}&=\sum_{k=1}^{N_{\mathrm{1st}}^{\ast}}{\left(\frac{\partial\Lambda_{N_{s}}}{\partial\Delta
q_{k}}\right)^{2}}\\\ &=\sum_{k=1}^{N_{\mathrm{1st}}^{\ast}}{\mu_{k}^{2}\Delta
q_{k}^{2}}.\end{split}$ (51)

With $\mu_{\mathrm{max}}:=\mu_{1}$ and
$\mu_{\mathrm{min}}:=\mu_{N_{\mathrm{1st}}^{\ast}}>0$, we thus have on the
ellipsoid surface, using Eq. (49):

$-\mu_{\mathrm{min}}\big{[}\varepsilon+b_{\varepsilon}\gamma^{2}\big{]}\leq\frac{1}{2\tilde{w}_{N_{s}}}\left|\frac{\partial\Lambda_{N_{s}}}{\partial\Delta\vec{q}}\right|^{2}\leq-\mu_{\mathrm{max}}\big{[}\varepsilon+b_{\varepsilon}\gamma^{2}\big{]}.$
(52)

Thus,
$\left|\partial\Lambda_{N_{s}}/\partial\Delta\vec{q}\right|\sim[-(\varepsilon+b_{\varepsilon}\gamma^{2})]^{1/2}$.
As a consequence, we obtain that the integral in Eq. (50) scales as the hyper-
surface area of the ellipsoid divided by
$[-(\varepsilon+b_{\varepsilon}\gamma^{2})]^{1/2}$:

$\Omega_{N_{s}}(\varepsilon,\gamma)\sim\Big{[}-(\varepsilon+b_{\varepsilon}\gamma^{2})\Big{]}^{(N_{\mathrm{1st}}^{\ast}-2)/2}.\\\
$ (53)

Hence, $\Omega\sim\Omega_{N_{s}}$ scales as a power law with the distance to
the transition point as long as $N_{\mathrm{1st}}^{\ast}>2$. In appendix F, we
show that the exponent in Eq. (53) can change if we allow for $K$th-order zero
modes with $K>1$ at the transition.

In appendix C, we discuss the more general case where the
$\tilde{B}_{p}^{(1)}$ do not vanish. The definition of $b_{\varepsilon}$
changes like in the athermal case, while the result Eq. (53) remains the same
up to an additional term $B_{N_{s}}^{(1)}\gamma$, which can be removed through
redefining $\gamma$.

###  IV.4 Free energy, tension, and shear modulus

The free energy of the system in the stiff-spring limit, $F_{S}$, is up to a
constant:

$\begin{split}F_{S}&=-\frac{1}{\beta}\log\Omega\\\
&=-\frac{N_{\mathrm{1st}}^{\ast}-2}{2\beta}\log{\Big{[}-(\varepsilon+b_{\varepsilon}\gamma^{2})\Big{]}},\end{split}$
(54)

where $\beta:=k_{B}T$ with $k_{B}$ being the Boltzmann constant. The index $S$
stresses that $F_{S}$ results entirely from entropic effects.

We then obtain for the tension $t_{S}=(\partial
F_{S}/\partial\varepsilon)/D\mathcal{V}$ and the shear modulus
$G_{S}=(\partial^{2}F_{S}/\partial\gamma^{2})/\mathcal{V}$, to lowest order:

$\displaystyle t_{S}$
$\displaystyle=-\frac{\kappa_{S}T}{\varepsilon+b_{\varepsilon}\gamma^{2}}$
(55) $\displaystyle G_{S}$
$\displaystyle=-2Db_{\varepsilon}\kappa_{S}T\frac{\varepsilon-
b_{\varepsilon}\gamma^{2}}{(\varepsilon+b_{\varepsilon}\gamma^{2})^{2}}$ (56)
$\displaystyle=2Db_{\varepsilon}t\left(1+\frac{2b_{\varepsilon}\gamma^{2}}{|\varepsilon+b_{\varepsilon}\gamma^{2}|}\right)$
(57)

with

$\kappa_{S}:=\frac{k_{B}(N_{\mathrm{1st}}^{\ast}-2)}{2D\mathcal{V}^{\ast}}.$
(58)

Note that only the regime with $\varepsilon+b_{\varepsilon}\gamma^{2}<0$ is
accessible here, which corresponds to the athermal floppy regime.

##  V General case

###  V.1 Partition sum

To study how entropic and energetic elasticity interact, we evaluate the
partition sum of the system

$Z=\int{\left(\prod_{n=1}^{N_{\mathrm{dof}}}{\mathrm{d}R_{n}}\right)e^{-\beta
E}}.$ (59)

Transforming the dofs $R_{n}$ to the dimensionless eigenmodes of the
compatibility matrix, $\Delta\tilde{r}_{n}$, we have:

$Z\sim\int{\left(\prod_{k=1}^{N_{\mathrm{dof}}}{\mathrm{d}\Delta\tilde{r}_{n}}\right)e^{-\beta
E}},$ (60)

where we again ignore prefactors proportional to $e^{\varepsilon}$.

We evaluate the partition sum for small $\varepsilon$, $\gamma$, and $T$.
Because $\varepsilon\ll 1$, we can use the energy expression from Eq. (25):

$E=\frac{E_{0}}{2}\sum_{p=1}^{N_{s}}{\Lambda_{p}^{2}}.$ (61)

Here we discuss the case where all $\tilde{B}_{p}^{(1)}=0$, while the more
general case of arbitrary $\tilde{B}_{p}^{(1)}$ is discussed in appendix C.

We focus on small temperatures, $k_{B}T\lesssim
E_{0}[\varepsilon+b_{\varepsilon}\gamma^{2}]^{2}$, such that states with
higher energies $E$ receive an exponentially smaller Boltzmann weight. Using
similar arguments as in the athermal and stiff-spring limits, we find again
that in all $\Lambda_{q}$ with $q\leq N_{s}$ all the $\Delta\tilde{r}_{p}$
with $p<N_{s}$ can be neglected in the terms
$\sim\Delta\tilde{r}_{m}\Delta\tilde{r}_{n}$ (details in appendix B.3). Thus,
to lowest order each of the $\Delta\tilde{r}_{p}$ with $p<N_{s}$ appears only
once, namely in the linear term of the respective $\Lambda_{p}$. As a
consequence, the integrals in the partition sum in Eq. (60) can be rearranged
as follows:

$Z\sim\int{\left(\prod_{n=N_{s}}^{N_{\mathrm{dof}}}{\mathrm{d}\Delta\tilde{r}_{n}}\right)\;e^{-\frac{\beta
E_{0}\Lambda_{N_{s}}^{2}}{2}}}\prod_{p=1}^{N_{s}-1}{\int{\mathrm{d}\Delta\tilde{r}_{p}}\;e^{-\frac{\beta
E_{0}\Lambda_{p}^{2}}{2}}}.$ (62)

The innermost first $N_{s}-1$ integrals are Gaussian integrals, each of which
evaluates to the constant $Z_{p}=\sqrt{2\pi/\beta E_{0}s_{p}^{2}}$.

In a system where each of the nodes also has a mass $m$, there is an
additional constant prefactor in the partition sum in Eq. (59), which arises
from the kinetic energy of the nodes [32]. This prefactor is given by
$\lambda^{-DN_{\mathrm{node}}}$, where $\lambda$ is the thermal wavelength
$\lambda:=\sqrt{\beta h^{2}/2\pi m}$ with Planck constant $h$. When combining
each of the $Z_{p}$ for $p<N_{s}$ with one of the thermal wavelength factors,
one obtains the partition sum of a harmonic oscillator,
$Z_{p}/\lambda=2\pi/\beta h\omega_{p}$, with eigen frequency
$\omega_{p}=s_{p}\sqrt{E_{0}/m}$. This means that each of the nonzero eigen
modes of the Hessian at the transition point just creates a harmonic
oscillator, like for isostatic and over-constrained systems. However, for
under-constrained systems, we additionally have the last integral, which stems
from the SSS that forms at the transition:

$Z_{N_{s}}:=\int{\left(\prod_{n=N_{s}}^{N_{\mathrm{dof}}}{\mathrm{d}\Delta\tilde{r}_{n}}\right)\;e^{-\frac{\beta
E_{0}\Lambda_{N_{s}}^{2}}{2}}}.$ (63)

Neglecting the factors $Z_{p}$ with $p<N_{s}$, which solely depend on
temperature, but not on strain, we have $Z\sim Z_{N_{s}}$.

###  V.2 Saddle point approximation

To simplify the remaining integral $Z_{N_{s}}$ in Eq. (63), we introduce in
its integrand an additional integral over some variable $\varepsilon_{S}$ that
trivially evaluates to one:

$\begin{split}Z_{N_{s}}&=\int{\left(\prod_{n=N_{s}}^{N_{\mathrm{dof}}}{\mathrm{d}\Delta\tilde{r}_{n}}\right)\;e^{-\frac{\beta
E_{0}\Lambda_{N_{s}}^{2}}{2}}}\\\
&\qquad\times\tilde{w}_{N_{s}}\int_{-\infty}^{\infty}{\mathrm{d}\varepsilon_{S}\;\delta\Big{(}\Lambda_{N_{s}}+\tilde{w}_{N_{s}}[\varepsilon_{S}-\varepsilon]\Big{)}}.\end{split}$
(64)

After exchanging the order of integration and substituting $\Omega_{N_{s}}$,
defined in Eq. (43), we have:

$\begin{split}Z_{N_{s}}&=\tilde{w}_{N_{s}}\int_{-\infty}^{\infty}{\mathrm{d}\varepsilon_{S}\;\Omega_{N_{s}}(\varepsilon_{S},\gamma)e^{-\beta
E_{N_{s}}(\varepsilon-\varepsilon_{S})}},\end{split}$ (65)

with:

$\begin{split}E_{N_{s}}(\varepsilon-\varepsilon_{S})&:=\frac{E_{0}\tilde{w}_{N_{s}}^{2}}{2}\big{(}\varepsilon-\varepsilon_{S}\big{)}^{2}\\\
&=\frac{D\mathcal{V}^{\ast}\kappa_{E}}{2}\big{(}\varepsilon-\varepsilon_{S}\big{)}^{2}.\end{split}$
(66)

Here, we used the definition for $\kappa_{E}$ in Eq. (36). From Eq. (53) in
the previous section, we have for
$\varepsilon_{S}+b_{\varepsilon}\gamma^{2}\leq 0$:

$\Omega_{N_{s}}(\varepsilon_{S},\gamma)\sim\Big{[}-(\varepsilon_{S}+b_{\varepsilon}\gamma^{2})\Big{]}^{(N_{\mathrm{1st}}^{\ast}-2)/2},$
(67)

while for $\varepsilon_{S}+b_{\varepsilon}\gamma^{2}>0$, no configurations are
possible, $\Omega_{N_{s}}(\varepsilon_{S},\gamma)=0$. These transformations
identify $\varepsilon_{S}$ as the entropic strain discussed in the companion
paper [26].

We use a saddle point approximation to simplify the integral in Eq. (65):

$Z_{N_{s}}=\tilde{w}_{N_{s}}\int_{-\infty}^{\infty}{\mathrm{d}\varepsilon_{S}\;e^{-\beta\bar{F}}},$
(68)

where $\bar{F}$ is the free energy of the system for imposed
$\varepsilon_{S}$:

$\begin{split}\bar{F}(\varepsilon_{S};\varepsilon,\gamma)&:=E_{N_{s}}(\varepsilon-\varepsilon_{S})-\frac{1}{\beta}\log{\Omega_{N_{s}}(\varepsilon_{S},\gamma)}\\\
&=D\mathcal{V}^{\ast}\bigg{(}\frac{\kappa_{E}}{2}\big{(}\varepsilon-\varepsilon_{S}\big{)}^{2}\\\
&\qquad\qquad-\kappa_{S}T\log{\Big{[}-(\varepsilon_{S}+b_{\varepsilon}\gamma^{2})\Big{]}}\bigg{)}\end{split}$
(69)

On the second line, we used the definition of $\kappa_{S}$, Eq. (58), and
ignored a constant offset.

For small temperatures, we can apply the saddle point approximation, where we
Taylor expand $\bar{F}(\varepsilon_{S})$ to second order around its minimum at
$\varepsilon_{S}=\hat{\varepsilon}_{S}$, which results in:

$Z_{N_{s}}\sim\int_{-\infty}^{\infty}{\mathrm{d}\varepsilon_{S}\;e^{-\beta\bar{F}}}=e^{-\beta\bar{F}(\hat{\varepsilon}_{S})}\sqrt{\frac{2\pi}{\beta\bar{F}^{\prime\prime}(\hat{\varepsilon}_{S})}},$
(70)

where $\bar{F}^{\prime\prime}:=\partial^{2}F/\partial\varepsilon_{S}^{2}$:

$\begin{split}\bar{F}^{\prime\prime}(\hat{\varepsilon}_{S})&=D\mathcal{V}^{\ast}\Bigg{(}\kappa_{E}+\frac{\kappa_{S}T}{[\hat{\varepsilon}_{S}+b_{\varepsilon}\gamma^{2}]^{2}}\Bigg{)}.\end{split}$
(71)

To find the minimum $\hat{\varepsilon}_{S}$, we transform
$\partial\bar{F}/\partial\varepsilon_{S}=0$, and obtain:

$\hat{\varepsilon}_{S}+b_{\varepsilon}\gamma^{2}=\frac{\varepsilon+b_{\varepsilon}\gamma^{2}}{2}-\sqrt{1+\theta}\left|\frac{\varepsilon+b_{\varepsilon}\gamma^{2}}{2}\right|,$
(72)

with

$\begin{split}\theta&:=\frac{4\kappa_{S}T}{\kappa_{E}(\varepsilon+b_{\varepsilon}\gamma^{2})^{2}}\\\
&=\frac{2(N_{\mathrm{1st}}^{\ast}-2)k_{B}T}{E_{0}\tilde{w}_{N_{s}}^{2}(\varepsilon+b_{\varepsilon}\gamma^{2})^{2}}.\end{split}$
(73)

Note that for $\gamma=0$, the equation
$\partial\bar{F}(\hat{\varepsilon}_{S})/\partial\varepsilon_{S}=0$, defining
the minimum $\hat{\varepsilon}_{S}$ of
$\bar{F}(\varepsilon_{S})=E(\varepsilon-\varepsilon_{S})+F_{S}(\varepsilon_{S})$,
corresponds to equating energetic and entropic tension, like in our intuitive
explanation in the companion paper [26].

###  V.3 Free energy, tension, and shear modulus

Up to terms that depend only on temperature and up to higher-order terms in
the strain variables, the free energy of the system, $F=-(\log{Z})/\beta$, is:

$F=\bar{F}(\hat{\varepsilon}_{S})+\frac{1}{2\beta}\log{\bar{F}^{\prime\prime}(\hat{\varepsilon}_{S})}.$
(74)

After some transformations, we find for the isotropic tension $t=(\partial
F/\partial\varepsilon)/D\mathcal{V}$ and the shear modulus
$G=(\partial^{2}F/\partial\gamma^{2})/\mathcal{V}$, to lowest order in
$\varepsilon$ and $\gamma$ (appendix G):

$t=\kappa_{E}\big{(}\varepsilon-\hat{\varepsilon}_{S}\big{)}\left[1+\frac{1}{2(N_{\mathrm{1st}}^{\ast}-2)}\;\frac{\theta}{1+\theta}\right]$
(75)

and

$\begin{split}G&=2Db_{\varepsilon}\kappa_{E}\big{(}\varepsilon-\hat{\varepsilon}_{S}\big{)}\Bigg{[}1+\frac{2b_{\varepsilon}\gamma^{2}}{|\varepsilon+b_{\varepsilon}\gamma^{2}|\sqrt{1+\theta}}\\\
&\quad+\frac{1}{2(N_{\mathrm{1st}}^{\ast}-2)}\;\frac{\theta}{1+\theta}\Bigg{(}1+\frac{2b_{\varepsilon}\gamma^{2}}{|\varepsilon+b_{\varepsilon}\gamma^{2}|\sqrt{1+\theta}}\\\
&\qquad\qquad\qquad\qquad\qquad\;\;-\frac{4b_{\varepsilon}\gamma^{2}}{(\varepsilon+b_{\varepsilon}\gamma^{2})(1+\theta)}\Bigg{)}\Bigg{]}.\end{split}$
(76)

Note that the last terms in both $t$ and $G$ scale as
$\sim(N_{\mathrm{1st}}^{\ast}-2)^{-1}$. These terms stem from the term
$\sim\log{\bar{F}^{\prime\prime}}$ in the free energy. For
$N_{\mathrm{1st}}^{\ast}\gg 1$, these contributions can be neglected, and we
have:

$\displaystyle t$
$\displaystyle=\kappa_{E}\big{(}\varepsilon-\hat{\varepsilon}_{S}\big{)}$ (77)
$\displaystyle G$
$\displaystyle=2Db_{\varepsilon}\kappa_{E}\big{(}\varepsilon-\hat{\varepsilon}_{S}\big{)}\Bigg{[}1+\frac{2b_{\varepsilon}\gamma^{2}}{|\varepsilon+b_{\varepsilon}\gamma^{2}|\sqrt{1+\theta}}\Bigg{]}$
(78)
$\displaystyle=2Db_{\varepsilon}t\Bigg{[}1+\frac{2b_{\varepsilon}\gamma^{2}}{|\varepsilon+b_{\varepsilon}\gamma^{2}|\sqrt{1+\theta}}\Bigg{]}.$
(79)

In these equations, the value of $\hat{\varepsilon}_{S}$ can be obtained using
Eq. (72).

###  V.4 Limits

####  V.4.1 Small temperature, $\theta\ll 1$

For $\theta\ll 1$, Eq. (72) becomes:

$\begin{split}\hat{\varepsilon}_{S}+b_{\varepsilon}\gamma^{2}&=\frac{\varepsilon+b_{\varepsilon}\gamma^{2}}{2}-\left(1+\frac{\theta}{2}\right)\left|\frac{\varepsilon+b_{\varepsilon}\gamma^{2}}{2}\right|\\\
&=(\varepsilon+b_{\varepsilon}\gamma^{2})\times\begin{cases}1+\theta/4&\text{for
$\varepsilon+b_{\varepsilon}\gamma^{2}<0$,}\\\ -\theta/4&\text{for
$\varepsilon+b_{\varepsilon}\gamma^{2}>0$.}\end{cases}\end{split}$ (80)

Thus, to lowest order in $\theta$:

$\varepsilon-\hat{\varepsilon}_{S}=(\varepsilon+b_{\varepsilon}\gamma^{2})\times\begin{cases}-\theta/4&\text{for
$\varepsilon+b_{\varepsilon}\gamma^{2}<0$,}\\\ 1&\text{for
$\varepsilon+b_{\varepsilon}\gamma^{2}>0$,}\end{cases}$ (81)

and we recover the purely entropic and energetic tensions and shear moduli,
respectively (compare Eqs. (33), (35), (55), and (57)):

$\displaystyle t$ $\displaystyle=\begin{cases}t_{S}&\text{for
$\varepsilon+b_{\varepsilon}\gamma^{2}<0$,}\\\ t_{E}&\text{for
$\varepsilon+b_{\varepsilon}\gamma^{2}>0$,}\end{cases}$ (82) $\displaystyle G$
$\displaystyle=\begin{cases}G_{S}&\text{for
$\varepsilon+b_{\varepsilon}\gamma^{2}<0$,}\\\ G_{E}&\text{for
$\varepsilon+b_{\varepsilon}\gamma^{2}>0$.}\end{cases}$ (83)

Note that for $\varepsilon+b_{\varepsilon}\gamma^{2}$ close to zero, we are
not in the limit $\theta\ll 1$ any more (see Eq. (73)).

####  V.4.2 Small strain, $\theta\gg 1$

For $\theta\gg 1$, Eq. (72) becomes:

$\hat{\varepsilon}_{S}+b_{\varepsilon}\gamma^{2}=-\sqrt{\frac{\kappa_{S}T}{\kappa_{E}}}+\frac{\varepsilon+b_{\varepsilon}\gamma^{2}}{2}$
(84)

And thus:

$\varepsilon-\hat{\varepsilon}_{S}=\sqrt{\frac{\kappa_{S}T}{\kappa_{E}}}+\frac{\varepsilon+b_{\varepsilon}\gamma^{2}}{2}.$
(85)

We then find:

$t=\frac{2N_{\mathrm{1st}}^{\ast}-3}{2N_{\mathrm{1st}}^{\ast}-4}\left(\sqrt{\kappa_{E}\kappa_{S}T}+\kappa_{E}\frac{\varepsilon+b_{\varepsilon}\gamma^{2}}{2}\right),$
(86)

and for the shear modulus:

$\begin{split}G&=Db_{\varepsilon}\bigg{[}2\sqrt{\kappa_{E}\kappa_{S}T}+\kappa_{E}(\varepsilon+b_{\varepsilon}\gamma^{2})\bigg{]}\\\
&\qquad\times\Bigg{[}\frac{2N_{\mathrm{1st}}^{\ast}-3}{2N_{\mathrm{1st}}^{\ast}-4}\Bigg{(}1+b_{\varepsilon}\gamma^{2}\sqrt{\frac{\kappa_{E}}{\kappa_{S}T}}\Bigg{)}\\\
&\qquad\qquad\qquad\qquad-\frac{b_{\varepsilon}\gamma^{2}}{2N_{\mathrm{1st}}^{\ast}-4}\;\frac{\kappa_{E}(\varepsilon+b_{\varepsilon}\gamma^{2})}{\kappa_{S}T}\Bigg{)}\Bigg{]}.\end{split}$
(87)

For $N_{\mathrm{1st}}^{\ast}\gg 1$ both expressions simplify to:

$\displaystyle t$
$\displaystyle=\sqrt{\kappa_{E}\kappa_{S}T}+\kappa_{E}\frac{\varepsilon+b_{\varepsilon}\gamma^{2}}{2}$
(88) $\displaystyle G$
$\displaystyle=2Db_{\varepsilon}t\Bigg{(}1+b_{\varepsilon}\gamma^{2}\sqrt{\frac{\kappa_{E}}{\kappa_{S}T}}\Bigg{)}.$
(89)

Note that we assumed in our derivation for the general case that
$k_{B}T\lesssim E_{0}[\varepsilon+b_{\varepsilon}\gamma^{2}]^{2}$. Using the
definitions of $\theta$ in Eq. (73), of $\kappa_{E}$ in Eq. (36), and of
$\kappa_{S}$ in Eq. (58), this condition corresponds to $\theta\lesssim
2(N_{\mathrm{1st}}^{\ast}-2)/\tilde{w}_{N_{s}}^{2}$. Thus, this limit
corresponds to $1\ll\theta\lesssim
2(N_{\mathrm{1st}}^{\ast}-2)/\tilde{w}_{N_{s}}^{2}$, which can be fulfilled
for large systems. Note that in our numerical results for randomly-cut
triangular networks, we also found that our prediction was accurate for finite
temperatures even for vanishing strain,
$\varepsilon+b_{\varepsilon}\gamma^{2}\rightarrow 0$ [26].

##  VI Discussion

In this article, we developed a theory for thermal, under-constrained systems
with fixed connectivity. We provide analytical expressions for the partition
sum $Z$, free energy $F$, and elastic material properties such as tension $t$
and shear modulus $G$ close to the athermal rigidity transition depending on
isotropic strain $\varepsilon$, shear strain $\gamma$, and temperature $T$.
The only three parameters appearing in our theory are energetic rigidity
$\kappa_{E}$, entropic rigidity $\kappa_{S}$, and an interaction parameter
between isotropic and shear strain, which we denote by $b_{\varepsilon}$.
These parameters are phenomenological coefficients that depend on the
microscopic structure of the system.

We first discussed the limit of zero temperature (athermal limit), which is
dominated by energetic rigidity, and the limit of infinitely stiff springs,
which is dominated by entropic rigidity. In the athermal limit, we
analytically derived the elastic properties before [3, 4], where the system is
floppy for $\varepsilon+b_{\varepsilon}\gamma^{2}\leq 0$ and rigid for
$\varepsilon+b_{\varepsilon}\gamma^{2}>0$. For the stiff-spring limit, we
showed that the accessible phase space volume corresponds to the surface of a
hyper-ellipsoid, whose dimension is given by the number of first-order zero
modes of the system, $N_{\mathrm{1st}}^{\ast}$ (subsection II.4). This
ultimately leads to a free energy that scales logarithmically with the
negative combined strain $-(\varepsilon+b_{\varepsilon}\gamma^{2})$. In both
limits, the mechanical properties are inherently linked to the properties of
the SSS that is created at the transition. Finally, we derive the partition
sum for the general case, where energetic and entropic rigidity are coupled
only through the SSS.

In the regime with $\epsilon<0$, previous numerical work has reported a shear
modulus scaling of $G\sim T^{0.8}$ [22, 33], while still reporting a linear
scaling of the isotropic tension, $t\sim T^{1}$ [33]. This partly contradicts
our results here, where in the regime with $\varepsilon<0$ and small
temperatures ($\theta\ll 1$), we obtain the scaling relations $G\sim t\sim
T^{1}$. We do not know the precise reason for the deviation of the shear
modulus scaling, but we see the following possibilities. First, all results
obtained here are based on the Taylor expansion in Eq. (18), which holds for
small deviations of the node positions from the transition point
configuration, and for small shear strain $\gamma$. A deviation of shear
modulus scaling but not of tension scaling could thus be explained by too high
values of the shear strain $\gamma$ for the configurations studied in Ref.
[33], such that terms of higher than quadratic order in $\gamma$ become
important in Eq. (18).

Second, Refs. [22, 33] seem to treat shear strain differently than we do here.
While we shear stabilize our configurations before the thermal simulations
[26], the authors of Ref. [22] do not seem to have done this, which can also
affect the apparent shear modulus scaling [4].

Third, in our simulations [26], we found that the numerically observed value
of the shear modulus $G$ can be more sensitive to insufficient thermal
equilibration than the system energy or the components of the stress tensor
are. Moreover, for the same number of simulations steps, insufficient
equilibration would more strongly affect the smaller temperatures. This might
be consistent with the results in Ref. [22], where a scaling of $G\sim
T^{0.8}$ was observed for small temperatures, while $G\sim T^{1}$ was observed
for larger temperatures.

Finally, Refs. [22, 33] suggested that the shear modulus scaling exponent
different from one was due to purely entropic effects in combination with the
disordered nature of the networks. However, in our approach, we fully take the
disordered network structure into account without any kind of mean-field-like
assumption, and we obtain an exponent of exactly 1. Indeed, in the infinitely-
stiff-spring limit, where elasticity is of purely entropic nature, an exponent
of exactly 1 is expected purely on dimensional grounds since the temperature
$T$ provides the only energy scale.

While not explicitly studied here, our results have immediate consequences
with respect to the scaling with network connectivity $z$. For $\gamma=0$,
earlier work by Zhang and Mao using EMT indicated three scaling regimes for
the shear modulus $G$ of randomly-cut triangular networks [24]. Remarkably,
these scaling regimes are very similar to what we report in subsection V.4,
where up to a prefactor, isotropic strain $\varepsilon$ is replaced by the
relative connectivity $\Delta z=z-z_{c}$, with $z_{c}$ being the connectivity
at isostaticity [24]. To understand how this could be, we discuss the $\Delta
z$ scaling that results from our findings by separately discussing the $\Delta
z$ scaling for $\kappa_{S}$ and $\kappa_{E}$. Eq. (58) predicts that
$\kappa_{S}$ essentially only depends on the number of first-order zero modes
$N_{\mathrm{1st}}^{\ast}$. Thus, if for large systems
$N_{\mathrm{1st}}^{\ast}$ scales to dominating order linearly with $\Delta
z=-4(N_{\mathrm{dof}}-N_{s})/N_{\mathrm{dof}}$, then this implies
$\kappa_{S}\sim|\Delta z|^{1}$. Meanwhile, how energetic elasticity
$\kappa_{E}$ scales with connectivity $\Delta z$ depends on the class of the
network [4]. For instance, for $|\Delta z|\ll 1$ our results for the
$a_{\ell}$ scaling in Refs. [3, 4] imply that $\kappa_{E}\sim|\Delta z|^{2}$
for Delaunay networks, $\kappa_{E}\sim|\Delta z|^{1}$ for randomly-cut,
packing-derived networks, and $\kappa_{E}\sim|\Delta z|^{0}$ for phantom
triangular networks. Because of Eq. (77), these scaling relations of
$\kappa_{S}$ and $\kappa_{E}$ with $\Delta z$ determine the scaling of tension
$t$. However, for the scaling of the shear modulus $G$, Eq. (79), the scalings
of $b_{\varepsilon}\kappa_{S}$ and $b_{\varepsilon}\kappa_{E}$ are the
relevant ones, and thus the scaling of $b_{\varepsilon}\sim b$ needs to be
taken into account. Based on earlier results for the scaling of $b$ [3, 4], we
have $b_{\varepsilon}\kappa_{S}\sim|\Delta z|^{-1}$ and
$b_{\varepsilon}\kappa_{E}\sim|\Delta z|^{0}$ for Delaunay networks,
$b_{\varepsilon}\kappa_{S}\sim|\Delta z|^{0}$ and
$b_{\varepsilon}\kappa_{E}\sim|\Delta z|^{0}$ for randomly-cut, packing-
derived networks, and $b_{\varepsilon}\kappa_{S}\sim|\Delta z|^{0.5}$ and
$b_{\varepsilon}\kappa_{E}\sim|\Delta z|^{-0.5}$ for phantom triangular
networks. Taken together, this suggests that the scaling of the elastic
properties with $\Delta z$ is expected to depend on the class of network
studied. For instance, the results by Zhang and Mao [24] can be explained from
the perspective of our results if the randomly-cut triangular networks have
$b_{\varepsilon}\kappa_{S}\sim|\Delta z|^{0}$ and
$b_{\varepsilon}\kappa_{E}\sim|\Delta z|^{0}$ like randomly-cut, packing-
derived networks. In this case, the scaling is entirely determined by the
scaling of the strain of the network state after network initialization,
$\varepsilon_{\mathrm{ini}}$, with connectivity $\Delta z$, which is generally
linear to lowest order: $\varepsilon_{\mathrm{ini}}\sim\Delta z$ (remember
that we define isotropic strain here with respect to the athermal transition
point) [3, 4]. This would be a simple explanation for why the shear modulus of
Zhang and Mao has the same form as our results when replacing $\varepsilon$ by
$\Delta z$. We leave a more thorough numerical test of all these predicitons
for the scaling with $\Delta z$, and a study of their consequences, for future
work.

We introduced under-constrained system as systems that have
$N_{\mathrm{dof}}>N_{s}$, and our results apply to a broad subset of such
systems. In particular, we focused on systems that fulfill the following
criteria:

  * •

_There is no SSS in the floppy regime._ Lifting the restriction on the number
of SSS in the floppy regime probably does not change our results.

  * •

_Only a single SSS forms at the transition._ In particular in disordered
networks, situations where more than a single SSS forms at the transition are
likely very rare. Also, we do not expect the number of SSS that form at the
transition to substantially affect our results.

  * •

_There are no $K$th-order zero modes at the transition,
$N_{K\mathrm{th}}^{\ast}=0$, for $K>1$._ We expect the occurrence of such
modes to be rare. Moreover, some preliminary work indicates that allowing for
$N_{K\mathrm{th}}^{\ast}>0$ for $K>1$ will mostly affect the prefactor for the
entropic rigidity $\kappa_{S}$, but may not change our overall results
(appendix F).

  * •

_There are more than two first-order zero modes at the transition,
$N_{\mathrm{1st}}^{\ast}>2$._ The marginal cases of
$N_{\mathrm{1st}}^{\ast}=1$ corresponds to the situation where there is no
floppy-rigid transition (compare Fig. 1), since no mechanisms would form in
the floppy regime. It will be interesting to study the marginal case
$N_{\mathrm{1st}}^{\ast}=2$ in future work.

  * •

_For our discussion inappendix B, we needed the additional assumption in Eq.
(94)._ Unfortunately, we currently needed to rely on this assumption. However,
we expect that it should be possible to obtain our result in appendix B quite
generally even without this assumption.

  * •

_There are no mechanisms at the transition._ The reason why we need this
assumption is that with a degeneracy $\vec{r}\,^{\ast}$ of the transition
point configurations, in principle, the Taylor expansion in Eq. (18) could
depend on $\vec{r}\,^{\ast}$. However, while properties like the singular
values, $s_{p}$ with $p<N_{s}$, of the compatibility matrix $C_{in}$ may
depend on $\vec{r}\,^{\ast}$ [34, 6], our results should be unaffected as long
as the SSS (more precisely $\tilde{w}_{N_{s}}$) and its behavior with respect
to shear (more precisely $b_{\varepsilon}$) are independent of the choice of
$\vec{r}\,^{\ast}$. While this seems to be true intuitively, preliminary
results suggest that stronger assumptions about the system than those listed
in subsection II.1 are required to rigorously show it. The tools required for
such a proof are outside the scope of the ideas presented here, and their
development is left for future work.

Thus, most of these criteria do not essentially restrict our results here
and/or correspond to marginal cases. Yet, including these cases provides
interesting avenues for future work.

###### Acknowledgements.

We thank Chris Santangelo, Jen Schwarz, and Manu Mannattil for fruitful
discussions. We thank the Centre Interdisciplinaire de Nanoscience de
Marseille (CINaM) for providing office space. The project leading to this
publication has received funding from France 2030, the French Government
program managed by the French National Research Agency (ANR-16-CONV-0001), and
from the Excellence Initiative of Aix-Marseille University - A*MIDEX.

##  Appendix A Implementation of periodic boundary conditions and shear

As mentioned in subsection II.1, our framework is valid both for simple and
pure shear. Here, we provide examples for the implementation of either, by
providing explicit expressions for the spring length functions
$L_{i}(\vec{R},\mathcal{V},\gamma)$ for a system of ordinary linear springs.

To implement simple shear in a $D=2$-dimensional system, we compute the length
of any spring $i$ that connects node $a_{i}$ to node $b_{i}$, we first define
the components of the spring length vector $\bm{L}_{i}=(L_{i,x},L_{i,y})$ as
follows [26]:

$\displaystyle L_{i,x}(\vec{R},\mathcal{V},\gamma)$
$\displaystyle=R_{b_{i},x}-R_{a_{i},x}+(q_{i,x}+\gamma
q_{i,y})\sqrt{\mathcal{V}}$ (90) $\displaystyle
L_{i,y}(\vec{R},\mathcal{V},\gamma)$
$\displaystyle=R_{b_{i},y}-R_{a_{i},y}+q_{i,y}\sqrt{\mathcal{V}},$

where $\bm{R}_{a}=(R_{a,x},R_{a,y})$ is the position of some node $a$. The
spring length is then given by:
$L_{i}(\vec{R},\mathcal{V},\gamma)=|\bm{L}_{i}(\vec{R},\mathcal{V},\gamma)|$.

In Eq. (90), the symbol $\bm{q}_{i}=(q_{i,x},q_{i,y})\in\mathbb{Z}^{D}$
denotes the _periodicity vector_ of spring $i$ [35, 13]. A value of
$q_{i,x}=0$ indicates that the spring does not cross the vertical periodic
boundary of the system, while $q_{i,x}\neq 0$ indicates that the spring _does_
cross the vertical periodic boundary. In particular, a positive value
$q_{i,x}>0$ indicates that when following the spring going from
$\bm{R}_{a_{i}}$ to $\bm{R}_{b_{i}}$ one leaves via the right vertical
boundary while entering via the left boundary exactly $q_{i,x}$ times.
Analogously, a negative value $q_{i,x}<0$ indicates that following the spring
from $\bm{R}_{a_{i}}$ to $\bm{R}_{b_{i}}$ one leaves the via _left_ boundary
while entering via the _right_ boundary exactly $-q_{i,x}$ times.

A similar way of implementing simple shear is given by:

$\displaystyle L_{i,x}(\vec{R},\mathcal{V},\gamma)$
$\displaystyle=R_{b_{i},x}-R_{a_{i},x}+q_{i,x}\sqrt{\mathcal{V}}$ (91)
$\displaystyle\qquad+\gamma\Big{(}R_{b_{i},y}-R_{a_{i},y}+q_{i,y}\sqrt{\mathcal{V}}\Big{)}$
$\displaystyle L_{i,y}(\vec{R},\mathcal{V},\gamma)$
$\displaystyle=R_{b_{i},y}-R_{a_{i},y}+q_{i,y}\sqrt{\mathcal{V}}.$

This corresponds to Lees-Edward boundary conditions [36].

To implement pure shear in a $D=3$-dimensional system, we define the
components of the length vector $\bm{L}_{i}=(L_{i,x},L_{i,y},L_{i,z})$ of a
spring $i$ that that connects node $a_{i}$ to node $b_{i}$ as follows:

$\displaystyle L_{i,x}(\vec{R},\mathcal{V},\gamma)$
$\displaystyle=R_{b_{i},x}-R_{a_{i},x}+q_{i,x}e^{\gamma}\sqrt[3]{\mathcal{V}}$
(92) $\displaystyle L_{i,y}(\vec{R},\mathcal{V},\gamma)$
$\displaystyle=R_{b_{i},y}-R_{a_{i},y}+q_{i,y}e^{-\gamma}\sqrt[3]{\mathcal{V}}$
$\displaystyle L_{i,z}(\vec{R},\mathcal{V},\gamma)$
$\displaystyle=R_{b_{i},z}-R_{a_{i},z}+q_{i,z}\sqrt[3]{\mathcal{V}}.$

Again, the spring length is then given by:
$L_{i}(\vec{R},\mathcal{V},\gamma)=|\bm{L}_{i}(\vec{R},\mathcal{V},\gamma)|$.
Different pure shear implementations can vary the $e^{\gamma}$ factors in the
periodicity terms, such that the product of all these terms remains one.

Note that in all example cases discussed here, the functions
$L_{i}(\vec{R},\mathcal{V},\gamma)$ are analytic as long as $L_{i}>0$ and
$\mathcal{V}>0$. Moreover, the relation Eq. (2) holds for all spring vector
definitions, Eqs. (90)–(92) with $d_{i}=1$. Finally, a similar approach as in
Eqs. (90)–(92) can be used to implement simple or pure shear in vertex models
instead of simple spring networks [35, 13].

##  Appendix B Orders of magnitude of the $\Lambda_{p}$ and
$\Delta\tilde{r}_{n}$

###  B.1 Athermal limit

We discuss the scaling of the $\Delta\tilde{r}_{n}$ and the $\Lambda_{p}$ with
$\varepsilon$ and $\gamma$ in the athermal limit at the energy minimum. While
in the floppy regime, where $E=0$, we have trivially $\Lambda_{p}=0$ for all
$p=1,\dots,N_{s}$, we discuss here the rigid regime, where $E>0$. For
simplicity, we ignore the existence of zero eigenvalues of the matrix
$\tilde{M}_{N_{s}mn}$ for $m,n=N_{s},\dots,N_{\mathrm{dof}}$ (see also
subsection IV.2).

We focus here for clarity on the case where $\tilde{B}_{p}^{(1)}=0$ for all
$p$, for which Eq. (24) becomes:

$\Lambda_{p}=s_{\underaccent{\dot}{p}}\Delta\tilde{r}_{\underaccent{\dot}{p}}+\frac{1}{2}\tilde{M}_{pmn}\Delta\tilde{r}_{m}\Delta\tilde{r}_{n}+\tilde{w}_{p}\varepsilon+\frac{1}{2}\tilde{B}_{p}^{(2)}\gamma^{2}.$
(93)

The general case of arbitrary $\tilde{B}_{p}^{(1)}$ is discussed in appendix
C.

We focus on the limit of small $\varepsilon$ and $\gamma$. More precisely, we
introduce a parameter $\eta$ and discuss the limit of $\eta\ll 1$ with
$\varepsilon\sim\eta$ and $\gamma\sim\eta^{1/2}$.

We now discuss the scaling of the $\Delta\tilde{r}_{n}$ and the $\Lambda_{p}$
with $\eta$. We assume that $\Delta\tilde{r}_{n}$ can be expanded in power
series of $\eta^{1/2}$, and it follows that the same is true for the
$\Lambda_{p}$. With coefficients $\rho_{n}^{(k/2)}$ and $\lambda_{p}^{(k/2)}$,
we respectively have:

$\displaystyle\Delta\tilde{r}_{n}$
$\displaystyle=\sum_{k=0}^{\infty}{\rho_{n}^{(k/2)}\eta^{k/2}},$ (94)
$\displaystyle\Lambda_{p}$
$\displaystyle=\sum_{k=0}^{\infty}{\lambda_{p}^{(k/2)}\eta^{k/2}}.$ (95)

First, we find that $\rho_{n}^{(0)}=0$ for all $n=1,\dots,N_{\mathrm{dof}}$,
because at the transition, where
$\eta\sim\varepsilon+b_{\varepsilon}\gamma^{2}\rightarrow 0$, all
$\Delta\tilde{r}_{n}=0$. Then, insertion of Eqs. (94) and (95) into Eq. (93)
yields:

$\displaystyle\lambda_{p}^{(0)}$ $\displaystyle=0$ (96)
$\displaystyle\lambda_{p}^{(1/2)}$
$\displaystyle=s_{\underaccent{\dot}{p}}\rho_{\underaccent{\dot}{p}}^{(1/2)}$
(97) $\displaystyle\lambda_{p}^{(1)}$
$\displaystyle=s_{\underaccent{\dot}{p}}\rho_{\underaccent{\dot}{p}}^{(1)}+\frac{1}{2}\tilde{M}_{pmn}\rho_{m}^{(1/2)}\rho_{n}^{(1/2)}+\mathrm{const.}$
(98)

To further discuss the values of $\rho_{n}^{(k/2)}$ and $\lambda_{p}^{(k/2)}$,
we use the energy minimum condition:

$0=\frac{1}{E_{0}}\frac{\partial
E}{\partial\Delta\tilde{r}_{n}}\qquad\text{for
{\hbox{n=1,\ldots,N_\mathrm{dof}}}.}$ (99)

From this we obtain using Eqs. (25) and (93):

$\displaystyle 0$
$\displaystyle=s_{\underaccent{\dot}{p}}\Lambda_{\underaccent{\dot}{p}}+\Lambda_{q}\tilde{M}_{qmp}\Delta\tilde{r}_{m}$
for $p=1,\dots,N_{s}-1$. (100) $\displaystyle 0$
$\displaystyle=\Lambda_{q}\tilde{M}_{qmn}\Delta\tilde{r}_{m}$ for
$n=N_{s},\dots,N_{\mathrm{dof}}$. (101)

On the right-hand side of both equations, sums over $q=1,\dots,N_{s}$ and
$m=1,\dots,N_{\mathrm{dof}}$ are implied.

To order $\eta^{1/2}$, Eq. (100) implies that $\lambda_{p}^{(1/2)}=0$ for
$p=1,\dots,N_{s}-1$. With Eq. (97), this implies $\rho_{p}^{(1/2)}=0$ for
$p=1,\dots,N_{s}-1$. Furthermore, since $s_{N_{s}}=0$, Eq. (97) also implies
that $\lambda_{N_{s}}^{(1/2)}=0$.

To order $\eta^{1}$, Eq. (100) also implies that $\lambda_{p}^{(1)}=0$ for
$p=1,\dots,N_{s}-1$. With Eq. (98) follows that

$0=s_{\underaccent{\dot}{p}}\rho_{\underaccent{\dot}{p}}^{(1)}+\frac{1}{2}\tilde{M}_{pmn}\rho_{m}^{(1/2)}\rho_{n}^{(1/2)}+\mathrm{const.}\quad~{}~{}\text{for
{\hbox{p<N_s}}.}$ (102)

Since $\rho_{n}^{(1/2)}=0$ for $n<N_{s}$, the implied sums in the second term
effectively run only over $m,n=N_{s},\dots,N_{\mathrm{dof}}$.

However, $\lambda_{N_{s}}^{(1)}$ is generally _non-zero_ with

$\lambda_{N_{s}}^{(1)}=\frac{1}{2}\tilde{M}_{N_{s}mn}\rho_{m}^{(1/2)}\rho_{n}^{(1/2)}+\mathrm{const.}$
(103)

Again, since $\rho_{n}^{(1/2)}=0$ for $n<N_{s}$, the implied sums in the
second term effectively run only over $m,n=N_{s},\dots,N_{\mathrm{dof}}$.

Consequentially, up to quadratic order in $\eta$, the system energy is
dominated by $\Lambda_{N_{s}}$, and Eq. (25) becomes
$E=E_{0}\Lambda_{N_{s}}^{2}/2$.

Finally, to order $\eta^{3/2}$, Eq. (101) together with the fact that in the
rigid regime $\lambda_{N_{s}}^{(1)}\neq 0$ implies that $\rho_{n}^{(1/2)}=0$
for $n\geq N_{s}$, where we ignore here zero eigen vectors of the matrix
$\tilde{M}_{N_{s}mn}$. Since we showed above that $\rho_{n}^{(1/2)}=0$ also
for $n<N_{s}$, we conclude that $\rho_{n}^{(1/2)}=0$ for all
$n=1,\dots,N_{\mathrm{dof}}.$

###  B.2 Stiff-spring limit

In the stiff-spring limit, the arguments are similar to the athermal case,
just that here, we already know from Eq. (40) that $\Lambda_{p}=0$ for all
$p=1,\dots,N_{s}$. With Eq. (97) follows that $\rho_{p}^{(1/2)}=0$ for
$p=1,\dots,N_{s}-1$. As a difference to the athermal limit, where also the
$\rho_{n}^{(1/2)}$ with $n=N_{s},\dots,N_{\mathrm{dof}}$ vanish, they take a
finite value in this case. Hence, $\Delta\tilde{r}_{p}\sim\eta^{1}$ for
$p<N_{s}$, and $\Delta\tilde{r}_{n}\sim\eta^{1/2}$ for $n\geq N_{s}$. As a
consequence, in all $\Lambda_{p}$ with $p=1,\dots,N_{s}$ all terms
$\sim\Delta\tilde{r}_{m}\Delta\tilde{r}_{n}$ with $m<N_{s}$ or $n<N_{s}$ can
be neglected.

###  B.3 General case

In the general case, we start from the partition sum, Eq. (60), which reads
with Eq. (25):

$Z\sim\int{\left(\prod_{n=1}^{N_{\mathrm{dof}}}{\mathrm{d}\Delta\tilde{r}_{n}}\right)\exp{\left[-\frac{\beta
E_{0}}{2}\sum_{p=1}^{N_{s}}{\Lambda_{p}^{2}}\right]}}.$ (104)

Using Eqs. (95) and (97), this transforms into:

$\begin{split}Z\sim&\int{\left(\prod_{n=1}^{N_{\mathrm{dof}}}{\mathrm{d}\Delta\tilde{r}_{n}}\right)}\\\
&\quad\times\exp{\left[-\frac{\beta
E_{0}\eta}{2}\sum_{p=1}^{N_{s}}{s_{p}^{2}[\rho_{p}^{(1/2)}]^{2}}\right]}\\\
&\quad\times\exp{\left[-\beta
E_{0}\eta^{3/2}\sum_{p=1}^{N_{s}}{s_{p}\rho_{p}^{(1/2)}\lambda_{p}^{(1)}}\right]}\\\
&\quad\times\exp{\left[-\frac{\beta
E_{0}\eta^{2}}{2}\sum_{p=1}^{N_{s}}{\Big{(}[\lambda_{p}^{(1)}]^{2}+2s_{p}\rho_{p}^{(1/2)}\lambda_{p}^{(3/2)}\Big{)}}\right]}\\\
&\quad\times\exp{\left[-\frac{\beta E_{0}}{2}\text{(higher-order terms in
$\eta$)}\right]}.\end{split}$ (105)

For $k_{B}T\lesssim E_{0}[\varepsilon+b_{\varepsilon}\gamma^{2}]^{2}$, we have
$\beta E_{0}\eta^{1}\gtrsim\eta^{-1}\gg 1$, which leads to an exponential
suppression of those terms in the integral where $\rho_{p}^{(1/2)}\neq 0$ for
any $p=1,\dots,N_{s}-1$. Consequentially, to lowest order,
$\Delta\tilde{r}_{p}=\rho_{p}^{(1)}\eta$ for $p=1,\dots,N_{s}-1$, and
$\Delta\tilde{r}_{n}=\rho_{p}^{(1/2)}\eta^{1/2}$ for
$n=1,\dots,N_{\mathrm{dof}}$. We thus obtain, using Eq. (98):

$\begin{split}Z&\sim\eta^{(N_{\mathrm{dof}}+N_{s}-1)/2}\\\
&\times\int{\left(\prod_{p=1}^{N_{s}-1}{\mathrm{d}\rho_{p}^{(1)}}\right)\left(\prod_{n=N_{s}}^{N_{\mathrm{dof}}}{\mathrm{d}\rho_{n}^{(1/2)}}\right)}\\\
&\quad\times\exp{\Bigg{[}-\frac{\beta
E_{0}\eta^{2}}{2}\sum_{p=1}^{N_{s}}{\bigg{(}s_{p}\rho_{p}^{(1)}}}\\\
&\qquad\quad~{}+\frac{1}{2}\sum_{m,n=N_{s}}^{N_{\mathrm{dof}}}{\tilde{M}_{pmn}\rho_{m}^{(1/2)}\rho_{n}^{(1/2)}}+\mathrm{const.}\bigg{)}^{2}\Bigg{]}.\end{split}$
(106)

Since the terms $\sim\eta^{2}$, are generally non-zero, in particular for
$p=N_{s}$, we neglect in the exponential terms of higher order than
$\eta^{2}$. Since in the exponential, any term with sum index $p$ does not
depend on any of the $\rho_{q}^{(1)}$ with $q\neq p$ (and $q<N_{s}$), we can
rearrange the integral into:

$\begin{split}Z&\sim\eta^{(N_{\mathrm{dof}}+N_{s}-1)/2}\\\
&\times\int{\left(\prod_{n=N_{s}}^{N_{\mathrm{dof}}}{\mathrm{d}\rho_{n}^{(1/2)}}\right)}\exp{\Bigg{[}-\frac{\beta
E_{0}\eta^{2}}{2}}\\\
&\qquad\quad\times\bigg{(}\frac{1}{2}\sum_{m,n=N_{s}}^{N_{\mathrm{dof}}}{\tilde{M}_{N_{s}mn}\rho_{m}^{(1/2)}\rho_{n}^{(1/2)}}+\mathrm{const.}\bigg{)}^{2}\Bigg{]}\\\
&\qquad\times\prod_{p=1}^{N_{s}-1}{\Bigg{(}}\int{\mathrm{d}\rho_{p}^{(1)}}\exp{\Bigg{[}}-\frac{\beta
E_{0}\eta^{2}}{2}\bigg{(}s_{p}\rho_{p}^{(1)}\\\
&\qquad\qquad+\frac{1}{2}\sum_{m,n=N_{s}}^{N_{\mathrm{dof}}}{\tilde{M}_{pmn}\rho_{m}^{(1/2)}\rho_{n}^{(1/2)}}+\mathrm{const.}\bigg{)}^{2}\Bigg{]}\Bigg{)}.\end{split}$
(107)

This corresponds to Eq. (62) in the main text.

##  Appendix C Non-zero first-order terms, $\tilde{B}^{(1)}_{p}$, in $\gamma$

Again to simplify our discussion here, we ignore the existence of mechanisms,
and more generally the possibility of zero modes of the matrix
$\tilde{M}_{N_{s}mn}$ for $m,n=N_{s},\dots,N_{\mathrm{dof}}$, i.e. this matrix
is strictly positive definite.

To discuss the general case of non-vanishing $\tilde{B}^{(1)}_{p}$, we first
simplify by applying shifts to both $\varepsilon$ and $\gamma$,

$\displaystyle\bar{\varepsilon}$
$\displaystyle:=\varepsilon-\Delta\varepsilon$ (108)
$\displaystyle\bar{\gamma}$ $\displaystyle:=\gamma-\Delta\gamma.$ (109)

We define $\Delta\varepsilon$ and $\Delta\gamma$ such that in the athermal
limit $(\bar{\varepsilon},\bar{\gamma})=(0,0)$ is on the transition line, and
such that after carrying out the Taylor expansion in Eq. (22) with respect to
$\bar{\gamma}$ instead of $\gamma$, the linear-order $\bar{\gamma}$ term in
$\Lambda_{N_{s}}$ vanishes at the energy minimum:

$\Lambda_{N_{s}}=\frac{1}{2}\tilde{M}_{N_{s}mn}\Delta\tilde{r}_{m}\Delta\tilde{r}_{n}+\tilde{w}_{N_{s}}\bar{\varepsilon}+\frac{1}{2}B_{N_{s}}^{(2)}\bar{\gamma}^{2}.$
(110)

Note that the $\Lambda_{p}$ with $p<N_{s}$ do in general still have a linear-
order term in $\bar{\gamma}$:

$\Lambda_{p}=s_{\underaccent{\dot}{p}}\Delta\tilde{r}_{\underaccent{\dot}{p}}+\frac{1}{2}\tilde{M}_{pmn}\Delta\tilde{r}_{m}\Delta\tilde{r}_{n}+\tilde{w}_{p}\hat{\varepsilon}+B_{p}^{(1)}\bar{\gamma}+\frac{1}{2}B_{p}^{(2)}\bar{\gamma}^{2}.$
(111)

Also, the $\Delta\tilde{r}_{n}$ may also be redefined due to the shifts in
$\varepsilon$ and $\gamma$ due to a redefinition according to subsection II.3
and a change in the compatibility matrix. If $\Delta\varepsilon$ and
$\Delta\gamma$ are sufficiently small such that to lowest order $C$, $M$, and
$B^{(2)}$ in Eq. (18) are unaffected, it can be shown that the shifts are
given by
$\Delta\varepsilon=(\tilde{B}^{(1)}_{N_{s}})^{2}/4\tilde{w}_{N_{s}}^{2}\\!b_{\varepsilon}$
and
$\Delta\gamma=-\tilde{B}^{(1)}_{N_{s}}/2\tilde{w}_{N_{s}}\\!b_{\varepsilon}$,
where $\tilde{B}^{(1)}_{N_{s}}$ is the linear-order term before the shift, and
$b_{\varepsilon}$ is given by Eq. (116) below.

We first discuss the athermal limit using arguments similar to appendix B.1.
Using $\bar{\eta}$ as small parameter and fixing
$\bar{\varepsilon}\sim\bar{\eta}$ and $\bar{\gamma}\sim\bar{\eta}^{1/2}$, we
use Eq. (100) to order $\bar{\eta}^{1/2}$ to conclude that for any $p<N_{s}$
up to order $\bar{\eta}^{1/2}$:

$0=\Lambda_{p}=s_{\underaccent{\dot}{p}}\Delta\tilde{r}_{\underaccent{\dot}{p}}+\tilde{B}^{(1)}_{p}\bar{\gamma}.$
(112)

Thus, to order $\bar{\eta}^{1/2}$:

$\Delta\tilde{r}_{p}=-\frac{\tilde{B}^{(1)}_{p}}{s_{p}}\bar{\gamma}\qquad\text{for
{\hbox{p<N_s}}.}$ (113)

Inserting this into Eq. (110) together with the diagonalization of the matrix
$\tilde{M}_{N_{s}mn}$ with $m,n\geq N_{s}$, Eq. (44), we obtain up to order
$\eta^{1}$:

$\Lambda_{N_{s}}=\frac{1}{2}\sum_{k=1}^{N_{\mathrm{dof}}-N_{s}+1}{\mu_{k}\Big{(}\Delta
q_{k}-\Delta
q^{0}_{k}\Big{)}^{2}}+\tilde{w}_{N_{s}}\big{[}\bar{\varepsilon}+b_{\varepsilon}\bar{\gamma}^{2}\big{]},$
(114)

where we introduced

$\Delta q^{0}_{k}:=\frac{\bar{\gamma}}{\mu_{k}}\sum_{\begin{subarray}{c}m\geq
N_{s}\\\
p<N_{s}\end{subarray}}{v^{k}_{m}\tilde{M}_{N_{s}mp}\frac{\tilde{B}^{(1)}_{p}}{s_{p}}}$
(115)

and

$\begin{split}b_{\varepsilon}&:=\frac{1}{2\tilde{w}_{N_{s}}}\bBigg@{5}[\tilde{B}^{(2)}_{N_{s}}+\sum_{p,q<N_{s}}{\frac{\tilde{B}^{(1)}_{p}}{s_{p}}\tilde{M}_{N_{s}pq}\frac{\tilde{B}^{(1)}_{q}}{s_{q}}}\\\
&\qquad\qquad\quad-\sum_{\begin{subarray}{c}p,q<N_{s}\\\ m,n\geq N_{s}\\\
<EMAIL_ADDRESS>(116)

Note that we assume here for simplicity that $\mu_{k}>0$.

The additional terms in $b_{\varepsilon}$ in Eq. (116) as compared to Eq. (29)
scale as $\sim\tilde{B}^{(1)}_{p}\tilde{B}^{(1)}_{q}$. The symmetry of the
system with respect to the sign of $\gamma$ combined with the central limit
theorem suggest that, in general, these linear terms scale with system size as
$\tilde{B}^{(1)}_{p}\sim N_{s}^{-1/2}$. Indeed, we numerically found earlier
for phanomized triangular networks that
$\tilde{B}^{(1)}_{N_{s}}\sim\Delta\gamma\sim N_{\mathrm{dof}}^{-1/2}\sim
N_{s}^{-1/2}$ [4]. We thus expect that the additional contributions in
$b_{\varepsilon}$ scale as $N_{s}^{-1}$.

For the stiff-spring limit and the general case of finite spring stiffness and
temperature, these arguments are essentially the same. In the stiff-spring
limit, for non-vanishing $\tilde{B}^{(1)}_{p}$, the ellipsoid that describes
the possible system configurations is now centered at $(\Delta q_{k}^{0})$
(compare Eqs. (46) and (114)). According to Eq. (115) this ellipsoid center
displaces linearly with shear strain $\bar{\gamma}$.

In the general case, a similar line of argument holds as in appendix B.3.
However, the quantities that become zero are not $\rho_{p}^{(1/2)}$ for
$p<N_{s}$, but instead
$\rho_{p}^{(1/2)}+(\tilde{B}^{(1)}_{p}/s_{p})\bar{\gamma}/\bar{\eta}$.

In our Taylor expansion, Eq. (18), we have left out terms that are linear in
both $\Delta\tilde{r}_{n}$ and $\gamma$. We expect such terms to also scale
with system size as $\sim N_{s}^{-1/2}$. They create additional contributions
to $\Delta q^{0}_{k}$ and, in interaction with the $\tilde{B}^{(1)}_{p}$
terms, contributions to $b_{\varepsilon}$. However, other than adding a
correction to the value of $b_{\varepsilon}$ we do not expect them to affect
the results presented in the main text.

##  Appendix D Relation between the Hessian and the matrix
$\tilde{M}_{N_{s}mn}$

We evaluate the Hessian of the system,

$\begin{split}H_{mn}&:=\frac{\partial^{2}E}{\partial\Delta\tilde{r}_{m}\partial\Delta\tilde{r}_{n}},\end{split}$
(117)

_in the rigid athermal regime_ , at an energy minimum close to the transition
point, where we assume for simplicity that $\tilde{B}_{p}^{(1)}=0$ for all
$p\leq N_{s}$.

Using Eq. (25):

$\begin{split}H_{mn}&=E_{0}\sum_{p=1}^{N_{s}}{\left(\frac{\partial\Lambda_{p}}{\partial\Delta\tilde{r}_{m}}\,\frac{\partial\Lambda_{p}}{\partial\Delta\tilde{r}_{n}}+\Lambda_{p}\frac{\partial^{2}\Lambda_{p}}{\partial\Delta\tilde{r}_{m}\partial\Delta\tilde{r}_{n}}\right)}.\end{split}$
(118)

Further, according to our discussion in appendix B.1, at an energy minimum,
all $\Delta\tilde{r}_{n}$ are at least of order $\sim\eta^{1}$. Also,
$\Lambda_{p}=0$ up to order $\sim\eta^{1}$ for $p<N_{s}$. Thus, we have up to
order $\sim\eta^{1}$:

$\begin{split}H_{mn}&=E_{0}\sum_{p=1}^{N_{s}-1}{s_{p}^{2}\delta_{pm}\delta_{pn}}\\\
&\quad+E_{0}\sum_{p=1}^{N_{s}-1}\sum_{l=1}^{N_{\mathrm{dof}}}{s_{p}\Delta\tilde{r}_{l}\Big{(}\delta_{pm}\tilde{M}_{pln}+\delta_{pn}\tilde{M}_{plm}\Big{)}}\\\
&\quad+E_{0}\Lambda_{N_{s}}\tilde{M}_{N_{s}mn}.\end{split}$ (119)

In this expression, the first term, which is of absolute order in
$\varepsilon$ and $\gamma$, corresponds to the Hessian at the transition
point. The second and third terms are of order $\sim\eta^{1}$.

Eq. (119) shows that to lowest order, the $N_{\mathrm{1st}}^{\ast}$ nonzero
eigen modes of $\tilde{M}_{N_{s}mn}$ for $m,n\geq N_{s}$, are those eigen
modes of $H_{mn}$ that have zero eigen value at the transition point, and
whose eigen values become nonzero in the rigid regime, scaling as
$\sim\Lambda_{N_{s}}\sim\eta^{1}$. This is because for $m,n\geq N_{s}$ only
the third term in Eq. (119) can be nonzero.

##  Appendix E Partition sum for stiff-spring limit

In Eq. (38), to define the configurational phase space volume $\Omega$, we use
a Dirac delta for each of the $N_{s}$ spring lengths in the system. One could
instead have the idea to use a single Dirac delta for the system energy $E$,
since $E=0$ is equivalent to all springs simultaneously attaining their rest
length. However, both definitions will lead to different results. This is
because the integration of a Dirac delta involves the evaluation of the
Jacobian of the argument of the Dirac delta with respect to the integration
variable. Since the Jacobian of the energy creates additional prefactors that
correspond to the prestresses in the system, such prefactors would
additionally appear when using a Dirac delta of the energy.

To see why the definition Eq. (38) is the correct one in the stiff-spring
limit, we start form the partition sum:

$Z=\int{\left(\prod_{n=1}^{N_{\mathrm{dof}}}{\mathrm{d}R_{n}}\right)\exp{\left(-\frac{\beta}{2}\sum_{i=1}^{N_{s}}{K_{i}\big{(}L_{i}-L_{0i}\big{)}^{2}}\right)}}.$
(120)

A Gaussian converges to a Dirac Delta in the limit of small standard
deviation:

$\lim_{K_{i}\rightarrow\infty}\sqrt{\frac{\beta
K_{i}}{2\pi}}\exp{\left(-\frac{\beta
K_{i}(L_{i}-L_{0i})^{2}}{2}\right)}=\delta(L_{i}-L_{0i}).$ (121)

Thus, whenever all $K_{i}$ are large enough such that for all
$i=1,\dots,N_{s}$ the $(\beta K_{i})^{-1/2}$ become much smaller than all
relevant length scales of the system, the partition sum can be expressed as:

$\begin{split}Z&=\left(\frac{2\pi}{\beta}\right)^{N_{s}/2}\prod_{i=1}^{N_{s}}{K_{i}^{-1/2}}\\\
&\qquad\qquad\times\int{\left(\prod_{n=1}^{N_{\mathrm{dof}}}{\mathrm{d}R_{n}}\right)\prod_{i=1}^{N_{s}}{\delta\big{(}L_{i}-L_{0i}\big{)}}}.\end{split}$
(122)

This indicates that the correct approach is to use Dirac deltas on the spring
lengths rather than the energy $E$. Note that in Eq. (38), we have ignored the
first prefactor that only depends on temperature, and so is not relevant for
the elastic system properties. Moreover, while the second prefactor does not
depend on strain or temperature, it converges to zero for
$K_{i}\rightarrow\infty$. However, to discuss the elastic system properties,
we only need derivatives of $\log Z$, where this term only creates a constant
albeit diverging offset.

##  Appendix F Higher-order dependencies of $\Delta\tilde{\ell}_{N_{s}}$ on
the dofs

Here we outline how some of our arguments change in the presence of $K$th-
order zero modes with $K>1$. In particular, we discuss how the phase space
volume scaling with strain changes.

We start from the following specific generalization of the expansion of
$\Delta\tilde{\ell}_{N_{s}}$ from Eqs. (22), (46):

$\Delta\tilde{\ell}_{N_{s}}=\sum_{K=1}^{\infty}\frac{1}{(K+1)!}\sum_{k=1}^{N_{K\mathrm{th}}^{\ast}}{\mu_{K,k}\Delta
q_{K,k}^{K+1}}+\frac{1}{2}\tilde{B}_{N_{s}}^{(2)}\gamma^{2}.$ (123)

Here, $(K+1)!$ denotes the faculty of $K+1$, the variables $\Delta q_{K,k}$
denote the amplitudes of the $K$th-order zero modes, and the associated
coefficients are $\mu_{K,k}>0$. For simplicity, we set again
$\tilde{B}_{p}^{(1)}=0$, and we ignore mixed terms between the different
$\Delta q_{K,k}$.

Assuming that the SSS mode $\Lambda_{N_{s}}$ still dominates the energy in
this case, we find that $N_{K\mathrm{th}}^{\ast}=0$ for even $K$. This is
because otherwise, a $K$th-order zero mode for even $K$ would allow for states
with $\Lambda_{N_{s}}=0$ for $\varepsilon>0$ and $\gamma=0$, even though the
system is supposed to be rigid with $E>0$ in this regime. $\Lambda_{N_{s}}$
can become zero for $N_{K\mathrm{th}}^{\ast}>0$ for even $K$, because a $K$th-
order zero mode would introduce a term $\sim\Delta q_{K,k}^{K+1}$ with an odd
exponent, and consequentially, $\Delta q_{K,k}$ can be adjusted by the energy
minimization such that $\Lambda_{N_{s}}=0$. Taken together, this implies that
the sum over $K$ in Eq. (123) actually only includes the odd integers.

Thus, the phase space volume is governed by an integral of the following kind:

$\begin{split}&\Omega_{N_{s}}\sim
I(\varepsilon,\gamma):=\int{\prod_{K=1}^{\infty}\prod_{k=1}^{N_{K\mathrm{th}}^{\ast}}{\mathrm{d}\Delta
q_{K,k}}}\\\
&\times\delta\left(\sum_{K=1}^{\infty}\frac{1}{(K+1)!}\sum_{k=1}^{N_{K\mathrm{th}}^{\ast}}{\mu_{K,k}\Delta
q_{K,k}^{K+1}}+\tilde{w}_{N_{s}}\Big{[}\varepsilon+b_{\varepsilon}\gamma^{2}\Big{]}\right),\end{split}$
(124)

To extract the dependency of $I$ on $\varepsilon$ and $\gamma$, we first
define:

$\begin{split}I_{1}&:=\int{\prod_{K=1}^{\infty}\prod_{k=1}^{N_{K\mathrm{th}}^{\ast}}{\mathrm{d}x_{K,k}}}\\\
&\qquad\quad\times\delta\left(\sum_{K=1}^{\infty}\frac{1}{(K+1)!}\sum_{k=1}^{N_{K\mathrm{th}}^{\ast}}{\mu_{K,k}x_{K,k}^{K+1}}+\tilde{w}_{N_{s}}\right),\end{split}$
(125)

_which depends neither on $\varepsilon$ nor on $\gamma$_. Now we can transform
the expression in Eq. (124) using the following substitutions for all
integration variables:

$\Delta
q_{K,k}=\big{[}\varepsilon+b_{\varepsilon}\gamma^{2}\big{]}^{1/(K+1)}x_{K,k},$
(126)

which yields:

$\Omega_{N_{s}}(\varepsilon,\gamma)\sim
I(\varepsilon,\gamma)=\big{[}\varepsilon+b_{\varepsilon}\gamma^{2}\big{]}^{N_{\Omega}}I_{1}$
(127)

with:

$N_{\Omega}=-1+\sum_{K=1}^{\infty}{\frac{N_{K\mathrm{th}}^{\ast}}{K+1}}.$
(128)

The first term in $N_{\Omega}$ arises from removing the factor
$[\varepsilon+b_{\varepsilon}\gamma^{2}]$ from the Dirac delta, and all the
remaining terms arise from the substitution of the integration variables. In
this argument, we have excluded modes that leave $\Delta\tilde{\ell}_{N_{s}}$
invariant. They do not contribute to $N_{\Omega}$. This is because any
substitution like in Eq. (126), while creating a prefactor, also rescales the
integral bounds by the same prefactor, both of which exactly cancel each other
out. Meanwhile, for any mode that _does_ affect $\Delta\tilde{\ell}_{N_{s}}$,
the rescaling of the integral bounds does not matter, since the integration
domain of these modes is effectively limited by the Dirac delta, not the
integral bounds.

Note that for the case discussed in the main text, i.e. for
$N_{\mathrm{1st}}^{\ast}>0$ and $N_{K\mathrm{th}}^{\ast}=0$ for $K>1$, we have
indeed $N_{\Omega}=(N_{\mathrm{1st}}^{\ast}-2)/2$, consistent with Eq. (128).
Taken together, as indicated by Eq. (128), also higher-order zero modes at the
transition may contribute to the exponent $N_{\Omega}$; they are just
attenuated according to their order $K$. If the SSS mode still dominates the
phase space volume scaling, we also have
$\kappa_{S}=N_{\Omega}k_{B}/D\mathcal{V}^{\ast}$.

##  Appendix G Tension and shear modulus in the general case

To obtain the isotropic tension $t:=(\partial
F/\partial\varepsilon)/D\mathcal{V}$, we first use the chain rule. To lowest
order, $\mathcal{V}=\mathcal{V}^{\ast}$, and so we obtain from Eq. (74):

$\begin{split}t=\frac{1}{D\mathcal{V}^{\ast}}\Bigg{[}&\frac{\partial\bar{F}}{\partial\hat{\varepsilon}_{S}}\frac{\partial\hat{\varepsilon}_{S}}{\partial\varepsilon}+\frac{\partial\bar{F}}{\partial\varepsilon}\\\
&\qquad+\frac{1}{2\beta\bar{F}^{\prime\prime}}\left(\frac{\partial\bar{F}^{\prime\prime}}{\partial\hat{\varepsilon}_{S}}\frac{\partial\hat{\varepsilon}_{S}}{\partial\varepsilon}+\frac{\partial\bar{F}^{\prime\prime}}{\partial\varepsilon}\right)\Bigg{]}.\end{split}$
(129)

We compute to lowest order:

$\displaystyle\frac{\partial\hat{\varepsilon}_{S}}{\partial\varepsilon}$
$\displaystyle=-\frac{\hat{\varepsilon}_{S}+b_{\varepsilon}\gamma^{2}}{\left|\varepsilon+b_{\varepsilon}\gamma^{2}\right|\sqrt{1+\theta}}$
(130)
$\displaystyle\frac{1}{D\mathcal{V}^{\ast}}\frac{\partial\bar{F}}{\partial\hat{\varepsilon}_{S}}$
$\displaystyle=-\kappa_{E}(\varepsilon-\hat{\varepsilon}_{S})-\frac{\kappa_{S}T}{\hat{\varepsilon}_{S}+b_{\varepsilon}\gamma^{2}}=0$
(131)
$\displaystyle\frac{1}{D\mathcal{V}^{\ast}}\frac{\partial\bar{F}}{\partial\varepsilon}$
$\displaystyle=\kappa_{E}(\varepsilon-\hat{\varepsilon}_{S})$ (132)
$\displaystyle\frac{1}{D\mathcal{V}^{\ast}}\frac{\partial\bar{F}^{\prime\prime}}{\partial\hat{\varepsilon}_{S}}$
$\displaystyle=-\frac{2\kappa_{S}T}{(\hat{\varepsilon}_{S}+b_{\varepsilon}\gamma^{2})^{3}}$
(133)
$\displaystyle\frac{1}{D\mathcal{V}^{\ast}}\frac{\partial\bar{F}^{\prime\prime}}{\partial\varepsilon}$
$\displaystyle=0.$ (134)

We thus obtain:

$\begin{split}t&=\kappa_{E}(\varepsilon-\hat{\varepsilon}_{S})\left[1+\frac{1}{2(N_{\mathrm{1st}}^{\ast}-2)}\,\frac{\theta}{1+\theta}\right],\end{split}$
(135)

where we have also used that

$\bar{F}^{\prime\prime}=-D\mathcal{V}^{\ast}\kappa_{E}\frac{\left|\varepsilon+b_{\varepsilon}\gamma^{2}\right|\sqrt{1+\theta}}{\hat{\varepsilon}_{S}+b_{\varepsilon}\gamma^{2}}$
(136)

and

$\kappa_{E}(\varepsilon-\hat{\varepsilon}_{S})=-\frac{\kappa_{S}T}{\hat{\varepsilon}_{S}+b_{\varepsilon}\gamma^{2}},$
(137)

which follows from Eq.(72). The term $\sim(N_{\mathrm{1st}}^{\ast}-2)^{-1}$ in
Eq. (135) corresponds to the terms involving $\bar{F}^{\prime\prime}$.

For the shear modulus $G:=(\partial^{2}F/\partial\gamma^{2})/\mathcal{V}$, we
have to lowest order:

$\begin{split}G=\frac{1}{\mathcal{V}^{\ast}}\Bigg{[}&\frac{\partial^{2}\bar{F}(\hat{\varepsilon}_{S}(\varepsilon,\gamma);\varepsilon,\gamma)}{\partial\gamma^{2}}\\\
&\qquad+\frac{1}{2\beta}\frac{\partial^{2}}{\partial\gamma^{2}}\log{\bar{F}^{\prime\prime}(\hat{\varepsilon}_{S}(\varepsilon,\gamma);\varepsilon,\gamma)}\Bigg{]}.\end{split}$
(138)

We compute with Eqs. (131) and (137):

$\begin{split}\frac{1}{\mathcal{V}^{\ast}}\frac{\partial\bar{F}(\hat{\varepsilon}_{S}(\varepsilon,\gamma);\varepsilon,\gamma)}{\partial\gamma}&=\frac{1}{\mathcal{V}^{\ast}}\frac{\partial\bar{F}(\hat{\varepsilon}_{S};\varepsilon,\gamma)}{\partial\gamma}\\\
&=-\frac{2D\kappa_{S}Tb_{\varepsilon}\gamma}{\hat{\varepsilon}_{S}+b_{\varepsilon}\gamma^{2}}\\\
&=2D\kappa_{E}(\varepsilon-\hat{\varepsilon}_{S})b_{\varepsilon}\gamma,\end{split}$
(139)

and further:

$\begin{split}&\frac{1}{\mathcal{V}^{\ast}}\frac{\partial^{2}\bar{F}(\hat{\varepsilon}_{S}(\varepsilon,\gamma);\varepsilon,\gamma)}{\partial\gamma^{2}}\\\
&=\frac{1}{\mathcal{V}^{\ast}}\Bigg{[}\frac{\partial^{2}\bar{F}(\hat{\varepsilon}_{S};\varepsilon,\gamma)}{\partial\gamma\partial\hat{\varepsilon}_{S}}\frac{\partial\hat{\varepsilon}_{S}}{\partial\gamma}+\frac{\partial^{2}\bar{F}(\hat{\varepsilon}_{S};\varepsilon,\gamma)}{\partial\gamma^{2}}\Bigg{]}\\\
\end{split}$ (140)

with

$\displaystyle\frac{1}{\mathcal{V}^{\ast}}\frac{\partial^{2}\bar{F}(\hat{\varepsilon}_{S};\varepsilon,\gamma)}{\partial\gamma\partial\hat{\varepsilon}_{S}}$
$\displaystyle=-2D\kappa_{E}b_{\varepsilon}\gamma$ (141)
$\displaystyle\frac{\partial\hat{\varepsilon}_{S}}{\partial\gamma}$
$\displaystyle=-b_{\varepsilon}\gamma\left(1+\frac{\varepsilon+b_{\varepsilon}\gamma^{2}}{|\varepsilon+b_{\varepsilon}\gamma^{2}|\sqrt{1+\theta}}\right)$
(142)
$\displaystyle=-2b_{\varepsilon}\gamma\frac{\varepsilon-\hat{\varepsilon}_{S}}{|\varepsilon+b_{\varepsilon}\gamma^{2}|\sqrt{1+\theta}}$
(143)
$\displaystyle\frac{1}{\mathcal{V}^{\ast}}\frac{\partial^{2}\bar{F}(\hat{\varepsilon}_{S};\varepsilon,\gamma)}{\partial\gamma^{2}}$
$\displaystyle=2D\kappa_{E}(\varepsilon-\hat{\varepsilon}_{S})b_{\varepsilon}.$
(144)

Thus:

$\begin{split}&\frac{1}{\mathcal{V}^{\ast}}\frac{\partial^{2}\bar{F}(\hat{\varepsilon}_{S}(\varepsilon,\gamma);\varepsilon,\gamma)}{\partial\gamma^{2}}\\\
&=2D\kappa_{E}(\varepsilon-\hat{\varepsilon}_{S})b_{\varepsilon}\Bigg{(}1+\frac{2b_{\varepsilon}\gamma^{2}}{|\varepsilon+b_{\varepsilon}\gamma^{2}|\sqrt{1+\theta}}\Bigg{)}.\end{split}$
(145)

Furthermore:

$\displaystyle\frac{\partial\bar{F}^{\prime\prime}(\hat{\varepsilon}_{S}(\varepsilon,\gamma);\varepsilon,\gamma)}{\partial\gamma}$
$\displaystyle=\frac{\partial\bar{F}^{\prime\prime}(\hat{\varepsilon}_{S};\varepsilon,\gamma)}{\partial\hat{\varepsilon}_{S}}\frac{\partial\hat{\varepsilon}_{S}}{\partial\gamma}+\frac{\partial\bar{F}^{\prime\prime}(\hat{\varepsilon}_{S};\varepsilon,\gamma)}{\partial\gamma}$
$\displaystyle=D\mathcal{V}^{\ast}\kappa_{E}\Bigg{[}\left(\frac{1}{\hat{\varepsilon}_{S}+b_{\varepsilon}\gamma^{2}}-\frac{\hat{\varepsilon}_{S}-\varepsilon}{(\hat{\varepsilon}_{S}+b_{\varepsilon}\gamma^{2})^{2}}\right)$
$\displaystyle\qquad\qquad\qquad\quad\times\left(-2b_{\varepsilon}\gamma\frac{\varepsilon-\hat{\varepsilon}_{S}}{|\varepsilon+b_{\varepsilon}\gamma^{2}|\sqrt{1+\theta}}\right)$
$\displaystyle\qquad\qquad\quad-2b_{\varepsilon}\gamma\frac{\hat{\varepsilon}_{S}-\varepsilon}{(\hat{\varepsilon}_{S}+b_{\varepsilon}\gamma^{2})^{2}}\Bigg{]}$
$\displaystyle=-2D\mathcal{V}^{\ast}\kappa_{E}b_{\varepsilon}\gamma\frac{\varepsilon-\hat{\varepsilon}_{S}}{(\hat{\varepsilon}_{S}+b_{\varepsilon}\gamma^{2})^{2}|\varepsilon+b_{\varepsilon}\gamma^{2}|\sqrt{1+\theta}}$
$\displaystyle\qquad\times\Big{[}\varepsilon+b_{\varepsilon}\gamma^{2}-|\varepsilon+b_{\varepsilon}\gamma^{2}|\sqrt{1+\theta}\Big{]}$
$\displaystyle=-4D\mathcal{V}^{\ast}\kappa_{E}b_{\varepsilon}\gamma\frac{\varepsilon-\hat{\varepsilon}_{S}}{(\hat{\varepsilon}_{S}+b_{\varepsilon}\gamma^{2})|\varepsilon+b_{\varepsilon}\gamma^{2}|\sqrt{1+\theta}},$
(146)

and thus:

$\begin{split}\frac{\partial}{\partial\gamma}\log{\bar{F}^{\prime\prime}(\hat{\varepsilon}_{S}(\varepsilon,\gamma);\varepsilon,\gamma)}=4b_{\varepsilon}\gamma\frac{\varepsilon-\hat{\varepsilon}_{S}}{(\varepsilon+b_{\varepsilon}\gamma^{2})^{2}(1+\theta)}.\end{split}$
(147)

The second derivative is:

$\begin{split}&\frac{\partial^{2}}{\partial\gamma^{2}}\log{\bar{F}^{\prime\prime}(\hat{\varepsilon}_{S}(\varepsilon,\gamma);\varepsilon,\gamma)}\\\
&=\frac{4b_{\varepsilon}(\varepsilon-\hat{\varepsilon}_{S})}{(\varepsilon+b_{\varepsilon}\gamma^{2})^{2}(1+\theta)}\\\
&\quad+\frac{4b_{\varepsilon}\gamma}{(\varepsilon+b_{\varepsilon}\gamma^{2})^{2}(1+\theta)}\times
2b_{\varepsilon}\gamma\frac{\varepsilon-\hat{\varepsilon}_{S}}{|\varepsilon+b_{\varepsilon}\gamma^{2}|\sqrt{1+\theta}}\\\
&\quad-\frac{16b_{\varepsilon}^{2}\gamma^{2}(\varepsilon-\hat{\varepsilon}_{S})}{(\varepsilon+b_{\varepsilon}\gamma^{2})^{3}(1+\theta)^{2}}\\\
&=\frac{4b_{\varepsilon}(\varepsilon-\hat{\varepsilon}_{S})}{(\varepsilon+b_{\varepsilon}\gamma^{2})^{2}(1+\theta)}\Bigg{[}1+\frac{2b_{\varepsilon}\gamma^{2}}{|\varepsilon+b_{\varepsilon}\gamma^{2}|\sqrt{1+\theta}}\\\
&\qquad\qquad\qquad\qquad\qquad\quad\;\;-\frac{4b_{\varepsilon}\gamma^{2}}{(\varepsilon+b_{\varepsilon}\gamma^{2})(1+\theta)}\Bigg{]}.\end{split}$
(148)

From Eqs. (58) and (137) follows:

$\begin{split}\frac{1}{2\beta\mathcal{V}^{\ast}}=\frac{D\kappa_{S}T}{N_{\mathrm{1st}}^{\ast}-2}=\frac{D\kappa_{E}}{N_{\mathrm{1st}}^{\ast}-2}\,\theta\frac{(\varepsilon+b_{\varepsilon}\gamma^{2})^{2}}{4}\end{split}$
(149)

Combining Eqs. (138), (145), (148), and (149), we obtain Eq. (76) in the main
text.

## References

  * Maxwell [1864] J. C. Maxwell, Philosophical Magazine Series 4 27, 294 (1864).
  * Calladine [1978] C. R. Calladine, International Journal of Solids and Structures 14, 161 (1978).
  * Merkel  _et al._ [2019] M. Merkel, K. Baumgarten, B. P. Tighe, and M. L. Manning, Proceedings of the National Academy of Sciences 116, 6560 (2019).
  * Lee and Merkel [2022] C.-T. Lee and M. Merkel, Soft Matter 18, 5410 (2022).
  * Lubensky  _et al._ [2015] T. C. Lubensky, C. L. Kane, X. Mao, A. Souslov, and K. Sun, Reports on Progress in Physics 78, 73901 (2015), arxiv:1503.01324 .
  * Mannattil  _et al._ [2022] M. Mannattil, J. M. Schwarz, and C. D. Santangelo, Physical Review Letters 128, 208005 (2022), arxiv:2112.04279 .
  * Damavandi  _et al._ [2022a] O. K. Damavandi, V. F. Hagh, C. D. Santangelo, and M. L. Manning, Physical Review E 105, 025004 (2022a).
  * Damavandi  _et al._ [2022b] O. K. Damavandi, V. F. Hagh, C. D. Santangelo, and M. L. Manning, Physical Review E 105, 025003 (2022b).
  * Alexander [1998] S. Alexander, Physics Report 296, 65 (1998), arxiv:1512.04885 .
  * Cui  _et al._ [2019] B. Cui, G. Ruocco, and A. Zaccone, Granular Matter 21, 69 (2019), arxiv:1901.09582 .
  * Zhang  _et al._ [2022] S. Zhang, E. Stanifer, V. V. Vasisht, L. Zhang, E. Del Gado, and X. Mao, Physical Review Research 4, 043181 (2022).
  * Wyart [2005] M. Wyart, Annales de Physique 30, 1 (2005).
  * Merkel and Manning [2018] M. Merkel and M. L. Manning, New Journal of Physics 20, 022002 (2018).
  * Onck  _et al._ [2005] P. R. Onck, T. Koeman, T. van Dillen, and E. van der Giessen, Physical Review Letters 95, 178102 (2005), arxiv:cond-mat/0502397 .
  * Broedersz and MacKintosh [2014] C. P. Broedersz and F. C. MacKintosh, Reviews of Modern Physics 86, 995 (2014), arxiv:1404.4332 .
  * Licup  _et al._ [2015] A. J. Licup, S. Münster, A. Sharma, M. Sheinman, L. M. Jawerth, B. Fabry, D. a. Weitz, and F. C. MacKintosh, Proceedings of the National Academy of Sciences 112, 9573 (2015), arxiv:1503.00924 .
  * Sharma  _et al._ [2016] A. Sharma, A. J. Licup, K. A. Jansen, R. Rens, M. Sheinman, G. H. Koenderink, and F. C. MacKintosh, Nature Physics 12, 584 (2016), arxiv:1506.07792 .
  * Moshe  _et al._ [2017] M. Moshe, M. J. Bowick, and M. C. Marchetti, Physical Review Letters 120, 268105 (2017), arxiv:1708.07848 .
  * Sussman and Merkel [2018] D. M. Sussman and M. Merkel, Soft Matter 14, 3397 (2018).
  * Wang  _et al._ [2020] X. Wang, M. Merkel, L. B. Sutter, G. Erdemci-Tandogan, M. L. Manning, and K. E. Kasza, Proceedings of the National Academy of Sciences 117, 13541 (2020), arxiv:2005.07283 .
  * Plischke and Joós [1998] M. Plischke and B. Joós, Physical Review Letters 80, 4907 (1998).
  * Dennison  _et al._ [2013] M. Dennison, M. Sheinman, C. Storm, and F. C. MacKintosh, Physical Review Letters 111, 095503 (2013).
  * Wigbers  _et al._ [2015] M. C. Wigbers, F. C. Mackintosh, and M. Dennison, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 92, 1 (2015), arxiv:1410.7860 .
  * Zhang and Mao [2016] L. Zhang and X. Mao, Physical Review E 93, 1 (2016).
  * Chen  _et al._ [2023] S. Chen, T. Markovich, and F. C. MacKintosh, Physical Review Letters 130, 088101 (2023).
  * Note [1] Reference to be inserted later.
  * Parker  _et al._ [2020] A. Parker, M. C. Marchetti, M. L. Manning, and J. M. Schwarz, How does the extracellular matrix affect the rigidity of an embedded spheroid? (2020), arxiv:2006.16203 [cond-mat, q-bio] .
  * Note [2] More precisely, we define the number of mechanisms $N_{\mathrm{mech}}$ of a system as the maximal integer $N_{m}$ for which there is a continuous function $\cc@accent{"717E}{x}:\cc@accent{"717E}{u}\mapstochar\rightarrow\cc@accent{"717E}{x}(\cc@accent{"717E}{u}),\mathbb{R}^{N_{m}}\rightarrow\mathbb{R}^{N_{\mathrm{dof}}}$ and some $\eta>0$ such that Eq. (13) holds within some open ball of radius $\eta$, i.e. for $|\cc@accent{"717E}{u}|<\eta$.
  * Note [3] In the stiff-spring limit, one might naively replace each spring by an exact constraint and compute the partition sum using Lagrangian mechanics. However, it has been shown that such an approach generally leads to a different result than taking the partition sum $Z$ of a system with finite spring stiffnesses $K_{i}$ and taking the limit $K_{i}\rightarrow\infty$. [37]. The latter is of course the physically more relevant approach, since there are no exact geometric constraints in nature.
  * Note [4] One only needs to show that for any given $k$ with $\mu_{k}=0$, there exists a two-times differentiable function $u\mapstochar\rightarrow\cc@accent{"717E}{y}(u)$ with $\cc@accent{"717E}{y}(0)=\cc@accent{"717E}{r}\,^{\ast}$, $\cc@accent{"717E}{y}^{\prime}(0)=\cc@accent{"717E}{z}_{k}$, and $\mathrm{d}^{2}\ell_{i}(\cc@accent{"717E}{y}(u),\gamma)/\mathrm{d}u^{2}=0$. Indeed, such a function is given by:  $\begin{split}y_{m}(u)&:=r_{m}^{\ast}+uV_{mn}v_{n}^{k}\\\ &\hskip 20.00003pt\hskip 20.00003pt-\frac{u^{2}}{2}\sum@\slimits@_{p=1}^{N_{s}-1}{\frac{V_{mp}}{s_{p}}\cc@accent{"707E}{M}_{pn^{\prime}n^{\prime\prime}}v_{n^{\prime}}^{k}v_{n^{\prime\prime}}^{k}}\end{split}$ with $m=1,\dots,N_{\mathrm{dof}}$.
  * Note [5] For any two-times differentiable function $u\mapstochar\rightarrow\cc@accent{"717E}{y}(u)$ with $\cc@accent{"717E}{y}(0)=\cc@accent{"717E}{r}\,^{\ast}$ and $\cc@accent{"717E}{y}^{\prime}(0)=\cc@accent{"717E}{z}_{k}$, it is impossible to have $\mathrm{d}^{2}\ell_{i}(\cc@accent{"717E}{y}(u),\gamma)/\mathrm{d}u^{2}=0$ for all $i=1,\dots,N_{s}$, because $\mathrm{d}^{2}\Delta\cc@accent{"707E}\ell_{N_{s}}(\cc@accent{"717E}{y}(u),\gamma)/\mathrm{d}u^{2}=\mu_{k}>0$.
  * Tuckerman [2010] M. E. Tuckerman,  _Statistical Mechanics: Theory and Molecular Simulation_ , illustrated edition ed. (Oxford University Press, Oxford, 2010).
  * Arzash  _et al._ [2023] S. Arzash, A. Gannavarapu, and F. C. MacKintosh, Mechanical criticality of fiber networks at a finite temperature (2023), arxiv:2304.12459 [cond-mat, physics:physics] .
  * Rocklin  _et al._ [2018] D. Z. Rocklin, V. Vitelli, and X. Mao, Folding mechanisms at finite temperature (2018), arxiv:1802.02704 [cond-mat] .
  * Merkel [2014] M. Merkel,  _From Cells to Tissues: Tissue Remodeling and Polarity Reorientation in Shearing Epithelia_ , Ph.D. thesis, Technische Universität Dresden (2014).
  * Lees and Edwards [1972] A. W. Lees and S. F. Edwards, Journal of Physics C: Solid State Physics 5, 1921 (1972).
  * van Kampen and Lodder [1984] N. G. van Kampen and J. J. Lodder, American Journal of Physics 52, 419 (1984).