# Partition sum of thermal, under-constrained systems

Cheng-Tai Lee, Matthias Merkel

► To cite this version:

Cheng-Tai Lee, Matthias Merkel. Partition sum of thermal, under-constrained systems. Physical Review E, 2024, 110 (6), pp.064147. 10.1103/PhysRevE.110.064147. hal-04251711v2

HAL Id: hal-04251711

https://hal.science/hal-04251711v2

Submitted on 29 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Partition sum of thermal, under-constrained systems

Cheng-Tai Lee and Matthias Merkel*

Aix Marseille Univ, Université de Toulon, CNRS, CPT (UMR 7332),

Turing Center for Living Systems, Marseille, France

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 ε, shear strain γ, and temperature T. These expressions contain only three parameters, entropic rigidity κS, energetic rigidity κE, and a parameter bε describing the interaction between isotropic and shear strain. We provide analytical expressions for these parameters based on the microscopic structure of the system. 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_{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_{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_{dof} - N_s + N_{SSS}, \quad (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, infinitesimal zero modes – while inducing no changes to the constraints to linear order – could induce changes to

higher order, corresponding to higher-order rigidity [7, 8]. Here, we refer to this third kind as Kth-order zero modes with some integer K > 0. Such modes do not affect the constraints to Kth order, but do so to (K + 1)th order (precise definitions in section II D). 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, 3, 7–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, 3, 5, 9, 12, 13]. 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–17] and vertex models for biological tissues [3, 13, 18–20].

Despite these advances for the athermal limit, there is still relatively little known about thermal, i.e. finite-temperature, under-constrained systems [21–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 small isotropic strain ε, shear strain γ, and temperature 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 foundation

<footnote>* matthias.merkel@cnrs.fr</footnote>
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 [26], we numerically verify the analytical results obtained here.

## II. MODEL

### A. 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 $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 the periodic boundary conditions with both kinds of shear.

A network consists of $N_{\text{node}}$ nodes whose position components combine into a $N_{\text{dof}} = DN_{\text{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_{\text{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. 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. 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$:

$$L_i(c\vec{R}, c^D\mathcal{V}, \gamma) = c^{d_i} L_i(\vec{R}, \mathcal{V}, \gamma). \quad (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 \left[ L_i(\vec{R}, \mathcal{V}, \gamma) - L_{0i} \right]^2, \quad (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 $V$, shear $\gamma$, and temperature $T$.

### B. 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 a typical scenario 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 $V^*$, where, consistent with our earlier findings [3, 4], we assume here that the system is floppy for $V < V^*$ and rigid for $V > V^*$.

We define the isotropic strain $\epsilon$ with respect to the critical volume $V^*$:

$$\epsilon := \frac{1}{D} \log \left( \frac{V}{V^*} \right). \quad (4)$$

Thus, for $\gamma = 0$, the transition occurs at $\epsilon = 0$. In other words, $\epsilon < 0$ corresponds to the athermal floppy regime with a vanishing shear modulus, and $\epsilon > 0$ corresponds to the athermal rigid regime with a finite shear modulus.

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, for simplicity of the argument, 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}^*$. For simplicity, we assume here that $\vec{R}^*$ is uniquely defined. A degeneracy in the athermal transition configuration $\vec{R}^*$ 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}^*$ degeneracy by mechanisms in the discussion section.

### C. Dimensionless quantities

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

$$\vec{r} := \frac{1}{\sqrt[\underline{D}}{\mathcal{V}}} \vec{R} = \frac{1}{e^{\epsilon} \sqrt[\underline{D}}{\mathcal{V}^*}} \vec{R}, \quad (5)$$

$$\ell_i := e^{-d_i \epsilon} \frac{L_i(\vec{R}, \mathcal{V}, \gamma)}{L_{0i}}. \quad (6)$$

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

$$\ell_i(\vec{r}, \gamma) = \frac{(\mathcal{V}^*)^{d_i/D}}{L_{0i}} L_i(\vec{r}, 1, \gamma). \quad (7)$$
In other words, when expressed using the dimensionless node positions $\vec{r}$, the dimensionless spring lengths do not depend on isotropic strain.

As discussed in section IIB, for $\gamma = 0$, the transition occurs at $\varepsilon = 0$ and $\vec{r} = \vec{r}^*$, 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(\vec{r} = \vec{r}^*, \gamma = 0) = 1. \qquad (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}}, \qquad (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 (e^{d_i \varepsilon} \ell_i - 1)^2. \qquad (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].

### D. 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_{\text{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. \qquad (11) $$

Here, $r_n$ and $z_n$ with $n = 1, \dots, N_{\text{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 Kth-order zero mode. In the absence of immobile nodes, the system shows translational invariance. We call an infinitesimal zero mode $\vec{z}$ trivial whenever it corresponds to one of the D global translations. Conversely, in the presence of immobile nodes, there are no trivial zero modes.

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} \to \mathbb{R}^{N_{\text{dof}}}$ with $\vec{x}(0) = \vec{r}$ and some range $\eta > 0$ such that for all $i = 1, \dots, N_s$:

$$ \vec{x}'(0) = \vec{z} \quad \text{and} \qquad (12) $$

$$ e^{d_i \varepsilon} \ell_i(\vec{x}(u), \gamma) = 1 \quad \text{for } |u| < \eta, \qquad (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_{\text{mech}}$ [28]. As discussed in section IIB, we focus for simplicity on a system where at the transition point, $\vec{r} = \vec{r}^*$ and $\varepsilon = \gamma = 0$, there are no mechanisms, $N_{\text{mech}}^* = 0$. Note however, that this does not preclude the existence of mechanisms away from the transition point.

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

$$ \vec{x}'(0) = \vec{z} \quad \text{and} \quad \left. \frac{d^k \ell_i(\vec{x}(u), \gamma)}{du^k} \right|_{u=0} = 0 \qquad (14) $$

for all integers $k \le K$, while this condition is not fulfilled for $K+1$. I.e. there is no $K+1$ times differentiable function $\vec{y} : u \mapsto \vec{y}(u), \mathbb{R} \to \mathbb{R}^{N_{\text{dof}}}$ with $\vec{y}(0) = \vec{r}$ such that for all $i = 1, \dots, N_s$:

$$ \vec{y}'(0) = \vec{z} \quad \text{and} \quad \left. \frac{d^k \ell_i(\vec{y}(u), \gamma)}{du^k} \right|_{u=0} = 0 \qquad (15) $$

for all integers $k \le K+1$. We stress that due to the last condition, a mechanism is never a Kth-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 Kth-order zero modes by $N_{1st}, N_{2nd}, \dots$. Specifically, at the transition point, $\vec{r} = \vec{r}^*$ and $\varepsilon = \gamma = 0$, we denote those numbers by $N_{1st}^*, N_{2nd}^*, \dots$.

While there is generally a finite number of first-order zero modes at the transition, $N_{1st}^*$ (see Fig. 1 for examples), we expect the occurrence of Kth-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_{Kth}^* = 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 \to 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.

### A. 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. \quad (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 (\Delta\ell_i + d_i\varepsilon)^2. \quad (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}^*$ and for small $\gamma$. Denoting the components of $\vec{r}$ by $r_n$ with $n = 1, \dots, N_{\text{dof}}$, we can write to second order in $\Delta r_n := r_n - r_n^*$ and $\gamma$:

$$ w_i \Delta \ell_i = C_{in}^{(0)} \Delta r_n + \frac{1}{2} M_{imn}^{(0)} \Delta r_m \Delta r_n + B_i^{(1)} \gamma + C_{in}^{(1)} \Delta r_n \gamma + \frac{1}{2} B_i^{(2)} \gamma^2. \quad (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_{\text{dof}}$). On the right-hand side, $C_{in}^{(0)} \equiv C_{in} := w_i(\partial\ell_i/\partial r_n)$ is a generalized compatibility matrix and the matrix $M_{imn}^{(0)} \equiv M_{imn} := w_i(\partial^2\ell_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_i(\partial\ell_i/\partial\gamma)$, $C_{in}^{(1)} := w_i(\partial^2\ell_i/\partial r_n\partial\gamma)$, and $B_i^{(2)} := w_i(\partial^2\ell_i/\partial\gamma^2)$. All derivatives are evaluated at the transition point, $\vec{r} = \vec{r}^*$, and $\gamma = 0$. Higher-order terms in $\Delta\vec{r}$ and $\gamma$ are not 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$ (section IIC).

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}^{(0)}$. To this end, we first perform a singular-value decomposition (SVD) on it:

$$ C_{in}^{(0)} = \sum_{p=1}^{N_s} U_{ip} s_p V_{np}, \quad (19) $$

where $U_{ip}$ and $V_{np}$ are orthogonal square matrices of dimensions $N_s \times N_s$ and $N_{\text{dof}} \times N_{\text{dof}}$, respectively. The $N_s$ singular values $s_p$ are sorted in decreasing order: $s_1 \ge \dots \ge 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_{\text{dof}}$ are infinitesimal zero modes.

Note that because $s_{N_s} = 0$, Eq. (19) still holds when replacing $U_{iN_s}$ by $-U_{iN_s}$ for all $i = 1, \dots, N_s$. The same is also true for the orthogonality of $U_{ip}$. In other words, the sign of the vector $(U_{iN_s})$ is not uniquely defined by Eq. (19). In the following, we choose this sign such that

$$ \sum_{i=1}^{N_s} U_{iN_s} d_i w_i \ge 0 \quad (20) $$

for later convenience.

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}^{(0)}$, respectively:

$$ \Delta \tilde{\ell}_p := \sum_i U_{ip} w_i \Delta \ell_i \quad \text{for } p = 1, \dots, N_s, \quad (21) $$

$$ \Delta \tilde{r}_p := V_{np} \Delta r_n \quad \text{for } p = 1, \dots, N_{\text{dof}}. \quad (22) $$

As a consequence of Eq. (19), $s_{N_s} = 0$, and Eq. (22), all $\Delta \tilde{r}_p$ with $p = N_s, \dots, N_{\text{dof}}$ are infinitesimal zero modes.

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

$$ \Delta \tilde{\ell}_p = s_p \Delta \tilde{r}_p + \frac{1}{2} \tilde{M}_{pqr} \Delta \tilde{r}_q \Delta \tilde{r}_r + \tilde{B}_p^{(1)} \gamma + \tilde{C}_{pq}^{(1)} \Delta \tilde{r}_q \gamma + \frac{1}{2} \tilde{B}_p^{(2)} \gamma^2, \quad (23) $$

where we have also transformed $M$, $C^{(1)}$, and $B^{(1/2)}$ into the eigenbases of $C^{(0)}$: $\tilde{M}_{pqr} := U_{ip}M_{imn}V_{mq}V_{nr}$, $\tilde{C}_{pq}^{(1)} := U_{ip}C_{in}^{(1)}V_{nq}$, 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^{(0)}$:

$$ 
\begin{aligned}
E &= \frac{E_0}{2} \sum_{i=1}^{N_s} (w_i \Delta \ell_i + d_i w_i \varepsilon)^2 \\
&= \frac{E_0}{2} \sum_{i,j=1}^{N_s} (w_i \Delta \ell_i + d_i w_i \varepsilon) \delta_{ij} (w_j \Delta \ell_j + d_j w_j \varepsilon) \\
&= \frac{E_0}{2} \sum_{i,j,p=1}^{N_s} (w_i \Delta \ell_i + d_i w_i \varepsilon) U_{ip} U_{jp} (w_j \Delta \ell_j + d_j w_j \varepsilon) \\
&= \frac{E_0}{2} \sum_{p=1}^{N_s} (\Delta \tilde{\ell}_p + \tilde{w}_p \varepsilon)^2.
\end{aligned} 
\quad (24) $$

In the second step of Eq. (24), 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. (21) and introduced $\tilde{w}_p := \sum_{i=1}^{N_s} U_{ip}d_i w_i$. Note that because of Eq. (20), we have $\tilde{w}_{N_s} \ge 0$.

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

$$ 
\begin{aligned}
\Lambda_p &:= \Delta \tilde{\ell}_p + \tilde{w}_p \varepsilon \\
&= s_p \Delta \tilde{r}_p + \frac{1}{2} \tilde{M}_{pqr} \Delta \tilde{r}_q \Delta \tilde{r}_r \\
&\quad + \tilde{w}_p \varepsilon + \tilde{B}_p^{(1)} \gamma + \tilde{C}_{pq}^{(1)} \Delta \tilde{r}_q \gamma + \frac{1}{2} \tilde{B}_p^{(2)} \gamma^2.
\end{aligned} 
\quad (25) $$

In the second step, we used Eq. (23). Inserting the definition of $\Lambda_p$ into the energy, Eq. (24), we obtain:

$$ E = \frac{E_0}{2} \sum_{p=1}^{N_s} \Lambda_p^2. \quad (26) $$

Thus, each $\Lambda_p$ indicates how much the corresponding left eigenmode $p$ of $C^{(0)}$ contributes to the system's energy.

## B. Energy minimum

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

At $\varepsilon = \gamma = 0$, the system is at its transition point $\vec{r} = \vec{r}^*$ (section II B), which implies $\Delta \tilde{r}_q = 0$ and thus, using Eq. (25), $\Lambda_p = 0$. However, as $\varepsilon$ and $\gamma$ become nonzero, also the $\Delta \tilde{r}_q$ and $\Lambda_p$ generally become nonzero. We provide the detailed arguments for the scaling of $\Lambda_p$ and $\Delta \tilde{r}_q$ with $\varepsilon$ and $\gamma$ in appendix B1, and provide here a brief sketch of these arguments. Minimizing the energy in Eq. (26) 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. (B10) in appendix B1):

$$ \Delta \tilde{r}_p = -\frac{1}{s_p} \left( \frac{1}{2} \sum_{q,r=N_s}^{N_{\text{dof}}} \tilde{M}_{pqr} \Delta \tilde{r}_q \Delta \tilde{r}_r + \tilde{w}_p \varepsilon + \frac{1}{2} \tilde{B}_p^{(2)} \gamma^2 \right). \quad (27) $$

However, $\Lambda_{N_s}$ cannot be set to zero in the same way, i.e. with $\Delta \tilde{r}_{N_s}$ given by Eq. (27), 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. (25) we have:

$$ \Lambda_{N_s} = \frac{1}{2} \sum_{q,r=N_s}^{N_{\text{dof}}} \tilde{M}_{N_s q r} \Delta \tilde{r}_q \Delta \tilde{r}_r + \tilde{w}_{N_s} \varepsilon + \frac{1}{2} \tilde{B}_{N_s}^{(2)} \gamma^2. \quad (28) $$

In the sum, the contributions by the terms with $q < N_s$ or $r < N_s$ can be neglected to lowest order, since they only generate higher order terms in $\varepsilon$ and $\gamma$ (Eq. (27) and appendix B1).

We know that for $\gamma = 0$ and $\varepsilon > 0$, where the system is in the athermal rigid regime, the shear modulus is positive, and thus 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 q r}$ (restricted to $q, r = N_s, \dots, N_{\text{dof}}$) is positive semi-definite; otherwise $\Lambda_{N_s}$ could become zero.

To further discuss the energy minimum, we introduce

$$ b_{\varepsilon} := \frac{\tilde{B}_{N_s}^{(2)}}{2\tilde{w}_{N_s}} \quad \text{if } \tilde{B}_p^{(1)} = \tilde{C}_{pq}^{(1)} = 0 \text{ for all } p, q, \quad (29) $$

so that $\Lambda_{N_s}$ becomes

$$ \Lambda_{N_s} = \frac{1}{2} \sum_{q,r=N_s}^{N_{\text{dof}}} \tilde{M}_{N_s q r} \Delta \tilde{r}_q \Delta \tilde{r}_r + \tilde{w}_{N_s} [\varepsilon + b_{\varepsilon} \gamma^2]. \quad (30) $$

Because $\tilde{M}_{N_s q r}$ is positive semi-definite, we obtain that for $\varepsilon + b_{\varepsilon} \gamma^2 > 0$, the energy is minimal if $\Delta \tilde{r}_q = 0$ for all $N_s \le n \le N_{\text{dof}}$, and it has the value

$$ E = \frac{E_0 \tilde{w}_{N_s}^2}{2} [\varepsilon + b_{\varepsilon} \gamma^2]^2. \quad (31) $$

Meanwhile, for $\varepsilon + b_{\varepsilon} \gamma^2 \le 0$, the minimal energy becomes zero since $\Lambda_{N_s}$ can become zero. Hence, $\varepsilon + b_{\varepsilon} \gamma^2 > 0$ corresponds to the athermal rigid regime and $\varepsilon + b_{\varepsilon} \gamma^2 < 0$ to the athermal floppy regime.

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

Note that the matrix $\tilde{M}_{N_s q r}$ is closely linked to the Hessian of the system and the emergence of rigidity in the athermal limit (appendix D).
### C. 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^*)^2 (1 + a_\ell^2)} [-\Delta L_0 + b\gamma^2]^2, \quad (32)$$

where $\Delta L_0 = -dL_0^*\epsilon$, and $L_0^*, 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. (31) and (32), we find:

$$\tilde{w}_{N_s} = d \left[ \frac{N_s}{1 + a_\ell^2} \right]^{1/2} \quad (33)$$

$$b_{\epsilon} = \frac{b}{dL_{0}^{*}}. \quad (34)$$

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)$.

### D. Elastic properties

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

$$t_E = \kappa_E [\epsilon + b_\epsilon \gamma^2] \quad (35)$$

$$G_E = 2Db_\epsilon\kappa_E [\epsilon + 3b_\epsilon\gamma^2] \quad (36)$$

$$= 2Db_\epsilon t \left( 1 + \frac{2b_\epsilon \gamma^2}{|\epsilon + b_\epsilon \gamma^2|} \right), \quad (37)$$

where we defined

$$\kappa_E := \frac{E_0 \tilde{w}_{N_s}^2}{D\mathcal{V}^*}. \quad (38)$$

The index $E$ indicates purely energetic rigidity.

## IV. STIFF-SPRING LIMIT ($K_i \to \infty$)

In the limit of infinitely stiff springs but at finite temperature $T$, elasticity is created purely by entropic effects. For $\epsilon+b_\epsilon\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 $\epsilon + b_\epsilon\gamma^2 \le 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 $\epsilon$ and $\gamma$.

## A. Configurational phase space volume

In appendix E we explicitly take the limit of the partition sum $Z$ for $K_i \to \infty$ to obtain an expression for the configurational phase space volume $\Omega$ [29]. Up to a constant factor, it is given by:

$$\Omega \sim \int \left( \prod_{n=1}^{N_{\text{dof}}} dR_n \right) \prod_{i=1}^{N_s} \delta(L_i - L_{0i}), \quad (39)$$

where $R_n$ for $n = 1, \dots, N_{\text{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_{\text{dof}}} dr_n \right) \prod_{i=1}^{N_s} \delta(\ell_i - e^{-d_i \epsilon}). \quad (40)$$

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

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

$$\Omega \sim \int \left( \prod_{n=1}^{N_{\text{dof}}} d\Delta r_n \right) \prod_{i=1}^{N_s} \delta(w_i \Delta \ell_i + d_i w_i \epsilon). \quad (41)$$

Here, we have included a constant factor of $w_i$ in each of the Dirac deltas, and we again ignored prefactors of $e^\epsilon$ in $\Omega$.

To better understand which set of points $\vec{r}$ is selected by the Dirac deltas in Eq. (41), 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{align}
\Omega &\sim \int \left( \prod_{q=1}^{N_{\text{dof}}} d\Delta \tilde{r}_q \right) \prod_{p=1}^{N_s} \delta(\Delta \tilde{\ell}_p + \tilde{w}_p \epsilon) \nonumber \\
&= \int \left( \prod_{q=1}^{N_{\text{dof}}} d\Delta \tilde{r}_q \right) \prod_{p=1}^{N_s} \delta(\Lambda_p). \tag{42}
\end{align}$$

In the second step, we have used the definition of the $\Lambda_p$, Eq. (25). 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. (25).

We now discuss which values of the dofs, $\Delta \tilde{r}_q$, satisfy the Dirac delta conditions in Eq. (42). Here, we focus again on the case where all $\tilde{B}_p^{(1)} = \tilde{C}_{pq}^{(1)} = 0$. The general case of non-zero $\tilde{B}_p^{(1)}$ and $\tilde{C}_{pq}^{(1)}$ is treated in appendix C. 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}_q\Delta\tilde{r}_r$ with $q < N_s$ or $r < N_s$ can be neglected (appendix B 2). Thus, we can split the integrations in Eq. (42) as follows:

$$ \Omega \sim \int \left( \prod_{q=N_s}^{N_{\text{dof}}} d\Delta\tilde{r}_q \right) \delta(\Lambda_{N_s}) \prod_{p=1}^{N_s-1} \int d\Delta\tilde{r}_p \delta(\Lambda_p). \quad (43) $$

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), \quad (44) $$

where

$$ \Omega_{N_s}(\varepsilon, \gamma) := \int \left( \prod_{q=N_s}^{N_{\text{dof}}} d\Delta\tilde{r}_q \right) \delta(\Lambda_{N_s}) \quad (45) $$

with $\Lambda_{N_s}$ given by Eq. (30):

$$ \Lambda_{N_s} = \frac{1}{2} \sum_{q,r=N_s}^{N_{\text{dof}}} \tilde{M}_{N_s qr} \Delta\tilde{r}_q \Delta\tilde{r}_r + \tilde{w}_{N_s} [\varepsilon + b_{\varepsilon} \gamma^2]. \quad (46) $$

To evaluate $\Omega_{N_s}(\varepsilon, \gamma)$, we first consider the symmetric matrix $\tilde{M}_{N_s qr}$ restricted to $q, r = N_s, \dots, N_{\text{dof}}$, which appears in $\Lambda_{N_s}$. We diagonalize this $(N_{\text{dof}} - N_s + 1) \times (N_{\text{dof}} - N_s + 1)$ sub matrix:

$$ \tilde{M}_{N_s qr} = \sum_{k=1}^{N_{\text{dof}}-N_s+1} \mu_k v_q^k v_r^k, \quad (47) $$

where $\mu_k$ are the eigenvalues and $v_q^k$ are the associated orthonormal eigenvectors. We sort the eigenvalues in decreasing order, $\mu_1 \ge \dots \ge \mu_{N_{\text{dof}}-N_s+1} \ge 0$. All eigenvalues are non-negative, because the matrix $\tilde{M}_{N_s qr}$ is positive semi-definite.

With

$$ \Delta q_k := \sum_{q=N_s}^{N_{\text{dof}}} v_q^k \Delta \tilde{r}_q, \quad (48) $$

the expression for $\Lambda_{N_s}$, Eq. (46), transforms into:

$$ \Lambda_{N_s} = \frac{1}{2} \sum_{k=1}^{N_{\text{dof}}-N_s+1} \mu_k \Delta q_k^2 + \tilde{w}_{N_s} [\varepsilon + b_{\varepsilon} \gamma^2]. \quad (49) $$

We now substitute $\Lambda_{N_s}$ into $\Omega_{N_s}$, Eq. (45):

$$ \Omega_{N_s} = \int \left( \prod_{k=1}^{N_{\text{dof}}-N_s+1} d\Delta q_k \right) \times \delta \left( \frac{1}{2} \sum_{k=1}^{N_{\text{dof}}-N_s+1} \mu_k \Delta q_k^2 + \tilde{w}_{N_s} [\varepsilon + b_{\varepsilon} \gamma^2] \right). \quad (50) $$

Because the $v_q^k$ represent a rotation, the differential transforms without prefactor.

## B. Eigenmodes of $\tilde{M}_{N_s qr}$

To evaluate the integral in Eq. (50), we first discuss how the matrix $\tilde{M}_{N_s qr}$ relates to the different kinds of zero modes in the system (see also section IID).

At the transition point, there are $N_{\text{dof}} - N_s + 1$ linearly independent infinitesimal zero modes. We consider here the zero modes $\tilde{z}_k$ with $k = 1, \dots, N_{\text{dof}} - N_s + 1$ that have the components $z_{k,m} := \sum_{q=N_s}^{N_{\text{dof}}} V_{mq} v_q^k$. Insertion in Eq. (11) with Eq. (19) confirms that they are infinitesimal zero modes. According to Eqs. (22) and (48) these modes contribute to $\Delta\tilde{r} = \tilde{r} - \tilde{r}^*$ with the respective amplitudes $\Delta q_k$.

Those of the $\Delta q_k$ modes that correspond to zero eigenvalues, $\mu_k = 0$, of the matrix $\tilde{M}_{N_s qr}$ do not contribute to the argument of the Dirac delta in Eq. (50). There are three possibilities for the corresponding infinitesimal zero modes $\tilde{z}_k$:

* $\tilde{z}_k$ could be one of the D global translations. In this case, the integral in Eq. (50) leads to a constant factor.

* $\tilde{z}_k$ could belong to a mechanism existing at the transition point, but here we focus on systems without such mechanisms (section II B). 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.

* $\tilde{z}_k$ belongs to a Kth-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 $\tilde{z}_k$ with $\mu_k = 0$ can be a first-order zero mode [30].

Conversely, all of the infinitesimal zero modes $\tilde{z}_k$ with $\mu_k > 0$ do contribute to the argument of the Dirac delta in Eq. (50). All of them are first-order zero modes. This is because these modes affect $\Delta\tilde{L}_{N_s}$ to second order [31]. 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 section IID that we denote the number of first-order zero modes at the transition point by $N_{1st}^*$.

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 in Fig. 1, where we show a few examples for $N_{1st}^* = 1, \dots, 3$. As soon as we apply isotropic compression, $\varepsilon < 0$ (while formally keeping $\gamma = 0$), exactly $N_{1st}^* - 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$. Together with this disappearing infinitesimal zero
<figure><img src="image_2.png" /><figcaption>FIG. 1. Zero mode counting for four example systems, each at the transition point ($\epsilon = 0$) and under compression ($\epsilon < 0$). In each case, we respectively indicate the number of spatial dimensions $D$, the number of degrees of freedom $N_{\text{dof}}$, the number of springs $N_s$, the number of states of self-stress $N_{\text{SSS}}$, the number of infinitesimal zero modes $N_0$, the number of 1st-order zero modes $N_{\text{1st}}$ (definition in section II D), and the number of mechanisms $N_{\text{mech}}$. At the transition, $\epsilon = 0$, a finite number $N_{\text{1st}}^*$ of 1st-order zero modes exists. Upon compression, $\epsilon < 0$, exactly $N_{\text{1st}}^* - 1$ of these modes turn into mechanisms. The remaining mode ceases to be an infinitesimal zero mode for $\epsilon < 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.</figcaption></figure>

mode, the SSS that existed at $\epsilon = 0$ disappears for $\epsilon < 0$ (compare Eq. (1)).

### C. Scaling of 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. (50) only leads to factors that are to lowest order finite constants, and we obtain:

$$ \Omega_{N_s} \sim \int \left( \prod_{k=1}^{N_{\text{1st}}^*} d\Delta q_k \right) \times \delta \left( \frac{1}{2} \sum_{k=1}^{N_{\text{1st}}^*} \mu_k \Delta q_k^2 + \tilde{w}_{N_s} [\varepsilon + b_{\varepsilon} \gamma^2] \right). \quad (51) $$

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

$$ \sum_{k=1}^{N_{\text{1st}}^*} \frac{\Delta q_k^2}{\rho_k^2} = 1 \quad (52) $$

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 (section IV B), as soon as compression is applied, $\varepsilon < 0$, there are $N_{\text{1st}}^* - 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}} dS \left| \frac{\partial \Lambda_{N_s}}{\partial \Delta \vec{q}} \right|^{-1}, \quad (53) $$

where $dS$ is an integration element of the hyper-ellipsoid surface $\partial\mathcal{E}$, and

$$ \begin{aligned} \left| \frac{\partial \Lambda_{N_s}}{\partial \Delta \vec{q}} \right|^2 &= \sum_{k=1}^{N_{\text{1st}}^*} \left( \frac{\partial \Lambda_{N_s}}{\partial \Delta q_k} \right)^2 \\ &= \sum_{k=1}^{N_{\text{1st}}^*} \mu_k^2 \Delta q_k^2. \end{aligned} \quad (54) $$

With $\mu_{\text{max}} := \mu_1$ and $\mu_{\text{min}} := \mu_{\text{1st}}^* > 0$, we thus have on the ellipsoid surface, using Eq. (52):

$$ -\mu_{\min} [\varepsilon + b_{\varepsilon} \gamma^2] \le \frac{1}{2\tilde{w}_{N_s}} \left| \frac{\partial \Lambda_{N_s}}{\partial \Delta \vec{q}} \right|^2 \le -\mu_{\max} [\varepsilon + b_{\varepsilon} \gamma^2]. \quad (55) $$

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

$$ \Omega_{N_s}(\varepsilon, \gamma) \sim \left[ - (\varepsilon + b_{\varepsilon} \gamma^2) \right]^{(N_{\text{1st}}^* - 2)/2}. \quad (56) $$

Hence, $\Omega \sim \Omega_{N_s}$ scales as a power law with the distance to the transition point as long as $N_{\text{1st}}^* > 2$. In appendix F, we prove this result in a more general context, where we show that the exponent in Eq. (56) can change if we allow for Kth-order zero modes with $K > 1$ at the transition.

In appendix C, we discuss the more general case where the $\tilde{B}_p^{(1)}$ and $\tilde{C}_{pq}^{(1)}$ do not vanish. The definition of $b_{\varepsilon}$ changes like in the athermal case, while the result Eq. (56) remains the same up to an additional term $B_{N_s}^{(1)}\gamma$, which can be removed through redefining $\gamma$.
### D. 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{aligned}
F_S &= -\frac{1}{\beta} \log \Omega \\
&= -\frac{N_{1\text{st}}^* - 2}{2\beta} \log \left[ -(\varepsilon + b_{\varepsilon} \gamma^2) \right],
\end{aligned} 
\quad (57) $$

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) / DV$ and the shear modulus $G_S = (\partial^2 F_S / \partial \gamma^2) / V$, to lowest order:

$$ t_S = -\frac{\kappa_S T}{\varepsilon + b_{\varepsilon} \gamma^2} \quad (58) $$

$$ G_S = -2Db_{\varepsilon}\kappa_S T \frac{\varepsilon - b_{\varepsilon}\gamma^2}{(\varepsilon + b_{\varepsilon}\gamma^2)^2} \quad (59) $$

$$ = 2Db_{\varepsilon}t \left( 1 + \frac{2b_{\varepsilon}\gamma^2}{|\varepsilon + b_{\varepsilon}\gamma^2|} \right) \quad (60) $$

with

$$ \kappa_S := \frac{k_B(N_{1\text{st}}^* - 2)}{2DV^*}. \quad (61) $$

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

### A. 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_{\text{dof}}} dR_n \right) e^{-\beta E}. \quad (62) $$

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

$$ Z \sim \int \left( \prod_{q=1}^{N_{\text{dof}}} d\Delta \tilde{r}_q \right) e^{-\beta E}, \quad (63) $$

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. (26):

$$ E = \frac{E_0}{2} \sum_{p=1}^{N_s} \Lambda_p^2. \quad (64) $$

Here we discuss the case where all $\tilde{B}_p^{(1)} = \tilde{C}_{pq}^{(1)} = 0$, while the more general case of arbitrary $\tilde{B}_p^{(1)}$ and $\tilde{C}_{pq}^{(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 \le N_s$ all the $\Delta \tilde{r}_p$ with $p < N_s$ can be neglected in the terms $\sim \Delta \tilde{r}_q \Delta \tilde{r}_p$ (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. (63) can be rearranged as follows:

$$ Z \sim \int \left( \prod_{q=N_s}^{N_{\text{dof}}} d\Delta \tilde{r}_q \right) e^{-\frac{\beta E_0 \Lambda_{N_s}^2}{2}} \prod_{p=1}^{N_s-1} \int d\Delta \tilde{r}_p e^{-\frac{\beta E_0 \Lambda_p^2}{2}}. \quad (65) $$

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. (62), which arises from the kinetic energy of the nodes [32]. This prefactor is given by $\lambda^{-DN_{\text{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_{q=N_s}^{N_{\text{dof}}} d\Delta \tilde{r}_q \right) e^{-\frac{\beta E_0 \Lambda_{N_s}^2}{2}}. \quad (66) $$

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}$.

### B. Saddle point approximation

To simplify the remaining integral $Z_{N_s}$ in Eq. (66), we introduce in its integrand an additional integral over some variable $\varepsilon_S$ that trivially evaluates to one:

$$ Z_{N_s} = \int \left( \prod_{q=N_s}^{N_{\text{dof}}} d\Delta \tilde{r}_q \right) e^{-\frac{\beta E_0 \Lambda_{N_s}^2}{2}} \times \tilde{w}_{N_s} \int_{-\infty}^{\infty} d\varepsilon_S \delta \left( \Lambda_{N_s} + \tilde{w}_{N_s} [\varepsilon_S - \varepsilon] \right). \quad (67) $$
After exchanging the order of integration and substituting $\Omega_{N_s}$, defined in Eq. (45), we have:

$$Z_{N_s} = \tilde{w}_{N_s} \int_{-\infty}^{\infty} d\epsilon_S \, \Omega_{N_s}(\epsilon_S, \gamma) e^{-\beta E_{N_s}(\epsilon - \epsilon_S)}, \quad (68)$$

with:

$$\begin{aligned} E_{N_s}(\epsilon - \epsilon_S) &:= \frac{E_0 \tilde{w}_{N_s}^2}{2} (\epsilon - \epsilon_S)^2 \\ &= \frac{D\mathcal{V}^* \kappa_E}{2} (\epsilon - \epsilon_S)^2. \end{aligned} \quad (69)$$

Here, we used the definition for $\kappa_E$ in Eq. (38). From the discussion in the previous section, we know that for $\epsilon_S + b_\epsilon \gamma^2 > 0$, no configurations are possible, i.e. $\Omega_{N_s}(\epsilon_S, \gamma) = 0$. Meanwhile, for $\epsilon_S + b_\epsilon \gamma^2 \le 0$, we have from Eq. (56) that:

$$\Omega_{N_s}(\epsilon_S, \gamma) \sim \left[ -(\epsilon_S + b_\epsilon \gamma^2) \right]^{(N_{1st}^* - 2)/2}. \quad (70)$$

These transformations identify $\epsilon_S$ as the entropic strain discussed in the companion paper [26].

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

$$Z_{N_s} = \tilde{w}_{N_s} \int_{-\infty}^{\infty} d\epsilon_S e^{-\beta \bar{F}}, \quad (71)$$

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

$$\begin{aligned} \bar{F}(\epsilon_S; \epsilon, \gamma) &:= E_{N_s}(\epsilon - \epsilon_S) - \frac{1}{\beta} \log \Omega_{N_s}(\epsilon_S, \gamma) \\ &= D\mathcal{V}^* \left( \frac{\kappa_E}{2} (\epsilon - \epsilon_S)^2 - \kappa_S T \log \left[ -(\epsilon_S + b_\epsilon \gamma^2) \right] \right) \end{aligned} \quad (72)$$

On the second line, we used the definition of $\kappa_S$, Eq. (61), and ignored a constant offset.

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

$$Z_{N_s} \sim \int_{-\infty}^{\infty} d\epsilon_S e^{-\beta \bar{F}} = e^{-\beta \bar{F}(\hat{\epsilon}_S)} \sqrt{\frac{2\pi}{\beta \bar{F}''(\hat{\epsilon}_S)}}, \quad (73)$$

where $\bar{F}'' := \partial^2 F / \partial \epsilon_S^2$:

$$\bar{F}''(\hat{\epsilon}_S) = D\mathcal{V}^* \left( \kappa_E + \frac{\kappa_S T}{[\hat{\epsilon}_S + b_\epsilon \gamma^2]^2} \right). \quad (74)$$

To find the minimum $\hat{\epsilon}_S$, we transform $\partial \bar{F} / \partial \epsilon_S = 0$, and obtain:

$$\hat{\epsilon}_S + b_\epsilon \gamma^2 = \frac{\epsilon + b_\epsilon \gamma^2}{2} - \sqrt{1 + \theta} \left| \frac{\epsilon + b_\epsilon \gamma^2}{2} \right|, \quad (75)$$

with

$$\begin{aligned} \theta &:= \frac{4\kappa_S T}{\kappa_E(\epsilon + b_\epsilon \gamma^2)^2} \\ &= \frac{2(N_{1st}^* - 2)k_B T}{E_0 \tilde{w}_{N_s}^2 (\epsilon + b_\epsilon \gamma^2)^2}. \end{aligned} \quad (76)$$

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

### C. 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{\epsilon}_S) + \frac{1}{2\beta} \log \bar{F}''(\hat{\epsilon}_S). \quad (77)$$

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

$$t = \kappa_E(\epsilon - \hat{\epsilon}_S) \left[ 1 + \frac{1}{2(N_{1st}^* - 2)} \frac{\theta}{1 + \theta} \right] \quad (78)$$

and

$$\begin{aligned} G = 2Db_\epsilon \kappa_E (\epsilon - \hat{\epsilon}_S) & \left[ 1 + \frac{2b_\epsilon \gamma^2}{|\epsilon + b_\epsilon \gamma^2|\sqrt{1+\theta}} \right. \\ & \left. + \frac{1}{2(N_{1st}^* - 2)} \frac{\theta}{1+\theta} \left( 1 + \frac{2b_\epsilon \gamma^2}{|\epsilon + b_\epsilon \gamma^2|\sqrt{1+\theta}} \right. \right. \\ & \left. \left. - \frac{4b_\epsilon \gamma^2}{(\epsilon + b_\epsilon \gamma^2)(1+\theta)} \right) \right]. \end{aligned} \quad (79)$$

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

$$t = \kappa_E(\epsilon - \hat{\epsilon}_S) \quad (80)$$

$$G = 2Db_\epsilon \kappa_E (\epsilon - \hat{\epsilon}_S) \left[ 1 + \frac{2b_\epsilon \gamma^2}{|\epsilon + b_\epsilon \gamma^2|\sqrt{1+\theta}} \right] \quad (81)$$

$$= 2Db_\epsilon t \left[ 1 + \frac{2b_\epsilon \gamma^2}{|\epsilon + b_\epsilon \gamma^2|\sqrt{1+\theta}} \right]. \quad (82)$$

In these equations, the value of $\hat{\epsilon}_S$ can be obtained using Eq. (75).
## D. Limits

1. Small temperature, $\theta \ll 1$

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

$$ 
\begin{aligned}
\hat{\epsilon}_S + b_\epsilon \gamma^2 &= \frac{\epsilon + b_\epsilon \gamma^2}{2} - \left(1 + \frac{\theta}{2}\right) \left|\frac{\epsilon + b_\epsilon \gamma^2}{2}\right| \\
&= (\epsilon + b_\epsilon \gamma^2) \times \begin{cases} 1 + \theta/4 & \text{for } \epsilon + b_\epsilon \gamma^2 < 0, \\ -\theta/4 & \text{for } \epsilon + b_\epsilon \gamma^2 > 0. \end{cases} \tag{83}
\end{aligned}
$$

Thus, to lowest order in $\theta$:

$$ \epsilon - \hat{\epsilon}_S = (\epsilon + b_\epsilon \gamma^2) \times \begin{cases} -\theta/4 & \text{for } \epsilon + b_\epsilon \gamma^2 < 0, \\ 1 & \text{for } \epsilon + b_\epsilon \gamma^2 > 0, \end{cases} \tag{84} $$

and we recover the purely entropic and energetic tensions and shear moduli, respectively (compare Eqs. (35), (37), (58), and (60)):

$$ t = \begin{cases} t_S & \text{for } \epsilon + b_\epsilon \gamma^2 < 0, \\ t_E & \text{for } \epsilon + b_\epsilon \gamma^2 > 0, \end{cases} \tag{85} $$

$$ G = \begin{cases} G_S & \text{for } \epsilon + b_\epsilon \gamma^2 < 0, \\ G_E & \text{for } \epsilon + b_\epsilon \gamma^2 > 0. \end{cases} \tag{86} $$

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

2. Small strain, $\theta \gg 1$

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

$$ \hat{\epsilon}_S + b_\epsilon \gamma^2 = -\sqrt{\frac{\kappa_S T}{\kappa_E}} + \frac{\epsilon + b_\epsilon \gamma^2}{2} \tag{87} $$

And thus:

$$ \epsilon - \hat{\epsilon}_S = \sqrt{\frac{\kappa_S T}{\kappa_E}} + \frac{\epsilon + b_\epsilon \gamma^2}{2}. \tag{88} $$

We then find:

$$ t = \frac{2N_{1st}^* - 3}{2N_{1st}^* - 4} \left( \sqrt{\kappa_E \kappa_S T} + \kappa_E \frac{\epsilon + b_\epsilon \gamma^2}{2} \right), \tag{89} $$

and for the shear modulus:

$$ G = Db_\epsilon \left[ 2\sqrt{\kappa_E \kappa_S T} + \kappa_E (\epsilon + b_\epsilon \gamma^2) \right] \times \left[ \frac{2N_{1st}^* - 3}{2N_{1st}^* - 4} \left( 1 + b_\epsilon \gamma^2 \sqrt{\frac{\kappa_E}{\kappa_S T}} \right) - \frac{b_\epsilon \gamma^2}{2N_{1st}^* - 4} \frac{\kappa_E (\epsilon + b_\epsilon \gamma^2)}{\kappa_S T} \right]. \tag{90} $$

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

$$ t = \sqrt{\kappa_E \kappa_S T} + \kappa_E \frac{\epsilon + b_\epsilon \gamma^2}{2} \tag{91} $$

$$ G = 2Db_\epsilon t \left( 1 + b_\epsilon \gamma^2 \sqrt{\frac{\kappa_E}{\kappa_S T}} \right). \tag{92} $$

Note that we assumed in our derivation for the general case that $k_B T \lesssim E_0[\epsilon + b_\epsilon \gamma^2]^2$. Using the definitions of $\theta$ in Eq. (76), of $\kappa_E$ in Eq. (38), and of $\kappa_S$ in Eq. (61), this condition corresponds to $\theta \lesssim 2(N_{1st}^* - 2)/\hat{w}_{N_\epsilon}^2$. Thus, this limit corresponds to $1 \ll \theta \lesssim 2(N_{1st}^* - 2)/\hat{w}_{N_\epsilon}^2$, which can be fulfilled for $N_{1st}^* \gg 1$. 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, $\epsilon + b_\epsilon \gamma^2 \to 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 $\epsilon$, 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_\epsilon$. These parameters depend on the microscopic structure of the system through Eqs. (38), (61), and (29) (or, more generally, Eq. (C10)), respectively.

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 $\epsilon + b_\epsilon \gamma^2 \le 0$ and rigid for $\epsilon + b_\epsilon \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_{1st}^*$ (section IIID). This ultimately leads to a free energy that scales logarithmically with the negative combined strain $-(\epsilon + b_\epsilon \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 $\epsilon < 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 as we have shown before [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 section V D, where up to a prefactor, isotropic strain $\epsilon$ 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. (61) predicts that $\kappa_S$ essentially only depends on the number of first-order zero modes $N_{1st}^*$. Thus, if for large systems $N_{1st}^*$ scales to dominating order linearly with $\Delta z = -4(N_{dof} - N_s)/N_{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. (80), 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. (82), the scalings of $b_\epsilon \kappa_S$ and $b_\epsilon \kappa_E$ are the relevant ones, and thus the scaling of $b_\epsilon \sim b$ needs to be taken into account. Based on earlier results for the scaling of $b$ [3, 4], we have $b_\epsilon \kappa_S \sim |\Delta z|^{-1}$ and $b_\epsilon \kappa_E \sim |\Delta z|^0$ for Delaunay networks, $b_\epsilon \kappa_S \sim |\Delta z|^0$ and $b_\epsilon \kappa_E \sim |\Delta z|^0$ for randomly-cut, packing-derived networks, and $b_\epsilon \kappa_S \sim |\Delta z|^{0.5}$ and $b_\epsilon \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_\epsilon \kappa_S \sim |\Delta z|^0$ and $b_\epsilon \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, $\epsilon_{ini}$, with connectivity $\Delta z$, which is generally linear to lowest order: $\epsilon_{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 $\epsilon$ by $\Delta z$. We leave a more thorough numerical test of all these predictons for the scaling with $\Delta z$, and a study of their consequences, for future work.

We introduced under-constrained systems as systems that have $N_{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.* It is well-known that at least one SSS forms at the transition [5]. Moreover, situations where more than a single SSS forms at the transition are likely very rare, in particular in disordered networks. Also, we do not expect the number of SSS that form at the transition to substantially affect our results.

* *There are no Kth-order zero modes at the transition, $N_{Kth}^* = 0$, for $K > 1$. We expect the occurrence of such modes to be rare. Moreover, some preliminary work indicates that allowing for $N_{Kth}^* > 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_{1st}^* > 2$. The marginal case of $N_{1st}^* = 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_{1st}^* = 2$ in future work.*
* For our discussion in appendix B, we needed the additional assumption in Eq. (B2). Intuitively, we believe that one can obtain our result in appendix B quite generally even without this assumption, and it will be interesting to test if this is true in future research.

* There are no mechanisms at the transition. The reason why we need this assumption is that with a degeneracy $\vec{r}^*$ of the transition point configurations, in principle, the Taylor expansion in Eq. (18) could depend on $\vec{r}^*$. However, while properties like the singular values, $s_p$ with $p < N_s$, of the compatibility matrix $C_{in}^{(0)}$ may depend on $\vec{r}^*$ [6, 34], 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}^*$. While this seems to be true intuitively, preliminary results suggest that stronger assumptions about the system than those listed in section II A 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, we expect that many 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.

In the future, it will be interesting to compare these analytical results to more numerical data [35]. Furthermore, as we have shown, all of the three parameters, $\kappa_E$, $\kappa_S$, and $b_\varepsilon$, are properties of the SSS that is created at the transition. Thus, better understanding the structure of this SSS will be important going forward. Recently, conformal invariance has been harnessed to derive scaling properties of rigidity percolation [36, 37]. It will be interesting to see whether similar approaches can be used to understand the structure of the SSS in under-constrained materials, at least for some classes of disordered networks. Finally, it will be interesting to understand how our work generalizes to systems that are constructed by adding to an under-constrained system a number of weak generalized springs to make it over-constrained [38, 39].

## ACKNOWLEDGMENTS

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 section II A, 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 $\mathbf{L}_i = (L_{i,x}, L_{i,y})$ as follows [26]:

$$ 
\begin{aligned}
L_{i,x}(\vec{R}, \mathcal{V}, \gamma) &= R_{b_i,x} - R_{a_i,x} + (q_{i,x} + \gamma q_{i,y})\sqrt{\mathcal{V}} \\
L_{i,y}(\vec{R}, \mathcal{V}, \gamma) &= R_{b_i,y} - R_{a_i,y} + q_{i,y}\sqrt{\mathcal{V}},
\end{aligned}
\quad (A1) $$

where $\mathbf{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) = |\mathbf{L}_i(\vec{R}, \mathcal{V}, \gamma)|$.

In Eq. (A1), the symbol $q_i = (q_{i,x}, q_{i,y}) \in \mathbb{Z}^D$ denotes the periodicity vector of spring $i$ [13, 40]. 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 $\mathbf{R}_{a_i}$ to $\mathbf{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 $\mathbf{R}_{a_i}$ to $\mathbf{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:

$$ 
\begin{aligned}
L_{i,x}(\vec{R}, \mathcal{V}, \gamma) &= R_{b_i,x} - R_{a_i,x} + q_{i,x}\sqrt{\mathcal{V}} \\
&\quad + \gamma(R_{b_i,y} - R_{a_i,y} + q_{i,y}\sqrt{\mathcal{V}}) \quad (A2) \\
L_{i,y}(\vec{R}, \mathcal{V}, \gamma) &= R_{b_i,y} - R_{a_i,y} + q_{i,y}\sqrt{\mathcal{V}}.
\end{aligned}
$$

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

To implement pure shear in a $D=3$-dimensional system, we define the components of the length vector $\mathbf{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:

$$ 
\begin{aligned}
L_{i,x}(\vec{R}, \mathcal{V}, \gamma) &= R_{b_i,x} - R_{a_i,x} + q_{i,x}e^{\gamma\sqrt{\mathcal{V}}} \\
L_{i,y}(\vec{R}, \mathcal{V}, \gamma) &= R_{b_i,y} - R_{a_i,y} + q_{i,y}e^{-\gamma\sqrt{\mathcal{V}}} \quad (A3) \\
L_{i,z}(\vec{R}, \mathcal{V}, \gamma) &= R_{b_i,z} - R_{a_i,z} + q_{i,z}\sqrt{\mathcal{V}}.
\end{aligned}
$$

Again, the spring length is then given by: $L_i(\vec{R}, \mathcal{V}, \gamma) = |\mathbf{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(\tilde{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. (A1)–(A3) with $d_i = 1$. Finally, a similar approach as in Eqs. (A1)–(A3) can be used to implement simple or pure shear in vertex models instead of simple spring networks [13, 40].

## Appendix B: Orders of magnitude of the $\Lambda_p$ and $\Delta\tilde{r}_q$

### 1. Athermal limit

We discuss the scaling of the $\Delta\tilde{r}_q$ 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 q r}$ for $q, r = N_s, \dots, N_{\text{dof}}$ (see also section IV B).

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

$$ \Lambda_p = s_p \Delta \tilde{r}_p + \frac{1}{2} \tilde{M}_{pqr} \Delta \tilde{r}_q \Delta \tilde{r}_r + \tilde{w}_p \varepsilon + \frac{1}{2} \tilde{B}_p^{(2)} \gamma^2. \quad (B1) $$

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

We focus on the limit of small $\varepsilon$ and $\gamma$. To simplify the subsequent calculations, we choose to approach this limit by Taylor expanding using the parameter $\eta \ll 1$ and setting $\varepsilon = \varepsilon_0\eta$ and $\gamma = \gamma_0\eta^{1/2}$. This way, for any given ratio $\varepsilon/\gamma^2 = \varepsilon_0/\gamma_0^2$, we can reach the small $\varepsilon$ and $\gamma$ limit by decreasing $\eta$ accordingly.

We now discuss the scaling of the $\Delta\tilde{r}_q$ and the $\Lambda_p$ with $\eta$. We assume that $\Delta\tilde{r}_q$ 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_q^{(k/2)}$ and $\lambda_p^{(k/2)}$, we respectively have:

$$ \Delta \tilde{r}_q = \sum_{k=0}^{\infty} \rho_q^{(k/2)} \eta^{k/2}, \quad (B2) $$

$$ \Lambda_p = \sum_{k=0}^{\infty} \lambda_p^{(k/2)} \eta^{k/2}. \quad (B3) $$

First, we find that $\rho_q^{(0)} = 0$ for all $q = 1, \dots, N_{\text{dof}}$, because at the transition, where $\eta \sim \varepsilon + b_\varepsilon \gamma^2 \to 0$, all $\Delta\tilde{r}_q = 0$. Then, insertion of Eqs. (B2) and (B3) into Eq. (B1) yields:

$$ \lambda_p^{(0)} = 0 \quad (B4) $$

$$ \lambda_p^{(1/2)} = s_p \rho_p^{(1/2)} \quad (B5) $$

$$ \lambda_p^{(1)} = s_p \rho_p^{(1)} + \frac{1}{2} \tilde{M}_{pqr} \rho_q^{(1/2)} \rho_r^{(1/2)} + \text{const.} \quad (B6) $$

To further discuss the values of $\rho_q^{(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}_q} \quad \text{for } q = 1, \dots, N_{\text{dof}}. \quad (B7) $$

From this we obtain using Eqs. (26) and (B1):

$$ 0 = s_q \Lambda_q + \Lambda_p \tilde{M}_{pqr} \Delta \tilde{r}_r \quad \text{for } q = 1, \dots, N_s - 1. \quad (B8) $$

$$ 0 = \Lambda_p \tilde{M}_{pqr} \Delta \tilde{r}_r \quad \text{for } q = N_s, \dots, N_{\text{dof}}. \quad (B9) $$

On the right-hand side of both equations, sums over $p = 1, \dots, N_s$ and $r = 1, \dots, N_{\text{dof}}$ are implied.

To order $\eta^{1/2}$, Eq. (B8) implies that $\lambda_p^{(1/2)} = 0$ for $p = 1, \dots, N_s - 1$. With Eq. (B5), this implies $\rho_p^{(1/2)} = 0$ for $p = 1, \dots, N_s - 1$. Furthermore, since $s_{N_s} = 0$, Eq. (B5) also implies that $\lambda_{N_s}^{(1/2)} = 0$.

To order $\eta^1$, Eq. (B8) also implies that $\lambda_p^{(1)} = 0$ for $p = 1, \dots, N_s - 1$. With Eq. (B6) follows that

$$ 0 = s_p \rho_p^{(1)} + \frac{1}{2} \tilde{M}_{pqr} \rho_q^{(1/2)} \rho_r^{(1/2)} + \text{const.} \quad \text{for } p < N_s. \quad (B10) $$

Since $\rho_q^{(1/2)} = 0$ for $q < N_s$, the implied sums in the second term effectively run only over $q, r = N_s, \dots, N_{\text{dof}}$. However, $\lambda_{N_s}^{(1)}$ is generally non-zero with

$$ \lambda_{N_s}^{(1)} = \frac{1}{2} \tilde{M}_{N_s q r} \rho_q^{(1/2)} \rho_r^{(1/2)} + \text{const.} \quad (B11) $$

Again, since $\rho_q^{(1/2)} = 0$ for $q < N_s$, the implied sums in the second term effectively run only over $q, r = N_s, \dots, N_{\text{dof}}$.

Consequently, up to quadratic order in $\eta$, the system energy is dominated by $\Lambda_{N_s}$, and Eq. (26) becomes $E = E_0 \Lambda_{N_s}^2 / 2$.

Finally, to order $\eta^{3/2}$, Eq. (B9) together with the fact that in the rigid regime $\lambda_{N_s}^{(1)} \neq 0$ implies that $\rho_q^{(1/2)} = 0$ for $q \ge N_s$, where we ignore here zero eigen vectors of the matrix $\tilde{M}_{N_s q r}$. Since we showed above that $\rho_q^{(1/2)} = 0$ also for $q < N_s$, we conclude that $\rho_q^{(1/2)} = 0$ for all $q = 1, \dots, N_{\text{dof}}$.

### 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. (42) that $\Lambda_p = 0$ for all $p = 1, \dots, N_s$. With Eq. (B5) 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_q^{(1/2)}$ with $q = N_s, \dots, N_{\text{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}_q \sim \eta^{1/2}$ for $q \ge N_s$. As a consequence, in all $\Lambda_p$ with $p = 1, \dots, N_s$ all terms $\sim \Delta\tilde{r}_q \Delta\tilde{r}_r$ with $q < N_s$ or $r < N_s$ can be neglected.
### 3. General case

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

$$Z \sim \int \left( \prod_{q=1}^{N_{\text{dof}}} d\Delta\tilde{r}_q \right) \exp \left[ -\frac{\beta E_0}{2} \sum_{p=1}^{N_s} \Lambda_p^2 \right]. \quad (B12)$$

Using Eqs. (B3) and (B5), this transforms into:

$$\begin{aligned} Z \sim & \int \left( \prod_{q=1}^{N_{\text{dof}}} d\Delta\tilde{r}_q \right) \\ & \times \exp \left[ -\frac{\beta E_0 \eta}{2} \sum_{p=1}^{N_s} s_p^2 [\rho_p^{(1/2)}]^2 \right] \\ & \times \exp \left[ -\beta E_0 \eta^{3/2} \sum_{p=1}^{N_s} s_p \rho_p^{(1/2)} \lambda_p^{(1)} \right] \\ & \times \exp \left[ -\frac{\beta E_0 \eta^2}{2} \sum_{p=1}^{N_s} \left( [\lambda_p^{(1)}]^2 + 2s_p \rho_p^{(1/2)} \lambda_p^{(3/2)} \right) \right] \\ & \times \exp \left[ -\frac{\beta E_0}{2} \text{(higher-order terms in } \eta \text{)} \right]. \end{aligned} \quad (B13)$$

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}_q = \rho_q^{(1/2)}\eta^{1/2}$ for $q = N_s, \dots, N_{\text{dof}}$. We thus obtain, using Eq. (B6):

$$\begin{aligned} Z \sim & \eta^{(N_{\text{dof}}+N_s-1)/2} \\ & \times \int \left( \prod_{p=1}^{N_s-1} d\rho_p^{(1)} \right) \left( \prod_{q=N_s}^{N_{\text{dof}}} d\rho_q^{(1/2)} \right) \\ & \times \exp \left[ -\frac{\beta E_0 \eta^2}{2} \sum_{p=1}^{N_s} \left( s_p \rho_p^{(1)} + \frac{1}{2} \sum_{q,r=N_s}^{N_{\text{dof}}} \tilde{M}_{pqr} \rho_q^{(1/2)} \rho_r^{(1/2)} + \text{const.} \right)^2 \right]. \end{aligned} \quad (B14)$$

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{aligned} Z \sim & \eta^{(N_{\text{dof}}+N_s-1)/2} \\ & \times \int \left( \prod_{q=N_s}^{N_{\text{dof}}} d\rho_q^{(1/2)} \right) \exp \left[ -\frac{\beta E_0 \eta^2}{2} \right. \\ & \qquad \left. \times \left( \frac{1}{2} \sum_{q,r=N_s}^{N_{\text{dof}}} \tilde{M}_{N_s q r} \rho_q^{(1/2)} \rho_r^{(1/2)} + \text{const.} \right)^2 \right] \\ & \times \prod_{p=1}^{N_s-1} \left( \int d\rho_p^{(1)} \exp \left[ -\frac{\beta E_0 \eta^2}{2} \left( s_p \rho_p^{(1)} + \frac{1}{2} \sum_{q,r=N_s}^{N_{\text{dof}}} \tilde{M}_{pqr} \rho_q^{(1/2)} \rho_r^{(1/2)} + \text{const.} \right)^2 \right] \right). \end{aligned} \quad (B15)$$

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

### Appendix C: Non-zero first-order terms in $\gamma$: $\tilde{B}_p^{(1)}$ and $\tilde{C}_{pq}^{(1)}$

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 q r}$ for $q, r = N_s, \dots, N_{\text{dof}}$, i.e. this matrix is strictly positive definite.

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

$$\tilde{\varepsilon} := \varepsilon - \Delta\varepsilon \quad (C1)$$

$$\tilde{\gamma} := \gamma - \Delta\gamma. \quad (C2)$$

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

$$\Lambda_{N_s} = \frac{1}{2} \tilde{M}_{N_s q r} \Delta \tilde{r}_q \Delta \tilde{r}_r + \tilde{w}_{N_s} \tilde{\varepsilon} + \tilde{C}_{N_s q}^{(1)} \Delta \tilde{r}_q \tilde{\gamma} + \frac{1}{2} \tilde{B}_{N_s}^{(2)} \tilde{\gamma}^2. \quad (C3)$$

Note that the $\Lambda_p$ with $p < N_s$ do in general still have a linear-order term in $\tilde{\gamma}$:

$$\begin{aligned} \Lambda_p = & s_p \Delta \tilde{r}_p + \frac{1}{2} \tilde{M}_{pqr} \Delta \tilde{r}_q \Delta \tilde{r}_r \\ & + \tilde{w}_p \tilde{\varepsilon} + \tilde{B}_p^{(1)} \tilde{\gamma} + \tilde{C}_{pq}^{(1)} \Delta \tilde{r}_q \tilde{\gamma} + \frac{1}{2} \tilde{B}_p^{(2)} \tilde{\gamma}^2. \end{aligned} \quad (C4)$$

Also, the $\Delta\tilde{r}_q$ may also be redefined due to the shifts in $\varepsilon$ and $\gamma$ due to a different non-dimensionalisation according to section IIC and a change in the compatibility matrix. If $\Delta\varepsilon$ and $\Delta\gamma$ are sufficiently small, it can be shown that they are given by $\Delta\varepsilon = (\tilde{B}_N^{(1)})^2/4\tilde{w}_N^2 b_\varepsilon$ and $\Delta\gamma = -\tilde{B}_N^{(1)}/2\tilde{w}_N b_\varepsilon$, where $\tilde{B}_N^{(1)}$ is the linear-order term before the shift, and $b_\varepsilon$ is given by Eq. (C10) below. 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}_{N_s}^{(1)} \sim N_s^{-1/2}$. Indeed, we numerically found earlier for phanomized triangular networks that $\tilde{B}_{N_s}^{(1)} \sim \Delta\gamma \sim N_{\text{dof}}^{-1/2} \sim N_s^{-1/2}$ [4].

We first discuss the athermal limit using arguments similar to appendix B 1. Using $\bar{\eta}$ as small parameter and fixing $\bar{\epsilon} \sim \bar{\eta}$ and $\bar{\gamma} \sim \bar{\eta}^{1/2}$, we use Eq. (B8) 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_p \Delta \tilde{r}_p + \tilde{B}_p^{(1)} \bar{\gamma}. \quad (C5)$$

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

$$\Delta \tilde{r}_p = - \frac{\tilde{B}_p^{(1)}}{s_p} \bar{\gamma} \quad \text{for } p < N_s. \quad (C6)$$

Inserting this into Eq. (C3) together with the diagonalization of the matrix $\tilde{M}_{N_sqr}$ with $q, r \ge N_s$, Eq. (47), we obtain up to order $\bar{\eta}^1$:

$$\Lambda_{N_s} = \frac{1}{2} \sum_{k=1}^{N_{\text{dof}}-N_s+1} \mu_k (\Delta q_k - \Delta q_k^0)^2 + \tilde{w}_{N_s} [\bar{\epsilon} + b_{\epsilon} \bar{\gamma}^2], \quad (C7)$$

where we introduced

$$\Delta q_k^0 := - \frac{\bar{\gamma}}{\mu_k} \sum_{r \ge N_s} v_r^k \zeta_r \quad (C8)$$

with

$$\zeta_r := \tilde{C}_{N_s r}^{(1)} - \sum_{p < N_s} \tilde{M}_{N_s rp} \frac{\tilde{B}_p^{(1)}}{s_p}, \quad (C9)$$

and

$$b_{\epsilon} := \frac{1}{2\tilde{w}_{N_s}} \left[ \tilde{B}_{N_s}^{(2)} + \sum_{p,q < N_s} \frac{\tilde{B}_p^{(1)}}{s_p} \tilde{M}_{N_s pq} \frac{\tilde{B}_q^{(1)}}{s_q} - 2 \sum_{p < N_s} \tilde{C}_{N_s p}^{(1)} \frac{\tilde{B}_p^{(1)}}{s_p} - \sum_{q,r \ge N_s} \zeta_q \frac{v_q^k v_r^k}{\mu_k} \zeta_r \right]. \quad (C10)$$

Remember that we assumed that the matrix $\tilde{M}_{N_sqr}$ for $q, r = N_s, \dots, N_{\text{dof}}$ is strictly positive definite, i.e. $\mu_k > 0$ for all $k$. In practice, to numerically compute $b_{\epsilon}$ we found that any modes $k$ with $\mu_k = 0$ can just be excluded from the last sum, where we numerically found that for these $k$, also $\sum_{q \ge N_s} \zeta_q v_q^k = 0$.

For the stiff-spring limit and the general case of finite spring stiffness and temperature, these arguments are essentially the same as for the athermal limit. In the stiff-spring limit, for non-vanishing $\tilde{B}_p^{(1)}$ and $\tilde{C}_{pq}^{(1)}$, the hyper-ellipsoid that describes the possible system configurations is now centered at $(\Delta q_k^0)$ (compare Eqs. (49) and (C7)). According to Eq. (C8) 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}_p^{(1)}/s_p)\bar{\gamma}/\bar{\eta}$.

Using the scaling arguments from appendix B, one can see that any further terms in the Taylor expansion, Eq. (18), will not change our results. Specifically, when allowing for finite $\tilde{B}_p^{(1)}$, we have that $\Delta \tilde{r}_q \sim \eta^{1/2}$ for all $q$, while $\Lambda_{N_s} \sim \eta^1$, and the $\Lambda_p$ with $p < N_s$ are of higher order than $\sim \eta^1$. Consequently, any higher-order term in Eq. (18) will only yield terms with order of at least $\sim \eta^{3/2}$.

## Appendix D: Relation between the Hessian and the matrix $\tilde{M}_{N_sqr}$

We evaluate the Hessian of the system,

$$H_{qr} := \frac{\partial^2 E}{\partial \Delta \tilde{r}_q \partial \Delta \tilde{r}_r}, \quad (D1)$$

in the rigid athermal regime, at an energy minimum close to the transition point, where we assume for simplicity that $\tilde{B}_p^{(1)} = \tilde{C}_{pq}^{(1)} = 0$ for all $p, q$.

Using Eq. (26):

$$H_{qr} = E_0 \sum_{p=1}^{N_s} \left( \frac{\partial \Lambda_p}{\partial \Delta \tilde{r}_q} \frac{\partial \Lambda_p}{\partial \Delta \tilde{r}_r} + \Lambda_p \frac{\partial^2 \Lambda_p}{\partial \Delta \tilde{r}_q \partial \Delta \tilde{r}_r} \right). \quad (D2)$$

Further, according to our discussion in appendix B 1, at an energy minimum, all $\Delta \tilde{r}_q$ 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$:

$$H_{qr} = E_0 \sum_{p=1}^{N_s-1} s_p^2 \delta_{pq} \delta_{pr} + E_0 \sum_{p=1}^{N_s-1} \sum_{s=1}^{N_{\text{dof}}} s_p \Delta \tilde{r}_s (\delta_{pq} \tilde{M}_{psr} + \tilde{M}_{psq} \delta_{pr}) + E_0 \Lambda_{N_s} \tilde{M}_{N_sqr}. \quad (D3)$$

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

For $q, r \ge N_s$ only the third term in Eq. (D3) is nonzero. Hence, to lowest order in $\eta$, the $N_{1st}^*$ nonzero eigen modes of $\tilde{M}_{N_sqr}$ for $q, r \ge N_s$, are those eigen modes of $H_{qr}$ 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$.

## Appendix E: Partition sum for the stiff-spring limit

Care needs to be taken when computing the partition sum for systems with fixed constraints [42]. In our
generalized spring networks, in the stiff-spring limit, all springs are imposed to be at their respective rest lengths. Yet, it is not immediately obvious how to impose these constraints to obtain the correct phase space volume. Naively, one may be tempted to explicitly impose the length of each spring, as we do in Eq. (39), to obtain a phase space volume:

$$ \Omega_{\text{def1}} = \int \left( \prod_{n=1}^{N_{\text{dof}}} dR_n \right) \prod_{i=1}^{N_s} \delta(L_i - L_{0i}). \quad (\text{E1}) $$

Yet, another possibility would be to just impose that the energy should vanish:

$$ \Omega_{\text{def2}} = \int \left( \prod_{n=1}^{N_{\text{dof}}} dR_n \right) \delta(E). \quad (\text{E2}) $$

After all, if the energy vanishes, all springs have to be at their respective rest length (compare Eq. (3)). However, both definitions, Eq. (E1) and Eq. (E2), 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, $K_i(L_i - L_{0i})$, the resulting phase space volumes $\Omega$ will differ between both definitions by a non-constant factor.

To obtain the correct expression for the phase space volume, we evaluate the partition sum

$$ Z = \int \left( \prod_{n=1}^{N_{\text{dof}}} dR_n \right) \exp \left( -\frac{\beta}{2} \sum_{i=1}^{N_s} K_i (L_i - L_{0i})^2 \right) \quad (\text{E3}) $$

in the stiff-spring limit. A Gaussian converges to a Dirac Delta in the limit of small standard deviation:

$$ \lim_{K_i \to \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}). \quad (\text{E4}) $$

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:

$$ Z = \left(\frac{2\pi}{\beta}\right)^{N_s/2} \prod_{i=1}^{N_s} K_i^{-1/2} \times \int \left(\prod_{n=1}^{N_{\text{dof}}} dR_n\right) \prod_{i=1}^{N_s} \delta(L_i - L_{0i}). \quad (\text{E5}) $$

This indicates that the correct approach is to use Dirac deltas on the spring lengths rather than the energy $E$. Note that in Eq. (39), 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 \to \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 Kth-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. (23), (49):

$$ \Delta\tilde{\ell}_{N_s} = \sum_{K=1}^{\infty} \frac{1}{(K+1)!} \sum_{k=1}^{N_{K\text{th}}^*} \mu_{K,k} \Delta q_{K,k}^{K+1} + \frac{1}{2} \tilde{B}_{N_s}^{(2)} \gamma^2. \quad (\text{F1}) $$

Here, $(K+1)!$ denotes the faculty of $K+1$, the variables $\Delta q_{K,k}$ denote the amplitudes of the Kth-order zero modes, and the associated coefficients are $\mu_{K,k} > 0$. For simplicity, we set again $\tilde{B}_p^{(1)} = \tilde{C}_{pq}^{(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\text{th}}^* = 0$ for even $K$. This is because otherwise, a Kth-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\text{th}}^* > 0$ for even $K$, because a Kth-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. (F1) actually only includes the odd integers.

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

$$ \Omega_{N_s} \sim I(\varepsilon, \gamma) := \int \prod_{K=1}^{\infty} \prod_{k=1}^{N_{K\text{th}}^*} d\Delta q_{K,k} \times \delta \left( \sum_{K=1}^{\infty} \frac{1}{(K+1)!} \sum_{k=1}^{N_{K\text{th}}^*} \mu_{K,k} \Delta q_{K,k}^{K+1} + \tilde{w}_{N_s} [\varepsilon + b_{\varepsilon} \gamma^2] \right), \quad (\text{F2}) $$

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

$$ I_1 := \int \prod_{K=1}^{\infty} \prod_{k=1}^{N_{K\text{th}}^*} dx_{K,k} \times \delta \left( \sum_{K=1}^{\infty} \frac{1}{(K+1)!} \sum_{k=1}^{N_{K\text{th}}^*} \mu_{K,k} x_{K,k}^{K+1} + \tilde{w}_{N_s} \right), \quad (\text{F3}) $$

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

$$ \Delta q_{K,k} = [\varepsilon + b_{\varepsilon} \gamma^2]^{1/(K+1)} x_{K,k}, \quad (\text{F4}) $$
which yields:

$$ \Omega_{N_s}(\varepsilon, \gamma) \sim I(\varepsilon, \gamma) = [\varepsilon + b_{\varepsilon}\gamma^2]^{N_{\Omega}} I_1 \quad (F5) $$

with:

$$ N_{\Omega} = -1 + \sum_{K=1}^{\infty} \frac{N_{K\text{th}}^*}{K+1}. \quad (F6) $$

The first term in $N_Ω$ arises from removing the factor $[\varepsilon + b_εγ^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\ell_{N_s}$ invariant. They do not contribute to $N_Ω$. This is because any substitution like in Eq. (F4), 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\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_{1st}^* > 0$ and $N_{Kth}^* = 0$ for $K > 1$, we have indeed $N_Ω = (N_{1st}^* - 2)/2$, consistent with Eq. (F6). Taken together, as indicated by Eq. (F6), also higher-order zero modes at the transition may contribute to the exponent $N_Ω$; 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_Ωk_B/DV^*$.

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

To obtain the isotropic tension $t := (\partial F/\partial \varepsilon)/DV$, we first use the chain rule. To lowest order, $V = V^*$, and so we obtain from Eq. (77):

$$ t = \frac{1}{DV^*} \left[ \frac{\partial \bar{F}}{\partial \hat{\varepsilon}_S} \frac{\partial \hat{\varepsilon}_S}{\partial \varepsilon} + \frac{\partial \bar{F}}{\partial \varepsilon} + \frac{1}{2\beta \bar{F}''} \left( \frac{\partial \bar{F}''}{\partial \hat{\varepsilon}_S} \frac{\partial \hat{\varepsilon}_S}{\partial \varepsilon} + \frac{\partial \bar{F}''}{\partial \varepsilon} \right) \right]. \quad (G1) $$

We compute to lowest order:

$$ \frac{\partial \hat{\varepsilon}_S}{\partial \varepsilon} = - \frac{\hat{\varepsilon}_S + b_{\varepsilon} \gamma^2}{|\varepsilon + b_{\varepsilon} \gamma^2| \sqrt{1 + \theta}} \quad (G2) $$

$$ \frac{1}{DV^*} \frac{\partial \bar{F}}{\partial \hat{\varepsilon}_S} = -\kappa_E(\varepsilon - \hat{\varepsilon}_S) - \frac{\kappa_S T}{\hat{\varepsilon}_S + b_{\varepsilon} \gamma^2} = 0 \quad (G3) $$

$$ \frac{1}{DV^*} \frac{\partial \bar{F}}{\partial \varepsilon} = \kappa_E(\varepsilon - \hat{\varepsilon}_S) \quad (G4) $$

$$ \frac{1}{DV^*} \frac{\partial \bar{F}''}{\partial \hat{\varepsilon}_S} = - \frac{2\kappa_S T}{(\hat{\varepsilon}_S + b_{\varepsilon} \gamma^2)^3} \quad (G5) $$

$$ \frac{1}{DV^*} \frac{\partial \bar{F}''}{\partial \varepsilon} = 0. \quad (G6) $$

We thus obtain:

$$ t = \kappa_E(\varepsilon - \hat{\varepsilon}_S) \left[ 1 + \frac{1}{2(N_{1st}^* - 2)} \frac{\theta}{1 + \theta} \right], \quad (G7) $$

where we have also used that

$$ \bar{F}'' = -DV^* \kappa_E \frac{|\varepsilon + b_\varepsilon \gamma^2| \sqrt{1 + \theta}}{\hat{\varepsilon}_S + b_\varepsilon \gamma^2} \quad (G8) $$

and

$$ \kappa_E(\varepsilon - \hat{\varepsilon}_S) = - \frac{\kappa_S T}{\hat{\varepsilon}_S + b_\varepsilon \gamma^2}, \quad (G9) $$

which follows from Eq. (75). The term $\sim (N_{1st}^* - 2)^{-1}$ in Eq. (G7) corresponds to the terms involving $\bar{F}''$.

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

$$ G = \frac{1}{V^*} \left[ \frac{\partial^2 \bar{F}(\hat{\varepsilon}_S(\varepsilon, \gamma); \varepsilon, \gamma)}{\partial \gamma^2} + \frac{1}{2\beta} \frac{\partial^2}{\partial \gamma^2} \log \bar{F}''(\hat{\varepsilon}_S(\varepsilon, \gamma); \varepsilon, \gamma) \right]. \quad (G10) $$

We compute with Eqs. (G3) and (G9):

$$ \begin{aligned} \frac{1}{V^*} \frac{\partial \bar{F}(\hat{\varepsilon}_S(\varepsilon, \gamma); \varepsilon, \gamma)}{\partial \gamma} &= \frac{1}{V^*} \frac{\partial \bar{F}(\hat{\varepsilon}_S; \varepsilon, \gamma)}{\partial \gamma} \\ &= -\frac{2D\kappa_S T b_\varepsilon \gamma}{\hat{\varepsilon}_S + b_\varepsilon \gamma^2} \\ &= 2D\kappa_E(\varepsilon - \hat{\varepsilon}_S)b_\varepsilon\gamma, \end{aligned} \quad (G11) $$

and further:

$$ \begin{aligned} \frac{1}{V^*} \frac{\partial^2 \bar{F}(\hat{\varepsilon}_S(\varepsilon, \gamma); \varepsilon, \gamma)}{\partial \gamma^2} &= \frac{1}{V^*} \left[ \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} \right] \end{aligned} \quad (G12) $$

with

$$ \frac{1}{V^*} \frac{\partial^2 \bar{F}(\hat{\varepsilon}_S; \varepsilon, \gamma)}{\partial \gamma \partial \hat{\varepsilon}_S} = -2D\kappa_E b_\varepsilon \gamma \quad (G13) $$

$$ \frac{\partial \hat{\varepsilon}_S}{\partial \gamma} = -b_\varepsilon \gamma \left( 1 + \frac{\varepsilon + b_\varepsilon \gamma^2}{|\varepsilon + b_\varepsilon \gamma^2|\sqrt{1 + \theta}} \right) \quad (G14) $$

$$ = -2b_\varepsilon \gamma \frac{\varepsilon - \hat{\varepsilon}_S}{|\varepsilon + b_\varepsilon \gamma^2|\sqrt{1 + \theta}} \quad (G15) $$

$$ \frac{1}{V^*} \frac{\partial^2 \bar{F}(\hat{\varepsilon}_S; \varepsilon, \gamma)}{\partial \gamma^2} = 2D\kappa_E(\varepsilon - \hat{\varepsilon}_S)b_\varepsilon. \quad (G16) $$

Thus:

$$ \begin{aligned} \frac{1}{V^*} \frac{\partial^2 \bar{F}(\hat{\varepsilon}_S(\varepsilon, \gamma); \varepsilon, \gamma)}{\partial \gamma^2} &= 2D\kappa_E(\varepsilon - \hat{\varepsilon}_S)b_\varepsilon \left( 1 + \frac{2b_\varepsilon \gamma^2}{|\varepsilon + b_\varepsilon \gamma^2|\sqrt{1 + \theta}} \right). \end{aligned} \quad (G17) $$
Furthermore:

$$
\begin{align*}
& \frac{\partial \bar{F}''(\hat{\varepsilon}_S(\varepsilon, \gamma); \varepsilon, \gamma)}{\partial \gamma} \\
&= \frac{\partial \bar{F}''(\hat{\varepsilon}_S; \varepsilon, \gamma)}{\partial \hat{\varepsilon}_S} \frac{\partial \hat{\varepsilon}_S}{\partial \gamma} + \frac{\partial \bar{F}''(\hat{\varepsilon}_S; \varepsilon, \gamma)}{\partial \gamma} \\
&= D\mathcal{V}^* \kappa_E \left[ \left( \frac{1}{(\hat{\varepsilon}_S + b_\varepsilon \gamma^2) - \frac{\hat{\varepsilon}_S - \varepsilon}{(\hat{\varepsilon}_S + b_\varepsilon \gamma^2)^2}} \right) \right. \\
&\qquad \left. \times \left( -2b_\varepsilon \gamma \frac{\varepsilon - \hat{\varepsilon}_S}{|\varepsilon + b_\varepsilon \gamma^2|\sqrt{1+\theta}} \right) \right. \\
&\qquad \left. - 2b_\varepsilon \gamma \frac{\hat{\varepsilon}_S - \varepsilon}{(\hat{\varepsilon}_S + b_\varepsilon \gamma^2)^2} \right] \\
&= -2D\mathcal{V}^* \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}} \\
&\qquad \times \left[ \varepsilon + b_\varepsilon \gamma^2 - |\varepsilon + b_\varepsilon \gamma^2| \sqrt{1+\theta} \right] \\
&= -4D\mathcal{V}^* \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}}, \quad (G18)
\end{align*}
$$

and thus:

$$
\frac{\partial}{\partial \gamma} \log \bar{F}''(\hat{\varepsilon}_S(\varepsilon, \gamma); \varepsilon, \gamma) = 4b_{\varepsilon}\gamma \frac{\varepsilon - \hat{\varepsilon}_S}{(\varepsilon + b_{\varepsilon}\gamma^2)^2(1 + \theta)}. \quad (G19)
$$

The second derivative is:

$$
\begin{align*}
& \frac{\partial^2}{\partial \gamma^2} \log \bar{F}''(\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)} \left[ 1 + \frac{2b_{\varepsilon}\gamma^2}{|\varepsilon + b_{\varepsilon}\gamma^2|\sqrt{1 + \theta}} \right. \\
&\qquad \left. - \frac{4b_{\varepsilon}\gamma^2}{(\varepsilon + b_{\varepsilon}\gamma^2)(1 + \theta)} \right]. \tag{G20}
\end{align*}
$$

From Eqs. (61) and (G9) follows:

$$
\frac{1}{2\beta V^*} = \frac{D\kappa_S T}{N_{1st}^* - 2} = \frac{D\kappa_E}{N_{1st}^* - 2} \theta \frac{(\varepsilon + b_\varepsilon \gamma^2)^2}{4} \quad (G21)
$$

Combining Eqs. (G10), (G17), (G20), and (G21), we ob-
tain Eq. (79) in the main text.

[1] J. C. Maxwell, Philosophical Magazine Series 4 **27**, 294 (1864).

[2] C. R. Calladine, International Journal of Solids and Structures **14**, 161 (1978).

[3] M. Merkel, K. Baumgarten, B. P. Tighe, and M. L. Manning, Proceedings of the National Academy of Sciences **116**, 6560 (2019).

[4] C.-T. Lee and M. Merkel, Soft Matter **18**, 5410 (2022).

[5] T. C. Lubensky, C. L. Kane, X. Mao, A. Souslov, and K. Sun, Reports on Progress in Physics **78**, 73901 (2015), arXiv:1503.01324.

[6] M. Mannattil, J. M. Schwarz, and C. D. Santangelo, Physical Review Letters **128**, 208005 (2022), arXiv:2112.04279.

[7] O. K. Damavandi, V. F. Hagh, C. D. Santangelo, and M. L. Manning, Physical Review E **105**, 025004 (2022).

[8] O. K. Damavandi, V. F. Hagh, C. D. Santangelo, and M. L. Manning, Physical Review E **105**, 025003 (2022).

[9] S. Alexander, Physics Report **296**, 65 (1998),
arXiv:1512.04885.

[10] B. Cui, G. Ruocco, and A. Zaccone, Granular Matter **21**, 69 (2019), arXiv:1901.09582.

[11] S. Zhang, E. Stanifer, V. V. Vasisht, L. Zhang, E. Del Gado, and X. Mao, Physical Review Research **4**, 043181 (2022).

[12] M. Wyart, Annales de Physique **30**, 1 (2005).

[13] M. Merkel and M. L. Manning, New Journal of Physics **20**, 022002 (2018).

[14] P. R. Onck, T. Koeman, T. van Dillen, and E. van der Giessen, Physical Review Letters **95**, 178102 (2005),
arXiv:cond-mat/0502397.

[15] C. P. Broedersz and F. C. MacKintosh, Reviews of Modern Physics **86**, 995 (2014), arXiv:1404.4332.

[16] 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.

[17] A. Sharma, A. J. Licup, K. A. Jansen, R. Rens, M. Shein-
man, G. H. Koenderink, and F. C. MacKintosh, Nature
Physics **12**, 584 (2016), arXiv:1506.07792.

[18] M. Moshe, M. J. Bowick, and M. C. Marchetti, Physical Review Letters **120**, 268105 (2017), arXiv:1708.07848.

[19] D. M. Sussman and M. Merkel, Soft Matter **14**, 3397 (2018).

[20] 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.

[21] M. Plischke and B. Joós, Physical Review Letters **80**, 4907 (1998).
[22] M. Dennison, M. Sheinman, C. Storm, and F. C. MacKintosh, Physical Review Letters 111, 095503 (2013).

[23] M. C. Wigbers, F. C. Mackintosh, and M. Dennison, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 92, 1 (2015), arXiv:1410.7860.

[24] L. Zhang and X. Mao, Physical Review E 93, 1 (2016).

[25] S. Chen, T. Markovich, and F. C. MacKintosh, Physical Review Letters 130, 088101 (2023).

[26] Reference to be inserted later.

[27] 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].

[28] More precisely, we define the number of mechanisms $N_{\text{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} \to \mathbb{R}^{N_{\text{dof}}}$ and some $\eta > 0$ such that Eq. (13) holds within some open ball of radius $\eta$, i.e. for $|\vec{u}| < \eta$.

[29] 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 \to \infty$. [42]. The latter is of course the physically more relevant approach, since there are no exact geometric constraints in nature.

[30] 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}^*$, $\vec{y}'(0) = \vec{z}_k$, and $d^2\ell_i(\vec{y}(u), \gamma)/du^2 = 0$. Indeed, such a function is given by:

$$y_m(u) := r_m^* + uV_{mq}v_q^k - \frac{u^2}{2} \sum_{p=1}^{N_s-1} \frac{V_{mp}}{s_p} \tilde{M}_{pqr}v_q^k v_r^k$$

with $m = 1, \dots, N_{\text{dof}}$.

[31] For any two-times differentiable function $u \mapsto \vec{y}(u)$ with $\vec{y}(0) = \vec{r}^*$ and $\vec{y}'(0) = \vec{z}_k$, it is impossible to have $d^2\ell_i(\vec{y}(u), \gamma)/du^2 = 0$ for all $i = 1, \dots, N_s$, because $d^2\Delta\ell_{N_s}(\vec{y}(u), \gamma)/du^2 = \mu_k > 0$.

[32] M. E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation, illustrated edition ed. (Oxford University Press, Oxford, 2010).

[33] S. Arzash, A. Gannavarapu, and F. C. MacKintosh, Phys. Rev. E 108, 054403 (2023).

[34] D. Z. Rocklin, V. Vitelli, and X. Mao, Folding mechanisms at finite temperature (2018), arXiv:1802.02704 [cond-mat].

[35] A. Ninarello, J. Ruiz-Franco, and E. Zaccarelli, Nature Communications 13, 527 (2022).

[36] N. Javerzat and M. Bouzid, Physical Review Letters 130, 268201 (2023).

[37] N. Javerzat, Physical Review Letters 132, 018201 (2024).

[38] R. Rens, C. Villarroel, G. Düring, and E. Lerner, Physical Review E 98, 062411 (2018).

[39] E. Lerner and E. Bouchbinder, Extreme Mechanics Letters 65, 102104 (2023).

[40] M. Merkel, From Cells to Tissues: Tissue Remodeling and Polarity Reorientation in Shearing Epithelia, Ph.D. thesis, Technische Universität Dresden (2014).

[41] A. W. Lees and S. F. Edwards, Journal of Physics C: Solid State Physics 5, 1921 (1972).

[42] N. G. van Kampen and J. J. Lodder, American Journal of Physics 52, 419 (1984).