Energetic rigidity I: A unifying theory of mechanical stability
Ojan Khatib Damavandi, Varda F. Hagh, Christian D. Santangelo, M. Lisa Manning
EEnergetic Rigidity: a Unifying Theory of Mechanical Stability
Ojan Khatib Damavandi, Varda F. Hagh, Christian D. Santangelo, ∗ and M. Lisa Manning † Department of Physics, Syracuse University, Syracuse, New York 13244, USA James Franck Institute, University of Chicago, Chicago, IL 60637, USA
Rigidity regulates the integrity and function of many physical and biological systems. In thiswork, we propose that “energetic rigidity,” in which all non-trivial deformations raise the energyof a structure, is a more useful notion of rigidity in practice than two more commonly used rigidiytests: Maxwell-Calladine constraint counting (first-order rigidity) and second-order rigidity. Wefind that constraint counting robustly predicts energetic rigidity only when the system has no statesof self stress. When the system has states of self stress, we show that second-order rigidity canimply energetic rigidity in systems that are not considered rigid based on constraint counting, andis even more reliable than shear modulus. We also show that there may be systems for which neitherfirst nor second-order rigidity imply energetic rigidity. We apply our formalism to examples in twodimensions: random regular spring networks, vertex models, and jammed packings of soft disks.Spring networks and vertex models are both highly under-constrained and first-order constraintcounting does not predict their rigidity, but second-order rigidity does. In contrast, jammed packingsare over-constrained and thus first-order rigid, meaning that constraint counting is equivalent toenergetic rigidity as long as prestresses in the system are sufficiently small. The formalism ofenergetic rigidity unifies our understanding of mechanical stability and also suggests new avenuesfor material design.
I. INTRODUCTION
How do we know if a material or structure is rigid?If we are holding it in our hands, we might choose topush on it to determine whether any small, applied dis-placement generates a proportional restoring force. Ifso, we say it is rigid. A structure that does not pushback is floppy. In this paper, we call this intuitive defini-tion of rigidity “energetic rigidity” by virtue of the factthat small deformations increase the elastic energy of thestructure. In many situations of interest, it is impossi-ble or impractical to push on a structure to measure therestoring force. In designing new mechanical metamate-rials, for example, we would like to sort through possibledesigns quickly, without having to push on every vari-ation of a structure. In biological tissues such as thecartilage of joints or the bodies of developing organisms,it is often difficult to develop non-disruptive experimen-tal rheological tools at the required scale. Or we maywish to understand how some tissues can tune their me-chanical rigidity in order to adapt such functionality intonew bio-inspired materials.To that end, we would like a theory that can predictwhether a given structure is energetically rigid rapidlyand without the need for large-scale simulations or ex-periments. Perhaps surprisingly, no such theory exists ineven simple cases. Indeed, determining whether even aplanar network of springs is rigid is NP-hard [1]. This hasinspired the search for proxies: simple tests that, whensatisfied, imply a structure is rigid [2–5]. Their existence ∗ To whom correspondence should be addressed:[email protected] † To whom correspondence should be addressed:[email protected] indicates that, in a practical sense, identifying energeti-cally rigid structures is possible in many, though not allcases.The standard proxy for rigidity in particulate systemscomes from Maxwell [6]. When two particles interact,for example through a contact, that interaction con-strains each particle’s motion. If a system has fewerconstraints than the particles have degrees of freedom,it is said to be under-constrained and therefore one ex-pects it to be floppy. This thinking has been successfullyapplied to many examples of athermal systems, such asjammed granular packings, randomly diluted spring net-works, and stress diluted networks [7–10]. A straight-forward extension of Maxwell’s argument, known as theMaxwell-Calladine index theorem [2, 11], shows that oneshould also subtract the number of states of self stress,equilibrium states of the system that can carry a load, be-cause they arise from redundant constraints. In hinge-barnetworks, these ideas can be exploited to design mechan-ical metamaterials with topologically protected mecha-nisms [11–15]. Indeed, the Maxwell-Calladine index the-orem applies generically to spring networks, that is foralmost every set of equilibrium edge lengths [5].A strong challenge to using this constraint countingproxy to determine rigidity is provided by recent ob-servations of rigidity transitions in experiments [16–20]and models [21–30] for biopolymer networks, and exper-iments [31–37] and models [38–44] for cellularized bi-ological tissues. These rigidity transitions do not in-volve changes to constraints or network topology, andare driven instead by tuning a continuous control pa-rameter – the applied strain in biopolymer networks andcell shape in biological tissues. Merkel et al. recentlysuggested that both types of rigidity transitions can beunderstood in terms of a geometric incompatibility be-tween the constraints and the physical space available a r X i v : . [ c ond - m a t . s o f t ] F e b to the system, and showed that the transition coincideswith the appearance of a system-spanning state of selfstress [29]. For finite-size systems, the state of self stresscan be used to predict scaling of elastic properties nearthe transition, though other numerical results suggestthat different scaling may arise in the thermodynamiclimit [45]. Materials that can be rigidified by tuninga continuous parameter provide an attractive path for-ward for designing metamaterials with programmed anddynamic properties [46, 47].So how and why does naive constraint counting fail?Seminal work by Yan and Bi [44] emphasizes that therank of the full Hessian matrix, the matrix of secondderivatives of the energy including effects beyond first-order perturbations, determines rigidity in vertex mod-els for cellularized biological tissues. In particulate sys-tems, weakly aspherical ellipsoids that are highly under-constrained show stable packings where finite-amplitudemotions cost energy [8, 48–50]. Even in over-constrainedelastic networks, prestresses have been shown to affectthe stability of the system [51]. Moreover, triangulatedorigami requires one to consider higher order perturba-tions to the constraints to correctly predict the numberof degrees of freedom [52]. Taken together, these analyseshighlight that rigidity is fundamentally nonlinear.The importance of nonlinear rigidity has already beenhighlighted by mathematicians and engineers in the con-text of bar-joint frameworks, origami, and tensegrities[3–5, 52–54]. In such systems, the focus has been on so-called ”structural rigidity” – whether the geometric con-straints are preserved – rather than on ”energetic rigid-ity” – whether the energy is preserved. In particular,Connelly and Whitely [4] demonstrate that in tensegri-ties there exist states where a different proxy, termed“second-order rigidity”, is sufficient to ensure that theconstraints are preserved. Importantly, that work focuseson the existence of such states, which do not have to beinitialized in a state at the local minimum of an energyfunctional.However, in many physical systems of interest, the dy-namics or boundary conditions that drive the system to-wards rigidity select specific, non-generic states at thetransition point [55]. Therefore, instead of demonstrat-ing the existence of states where second-order rigidity im-plies structural rigidity, we instead ask a different ques-tion: what can we say about energetic rigidity for sys-tems that are at an energy minimum and correspond tohighly non-generic states selected by physical dynamics?In particular, is it possible to find or design structureswhere motions preserve the energy but not the individ-ual constraints? In an important sense, such a structurewould still be floppy.To answer this question we develop a generalized for-malism for understanding the rigidity of energetically sta-ble physical materials, extending previous work on thetheory of second-order rigidity. Specifically, we demon-strate that the onset of rigidity upon tuning a contin-uous parameter emerges from the effects of geometric incompatibility arising from higher-order corrections toMaxwell-Calladine constraint counting. Depending onthe prestresses in the system and features of the eigen-value spectrum, we identify different cases where first-order or second-order rigidity imply energetic rigidity.We then apply the formalism to three systems: under-constrained spring networks, vertex model, and jammedpackings of 2D disks. We show that our formalism canexplain the origin of rigidity in these systems and discussnew research directions suggested by this understanding. II. FORMALISM
In this section, we will discuss what makes a phys-ical system rigid. We assume the state of the systemis described by N dof generalized coordinates, x n . Forexample, the coordinates { x n } might represent the com-ponents of the positions of all vertices in a spring net-work. We also introduce M strains of the form f α ( { x n } )and assume the physical system is characterized by theHooke-like energy, E , of the form E = 12 M (cid:88) α =1 k α f α ( { x n } ) , (1)where k α > f α then serve as theconstraints in Maxwell-Calladine counting arguments.Without loss of generality, we absorb k α into f α by rescal-ing it by √ k α and writing E = (cid:80) Mα =1 f α /
2. Throughoutthis paper, we will assume that the system is at a localminimum of energy, meaning that no small perturbationscan lower the energy.As an example, for a d − dimensional spring networkof N vertices connected via M springs with rest length L in a fixed periodic box, N dof = dN and the strainassociated with spring α connecting vertices i and j atpositions X i and X j is simply the strain of the spring, f α = L α − L , where L α = | X i − X j | is the actual lengthof the spring.We can capture the intuitive notion of rigidity or flop-piness by considering the behavior of Eq. (1) under de-formations. A system is energetically rigid if any global motion that is not a trivial translation or rotation in-creases the energy. A global motion is one that extendsthrough the entire system so as to exclude rattlers ordanglers. If there exists a nontrivial, global motion thatpreserves the energy, we call the system floppy . If, fora given system at an energy minimum, all the strainsvanish, f α = 0 for all α , and the system is unstressed .Otherwise, we say the system is prestressed . A. Standard proxies of energetic rigidity
Experimentally, the standard proxy used to determinewhether the system is energetically rigid is the shearmodulus, G , defined as the second derivative of energywith respect to a shear variable γ in the limit of zeroshear [37, 43]: G = 1 V d E d γ = 1 V (cid:32) ∂ E∂γ − (cid:88) l λ l (cid:34)(cid:88) n ∂ E∂γ∂x n u ( l ) n (cid:35)(cid:33) , (2)where V is the volume of the system while λ l and u ( l ) n arerespectively the eigenvalues and eigenvectors of the Hes-sian matrix, H nm = ∂ E/∂x n ∂x m , and the sum excludeseigenmodes with λ l = 0. When G (cid:54) = 0, the system is cer-tainly energetically rigid. Note that this is closely alliedwith the mathematical notion of prestress stability [4](see Appendix A). On the other hand, if H nm has global,nontrivial zero eigenmodes (or more precisely, zero eigen-modes that overlap with the shear degree of freedom), G = 0 [43].Importantly, defining rigidity based on G is not equiv-alent to energetic rigidity. Specifically, G (cid:54) = 0 implies thesystem is energetically rigid, but G = 0 does not implyfloppiness. As highlighted in Appendix A there may bequartic corrections in δx n that increase the energy evenwith vanishing shear modulus. Moreover, in many casesof interest these quartic corrections are expected to dom-inate precisely at the onset of rigidity.A definition of rigidity based on G is equivalent toexamining the Hessian matrix H directly: if H is positivedefinite on the global, non-trivial deformations, then thesystem is also energetically rigid. Writing out the Hessianmatrix in terms of the constraints, we find H nm = ∂ E∂x n ∂x m = (cid:88) α (cid:20) ∂f α ∂x n ∂f α ∂x m + f α ∂ f α ∂x n ∂x m (cid:21) = ( R T R ) nm + P nm , (3)where R αn = ∂f α ∂x n (4)is known as the rigidity matrix. We call ( R T R ) nm theGram term (as it is the Gramian of rigidity matrix), and P nm the prestress matrix because it is only nonzero if f α (cid:54) = 0 (Gram term and prestress matrix are sometimescalled stiffness matrix and geometric stiffness matrix re-spectively in structural engineering [4, 53]). If the Hes-sian has at least one global nontrivial zero direction, we obtain the necessary (but not sufficient) condition forfloppiness, (cid:88) nm P nm δx n δx m = − (cid:88) nm ( R T R ) nm δx n δx m = − (cid:88) α (cid:32)(cid:88) n ∂f α ∂x n δx n (cid:33) , (5)where the sum over α is over all constraints and, again,trivial Euclidean modes have been excluded. Analogousto our discussion of G above, a definition of rigidity basedon H is also not equivalent to energetic rigidity, due tothe importance of quartic terms in cases of interest (in-cluding at the transition point). B. First- and second-order rigidity tests
The existence of any global, non-trivial, and continu-ous motion of the system x n ( t ) that preserves the con-straints f α ( { x n ( t ) } ) implies the system is floppy. A sys-tem is structurally rigid when no such motions exist, adefinition highlighted in Table I. Energetic rigidity is notnecessarily equivalent to structural rigidity when the sys-tem is prestressed ( E > E = 0, as discussed in more detail later.Though determining whether a system is structurallyrigid is NP-hard [1], there are several simpler conditionsthat, if they hold true, imply that a system is structurallyrigid [2–5]. These tests (and more) are reviewed in moredetail in Appendix A and summarized in Table I. In thissection, we turn to the question: under what conditionsdo these structural rigidity tests, valid at first and secondorder in the vertex displacements, also imply energeticrigidity? We will conclude with classification of systemsinto those for which structural rigidity at first or secondorder is sufficient to imply energetic rigidity and thosefor which the question is still not easily determinable.We first consider a first-order perturbation to the con-straints, δf α = (cid:80) n ∂f α /∂x n δx n . We define a linear(first-order) zero mode (LZM) δx (0) n as one that preserves f α to linear order, (cid:88) n ∂f α ∂x n δx (0) n = (cid:88) n R αn δx (0) n = 0 . (6)We can see that LZMs are in the right nullspace ofthe rigidity matrix. Excluding Euclidean motions, anontrivial LZM is often called floppy mode (FM) inphysics [11]. The Maxwell-Calladine index theorem (alsoknown as the rigidity rank-nullity theorem) states that N dof − M = N − N s , where N is the number of LZMsand N s is the number of states of self stress [2]. A stateof self stress is a vector σ α in the left nullspace of therigidity matrix: (cid:88) α σ α R αn = 0 . (7) A system is . . . when . . .Energetically rigid any nontrivial global motion increases the energyStructurally rigid no nontrivial global motion preserves the constraints f α First-order rigid no nontrivial global motion preserves the constraints f α to first orderSecond-order rigid no nontrivial global motion preserves the constraints f α to second orderTABLE I. Different definitions of rigidity. For a system to be prestressed, it must have a state of selfstress. In fact, if we write the Euler-Lagrange equationsfor the system at the energy minimum: (cid:88) α f α ∂f α ∂x n = (cid:88) α f α R αn = 0 , ∀ n. (8)Therefore, if f α (cid:54) = 0, f α has to be a state of self stress.Note, however, the converse is not true. The existenceof states of self stress only depends on the geometry ofthe system and does not imply that the system has to beprestressed. For example, take a system with constraints f α ( { x n } ) = F α ( { x n } ) − F α at a particular mechanicallystable configuration { ¯ x n } that has a state of self stressand choose F α = F α ( { ¯ x n } ). The system will be un-stressed at { ¯ x n } but still has a state of self stress.Maxwell constraint counting suggests that an over-constrained system ( N dof < M ) must be rigid while anunder-constrained system ( N dof > M ) must be floppy.Generic states in over-constrained systems have no FMs.Such a system is called first-order rigid and, indeed, insuch systems first-order rigidity implies structural rigid-ity as defined in Table I [3, 4]. However, Eq. (5) leavesopen the possibility that a system with negative pre-stresses could be first-order rigid with a zero shear mod-ulus. Conversely, when the system is prestressed, thesystem may be energetically rigid even if it is not first-order rigid.We thus study the variations to second order. If thereare motions that preserve f α to second order in δx n , Tay-lor expansion of f α results in: δf α ≈ (cid:88) n R αn δx n + 12 (cid:88) nm ∂ f α ∂x n ∂x m δx n δx m = 0 , (9)where we used Eq. 4 for the linear term in the expansion.If the only LZMs that satisfy Eq. 9 are trivial ones, thesystem is called second-order rigid and, consequently, isstructurally rigid [3, 4]. It can be shown that a LZM, δx (0) n , must satisfy (cid:88) α (cid:88) nm σ α,I ∂ f α ∂x n ∂x m δx (0) n δx (0) m = 0 , (10)for all states of self stress σ α,I and solutions to Eq. 7 tobe a second-order zero mode ([4, 5]; Appendix A).Finally, we note that going beyond second order is lesshelpful than one might suppose. There are examples ofsystems that are rigid only at third order or beyond yetremain floppy [56]. We now ask the question: when does second-orderrigidity imply energetic rigidity? We identify three cases,which encompass several examples of physical interest,where second-order rigidity implies energetic rigidity, anddemonstrate that second-order rigidity is a better proxyfor energetic rigidity than the shear modulus. We iden-tify a fourth case where second-order rigidity does notimply energetic rigidity – for example there may be sys-tems with large prestresses that do not preserve f α tosecond-order but preserve energy. We classify these dis-tinct cases using the eigenspectrum of P nm and the statesof self stress. In all the cases, we will assume that if thesystem has FMs, at least one is global. Case 1: The system has no self stresses ( N s = 0 ) In this case, constraint counting alone provides a goodproxy for energetic rigidity. From Eq. (8), the system isalso unstressed and Eq. (5) reduces to (cid:88) α ( (cid:88) n ∂ n f α δx n ) = 0 . (11)The solutions are LZMs, δx (0) n (Eq. (6)). If a systemdoes not have any FMs, it is energetically rigid. An en-ergetically rigid system with no states of self stress isalso called isostatic. This also means that there are nomotions that preserve f α even to first order, thus thesystem is first-order rigid. Examples of systems belong-ing to Case 1 include under-constrained spring networks,unstressed vertex models with no area terms (which wewill study below), and the special, non-generic framesdescribed in Fig. 4(a)-(c) of [11]. Case 2: The system has at least one self stress( N s ≥ ) While configurations with self stresses are, in princi-ple, rare, they still occur quite commonly in physical sys-tems. Indeed, any local extremum of the energy necessar-ily leads to a self stress (Appendix A) and, consequently,any system whose configurations are determined by min-imizing an energy can easily find itself in a state with oneor more self stresses. As seen in Eq. (10), the existenceof states of self stress puts additional constraints on zeromodes. Here, we explore the implications of these con-straints on energetic rigidity. To do so, we look at severalsub-cases:
Case 2A: The system is unstressed ( P nm = 0 ) This case includes systems with either no prestress, f α = 0, or systems for which the prestress is perpen-dicular to its second order expansion such that P nm = (cid:80) α f α ∂ n ∂ m f α = 0. If the system is first-order rigid, it isagain energetically rigid. If there are global FMs, G = 0;however, it can be shown (Appendix A) that the fourthorder expansion of energy for these modes will be δE ≈ N s (cid:88) I =1 (cid:34) (cid:88) α,nm σ α,I ∂ n ∂ m f α δx (0) n δx (0) m (cid:35) (12)Therefore, if the system is second-order rigid in the spaceof its global FMs, it is energetically rigid even though G = 0. Examples include random regular spring net-works with coordination number z = 3 and vertex modelsexactly at the rigidity transition. Case 2B: P nm is positive semi-definite For a system with a positive semi-definite P nm , theHessian has a zero eigenmode if and only if both LHSand RHS of Eq. (5) are zero for δx n . The RHS is zeroonly for LZMs. Then if the system is first-order rigid,it is again energetically rigid. For a system with globalFMs, we reduce Eq. (5) to (cid:88) nm P nm δx (0) n δx (0) m = (cid:88) nm (cid:88) α f α ∂ n ∂ m f α δx (0) n δx (0) m = 0 , (13)where x (0) n is now a global FM. We show below thatsecond-order rigidity implies energetic rigidity, but de-pending on N s , G may be zero. If the system has a single self stress : Calling thisstate of self stress σ α , we conclude from Eq. (8) that f α ∝ σ α , meaning Eq. (13) is identical to Eq. (10) in thiscase. This means that if this system is second-order rigid,it is energetically rigid and G >
0. We will demonstratein Section III that both spring networks under tensionand vertex models with only the perimeter term fall intothis category.
If the system has multiple self stresses : In Ap-pendix A we show that if the system is second-orderrigid in the space of global FMs, it is energetically rigid(Eq. (12)). However, the Hessian may still have zeroeigenmodes if in the minimized state f α is a linear com-bination of self stresses that satisfies Eq. (13). This sug-gests that the system may be energetically rigid but with G = 0. We have not been able to identify an exampleof a second-order rigid system with multiple self stressesand G = 0, but if one exists, it may lead to interestingideas for material design. Case 2C: P nm has negative eigenvalues In this case, we have been unable to derive analytic re-sults for whether first-order or second-order rigidity im-plies energetic rigidity. As the models that fall into thisclass are quite diverse, it is likely that more restrictiveconditions are necessary in specific cases to develop ana-lytic results.One example in this category is vertex models withan area term in addition to a perimeter term when pre-stressed. In Section III, we demonstrate numerically thatin such models there is always only one state of self stressthat is non-trivial, and that P nm has negative eigenval-ues. However, the Hessian itself is still positive-definite(excluding trivial LZMs) and therefore the system is en-ergetically rigid. Another example is a rigid jammedpacking, which exhibits quite different behavior for theeigenspectra of P nm .More generally, we cannot rule out the possibility thatthere may be examples where the Hessian of a first-orderor second-order rigid system could have global zero di-rections for non-zero modes. Such a system would bemarginally stable because if any negative eigenmode of P nm becomes too negative, the Hessian would have anegative direction and the system would not be at anenergy minimum anymore. Furthermore, states of selfstress place the same constraints as in Eq. (10) on thesenon-zero modes. If those constraints are not satisfied, theenergy would increase at fourth order (Appendix A), sug-gesting that again the shear modulus could be zero whilethe energy is not preserved. Even though it is highlynongeneric, this case could aid in the design of struc-tures that become unstable by varying the prestress [51]or new materials that are flexible even though individualconstraints are not preserved.Fig. 1 summarizes the cases describing when eitherfirst-order or second-order rigidity imply energetic rigid-ity. In Appendix A, we provide another flowchart (Fig.7) to clearly establish the connection between energeticrigidity and structural rigidity as understood by mathe-maticians. We also provide several propositions to showthat energetic rigidity and structural rigidity are inter-changeable when E = 0 but not necessarily otherwise.For instance, it can be shown that first-order and second-order rigidity both imply structural rigidity [5], but wesaw that they do not always imply energetic rigidity. Thisis because for a system which possesses self stress at anenergy minimum, mathematicians only require the exis-tence of a linear combination of self stresses that wouldmake the system rigid [4], however, that particular selfstress may not be the linear combination of self stressesthat the system chooses as its prestress based on externalforces [55].In the next section, we apply this formalism to studyrigidity in three important examples previously discussedin the literature: 2D spring networks, vertex model, andjamming. 𝑁 " = 0 𝑁 " > 0𝒫 ’( ≥ 0 𝒫 ’( has negative eigenvalues 𝑁 *(,,’.) > 0 st order rigidityenergetic rigidity1 st order rigidityenergetic rigidity 2 nd order rigidityenergetic rigidity System dependentDoes either first-order or second-order rigidity imply energetic rigidity? 𝑁 *(,,’.) = 0 FIG. 1.
Flowchart of Cases summarizing the classificationof systems based on the findings of second-order rigidity for-malism. N ( g,nt )0 refers to the number of global non-trivialLZMs (i.e. global FMs). III. EXAMPLES
We use the formalism introduced in the previous sec-tion to study two examples of under-constrained systems,2D random regular spring networks and vertex models,from analytical and numerical perspectives and examinewhether their linear zero modes make them floppy. Wethen briefly examine rigidity of jammed packings in thecontext of this new formalism.
A. 2D Spring Network
Our first example is a 2D spring network comprised of N vertices connected by springs with coordination num-ber z = 3 in a periodic square of fixed length. We choose z = 3 so that the system is under-constrained and thenetwork is regular, and there are no dangling vertices, butthe results are valid for any z < z > L .For this system, N dof = 2 N and M = 3 N/
2. From con-straint counting, the number of LZMs, N ≥ N/
2, sothere are many system spanning FMs and one might beinclined to conclude that the system is floppy. However,we will see that for a range of spring rest lengths, thesystem possesses a self stress and becomes energeticallyrigid.
1. Analytical Results
Calling the length of the α th spring L α , the energy ofthis spring network is E = K L N/ (cid:88) α =1 ( L α − L ) , (14)which defines the constraints f α = L α − L . Here, springconstants are identical and equal to 2 K L . In Appendix B,we show that the condition for FMs to be second-order(Eq. (10)) can be written as (cid:88) α σ α,I ( δ L α ) L α = 0 . (15)Assuming the edge α connects vertices n and m , δ L α = δ X (0) n − δ X (0) m is the vectorial displacement of the edgein response to FM { δ X (0) n } . Since trivial LZMs are ex-cluded, this displacement must be perpendicular to theedge, δ L α = δ L ⊥ α . An important observation is that ifthere exists a positive self stress ( σ α,I > G = 0condition for spring networks is (Eq. 5) (cid:88) α ( L α − L ) (cid:0) δ L ⊥ α (cid:1) L α = − (cid:88) α (cid:16) δL (cid:107) α (cid:17) , (16)for any global mode of motion, where δL (cid:107) α is the com-ponent of δ L α parallel to the edge. The left-hand sideof Eq. (16) is equivalent to (cid:80) mn P mn δx n δx m from theprevious section.Next, using numerics we show that even though thesystem has at least N/ L ≤ L ∗ which implies energetic rigidity.
2. Numerical Results
We simulate the spring network by a random Voronoitessellation of space in two dimensions which is guar-anteed to produce networks with z = 3. Details aregiven in Appendix D. Defining the rigidity matrix as R αm = ∂f α /∂ X m , where X m is the position vector ofvertex m , the LZMs are found by solving for the zeromodes of R T R . The states of self stress similarly arezero modes of RR T . Fixing the box size, we vary L andminimize the network, observing that the system goesthrough a transition from f α = 0 for L > L ∗ to f α > L < L ∗ .From constraint counting, the number of LZMs N mi-nus the number of states of self stress N s is N − N s = N/
2, which is confirmed numerically: we find N = N/ N s = 0 for L > L ∗ , and N = N/ N s = 1for L < L ∗ (Fig. 2). There are two trivial LZMs (no (a)(b) FIG. 2.
Spring Network Density of states D as a functionof the eigenvalue λ for the Hessian matrix (solid black line),Gram matrix (dashed yellow line) and prestress matrix (dash-dotted red line), averaged over 10 samples with N = 1000 and K L = 2. (a) In the Floppy regime ( L = 0 . N/ L = 0 . N/ N s = 1, while the Hessian onlypossesses two trivial (translational) zero modes. The HessianDOS is dominated by the Gram term at high frequency andby the prestress matrix at low frequencies. rotation is allowed because of periodic boundary condi-tions).Consistent with previous work [37, 43], we use theshear modulus to quantify the energetic rigidity of thesystem. Specifically, we first calculate the eigenmodesof the Hessian, which then allows us to calculate theshear modulus G . From the observation that f α = 0for L > L ∗ (Case 1) and that there are many FMs,we would expect the system to be floppy in this regime.We find that, indeed, G = 0 for L > L ∗ , agreeing withour prediction (Fig. 3a), and this is also consistent withpreviously reported results [29].For L < L ∗ , the Hessian has no nontrivial zero eigen-modes. We can see that by examining the eigenvaluesof P mn given by the LHS of Eq. 16, which is positivesemi-definite since f α >
0. The high-frequency spectrumremains nearly identical to that in the unstressed case,which makes sense because the Gram term that domi-nates at high frequencies depends only on geometry, anddoes not change significantly as L is lowered in the rigidregime. In contrast, in the rigid phase a new band of low-frequency modes appears in the Hessian, clearly comingfrom the prestress eigenspectrum. Moreover, the average magnitude of this low-frequency band is very sensitiveto the control parameter L – eigenvalues shift to higherfrequencies as L is lowered (Fig. 3b).Numerical results indicate that the system only has onestate of self stress, consistent with previous work [29].Therefore it falls under Case 2B, where second-orderrigidity would imply energetic rigidity, and Eq. (5) re-duces to Eq. (15) with σ α,I → f α . Due to the exis-tence of the positive self stress, the system is second-orderrigid and energetically rigid. Exactly at the transition L = L ∗ , the system is unstressed and G = 0 (see [29]).However, at L → L ∗− , since the system has a positiveself stress and is second-order rigid, it is energeticallyrigid as well (Case 2A). B. 2D Vertex Model
Here we discuss rigidity of another highly under-constrained system, the 2D vertex model. The 2D vertexmodel consists of N cell polygonal cells tiling a 2D periodicsquare. The energy of the system is E = N cell (cid:88) α =1 (cid:2) K A ( A α − A ) + K P ( P α − P ) (cid:3) , (17)where A α and P α are the area and perimeter of the α th cell, respectively. We have assumed that A α has a pre-ferred value A and P α has a preferred value P . Theenergy is still of the form Eq. (1) only with two sets ofconstraints f α, = A α − A and f α, = P α − P . The totalnumber of constraints is thus M = 2 N cell . In the vertexmodel, DOFs are the vertices. Thus, in a periodic box, N dof = 4 N cell and we have at least 2 N cell LZMs fromconstraint counting. These constraints are not all inde-pendent, however, because they act on the same vertices.In the numerical section, we will show by looking at therank of the rigidity matrix that there is in fact one re-dundant constraint and there is a state of self stress andan extra zero mode because of it. We will also show thatfor P < P ∗ , a second state of self stress appears due togeometric incompatibility, similar to spring networks.In the analytical results section, we limit ourselves tothe more tractable version of the model with no areaconstraints so that K A = 0, and in the numerical sectionwe will study the general version given by Eq. (17).
1. Analytical Results
We look the equations governing second-order zeromodes and the G = 0 condition (Eq. (5)) for the ver-tex model with no area term ( K A = 0), thus M = N cell .The constraints on the vertices are given by f α = P α − P .The self stresses of R αm impose the following quadraticconstraints on LZMs, δ X n (Appendix C): (a)(b)(c)(d) FIG. 3.
Spring Network Comparison of Density ofStates as a function of the eigenvalue λ for the Hessian ma-trix ( D H ), Gram matrix ( D G ) and prestress matrix ( D P ),averaged over 10 samples with N = 1000 and K L = 2 for dif-ferent values of L . (a) Shear modulus ( G ) of spring networksas a function of L shows a transition from fluid to solid at L ∗ ≈ .
63. The error bars represent the standard deviationover 10 samples. (b) Transitioning from floppy to rigid coin-cides with the appearance of low frequency eigenmodes in theHessian DOS that shift to higher frequencies as the system be-comes stiffer. (c) Even though in the rigid regime there is onlyone new zero mode, the Gram matrix DOS looks significantlydifferent when compared to systems with L < L ∗ becausethe geometry of the system does not change significantly inthe rigid regime. (d) Unlike the Gram matrix, the prestressincreases linearly as L is decreased in the rigid regime, asreflected by the shift in the prestress matrix DOS towardshigher frequencies. (cid:88) α σ α,I (cid:88) edge j ∈ cell α (cid:0) δ L ⊥ j (cid:1) L j = 0 . (18)Since all vertices are connected to three edges and thebox size is fixed, for a generic system, there are no non-trivial motions that do not introduce a δ L ⊥ j and thus theinner sum is positive definite. Hence, similar to springnetworks, if a self stress σ α,I > (cid:88) α f α (cid:88) edge j ∈ cell α (cid:0) δ L ⊥ j (cid:1) L j = − (cid:88) α (cid:32)(cid:88) m R αm · δ X m (cid:33) . (19)We will see in the numerical results section that similarto spring networks, vertex model with K A = 0 has apositive state of self stress for P ≤ P ∗ and thus is bothsecond-order rigid and energetically rigid. For P > P ∗ however, all constraints are satisfied and the system hasno state of self stress (Case 1). Therefore the system isfloppy because it has many (3 N cell ) LZMs. We will alsosee that vertex model with K A = 1 belongs to Case 2Cwhen prestressed ( P < P ∗ ), but it is still second-orderrigid and energetically rigid.
2. Numerical Results
We simulate the vertex model (with both area andperimeter terms) using the same algorithm as springnetworks but with the energy function given in Eq. 17with A = 1 and varying P in a fixed periodic box( A box = N cell ). Details, along with numerical resultsfor systems with no area constraints ( K A = 0), are givenin Appendix C and D. Each component of the rigiditymatrix R αm now is a 2 × α and two spatial components for each vertex m .We find that for P > P ∗ , all constraints are satisfiedand the Hessian zero eigenmodes are the same as theLZMs ( N = 2 N cell + 1) (Fig. 4a). Even though thereare no prestresses, we find N s = 1 (satisfying the rank-nullity theorem). To understand the source of this selfstress, note that the sum of all the cell areas must beequal to the fixed box size. Thus A merely changes theoverall pressure of the system [57] and we can rewrite theenergy function as: E = N cell (cid:88) α =1 (cid:20) K A ( A α − A box N cell ) + K P ( P α − P ) (cid:21) + K A N cell ( A box N cell − A ) . (20) (a)(b) FIG. 4.
Vertex Model Density of states D as a functionof the eigenvalue λ for the Hessian matrix (solid black line),Gram matrix (dashed yellow line) and prestress matrix (dash-dotted red line), averaged over 10 samples with N cell = 500and K A = K P = 1. (a) In the Floppy regime ( P = 3 . N cell + 1 zeromodes and their DOS curves overlap as there is no prestressin the network. N s = 1 which is the trivial area self stress.(b) In the rigid regime ( P = 3 . N cell + 2 zero modes and N s = 2, while the Hessian onlypossesses two trivial (translational) zero modes. The HessianDOS is dominated by the Gram term at high frequency and bythe prestress matrix at low frequencies. While the prestressmatrix has some negative eigenvalues, the vast majority arepositive. Therefore, changing A amounts to increasing the over-all tissue pressure without breaking force balance andby definition there must be a self stress associated withit. Indeed, we see that the area components of selfstress are 1 while its perimeter components are 0: σ =(1 , , . . . , , N cell − P < P ∗ , none of the constraints are satisfied dueto geometric incompatibility (Fig. 4b). The perimeterprestresses are all positive ( f α, > A = 1 and A box = N cell , theaverage, ¯ f α, = 0. We find N s = 2 and N = 2 N cell + 2.One of the self stresses is due to the geometric incompat-ibility we already saw in the spring network, σ = ( A − A box /N cell , P − P , . . . , A N cell − A box /N cell , P N cell − P ).The other one arises from the area constraint redundancy, σ = (1 , , . . . , , G = 0 for P > P ∗ , and G > P < P ∗ (Fig. 5a) suggesting that the system is second-order rigid when σ appears. Decreasing P further shifts the Hes-sian eigenmodes to higher frequency and increases G sim-ilar to spring networks (Fig. 5b). This is again due to ashift in P nm and not the Gram term (Fig. 5c, d).Looking at the spectrum of P nm , we can see that ithas negative eigenvalues (Fig. 5d) and thus the systemfalls under Case 2C. It is empirically still rigid because P nm does not have directions negative enough to satisfyEq. (5). Even though its negative eigenvalues becomemore negative as P is decreased, its positive ones domi-nate (Fig. 5c).One way to identify whether it is the Gram or prestressmatrix that is responsible for rigidity is to multiply theprestress matrix P nm by an arbitrary (cid:15) >
1. This sug-gests that the Hessian is dominated by the positive eigen-values of the prestress matrix and not from a competitionwith the Gram term. For K A = 0 with P < P ∗ , P nm ispositive semi-definite and N s = 1, so it falls under Case2B similar to spring networks. At the onset of rigidity, P → P ∗− , both vertex models will show G →
0, butsince they both have a nontrivial self stress (Case 2A)and are second-order rigid, our formalism indicates theyare energetically rigid as well.
C. 2D Jammed Packings
Athermal packings of soft or hard spheres are a usefulmodel for studying granular matter and glasses at zerotemperature. A 2D disk packing with one state of selfstress is in a way very similar to the spring networksthat we studied above. States of self stress in jammedsystems are extended over the entire system with com-pressive forces everywhere [9] which resembles the case ofspring networks under tension. Here, the energy is givenby: E = k N − (cid:88) α =1 h α Θ( h α ) , (21)where h α = (1 − ρ α σ α ) is the dimensionless overlap betweenparticle pair α ≡ i, j , with ρ α being the distance betweentwo disk centers and σ α being the sum of their radii. TheHeaviside step function is used to count contributionsfrom positive overlaps only.
1. Analytical Results
The analysis presented in section III A can also be usedto describe the energetic rigidity in jammed packings ofsoft harmonic disks/spheres. One difference between aspring network under tension and a critically jammedpacking under compression (with one state of self stressonly) is that all of the prestress forces in a jammed pack-ing are negative and therefore any terms in the expansion0 (a) (b)(c)(d)
FIG. 5.
Vertex Model Comparison of Density of States as a function of the eigenvalue λ for the Hessian matrix ( D H ),Gram matrix ( D G ) and prestress matrix ( D P ), averaged over10 samples with N cell = 500 and K A = K P = 1 for dif-ferent values of P . (a) Shear modulus ( G ) of vertex modelwith as a function of P shows a transition from fluid to solidat P ∗ ≈ .
84. The error bars represent standard deviationover 10 samples. (b-d) Comparison between DOS of Hessianand its components for different values of P . (b) Transi-tioning from floppy to rigid coincides with the appearance oflow frequency eigenmodes in the Hessian DOS that shift tohigher frequencies as the system becomes more rigid. (c) Eventhough in the rigid regime there is only one new zero mode,the Gram matrix DOS looks significantly different when com-pared to systems with P < P ∗ because the geometry of thesystem does not change significantly in the rigid regime. (d)Unlike the Gram matrix, the prestress increases linearly as L is decreased in the rigid regime, as reflected by the shiftin the prestress matrix DOS towards higher frequencies (bothpositive and negative). of the energy that are proportional to the first deriva-tives will be negative. Jammed packings thus belong toCase 2C. We do not present the analytical work for thesesystems here because it will be a repetition of the calcu-lations in section III A 1.
2. Numerical Results
We create an ensemble of ten 2D disk packings veryclose to the critical jamming using standard methods asdescribed in Appendix D, so that there is one state ofself stress in each packing. We calculate the eigenspectraof the Hessian matrix as well as the Gram and prestressterms; the results are presented in Fig. 6.Unlike spring networks that have many non-zero modesin their floppy regime, a critically jammed packing willcompletely unjam and reach a global minimum with zeroenergy and zero eigenmodes everywhere if the densityis lowered. Therefore, we do not report data from thefloppy side of the transition. As can be seen from Fig. 6,the density of states of the full Hessian almost matchesthe density of states of Gram term in the Hessian. Thisis because at critical jamming, the prestress forces areinfinitesimal and the magnitude of the eigenvalues of pre-stress term are several orders of magnitude smaller thantheir equivalent eigenvalues in the Gram term. At thispoint, the rigidity of the system is mainly determined bythe Gram term in the Hessian. Since both Gram andprestress terms have two zero eigenmodes, the full Hes-sian also has two zero eigenmodes which is typical of arigid body in 2D.A consequence of this is that the energetic rigidityof jammed systems can be fully described using theMaxwell-Calladine count, since even at the jammingpoint where the pressure is zero and the prestress forcesare infinitesimal, the system is first-order rigid. The pre-stress forces can only play a role in the energetic rigidityof the system when the pressure is large enough to pushthe system to an instability. This marks another differ-ence between the spring networks under tension and softharmonic particles under compression at one state of selfstress.
IV. DISCUSSION AND CONCLUSIONS
We term an “energetically rigid” structure as onewhere any sufficiently small applied displacement in-creases the structure’s energy. Our focus on motionsthat preserve energy contrasts with previous work onstructural rigidity that has focused on motions that pre-serve constraints. There are interesting differences be-tween the two approaches. Unlike structural rigidity, en-ergetic rigidity is not defined solely by the geometry –predictions also depend on the energy functional. Herewe studied a Hooke-like energy that is quadratic in theconstraints, which is the simplest nontrivial energy func-1
FIG. 6.
Jammed packing Density of states D as a func-tion of the eigenvalue λ for the Hessian matrix (solid blackline), Gram matrix (dashed yellow line) and prestress matrix(dash-dotted red line), in the rigid regime (at one state of selfstress) averaged over 10 samples with N = 1000 and k = 1.Each of these three data sets only have 2 zero modes persample. tional that encompasses a large number of physical sys-tems, but other choices are possible. On the other hand,this choice opens the possibility that in some structuresthere may exist motions that preserve the energy with-out preserving individual constraints. Importantly, theframework developed here would allow us to identify suchsystems as floppy.Specifically, we want to understand under which pre-cise circumstances structural rigidity implies energeticrigidity, and in the process identify underlying geomet-ric mechanisms that are responsible for rigidity in spe-cific materials. It is understood that predicting whethera planar graph is structurally rigid is already an NP-hard problem, and so previous work has proposed several“quick” tests for rigidity, which work in limited circum-stances. One test is the Maxwell-Calladine index theo-rem, also called first-order rigidity, which tests whetherthe constraints f α that define the energy functional canbe satisfied to first order. Another test is second-orderrigidity, which checks whether constraints can be satisfiedto second order.In this work we have developed a systematic frameworkthat clarifies the relationship between energetic rigidityand these other previously proposed rigidity tests. Wedemonstrate that first-order rigidity always implies ener-getic rigidity when there are no states of self stress. How-ever, when the system does possess states of self stress,the eigenvalue spectrum of the prestress matrix P nm con-trols whether first- or second-order rigidity (or neither)implies energetic rigidity. For several physical systemsof interest, we show second-order rigidity is sufficient toguarantee energetic rigidity, while in others it is not.Specifically, in both spring networks and vertex modelswithout an area term, where the prestress matrix is posi-tive semi-definite, second-order rigidity implies energeticrigidity precisely at the rigidity transition (even thoughthe shear modulus is zero) and in the rigid regime (when the shear modulus is finite).When the prestress matrix is indefinite or negativesemi-definite, we can still show analytically that at therigidity transition, second-order rigidity implies energeticrigidity. But away from the transition point neither first-order nor second-order rigidity guarantee energetic rigid-ity, and indeed we observe widely divergent examples inthis category: in simulations of vertex models with anarea term, second-order rigidity is always seen to implyenergetic rigidity, while in simulations of jammed pack-ings of soft spheres, first-order rigidity is sufficient to pre-dict the onset of energetic rigidity.Moving forward, it would be useful to identify featuresthat distinguish examples in this category, dividing it intosub-cases that are at least partially analytically tractable.One intriguing possibility is to classify a structure’s re-sponse to applied loads. For example, one could artifi-cially increase the prestresses in a structure, multiply-ing P nm by a coefficient (cid:15) >
1, which will only increasethe overall magnitude of the state of self stress but notchange the geometry of the network or the Gram termin the Hessian. Such a perturbation does not destabi-lize vertex models with an area term, where second-orderrigidity is observed to guarantee energetic rigidity, but itdoes destabilize jammed packings where first-order rigid-ity governs energetic rigidity [58].This also suggests that it may be possible to programtransitions between minima in the potential energy land-scape via careful design of applied load. For example,while the type of spring network we studied was com-pletely tensile for L < L ∗ , one could create spring net-works with both tensile and compressed edges [51] or atensegrity with tensile cables and compressed rods. Itwill be interesting to see if we can design such systemsto have a negative-definite prestress matrix. If so, ap-plied loads may destabilize the structure along a specifiedmode towards a new stable configuration.A related question is whether we can move such a sys-tem from one energy minimum to another in a more effi-cient manner. Traditionally, to push a system out of itslocal minimum into a nearby minimum, one rearrangesthe internal components of the system locally or globally,while it is rigid, by finding a saddle point on the energylandscape. An alternate design could be to (1) apply aglobal perturbation that makes the system floppy, (2) re-arrange its components at no energy cost, and (3) applya reverse global perturbation to make it rigid again. Inother words, the fact that the system can transition fromrigid to floppy using very small external forces withoutadding or removing constraints could help us generatereconfigurable materials with very low energy cost.Another interesting avenue for design is to perturb theenergy functional itself. In this work we focused on anenergy that is Hookean in the constraints, but it wouldbe interesting to explore whether different choices of en-ergy functional still generate the same relationships be-tween energetic rigidity and first- or second-order rigidityidentified in Fig 1. If not, such functionals may enable2structures with interesting floppy modes.Taken together, this highlights that the subtleties in-volved in determining energetic rigidity could be ex-ploited to drive new ideas in material design. Withthe framework described here, we now fully understandwhen we can use principles based on first-order con-straint counting or second-order rigidity to ensure en-ergetic rigidity in designed materials. Moreover, theremay be some new design principles available, especiallyfor dynamic and activated structures, if we focus on caseswhere these standard proxies fail. ACKNOWLEDGMENTS
We are grateful to Z. Rocklin for an inspiring initialconversation pointing out the connection between rigid-ity and origami, and to M. Holmes-Cerfon for substantialcomments on the manuscript. This work is partially sup-ported by grants from the Simons Foundation No 348126to Sid Nagel (VH), No 454947 to MLM (OKD and MLM)and No 446222 (MLM). CDS acknowledges funding fromthe NSF through grant DMR-1822638, and MLM ac-knowledges support from NSF-DMR-1951921.
Appendix A: Derivation of second-order rigidity condition and implications for energetic rigidity
Here we first provide a careful connection between our work and the structural rigidity of bar-joint frameworks. Westart with definitions and important theorems and then establish relationships between energetic rigidity and variousrelevant results from the rigidity literature (summarized in Fig. 7). Several of these theorems are adapted from [4].We also provide derivations of second-order rigidity and energetic rigidity that we have omitted from the main text.
1. Basic results on structural rigidity
Let x n be a point in a space of configurations and let F α ( { x n } ) be a set of measures (for example, in a fibernetwork F α ( { x n } ) might give the length of the fibers). From now on we denote the configuration { x n } as x forsimplicity. We start with some basic definitions: Definition:
A nontrivial isometry (or, sometimes, flex) is a one-parameter family of deformations, x ( t ), such that F α ( x ( t )) = F α (for some F α ) and x ( t ) is not a translation or rotation. We similarly refer to a nontrivial deformationas any deformation δx ( t ) that is not a translation or rotation. Definition:
A linear zero mode, also known as a first-order isometry or a first-order flex, at a configuration ¯ x , ˙ x ,is a solution to the equation (cid:80) n ∂ n F α (¯ x ) ˙ x n = 0. A system is first-order rigid if there are no solutions to this equation. Definition:
A self stress, σ α , at ¯ x is a solution to (cid:80) α σ α ∂ n F α (¯ x ) = 0. Definition:
A second-order isometry (or a second-order flex) at ¯ x is a first-order isometry such that (cid:80) α (cid:80) nm σ α,I ∂ n ∂ m F α (¯ x ) ˙ x n ˙ x m = 0 has a solution where { σ α, , σ α, , · · · , σ α,N s } is a basis of self stresses at ¯ x .A system is second-order rigid if it has nontrivial zero modes but no nontrivial second-order isometries.We finally have the main result that we want: a system that is either first-order or second-order rigid, is structurallyrigid [4]. It can be hard – still – to test for structural rigidity at second order because it involves solving a system ofquadratic equations. It is, therefore, convenient to introduce a stronger condition: Definition:
A system is prestress stable at ¯ x if there is a self stress at ¯ x , σ α , such that (cid:80) α σ α ∂ n ∂ m F α (¯ x ) is positivedefinite on every nontrivial zero mode.With this definition, we prove that a system that is prestress stable at ¯ x is also second-order rigid at ¯ x (and hence,structurally rigid) . This follows because there is a self stress σ α such that (cid:80) α σ α ∂ i ∂ j F α (¯ x ) is positive definite onnontrivial first-order flexes. We can construct a basis for the self stresses with σ α as one of its elements. Therefore,it is second-order rigid as well.According to Connelly and Whitely [4], there are examples of second-order rigid structures that are not prestressstable in 2D and, especially, 3D. The notion of prestress stability is related to notions of an energy.3 Energetic rigidity E (¯ x ) = 0 structural rigidityEnergetic rigidity E (¯ x ) > 0 Prestress stability1st order rigidityEnergetic critical point E (¯ x ) > 0 has a self stress ¯ x N S = 1 ∃ F α ∃ F α E ′ ′ (¯ x ) > 0 FIG. 7.
Relations between various definitions for a given configuration ¯ x . The numbers on arrows refers to propositionswith the same numbers. We can see that only when the system is unstressed ( E (¯ x ) = 0), energetic rigidity and structuralrigidity are equivalent (one is always guaranteed to imply the other). Dotted arrows labeled with ∃ F α mean that the implicationis only valid for specific choices of F α and thus prestress. E (cid:48)(cid:48) (¯ x ) > N s = 1 means that the implication is guaranteed when there is only one state of selfstress. a. Energetic rigidity We define an energy functional E ( x ) = (cid:80) α f α ( x ) / (cid:80) α ( F α ( x ) − F α ) /
2; this energy functional implicitlydepends on our choice of F α . Any stiffnesses can be absorbed into the definitions of F α and F α . We say that a systemis energetically rigid at ¯ x if there exists a c such that E (¯ x + (cid:15)δx ) > E (¯ x ) for any nontrivial deformation δx and any < (cid:15) < c . In other words, it is energetically rigid if all sufficiently small, finite deformations increase the energy.Similarly, a system is energetically rigid at n th order at the configuration ¯ x if (cid:80) i ··· i n ∂ i · · · ∂ i n E (¯ x ) δx i · · · δx i n > for any nontrivial deformation, δx .With the definitions above, we can start relating energetic rigidity to other notions of rigidity. These relationshipsare summarized in Fig. 7. Important to note is that the dashed arrows signify that while the implication can beproved for some choice of self stress, it is not guaranteed that a given system has picked that particular self stress atthe energy minimum (i.e. the actual prestress may be a different linear combination of self stresses). The numberslabeling the propositions below refer to the arrows in Fig. 7 labeled with the same number. Proposition: (1) Energetic rigidity at ¯ x with E (¯ x ) > x is a critical point of the energy.Let ¯ x be a point that is energetically rigid. This means that E (¯ x + (cid:15)δx ) > E (¯ x ) for all nontrivial δx and for all0 < (cid:15) < c . Taking the derivative with respect to (cid:15) giveslim (cid:15) → ∂ (cid:15) E (¯ x + (cid:15)δx ) = (cid:88) n ∂ n E (¯ x ) δx n . (A1)If this were not a critical point then taking δx → − δx would give us a nontrivial deformation that decreases theenergy for some (cid:15) that was small enough. Therefore, it must be a critical point. Proposition: (2) The point ¯ x is a critical point of some energy with E (¯ x ) > x . Theconverse is also true for specific choices of F α .4We first assume ¯ x is a critical point with E (¯ x ) >
0. Then ∂ n E (¯ x ) = 0, which requires0 = (cid:88) α [ F α (¯ x ) − F α ] ∂ n F α (¯ x ) . (A2)Since E (¯ x ) (cid:54) = 0, F α (¯ x ) (cid:54) = F α . Therefore, F α (¯ x ) − F α is a self stress.Now assume that we have a point ¯ x where σ α is a self stress. Then choose F α = F α (¯ x ) + cσ α . We can now verifythat ¯ x is a critical point of E ( x ) = (cid:80) α [ F α ( x ) − F α (¯ x ) + cσ α ] for any c . Proposition: (3) The configuration ¯ x is energetically rigid at E ( x ) with E (¯ x ) = 0 if and only if ¯ x is rigid.We first assume that ¯ x is rigid. Then let F α = F α (¯ x ). We get E (¯ x ) = 0. Let δx be any nontrivial deformation. Since F α (¯ x + cδx ) (cid:54) = F α for sufficiently small c we must have E (¯ x + cδx ) > x is energetically rigid with E (¯ x ) = 0. Then F α (¯ x ) = F α . Since E (¯ x + cδx ) > c , we must have F α (¯ x + cδu ) (cid:54) = F α . Proposition: (4) Let ¯ x be an extremum of E ( x ) such that E (¯ x ) (cid:54) = 0 and suppose that ¯ x is energetically rigid. Thenthe system is structurally rigid at ¯ x as well.Suppose that ¯ x is an extremum of E ( x ) such that E (¯ x ) (cid:54) = 0 but such that ¯ x is energetically rigid. That is, allnontrivial directions raise the energy further. Then there cannot be any nontrivial isometries x ( t ) passing through ¯ x since if there were E would have to be constant along them and this contradicts the assumption.Note that this can be extended to energy maxima as well. The converse need not be true though. If a system isrigid at ¯ x , choosing F α so that ¯ x is an extremum does not mean that it will be energetically rigid. Let’s suppose that x ( t ) is a one-parameter family of constant energy trajectories. Then ∂ t E [ x ( t )] = 0 = (cid:88) α (cid:88) n [ F α ( x ( t )) − F α ] ∂ n F ( x ( t )) ˙ x n . (A3)This can only be true if x ( t ) are all extrema of E with E ( x ( t )) (cid:54) = 0. In addition, there must be at least one self stressalong the entire trajectory x ( t ).The notion of prestress stability is intimately related to energetic rigidity at quadratic order. The next propositionestablishes the equivalence of prestress stability (as defined above) and energetic rigidity to quadratic order: Proposition: (5) A system is prestress stable at ¯ x if and only if there is a choice F α such that it is an extremum ofthe energy with E (¯ x ) (cid:54) = 0 and is energetically rigid at quadratic order.To prove this we first assume that the system is prestress stable and let σ α be the self stress such that (cid:80) α σ α ∂ n ∂ m F α (¯ x ) is positive definite on nontrivial first-order flexes. Then define an energy functional E ( x ) = (cid:88) α [ F α ( x ) − F α (¯ x ) + cσ α ] , (A4)where c > x is an extremum, ∂ n E (¯ x ) = c (cid:80) α σ α ∂ n F α (¯ x ) = 0.Computing the Hessian, we find H nm = (cid:88) α ∂ n F α (¯ x ) ∂ m F α (¯ x ) + c (cid:88) α σ α ∂ n ∂ m F α (¯ x ) . (A5)This is positive definite on nontrivial first-order flexes by the assumption of prestress stability, for any c . On modesthat are not nontrivial first-order flexes, we can always choose c > c to be smaller than the smallest eigenvalue of the Gram term). Therefore, ¯ x is an energetically stableextremum of E ( x ) when F α = f α (¯ x ) − cσ α .Going the other way, let’s assume that our system is energetically rigid at quadratic order at an extremum ¯ x . Thenlet ˙ x n be any nontrivial, first-order flex. We have (cid:88) nm H nm ˙ x n ˙ x m = (cid:88) nm (cid:88) α [ F α (¯ x ) − F α ] ∂ n F (¯ x ) ˙ x n ˙ x m > . (A6)That implies that F α (¯ x ) − F α is a self stress and that it is prestress stable.It is worth noting that prestress stability at ¯ x does not mean that a system is energetically rigid at ¯ x for a particularchoice of F α , only for some choice.However, we also have the following two examples of showing to what extent these arrows are reversible.51. A system that is second-order rigid is not necessarily prestress stable. Examples appear in Connelly and Whitely.However, Proposition:
A system that is second-order rigid but has one self stress is prestress stable. This is also in [4].We must have cσ α ∂ n ∂ m f (¯ x ) positive definite for some, potentially negative, c . Then choosing F α = F α (¯ x ) − cσ α is energetically rigid to quadratic order and, hence, prestress stable.2. A system that is prestress stable may not be energetically rigid for a particular choice of F α . Suppose that asystem is prestress stable but has a self stress σ α for which the prestress matrix is not positive definite on thenontrivial first-order flexes. Choose F α = F α (¯ x ) − cσ α . This shows that the system with this choice is notenergetically rigid at quadratic order. In other words, the prestress that the system picks at ¯ x may not be onethat makes the system prestress stable.Finally, the following proposition deals with the nonlinear nature of rigidity: Proposition:
A system is energetically rigid at ¯ x with E (¯ x ) = 0 to fourth order if it is second-order rigid.This proposition shows that even if the standard checks of energetic rigidity (e.g. shear modulus) suggest floppiness,the system may still be energetically rigid to finite deformations. We will prove this proposition in the following section,where we also show a more detailed derivation of the equations in section II.
2. Second-order rigidity and energetic rigidity
Our goal here is to derive conditions for second-order zero modes and study the energy of systems that are second-order rigid. We will show that a system that has no prestress (Case 2A) but is second-order rigid is energetically rigidas well at fourth order in deformations. For prestressed systems, we show derivations of our claims for Case 2B and2C.Take constraints f α on a given system, e.g., f α ( { x n } ) may be the displacements of edges of a graph from theirequilibrium lengths. The energy functional is E = k (cid:80) Mα =1 f α / M is the number of constraints. We set k = 1without loss of generality. For a more general case with constraint dependent stiffnesses k α , we can simply rescale theconstraints to f (cid:48) α = √ k α f α . Imagine that ¯ x n is at a critical point of E .At a critical point, (cid:80) α f α ( { ¯ x n } ) ∂ m f α ( { ¯ x n } ) = 0. Let { σ α, , · · · , σ α,N s , e α, , · · · , e α,M − N s } be an orthonormal basisin R M where (cid:80) α σ α,I · ∂ n f α ( { ¯ x n } ) = 0 (so σ α,I are self stresses). Let us further assume f α ( { ¯ x n } ) = Cσ α, with C >
0, which we can do without loss of any generality.To find zero modes, we Taylor expand f α for small perturbations around ¯ x n . To easily keep track of the order ofexpansion, we parametrize deformations in time so that at an infinitesimal time δt we have a deformation x n ( δt ) suchthat x n (0) = ¯ x n . We then have f α ( { x n ( δt ) } ) ≈ Cσ α, + (cid:88) n ∂ n f α ˙ x n δt + 12 (cid:34)(cid:88) n ∂ n f α ¨ x n + (cid:88) nm ∂ n ∂ m f α ˙ x n ˙ x m (cid:35) δt + O ( δt ) , (A7)where partial derivatives are evaluated at ¯ x n . Also, ˙ x n is short hand for ˙ x n (0) and ¨ x n is short hand for ¨ x n (0). Thatis, these are explicitly independent vectors that determine the first two terms in a Taylor expansion of x n ( t ) around t = 0.It is useful to project f α ( { x n ( δt ) } ) along the orthonormal basis vectors (cid:88) α σ α,I f α ( { x n ( δt ) } ) ≈ Cδ I + (cid:88) α (cid:88) nm σ α,I ∂ n ∂ m f α ˙ x n ˙ x m δt , (A8) (cid:88) α e α,I f α ( { x n ( δt ) } ) ≈ (cid:88) α e α,I (cid:88) n ∂ n f α ˙ x n δt + 12 (cid:88) α e α,I (cid:34)(cid:88) n ∂ n f α ¨ x n + (cid:88) nm ∂ n ∂ m f α ˙ x n ˙ x m (cid:35) δt . (A9)6To find second-order zero modes, modes that preserve f α to second-order, Eqs. (A8-A9) imply the system (cid:88) α e α,I (cid:88) n ∂ n f α ˙ x n = 0 (cid:88) α e α,I (cid:34)(cid:88) n ∂ n f α ¨ x n + (cid:88) nm ∂ n ∂ m f α ˙ x n ˙ x m (cid:35) = 0 (cid:88) α (cid:88) nm σ α,I ∂ n ∂ m f α ˙ x n ˙ x m = 0where the first equation implies ˙ x n is along a linear zero mode (note that (cid:80) n ∂ n f α ˙ x n must have a nonzero projectionon at least one e α,I since it is perpendicular to all self stresses σ α,I by definition), the middle equation is associatedto the curvature of the linear zero mode as we proceed along t , and the last equation gives an additional quadraticconstraint that these tangents must satisfy to be second-order zero modes. Multiplying the last equation by δt , werecover Eq. 10.Notice that the middle equation always has a solution. To see this, we note that it is a linear equation of the form A ¨ x − b = 0. Since b is explicitly in the image of A , ¨ x has a solution that is unique up to zero modes. Since the linearzero modes are already included in ˙ x n , we can choose ¨ x n to be orthogonal to them without loss of generality. Withtaht choice, the matrix (cid:80) α e α,I ∂ n f α is invertible.Putting all of this into the energy, we find that E ≈ M − N s (cid:88) I =1 (cid:34)(cid:88) α e α,I (cid:88) n ∂ n f α ˙ x n + 12 (cid:88) α e α,I (cid:34)(cid:88) n ∂ n f α ¨ x n + (cid:88) nm ∂ n ∂ m f α ˙ x n ˙ x m (cid:35) δt (cid:35) δt (A10)+ 12 (cid:34) C + 12 (cid:88) α (cid:88) nm σ α, ∂ n ∂ m f α ˙ x n ˙ x m δt (cid:35) + 18 N s (cid:88) I =2 (cid:34)(cid:88) α (cid:88) nm σ α,I ∂ n ∂ m f α ˙ x n ˙ x m (cid:35) δt . What we are interested in is whether we can find a solution x n ( t ) such that E ( t ) increases, decreases, or staysconstant to a particular order in δt .Let us consider what happens when C → C = 0 because of the way they are prepared. Here, we assume that the energy can be continuously modulated tozero. Such a system is not prestressed, but can still possess self stresses (e.g. the onset of geometric incompatibility[29]). In that case, E ≈ M − N s (cid:88) I =1 (cid:34)(cid:88) α e α,I (cid:88) n ∂ n f α ˙ x n + 12 (cid:88) α e α,I (cid:34)(cid:88) n ∂ n f α ¨ x n + (cid:88) nm ∂ n ∂ m f α ˙ x n ˙ x m (cid:35) δt (cid:35) δt (A11)+ 18 N s (cid:88) I =1 (cid:34)(cid:88) α (cid:88) nm σ α,I ∂ n ∂ m f α ˙ x n ˙ x m (cid:35) δt . The energy is constant as long as the coefficients of δt , δt , and so on vanish. These lead to (cid:88) α e α,I (cid:88) n ∂ n f α ˙ x n = 0 , (A12)to second order, and we have the two equations (cid:88) α e α,I (cid:34)(cid:88) n ∂ n f α ¨ x n + (cid:88) nm ∂ n ∂ m f α ˙ x n ˙ x m (cid:35) = 0 , (A13)and (cid:88) α (cid:88) nm σ α,I ∂ n ∂ m f α ˙ x n ˙ x m = 0 , (A14)to fourth order. The third order term already vanishes if the quadratic term vanishes. These are the three equationsthat defined a quadratic isometry previously. Hence, E is constant along any quadratic isometry. Similarly, if E is7constant along a direction, the trajectory must be along a quadratic isometry. So at the critical point, second-orderrigidity implies energetic rigidity to this order in δt . This also proves the last proposition in the previous section.Now, one might wonder what happens as C increases. We then have E = C δt M − N s (cid:88) I =1 (cid:32)(cid:88) α e α,I (cid:88) n ∂ n f α ˙ x n (cid:33) + (cid:88) α (cid:88) nm Cσ α, ∂ n ∂ m f α ˙ x n ˙ x m + 12 δt M − N s (cid:88) I =1 (cid:32)(cid:88) α e α,I (cid:88) n ∂ n f α ˙ x n (cid:33) (cid:32)(cid:88) α (cid:48) e α (cid:48) ,I (cid:34)(cid:88) n ∂ n f α (cid:48) ¨ x n + (cid:88) nm ∂ n ∂ m f α (cid:48) ˙ x n ˙ x m (cid:35)(cid:33) (A15)+ 18 δt M − N s (cid:88) I =1 (cid:32)(cid:88) α e α,I (cid:34)(cid:88) n ∂ n f α ¨ x n + (cid:88) nm ∂ n ∂ m f α ˙ x n ˙ x m (cid:35)(cid:33) + 18 δt N s (cid:88) I =1 (cid:34)(cid:88) α (cid:88) nm σ α,I ∂ n ∂ m f α ˙ x n ˙ x m (cid:35) . The second-order term is the Hessian. If that has a direction that is negative, then we have not expanded arounda local minimum. However, one can ask whether or not zero directions might arise even if the system is second-orderrigid. For that to happen, however, ˙ x n cannot be along a zero mode. If it was along a zero mode and the Hessianwas zero, the fact that the system is second-order rigid would imply that the energy increases to fourth order. If ˙ x n was not along a zero mode and the Hessian was zero, for it to not increase the energy to the fourth order, it has tosatisfy Eq. A14, similar to second-order zero modes (this system would belong to Case 2C). Appendix B: Analytical calculations for springnetworks
Here we provide the details of our spring network cal-culations discussed in Section III A 1.It is useful to express f α explicitly in terms of theDOFs, i.e. the vertex positions X n (note that here weare using a vectorial notation so n ∈ { , . . . , N } ). To doso, we define the incidence matrix of the network, I αn I αn = , α leaves vertex n − , α enters vertex n . otherwise (B1)With this definition, we can rewrite the spring length as L α = (cid:12)(cid:12)(cid:12) (cid:88) n I αn X n (cid:12)(cid:12)(cid:12) = (cid:34)(cid:88) n,m I αn I αm X n · X m (cid:35) / . (B2)and the constraints as f α = (cid:34)(cid:88) n,m I αn I αm X n · X m (cid:35) / − L (B3)Now we study behavior of the zero modes. To do that,we first perturb the network X n → X n + δ X n . Taylorexpanding Eq. B3 to second order in δ X n , we find δf α = ( (cid:80) n I αn X n ) · ( (cid:80) m I αm δ X m ) L α +( (cid:80) n I αn δ X n ) L α − [( (cid:80) n I αn X n ) · ( (cid:80) m I αm δ X m )] L α . (B4) By comparing with Eq. 9, we determine the rigidity ma-trix of the system is defined as R αm = (cid:80) n I αn I αm X n L α . (B5)We can then simplify Eq. B4 δf α = (cid:88) m R αm · δ X m + ( (cid:80) n I αn δ X n ) − ( (cid:80) m R αm · δ X m ) L α . (B6)Linear zero modes are the solutions to (cid:80) m R αm · δ X (0) m = 0. If a linear zero mode is also a second-orderzero mode, it must additionally satisfy (Eq. 10) (cid:88) α σ α,I (cid:16)(cid:80) n I αn δ X (0) n (cid:17) L α = 0 . (B7)( (cid:80) n I αn δ X n ) > σ α,I >
0) exists, no zeromode will satisfy Eq. 15. The G = 0 condition is (Eq. 5) (cid:88) α f α ( (cid:80) n I αn δ X n ) − ( (cid:80) m R αm · δ X m ) L α = − (cid:88) α (cid:32)(cid:88) m R αm · δ X m (cid:33) (B8)for any mode δ X n . To retrieve Eqs. 15 and 16, we define L α = (cid:80) n I αn X n to be the vector along the edge α , and δ L α = (cid:80) n I αn δ X n to be its change due to perturbation δ X n . The component of δ L α parallel to L α is δL (cid:107) α = δL α = (cid:80) m R αm · δ X m + O ( δ X n ).8 Appendix C: Analytic calculations for Vertex models1. Second-order rigidity of vertex model with K A = 0 Here, we look the equations governing second-orderzero modes and the G = 0 condition for the vertex modelwith no area term ( K A = 0) discussed in section III B 1,thus M = N cell . The constraints on the vertices are givenby f α = P α − P = (cid:88) edge j ∈ cell α L j − P = (cid:88) edge j ∈ cell α (cid:12)(cid:12)(cid:12)(cid:12) (cid:88) n I jn X n (cid:12)(cid:12)(cid:12)(cid:12) − P . (C1)If we define a cell-edge adjacency matrix A αj by: A αj = (cid:26) , edge j ∈ cell α , otherwise (C2)it allows us to rewrite f α to show the dependence on α more explicitly: f α = (cid:88) j A αj (cid:12)(cid:12)(cid:12)(cid:12) (cid:88) n I jn X n (cid:12)(cid:12)(cid:12)(cid:12) − P , (C3)where j now runs through all the edges and n throughall vertices.Now, we perturb vertex positions X n → X n + δ X n .We get for δf α an expression that is similar to Eq. B4: δf α = (cid:88) m R αm · δ X m + 12 (cid:88) j A αj (cid:20) ( (cid:80) n I jn δ X n ) L j − ( (cid:80) nm I jn I jm X n · δ X m ) L j (cid:21) . (C4)Where we have defined the rigidity matrix as R αm = (cid:88) j (cid:88) n A αj I jn I jm X n L j . (C5)Note that the last term in Eq. C4 cannot be written interms of R αm anymore.The self stresses of R αm impose the following quadraticconstraints on zero modes δ X n : (cid:88) α σ α,I (cid:88) j A αj (cid:20) (cid:16)(cid:80) n I jn δ X (0) n (cid:17) L j − (cid:16)(cid:80) nm I jn I jm X n · δ X (0) m (cid:17) L j (cid:21) = 0 , (C6)which simplifies to (cid:88) α σ α,I (cid:88) j A αj (cid:0) δ L ⊥ j (cid:1) L j = 0 . (C7) (a)(b) FIG. 8.
Vertex Model ( K A = 0 ) Density of states D asa function of the eigenvalue λ for the Hessian matrix (solidblack line), Gram matrix (dashed yellow line) and prestressmatrix (dash-dotted red line), averaged over 10 samples with N cell = 500 and and K P = 1. (a) In the Floppy regime( P = 3 . N cell zero modes and their DOS curves overlap as there is noprestress in the network. (b) In the rigid regime ( P = 3 . N cell + 1 zero modes and N s = 1,while the Hessian only possesses two trivial (translational)zero modes. The Hessian DOS is dominated by the Gramterm at high frequency and by the prestress matrix at lowfrequencies. Since all vertices are connected to three edges and thebox size is fixed, there are no nontrivial motions that donot introduce a δ L ⊥ j and thus the inner sum is positivedefinite. Thus, similar to spring networks, if a self stress σ α,I > G = 0condition (Eq. (5)): (cid:88) α f α (cid:88) j A αj (cid:0) δ L ⊥ j (cid:1) L j = − (cid:88) α (cid:32)(cid:88) m R αm · δ X m (cid:33) , (C8)which cannot be satisfied for any nontrivial zero mode if f α >
0. In Fig. 8 and 9 we show plots for vertex modelsimulations with K A = 0 similar to Fig. 4 and 5.
2. Prestresses in the floppy regime of vertex modelwith K A (cid:54) = 0 It is numerically possible for vertex model configura-tions in the P > P ∗ regime to be prestressed locally.9 (a)(b)(c)(d) FIG. 9.
Vertex Model ( K A = 0 ) Comparison of Densityof States as a function of the eigenvalue λ for the Hessianmatrix ( D H ), Gram matrix ( D G ) and prestress matrix ( D P ),averaged over 10 samples with N cell = 500 for different valuesof P . (a) Shear modulus ( G ) with as a function of P shows atransition from fluid to solid at P ∗ ≈ .
78. The error bars rep-resent standard deviation over 10 samples. (b-d) Comparisonbetween DOS of Hessian and its components for different val-ues of P . (b) Transitioning from floppy to rigid coincides withthe appearance of low frequency eigenmodes in the HessianDOS that shift to higher frequencies as the system becomesmore rigid. (c) Even though in the rigid regime there is onlyone new zero mode, the Gram matrix DOS looks significantlydifferent when compared to systems with P < P ∗ becausethe geometry of the system does not change significantly inthe rigid regime. (d) Unlike the Gram matrix, the prestressincreases linearly as L is decreased in the rigid regime, asreflected by the shift in the prestress matrix DOS towardshigher frequencies. Unlike the vertex model with K A (cid:54) = 0,the prestress matrix has no negative eigenvalues. This phenomenon has been reported before [29]. Like-wise, we have encountered some cases with four-sidedpolygons that were prestressed at P = 3 .
90. This is be-cause those four-sided polygons could not achieve boththeir preferred area and perimeters, A and P , even witha zero shear modulus as the prestress is localized. Fig. 4and 5 in the main text excludes such cases. Appendix D: Numerical methods1. Structure initialization for spring networks andvertex model
For both spring networks and vertex model, we use cell-GPU [59] to initialize N cell cell centers randomly in a pe-riodic box of size L x L y = N cell with L x = L y = √ N cell .A Voronoi tessellation is applied to get N cell polygon cellswith 2 N cell vertices with coordination number z = 3.The final step in the initialization process involves mov-ing the cell centers for a few time steps using a self-propelled Voronoi model [41] to make cell areas more uni-form. After the initialization process, the energy (Eq. 14for spring networks and Eq. 17 for vertex model) is min-imized by moving the vertices using the FIRE minimizer[60] with a force cutoff of 10 − . For vertex model, a T t = 0 . t max = 0 .
2. Structure initialization for jammed packings
We create 2D disk packings using a quad-precisionGPU implementation of the FIRE algorithm [61, 62].First, N particles with a polydispersity of 20% are ran-domly distributed in a periodic box of size L = 1 andthen radii are re-scaled to a packing density well abovejamming transition which is typically φ J (cid:39) .
84 for 2Dsystems. Finally, the system is minimized to its inherentstructure using the FIRE minimizer with a force cut-off of 10 − . At densities far from jamming, a packingwill have many states of self stress. To bring the systemto the critical jamming with one state of self stress, wesuccessively re-scale the density to smaller values and re-minimize the energy. Once the system reaches one stateof self stress, the pressure will be in order of p ≈ − and the initialization process is halted.
3. Density of states and shear modulus
To find the eigenvalue spectrum of the Hessian, Gramterm and prestress matrix, we calculate the Hessian ma-trix and its components for a given system at an energyminimum and consider eigenvalues with an absolute value0smaller than 10 − as zero eigenmodes. For the density ofstates plot, we sort all the eigenvalues and use equiprob-able (Dirichlet) binning with 150 bins such that there isan equal number of eigenvalues in each bin, from whichwe can plot a normalized histogram representing the den-sity of states. We also apply a centered moving averagewith window size 3 to smooth the curves. Zero modeswould represent a peak at 0 which are not plotted. Forthe shear modulus plots, we modify the periodic bound-ary conditions to accommodate a skew (i.e. Lees-Edwardsboundary conditions) with a simple shear parameter γ ,which allows us to write the energy function as a func- tion of γ [43]. We then use Eq. (2) to calculate the shearmodulus.
4. Numerical results for vertex models with noarea constraints ( K A = 0 ) Using the same analyses as discussed for Figures 4 and5, we calculate the contributions to the Hessian densityof vibrational states for vertex models where there areno area constraints, ( K A = 0). These data are shown inFig 8 and 9, respectively. [1] T. G. Abbott, Generalizations of Kempe’s universal-ity theorem , Master’s thesis, Massachusetts Institute ofTechnology (2008).[2] C. Calladine, Buckminster fuller’s “tensegrity” structuresand clerk maxwell’s rules for the construction of stiffframes, International Journal of Solids and Structures ,161 (1978).[3] C. R. Calladine and S. Pellegrino, First-order infinites-imal mechanisms, International Journal of Solids andStructures , 10.1016/0020-7683(91)90137-5 (1991).[4] R. Connelly and W. Whiteley, Second-order rigid-ity and prestress stability for tensegrity frame-works, SIAM Journal on Discrete Mathematics ,10.1137/S0895480192229236 (1996).[5] R. Connelly and S. Guest, Frameworks, tensegritiesand symmetry: understanding stable structures (CornellUniv., College Arts Sci., 2015).[6] J. C. Maxwell, On the calculation of the equilibrium andstiffness of frames, The London, Edinburgh, and DublinPhilosophical Magazine and Journal of Science , 294(1864), https://doi.org/10.1080/14786446408643668.[7] F. Lechenault, O. Dauchot, G. Biroli, and J. P.Bouchaud, Critical scaling and heterogeneous superdiffu-sion across the jamming/rigidity transition of a granularglass, EPL , 10.1209/0295-5075/83/46003 (2008).[8] M. V. Hecke, Jamming of soft particles: Geometry, me-chanics, scaling and isostaticity, Journal of Physics Con-densed Matter , 033101 (2009).[9] W. G. Ellenbroek, V. F. Hagh, A. Kumar, M. F.Thorpe, and M. V. Hecke, Rigidity loss in disorderedsystems: Three scenarios, Physical Review Letters ,10.1103/PhysRevLett.114.135501 (2015).[10] D. B. Liarte, X. Mao, O. Stenull, and T. C. Lubensky,Jamming as a multicritical point, Physical Review Let-ters , 10.1103/PhysRevLett.122.128006 (2019).[11] T. C. Lubensky, C. L. Kane, X. Mao, A. Souslov,and K. Sun, Phonons and elasticity in critically coordi-nated lattices, Reports on Progress in Physics , 073901(2015).[12] D. Z. Rocklin, B. G.-g. Chen, M. Falk, V. Vitelli, andT. C. Lubensky, Mechanical weyl modes in topologicalmaxwell lattices, Phys. Rev. Lett. , 135503 (2016).[13] K. Bertoldi, V. Vitelli, J. Christensen, and M. Hecke,Flexible mechanical metamaterials, Nature Reviews Ma-terials , 17066 (2017). [14] X. Mao and T. C. Lubensky, Maxwell lattices and topo-logical mechanics, Annual Review of Condensed MatterPhysics , 413 (2018), https://doi.org/10.1146/annurev-conmatphys-033117-054235.[15] D. Zhou, L. Zhang, and X. Mao, Topological edge floppymodes in disordered fiber networks, Phys. Rev. Lett. ,068003 (2018).[16] J. Xu, Y. Tseng, and D. Wirtz, Strain hardening of actinfilament networks: Regulation by the dynamic cross-linking protein α -actinin, Journal of Biological Chemistry , 10.1074/jbc.M002377200 (2000).[17] S. Rammensee, P. A. Janmey, and A. R. Bausch, Me-chanical and structural properties of in vitro neuro-filament hydrogels, European Biophysics Journal ,10.1007/s00249-007-0141-7 (2007).[18] G. H. Koenderink, Z. Dogic, F. Nakamura, P. M.Bendix, F. C. MacKintosh, J. H. Hartwig, T. P.Stossel, and D. A. Weitz, An active biopolymer net-work controlled by molecular motors, Proceedings ofthe National Academy of Sciences , 1358 (2010), pMID: 20392048,https://doi.org/10.1021/bm100136y.[20] F. Burla, J. Tauber, S. Dussi, J. Gucht, and G. Koen-derink, Stress management in composite biopolymer net-works, Nature Physics , 549–553 (2019).[21] W. Tang and M. F. Thorpe, Percolation of elas-tic networks under tension, Physical Review B ,10.1103/PhysRevB.37.5539 (1988).[22] C. Storm, J. J. Pastore, F. C. MacKintosh, T. C. Luben-sky, and P. A. Janmey, Nonlinear elasticity in biologicalgels, Nature , 191 (2005).[23] M. Wyart, H. Liang, A. Kabla, and L. Mahadevan, Elas-ticity of floppy and stiff random networks, Phys. Rev.Lett. , 215501 (2008).[24] M. Sheinman, C. P. Broedersz, and F. C. MacKintosh,Nonlinear effective-medium theory of disordered springnetworks, Physical Review E - Statistical, Nonlinear, andSoft Matter Physics , 10.1103/PhysRevE.85.021801(2012).[25] J. Feng, H. Levine, X. Mao, and L. M. Sander, Nonlinearelasticity of disordered fiber networks, Soft Matter ,10.1039/c5sm01856k (2016). [26] A. Sharma, A. Licup, K. Jansen, R. Rens, M. Sheinman,G. Koenderink, and F. MacKintosh, Strain-controlledcriticality governs the nonlinear mechanics of fibre net-works, Nature Physics , 584 (2016).[27] M. F. J. Vermeulen, A. Bose, C. Storm, and W. G. Ellen-broek, Geometry and the onset of rigidity in a disorderednetwork, Phys. Rev. E , 053003 (2017).[28] R. Rens, C. Villarroel, G. D¨uring, and E. Lerner, Mi-cromechanical theory of strain stiffening of biopolymernetworks, Phys. Rev. E , 062411 (2018).[29] M. Merkel, K. Baumgarten, B. P. Tighe, and M. L.Manning, A minimal-length approach unifies rigid-ity in underconstrained materials, Proceedings of theNational Academy of Sciences ,10.1103/PhysRevLett.122.188003 (2019).[31] T. E. Angelini, E. Hannezo, X. Trepatc, M. Marquez,J. J. Fredberg, and D. A. Weitz, Glass-like dynamicsof collective cell migration, Proceedings of the NationalAcademy of Sciences of the United States of America , 10.1073/pnas.1010059108 (2011).[32] M. Sadati, N. T. Qazvini, R. Krishnan, C. Y. Park, andJ. J. Fredberg, Collective migration and cell jamming,Differentiation , 121 (2013).[33] K. E. Kasza, D. L. Farrell, and J. A. Zallen, Spa-tiotemporal control of epithelial remodeling by regu-lated myosin phosphorylation, Proceedings of the Na-tional Academy of Sciences of the United States of Amer-ica , 10.1073/pnas.1400520111 (2014).[34] S. Garcia, E. Hannezo, J. Elgeti, J. F. Joanny, P. Sil-berzan, and N. S. Gov, Physics of active jamming duringcollective cellular motion in a monolayer, Proceedings ofthe National Academy of Sciences of the United Statesof America , 10.1073/pnas.1510973112 (2015).[35] J.-A. Park, J. H. Kim, D. Bi, J. A. Mitchel, N. T.Qazvini, K. Tantisira, C. Y. Park, M. McGill, S.-H. Kim,B. Gweon, J. Notbohm, R. S. Jr, S. Burger, S. H. Ran-dell, A. T. Kho, D. T. Tambe, C. Hardin, S. A. Shore,E. Israel, D. A. Weitz, D. J. Tschumperlin, E. Henske,S. T. Weiss, M. L. Manning, J. P. Butler, J. M. Drazen,and J. J. Fredberg, Unjamming and cell shape in theasthmatic airway epithelium, Nature Materials , 1040(2015).[36] J. A. Park, L. Atia, J. A. Mitchel, J. J. Fredberg, andJ. P. Butler, Collective migration and cell jamming inasthma, cancer and development, Journal of Cell Science , 10.1242/jcs.187922 (2016).[37] X. Wang, M. Merkel, L. B. Sutter, G. Erdemci-Tandogan,M. L. Manning, and K. E. Kasza, Anisotropy links cellshapes to tissue flow during convergent extension, Pro-ceedings of the National Academy of Sciences of theUnited States of America , 10.1073/pnas.1916418117(2020).[38] R. Farhadifar, J. C. R¨oper, B. Aigouy, S. Eaton, andF. J¨ulicher, The influence of cell mechanics, cell-cell inter-actions, and proliferation on epithelial packing, CurrentBiology , 10.1016/j.cub.2007.11.049 (2007).[39] D. B. Staple, R. Farhadifar, J. C. R¨oper, B. Aigouy,S. Eaton, and F. J¨ulicher, Mechanics and remodellingof cell packings in epithelia, European Physical JournalE , 10.1140/epje/i2010-10677-0 (2010). [40] D. Bi, J. H. Lopez, J. M. Schwarz, and M. L. Manning, Adensity-independent rigidity transition in biological tis-sues, Nature Physics , 10.1038/nphys3471 (2015).[41] D. Bi, X. Yang, M. C. Marchetti, and M. L. Man-ning, Motility-driven glass and jamming transitions inbiological tissues, Physical Review X , 10.1103/Phys-RevX.6.021011 (2016).[42] M. Moshe, M. J. Bowick, and M. C. Marchetti, Geo-metric frustration and solid-solid transitions in model2d tissue, Physical Review Letters , 10.1103/Phys-RevLett.120.268105 (2018).[43] M. Merkel and M. L. Manning, A geometrically con-trolled rigidity transition in a model for confluent3d tissues, New Journal of Physics , 10.1088/1367-2630/aaaa13 (2018).[44] L. Yan and D. Bi, Multicellular rosettes drive fluid-solidtransition in epithelial tissues, Physical Review X ,10.1103/PhysRevX.9.011029 (2019).[45] S. Arzash, J. L. Shivers, and F. C. MacKintosh, Finitesize effects in critical fiber networks, Soft Matter ,10.1039/d0sm00764a (2020).[46] J. A. Jackson, M. C. Messner, N. A. Dudukovic, W. L.Smith, L. Bekker, B. Moran, A. M. Golobic, A. J. Pascall,E. B. Duoss, K. J. Loh, and C. M. Spadaccini, Fieldresponsive mechanical metamaterials, Science Advances , 10.1126/sciadv.aau6419 (2018).[47] C. Merrigan, C. Nisoli, and Y. Shokef, Topological mem-ory and hysteresis in ice-like mechanical metamaterials(2020), arXiv:2010.05227 [cond-mat.soft].[48] A. Donev, R. Connelly, F. H. Stillinger, and S. Torquato,Underconstrained jammed packings of nonspherical hardparticles: Ellipses and ellipsoids, Physical Review E- Statistical, Nonlinear, and Soft Matter Physics ,10.1103/PhysRevE.75.051304 (2007).[49] M. Mailman, C. F. Schreck, C. S. O’Hern, andB. Chakraborty, Jamming in systems composed of fric-tionless ellipse-shaped particles, Physical Review Letters , 10.1103/PhysRevLett.102.255501 (2009).[50] Z. Zeravcic, N. Xu, A. J. Liu, S. R. Nagel, and W. V.Saarloos, Excitations of ellipsoid packings near jamming,EPL , 10.1209/0295-5075/87/26001 (2009).[51] A. Bose, M. F. Vermeulen, C. Storm, and W. G. Ellen-broek, Self-stresses control stiffness and stability in over-constrained disordered networks, Physical Review E ,10.1103/PhysRevE.99.023001 (2019).[52] B. G. G. Chen and C. D. Santangelo, Branches of trian-gulated origami near the unfolded state, Physical ReviewX , 10.1103/PhysRevX.8.011034 (2018).[53] S. Guest, The stiffness of prestressed frameworks: A uni-fying approach, International Journal of Solids and Struc-tures , 10.1016/j.ijsolstr.2005.03.008 (2006).[54] M. Holmes-Cerfon, L. Theran, and S. J. Gortler,Almost-rigidity of frameworks (2020), arXiv:1908.03802[math.MG].[55] J. D. Sartor and E. I. Corwin, Direct measurement offorce configurational entropy in jamming, Physical Re-view E , 10.1103/PhysRevE.101.050902 (2020).[56] R. Connelly and H. Servatius, Higher-order rigid-ity—what is the proper definition?, Discrete and Com-putational Geometry , 193 (1994).[57] X. Yang, D. Bi, M. Czajkowski, M. Merkel, M. L. Man-ning, and M. C. Marchetti, Correlating cell shape andcellular stress in motile confluent tissues, Proceedings ofthe National Academy of Sciences of the United States of America , 10.1073/pnas.1705921114 (2017).[58] E. Degiuli, A. Laversanne-Finot, G. D¨uring, E. Lerner,and M. Wyart, Effects of coordination and pressureon sound attenuation, boson peak and elasticity inamorphous solids, Soft Matter , 10.1039/c4sm00561a(2014).[59] D. M. Sussman, cellgpu: Massively parallel simulationsof dynamic vertex models, Computer Physics Communi-cations , 10.1016/j.cpc.2017.06.001 (2017). [60] E. Bitzek, P. Koskinen, F. G¨ahler, M. Moseler, andP. Gumbsch, Structural relaxation made simple, Phys-ical Review Letters , 10.1103/PhysRevLett.97.170201(2006).[61] P. Charbonneau, E. I. Corwin, G. Parisi, and F. Zam-poni, Universal microstructure and mechanical stabil-ity of jammed packings, Physical Review Letters ,10.1103/PhysRevLett.109.205501 (2012).[62] P. K. Morse and E. I. Corwin, Geometric signatures ofjamming in the mechanical vacuum, Physical Review Let-ters112