Cheng-Tai Lee

Matthias Merkel

Generic elasticity of thermal, under-constrained systems

Athermal (i.e. zero-temperature) under-constrained systems are typically floppy, but they can be rigidified by the application of external strain, which is theoretically well understood. Here and in the companion paper, we extend this theory to finite temperatures for a very broad class of underconstrained systems. In the vicinity of the athermal transition point, we derive from first principles expressions for elastic properties such as isotropic tension t and shear modulus G on temperature T , isotropic strain ε, and shear strain γ, which we confirm numerically. These expressions contain only three parameters, which depend on the microscopic structure of the system. These respectively describe entropic rigidity, energetic rigidity, and an interaction between isotropic and shear strain. Our results imply that in under-constrained systems, entropic and energetic rigidity interact like two springs in series. This also allows for a simple explanation of the previously observed scaling relation t ∼ G ∼ T 1/2 at ε = γ = 0. Our work unifies the physics of systems as diverse as polymer fibers & networks, membranes, and vertex models for biological tissues.

Understanding rigidity of amorphous materials still represents a fundamental challenge, which is relevant for questions ranging from the glass transition in inert materials [1] to solid-fluid transitions in living systems [2].

In the athermal limit, i.e. at zero temperature, the rigidity of many systems, including spring networks and particulate matter, can be predicted using Maxwell-Calladine constraint counting [3][4][5][6]. In this approach, one compares the number of constraints, N s (e.g. the number of springs in a spring network), to the number of degrees of freedom, N dof (e.g. the node positions in a spring network). Roughly, if N s > N dof , the system is expected to be rigid. Conversely, if N s < N dof , the system is under-constrained and expected to be floppy.

However, it is known that even under-constrained systems can be rigidified [4,[7][8][START_REF] Merkel | Manning[END_REF][10][11][12]. This occurs for instance when there is geometric incompatibility between spring FIG. 1.

Analytical predictions for the isotropic tension t depending on isotropic strain ε, shown for the athermal limit (red dashed line), the infinitely-stiff-spring limit (blue dashed line), and the general case of finite spring stiffness and temperature (black solid line). (inset) Same plot in log-log scaling.

rest lengths and externally applied strain. For example, in the simple spring network in Fig. 2b, we illustrate what happens when progressively increasing the externally applied isotropic strain ε while keeping the spring rest lengths fixed. Eventually, as ε is increased above a critical value, here for convenience defined as the point of zero strain, ε = 0, the springs necessarily need to be stretched. At the boundary of geometric incompatibility (ε = 0), a so-called state of self-stress (SSS) emerges [6]. This means that specific combinations of virtual tensions can be put on the springs without resulting in any net forces on any of the (movable) nodes. In previous work, we used the properties of the SSS created at the (a) Schematic: Variation of the accessible phase space volume Ω with isotropic strain ε in the limit of infinitely stiff springs. Starting at the athermal transition point at ε = 0, upon isotropic compression (ε < 0), Ω is expected to increase with the distance |ε| to the athermal transition point. The abscissa represents the N dof -dimensional configuration space. (b) Illustration of a simple two-spring system in 3D of linear dimension 3 √ V = e ε 3 √ V * , where the black nodes have fixed positions and the blue node is freely movable. In panels a and b, blue dashed lines represent the accessible parts of the phase space.

rigidity transition to analytically derive the elastic properties in the athermal limit [START_REF] Merkel | Manning[END_REF]13]. Our approach applies not only to spring networks, but also to other systems such as polymer networks [14][15][16][17] and vertex models for biological tissues [START_REF] Merkel | Manning[END_REF]10,[18][19][20].

While we analytically understand rigidity in the athermal limit, there are mostly only numerical results for the case of a finite temperature, i.e. for thermal underconstrained systems [21][22][23][24][25]. For instance, previous studies showed numerically that upon increasing temperature T from the athermal rigidity transition point, the shear modulus G scales as G ∼ T 1/2 [22,23]. For specific spring networks, this scaling could be derived analytically using effective-medium theory (EMT) before [26], but a derivation from first principles independent from the assumptions inherent to EMT has been lacking so far. There is also work treating a related abstract constraint satisfaction problem using the replica approach [27], work studying these questions in isostatic or overconstrained systems [26,[28][START_REF] Rocklin | Folding mechanisms at finite temperature[END_REF][START_REF] Chen | [END_REF], and work discussing the effect of singularities on under-constrained systems [31]. However, despite these advances, we still miss a general understanding for the elastic properties of thermal, under-constrained systems, in particular their response to isotropic and shear strain.

Here and in the companion paper [32], we develop a generic analytical theory for thermal, under-constrained systems in the vicinity of the athermal transition point. We first derive the elastic system properties in the limit of infinitely stiff springs, where only entropic elasticity plays a role. To then discuss the general case of finite spring stiffness and finite temperature (black solid line in Fig. 1), we then combine the stiff-spring results (blue dashed line) with our earlier results for the athermal limit, where only energetic elasticity plays a role (red dashed line) [START_REF] Merkel | Manning[END_REF]13]. We provide an intuition for why the system behaves as if the two limiting cases were "put in series", i.e. entropic and energetic strains add up (illustration in Fig. 1). Finally, we test our analytical results numerically, using simulations of randomly cut triangular networks. With only three fit parameters, our analytical results reproduce the observed behavior of tension t and shear modulus G over many orders of magnitude of varying isotropic strain ε, shear strain γ, and temperature T .

To discuss the key ideas of our approach, we focus here on an under-constrained network of linear springs with equal spring constants K and rest lengths L 0 . The more general case of varying spring constants and rest lengths is discussed in the companion paper [32]. The system energy is:

E = K 2 Ns i=1 (L i -L 0 ) 2 , (1)
where the sum is over all springs i. The length L i of spring i depends on isotropic strain ε, shear strain γ, and the positions R a of the nodes it connects, where a is the node index. We use periodic boundary conditions in D spatial dimensions with periodic box volume V . Isotropic strain ε is defined as linear strain, i.e. V = e Dε V * , where V * is the system volume at the athermal transition point. Shear strain γ can be either pure shear strain, or simple shear strain using Lees-Edwards-like boundary conditions [33].

It will be convenient to consider geometric quantities that are rescaled with respect to system size. For instance, we consider the rescaled node positions r a :=

V -1/D R a = e -ε (V * ) -1/D R a .
Because the springs in our network are linear, the spring lengths are homogeneous with respect to isotropic strain. Thus, the rescaled spring lengths i ({r a }, γ) := e -ε L i ({R a }, ε, γ)/L 0 do not depend on isotropic strain ε any more.

In the athermal limit, the elastic properties arise from purely energetic interactions, and their analytical expressions have been derived before [START_REF] Merkel | Manning[END_REF]13]. Briefly, for small isotropic and shear strain, the system is rigid for ε + b ε γ 2 > 0, where b ε is a constant that depends on the microscopic network structure. In the rigid regime close to the transition, energetically created isotropic tension t E and shear modulus G E are given by (red dashed line in Fig. 1) [32]:

t E = κ E ε + b ε γ 2 (2)
G E = 2Db ε κ E ε + 3b ε γ 2 , (3)
where κ E is a constant that depends on the microscopic network structure.

In the thermal case, the configurational degeneracy of the system becomes important (Fig. 2). To discuss its influence, we first consider the stiff-spring limit, K → ∞, which implies a hard constraint for all spring lengths: L i = L 0 . To sketch the derivation of the elastic system properties in this limit we set for simplicity γ = 0 here. The general case including γ is described in the companion paper [32]. To derive how the accessible phase space volume Ω scales with ε, we first note that the spring length constraints L i = L 0 link the rescaled spring lengths i to isotropic strain ε: To linear order in ε, the definition of i implies: ∆ i := i -1 = -ε. We then Taylor-expand all rescaled spring lengths around the transition:

-ε ≡ ∆ i = C in ∆r n + 1 2 M imn ∆r m ∆r n , (4)
where we use Einstein notation (implied sum over equal indices). We introduced ∆r n := r n -r * n , where the r n with n = 1, . . . , N dof correspond to all components of the rescaled node positions r a , and r * n denotes their respective values at the transition point. The matrix C in is the compatibility matrix at the transition point, and M imn is the second spring length derivative at the transition point. Performing a singular-value decomposition on C in , we can express Eq. ( 4) in terms of the "eigen modes" of C in :

-ε wp = s ṗ∆r ṗ + 1 2 Mpmn ∆r m ∆r n , (5)
where the s p with p = 1, . . . , N s are the singular values of C in , indices with underdot are not summed over, and the quantities with tilde were transformed into the respective "eigen space" of C in . In particular, the N s -dimensional vector ( wp ) corresponds to the N s -dimensional vector (w i ) = (1, . . . , 1) transformed into the eigen space of C in .

In general, in disordered networks, a single SSS is created at the transition point [6,[START_REF] Merkel | Manning[END_REF]. Sorting the s p in decreasing order, we thus have that s p > 0 for p < N s and s Ns = 0.

Retaining only the lowest-order terms in ε, we obtain from Eq. ( 5) after some discussion [32]:

-ε wp = s ṗ∆r ṗ+ 1 2 N dof m,n=Ns
Mpmn ∆r m ∆r n for p = 1, . . . , N s . ( 6)

Each of these N s equations includes all the ∆r n with n ≥ N s . However, each of the ∆r p with p < N s appears only once in total; as linear terms in the equation with index p. Thus, for a given value of ε, the first N s -1 equations only fix the values of the ∆r p with p < N s . Meanwhile, the last equation, p = N s , puts a constraint on the ∆r n with n ≥ N s . Because s Ns = 0, this equation describes the surface of a hyper-ellipsoid (compare Fig. 2b) whose half axes scale as ∼ √ -ε. Some of the eigen values of the matrix MNsmn are zero; for instance those that correspond to a global translation of all nodes. The ellipsoid is thus N e -dimensional, where N e is the number of non-zero eigen values of MNsmn , which fulfills N e ≤ N dof -N s -D + 1. Thus, the surface of the ellipsoid scales as ∼ (-ε) (Ne-1)/2 . Together with an additional prefactor of ∼ 1/ √ -ε due to the constraint on the spring lengths [32], we obtain for the accessible phase space volume:

Ω ∼ (-ε) Ne -2 2 . ( 7
)
The free energy consists only of an entropic contribution

F S = -k B T log Ω = -k B T ([N e -2]/2) log [-(ε + b ε γ 2 )],
where k B denotes the Boltzmann constant, and we neglected constant offsets. Also, we have included the dependency on γ, which we derive in the companion paper [32]. We thus obtain for the purely entropic tension t S = (∂F S /∂ε)/V and shear modulus G S = (∂ 2 F S /∂γ 2 )/V :

t S = - κ S T ε + b ε γ 2 (8)
G S = -2Db ε κ S T ε -b ε γ 2 (ε + b ε γ 2 ) 2 , (9)
with κ S = k B (N e -2)/2DV * . Hence, not only energetic rigidity [START_REF] Merkel | Manning[END_REF]13], but also entropic rigidity is determined by 11) and ( 12) using a single parameter set (κE, κS, bε) [START_REF]We extracted κE ≈ 0.161 from the tension data for the athermal limit using a linear fit. The value of κS/kB ≈ 0.0629 was then obtained from the fit in the inset to Fig. 3a. The parameter bε ≈ 0.835 describing the interaction between isotropic and shear strain was ob[END_REF].

the properties of the SSS that is created at the athermal transition point.

For γ = 0, the scaling of the tension is plotted in Fig. 1 (blue dashed line). This scaling is consistent with known results for some under-constrained systems. For instance, 11) and ( 12) using the same parameter set as in Fig. 3.

the tension of a freely jointed chain scales as ∼ -T /ε close to the fully stretched chain configuration [START_REF] Rubinstein | Polymer Physics[END_REF], and so does the tension of an incompressible membrane without bending rigidity close to the fully stretched state [36].

We now discuss the general case of finite spring stiffness K and finite temperature T , where entropic and energetic elasticity compete with each other. Entropic effects generally tend to elongate individual springs in order to increase the phase space volume, but strongly elongated springs are less favorable energetically. To get a more quantitative intuition for the interplay between entropic and energetic elasticity, we consider boundary conditions that ensure γ = 0 and constant tension t > 0 (horizontal gray line in Fig. 1). In the thermal stiff-spring limit, K → ∞, this leads to a negative strain ε S < 0 (blue arrow), which is according to Eq. ( 8): ε S = -κ S T /t (blue dashed line). When we now allow for a finite spring constant K, the springs additionally stretch, increasing the total strain by some amount ε E to the total strain ε = ε S + ε E (red and black arrows, respectively). For sufficiently small temperatures, the network geometry is on average not very different from the athermal network. So, using Eq. ( 2), we have ε E = t/κ E (red dashed line). Thus, the total strain depends on tension as (black solid line)

ε = - κ S T t + t κ E . (10)
In other words, the interaction between entropic and energetic elasticity corresponds to an entropic and an energetic spring acting in series (Fig. 1). The rigorous analysis in our companion paper confirms that this intuitive explanation leads to the exact result for N e 1 [32]. Inverting Eq. ( 10) we find:

t = κ E 2 ε + b ε γ 2 + ε + b ε γ 2 √ 1 + θ , (11)
where

θ := 4κ S T /[κ E (ε + b ε γ 2 ) 2 ]
. The γ dependency that we have included here is derived in the companion paper [32]. We find for the shear modulus [32]:

G = 2Db ε t 1 + 2b ε γ 2 |ε + b ε γ 2 | √ 1 + θ . (12)
Thus, the shear modulus is proportional to tension for γ = 0, but not for varying γ.

Our approach also provides a simple explanation for the previously observed G ∼ T 1/2 scaling for ε = 0 and γ = 0. In this limit, we find from Eqs. (11) and ( 12) that t = √ κ E κ S T and G = 2Db ε t. The relation for t can also be obtained immediately from Eq. ( 10) with ε = 0.

More generally, we find three scaling regimes, which we discuss here for γ = 0 (Fig. 3a and inset to Fig. 1). (i) For θ 1 and ε < 0, the system behaves like in the stiff-spring limit with G ∼ t ∼ T 1 |ε| -1 . (ii) For θ 1 and ε > 0, the system behaves like in the athermal regime with G ∼ t ∼ T 0 ε 1 . (iii) For θ 1, the system shows the scaling G ∼ t ∼ T 1/2 ε 0 . This is remarkably analogous to previous results [26], where strain is replaced by the network connectivity. We discuss a possible reason for this in the companion paper [32].

We tested our analytical results using Monte-Carlo (MC) simulations of a 2D under-constrained randomly cut triangular network in periodic boundary conditions [37]. We first set shear strain to zero, γ = 0, varying only isotropic strain ε and temperature T (Fig. 3). We confirmed the t ∼ T 1/2 scaling for ε = 0 (inset to Fig. 3a). Furthermore, Eqs. (11) and (12) predict that for γ = 0 the tension-strain and shear-modulus-strain curves respectively collapse when rescaling each of t, G, and ε by √ T , which we confirm numerically in Fig. 3b,c. The theoretical predictions are shown as dashed and solid curves, where the values of κ E , κ S , and b ε were obtained independently, from data at T = 0 and at ε = 0 [START_REF]We extracted κE ≈ 0.161 from the tension data for the athermal limit using a linear fit. The value of κS/kB ≈ 0.0629 was then obtained from the fit in the inset to Fig. 3a. The parameter bε ≈ 0.835 describing the interaction between isotropic and shear strain was ob[END_REF]. Note that the value of κ S can also directly be obtained from N e , by counting the non-zero eigen values of MNsmn with m, n ≥ N s . Indeed, for our network studied here, we find N e = 433, which is consistent with the observed value for κ S [START_REF]We extracted κE ≈ 0.161 from the tension data for the athermal limit using a linear fit. The value of κS/kB ≈ 0.0629 was then obtained from the fit in the inset to Fig. 3a. The parameter bε ≈ 0.835 describing the interaction between isotropic and shear strain was ob[END_REF].

We further tested our analytical results for varying shear strain, γ > 0, fixing isotropic strain ε (Fig. 4). For the tension t, Eq. ( 11) again implies a scaling collapse onto a theoretical curve, which we confirm numerically in Fig. 4a. The solid and dashed lines in Fig. 4a and inset indicate again the prediction according to Eq. ( 11), using the same parameter values as before [START_REF]We extracted κE ≈ 0.161 from the tension data for the athermal limit using a linear fit. The value of κS/kB ≈ 0.0629 was then obtained from the fit in the inset to Fig. 3a. The parameter bε ≈ 0.835 describing the interaction between isotropic and shear strain was ob[END_REF]. Similarly, the shear modulus data follows the analytical prediction in Eq. ( 12) (Fig. 4b). Deviations in the shear modulus at the lowest temperature are consistent with slightly insufficient equilibration, and deviations at larger shear strain are likely due to higher-order terms.

In conclusion, we have developed a generic analyical theory for the elastic properties of thermal, underconstrained systems, and their behavior under isotropic and shear strain. We expect this work to unify the physics of a broad class of materials, including polymer fibers, polymer networks, membranes, and vertex models for biological tissues.

[36] M. Durand, Soft Matter 18, 3891 (2022).

[37] The network was initialized with 40 × 40 nodes with connectivity z = 3.2 [13,[START_REF] Broedersz | [END_REF][39][40], were we used homogeneous spring constant K = 1 and rest length L0 = 1. Dangling springs and isolated islands were removed and were not counted towards the connectivity z [13]. Before running the actual MC simulations, we removed any asymmetry with respect to shear strain γ by shear-stabilizing the network [START_REF] Merkel | Manning[END_REF]13]. For each strain pair (ε, γ), we first minimized the network energy to find the corresponding athermal state before running for each temperature 10 different realizations of the MC simulation. In a slight deviation from the analytical part, in the simulations, we defined the isotropic strain such that V = (1+ε) D V * . The difference is only of higher-order in ε, and thus not relevant here.

For the MC simulation, we used the Metropolis-Hastings algorithm [START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications[END_REF] with 10 10 MC steps, where in each step we moved a single, randomly selected node. For each temperature and strain, we chose the constant node displacement such that the transition probability was ≈ 0.5. We computed the network stress tensor as the time-averaged force dipole density of the individual springs [42] and subtracted the ideal gas contribution [22,43,[START_REF] Tuckerman | Statistical Mechanics: Theory and Molecular Simulation[END_REF]. The shear modulus was computed from the shear stress using the central finite difference quotient with a shear strain step of ∆γ = 2.5 × 10 -3 .

FIG. 2 .

 2 FIG. 2.(a) Schematic: Variation of the accessible phase space volume Ω with isotropic strain ε in the limit of infinitely stiff springs. Starting at the athermal transition point at ε = 0, upon isotropic compression (ε < 0), Ω is expected to increase with the distance |ε| to the athermal transition point. The abscissa represents the N dof -dimensional configuration space. (b) Illustration of a simple two-spring system in 3D of linear dimension3 √ V = e ε 3 √ V * , where the black nodes have fixed positions and the blue node is freely movable. In panels a and b, blue dashed lines represent the accessible parts of the phase space.

FIG. 3 .

 3 FIG. 3. Numerical results for zero shear strain, γ = 0. (a) We find three scaling regimes. (a inset) For ε = 0, we numerically confirm the scaling t ∼ T 1/2 . The red dashed line is a fit to the curve t = √ κEκST for T ≤ 10 -6 to extract the value of √ κEκS. (b,c) When plotting isotropic tension t (b) or shear modulus G (c) over isotropic strain ε, a rescaling by √ T leads to a collapse of the curves for the different temperatures. (respective inset) Same data without rescaling. All curves in panels b, c and their insets are predictions from Eqs. (11) and (12) using a single parameter set (κE, κS, bε)[START_REF]We extracted κE ≈ 0.161 from the tension data for the athermal limit using a linear fit. The value of κS/kB ≈ 0.0629 was then obtained from the fit in the inset to Fig. 3a. The parameter bε ≈ 0.835 describing the interaction between isotropic and shear strain was ob[END_REF].

4 FIG. 4 .

 44 FIG.4. Numerical results for varying shear strain, γ > 0, with fixed ε = -9.15 × 10 -5 . (a) When plotting isotropic tension t over the combined strain, ε + bεγ 2 , a rescaling by √ T leads to a collapse of the curves for the different temperatures. (a inset) Same data without rescaling, and on the x axis shear strain is plotted directly. (b) Shear modulus G over shear strain γ. All curves are predictions from Eqs. (11) and (12) using the same parameter set as in Fig.3.

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.


Cheng-Tai Lee
email: matthias.merkel@cnrs.fr

published or not. The documents may come

Partition sum of thermal, under-constrained systems

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 (Dated: April 17, 2023) 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 . 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, underconstrained 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 [START_REF] Merkel | Manning[END_REF]4]. However, such a framework does not yet exist for thermal, i.e. finite-temperature, under-constrained systems.

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 , (1)
where N SSS is the number of linearly independent states of self-stress (SSS). A SSS is a combination of virtual tensions put on the springs of a network that does not result in any net force on the dofs. According to Eq. ( 1), under-constrained systems always have a positive number of infinitesimal zero modes.

There are three possibilities regarding the nature of infinitesimal zero modes. First, they could be trivial zero modes, i.e. global translations or rotations of the system, which leave the "shape" of the system unchanged. Second, they could correspond to mechanisms, i.e. finite displacements that change the "shape" of the system, but do not affect any of the constraints [5,6]. Third, while inducing no changes to the constraints to linear order, they could induce changes to higher order, corresponding to higher-order rigidity [7,8]. Here, we refer to this third kind as lower-order zero mode. At zero temperature, under-constrained are typically floppy because many of the infinitesimal zero modes often correspond to mechanisms.

Still, under-constrained systems can be rigidified by global strain, such as isotropic expansion or shear [2,[START_REF] Merkel | Manning[END_REF][7][8][9][10][11]. This athermal rigidity transition is due to the shrinkage of the hyper-volume of dof configurations that are geometrically compatible with the constraints. At this rigidity transition, a SSS is created [5], and a positive number of mechanisms turns into lower-order zero modes. Beyond this transition point, the constraints cannot be fulfilled any more, i.e. the external strain is geometrically incompatible with the constraints, and in the case of a spring network, any further strain creates prestresses in the springs. These prestresses, together with the SSS created at the transition, rigidify the system [2,[START_REF] Merkel | Manning[END_REF]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 [START_REF] Merkel | Manning[END_REF]4]. These ideas apply not only to spring networks, but also to other systems such as polymer networks [14][15][16][17] and vertex models for biological tissues [START_REF] Merkel | Manning[END_REF]13,[18][19][20].

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

Here, we analytically derive from first principles the partition sum and elastic properties of thermal, underconstrained systems close to the athermal transition for varying isotropic strain ε, shear strain γ, and small temperatures T . After defining the class of underconstrained systems that we discuss here (section II), we first revisit the athermal limit (section III). This will provide the foundations for discussing at finite temperature both the limit of infinitely stiff springs (section IV) and the general case of finite spring stiffness (section V). In the companion paper [26], we numerically verify the analytical results obtained here.

II. MODEL

A. Under-constrained system

We consider a Hookean network with fixed connectivity of N s generalized springs with length L i , labeled by i = 1, . . . , N s . The system energy is:

E = 1 2 Ns i=1 K i L i -L 0i 2 , (2)
where the sum is over all springs i with spring constants K i and rest lengths L 0i . Each generalized spring length L i has dimension d i . For example, angular springs have d i = 0, ordinary linear springs or perimeter springs in the vertex model have d i = 1, and area and volume springs in the vertex model have d i = 2 and d i = 3, respectively. Thus, also models that, for example, combine 2D or 3D vertex models with spring networks are accounted for [START_REF] Parker | How does the extracellular matrix affect the rigidity of an embedded spheroid?[END_REF]. The notion of the dimension d i of a generalized spring length L i is further formalized in Eq. ( 5) below.

The spring network is embedded in D spatial dimensions with periodic boundary conditions with total system volume V . The periodic box can be sheared with shear strain variable γ. Our framework applies both to simple and pure shear [START_REF] Merkel | Manning[END_REF].

The generalized spring lengths depend on N node node positions R a with a = 1, . . . , N node , which correspond to the N dof = DN node degrees of freedom. We denote the collection of all components of all node positions by the N dof -dimensional vector #-R. Here we consider an underconstrained spring network, where N dof > N s . To be precise, we discuss here a slightly stricter class of underconstrained models that will be specified later.

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 γ, and temperature T .

B. Athermal transition

We consider a network that in the athermal limit transitions from floppy (zero shear modulus) to rigid (finite shear modulus) at some critical volume. For spring networks or vertex models, this is typically the case when the network percolates the periodic box (in the sense of connectivity percolation). For zero shear strain, γ = 0, we denote this critical volume by V * , where, consistent with earlier findings [START_REF] Merkel | Manning[END_REF]4], we assume here that the system is floppy for V < V * and rigid for V > V * .

For simplicity, we assume here and in the following that there are no states of self-stress (SSS) in the floppy regime. Also, we assume that only a single SSS forms at the transition, which is generally true for disordered networks.

C. Isotropic strain

We define isotropic strain ε with respect to the critical volume V * :

ε := 1 D log V V * . (3)
Thus, for γ = 0, the transition occurs at ε = 0. We introduce the dimensionless degrees of freedom (ddofs) #r as:

#-r := V -1/D #- R = e -ε (V * ) -1/D #- R. (4)
This allows us to formalize the notion of spring dimension. We say that a spring i has dimension d i if the spring length L i ( #-R, ε, γ) is homogeneous of degree d i with respect to an isotropic rescaling. In other words, L i can be written as:

L i ( #- R, ε, γ) = L 0i e diε i ( #-r , γ) (5)
with a function i ( #r , γ) that depends on the ddofs #r and shear strain γ, but not on isotropic strain ε. We call i the dimensionless spring length. At the transition, where

L i = L 0i , we have i = 1.
Finally, we define the dimensionless weights

w i := K i L 2 0i E 0 , (6)
where we have introduced the energy scale E 0 := ( i K i L 2 0i )/N s such that i w 2 i = N s . This allows us to rewrite the system energy as

E = E 0 2 Ns i=1 w 2 i e diε i -1 2 . (7)
To lowest order in ε, and with ∆ i = i -1, we thus have:

E = E 0 2 Ns i=1 w 2 i ∆ i + d i ε 2 . (8)
To obtain Eq. ( 8) from Eq. ( 7), we have factored out the prefactor e diε in front of the parentheses, which scale to lowest order as ∼ ε. The prefactor of e diε only adds higher-order terms and is thus neglected.

III. ATHERMAL LIMIT (T → 0)

While in previous work we derived the elastic system properties in the athermal limit based on a minimallength function [START_REF] Merkel | Manning[END_REF]4], we present here an alternative derivation, which will prepare our discussions in the subsequent sections.

We denote the values of the ddofs #r at the transition by #r * . For simplicity, we assume here that #r * is uniquely defined. Global translations do not change any of the arguments below, and we focus on networks where mechanisms do not have any influence on the SSS that forms at the transition.

We now expand the dimensionless spring lengths i around the transition point at #r * and for small γ. This is possible because each dimensionless spring length i is an analytical function of #r and γ [START_REF] See | the appendices of Refs[END_REF]. Denoting the components of #r by r n , we can write to second order in ∆r n := r n -r * n and γ:

w ˙i∆ ˙i = C in ∆r n + 1 2 M imn ∆r m ∆r n + B (1)
i γ + 1 2 B (2) i γ 2 . ( 9
)
Here and in the following, we imply summation over equal indices (except for indices with underdots, and in the absence of explicit sums). On the right-hand side, C in is a generalized compatibility matrix, and the matrix M imn is symmetric in the last two indices. We introduce the shear strain dependence through the coefficients B

(1/2) i . In general, matrices C in and M imn may depend on γ. However, including this dependency turns out not to be important to lowest order in our final results, both here and in the following sections (appendix B). There are no ε-dependent terms in Eq. ( 9), because the i functions are independent of ε.

We further simplify the problem of minimizing the system energy by performing a singular-value decomposition (SVD) on the compatibility matrix:

C in = Ns p=1 U ip s p V np , (10)
where each of U ip and V np is an orthogonal square matrix, and the singular values s p are sorted in decreasing order:

s 1 ≥ • • • ≥ s Ns .
Because at the transition a single SSS is created, we have s Ns-1 > 0 and s Ns = 0, where U iNs is the corresponding SSS. Furthermore, all V np with p = N s , . . . , N dof are infinitesimal zero modes.

Based on the SVD of the compatibility matrix, Eq. ( 10), we define changes of spring lengths and of degrees of freedom in the "eigen bases" of C in :

∆ ˜ p := i U ip w i ∆ i ( 11
)
∆r p := V np ∆r n . (12)
Thus, the expansion in Eq. ( 9) transforms into:

∆ ˜ p = s ṗ∆r ṗ + 1 2 Mpmn ∆r m ∆r n + B(1) p γ+ 1 2 B(2) p γ 2 , ( 13
)
where we have also transformed M and B (1/2) into the eigenbases of C:

Mpmn := U ip M im n V m m V n n and B(1/2) p := U ip B (1/2) i
.

We now insert the Taylor expansion of w ˙i∆ ˙i into the system energy, Eq. ( 8). First, using Parseval's theorem, we have:

E = E 0 2 Ns i=1 w ˙i∆ ˙i + d ˙iw ˙iε 2 = E 0 2 Ns p=1 ∆ ˜ p + wp ε 2 . (14)
Here, we have introduced wp := i U ip d i w i . Note that Eq. ( 10) is invariant with respect to a sign flip in the SSS, whose components are given by U iNs for each spring i.

For convenience, we choose this sign such that wNs = i U iNs d i w i ≥ 0. Inserting Eq. ( 13) into the energy, Eq. ( 14), we obtain:

E = E 0 2 Ns p=1 G 2 p . (15)
with

G p := ∆ ˜ p + wp ε = s ṗ∆r ṗ + 1 2 Mpmn ∆r m ∆r n + wp ε + B(1) p γ + 1 2 B(2) p γ 2 . (16)
To discuss the energy minimum we focus here for clarity on the case where B(1) p = 0 for all p. The general case of arbitrary B(1) p is discussed in appendix B. As we discuss in appendix A, at the energy minimum, the ∆r n do not scale in the same way with the strain variables ε and γ for different n. For p < N s we find that ∆r p is of order of at least O(ε) + O(γ 2 ), and for n ≥ N s we find that ∆r n is of order of at least O(|ε| 1/2 ) + O(γ). Consequentially, all of the M terms involving any of the ∆r p with p < N s can be neglected in all the G q with q ≤ N s . Thus, each of the ∆r p with p < N s only appears in total once, namely as linear term in the respective G p . Thus, to the lowest relevant order, when minimizing the energy in Eq. ( 15), each of the first N s -1 terms G p for p < N s can become zero by adjusting ∆r p accordingly:

∆r p = - 1 s ṗ 1 2 N dof m,n=Ns M˙p mn ∆r m ∆r n + wp ε + 1 2 B(2) ṗ γ 2 . ( 17
)
The remaining term is:

G Ns = 1 2 N dof m,n=Ns MNsmn ∆r m ∆r n + wNs ε + 1 2 B(2) Ns γ 2 . ( 18
)
We know that for γ = 0 and ε > 0, where the system is in the rigid regime, the energy minimum is strictly positive, E = E 0 G 2 Ns /2 > 0. This is only possible if the matrix MNsmn is positive semi-definite; otherwise G Ns could become zero. As a consequence, for the athermal rigid regime, ε+b ε γ 2 > 0, the energy is minimal if ∆r n = 0 for all N s ≤ n ≤ N dof , and it has the value

E = E 0 w2 Ns 2 ε + b ε γ 2 2 , (19)
where we introduced

b ε := B (2) Ns 2 wNs if B(1) p = 0 for all p. (20)
In the more general case where the B(1)

p
do not vanish, the definition of b ε includes more terms, while Eq. ( 19) remains the same up to a shift in γ (appendix B).

Note that the matrix MNsmn is closely linked to the Hessian of the system and the emergence of rigidity in the athermal limit (appendix C).

B. Comparison to earlier results

The published expression for the energy of a network with only d-dimensional springs in the rigid regime corresponds to [START_REF] Merkel | Manning[END_REF]4]:

E = E 0 N s 2(L * 0 ) 2 (1 + a 2 ) -∆L 0 + bγ 2 2 , (21)
where ∆L 0 = -dL * 0 ε, and L * 0 , a , 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 [START_REF] Merkel | Manning[END_REF]4]. Comparing Eqs. ( 19) and ( 21), we find:

wNs = d N s 1 + a 2 1/2 (22) b ε = b dL * 0 . (23)
In addition to the earlier work [START_REF] Merkel | Manning[END_REF]4], our approach here also clarifies why the result is analytic in γ; it is inherited from the analytic nature of the i functions (Eq. ( 9)).

C. Elastic properties

In the rigid regime of the athermal limit, i.e. for ε + b ε γ 2 > 0, tension t E = ∂E/∂V = (∂E/∂ε)/DV and shear modulus G E = (∂ 2 E/∂γ 2 )/V are to lowest order in ε and γ:

t E = κ E ε + b ε γ 2 (24)
G E = 2Db ε κ E ε + 3b ε γ 2 (25) = 2Db ε t 1 + 2b ε γ 2 |ε + b ε γ 2 | , (26)
where we defined

κ E := E 0 w2 Ns DV * . (27)
The index E indicates purely energetic rigidity.

IV. STIFF-SPRING LIMIT (Ki → ∞)

In the limit of infinitely stiff springs but at finite temperature T , elasticity is created purely by entropic effects. For ε+b ε γ 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 ε + b ε γ 2 ≤ 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 Ω, from which we can directly derive the elastic network properties. We focus again on the limit of small ε and γ, in which the network is almost rigid [START_REF] Holmes-Cerfon | [END_REF].

A. Configurational phase space volume

Up to a constant factor, the configurational phase space volume is given by (appendix D)

Ω ∼ N node a=1 d D R a Ns i=1 δ(L i -L 0i ), (28)
which, using the rescaling from Eqs. ( 4) and ( 5), transforms into

Ω ∼ N dof n=1 dr n Ns i=1 δ(e diε i -1), (29)
where we ignore any prefactors in Ω that are powers of e ε , because these factors only create higher-order terms in the final result. We aim to obtain an analytical expression for Ω in the vicinity of the transition point, i.e. for |ε| 1 and |γ| 1. We can thus rewrite Ω as:

Ω ∼ N dof n=1 d∆r n Ns i=1 δ w i ∆ i + d i w i ε . (30)
Here, we have included a constant factor of w i in each of the Dirac deltas, and we again ignored prefactors of e ε in Ω.

To insert the Taylor expansion, Eq. ( 13), into Eq. ( 30), we first apply a rotation by V np to the integration variables ∆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 Ω:

Ω ∼ N dof n=1 d∆r n Ns p=1 δ ∆ ˜ p + wp ε = N dof n=1 d∆r n Ns p=1 δ G p . (31)
The constraints by the Dirac deltas in this equations imply G p = 0 for all p ≤ N s .

In the case where all B

(1) p = 0, we can again use the same arguments as in appendix A to show that all terms involving the ∆r p with p < N s can be neglected in all G q with q ≤ N s . Thus, like in the athermal limit, the value of each of the ∆r p with p < N s is adjusted by the respective equation G p = 0, and the remaining G Ns simplifies to Eq. (18). In other words, we can split the integrations in Eq. ( 31) as follows:

Ω ∼ N dof n=Ns d∆r n δ(G Ns ) Ns-1 p=1 d∆r p δ(G p ). (32)
Here, each of the N -1 inner integrals evaluates to a constant factor of 1/s p . We are thus left with:

Ω ∼ Ω Ns (ε, γ) (33)
with

Ω Ns (ε, γ) := N dof n=Ns d∆r n δ G Ns = N dof n=Ns d∆r n × δ 1 2 N dof m,n=Ns
MNsmn ∆r m ∆r n

+ wNs ε + b ε γ 2 . ( 34
)
To evaluate Ω Ns (ε, γ), we first diagonalize the matrix MNsmn :

MNsmn = N dof -Ns+1 k=1 µ k v k m v k n , (35)
with m, n = N s , . . . , N dof , and where µ k are the eigenvalues and v k m are the associated orthonormal eigenvectors. We sort the eigenvalues in decreasing order,

µ 1 ≥ • • • ≥ µ N dof -Ns+1
, where µ k ≥ 0, since the matrix MNsmn is positive semi-definite.

With ∆q k := v k m ∆r m , the equation G Ns = 0 transforms into:

-wNs ε + b ε γ 2 = 1 2 N dof -Ns+1 k=1 µ k ∆q 2 k . ( 36
)
Because the v k m represent a rotation, the differential transforms without prefactor:

Ω Ns = N dof -Ns+1 k=1 d∆q k × δ 1 2 N dof -Ns+1 k=1 µ k ∆q 2 k + wNs ε + b ε γ 2 .
(37)

To evaluate this integral, we first discuss the zero modes of MNsmn , i.e. modes k with µ k = 0, for which there are two possibilities: First, such a mode k could be one of the D global translations, or it could be related to a mechanism. In the first case, the integral over a global translation leads to a factor of e ε , which can again be neglected in the end result. Similarly, the integral over the manifold corresponding to the mechanisms that exists at the transition yields a constant (describing the hypervolume of that manifold) plus, possibly, a higher-order term in ε and/or γ (describing the change of the manifold volume upon strain), which can be neglected. Another possibility is however that any change of the zero mode ∆q k does affect ∆ Ns , but only to higher than quadratic order. We believe this is an exceptional case and do not consider it further in the main text, but we discuss it in appendix E.

We denote the number of non-zero eigen values of MNsmn by N e . Integrating out the zero modes in Eq. (37), we have:

Ω Ns ∼ Ne k=1 d∆q k × δ 1 2 Ne k=1 µ k ∆q 2 k + wNs ε + b ε γ 2 . ( 38
)
The Dirac delta ensures that the integral is over the surface ∂E of an ellispoid E in N e -dimensional space, defined by the equation:

Ne k=1 ∆q 2 k ρ 2 k = 1 ( 39
)
with half axes

ρ k := (-2 wNs [ε + b ε γ 2 ]/µ k ) 1/2 ∼ [-(ε + b ε γ 2 )] 1/2 .
Using standard transformation rules of the Dirac delta, we get:

Ω Ns ∼ ∂E dS ∂G Ns ∂∆q -1 , (40)
where

∂G Ns ∂∆q 2 = Ne k=1 ∂G Ns ∂∆q k 2 = Ne k=1 µ 2 k ∆q 2 k . (41)
With µ max := µ 1 and µ min := µ Ne > 0, we thus have on the ellipsoid surface:

-µ min ε + b ε γ 2 ≤ 1 2 wNs ∂G Ns ∂∆q 2 ≤ -µ max ε + b ε γ 2 .
(42) Thus, |∂G Ns /∂∆q| ∼ [-(ε+b ε γ 2 )] 1/2 . As a consequence, we obtain that the integral in Eq. ( 40) scales as the hypersurface area of the ellipsoid divided by

[-(ε + b ε γ 2 )] 1/2 : Ω Ns (ε, γ) ∼ -(ε + b ε γ 2 ) (Ne-2)/2 . ( 43
)
Hence, Ω ∼ Ω Ns scales as a power law with the distance to the transition point as long as N e > 2.

In appendix B, we discuss the more general case where the B(1) p do not vanish. The definition of b ε changes like in the athermal case, while the result Eq. ( 43) remains the same up to an additional term B

(1) Ns γ, which can be removed through redefining γ.

B. Free energy, tension, and shear modulus

The free energy of the system in the stiff-spring limit, F S , is up to a constant:

F S = - 1 β log Ω = - N e -2 2β log -(ε + b ε γ 2 ) , (44)
where β := 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 = (∂F S /∂ε)/DV and the shear modulus G S = (∂ 2 F S /∂γ 2 )/V , to lowest order:

t S = - κ S T ε + b ε γ 2 (45) G S = -2Db ε κ S T ε -b ε γ 2 (ε + b ε γ 2 ) 2 (46) = 2Db ε t 1 + 2b ε γ 2 |ε + b ε γ 2 | ( 47
)
with

κ S := k B (N e -2) 2DV * . ( 48
)
Note that only the regime with ε + b ε γ 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 = N node a=1 d D R a e -βE . (49)
In terms of rotated dimensionless node positions relative to the transition point, ∆r n , we have:

Z ∼ N dof k=1 d∆r n e -βE , (50)
where we again ignore prefactors proportional to e ε . We evaluate the partition sum for small ε, γ, and T . Because ε 1, we can use the energy expression from Eq. ( 15):

E = E 0 2 Ns p=1 G 2 p . (51)
Here we discuss the case where all B(1) p = 0, while the more general case of arbitrary B(1)

p is discussed in ap- pendix B.
We focus on small temperatures, k B T E 0 [ε+b ε γ 2 ] 2 , such that states with smaller energies E receive an exponentially higher Boltzmann weight. Using our discussion in appendix A, we find again that in all G q with q ≤ N s all the ∆r p with p < N s can be neglected in the quadratic terms. Thus, to lowest order each of the ∆r p with p < N s appears only once, namely in the linear term of in the respective G p . As a consequence, the integrals in the partition sum in Eq. (50) can be rearranged as follows:

Z ∼ N dof n=Ns d∆r n e - βE 0 G 2 Ns 2 Ns-1 p=1 d∆r p e - βE 0 G 2 p 2 . ( 52
)
The innermost first N s -1 integrals are Gaussian integrals, each of which evaluates to the constant Z p = 2π/βE 0 s 2 p . When combined with the thermal wavelength of the nodes, each of the factors Z p corresponds to the partition sum of a harmonic oscillator with eigen frequency ω p = s p E 0 /m with m being the mass of the nodes. 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 Ns := N dof n=Ns d∆r n e - βE 0 G 2 Ns 2 . ( 53
)
Neglecting the factors Z p with p < N s , which solely depend on temperature, but not on strain, we have Z ∼ Z Ns .

B. Saddle point approximation

To simplify the remaining integral Z Ns in Eq. ( 53), we introduce in its integrand an additional integral over some variable ε S that trivially evaluates to one:

Z Ns = N dof n=Ns d∆r n e - βE 0 G 2 Ns 2 × wNs ∞ -∞ dε S δ G Ns + wNs [ε S -ε] . (54)
After exchanging the order of integration and substituting Ω Ns , defined in Eq. (34), we have:

Z Ns = wNs ∞ -∞ dε S Ω Ns (ε S , γ)e -βE Ns (ε-ε S ) , (55)
with:

E Ns (ε -ε S ) := E 0 w2 Ns 2 ε -ε S 2 = DV * κ E 2 ε -ε S 2 . (56)
Here, we used the definition for κ E in Eq. ( 27). From Eq. (43) in the previous section, we have for ε S +b ε γ 2 ≤ 0:

Ω Ns (ε S , γ) ∼ -(ε S + b ε γ 2 ) (Ne-2)/2 , (57)
while for ε S + b ε γ 2 > 0, no configurations are possible, Ω Ns (ε S , γ) = 0. These transformations identify ε S as the entropic strain discussed in the companion paper [26]. We use a saddle point approximation to simplify the integral in Eq. (55):

Z Ns = wNs ∞ -∞ dε S e -β F , ( 58
)
where F is the free energy of the system for imposed ε S :

F (ε S ; ε, γ) := E Ns (ε -ε S ) - 1 β log Ω Ns (ε S , γ) = DV * κ E 2 ε -ε S 2 -κ S T log -(ε S + b ε γ 2 ) (59)
On the second line, we used the definition of κ S , Eq. ( 48), and ignored a constant offset.

For small temperatures, we can apply the saddle point approximation, where we Taylor expand F (ε S ) to second order around its minimum at ε S = εS , which results in:

Z Ns ∼ ∞ -∞ dε S e -β F = e -β F (ε S ) 2π β F (ε S ) , (60)
where F := ∂ 2 F/∂ε 2 S :

F (ε S ) = DV * κ E + κ S T [ε S + b ε γ 2 ] 2 . ( 61
)
To find the minimum εS , we transform ∂ F/∂ε S = 0, and obtain:

εS + b ε γ 2 = ε + b ε γ 2 2 - √ 1 + θ ε + b ε γ 2 2 , (62)
with

θ := 4κ S T κ E (ε + b ε γ 2 ) 2 = 2(N e -2)k B T E 0 w2 Ns (ε + b ε γ 2 ) 2 . ( 63
)
Note that for γ = 0, the equation ∂ F (ε S )/∂ε S = 0, defining the minimum εS of F (ε S ) = E(ε -ε S ) + F S (ε 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 higher-order terms in the strain variables, the free energy of the system, F = -(log Z)/β, is:

F = F (ε S ) + 1 2β log F (ε S ). ( 64
)
After some transformations, we find for the isotropic tension t = (∂F/∂ε)/DV and the shear modulus G = (∂ 2 F/∂γ 2 )/V , to lowest order in ε and γ (appendix F):

t = κ E ε -εS 1 + 1 2(N e -2) θ 1 + θ (65)
and

G = 2Db ε κ E ε -εS 1 + 2b ε γ 2 |ε + b ε γ 2 | √ 1 + θ + 1 2(N e -2) θ 1 + θ 1 + 2b ε γ 2 |ε + b ε γ 2 | √ 1 + θ - 4b ε γ 2 (ε + b ε γ 2 )(1 + θ) . (66)
Note that the last terms in both t and G scale as ∼ (N e -2) -1 . These terms stem from the term ∼ log F in the free energy. For N e 1, these contributions can be neglected, and we have:

t = κ E ε -εS (67) G = 2Db ε κ E ε -εS 1 + 2b ε γ 2 |ε + b ε γ 2 | √ 1 + θ (68) = 2Db ε t 1 + 2b ε γ 2 |ε + b ε γ 2 | √ 1 + θ . ( 69
) D. Limits 1. Small temperature, θ 1
For θ 1, Eq. ( 62) becomes:

εS + b ε γ 2 = ε + b ε γ 2 2 -1 + θ 2 ε + b ε γ 2 2 = (ε + b ε γ 2 ) × 1 + θ/4 for ε + b ε γ 2 < 0, -θ/4 for ε + b ε γ 2 > 0. ( 70
)
Thus, to lowest order in θ:

ε -εS = (ε + b ε γ 2 ) × -θ/4 for ε + b ε γ 2 < 0, 1 for ε + b ε γ 2 > 0, (71)
and we recover the purely entropic and energetic tensions and shear moduli, respectively (compare Eqs. ( 24), ( 26), (45), and (47)):

t = t S for ε + b ε γ 2 < 0, t E for ε + b ε γ 2 > 0, (72)
G = G S for ε + b ε γ 2 < 0, G E for ε + b ε γ 2 > 0. ( 73
)
2. Small strain, θ 1

For θ 1, Eq. ( 62) becomes:

εS + b ε γ 2 = - κ S T κ E + ε + b ε γ 2 2 (74)
And thus:

ε -εS = κ S T κ E + ε + b ε γ 2 2 . ( 75
)
We then find:

t = 2N e -3 2N e -4 κ E κ S T + κ E ε + b ε γ 2 2 , ( 76
)
and for the shear modulus:

G = Db ε 2 κ E κ S T + κ E (ε + b ε γ 2 ) × 2N e -3 2N e -4 1 + b ε γ 2 κ E κ S T - b ε γ 2 2N e -4 κ E (ε + b ε γ 2 ) κ S T . (77)
For N e 1 both expressions simplify to:

t = κ E κ S T + κ E ε + b ε γ 2 2 (78) G = 2Db ε t 1 + b ε γ 2 κ E κ S T . ( 79
)
Note that we assumed in our derivation for the general case that

k B T E 0 [ε+b ε γ 2 ] 2 .
Using the definitions of θ in Eq. ( 63), of κ E in Eq. ( 27), and of κ S in Eq. ( 48), this condition corresponds to θ 2(N e -2)/ w2 Ns . Thus, this limit corresponds to 1 θ 2(N e -2)/ w2 Ns , which can be fulfilled for large systems. Note that in our numerical results for randomly-cut triangular networks, we also found that our prediction was accurate for finite temperatures even for vanishing strain, ε + b ε γ 2 → 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 elastic material properties such as tension t and shear modulus G close to the athermal rigidity transition depending on isotropic strain ε, shear strain γ, and temperature T . The only three parameters are κ E , which characterizes energetic rigidity, κ S , which characterizes entropic rigidity, and b ε , which characterizes the interaction between isotropic and shear strain. These parameters are phenomenological coefficients that depend on the microscopic structure of the system.

We first discussed the limits of zero temperature (athermal limit), and infinitely stiff springs, before discussing the general case of finite temperature and finite spring stiffness. In the athermal limit, the elastic properties were analytically derived before [START_REF] Merkel | Manning[END_REF]4], where the system is floppy for ε+b ε γ 2 ≤ 0 and rigid for ε+b ε γ 2 > 0. We showed that in the stiff-spring limit, the phase space volume is given by the surface of a hyper-ellipsoid, which leads to a free energy that scales logarithmically with the combined strain ε + b ε γ 2 . In both cases, 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.

We introduced under-constrained system as systems that have N dof > N s , and our results apply to a subset of such systems. In particular, we consider systems that fulfill the following three criteria: (i) there is no SSS in the floppy regime, (ii) only a single SSS forms at the transition, and (iii) N e > 2. Here, N e is the number of non-zero eigen values of the MNsmn matrix. In other words, in the athermal limit N e is the number of eigen values of the Hessian that are zero in the floppy regime and become finite in the rigid regime (appendix C). Lifting the restriction on the number of SSS in the floppy regime probably does not change our results. Similarly, we do not expect the number of SSS that form at the transition to substantially affect our results. Finally, it will be interesting to study the marginal cases of N e = 1 and N e = 2 in future work.

While not explicitly studied here, our results have immediate consequences with respect to the scaling with network connectivity z. For γ = 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 ε is replaced by the relative connectivity ∆z = z -z c , with z c being the connectivity at isostaticity [24]. To understand how this could be, we discuss the ∆z scaling that results from our findings by separately discussing the ∆z scaling for κ S and κ E . Eq. ( 48) predicts that κ S essentially only depends on N e , i.e. the number of non-zero eigen values of the matrix MNsmn with m, n ≥ N s . Thus, if for large systems N e scales to dominating order linearly with ∆z = -4(N dof -N s )/N dof , then this implies κ S ∼ |∆z| 1 . Meanwhile, how energetic elasticity κ E scales with connectivity ∆z depends on the class of the network [4]. For instance, for |∆z| 1 our results for the a scaling in Refs. [START_REF] Merkel | Manning[END_REF]4] imply that κ E ∼ |∆z| 2 for Delaunay networks, κ E ∼ |∆z| 1 for randomly-cut, packing-derived networks, and κ E ∼ |∆z| 0 for phantom triangular networks. Because of Eq. (67), these scaling relations of κ S and κ E with ∆z determine the scaling of tension t. However, for the scaling of the shear modulus G, Eq. ( 69 and b ε κ E ∼ |∆z| -0.5 for phantom triangular networks. Taken together, this suggests that the scaling of the elastic properties with ∆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 ε κ S ∼ |∆z| 0 and b ε κ E ∼ |∆z| 0 like randomlycut, packing-derived networks. In this case, the scaling is entirely determined by the scaling of the strain of the network state after network initialization, ε ini , with con-nectivity ∆z, which is generally linear to lowest order: ε ini ∼ ∆z (remember that we define isotropic strain here with respect to the athermal transition point) [START_REF] Merkel | Manning[END_REF]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 ε by ∆z. We leave a more thorough numerical test of all these predicitons for the scaling with ∆z, and a study of their consequences, for future work.

In general, there can be a degeneracy of the possible network configurations #r * at the transition, and here we have assumed that this does not affect our result. To this end, we needed to explicitly restrict ourselves to systems where the choice of the transition configuration #r * does not affect the properties of the SSS. In particular, while properties like the other singular values, s p with p < N s , of the compatibility matrix C in may depend on #r * [6,[START_REF] Rocklin | Folding mechanisms at finite temperature[END_REF], our results should be unaffected as long as the SSS (more precisely wNs ) and its behavior with respect to shear (more precisely b ε ) are independent of the choice of #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.

Previous work discussing the athermal limit was based on a minimal-length function [START_REF] Merkel | Manning[END_REF]4], where a Taylor expansion with a quadratic dependence on γ relied on the assumption that the minimal-length function was analytic in γ. With the approach here we now that the dependence on γ is analytic, because we know that the functions i ( #r , γ) are analytic. However, it is in principle still possible that higher-order terms could create isolated discontinuities in the emergent behavior of the system.

Finally, we note that all results obtained here are based on Taylor expansions for systems of a given finite size. As discussed in Ref. [4], the strain ranges over which these predictions hold may depend on system size. In particular, depending on the type of network studied, the extent of these ranges may decrease with increasing system size, and beyond this strain range, different scaling regimes may appear [START_REF] Arzash | [END_REF]. B(1)

+ wNs ε + b ε γ2 , (
p s p MNspm v k m v k n µ k MNsnq B(1) q s q       . ( B8
)
We assumed here that for all k with µ k = 0 the relation n≥Ns,p<Ns v k n MNsnp B(1) p /s p = 0. If the second-order zero mode k stems from a mechanism or a trivial zero mode (i.e. linked to a global translation), this relation is necessarily true; otherwise there would be no energy minimum. If the mode k is neither mechanism nor trivial zero mode, but a higher-than-second-order zero mode, then higher-order terms need to be taken into account. However, we regard this as an exceptional case.

The additional terms in b ε in Eq. (B8) as compared to Eq. ( 20) scale as ∼ B(1) p B(1) q . The symmetry of the system with respect to the sign of γ combined with the central limit theorem suggest that, in general, these linear terms scale with system size as B( 1)

p ∼ N -1/2 s
. Indeed, we numerically found earlier for phanomized triangular networks that B(1)

Ns ∼ ∆γ ∼ N -1/2 dof ∼ N -1/2 s [4]
. We thus expect that the additional contributions in b ε scale as N -1 s . For the stiff-spring limit and the general case of finite spring stiffness and temperature, these arguments are essentially the same. The only difference in the stiff-spring limit is that we immediately have G p = 0. Moreover, for non-vanishing B(1) p , the ellipsoid that describes the possible system configurations is now centered at (∆q 0 k ) (compare Eqs. ( 36) and (B6)). According to Eq. (B7) this ellipsoid center displaces linearly with shear strain γ. In the general case, we can follow the same arguments as in the athermal limit, using that the terms where the G p with p < N s are of order O(ε) + O(γ 2 ) or lower have an exponentially reduced Boltzmann weight and can thus be neglected in the partition sum.

In our Taylor expansion, Eq. ( 9), we have left out terms that are linear in both ∆r n and γ. We expect such terms to also scale with system size as ∼ N -1/2 s . They create additional contributions to ∆q 0 k and, in interaction with the B(1) p terms, contributions to b ε . However, other than adding a correction to the value of b ε they do not affect the results presented in the main text. We evaluate the Hessian of the system,

H mn := ∂ 2 E ∂∆r m ∂∆r n , ( C1
)
in the rigid athermal regime, at an energy minimum close to the transition point, where we assume for simplicity that B(1) p = 0 for all p ≤ N s . Using Eq. ( 15):

H mn = E 0 Ns p=1 ∂G p ∂∆r m ∂G p ∂∆r n + G p ∂ 2 G p ∂∆r m ∂∆r n . (C2)
Further, according to our discussion in section III A, at an energy minimum, ∆r p is at least of order O(ε)+O(γ 2 ) for p < N s , and ∆r n = 0 for n ≥ N s . Also, G p = 0 for p < N s up to order O(ε) + O(γ 2 ). Thus, we have up to order O(ε) + O(γ 2 ):

H mn = E 0 Ns-1 p=1 s 2 p δ pm δ pn + E 0 Ns-1 p,q=1
s p ∆r q δ pm Mpqn + δ pn Mpqm + E 0 G Ns MNsmn .

(C3)

In this expression, the first term, which is of absolute order in ε and γ, corresponds to the Hessian at the transition point. The second and third term are of order O(ε) + O(γ 2 ). Eq. (C3) shows that to lowest order, the N e nonzero eigen modes of MNsmn for m, n ≥ N s , are those eigen modes of H mn that have zero eigen value at the transition point, and whose eigen values become nonzero in the rigid regime, scaling as ∼ O(ε) + O(γ 2 ). This is because for m, n ≥ N s only the third term in Eq. (C3) can be nonzero.

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

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

Z = N node a=1 d D R a exp - β 2 Ns i=1 K i L i -L 0i 2 .
(D1) A Gaussian converges to a Dirac Delta in the limit of small standard deviation:

lim Ki→∞ βK i 2π exp - βK i (L i -L 0i ) 2 2 = δ(L i -L 0i ).
(D2) Thus, whenever all K i are large enough such that the (βK i ) -1/2 becomes much smaller than all relevant length scales of the system, the partition sum can be expressed as:

Z = 2π β Ns/2 Ns i=1 K i -1/2 × N node a=1 d D R a Ns i=1 δ L i -L 0i . (D3)
This indicates that the correct approach is to use Dirac deltas on the spring lengths rather than the energy E.

Note that in Eq. ( 29), 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 → ∞. However, to discuss the elastic system properties, we only need derivatives of log Ω, where this term only creates a constant albeit diverging offset.

Appendix E: Higher-order dependencies of ∆ Ns on the dofs

Here we outline how our arguments change in the presence of eigen vectors with eigen value zero of the matrix MNsmn with m, n ≥ N s , but which affect ∆ ˜ Ns to some order α > 2. We first note that the exponent α needs to be even, because otherwise, such a mode would allow for states with G Ns = 0 for ε > 0 and γ = 0, even though the system is supposed ot be rigid with E > 0 in this regime.

Setting again for simplicity B(1) p = 0, we discuss the following generalization of the expansion of ∆ ˜ Ns from Eqs. ( 13), (36):

∆ ˜ Ns = α 1 α! N (α) e k=1 µ α,k ∆q α α,k + 1 2 B(2) Ns γ 2 , (E1)
where the sum over α is over all positive even integers, α! denotes the faculty of α, the variable N (α) e denotes the number of modes ∆q α,k that enter ∆ ˜ Ns to lowest order as ∼ ∆q α α,k , and µ α,k are the associated coefficients. Thus, the phase space volume is governed by an integral of the following kind:

Ω Ns ∼ I(ε, γ) :=   α N (α) e k=1 d∆q α,k   × δ   α 1 α! N (α) e k=1 µ α,k ∆q α α,k + wNs ε + b ε γ 2   , (E2)
To extract the dependency of I on ε and γ, we first define:

I 1 :=   α N (α) e k=1 dx α,k   × δ   α 1 α! N (α) e k=1 µ α,k x α α,k + wNs   , (E3)
which depends neither on ε nor on γ. Now we can transform the expression in Eq. (E2) using the following substitutions for all integration variables:

∆q α,k = ε + b ε γ 2 1/α x α,k , (E4)
which yields:

Ω Ns (ε, γ) ∼ I(ε, γ) = ε + b ε γ 2 NΩ I 1 (E5)
with:

N Ω = -1 + α N (α) e α . ( E6
)
The first term in N Ω arises from removing the factor [ε + b ε γ 2 ] from the Dirac delta, and all the remaining terms arise from the substitution of the integration variables. Note that for the case discussed in the main text, i.e. for N e := N = 0 for α > 2, we have indeed N Ω = (N e -2)/2. Hence, also higher-order dependencies of ∆ ˜ Ns on the dofs also contribute to the exponent N Ω , and thus to κ S = N Ω k B /DV * ; they are just attenuated according to their order α.

In this argument, we have excluded modes that leave ∆ ˜ Ns invariant. They do not contribute to N Ω . This is because any substitution like in Eq. (E4), 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 ∆ ˜ Ns , 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.

Appendix F: Tension and shear modulus in the general case

To obtain the isotropic tension t := (∂F/∂ε)/DV , we first use the chain rule. To lowest order, V = V * , and so we obtain from Eq. (64):

t = 1 DV * ∂ F ∂ εS ∂ εS ∂ε + ∂ F ∂ε + 1 2β F ∂ F ∂ εS ∂ εS ∂ε + ∂ F ∂ε .
(F1)

We compute to lowest order:

∂ εS ∂ε = - εS + b ε γ 2 |ε + b ε γ 2 | √ 1 + θ (F2) 1 DV * ∂ F ∂ εS = -κ E (ε -εS ) - κ S T εS + b ε γ 2 = 0 (F3) 1 DV * ∂ F ∂ε = κ E (ε -εS ) (F4) 1 DV * ∂ F ∂ εS = - 2κ S T (ε S + b ε γ 2 ) 3 (F5) 1 DV * ∂ F ∂ε = 0. (F6)
We thus obtain:

t = κ E (ε -εS ) 1 + 1 2(N e -2) θ 1 + θ , ( F7
)
where we have also used that

F = -DV * κ E ε + b ε γ 2 √ 1 + θ εS + b ε γ 2 (F8)
and

κ E (ε -εS ) = - κ S T εS + b ε γ 2 , ( F9
)
which follows from Eq.( 62). The term ∼ (N e -2) -1 in Eq. (F7) corresponds to the terms involving F . For the shear modulus G := (∂ 2 F/∂γ 2 )/V , we have to lowest order:

G = 1 V * ∂ 2 F (ε S (ε, γ); ε, γ) ∂γ 2 + 1 2β
∂ 2 ∂γ 2 log F (ε S (ε, γ); ε, γ) .

(F10)

We compute with Eqs. (F3) and (F9):

∂ εS ∂γ = -b ε γ 1 + ε + b ε γ 2 |ε + b ε γ 2 | √ 1 + θ (F14) = -2b ε γ ε -εS |ε + b ε γ 2 | √ 1 + θ (F15) 1 V * ∂ 2 F (ε S ; ε, γ) ∂γ 2 = 2Dκ E (ε -εS )b ε . (F16)
Thus:

1 V * ∂ 2 F (ε S (ε, γ); ε, γ) ∂γ 2 = 2Dκ E (ε -εS )b ε 1 + 2b ε γ 2 |ε + b ε γ 2 | √ 1 + θ .
(F17) Furthermore:

∂ F (ε S (ε, γ); ε, γ) ∂γ

= ∂ F (ε S ; ε, γ) ∂ εS ∂ εS ∂γ + ∂ F (ε S ; ε, γ) ∂γ = DV * κ E 1 εS + b ε γ 2 - εS -ε (ε S + b ε γ 2 ) 2 × -2b ε γ ε -εS |ε + b ε γ 2 | √ 1 + θ -2b ε γ εS -ε (ε S + b ε γ 2 ) 2 = -2DV * κ E b ε γ ε -εS (ε S + b ε γ 2 ) 2 |ε + b ε γ 2 | √ 1 + θ × ε + b ε γ 2 -|ε + b ε γ 2 | √ 1 + θ = -4DV * κ E b ε γ ε -εS (ε S + b ε γ 2 )|ε + b ε γ 2 | √ 1 + θ , (F18)
and thus:

∂ ∂γ log F (ε S (ε, γ); ε, γ) = 4b ε γ ε -εS (ε + b ε γ 2 ) 2 (1 + θ) . (F19)
The second derivative is:

∂ 2 ∂γ 2 log F (ε S (ε, γ); ε, γ) = 4b ε (ε -εS ) (ε + b ε γ 2 ) 2 (1 + θ) + 4b ε γ (ε + b ε γ 2 ) 2 (1 + θ) × 2b ε γ ε -εS |ε + b ε γ 2 | √ 1 + θ - 16b 2 ε γ 2 (ε -εS ) (ε + b ε γ 2 ) 3 (1 + θ) 2 = 4b ε (ε -εS ) (ε + b ε γ 2 ) 2 (1 + θ) 1 + 2b ε γ 2 |ε + b ε γ 2 | √ 1 + θ - 4b ε γ 2 (ε + b ε γ 2 )(1 + θ) .
(F20)

From Eqs. (48) and (F9) follows:

1 2βV * = Dκ S T N e -2 = Dκ E N e -2 θ (ε + b ε γ 2 ) 2 4 (F21)
Combining Eqs. (F10), (F17), (F20), and (F21), we obtain Eq. (66) in the main text.

  ), the scalings of b ε κ S and b ε κ E are the relevant ones, and thus the scaling of b ε ∼ b needs to be taken into account. Based on earlier results for the scaling of b [3, 4], we have b ε κ S ∼ |∆z| -1 and b ε κ E ∼ |∆z| 0 for Delaunay networks, b ε κ S ∼ |∆z| 0 and b ε κ E ∼ |∆z| 0 for randomlycut, packing-derived networks, and b ε κ S ∼ |∆z| 0.

5

 5

  k:µ k >0

  Appendix C: Relation between the Hessian and the matrix MNsmn

∂∂ 2 F∂ 2 F∂ 2 F∂ 2 F

 2222 F (ε S (ε, γ); ε, γ) ∂γ = 1 V * ∂ F (ε S ; ε, γ) ∂γ = -2Dκ S T b ε γ εS + b ε γ 2 = 2Dκ E (ε -εS )b ε γ, (ε S (ε, γ); ε, γ) ∂γ 2 (ε S ; ε, γ) ∂γ∂ εS ∂ εS ∂γ + (ε S ; ε, γ) ∂γ 2 (ε S ; ε, γ) ∂γ∂ εS = -2Dκ E b ε γ (F13)

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. 1

Appendix A: Orders of magnitude of the ∆rn We discuss the scaling of the ∆r n with ε and γ in the athermal limit at the energy minimum, for which we will also have to discuss the scaling of the G p . We focus here for clarity on the case where B(1) p = 0 for all p. The general case of arbitrary B (1) p is discussed in appendix B. Using Eq. ( 15), the energy minimum must fulfill:

At the transition all ∆r n = 0, and so they must be of higher than absolute order in ε or γ. Moreover, let us choose p < N s such that none of the of G p with p < N s is of lower order in ε than G p. Because of Eq. (A1), the order of G p in ε is higher than that of the minimum order of all G q with q ≤ N s . Because the order of G p is minimal among the G p with p < N s , this implies that the order of G p is higher than that of G Ns . According to Eq. ( 18), G Ns is of order O(ε) or higher [32]. Thus, because the order of G p is minimal, all G p with p < N s are of higher order than O(ε). Analogously, one can argue that all G p with p < N s are of higher order than O(γ 2 ). Thus, up to order [O(ε) + O(γ 2 )] 2 , the energy in Eq. ( 15) becomes E = E 0 G 2 Ns /2. We know that for γ = 0 and for positive ε, where the system is in the rigid regime, the minimum of the energy is strictly positive, E > 0. Using Eq. (18), this is only possible if the matrix MNsmn is positive semi-definite. Minimization of the energy E = E 0 G 2 Ns /2 implies that G Ns = 0 or else G Ns = 0 is minimal. As a consequence, using Eq. ( 18), none of the ∆r n with 1 ≤ n ≤ N dof can be of lower order than O(|ε| 1/2 ) + O(γ). Thus, because up to order O(ε) + O(γ 2 ) we have G p = 0 for p < N s , Eq. ( 16) implies that none of the ∆r p with p < N s is of lower order than O(ε) + O(γ 2 ). Further, this implies that any terms involving the ∆r p with p < N s can be neglected in the G p with p ≤ N s . Thus, each of the ∆r p with p < N s appears in total only once, namely as linear term in the respective G p . As a consequence, each of the equations G p = 0 for p < N s can be fulfilled by adjusting ∆r p accordingly. To discuss the general case of non-vanishing B(1) p , we first simplify by applying shifts to both ε and γ,

such that in the athermal limit (ε, γ) = (0, 0) is on the transition line and, after again carrying out the Taylor expansion, Eq. ( 13), the linear-order γ term in G Ns vanishes at the energy minimum:

Note that the G p with p < N s do in general still have a linear-order term in γ:

p γ2 , (B4) and that the ∆r n are also redefined due to the shifts in ε and γ. If ∆ε and ∆γ are sufficiently small such that to lowest order C, M , and B (2) in Eq. ( 9) are unaffected, it can be shown that the shifts are given by ∆ε = ( B(1) Ns ) 2 /4 w2 Ns b ε and ∆γ = -B(1) Ns /2 wNs b ε , where B (1) Ns is the linear-order term before the shift, and b ε is given by Eq. (B8) below.

We first discuss the athermal limit. Using the arguments from appendix A, one finds again that at the energy minimum, all G p = 0 for p < N s up to order O(ε) + O(γ 2 ), and none of the ∆r n with 1 ≤ n ≤ N dof can be of lower order than O(|ε| 1/2 ) + O(γ). As a consequence, up to linear order in γ, we have from Eq. (B4):