Stiffening of under-constrained spring networks under isotropic strain

Abstract

Disordered spring networks are a useful paradigm to examine macroscopic mechanical properties of amorphous materials. Here, we study the elastic behavior of under-constrained spring networks, i.e. networks with more degrees of freedom than springs. While such networks are usually floppy, they can be rigidified by applying external strain. Recently, an analytical formalism has been developed to predict the scaling behavior of the elastic network properties close to this rigidity transition. Here we numerically show that these predictions apply to many different classes of spring networks, including phantom triangular, Delaunay, Voronoi, and honeycomb networks. The analytical predictions further imply that the shear modulus G scales linearly with isotropic stress T close to the rigidity transition. However, this seems to be at odds with recent numerical studies suggesting an exponent between G and T that is smaller than one for some network classes. Using increased numerical precision and shear stabilization, we demonstrate here that close to the transition a linear scaling, G ∼ T, holds independent of the network class. Finally, we show that our results are not or only weakly affected by finite-size effects, depending on the network class.

---

Introduction

Understanding macroscopic rigidity and how it depends on the microscopic structure in amorphous materials such as fibrous networks, glasses, jammed colloids, and granular materials has been a long-standing challenge in the field. While the macroscopic mechanics of crystalline materials can be computed explicitly by exploiting their spatially periodic microscopic structure, this is not possible for disordered materials. In particular, upon deformation disordered materials generally display non-affine microscopic displacements, which are hard to predict.^1–11

A classical way to predict the onset of rigidity in many systems is to use Maxwell's constraint counting, which states that rigidity emerges whenever the constraints in a system outnumber its degrees of freedom.^12–14 In systems with pair interactions, this is equivalent to comparing the average connectivity z, i.e. the average number of pair interactions each particle is involved in, to the number of degrees of freedom per particle, which is given by the dimension of space, D. Such systems are predicted to be rigid if z exceeds the isostatic point, z > z_c := 2D. In this case the system is called over-constrained. Otherwise, for z < z_c, the system is called under-constrained or sub-isostatic, and is predicted to be floppy.

While Maxwell's constraint counting predicts under-constrained systems to be floppy, these systems can still be rigidified, either through the application of external strain or the presence of residual stresses.^{3,6,7,11,15–21} As a simple model to study such strain-induced rigidity, we discuss here strain-induced rigidity in athermal, under-constrained disordered spring networks.^{16,18,19,22–24} Strain-induced rigidification is illustrated in Fig. 1a for a network to which isotropic and shear strain has been applied.

Fig. 1 Under-constrained spring networks can be rigidified by the external application of strain. (a) Network configurations and spring tensions for an example network under either isotropic or shear strain. (b) Change of minimal network energy E, bulk modulus B, and shear modulus G when isotropically deforming a network across the rigidity transition, which occurs at the critical strain value ε*. Below ε*, all springs can attain their rest lengths, the system is floppy, and E = B = G = 0. Beyond ε*, springs start to deviate from their rest lengths. According to ref. 18, the bulk modulus B shows a discrete jump at ε*, while G increases linearly and E increases quadratically with the distance from the transition point ε*. (c) Two types of spring potentials are used in our simulations, harmonic (left) and rope-like (right). [small script l] and [small script l]₀ denote spring length and rest length, respectively.

The mechanism creating strain-induced rigidity has been discussed in the literature before.^15,16,25,26 When approaching the transition from the floppy side, a state of self-stress (SSS) forms right at the transition. A SSS is a set of tensions that could be put on the springs without any net forces on the nodes. The SSS that appears at the rigidity transition couples to isotropic strain, and using known approaches it can be shown that this induces a jump in the bulk modulus right at the transition (Fig. 1b).^14,18 Meanwhile, the shear modulus shows a continuous transition, whenever the SSS that appears at the transition has no net overlap with shear strain. Previously, the floppy side of the strain-stiffening transition was discussed in the limit where the springs are infinitely rigid.^25,26 Here, we are interested in the network mechanics of the rigid side of the transition when spring constants are finite.

Recent work involving one of us proposed a theoretical approach that allows to analytically predict the behavior of the elastic properties of under-constrained materials close to the strain-induced rigidity transition.¹⁸ This approach is based on a minimal-length function that formalizes the relationship between spring lengths and the applied global strain. This minimal-length function both reflects the critical point where the network starts to rigidify and allows to predict the behavior of the elastic network properties in the rigid regime. In ref. 18, this approach was numerically verified both on models for disordered cellular materials and for packing-derived spring networks. However, it has never been explicitly tested for other classes of under-constrained spring networks.

The approach in ref. 18 allows to predict the behavior of the elastic moduli close to the transition, where the bulk modulus B shows a discontinuity, while the shear modulus G increases linearly with isotropic strain ε (Fig. 1b). One can show that as a consequence of both, one would expect the shear modulus G to linearly increase also with isotropic stress T close to the transition. This is also consistent with earlier work on stress-induced rigidity.^3,15,27 However, more recent numerical work¹⁹ on under-constrained disordered spring networks suggested that the value for the scaling exponent between G and T can differ from one, depending on the class of network studied. The reason for this deviation from the analytical predictions is so far unclear. Other recent work proposed that the numerical results in ref. 18 could potentially be affected by finite-size effects caused by a diverging length scale when shearing the system.¹¹ Could similar finite-size effects be the reason for this contradiction between predicted and numerically obtained exponents between G and T?

Here, we numerically test the predictions from ref. 18 on several different classes of athermal spring networks. These include phantom triangular, Delaunay, Voronoi, and honeycomb networks, where we study two types of spring potentials, harmonic and rope-like (Fig. 1c). In the following, we first summarize the analytical approach from ref. 18 in Section 1. We then test the analytical predictions on the four different network classes in Section 2, and show that they follow the predicted behavior (Section 2.2). In Section 2.3, we furthermore show that the scaling behavior of the coefficients appearing in the minimal-length function with connectivity z depends on the network class. We then numerically explore the scaling behavior of the shear modulus G over isotropic stress T with increased numerical precision and find a scaling exponent of one, independent of network class (Section 2.4). Finally, we show that depending on the network class, there is no or a weak system-size dependence affecting these results (Section 2.5).

1 Theoretical predictions

We start by summarizing the approach of ref. 18, which allows to predict the scaling behavior of the elastic properties of under-constrained spring networks close to the rigidity transition.

In general, the formalism of ref. 18 applies to any disordered Hookean spring network of N springs, where each spring i has a different spring constant k_i and rest length [small script l]_0i. The energy of such a network is:

where [small script l]_i is the length of spring i. The springs are connected at movable nodes, around which they can freely rotate. While the approach can be applied largely independently of the precise boundary conditions, we focus here on periodic boundary conditions with fixed system size. Unless stated otherwise, we use dimensionless quantities, where the length unit is L_c := (V/N)^1/D with D being the dimension of space and V the system volume. We define the energy unit such that . Using dimensionless lengths will later allow us to describe the effect of isotropic strain (Section 1.3).

Here, to explain just the key ideas of the approach, we focus for simplicity on the special case of a network with homogeneous spring constants k_i = 1 and rest lengths [small script l]_0i = [small script l]₀:

The behavior of networks with heterogeneous spring properties can be predicted by formally mapping them onto eqn (2) as discussed in Appendix A.

The elastic properties of disordered networks are in general difficult to predict analytically. Formally, these elastic properties can be computed from derivatives of a minimal energy function e_min(ε,γ) with respect to external isotropic strain ε or shear strain γ. This function corresponds to the minimized system energy e({r_n},ε,γ) with respect to the node positions r_n at constant strain variables ε,γ. However, applying strain to a disordered network generally induces non-affine displacements of the node positions, which are typically hard to predict without numerical energy minimization. To nevertheless make non-trivial predictions about the scaling of the elastic properties, ref. 18 introduced a different approach. Instead of explicitly following the node motion, progress can already be made by focusing on the relation between spring lengths and external strain.

Note that while we focus in this section on harmonic springs, the formalism can also be applied to networks with rope-like pair interactions (Fig. 1c). This is because a rope-like pair interaction can be perfectly mimicked by a chain of two or more harmonic springs¹⁸ (see also Appendix D).

1.1 Key idea

To obtain an explicit expression for e_min in terms of external strain, we first transform the expression in eqn (2) into a sum of two squares:Here, and are average and variance of the spring lengths, respectively.

The expression in eqn (3) allows us to more conveniently discuss the minimal network energy e_min and its behavior once we strain the system. Because e is the sum of two squares, an energy minimum is attained whenever both and σ_{[small script l]} are as small as possible. There are two possibilities. First, if there is a set of node positions such that both squares can simultaneously attain zero, then the minimal energy is zero e_min = 0. Because elastic stresses and moduli correspond to derivatives of e_min, the system is floppy in this parameter regime. Second, there might be no set of node positions such that both terms and σ_{[small script l]} can simultaneously vanish. In this regime, the system is typically rigid.

To access the value of e_min in the rigid regime, we need to understand how the system compromises between minimizing and σ_{[small script l]} in eqn (3). To this end, we first need a way to express which combinations of and σ_{[small script l]} are geometrically possible. As shown in ref. 18, this can be done using a minimal-length function which returns the minimally possible for a given σ_{[small script l]}. In other words, a combination of and σ_{[small script l]} is geometrically possible only if:

For instance, for σ_{[small script l]} = 0 it is possible to find only network configurations with . Thus, for the network will be floppy, because both squares in eqn (3) can simultaneously vanish, which implies that e_min and its derivatives vanish. Conversely, for the first term in eqn (3) can not vanish with σ_{[small script l]} = 0, because only configurations with are possible. Thus, is the transition point between floppy and rigid regime.

In general, the precise functional form of depends on the network structure. However, we showed in ref. 18 that at the transition point to first order in σ_{[small script l]} this function can be expanded as

where and a_{[small script l]} are constants that encode the network structure. Eqn (5) holds in the limit of small σ_{[small script l]}, which means that the system is close to the transition point, where σ_{[small script l]} = 0. We expect that deriving expressions for and a_{[small script l]} from first principles is very hard for disordered networks. Besides some exceptions, and a_{[small script l]} will need to be determined numerically.

We briefly outline how eqn (5) is a consequence of the SSS that is created at the transition point (details in ref. 18). At the transition point the N-dimensional vector of all spring lengths [small script l]_i is given by with e = (1,…,1). We consider a state of the network that is slightly shifted into the rigid regime, where . One can then show that the definition of an SSS implies that is perpendicular to it:¹⁸

Here, t is the SSS created at the transition. To find the minimal possible , we first decompose both and t into parts parallel and perpendicular to and t ∼ e + a_{[small script l]}m_t, where we use the normalization m_{[small script l]}² = m_t² = N. Note that this defines a_{[small script l]} as the coefficient of variation of the SSS components. We then insert both decompositions into eqn (6) and obtain:Thus, to minimize for fixed σ_{[small script l]}, the scalar product m_{[small script l]}·m_t needs to be maximized. In the presence of only a single SSS, m_{[small script l]} can be any normalized vector perpendicular to e, and thus m_{[small script l]}·m_t is maximized by m_{[small script l]} = m_t, for which eqn (7) becomes eqn (5).

To derive an expression for the minimal energy e_min in the solid regime, we combine two parts: the energy in eqn (3) and the condition of geometrically possible combinations in eqn (4) and (5). First, eqn (4) implies that for fixed σ_{[small script l]}, the energy in eqn (3) is minimized when . Combining this with eqn (5), insertion into eqn (3), and minimization with respect to σ_{[small script l]}, yields:

This expression only depends on the spring number N, the rest length [small script l]₀, and the two parameters and a_{[small script l]} that encode the network structure. Note that from eqn (8) we see that the system energy is that of a single effective spring with rest length [small script l]₀.

1.2 Simple example network

To illustrate the ideas of the previous section, we discuss a simple example network (Fig. 2a left). The network consists of four springs with equal dimensionless spring constants k = 1 and rest lengths [small script l]₀. Two of the springs are connected to fixed points (black dots) located at positions (−1/2,0) and (1/2,0), respectively. The two internal nodes (red dots) at positions r_n with n = 1, 2 are movable. We will use the ideas of the previous section to derive an expression for the minimal energy e_min.

Fig. 2 Illustration of the analytical formalism using a 4-spring example network. (a) The network can be rigidified by either decreasing the dimensional spring rest length or increasing the system size, both of which has the effect of decreasing the dimensionless parameter [small script l]₀. Black and red dots indicate fixed and movable nodes, respectively. (b) Dependence of average and standard deviation σ_{[small script l]} of the four spring lengths on the internal node positions r_n = (x_n,y_n) with i ∈ {1,2} and y₁ = y₂ = 0. The axes are x̃₁ = x₁ + 1/6 and x̃₂ = x₂ − 1/6. Curves of constant are diagonal lines (with increasing : blue solid, black dashed, red dotted lines), while curves of constant σ_{[small script l]} are ellipses. The configuration of minimal for given σ_{[small script l]} is indicated by the red dot. Because the linear size of the ellipse scales with σ_{[small script l]}, the minimal for given σ_{[small script l]} decreases linearly with σ_{[small script l]}.

For [small script l]₀ ≥ 1/3, there are always configurations where all springs can attain their rest lengths [small script l]_i = [small script l]₀ (Fig. 2a top). This implies that and σ_{[small script l]} = 0, i.e. both terms in eqn (3) can simultaneously vanish.

Conversely, for [small script l]₀ < 1/3, the springs will be under tension (Fig. 2a bottom). Our 4-spring example network is simple enough so that we can explicitly minimize the energy with respect to the inner node positions r_n = (x_n,y_n) with n = 1, 2. This will allow us to first directly test whether the minimal energy has the form predicted by eqn (8) in the previous section. The energy of our example network is

Here, to simplify the following discussion, we have set y₁ = y₂ = 0. The energy e has a global minimum at x₁ = −(1 + 2[small script l]₀)/10, x₂ = (1 + 2[small script l]₀)/10, where its value is

This expression is indeed of the predicted form eqn (8).

We now demonstrate how the minimal energy e_min can instead be obtained using the ideas of the previous section. We first discuss which pairs of and σ_{[small script l]} are geometrically possible. To this end, we express and σ_{[small script l]} in terms of the internal degrees of freedom x₁,x₂:

where we defined x̃₁ = x₁ + 1/6 and x̃₂ = x₂ − 1/6. Both eqn (11) and (12) are illustrated in Fig. 2b. Curves of constant correspond to lines inclined by 45, where increases as x̃₁ decreases and x̃₂ increases. Meanwhile, curves of constant σ_{[small script l]} correspond to ellipses centered at x̃₁ = x̃₂ = 0, whose main axes scale with σ_{[small script l]} and are oriented at 45 angles with respect to the x̃₁ and x̃₂ axes.† Thus, for a given value of σ_{[small script l]}, any combination of x₁ and x₂ can give rise to values for only in an interval between (blue solid line) and (red dashed line). The upper bound only exists because we set y₁ = y₂ = 0 before; without this constraint, can become arbitrarily large for a given σ_{[small script l]}.‡ Meanwhile, the lower bound decreases linearly with the distance between origin (black dot) and the intersection point (red dot) in Fig. 2b, which is proportional to σ_{[small script l]}. As a consequence, using eqn (11) and (12):This is of the form of eqn (5), where we identify and a_{[small script l]} = 1/3. Inserting this into eqn (8), we obtain indeed eqn (10).

In our discussion here we included the internal degrees of freedom x₁, x₂ to demonstrate their connection to geometrically possible combinations of and σ_{[small script l]}, and to obtain explicit values for and a_{[small script l]}. In general, however, the approach from ref. 18 does not require a discussion of internal degrees of freedom. Eqn (3)–(5) are sufficient to understand the overall system behavior close to the rigidity transition, unless one wants to derive the values of the coefficients and a_{[small script l]} from first principles.

1.3 Effect of isotropic strain

We now discuss how the effect of isotropic strain ε is incorporated into the formalism. The 4-spring system in Fig. 2a transitions from floppy to rigid when decreasing the dimensionless parameter [small script l]₀. Such a decrease in [small script l]₀ can correspond either to a decrease in the dimensional spring rest length while keeping the system size constant (Fig. 2a left), or to an increase in system size while keeping the dimensional rest length constant (Fig. 2a right). Thus, [small script l]₀ is a control parameter combining both dimensional spring rest length and isotropic strain.

Let us consider simulations where the dimensional spring rest length L₀ is kept constant, but the system size V is changing. In this case, the combined control parameter [small script l]₀ encodes isotropic strain. We define (linear) isotropic strain as ε := (V/V_ref)^1/D − 1, where the V_ref is the system volume right after creation of the network. From our length non-dimensionalization follows that we can convert between [small script l]₀ and bulk strain ε via:

Inserting this equation into eqn (8) provides an explicit expression of the dimensionless system energy on isotropic strain ε.

1.4 Effect of shear strain

To understand how shear strain enters the formalism, we first note that shearing the system does not change the energy formula eqn (3). However, shearing the system will change the set of geometrically possible combinations . Thus, shear strain needs to be included as a parameter in the minimal-length function . In ref. 18 this function is Taylor expanded to second order in shear strain, so that eqn (5) becomes most generally:For later compactness of notation, here we also substituted the notation of parameter by .

Note that the linear order term in [small gamma, Greek, circumflex] appears only because disordered systems with a finite size generally display a small but finite anisotropy. Eqn (16) can be simplified by removing this anisotropy through defining a new shear variable γ = [small gamma, Greek, circumflex] − Δγ₀, where [small gamma, Greek, circumflex] = Δγ₀ is defined as the shear where the function is minimal: Δγ₀ = − b₁/2b. Defining this leads to the minimal-length function:

The anisotropy Δγ₀ is expected to disappear in the limit of a large network size.

1.5 Elastic properties near the rigidity transition

Substituting eqn (17) into eqn (3) and minimizing with respect to σ_{[small script l]}, we obtain the following explicit energy expression in terms of the control parameters [small script l]₀ and γ:Derivatives of this expression with respect to [small script l]₀ (which is related to isotropic strain ε viaeqn (15)) and shear strain γ allow to derive the following quantities, here for the 2D case:¹⁸Here, T, σ, ΔB and G are isotropic stress, shear stress, bulk modulus discontinuity, and shear modulus, respectively. These formulas hold close to the rigidity transition in the region where eqn (17) is accurate. As shown by eqn (19)–(22), the three parameters , a_{[small script l]} and b fully describe the macroscopic elastic behavior in this regime.

2 Numerical results

While in ref. 18 the analytical predictions in eqn (19)–(22) were numerically tested on packing-derived networks only, we test these predictions here on a set of additional network classes. These include phantom triangular and Delaunay networks (both with varying connectivity z), as well as honeycomb and Voronoi networks (which both have fixed connectivity z = 3). We probe the elastic properties of these networks under isotropic (i.e. bulk) strain.

2.1 Network generation and energy minimization

Networks of freely hinging nodes are created in a periodic box following existing protocols^19,28 (details in Appendix B). We probe the system by varying isotropic strain ε. Each time, we first use bisection to detect the transition point ε*, before we carry out exponential and/or linear sweeps in isotropic strain ε (details in Appendix C). To ensure high precision in our energy minimization, we use an optimized conjugate gradient scheme that allows to reduce the average residual force per degree of freedom to less than 10⁻¹².¹⁸

Right after creation, where [small gamma, Greek, circumflex] = 0, the disordered networks will generally display an anisotropy. To remove this anisotropy, we need to shear the system to the state [small gamma, Greek, circumflex] = Δγ₀ (i.e. γ = 0, see Section 1.4). At this point, according to eqn (20), shear stress vanishes, σ = 0. Thus, the anisotropy in the networks can be numerically removed using shear stabilization.²⁹ Shear stabilization means that shear strain is treated as an additional degree of freedom during the energy minimization. Unless stated otherwise, we always apply this method during the bisection phase to search for the transition point, so that our system right after the bisection phase is at (ε,γ) = (ε*,0). During the subsequent ε sweeps, we keep shear strain γ fixed (details in Appendix C).

2.2 Elastic moduli close to the transition

To numerically characterize the nature of the transition, we first carry out a combination of exponential and linear sweeps around the transition point ε* (details in Appendix C). In Fig. 3, we plot bulk modulus B and shear modulus G against isotropic strain ε for single network realizations with varying connectivity z, where we use harmonic spring potentials.

Fig. 3 Behavior of bulk and shear moduli across the transition for four different classes of networks and comparison to analytical predictions. (a)–(c) Sketches of network structures of phantom triangular (a, z = 3.2, W = 10), Delaunay (b, z = 3.2, W = 12), Voronoi (c, z = 3, W = 12) and honeycomb (c, z = 3, W = 12) networks. In (a), the solid lines are springs of length one, while the dashed ones are phantom springs, i.e. springs that cross one or more nodes without being connected to them. We show phantom springs, which are actually straight, slightly curved only for better visualization. In (b), gray lines indicate the removed springs from the initial full Delaunay network, leaving the black springs in the actual network. (d)–(i) Numerically obtained bulk modulus B and shear modulus G at γ = 0 against increasing isotropic strain ε for different network classes: phantom triangular (d, g, W = 40) and Delaunay (e, h, W = 20) networks with variable connectivity z, as well as Voronoi (f, i, W = 70) and honeycomb (f, i, W = 60). Here, we use harmonic spring potentials, and we shear stabilized the networks before the ε sweeps (Appendix C). The discontinuity ΔB in the bulk modulus at the transition point and linear scaling of G predicted from eqn (21) and (22) are indicated as solid bars and solid lines, respectively. (h inset) Log-log plot of G against strain difference ε − ε*, to resolve the vicinity of the transition point ε*. Symbols are numerical data and lines are analytical predictions.

At the transition, all networks show a discontinuity ΔB in the bulk modulus, while the transition is continuous in the shear modulus G. This is qualitatively consistent with our analytical predictions (Section 1.5) and the behavior of packing-derived networks.¹⁸ We also observe that for both phantom triangular and Delaunay networks the transition point ε* decreases with the average connectivity z.

To compare these data to the prediction for the bulk modulus discontinuity ΔB according to eqn (21), we need the values of and a_{[small script l]} for our simulations. To extract , we insert the transition point strain value ε* into eqn (15). To extract a_{[small script l]}, we plot σ_{[small script l]} over (inset of Fig. 4d) and perform a linear fit whose slope is a_{[small script l]} for small (see Appendix D). Note that for symmetry reasons, the honeycomb lattice has σ_{[small script l]} = 0 and thus a_{[small script l]} = 0. The resulting predictions for the bulk modulus discontinuities ΔB are respectively indicated as horizontal bars in Fig. 3d–f. Indeed, our predictions match well the discontinuities present in the simulation data for all four network classes and all connectivities z.

Fig. 4 Scaling of the parameters , a_{[small script l]} and b with network connectivity z for phantom triangular (system size W = 40), Delaunay (W = 20), Voronoi (W = 70; z = 3) and honeycomb (W = 60; z = 3) networks with harmonic spring potentials. The parameters , a_{[small script l]} and b are extracted for each network realization according to the protocols in Section 2 and Appendix D. For instance, (d inset) a_{[small script l]} is extracted using a linear fit of σ_{[small script l]} over , and (g inset) b is extracted using a linear fit of G over . Error bars in all panels indicate the standard error of the mean. We use the first 5 data points away from Δz = 0 to fit and obtain 2.07–0.31Δz for phantom triangular and 1.30–0.30Δz for Delaunay networks (red dashed lines in panels (a) and (b)). For the honeycomb network, the numerically obtained value for is consistent with its theoretical value , and the parameter a_{[small script l]} is exactly zero, because σ_{[small script l]} = 0 due to symmetry.

Some of the data points right at the transition fall below the analytically predicted value for ΔB. These deviations occur in our data for the strain values closest to the transition point, while ΔB values of the same network at similar strain values match closely with our analytical prediction. These deviations are likely due to numerical residues, an effect that we observed before.¹⁸

For the shear modulus, eqn (22) predicts a continuous transition with a linear scaling . Using eqn (15), this implies also a linear scaling to lowest order in Δε. Indeed, this is what we observed close to the transition (inset of Fig. 4g). We indicate this linear scaling also in Fig. 3g–i. Note that for larger Δε, non-linearities in and in eqn (15) create deviations from this prediction.

Note that both honeycomb and Voronoi networks have their transition points at ε* = 0. This means that these networks have a SSS already right at creation. While this is clearly the case for the honeycomb lattice, we show in Appendix E that it is also true for any Voronoi network.

Taken together, the elastic properties of the system close to the transition, such as the transition point ε*, the magnitude of the discontinuity ΔB in the bulk modulus, and the linear scaling coefficient for the shear modulus G, can be predicted from the coefficients , a_{[small script l]} and b.

2.3 Scaling of , a_{[small script l]} and b with connectivity z

In Fig. 4, we plot the parameters , a_{[small script l]} and b for phantom triangular, Delaunay, Voronoi, and honeycomb networks with harmonic spring potentials. For phantom triangular and Delaunay networks, we show the dependency on the connectivity z. For the disordered networks (i.e. phantom triangular, Delaunay, and Voronoi) we average each time over 50 random realizations.

In both phantom triangular and Delaunay networks, close to the isostatic point the parameter exhibits a linear dependence on Δz with a negative coefficient (Fig. 4a and b), which has also been observed in 2D packing-derived networks.¹⁸

We also examined the z-scaling exponents of a_{[small script l]} and b close to isostaticity (Fig. 4d, e, g and h). For a_{[small script l]}, we find for phantom triangular networks a scaling exponent of ≈0.2, for Delaunay networks, we find an exponent close to −1, while we found an exponent of −0.5 in earlier work for packing-derived networks. Meanwhile for b, we find for phantom triangular networks an exponent of ≈−0.5 or smaller, for Delaunay networks an exponent of roughly −2, while we have found before for packing-derived networks an exponent of −1. Hence, the scaling exponent of both parameters strongly depends on network class.

Note that for a small fraction of the Delaunay networks, we did not observe a linear scaling between σ_{[small script l]} and , suggesting that the linear relation between and σ_{[small script l]} might be violated for these networks (Appendix D). A more detailed examination suggests that this could quite possibly be due to finite numerical cutoffs required to identify the transition point, which would make us miss the regime where this scaling is linear (Appendix F). We excluded these networks from the averages shown in Fig. 4. We stress that we only found deviations from the linear scaling for Delaunay networks with harmonic springs, while we could numerically confirm the predicted linear scaling for all phantom triangular and Voronoi networks, as well as the honeycomb network.

2.4 The shear modulus scales linearly with isotropic stress.

In the previous Sections (1.5 and 2.2), we showed that the shear modulus G scales linearly with the isotropic strain beyond the transition point, Δε = ε −ε*. Moreover, a finite bulk modulus discontinuity at ε* implies that the isotropic stress T also scales linearly with Δε to lowest order, which can be derived form eqn (15) and (19). Hence, we would expect from the analytical predictions in Section 1 that the shear modulus scales linearly with the isotropic stress:

G ∼ T^α with α = 1. (23)

However, recent numerical work has suggested different values for α. For instance, ref. 19 studied networks with rope-like potentials, and for z = 3.2 the results suggested an exponent of α ≈ 0.85 for phantom triangular and α ≈ 0.9 for Delaunay networks, while α = 1 was found for honeycomb and Voronoi networks.

To resolve this contradiction between the numerical results from ref. 19 and our analytical results from Section 1 and ref. 18, we simulate here different kinds of rope-like networks with a high numerical precision, where we vary linear system size W by more than an order of magnitude. Fig. 5a shows the scaling of the shear modulus G against the isotropic stress T, both averaged over 50 realizations, for phantom triangular networks, where we used two protocols. The open symbols correspond to a protocol without any shear stabilization. This means that no shear strain was applied after the creation of the network, and the ε sweep was carried out at [small gamma, Greek, circumflex] = 0. The closed symbols correspond to a protocol where we used shear stabilization when searching for the transition point, and as a consequence the ε sweep was carried out at γ = 0 (see Sections 1.4).

Fig. 5 The shear modulus G scales linearly with isotropic stress T for shear-stabilized networks. Without shear stabilization an additional plateau appears at small T. (a) Dependence of G on isotropic stress T for phantom triangular networks with z = 3.2 and rope-like spring potentials for varying system size W, where we either use shear stabilization (closed circles) or not (open circles). All shear-stabilized networks show a linear scaling close to the transition point (at small T). The gray vertical line indicates the lowest T value probed in previous work.¹⁹ (a inset) The linear scaling in G(T) also appears for shear-stabilized honeycomb (W = 60), Voronoi (W = 70) and Delaunay (W = 20) networks. Error bars in panel a & inset indicate the standard error of the mean. (b) The G(T) curves for individual phantom triangular networks with z = 3.2 and W = 40 without shear stabilization also exhibit a plateau, confirming that the plateaus in panel a are not due to an averaging effect. (b inset) The variance of the network anisotropy Δγ₀ (defined below eqn (16)) across different randomly generated networks for a given system size W scales inversely proportional with W².

We find that indeed, for the protocol with shear stabilization (closed symbols), the shear modulus G scales linearly with isotropic stress T over many orders of magnitude (Fig. 5a for phantom triangular networks & inset for the other network classes). This observation is independent of system size. However, without shear stabilization (open symbols), at small stress T we observe a plateau, whose value depends on system size. Away from the plateau regions the curves largely collapse for different system sizes.

The appearance of a plateau in G(T) in the protocol without shear stabilization can be readily understood from our analytical results. eqn (22) states that G is proportional to , where and γ = [small gamma, Greek, circumflex] − Δγ₀. Without shear stabilization, [small gamma, Greek, circumflex] = 0 and so γ = −Δγ₀. This implies a plateau in G that is proportional to Δγ₀². In other words, the plateau in G(T) is related to the small anisotropy in the disordered networks. Shear stabilization removes this anisotropy and thus also the plateau in G(T).

To test whether the plateaus that we find in Fig. 5a do not result form an averaging effect, we plot in Fig. 5b the same curves for individual realizations for a given system size. We find that the plateau is also present in individual simulations, and that its height fluctuates across realizations. This makes sense, because the network anisotropy Δγ₀ also fluctuates across realizations. Moreover, we find that the variance of Δγ₀ decreases inversely proportional to the number of springs in the system (Fig. 5b), which scales as ∼W². Hence, the plateau in G(T) corresponds to a finite-size effect. A similar conclusion was drawn also in ref. 30 following a different line of argument.

2.5 Linear range of shows no or weak system-size dependence

In recent work, it was pointed out that scaling exponents in spring networks under shear strain may be affected by finite-size effects.¹¹ In particular, it was suggested that for networks of size W, finite-size effects could affect scaling exponents when shearing the systems by Δγ ≲ W^−1/ν beyond the transition point γ*, where ν > 0. This would correspond to a diverging length scale ξ ∼ |Δγ|^−ν. While in this article, we probe the system mechanics with respect to isotropic strain ε instead, we still wanted to check whether such finite-size effects could affect our results.

The system mechanics with respect to isotropic strain is crucially determined by how the minimal-length function scales with σ_{[small script l]} in eqn (15) (Section 1, Appendix H). Hence, we were wondering whether the linear scaling of with σ_{[small script l]} is only valid close to the transition point with strains Δε < Δε_max ∼ W^{−1/ν_{[small script l]}} for some ν_{[small script l]} > 0. In other words, we wondered whether the range Δε_max of linear scaling would algebraically decrease to zero with increasing system size W.

In Fig. 6 we show the resulting dependency of Δε_max on system size W for both phantom triangular networks and packing-derived networks, both with rope-like spring potentials, where Δε_max is quantified as described in Appendix G. For the phantom networks, beyond an initial quick decrease in Δε_max for small W, we find that for W ≥ 60 our data indicates a finite-size scaling exponent in the range 1/ν_{[small script l]} ∈ [0.04,0.35]. Thus, the range of linear scaling in the function slowly decreases with system size. In contrast, for the packing-derived networks discussed in ref. 18, we find a range of 1/ν_{[small script l]} ∈ [−0.03,0.08]. This means that the linear scaling range of is subject to none, or at most a weak finite-size effect. Thus, intriguingly, the effect of finite-size effects on the linear scaling range of appears to depend on the class of network studied.

Fig. 6 Plot of the isotropic strain range Δε_max = ε_max − ε* within which the function scales linearly with σ_{[small script l]}, shown here for increasing system size W for phantom triangular networks, and for shear-stabilized packing-derived networks (inset). In both cases, we use z = 3.2 and rope-like potentials. The range Δε_max is quantified as explained in Appendix G and Fig. 9. While for W ≳ 60, the linear range of the phantom triangular networks appear to show a weak power-law dependence on W, there is no significant dependence on system size for packing-derived networks. Error bars indicate the standard error of the mean.

Given that strain-controlled transitions in spring networks have been shown to be critical transitions,^{4,6,7,10,11,19} we wondered whether we would also observe a divergence in the fluctuations close to the transition. Focusing on the scaling behavior of a non-affine motion parameter Γ with system size W and distance to the transition point, ε − ε*, we find indeed such a divergence (Appendix G). Moreover, we also find finite-size effects with exponents of 1/ν_Γ ≈ 0.75 for phantom triangular networks and 1/ν_Γ ≈ 0.6 for packing-derived networks. Intriguingly, these exponents are very different from what we observe for the linear scaling regime of , i.e. ν_Γ ≠ ν_{[small script l]}. This suggests that the linear range of is not controlled by the diverging length scale that controls non-affine motion. A possible reason for this is that close to the transition the non-affinity parameter Γ mostly captures motions that are (to first-order) unconstrained by spring lengths, while characterizes spring length behavior.

3 Discussion

We studied the elastic behavior of sub-isostatic spring networks that are rigidified by isotropic expansion, comparing numerical simulation results with analytical predictions from ref. 18. We first summarized the approach from ref. 18, which proposed an analytical framework to predict the behavior of the elastic network properties using a minimal-length function (eqn (17)). This minimal-length function allows to map the physical problem of the strain-induced stiffening transition to the purely geometric problem of finding a minimal length. Because of this reduction to a geometric problem, we expect these results to hold quite generally for any athermal, under-constrained material with a fixed connectivity.

The formalism allows to make several predictions of the elastic network behavior close to the transition.¹⁸ These predictions include the rigid-floppy boundary with respect to shear and isotropic strain, the value of the bulk modulus discontinuity at the transition, the linear scaling coefficient of shear modulus with isotropic tension, the value of the shear modulus discontinuity for networks under shear strain, the coefficient of the linear shear modulus scaling beyond this transition, and the coefficient describing the anomalous Poynting effect. All of these predictions are based on the three parameters , a_{[small script l]} and b. In other words, ref. 18 reduced the question of quantitatively predicting the elastic properties of disordered spring networks close to the transition to the question of determining these three parameters a priori, which we believe is a very hard problem that remains for future work. Currently, these parameters need to be determined numerically. However, by combining our predictions one also obtains non-trivial parameter-free predictions that apply to any athermal under-constrained material.¹⁸

A key result required to derive all our predictions is that close to the transition the minimal average spring length depends linearly on the standard deviation of the spring lengths σ_{[small script l]}. We have shown that this scaling results from the appearance of a state of self-stress (SSS) at the onset of rigidity.¹⁴ We note that an equivalent linear scaling has been reported before, in a system of spheres close to the jamming point.³¹ However, in such a system each contact change will also change the set of SSSs. Contact changes occur also in a host of other systems, such as vitrimers.³² Hence, investigating in depth how varying connectivity affects the relation between and σ_{[small script l]} is an interesting question for future research.

To numerically test the predictions, we first verified the linear scaling of the minimal-length function near the transition and extracted the three parameters , a_{[small script l]} and b for four different network classes, including phantom triangular, Delaunay, honeycomb, and Voronoi networks. Based on these parameters, we compute the bulk modulus discontinuity ΔB, which predicts well our numerical results for all network classes (Fig. 3). Moreover, we also recovered the predicted linear scaling of the shear modulus G with isotropic tension T close to the transition.

Next we explored the scaling of the parameters , a_{[small script l]} and b with respect to connectivity z. We found that the scaling of the parameters a_{[small script l]} and b with the distance to isostaticity Δz = 4 − z strongly depends on the network class. The scaling exponent for a_{[small script l]} can even change sign, varying from −1 for Delaunay networks to ≈0.2 for phantom triangular networks (Fig. 4, with an exponent of −0.5 for packing-derived networks¹⁸). The scaling exponent for b varies from −2 for Delaunay networks to −0.5 for phantom triangular networks (Fig. 4, with an exponent of −1 for packing-derived networks¹⁸). This dependency on network class is not too surprising, since the parameters , a_{[small script l]} and b depend on the microscopic network structure, which varies with network class. In contrast, the value of always showed a linear dependency on Δz, where intercept and slope depend on network class.

One prediction of the formalism in ref. 18 is a linear scaling of the shear modulus G with the isotropic stress T close to the transition point: G ∼ T^α with α = 1. This is a direct consequence of the discontinuity ΔB in the bulk modulus and of the linear scaling of the shear modulus G with strain Δε. There are different ways to derive such a linear scaling (e.g. ref. 15 and 27); ref. 18 showed that it can be understood as a consequence of the linear scaling of the function with σ_{[small script l]}. While in Section 1.4 we make the assumption that at its minimum is analytical in γ, even certain non-analytic behavior in γ would still lead to an integer exponent α (Appendix H). However, the prediction of integer α seems to be at odds with recent numerical work, which suggested that the value of α can be different from one for networks with a rope-like interaction potential, depending on the disordered nature of the network.¹⁹ In particular, ref. 19 found an exponent of α ≈ 0.85 for phantom triangular networks and α ≈ 0.9 for Delaunay networks with connectivity of z = 3.2.

To reconcile the two results from ref. 18 and 19, we numerically studied the G(T) scaling with an increased numerical precision, and our results confirmed the analytically predicted scaling exponent of α = 1 in both phantom triangular and Delaunay networks (Fig. 5a and inset). We show that the result also depends on a small random anisotropy in the generated network. In the presence of such a finite anisotropy, we observed a plateau in the shear modulus G(T) for small isotropic stress T, consistent with the analytic prediction, eqn (22). This plateau disappears when using shear stabilization,²⁹ which removes the network anisotropy by shearing the network by a shear strain Δγ₀ (Section 1.4). We moreover show that the plateau disappears for larger system sizes (Fig. 5b inset). Hence, while without shear stabilization large system sizes are required to probe the behavior close to the transition, shear stabilization allows to explore this regime already for smaller systems.

We see two possible reasons for the discrepancy in the G(T) scaling between ref. 18 and 19. First, we used the conjugate gradient minimizer code developed in ref. 18, which allows us to probe the system at least two orders of magnitude closer to the transition point than ref. 19 (see gray vertical line in Fig. 5a). For instance in phantom triangular networks we observe an exponent of α < 1 for larger isotropic stress T ≳ 10⁻³, which seems consistent with the value of 0.85 given by ref. 19, while we observe an exponent of α = 1 for stress T smaller than that. Second, we show that a small anisotropy in the generated network can lead to a plateau in the shear modulus curve G(T), which could in turn affect the inferred scaling exponent.

Previous work suggested that finite-size effects could affect scaling exponents in spring networks.¹¹ This can occur whenever the system size is on the order of or smaller than a length scale that diverges close to the transition point. While ref. 11 focused on shear simulations, we wanted to test whether such an effect could also arise in our isotropic-strain simulations. To this end, we numerically tested in what range around the transition point the linear scaling of the function holds. Our results suggest that this potentially depends on the class of network studied. While in phantom triangular networks, this range decreases weakly with system size, we did not find a significant decrease in packing-derived networks. This is also consistent with Fig. 5a, which suggests that G(T) is largely independent of system size W for the range of W probed.

Intriguingly, we found that non-affine motion Γ shows a much stronger system-size dependence, suggesting that it is controlled by a length scale that does not affect the linear range of . One possible reason for this is that close to the transition point non-affine motions are to linear order unconstrained in under-constrained networks. Better understanding this difference in the finite-size scaling behavior of Γ and is an interesting avenue for future research.

Conflicts of interest

The authors declare no conflict of interest.

Appendix

A Generalization to networks with heterogeneous spring constants and rest lengths

In the main text, we focused on the case where all springs share the same rest length [small script l]₀ and the same spring constant k. Here, we generalize this to networks where rest length and spring constant may differ among the springs, as in eqn (1). Similar to ref. 18, we introduce re-scaled spring lengths , re-scaled spring constants k̃_i, and an average spring rest length [small script l]₀ in a way that allows us to rewrite eqn (1) in the form:For this to work, we need to define the re-scaled spring lengths asThis will accordingly give rise to a new re-scaled spring constantsFinally, we choose to define [small script l]₀ as the quadratic mean of the [small script l]_0i, weighted by the k_i:Using the re-scaling eqn (25)–(27), eqn (1) can be exactly re-expressed as eqn (24). This network energy can be transformed intoHere, and σ_{[small script l]} are defined as the average and standard deviation of the re-scaled spring length with the weighting factors k̃_i:The subsequent discussion in Sections 1.1–1.5 remains unchanged.

B Network generation

Networks were created using the following protocols.

Phantom triangular (Fig. 3a)²⁸

Following ref. 19, a 2D triangular lattice of spacing 1 is first constructed by depositing three sets of W parallel filaments each at angles of 0°, 60° and 120° with the x-axis, respectively. To reduce the connectivity from z = 6 to values observed in e.g. collagen networks³³ of 3 ≤ z < 4, we first detach at each node one filament, which is randomly chosen among the three crossing filaments. This creates a network of homogeneous connectivity z = 4. To avoid system-spanning filaments, one spring is removed at a random position on each filament, giving the average connectivity z = 4 − 6/W. To further reduce the connectivity to a defined value z, we implement an iterative procedure. At each iteration, we randomly remove only a few of the springs and then clear off all of the dangling springs and isolated islands. This is repeated until the desired connectivity z is reached.

Delaunay (Fig. 3b)

Delaunay networks are constructed from W² nodes that are placed at uncorrelated random positions in a square box of side W. The connectivity of initially z = 6 is decreased to the desired value z by employing the same protocol using random cuts as for the phantom triangular networks.

Honeycomb (Fig. 3c)

We construct a network of W²/3 regular hexagons with side length 1.

Voronoi (Fig. 3c)

Voronoi networks correspond to the Voronoi tessellation of W²/2 nodes at uncorrelated random positions in a square box of side W.

In all four network classes, we set the dimensionless spring rest lengths [small script l]_0i to the respective initial spring lengths before any deformation is applied, i.e. at (ε,[small gamma, Greek, circumflex]) = (0,0). We set the dimensionless spring constants as the inverse of the respective rest length at zero strain, k_i = 1/[small script l]_0i.

C Details of numerical strain sweeps and computation of the elastic moduli

In this paper, we exclusively carry out sweeps of isotropic strain ε. Before each sweep, we first identified the transition point ε*. To this end, we implemented a bisection scheme, which we optimized by linearly interpolating the transition point in each step. We defined networks as rigid whenever their isotropic stress T is above a cutoff value of 10⁻¹⁰ (two orders of magnitude above the tolerance for the residual force cutoff per degree of freedom, 10⁻¹²). We use isotropic stress T as a criterion for network rigidity, because it is much faster to compute than an elastic modulus.

In the bisection to identify the transition point, we also implemented the option to perform shear stabilization to remove network anisotropy (Sections 1.4 and 2.1). This is done by treating the shear strain [small gamma, Greek, circumflex] as an additional degree of freedom during each energy minimization of the bisection process. In any case, shear stabilization was always turned off (i.e. shear remains constant) after the transition point ε* has been identified.

We apply an exponential sweep of isotropic strain to probe the scaling behavior of network mechanics close to the transition point ε*. In particular, we probed strain values ε − ε* = 10^−10+0.2k, where the step index k ranged from 0 to 51 by default, with only two exceptions. First, in Fig. 3, we apply the same exponential sweep, yet with k ranging from 0 to 7 only, which is then followed by a linear sweep. Second, for the Voronoi networks of size W = 70 (Fig. 4 and Fig. 5a inset) we needed to increase the residual force cutoff for the energy minimization to 10⁻¹⁰, and so we also raised the cutoff in isotropic stress T to identify the transition point to 10⁻⁸. Accordingly, we changed the sweep to the values ε − ε* = 10^−8+0.2k with k ranging from 0 to 41.

We computed the elastic moduli using two different methods. For not too big networks, we diagonalized the Hessian of the system energy and used the resultant eigenvalues to compute elastic moduli.^{18,27,34–36} This approach produces a higher numerical precision and was suitable for typical system sizes W < 100 (Fig. 3 and 4). However, in Fig. 5 we studied networks with a large system size, and so we used a less time-intensive way of computing the shear modulus G. We computed G through a difference quotient of the shear stress over the shear strain: G(ε, γ = 0) = [σ(ε, Δ) − σ(ε, −Δ)]/2Δ, where we numerically tuned and found the optimized shear strain Δ = 5 × 10⁻⁵. We also noticed that for ε − ε* < 10⁻⁷ the shear modulus computed with this method could deviate significantly from the true value. We hence excluded these data points in Fig. 5 and the lowest isotropic stress there is accordingly T ≈ 10⁻⁷.

D Extraction of the parameters a_{[small script l]} and b of the minimal-length function

To extract a_{[small script l]} from numerical data, one could just directly use the function (eqn (5)). However, this approach depends on the correct identification of the transition point . While we can identify with a relatively high precision of ∼10⁻¹⁰, we could even remove the dependency on entirely when determining a_{[small script l]}. To this end, we note that in an energy-minimized state, the energy is also minimal with respect to variation of σ_{[small script l]}, i.e. de/dσ_{[small script l]} = 0. From eqn (3), and using the insight that in the rigid regime, the minimization condition reads:Using eqn (5), the derivative of the minimal-length function is . Taken together, we thus obtain the linear relation:Based on eqn (31), examining the relation between σ_{[small script l]} and (e.g.Fig. 4d inset) allows both to effectively verify the scaling of the minimal-length function (eqn (5)), and to extract the value of a_{[small script l]}. This approach does not involve the critical value which we obtain with a lower precision as compared to and σ_{[small script l]} (as precise as 10⁻¹²).

To extract the parameter b, we use the derived shear modulus formula (eqn (22)), instead of directly using the original minimal-length function (eqn (17)) since we do not shear the networks (i.e. γ = 0). As before, we intend not to use the critical value . Thus, to replace the term that appears in the shear modulus formula (eqn (22)), we insert eqn (31) back into the minimal-length function (eqn (5)) and obtain

Combining this equation with the shear modulus formula (eqn (22)) yieldsWe used this equation to extract b from the plots of G over (Fig. 4g inset).

By default we use the first 25 data points from an exponential sweep to numerically fit eqn (31) and (33) and extract the parameters a_{[small script l]} and b. Note that a_{[small script l]} in Fig. 6 is defined in the very same way, based on the first 25 data points of an exponential sweep. Meanwhile in Fig. 3 we use only the first 5 data points due to a decreased step number n in the exponential sweep of isotropic strain (Appendix C).

We note that rope-like spring potentials can be treated as well with the analytical framework in Section 1, which we took into account when computing and σ_{[small script l]} here. While one way to treat rope-like spring potentials was discussed in ref. 18, where each spring is subdivided into a series of shorter springs, we chose here an alternative approach. We used the fact that for rope-like spring potentials, a spring i only affects the mechanics when [small script l]_i > [small script l]_0i, while springs with [small script l]_i < [small script l]_0i do not contribute. Thus, to compute and σ_{[small script l]}, whenever for any spring i the distance of the two connected nodes is smaller than the rest length [small script l]_0i, we set the spring length to [small script l]_i = [small script l]_0i. This redefinition of [small script l]_i does not affect the computation of shear modulus G and tension T.

E Any Voronoi network at creation has a state of self stress

We numerically found that Voronoi networks have a critical isotropic strain very close to zero, ε* ≈ 0. Here we show that the critical strain is indeed exactly zero, by proving that there is a state of self stress right at creation of these networks. In other words, at creation (ε = 0) these networks can sustain finite tensions in a subset of springs, while force balance is maintained at the internal nodes.

The geometric structure of a Voronoi network allows for the following set of spring tensions t_i (with i being a spring index) to be a state of self stress:

t_i = η |C²_i− C¹_i|. (34)

Here, η is some constant factor, the vectors C¹_i and C²_i refer to the two Voronoi seed points that are closest to spring i (Fig. 7; i.e.C¹_i and C²_i are the two points that generated the line that defines spring i), and |·| denotes the length of a vector.

Fig. 7 Illustration of local force balance at any node n of a Voronoi network, demonstrating that Voronoi networks have a SSS at creation. Red dots are the internal nodes, while blue dots are the neighboring random seeds used for the Voronoi tessellation. Geometrically, node n is created as the circumcenter of the local triangle (gray dashed lines) formed by these seeds, and the three local springs i, j and k (black segments) are the vertical bisectors of the respective sides. The vectors C¹_i and C²_i refer to the two seed points at the side perpendicular to spring i (and similarly for springs j and k). The node n is force-balanced when the magnitude of the spring tensile forces (f_i, f_j and f_k, in black arrows) follow the proportionality relation, eqn (34). These forces will form a closed triangle that is similar to the local triangle by a rotation of 90°, thus giving zero net force, see also eqn (35) and (36).

To show that the set of spring tensions t_i form a state of self stress, we demonstrate that they satisfy local force balance at each node. To this end, we focus here on a node n that is connected to springs i, j, k as shown in Fig. 7. The force that spring i exerts on the node n is f_i = t_ie_i, where e_i is the unit tangent vector of spring i pointing away from node n. Furthermore, we have

Here, in the second line, we inserted the spring tensions eqn (34). In the third line, we used the fact that spring i is perpendicular to the segment connected by the two seed points C¹_i and C²_i, while the operator R(π/2) performs a counter-clockwise rotation by an angle of π/2. An analogous equation to eqn (35) holds also for the forces by springs j and k. As a consequence, the sum of these three forces is zero:

f_i + f_j + f_k = 0. (36)

In other words, force balance on node n holds. This proof is also illustrated at the bottom of Fig. 7; up to the factor of proportionality η, the three forces f_i, f_j, f_k correspond to the three triangle sides rotated by π/2, which is why they add up to zero. Hence, Voronoi networks at creation have a state of self stress given by eqn (34).

F Apparent non-linear scaling of in some Delaunay networks

For Delaunay networks with harmonic spring potentials, we observed that a fraction of the networks did not seem to follow the linear relation (31) between σ_{[small script l]} and (blue and red data points in Fig. 8a inset). This is also apparent from the absence of a plateau in (compare blue and red with black data points in Fig. 8a). From our arguments in Appendix D, it follows that this non-linearity also implies a non-linear scaling of the minimal-length function with σ_{[small script l]}, which would also affect the elastic network properties, eqn (19)–(22).

Fig. 8 (a) The quotient plotted versus for Delaunay networks with harmonic spring potentials, shown here for three networks. These networks show either a linear (black diamonds) or non-linear (blue dots and red squares) scaling between σ_{[small script l]} and close to the transition point (inset), where a linear scaling is reflected by a plateau in the quotient . (b) Histograms of the extent T_max of the plateau in the quotient for Delaunay networks with harmonic spring potentials of different connectivity z. We compute T_max as the isotropic stress value where the quotient starts to show a deviation of more than 10% from the one computed at the detected transition point.

We wondered whether this non-linear scaling between σ_{[small script l]} and was just due to finite numerical cutoffs, or whether it reflects the real scaling behavior infinitesimally close to the transition point. Numerical limitations arise because we cannot probe the networks arbitrarily close to the true transition point. Indeed, we used a cutoff value of T_cutoff = 10⁻¹⁰ for the isotropic stress T to numerically identify the transition point. In other words, at the detected transition point we are already in the rigid regime by some small extent beyond the true transition point. If the plateau in exists only close to the true transition point until some isotropic stress value T_max < T_cutoff, we will not detect it since we missed that regime. To test if this could be the case, we created histograms of the extent T_max of the plateau for different connectivity z (Fig. 8b). For a given network, we define T_max as the isotropic stress of the data point at which the quotient first deviates by more than 10% from the value of this quotient at the detected transition point. For networks where the plateau ends below T_cutoff, we would find with this approach T_max ≈ T_cutoff = 10⁻¹⁰. If there is a significant excess of networks where we numerically do not observe a plateau, this could be an indication that there is indeed no plateau.

Fig. 8b shows that for Delaunay networks with harmonic springs, T_max generally decreases with connectivity z, and that we only observe a peak around T_max ≈ T_cutoff occur mostly for the two largest values of z. Even in these cases, the peak is not very pronounced and may very well arise from the integral of the real T_max distribution from 0 to T_cutoff. In other words, these networks may possibly have a plateau which ends just too close to the transition point for us to detect it.

This is also consistent with the observation that most of these curves appear to collapse with the curves that do show a plateau beyond the end of the plateau (Fig. 8a). This suggests that the non-linear scaling regime just corresponds to a regime governed by higher-order terms. In future work, it will be interesting to see if these higher-order terms could also be predicted from first principles.

G Finite-size effects

In the main text, we examined the range of validity of the linear scaling of with σ_{[small script l]}. In Fig. 9 we show how we determined this range using a 10% cutoff on the ratio . We find in Section 2.5 that this range does not or only weakly change with system size.

Fig. 9 Δε_max in Fig. 6 is defined as the strain range where the ratio shows a deviation of less than 10% from its value a_{[small script l]} defined very close to the transition point (Appendix C).

Strain-stiffening of spring networks has been shown to be a critical transition when using shear strain as control parameter, which includes diverging fluctuations when approaching the transition point.^6,10,11,19 We wondered whether we would also observe diverging fluctuations when using isotropic strain as control parameter. Analogous to previous work,^6,10,11,19 we quantify fluctuations using a non-affinity parameter Γ, which we define as:

Here, δu^NA_n is the (dimensionful) non-affine displacement of node n during an isotropic expansion by strain δε. The factor L_c is the length unit defined below eqn (1). Because the affine transformation corresponds in our case to uniform isotropic inflation, eqn (37) can be simplified using dimensionless node positions r, so that in practise we compute Γ as:Here, ε_k is the strain step with index k within a sweep, and r_n,k is the corresponding dimensionless position of node n.

Using the same networks as in Fig. 6, we numerically studied Γ(ε) and its dependence on system size. For all system sizes we observed a plateau in Γ(ε) (Fig. 10a), for both phantom triangular and packing-derived networks.

Fig. 10 Finite-size scaling of the non-affinity parameter Γ defined in eqn (37). (a) A plateau regime is observed in Γ at strain values close to the transition point ε*. Here we show only three example networks with W = 140 of different class. (b and c) As system size W increases, the range Δε_pl of the plateau decreases as a power-law, while the plateau value Γ_pl increases as a power-law. The plateau range Δε_pl here is defined as the strain difference between the transition point ε* and the point where Γ shows a 10% deviation from its plateau value Γ_pl. Error bars indicate the standard error of the mean (≈50 networks for each system size). Here we studied the same set of networks as in Fig. 6.

We examined how both height Γ_pl and extent Δε_pl of the plateau depend on system size (Fig. 10b and c). We quantified the extent Δε_pl as the value of ε where Γ deviates by 10% from the plateau value, where we also performed linear interpolation between neighboring ε_k values. We found power law scaling with system size of both plateau value Γ_pl ∼ W^λ_Γ/ν_Γ and plateau extent Δε_pl ∼ W^−1/ν_Γ. For the phantom triangular networks we found the plateau height exponent λ_Γ/ν_Γ ≈ 1.1 and for plateau extent 1/ν_Γ ≈ 0.75. For the packing-derived networks, we found the plateau height exponent λ_Γ/ν_Γ ≈ 0.5 and for plateau extent 1/ν_Γ ≈ 0.6.

These findings are consistent with the finite-size scaling behavior of the non-affinity parameter Γ with respect to shear strain γ:^6,10,11,19 The non-affinity parameter generally diverges when approaching the transition as Γ ∼ |Δγ|^−λ_Γ with λ_Γ > 0, but a diverging length scale ξ_Γ ∼ |Δγ|^−ν_Γ with ν_Γ > 0 changes this behavior for system sizes W ≲ ξ_Γ. As a consequence, the non-affinity parameter Γ (γ) has a plateau whose height scales as Γ_pl ∼ W^λ_Γ/ν_Γ, and whose extent scales as Δγ_pl ∼ W^−1/ν_Γ. Here we demonstrated that this behavior also appears when using isotropic strain instead of shear strain as control parameter.

Noticeably, the 1/ν_Γ values are much larger than the power-law exponent 1/ν_{[small script l]} that characterizes the range of the linear scaling of (Fig. 6; for phantom triangular networks 1/ν_{[small script l]} ∈ [0.04,0.35], while for packing-derived ones 1/ν_{[small script l]} ∈ [−0.03,0.08]). This indicates that the linear scaling regime of is not controlled by the diverging length scale ξ_Γ that governs the apparent divergence of the non-affine motions.

H General form of the minimal length function

In eqn (5), we Taylor-expanded the minimal length function to the second order in shear strain γ, while treating the coefficient a_{[small script l]} as independent of γ. Here we discuss a more general form of that can include potentially non-analytic dependencies on γ:where the coefficient a_{[small script l]}(γ) is a function of γ. We also newly introduced the function g(γ), where we choose the convention g(0) = 0; any offset can be absorbed into . Note that eqn (39) reflects an arbitrary dependency of on γ, while we keep the linear dependency on σ_{[small script l]}. Both a_{[small script l]} and g do not need to be analytic, but they need to be twice differentiable at γ = 0, or at least arbitrarily close to γ = 0.

After minimizing the energy with respect to inner degrees of freedom and the standard deviation σ_{[small script l]}, the resultant energy e is:

Here we defined . Using G = (d²e/dγ²)|_γ=0/N we then obtain for the shear modulus:where Q₀, Q₁, and Q₂ are coefficients that depend only on a_{[small script l]}, g, and their derivatives with respect to γ.

Note that according to eqn (19), Δ[small script l] is proportional to isotropic tension: T ∼ Δ[small script l], which can be understood as a consequence of the bulk modulus discontinuity. Since the Qs in eqn (41) do not depend on Δ[small script l], one already observes from this equation that any scaling G ∼ T^α needs to have α ∈ {0,1,2}. In this sense, the integer scaling exponent between G and T is inherited from the linear scaling of the function with σ_{[small script l]} in eqn (39).

Which of the three values α ∈ {0,1,2} is attained depends on the coefficients Q₀ to Q₂ in eqn (41), which are:

Q₀ = (1 + a_{[small script l]}²)²g′² (42)

Here, for simplicity we used the superscripts ′ and ′′ for the first and second derivatives with respect to γ, respectively. From eqn (41)–(44) follows that the system has finite shear modulus only if at least one of a_{[small script l]} or g has a finite first or second derivative with respect to γ at γ = 0.

With respect to the scaling exponent α we can say that first, α = 0 only if g has a finite first derivative. This corresponds to the case where there is a discontinuity in the shear modulus at the transition point. In these cases, the SSS that appears at the transition must have finite overlap with the shear deformation, i.e. the network is asymmetric (Section 1.4, e.g. the non-shear stabilized simulations in Fig. 5). Second, α = 1 only if g′ = 0 (i.e. the network is symmetric) and g has a finite second derivative. This is the typical case that we observe for shear-stabilized networks. Third, α = 2 would appear if the network is symmetric (g′ = 0), the second derivative of g vanishes, and a_{[small script l]} has finite first or second derivative. This situation might appear at a bifurcation (where g′′ as bifurcation parameter crosses zero), possibly related to a structural transition in the network.

Eqn (42)–(44) include first and second derivatives of g and a_{[small script l]}, which may formally not exist at γ = 0. An important possibility is that at least one of these four derivatives might diverge as γ → 0. However, with eqn (41) this would imply that generally also the shear modulus G diverges as γ → 0, which to our knowledge has not yet been observed in spring networks.

Acknowledgements

We thank Martin Lenz 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. The Centre de Calcul Intensif dAix-Marseille is acknowledged for granting access to its high-performance computing resources.

References

1. W. G. Ellenbroek, Z. Zeravcic, W. Van Saarloos and M. Van Hecke, EPL, 2009, 87, 6.
2. J. L. Silverberg, A. R. Barrett, M. Das, P. B. Petersen, L. J. Bonassar and I. Cohen, Biophys. J., 2014, 107, 1721–1730.
3. A. J. Licup, S. Münster, A. Sharma, M. Sheinman, L. M. Jawerth, B. Fabry, D. A. Weitz and F. C. MacKintosh, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 9573–9578.
4. J. Feng, H. Levine, X. Mao and L. M. Sander, Soft Matter, 2016, 12, 1419–1424.
5. A. S. V. Oosten, M. Vahabi, A. J. Licup, A. Sharma, P. A. Galie, F. C. MacKintosh and P. A. Janmey, Sci. Rep., 2016, 6, 19270.
6. A. Sharma, A. J. Licup, R. Rens, M. Vahabi, K. A. Jansen, G. H. Koenderink and F. C. MacKintosh, Phys. Rev. E, 2016, 94, 042407.
7. A. Sharma, A. J. Licup, K. A. Jansen, R. Rens, M. Sheinman, G. H. Koenderink and F. C. Mackintosh, Nat. Phys., 2016, 12, 584–587.
8. A. J. Licup, A. Sharma and F. C. Mackintosh, Phys. Rev. E, 2016, 93, 012407.
9. K. A. Jansen, A. J. Licup, A. Sharma, R. Rens, F. C. MacKintosh and G. H. Koenderink, Biophys. J., 2018, 114, 2665–2678.
10. J. L. Shivers, S. Arzash, A. Sharma and F. C. MacKintosh, Phys. Rev. Lett., 2019, 122, 188003.
11. S. Arzash, J. L. Shivers and F. C. MacKintosh, Soft Matter, 2020, 169, 6784–6793.
12. J. C. Maxwell, Lond. Edinb. Dublin Philos. Mag. J. Sci., 1864, 27, 294–299.
13. C. Calladine, Int. J. Solids Struct., 1978, 14, 161–172.
14. T. C. Lubensky, C. L. Kane, X. Mao, A. Souslov and K. Sun, Rep. Prog. Phys., 2015, 78, 73901.
15. S. Alexander, Phys. Rep., 1998, 296, 65–236.
16. M. Wyart, H. Liang, A. Kabla and L. Mahadevan, Phys. Rev. Lett., 2008, 101, 1–4.
17. D. E. Ingber, N. Wang and D. Stamenović, Rep. Prog. Phys., 2014, 77, 046603.
18. M. Merkel, K. Baumgarten, B. P. Tighe and M. L. Manning, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 6560–6568.
19. S. Arzash, J. L. Shivers, A. J. Licup, A. Sharma and F. C. MacKintosh, Phys. Rev. E, 2019, 99, 042412.
20. B. Cui, G. Ruocco and A. Zaccone, Granul. Matter, 2019, 21, 69.
21. O. K. Damavandi, V. F. Hagh, C. D. Santangelo and M. L. Manning, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2022, 105, 025003.
22. P. R. Onck, T. Koeman, T. van Dillen and E. van der Giessen, Phys. Rev. Lett., 2005, 95, 178102.
23. M. Sheinman, C. P. Broedersz and F. C. MacKintosh, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 021801.
24. M. F. J. Vermeulen, A. Bose, C. Storm and W. G. Ellenbroek, Phys. Rev. E, 2017, 96, 053003.
25. G. Düring, E. Lerner and M. Wyart, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 1–12.
26. R. Rens, C. Villarroel, G. Düring and E. Lerner, Phys. Rev. E, 2018, 98, 062411.
27. M. Merkel and M. L. Manning, New J. Phys., 2018, 20, 022002.
28. C. P. Broedersz and F. C. MacKintosh, Soft Matter, 2011, 7, 3186–3191.
29. S. Dagois-Bohy, B. P. Tighe, J. Simon, S. Henkes and M. Van Hecke, Phys. Rev. Lett., 2012, 109, 1–5.
30. O. K. Damavandi, M. L. Manning and J. M. Schwarz, EPL, 2022, 138, 27001.
31. S. Kooij and E. Lerner, Phys. Rev. E, 2019, 100, 042609.
32. S. Ciarella, F. Sciortino and W. G. Ellenbroek, Phys. Rev. Lett., 2018, 121, 58003.
33. S. B. Lindström, D. A. Vader, A. Kulachenko and D. A. Weitz, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 051905.
34. K. Huang and M. Born, Proc. R. Soc. London, Ser. A, 1950, 203, 178–194.
35. M. Born and K. Huang, Am. J. Phys., 1955, 23, 474.
36. A. Lemaître and C. Maloney, J. Stat. Phys., 2006, 123, 415–453.

---

Footnotes

† This is because eqn (12) can be transformed into
where u = x̃₁ + x̃₂ and w = x̃₁ − x̃₂. This is the equation of an ellipse whose main axes are diagonally oriented and scale with σ_{[small script l]}.

‡ That can become arbitrarily large for given σ_{[small script l]} < 1/4 can be shown explicitly by considering a subset of configurations parameterized by two scalars w and h as r₁ = (−w/2,h) and r₂ = (−w/2,h). Then one can show that the choice
leads to the correct value for the standard deviation of the spring lengths σ_{[small script l]}. Moreover, one can show that for this choice, the relation
holds, and that with given by eqn (14). Finally, for fixed σ_{[small script l]}, the function increases monotonically with h without upper bound.

---