Rigidity and auxeticity transitions in networks with strong bond-bending interactions

A widely-studied model for gels or biopolymeric fibrous materials are networks with central force interactions, such as Hookean springs. Less commonly studied are materials whose mechanics are dominated by non-central force interactions such as bond-bending potentials. Inspired by recent experimental advancements in designing colloidal gels with tunable interactions, we study the micro- and macroscopic elasticity of two-dimensional planar graphs with strong bond bending potentials, in addition to weak central forces. We introduce a theoretical framework that allows us to directly investigate the limit in which the ratio of characteristic central-force to bending stiffnesses vanishes. In this limit we show that a generic isostatic point exists at $z_c=4$, coinciding with the isostatic point of frames with central force interactions in two dimensions. We further demonstrate the emergence of a stiffening transition when the coordination is increased towards the isostatic point, which shares similarities with the strain-induced stiffening transition observed in biopolymeric fibrous materials, and coincides with an auxeticity transition above which the material's Poisson's ratio approaches -1 when bond-bending interactions dominate.

A widely-studied model for gels or biopolymeric fibrous materials are networks with central force interactions, such as Hookean springs. Less commonly studied are materials whose mechanics are dominated by non-central force interactions such as bond-bending potentials. Inspired by recent experimental advancements in designing colloidal gels with tunable interactions, we study the microand macroscopic elasticity of two-dimensional planar graphs with strong bond bending potentials, in addition to weak central forces. We introduce a theoretical framework that allows us to directly investigate the limit in which the ratio of characteristic central-force to bending stiffnesses vanishes. In this limit we show that a generic isostatic point exists at zc = 4, coinciding with the isostatic point of frames with central force interactions in two dimensions. We further demonstrate the emergence of a stiffening transition when the coordination is increased towards the isostatic point, which shares similarities with the strain-induced stiffening transition observed in biopolymeric fibrous materials, and coincides with a auxeticity transition above which the material's Poisson's ratio approaches -1 when bond-bending interactions dominate.

I. INTRODUCTION
In 1864 Maxwell spelled out a criterion that frames of freely-hinged struts need to satisfy in order to be mechanically stable [1]: if the average connectivity z is higher than a threshold value z c ≡ 2d ind spatial dimensions, rigidity of the frame is guaranteed, regardless of the exact way the elements are connected (as long as fluctuations in connectivity are limited [2]). In frames with z < z c collective modes exist that are floppy [3][4][5], which means that motion associated with these modes will respect the perfect rigidity of the struts. The gradual disappearance of such floppy modes as z → z c is known as the jamming transition [6][7][8], and has been related to various mechanical and dynamical phenomena such as the divergence of viscosity in non-Brownian suspensions [9] and the fragility of chalcogenide glass formers [10].
Most studies of jamming phenomena focus on the role of steric interactions, often modeled by some form of central forces e.g. Hookean springs [11] or hard-sphere repulsions [12,13]. In this work we explore the mechanical properties of a different class of materials: disordered networks in which the dominant interaction takes the form of bond bending. Our focus is motivated by recent advancements in the fabrication of colloidal gells with tunable interactions; in particular, we were inspired by the work of Schall et al. [14], who built and controlled nano-and micrometer size superstructures using critical Casimir forces on patchy colloidal particles. By measuring fluctuations of the constituent colloids, the authors of [14] established that the stiffness of bond bending in their superstructures is much larger than the stiffness of radial interactions, by up to two orders of magnitude or more [15].
Inspired by these outstanding experiments mentioned above, we set out to address the following questions: (i) what are the elastic properties of materials whose mechanics are dominated by bond-bending interactions, (ii) how should the micromechanics of this class of materials be understood from a geometric perspective, and (iii) what sort of jamming phenomenology emerges in this class of systems.
In this work we consider disordered networks in two dimensions (2D) of mean coordination z, and introduce both radial and bond bending interactions, characterized by stiffnesses k r and k θ , respectively. We consider the ratio µ ≡ k r /k θ between these stiffnesses as a key tunable parameter of the material (in addition to the coordination z), and investigate numerically and theoretically the behavior of elastic moduli under variations of z and µ. We show that as the limit µ → 0 is approached, scaling behavior of elastic moduli emerges, as a function of the distance between the mean connectivity and the system's jamming point, shown in what follows to coincide with the the Maxwell threshold [1] z c = 4 at which the generic, 'central-force' isostatic point occurs.
We take two complementary routes in order to study theoretically the limit µ → 0 at which the scaling behavior of elastic moduli emerges. First, we fix k θ and set k r = 0; this arrangement allows us to understand the scaling behavior of elastic moduli in the hyperstatic regime z > z c . We further put forward a scaling argument supported by numerical tests that the hyperstatic regime is characterized by a diverging length c ∼ (z − z c ) −1/2 . Even more intruiging is the limit µ → 0 obtained by fixing k r and sending k θ → ∞; this is achieved by considering the angles formed between bonds that share a common node to be entirely fixed, i.e. they serve as geometric constraints. In this limit, and away from the isostatic point, hypostatic networks feature elastic moduli ∼ k r . Interestly, we find that the shear modulus G diverges as z is made to approach z c from below, while the bulk modulus K remains regular. We present a theoretical framework that allows us to predict the divergence of G with z c − z in this limit, and find good agreement with our numerical calculations.
Finally, we consider the auxeticity of our model material; we find that far into the hypostatic regime the material features a positive Poisson's ratio, of around 0.3, a value characteristic to many disordered materials [16]. However, ν rapidly decreases as z is increased. Interestingly, in the limit µ → 0 the material approaches perfect auxeticity ν → −1 as z → z c , and remains perfectly auxetic in the entire hyperstatic regime.
Our work is structured as follows; in Sect. II we spell out the model ingredients and observables considered in our study. In Sect. III we present a numerical investigation of the elastic properties of our model, as a function of the two key control parameters, namely the ratio of stiffnesses µ ≡ k r /k θ and the mean coordination z. In Sect. IV we explain the occurrence of an isostatic point at z c = 4 in our model. In Sect. V we consider the hyperstatic regime, and study theoretically the limit µ → 0, while Sect. VI presents a theoretical framework that allows us to study the hypostatic regime in the limit µ → 0. In Sect. VII we provide scaling arguments and show numerically that a characteristic lengthscale diverges as the isostatic point is approached. Sect. VIII discusses the auxeticity of our model, and our work is summarized in Sect. IX, where we discuss future research directions.

II. MODEL AND KEY OBSERVABLES
We consider 2D disordered planar graphs (networks) with periodic boundary considitions, whose topology is characterized by the mean number of edges per node, denoted by z. We build these disordered networks by adopting the contact network of dense packings of soft spheres, and pruning edges according to a protocol that maintains low node-to-node fluctuations of connectivity. Our network generation protocol is described in Appendix A, and an example of a network with z = 3.95 is shown in Fig. 1a. We introduce the following potential energy U for our disordered networks where k θ and k r denote the bond-bending and centralforce stiffnesses, respectively, and¯ denotes the microscopic units of length. The first term on the RHS of Eq. (1) represents a sum over all angles θ ijk formed between pairs of edges that share a common node, with no other edges found in between the said pair, as illustrated in Fig. 1b. We define the deviations from the rest-angles ijk denoting the initial (ground state) rest-angles. The second term on the RHS of Eq. (1) represents a sum over the network's edges, and we define the deviation from the rest-lengths ∆r ij ≡ r ij − ij where ij denotes the rest-length of the edge connecting the i, j pair of nodes. Lengths are expressed in terms of¯ which denotes the mean rest-length. Below we will always assume that the networks reside at their respective ground states, i.e. all of the angles are equal to their rest-angles, and all edges reside at their rest-lengths, implying that U = 0 and that there are no stresses in the material.
In what follows we consider simple shear and expansive strains, described by the strain tensor = (γ, η) most conveniently parameterized by simple shear and expansive strain parameters γ and η, respectively. In 2D the strain tensor assumes the form Given a strain tensor that describes an imposed deformation mode, distances r ij = √ x ij · x ij between the coordinates of any two nodes x i , x j vary under such imposed deformations according to where We focus on macroscopic elastic properties as seen in the athermal shear and bulk moduli, denoted as G and K, respectively. The athermal shear modulus is defined as whereas the athermal bulk modulus is given by Total dervatives e.g. d/dγ are understood as taken in the athermal limit, i.e. under the constraints dictated by mechanical equilibrium [17]. In this representation we see that (i) a jamming transition occurs at zc = 4, (ii) that G ∼ k θ deep in the hyperstatic regime z > zc, and (iii) that G ∼ kr deep in the hypostatic regime z < zc. The black circles correspond to the limit µ → 0, discussed in detail separately for the hyperstatic regime in Sect. V and for the hypostatic regime in Sect. VI.

III. ELASTIC PROPERTIES OF BOND-BENDING-DOMINATED NETWORKS
We start the presentation of our results with a numerical investigation of the shear modulus variation under changes of the mean coordination z and the ratio µ of bond-bending to central-force stiffnesses. In Fig. 2 we plot the sample-to-sample means of the shear modulus averaged over 20 random networks of N = 25600 nodes. The left panel plots the ratio G/k θ vs. the coordination z; noticeably, z = 4 marks the onset of an underlying jamming transition, that becomes more pronounced as µ → 0. For small µ the ratio G/k θ grows by several orders of magnitude as z approaches the critical coordination z c = 4, in a fashion reminiscent to the strain-stiffening transition observed upon deformation of biopolymeric fibrous materials [18][19][20][21]. Far above z c the ratio G/k θ becomes roughly independent of µ, indicating that in this regime G ∼ k θ .
In the right panel of Fig. 2 we show the same sampleto-sample means of the shear modulus, this time rescaled by k r , to find that deep in the hypostatic regime z < 4, G ∼ k r . The data pertaining to µ = 0 in Fig. 2 will be discussed in detail in what follows.
In Fig. 3 we show a scaling plot for G; here G/k θ is rescaled by |δz| f and plotted against the ratio k r /|δz| φ where δz ≡ z − z c denotes the distance to the critical coordination z c ≡ 4. The best collapse is found using the exponents f = 1.25 and φ = 2. 25. In what follows we will argue that the mean field exponents are f = 1 and φ = 2, and discuss the discrepancy we find with the mean field predictions.
To understand the scaling behavior of the shear mod-ulus as seen in Figs. 2 and 3, in the next Sections we explain the occurrence of an isostatic point at z c = 4 in our model, and consider the limit µ → 0 separately in the hyperstatic z > z c and the hypostatic z < z c regimes.

IV. THE ANGLE-PRESERVING ISOSTATIC POINT OF 2D PLANAR NETWORKS
The emergence of an isostatic point becomes apparent when the limit µ → 0 is considered, as seen in Fig. 2. In the hypostatic regime z ≤ z c the limit µ → 0 can be obtained by fixing k r and sending k θ → ∞. Under these cir-cumstances the bond-bending interactions can be treated as geometric constraints; a displacement field u on the network's nodes will leave the angles θ ijk invariant (to leading order in u ) if it satisfies for every angle θ ijk , where here and in what follows repeated node indices are understood to be summed over. An example of the field ∂θ ijk ∂x is shown in Fig. 1b. It is useful to define the linear operator which takes vectors from the space of the nodes' coordinates to the space defined by the entire set of angles. The number of rows Q features is equal to the total number of angles N z, whereas the number of columns is 2N , each represented a spatial coordinate of a node. Nontrivial solutions to Eq. 6 are expected to exist if the the rank of the operator Q is smaller than the number of degrees of freedom 2N available to the displacement field u. We next argue that the rank of Q becomes exactly equal to 2N at z = 4; to this aim we define z i to be the number of edges connected to the i th node, and we note that the number of angles that surround each node is equal to z i . A displacement field that preserves z i −1 angles that surround a single node must also preserve the z th i angle as well, meaning that a single angle for each node can be selected, and the row of Q corresponding to that particular angle can be eliminated (it will be expressible as a linear combination of other rows). This amounts to eliminating N rows out of the N z rows of Q.
We finally note that each face of our planar network is a polygon with m edges, built from m angles; a displacement field that preserves m−1 angles of a face must also preserve the m th angle, further reducing the number of independent rows of Q by the number F of faces of the network. The latter is related to the number of nodes N and the number of edges 1 2 N z via Euler's formula for a planar graph embedded on a torus The total number of independent rows left after the eliminations described above is thus which becomes equal to the dimension of configuration space 2N when the coordination reaches z c = 4. In these considerations we neglect corrections of order N −1 that arise from collective translations or deformations of space [22]. We conclude that in our two-dimensional disordered networks nontrivial displacement fields that preserve the entire set of angles will exist if z < 4, whereas when z > 4 no such displacements exist. Interestingly, the critical coordination that separates these two regimes coincides with the Maxwell threshold z c = 2d of networks of rigid struts [1] discussed intensively in the jamming literature, see e.g. [7,8].
Having established that an underlying isostatic point exists at z = 4, we next turn to discussing the scaling behavior of elastic moduli as the critical coordination is approached from above and from below.

V. THE HYPERSTATIC REGIME z > zc
It is illuminating to study the elasticity of our model in the hyperstatic regime by fixing k θ = 1 and sending k r → 0, resulting in µ → 0 (see black circles in Fig. 2a). These circumstances correspond to eliminating the second sum on the RHS of Eq. (1) for the potential energy, leaving us with It is illuminating to study the elasticity of our model in the hyperstatic regime by fixing k θ = 1 and sending k r → 0, resulting in µ → 0 (see black circles in Fig. 2a). These circumstances correspond to eliminating the second sum on the RHS of Eq. (1) for the potential energy, leaving us with In the remainder of this Section we show how in this limit elastic moduli expressions can be simply expressed via the geometry of the network, and provide scaling arguments to moduli's dependence on coordination for z > z c .

A. Elastic moduli in the hyperstatic regime
The microscopic expression for the athermal shear modulus G for systems governed by a potential energy U reads [17] where is known as the dynamical matrix. Recall that our systems are assumed to have U = 0 (i.e. all of the angles reside exactly at their rest-angles) before any deformation is imposed, leading to a simple form for the dynamical matrix; it reads Working in units such that¯ = 1 and k θ = 1, employing Eq. 7 and adopting bra-ket and matrix notations, the dynamical matrix takes the form where here and in what follows we assume that the redundent rows of Q have been eliminated as described in Sect. IV. We next see that under the circumstances of fixing k θ = 1 and sending k r to zero, which can be expressed using our bra-ket notation as where we denoted ∂ γ ≡ ∂/∂γ and ∂ 2 x,γ ≡ ∂ 2 /∂γ∂x. We note that the space of angles θ ijk is assumed here to only consist of those angles that were not eliminated from the corresponding rows of Q.
Combining now Eqs. (11), (14) and (16), we arrive at a simple expression for the shear modulus: A similar expression was first put forward in [23] for the case of random networks of relaxed Hookean springs.

B. Scaling argument for hyperstatic moduli
We follow a similar line of argumentation as presented in [23] for elastic moduli of random networks of relaxed Hookean springs. From general considerations (see e.g. [24] for a detailed discussion) it can be shown that for z > z c where |φ are the zero-modes of the operator QQ T , namely they satisfy These objects are akin to the so-called states of self stress studied intensively in the context of the jamming transition [24][25][26][27] and the physics of topological metamaterials [28,29]. Simple counting arguments [23,25,27] suggest that in a system of size N with coordination z there are N (z − z c ) orthonormal modes |φ . Assuming these modes are extended and random objects one expects ∂ γ θ|φ ∼ O(1), and therefore we predict in the hyperstatic regime, and in the limits z → z + c and µ → 0. This prediction is also in agreement with Effective Medium calculations, see e.g. [30,31]; it suggests that the scaling exponent f = 1, c.f. Fig. 3.
In Fig. 4b we show our measurements of the shear modulus for k θ = 1 and k r = 0. We do not find perfect agreement with the theoretical prediction. However, this disagreement with the scaling argument and mean field theory seems not to be unique to our system in which bending interactions dominate; the shear modulus of randomly diluted hyperstatic spring networks in 2D shows a similar super linear scaling with δz too, as demonstrated in previous works [32,33]. We speculate [34] that dimensionality may be playing a role. Uncovering the origin of these observed disagreements is left for future work.
As for the bulk modulus, we note that isotropic expansions leave the angles θ ijk unchanged. As a consequence, we expect the bulk modulus to trivially scale with k r throughout the entire coordination range, and to trivially grow with increasing coordination. In Fig. 5 we show that this is indeed the case.

VI. THE HYPOSTATIC REGIME z < zc
To study the elasticity of our model in the hypostatic regime, we fix k r and send k θ → ∞, resulting in µ → 0. In this route of taking the µ → 0 limit the bond-bending interactions can be treated as geometrical constraints, as mentioned in Sect. IV.
In what follows we assume z < z c , we fix k r and send k θ → ∞, resulting in µ → 0. We aim at incorporating the geometric constraints of fixed angles while deriving microscopic expressions for elastic moduli; a similar derivation was recently presented in [20] in the context of the nonlinear mechanics of biopolymeric fibrous materials. In what follows we will consider the angles θ ijk as being entirely fixed; the potential energy of the material then reduces to i.e. we only consider the central-force part of the energy. We aim at taking the derivative of the energy with respect to the imposed deformation, under two sets of constraints: (i) the system must remain in mechanical equilibrium along the deformation, and (ii) the geometric constraints of fixed angles must be satisfied by the total displacements of the network's nodes. In order to satisfy condition (i), we introduce Lagrange multipliers τ ijk that correspond to clamping torques associated with each fixed angle θ ijk ; in terms of these torques and the potential energy given by Eq. (21), the net force experienced by the i th node must vanish, namely Condition (i) is satisfied by demanding that not only is the net force zero, but it also does not change under imposed deformations, namely where the nonaffine displacements y are additional displacements that the nodes must perform on top of the imposed deformation, in order to satisfy the mechanical equilibrium constraints, and leave the angles invariant under the external deformation. The latter requirement can be expressed as a constraint equation for each angle θ ijk that reads Setting k r = 1, and defining the linear operator [3] we show in Appendix B that Eq. (23) and (24) can be incorporated into a single relation as where we have used the notations ∂ γ r ij ≡ ∂r ij /∂γ and ∂ γ θ ijk ≡ ∂θ ijk /∂γ. Eq. (26) forms a closed linear system that can be inverted in favor of the nonaffine displacement field y and the torque variationsτ .

A. Elastic moduli in the hypostatic regime
We are now in the position to derive expressions for elastic moduli, in the limit µ → 0 obtained by taking k θ → ∞ and k r = 1; the second derivatives of the energy with respect to deformation reads Assuming again an unstressed material, namely that all torques τ ijk = 0 vanish then following Eq. (22) Notice further that for the simple form of the potential Eq. (21), and by setting k r = 1, one finds and Combining Eqs. (28), (29) and (30) with Eq. (27), we arrive at an expression for elastic moduli in the hypostatic regime z < z c in the limit µ → 0: B. Scaling arguments for hypostatic moduli Examining Eq. (31) it is clear that if the characteristic scale of nonaffine displacements y ≡ y|y /N diverges as z → z c , one expects to observe scaling laws of elastic moduli with respect to z−z c . To understand the behavior of the nonaffine velocities, we rearrange Eq. (26) in favor of y as We now perform a mean-field approximation, and consider a potential energy that consists of Hookean springs of unit stiffness that connect each node to its absolute initial position x (0) , namely In this case the dynamical matrix ∂ 2 U mf /∂x∂x reduces to the identity tensor I (instead of S T S, as given by Eq. (29)), and the nonaffine displacements assume the simple form Recall that when z = z c a single solution |φ to the equation Q T |φ = 0 appears. We therefore expect the operator QQ T to posses lower and lower frequency modes as z approaches z c , that, in turn, give rise to diverging nonaffine displacements. This is akin to the behavior observed for the operator SS T in floppy networks [5]. We next note that the first term in the right hand side (RHS) of Eq. (34) is regular, and in Appendix C we argue that the first term in the brackets on the RHS of Eq. (34) (the term involving Q) can be neglected compared to the second term as z → z c . Therefore, as z → z c , the characteristic scale of nonaffine displacements follows Assuming that the operator Q is random, it has been shown in [35] that the spectrum of concatenations of the form QQ T should depend on the dimensions of the operator Q; in particular, one expects the density of states (i.e. the distribution of the square root of the eigenvalues) of QQ T to follow the Marchenko-Pastur distribution [35], which in the small frequency and z → z c limits takes the form with ω ∼ z c − z. Assuming next that the eigenmodes of QQ T are random, extended objects (similarly to the arguments made before Eq. (20)), we estimate Finally, following Eq. (31) we expect that to leading order in y, G ∼ y 2 . We therefore conclude this Section with the prediction in the hypostatic regime, in the limits z → z − c and µ → 0. In Fig. 4a we plot G vs. δz for our hypostatic systems, and find excellent agreement with the theoretical prediction Eq. (38). We note that similar results were shown for an elastic system subjected to radial constraints in [20]. Eq. (38) implies the scaling relation f − φ = −1; using the mean-field exponent f = 1 (see Sect. V B), one expects φ = 2. However, the best collapse in Fig. 3 is achieved using f = 1.25 and φ = 2.25, consistent with our predicted scaling relation, and with the direct measurement in the hyperstatic regime shown in Fig. 4.   [24,36].

VII. DIVERGING LENGTHSCALE IN THE HYPERSTATIC REGIME
In Sections V and VI we rationalize the scaling of the shear modulus in the limit µ → 0 both in the hyperstatic and hypostatic regimes, respectively. Building on the discussions held in both of these Sections, we now make a scaling argument that predicts a diverging lengthscale in the hyperstatic regime, that is expected to follow The argument is made as follows; consider the vibrational spectra of our bending dominated disordered networks with µ = 0 and z larger than but close to the critical coordination z c . As shown in Sect. V A, for µ = 0 the dynamical matrix assumes the form We note that, apart from possible zero modes, the spectra of the concatenations QQ T and of Q T Q are identical [35]. We therefore expect the occurrence of a plateau of disordered vibrational modes above the frequency scale ω ∼ z−z c [35], as discussed in Sect. VI B.
On the other hand, our system's potential energy is invariant to global translations, and therefore Goldstone modes are expected to be present at small frequencies. In particular, the frequencyω min of the longest-wavelength Goldstone mode depends on the system size L and the shear modulus G asω Consider now the response to a localized force perturbation in a random bending-dominated network of coordination z and linear size L; ifω min ω , we expect the far field of the response to the local perturbation to exhibit a continuum-elastic-like structure, whilst ifω min ω no continuum-elastic-like response is expected since low-frequency disordered modes will overwhelm the response. Since the two frequency scalesω min and ω become comparable when L is of the order of 1/ √ z−z c , we expect to see a signature of a diverging length c ∼ 1/ √ z−z c in the spatial structure of the response to point perturbations.
We stress that both in floppy (hypostatic) [5] and in hyperstatic random networks [36] the length c was observed, as well in several other previous works. Ref. [24] provides a comprehensive summary of additional observations of the length c in the existing literature.
To test our argument, we select randomly an angle θ ijk impose a localized force perturbation of the form An example of this force can be seen in Fig. 1b. The linear displacement response to the force f m reads We denote any two angles in the system θ 1 and θ 2 , and by r 12 the distance between the nodes associated with these angles. In Fig. 6 an example of the response field is visualized for two different coordinations. The growing length scale is clearly showing. We next define where r12 denotes the average over all pairs of angles θ 1 , θ 2 separated by a distance r 12 . Fig. 7 shows the results of our numerical calculations of C(r): in panel (a) we show the mean of 50 different force perturbations. In panel (b) we plot the same data shown in panel (a), this time scaled by r 4 δz 2.5 , and plotted against r √ δz ∼ r/ c , clearly revealing that the lengthscale governing the transition to a continuum-like response is c ∼ 1/ √ z−z c , consistent with our scaling argument.

VIII. AUXETICITY
In this final Section, before summarizing our work, we show the auxetic behavior of our model material. Auxeticity is quantified via the Poisson's ratio ν, defined in 2D as Materials possessing ν = 1/2 are known as incompressible, whereas materials with ν < 0 are termed auxetic; the latter are nongeneric, and as such draw attention in the field of architected metamaterials [37][38][39].
In Fig. 8 we plot the Poisson's ratio ν averaged over 20 realizations as a function of the coordination, for various values of the stiffness ratio µ. We find that our system becomes auxetic in the entire range of µ explored (µ ≤ 10 −2 ), for z roughly larger than 3.5. Interestingly, in the limit µ → 0 we find a transition to perfect auxeticity µ = −1 at the isostatic point z c = 4. From the scaling behavior of G discussed in Sect. VI we expect ν+1 ∼ z c −z as z → z c .

IX. SUMMARY AND DISCUSSION
In this work we explored the elastic properties of a model system that represents a material whose mechanics is dominated by bond-bending interactions, rather than the commonly-studied steric or radial interactions. Our study was motivated by recent intriguing experimental work by Schall and coworkers [14], who fabricated colloidal superstructures using critical Casimir forces [40], in which bond-bending interactions between the constituent patchy colloidal particles were shown to be much stiffer than radial ones. In our model system we observe numerically and rationalize theoretically the existence of an angle-preserving isostatic point z c = 4 that governs an underlying jamming transition observed when the ratio of stretching to bending stiffnesses of the interactions vanishes. The jamming behavior we observe as the coordination is made to approach the critical coordination from below is reminiscent of the strain stiffening transition observed in fibrous biological materials [18][19][20][21].
Our theoretical arguments for the scaling of elastic moduli in the limits µ → 0, z → z + c (the hyperstatic regime) and z → z − c (the hypostatic regime) follow closely two previous theoretical approaches to the jamming [23] and the strain-stiffening [20] transitions in random networks, respectively. In the hyperstatic regime, the arguments put forward by Wyart [23] using the operator S (see Eq. (25)) -that represents the constraints assocaited with radial interactions -were applied here using the operator Q that represents constraints associated with conserving the angles of our networks. In the hypostatic regime, our argumentation echos the framework and reasoning presented in [20], where again the role of the operator S in [20] is played by the operator Q in the present work. The applicability of these approaches establishes their generality, and highlights the key physical ingredient governing the mechanics near jamming and stiffening transitions -the interplay between interactions-induced constraints and configurational degrees of freedom.
Here we argue that the vector |b = QS T |∂ γ r − |∂ γ θ can be approximated as |b ∼ −|∂ γ θ near the critical point. The reason is that the first of these two terms features an additional factor Q; since this term is expected to be regular when contracted with Q T (QQ T ) −1 (and see similar discussion in [20]), it can be neglected in our scaling analysis.