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

TL;DR

This paper investigates the elasticity and phase transitions of 2D networks with strong bond-bending interactions, revealing a stiffening and auxeticity transition near the isostatic point, with implications for designing tunable materials.

Contribution

It introduces a theoretical framework for analyzing networks with dominant bond-bending forces and identifies a new auxeticity transition in such systems.

Findings

Existence of a generic isostatic point at z_c=4 with bond-bending interactions

Emergence of a strain-induced stiffening transition near the isostatic point

Identification of an auxeticity transition with Poisson's ratio approaching -1

Abstract

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.

Figures

Equations

$$U=\frac{k_{\theta }{\bar{\ell }}^2}{2}\sum _{⟨i,j,k⟩}\Delta {{\theta }_{ijk}}^2+\frac{k_r}{2}\sum _{⟨i,j⟩}\Delta {r_{ij}}^2,$$

$$H=\left(\begin{array}[]{cc}1+\eta&\gamma\\ 0&1+\eta\end{array}\right)\,.$$

$$\epsilon(\gamma,\eta)=\frac{1}{2}\!\left(H^{T}\cdot H-{\cal I}\right)=\frac{1}{2}\!\left(\begin{array}[]{cc}2\eta\!+\!\eta^{2}&\gamma\!+\!\gamma\eta\\ \gamma\!+\!\gamma\eta&2\eta\!+\!\eta^{2}\!+\!\gamma^{2}\end{array}\right)\,,$$

$$\delta r_{ij}\simeq \frac{\text{mathitbf}{x}_{ij}\cdot ϵ\cdot \text{mathitbf}{x}_{ij}}{r_{ij}}-\frac{1}{2}\frac{(\text{mathitbf}{x}_{ij}\cdot ϵ\cdot \text{mathitbf}{x}_{ij}{)}^2}{r_{ij}^3},$$

$$G\equiv \frac{1}{V}\frac{d^{2}U}{d{\gamma }^2},$$

$$K\equiv \frac{1}{V}\left(\frac{1}{4}\frac{d^{2}U}{d{\eta }^2}-\frac{1}{2}\frac{dU}{d\eta }\right).$$

$$\frac{\partial {\theta }_{ijk}}{\partial \text{mathitbf}{x}_{\ell }}\cdot \text{mathitbf}{u}_{\ell }=0,$$

$${\mathcal{Q}}_{ijk,\ell }\equiv \frac{\partial {\theta }_{ijk}}{\partial \text{mathitbf}{x}_{\ell }},$$

$$\mathcal{F}=\frac{1}{2}Nz-N.$$

$$\mbox{rank}(Q)=Nz-N-\big{(}{\textstyle\frac{1}{2}}Nz-N\big{)}={\textstyle\frac{1}{2}}Nz\,,$$

$$U_{k_r=0}=\frac{k_{\theta }{\bar{\ell }}^2}{2}\sum _{⟨i,j,k⟩}\Delta {{\theta }_{ijk}}^2.$$

$$G=\frac{1}{V}\left(\frac{{\partial }^{2}U}{\partial {\gamma }^2}-\frac{{\partial }^{2}U}{\partial \gamma \partial \text{mathitbf}{x}_{\ell }}\cdot {\mathcal{M}}_{\ell m}^{-1}\cdot \frac{{\partial }^{2}U}{\partial \text{mathitbf}{x}_m\partial \gamma }\right),$$

$$\mathcal{M}\equiv \frac{{\partial }^{2}U}{\partial \text{mathitbf}x\partial \text{mathitbf}x}$$

$${\mathcal{M}}_{\ell m}=k_{\theta }{\bar{\ell }}^2\sum _{⟨i,j,k⟩}\frac{\partial {\theta }_{ijk}}{\partial \text{mathitbf}{x}_{\ell }}\frac{\partial {\theta }_{ijk}}{\partial \text{mathitbf}{x}_m}.$$

$$\mathcal{M}={\mathcal{Q}}^T\mathcal{Q},$$

$$\frac{{\partial }^{2}U}{\partial \text{mathitbf}{x}_{\ell }\partial \gamma }=\sum _{⟨i,j,k⟩}\frac{\partial {\theta }_{ijk}}{\partial \text{mathitbf}{x}_{\ell }}\frac{\partial {\theta }_{ijk}}{\partial \gamma },$$

$$|{\partial }_{\text{mathitbf}x,\gamma }^{2}U⟩={\mathcal{Q}}^T|{\partial }_{\gamma }\theta ⟩,$$

$$G=\frac{⟨{\partial }_{\gamma }\theta |\mathcal{I}-\mathcal{Q}({\mathcal{Q}}^T\mathcal{Q}{)}^{-1}{\mathcal{Q}}^T|{\partial }_{\gamma }\theta ⟩}{V}.$$

$$\mathcal{I}-\mathcal{Q}({\mathcal{Q}}^T\mathcal{Q}{)}^{-1}{\mathcal{Q}}^T=\sum _{\ell }|{\varphi }_{\ell }⟩⟨{\varphi }_{\ell }|,$$

$${\cal Q}^{T}|\phi_{\ell}\rangle=\mbox{\boldmath$0$}\,.$$

$$G=\frac{1}{V}\sum _{\ell }⟨{\partial }_{\gamma }\theta |{\varphi }_{\ell }{⟩}^2\sim z-z_c,$$

$$U_{k_{\theta }=0}=\frac{k_r}{2}\sum _{⟨i,j⟩}\Delta {r_{ij}}^2,$$

$$\text{mathitbf}{f}_m=\sum _{⟨i,j,k⟩}{\tau }_{ijk}\frac{\partial {\theta }_{ijk}}{\partial \text{mathitbf}{x}_m}-\frac{\partial U}{\partial \text{mathitbf}{x}_m}=0.$$

$$\frac{d\text{mathitbf}{f}_m}{d\gamma }=\frac{\partial \text{mathitbf}{f}_m}{\partial \gamma }+\frac{\partial \text{mathitbf}{f}_m}{\partial \text{mathitbf}{x}_n}\cdot \text{mathitbf}{y}_n+\sum _{⟨i,j,k⟩}\frac{\partial \text{mathitbf}{f}_m}{\partial {\tau }_{ijk}}\frac{d{\tau }_{ijk}}{d\gamma }=0,$$

$$\frac{d{\theta }_{ijk}}{d\gamma }=\frac{\partial {\theta }_{ijk}}{\partial \gamma }+\frac{\partial {\theta }_{ijk}}{\partial \text{mathitbf}{x}_m}\cdot \text{mathitbf}{y}_m=0.$$

$${\mathcal{S}}_{ij,k}\equiv \frac{\partial r_{ij}}{\partial \text{mathitbf}{x}_k},$$

$$\left(\begin{array}[]{cc}{\cal S}^{T}{\cal S}&-{\cal Q}^{T}\\ -{\cal Q}&0\end{array}\right)\left(\begin{array}[]{c}|\mathitbf y\rangle\\ |\dot{\tau}\rangle\end{array}\right)=\left(\begin{array}[]{c}-{\cal S}^{T}|\partial_{\gamma}r\rangle\\ |\partial_{\gamma}\theta\rangle\end{array}\right)\,,$$

$$\frac{d^{2}U}{d{\gamma }^2}$$

$$\frac{\partial U}{\partial\mathitbf x_{k}}=\mbox{\boldmath$0$}\,.$$

$$\frac{{\partial }^{2}U}{\partial \text{mathitbf}{x}_k\partial \text{mathitbf}{x}_{\ell }}=\sum _{⟨i,j⟩}\frac{\partial r_{ij}}{\partial \text{mathitbf}{x}_k}\frac{\partial r_{ij}}{\partial \text{mathitbf}{x}_{\ell }}\leftrightarrow {\mathcal{S}}^T\mathcal{S},$$

Full text

I Introduction

In 1864 Maxwell spelled out a criterion that frames of freely-hinged struts need to satisfy in order to be mechanically stable Maxwell (1864): if the average connectivity $z$ is higher than a threshold value $z_{c}\!\equiv 2{\,\mathchar 22\relax\mkern-12.0mud}$ in ${\,\mathchar 22\relax\mkern-12.0mud}$ 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, and in the absence of over-constrained clusters Ellenbroek et al. (2015)). In frames with $z\!<\!z_{c}$ collective modes exist that are floppy Calladine (1978); Heussinger and Frey (2006); Düring et al. (2013), 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\!\to\!z_{c}$ is known as the jamming transition O’Hern et al. (2003); van Hecke (2010); Liu and Nagel (2010), and has been related to various mechanical and dynamical phenomena such as the divergence of viscosity in non-Brownian suspensions DeGiuli et al. (2015) and the fragility of chalcogenide glass formers Yan and Wyart (2014).

Most studies of jamming phenomena focus on the role of steric interactions, often modeled by some form of central forces e.g. Hookean springs Ellenbroek et al. (2009) or hard-sphere repulsions Charbonneau et al. (2014); DeGiuli et al. (2014a). 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 gels with tunable interactions; in particular, we were inspired by the work of Schall et al. Nguyen et al. (2017), 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 Nguyen et al. (2017) 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 pri (a).

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_{\theta}$, respectively. We consider the ratio $\mu\!\equiv\!k_{r}/k_{\theta}$ 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 $\mu$. We show that as the limit $\mu\!\to\!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 Maxwell threshold Maxwell (1864) $z_{c}\!=\!4$ at which the generic, ‘central-force’ isostatic point occurs.

We take two complementary routes in order to study theoretically the limit $\mu\!\to\!0$ at which the scaling behavior of elastic moduli emerges. First, we fix $k_{\theta}$ 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 $\ell_{c}\!\sim\!(z-z_{c})^{-1/2}$.

Even more intruiging is the limit $\mu\!\to\!0$ obtained by fixing $k_{r}$ and sending $k_{\theta}\!\to\!\infty$; 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 $\sim 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 Wang (2011). However, $\nu$ rapidly decreases as $z$ is increased. Interestingly, in the limit $\mu\!\to\!0$ the material approaches perfect auxeticity $\nu\!\to\!-1$ as $z\!\to\!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 $\mu\!\equiv\!k_{r}/k_{\theta}$ 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 $\mu\!\to\!0$, while Sect. VI presents a theoretical framework that allows us to study the hypostatic regime in the limit $\mu\!\to\!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

[TABLE]

where $k_{\theta}$ and $k_{r}$ denote the bond-bending and central-force stiffnesses, respectively, and $\bar{\ell}$ denotes the microscopic units of length. The first term on the RHS of Eq. (1) represents a sum over all angles $\theta_{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 $\Delta\theta_{ijk}\!\equiv\!\theta_{ijk}\!-\!\theta_{ijk}^{(0)}$ with $\theta_{ijk}^{(0)}$ 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 $\Delta r_{ij}\!\equiv\!r_{ij}\!-\!\ell_{ij}$ where $\ell_{ij}$ denotes the rest-length of the edge connecting the $i,j$ pair of nodes. Lengths are expressed in terms of $\bar{\ell}$ 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 that result from the application of the affine transformation $H(\gamma,\eta)$ on coordinates $\mathitbf x$ as $\mathitbf x\!\to\!H\!\cdot\!\mathitbf x$. The transformation $H(\gamma,\eta)$ is most conveniently parameterized by simple shear and expansive strain parameters $\gamma$ and $\eta$, respectively, and has the following form in 2D

[TABLE]

Using this transformation, the strain tensor $\epsilon$ assumes the form

[TABLE]

where ${\cal I}$ represents the identity tensor. We note that quadratic terms in the strain parameters $\eta$ and $\gamma$ are kept in the entries of $\epsilon$ since we will be interested in second-order derivatives with respect to those parameters, as shown below. Given a strain tensor that describes an imposed deformation mode, distances $r_{ij}\!=\!\sqrt{\mathitbf x_{ij}\cdot\mathitbf x_{ij}}$ between the coordinates of any two nodes $\mathitbf x_{i},\mathitbf x_{j}$ vary under such imposed deformations according to

[TABLE]

where $\mathitbf x_{ij}\!\equiv\!\mathitbf x_{j}\!-\!\mathitbf x_{i}$. We note importantly that in what follows we only consider either simple shear or expansive strains, and not combinations of these.

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

[TABLE]

whereas the athermal bulk modulus is given by

[TABLE]

Total dervatives e.g. $d/d\gamma$ are understood as taken in the athermal limit, i.e. under the constraints dictated by mechanical equilibrium Lutsko (1989).

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 $\mu$ 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_{\theta}$ vs. the coordination $z$; noticeably, $z\!=\!4$ marks the onset of an underlying jamming transition, that becomes more pronounced as $\mu\!\to\!0$. For small $\mu$ the ratio $G/k_{\theta}$ grows by several orders of magnitude as $z$ approaches the critical coordination $z_{c}\!=\!4$, in a fashion reminiscent of the strain-stiffening transition observed upon deformation of biopolymeric fibrous materials Broedersz et al. (2011); Rens et al. (2016); Jansen et al. (2018); Rens et al. (2018); Merkel et al. (2019). Far above $z_{c}$ the ratio $G/k_{\theta}$ becomes roughly independent of $\mu$, indicating that in this regime $G\!\sim\!k_{\theta}$.

In the right panel of Fig. 2 we show the same sample-to-sample means of the shear modulus, this time rescaled by $k_{r}$, to find that deep in the hypostatic regime $z\!<\!4$, $G\!\sim k_{r}$. The data pertaining to $\mu\!=\!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_{\theta}$ is rescaled by $|\delta z|^{f}$ and plotted against the ratio $k_{r}/|\delta z|^{\phi}$ where $\delta z\!\equiv\!z\!-\!z_{c}$ denotes the distance to the critical coordination $z_{c}\!\equiv\!4$. The best collapse is found using the exponents $f\!=\!1.25$ and $\phi\!=\!2.25$. In what follows we will argue that the mean field exponents are $f\!=\!1$ and $\phi\!=\!2$, and discuss the discrepancy we find with the mean field predictions.

To understand the scaling behavior of the shear modulus 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 $\mu\!\to\!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

In this Section we extend the Maxwell-Calladine linear-algebraic constructions for the rigidity of frames of struts Maxwell (1864); Calladine (1978) to 3-body (bond-bending) geometries, in order to establish that the angle-preserving isostatic point of 2D planar networks is $z_{c}=4$.

The emergence of an isostatic point in 2D planar networks becomes apparent when the limit $\mu\!\to\!0$ is considered, as seen in Fig. 2. In the hypostatic regime $z\!\leq\!z_{c}$ the limit $\mu\!\to\!0$ can be obtained by fixing $k_{r}$ and sending $k_{\theta}\!\to\!\infty$. Under these circumstances the bond-bending interactions can be treated as geometric constraints; a displacement field $\mathitbf u_{\ell}$ on the network’s nodes will leave the angles $\theta_{ijk}$ invariant (to leading order in $\mathitbf u_{\ell}$) if it satisfies

[TABLE]

for every angle $\theta_{ijk}$, where here and in what follows repeated node indices are understood to be summed over. An example of the field $\frac{\partial\theta_{ijk}}{\partial\mathitbf x_{\ell}}$ is shown in Fig. 1b.

It is useful to define the linear operator

[TABLE]

which takes vectors from the space of the nodes’ coordinates to the space defined by the entire set of angles. The number of rows ${\cal Q}$ features is equal to the total number of angles $Nz$, whereas the number of columns is $2N$, each represented a spatial coordinate of a node. Nontrivial solutions to Eq. 7 are expected to exist if the rank of the operator ${\cal Q}$ is smaller than the number of degrees of freedom $2N$ available to the displacement field $\mathitbf u$.

We next argue that the rank of ${\cal 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^{\mbox{\tiny 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_{i}^{\mbox{\tiny th}}$ angle as well, meaning that a single angle for each node can be selected, and the row of ${\cal 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 $Nz$ rows of ${\cal 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^{\mbox{\tiny th}}$ angle, further reducing the number of independent rows of ${\cal Q}$ by the number ${\cal F}$ of faces of the network. The latter is related to the number of nodes $N$ and the number of edges $\frac{1}{2}Nz$ via Euler’s formula for a planar graph embedded on a torus

[TABLE]

The total number of independent rows left after the eliminations described above is thus

[TABLE]

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 Goodrich et al. (2014).

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}\!=\!2{\,\mathchar 22\relax\mkern-12.0mud}$ of networks of rigid struts Maxwell (1864); Calladine (1978) discussed intensively in the jamming literature, see e.g. van Hecke (2010); Liu and Nagel (2010).

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\!>\!z_{c}$

It is illuminating to study the elasticity of our model in the hyperstatic regime by fixing $k_{\theta}\!=\!1$ and sending $k_{r}\!\to\!0$, resulting in $\mu\!\to\!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

[TABLE]

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

V.1 Elastic moduli in the hyperstatic regime

The microscopic expression for the athermal shear modulus $G$ for systems governed by a potential energy $U$ reads Lutsko (1989)

[TABLE]

where

[TABLE]

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

[TABLE]

Working in units such that $\bar{\ell}\!=\!1$ and $k_{\theta}\!=\!1$, employing Eq. 8 and adopting bra-ket and matrix notations, the dynamical matrix takes the form

[TABLE]

where here and in what follows we assume that the redundent rows of ${\cal Q}$ have been eliminated as described in Sect. IV.

We next see that under the circumstances of fixing $k_{\theta}\!=\!1$ and sending $k_{r}$ to zero,

[TABLE]

which can be expressed using our bra-ket notation as

[TABLE]

where we denoted $\partial_{\gamma}\!\equiv\!\partial/\partial\gamma$ and $\partial^{2}_{\mathitbf x,\gamma}\!\equiv\!\partial^{2}/\partial\gamma\partial\mathitbf x$. We note that the space of angles $\theta_{ijk}$ is assumed here to only consist of those angles that were not eliminated from the corresponding rows of ${\cal Q}$.

Combining now Eqs. (12),(15) and (17), we arrive at a simple expression for the shear modulus:

[TABLE]

A similar expression was first put forward in Wyart, M. (2005) for the case of random networks of relaxed Hookean springs.

V.2 Scaling argument for hyperstatic moduli

We follow a similar line of argumentation as presented in Wyart, M. (2005) for elastic moduli of random networks of relaxed Hookean springs. From general considerations (see e.g. Lerner (2018) for a detailed discussion) it can be shown that for $z\!>\!z_{c}$

[TABLE]

where $|\phi_{\ell}\rangle$ are the zero-modes of the operator ${\cal Q}{\cal Q}^{T}$, namely they satisfy

[TABLE]

These objects are akin to the so-called states of self stress studied intensively in the context of the jamming transition Ellenbroek et al. (2009, 2015); Sussman et al. (2016); Lerner (2018) and the physics of topological metamaterials Kane and Lubensky (2014); Paulose et al. (2015). Simple counting arguments Wyart, M. (2005); Ellenbroek et al. (2009); Sussman et al. (2016) suggest that in a system of size $N$ with coordination $z$ there are $N(z\!-\!z_{c})$ orthonormal modes $|\phi_{\ell}\rangle$. Assuming these modes are extended and random objects one expects $\langle\partial_{\gamma}\theta|\phi_{\ell}\rangle\!\sim\!{\cal O}(1)$, and therefore we predict

[TABLE]

in the hyperstatic regime, and in the limits $z\!\to\!z_{c}^{+}$ and $\mu\!\to\!0$. This prediction is also in agreement with Effective Medium calculations, see e.g. Wyart (2010); DeGiuli et al. (2014b); 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_{\theta}\!=\!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 $\delta z$ too, as demonstrated in previous works Tighe (2011); Baumgarten and Tighe (2018). We speculate pri (b) 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 $\theta_{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\!<\!z_{c}$

To study the elasticity of our model in the hypostatic regime, we fix $k_{r}$ and send $k_{\theta}\!\to\!\infty$, resulting in $\mu\!\to\!0$. In this route of taking the $\mu\!\to\!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_{\theta}\!\to\!\infty$, resulting in $\mu\!\to\!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 Rens et al. (2018) in the context of the nonlinear mechanics of biopolymeric fibrous materials. In what follows we will consider the angles $\theta_{ijk}$ as being entirely fixed; the potential energy of the material then reduces to

[TABLE]

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 $\tau_{ijk}$ that correspond to clamping torques associated with each fixed angle $\theta_{ijk}$; in terms of these torques and the potential energy given by Eq. (22), the net force experienced by the $i^{\mbox{\tiny th}}$ node must vanish, namely

[TABLE]

Condition $(i)$ is satisfied by demanding that not only is the net force zero, but it also does not change under imposed deformations, namely

[TABLE]

where the nonaffine displacements $\mathitbf 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 $\theta_{ijk}$ that reads

[TABLE]

Setting $k_{r}\!=\!1$, and defining the linear operator Calladine (1978)

[TABLE]

we show in Appendix B that Eq. (24) and (25) can be incorporated into a single relation as

[TABLE]

where we have used the notations $\partial_{\gamma}r_{ij}\!\equiv\!\partial r_{ij}/\partial\gamma$ and $\partial_{\gamma}\theta_{ijk}\!\equiv\!\partial\theta_{ijk}/\partial\gamma$. Eq. (27) forms a closed linear system that can be inverted in favor of the nonaffine displacement field $\mathitbf y$ and the torque variations $\dot{\tau}$.

VI.1 Elastic moduli in the hypostatic regime

We are now in the position to derive expressions for elastic moduli, in the limit $\mu\!\to\!0$ obtained by taking $k_{\theta}\!\to\!\infty$ and $k_{r}\!=\!1$; the second derivatives of the energy with respect to deformation reads

[TABLE]

Assuming again an unstressed material, namely that all torques $\tau_{ijk}\!=\!0$ vanish then following Eq. (23)

[TABLE]

Notice further that for the simple form of the potential Eq. (22), and by setting $k_{r}\!=\!1$, one finds

[TABLE]

and

[TABLE]

Combining Eqs. (29), (30) and (31) with Eq. (28), we arrive at an expression for elastic moduli in the hypostatic regime $z\!<\!z_{c}$ in the limit $\mu\!\to\!0$:

[TABLE]

VI.2 Scaling arguments for hypostatic moduli

Examining Eq. (32) it is clear that if the characteristic scale of nonaffine displacements $y\!\equiv\!\sqrt{\langle\mathitbf y|\mathitbf y\rangle/N}$ diverges as $z\!\to\!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. (27) in favor of $\mathitbf y$ as

[TABLE]

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 $\mathitbf x^{(0)}$, namely

[TABLE]

In this case the dynamical matrix $\partial^{2}U_{\mbox{\tiny mf}}/\partial\mathitbf x\partial\mathitbf x$ reduces to the identity tensor ${\cal I}$ (instead of ${\cal S}^{T}{\cal S}$, as given by Eq. (30)), and the nonaffine displacements assume the simple form

[TABLE]

Recall that when $z\!=\!z_{c}$ a single solution $|\phi\rangle$ to the equation {\cal Q}^{T}|\phi\rangle\!=\!\mbox{\boldmath0} appears. We therefore expect the operator ${\cal Q}{\cal Q}^{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 ${\cal S}{\cal S}^{T}$ in floppy networks Düring et al. (2013).

We next note that the first term in the right hand side (RHS) of Eq. (35) is regular, and in Appendix C we argue that the first term in the brackets on the RHS of Eq. (35) (the term involving ${\cal Q}$) can be neglected compared to the second term as $z\!\to\!z_{c}$. Therefore, as $z\!\to\!z_{c}$, the characteristic scale of nonaffine displacements follows

[TABLE]

Assuming that the operator ${\cal Q}$ is random, it has been shown in Beltukov (2015) that the spectrum of concatenations of the form ${\cal Q}{\cal Q}^{T}$ should depend on the dimensions of the operator ${\cal Q}$; in particular, one expects the density of states (i.e. the distribution of the square root of the eigenvalues) of ${\cal Q}{\cal Q}^{T}$ to follow the Marchenko-Pastur distribution Beltukov (2015), which in the small frequency and $z\!\to\!z_{c}$ limits takes the form

[TABLE]

with $\omega_{\star}\!\sim\!z_{c}\!-\!z$. Assuming next that the eigenmodes of ${\cal Q}{\cal Q}^{T}$ are random, extended objects (similarly to the arguments made before Eq. (21)), we estimate

[TABLE]

Finally, following Eq. (32) we expect that to leading order in $y$, $G\!\sim\!y^{2}$. We therefore conclude this Section with the prediction

[TABLE]

in the hypostatic regime, in the limits $z\!\to\!z_{c}^{-}$ and $\mu\!\to\!0$.

In Fig. 4a we plot $G$ vs. $\delta z$ for our hypostatic systems, and find excellent agreement with the theoretical prediction Eq. (39). We note that similar results were shown for an elastic system subjected to radial constraints in Rens et al. (2018). Eq. (39) implies the scaling relation $f\!-\!\phi=-1$; using the mean-field exponent $f\!=\!1$ (see Sect. V.2), one expects $\phi\!=\!2$. However, the best collapse in Fig. 3 is achieved using $f\!=\!1.25$ and $\phi\!=\!2.25$, consistent with our predicted scaling relation, and with the direct measurement in the hyperstatic regime shown in Fig. 4.

VII Diverging lengthscale in the hyperstatic regime

In Sections V and VI we rationalize the scaling of the shear modulus in the limit $\mu\!\to\!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

[TABLE]

The argument is made as follows; consider the vibrational spectra of our bending dominated disordered networks with $\mu\!=\!0$ and $z$ larger than but close to the critical coordination $z_{c}$. As shown in Sect. V.1, for $\mu\!=\!0$ the dynamical matrix assumes the form

[TABLE]

We note that, apart from possible zero modes, the spectra of the concatenations ${\cal Q}{\cal Q}^{T}$ and of ${\cal Q}^{T}{\cal Q}$ are identical Beltukov (2015). We therefore expect the occurrence of a plateau of disordered vibrational modes above the frequency scale $\omega_{\star}\!\sim\!z\!-\!z_{c}$ Beltukov (2015), as discussed in Sect. VI.2.

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 $\tilde{\omega}_{\mbox{\tiny min}}$ of the longest-wavelength Goldstone mode depends on the system size $L$ and the shear modulus $G$ as

[TABLE]

Consider now the response to a localized force perturbation in a random bending-dominated network of coordination $z$ and linear size $L$; if $\tilde{\omega}_{\mbox{\tiny min}}\!\ll\!\omega_{\star}$, we expect the far field of the response to the local perturbation to exhibit a continuum-elastic-like structure, whilst if $\tilde{\omega}_{\mbox{\tiny min}}\!\gg\!\omega_{\star}$ no continuum-elastic-like response is expected since low-frequency disordered modes will overwhelm the response. Since the two frequency scales $\tilde{\omega}_{\mbox{\tiny min}}$ and $\omega_{\star}$ become comparable when $L$ is of the order of $1/\sqrt{z\!-\!z_{c}}$, we expect to see a signature of a diverging length $\ell_{c}\!\sim\!1/\sqrt{z\!-\!z_{c}}$ in the spatial structure of the response to point perturbations.

We stress that both in floppy (hypostatic) Düring et al. (2013) and in hyperstatic random networks Lerner et al. (2014) the length $\ell_{c}$ was observed, as well in several other previous works. Ref. Lerner (2018) provides a comprehensive summary of additional observations of the length $\ell_{c}$ in the existing literature.

To test our argument, we select randomly an angle $\theta_{ijk}$, and impose a localized force perturbation of the form

[TABLE]

An example of this force can be seen in Fig. 1b. The linear displacement response to the force $\mathitbf f_{m}$ reads

[TABLE]

We denote any two angles in the system $\theta_{1}$ and $\theta_{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

[TABLE]

where $\langle\rangle_{r_{12}}$ denotes the average over all pairs of angles $\theta_{1},\theta_{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}\delta z^{2.5}$, and plotted against $r\sqrt{\delta z}\!\sim\!r/\ell_{c}$, clearly revealing that the lengthscale governing the transition to a continuum-like response is $\ell_{c}\!\sim\!1/\sqrt{z\!-\!z_{c}}$, consistent with our scaling argument.

VIII Auxeticity

In this penultimate Section we discuss the auxetic behavior of our model material. Auxeticity is quantified via the Poisson’s ratio $\nu$, defined in 2D as

[TABLE]

Materials possessing $\nu\!=\!1/2$ are known as incompressible, whereas materials with $\nu\!<\!0$ are termed auxetic; the latter are nongeneric, and as such draw attention in the field of architected metamaterials Larsen et al. (1997); Yang et al. (2004); Reid et al. (2018).

In Fig. 8 we plot the Poisson’s ratio $\nu$ averaged over 20 realizations as a function of the coordination, for various values of the stiffness ratio $\mu$. We find that our system becomes auxetic in the entire range of $\mu$ explored ($\mu\!\leq\!10^{-2}$), for $z$ roughly larger than 3.5. As expected from the scaling of $G$ and $K$ discussed previously, in the limit $\mu\!\to\!0$ we find a transition to perfect auxeticity $\mu\!=\!-1$ at the isostatic point $z_{c}\!=\!4$. From the scaling behavior of $G$ discussed in Sect. VI we expect $\nu\!+\!1\sim z_{c}\!-\!z$ as $z\!\to\!z_{c}$, as supported by our data shown in Fig. 8.

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 Nguyen et al. (2017), who fabricated colloidal superstructures using critical Casimir forces Stuij et al. (2017), 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 Rens et al. (2016); Jansen et al. (2018); Rens et al. (2018); Merkel et al. (2019).

Our theoretical arguments for the scaling of elastic moduli in the limits $\mu\!\to\!0$, $z\!\to\!z_{c}^{+}$ (the hyperstatic regime) and $z\!\to\!z_{c}^{-}$ (the hypostatic regime) follow closely two previous theoretical approaches to the jamming Wyart, M. (2005) and the strain-stiffening Rens et al. (2018) transitions in random networks, respectively. In the hyperstatic regime, the arguments put forward by Wyart Wyart, M. (2005) using the operator ${\cal S}$ (see Eq. (26)) — that represents the constraints associated with radial interactions — were applied here using the operator ${\cal 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 Rens et al. (2018), where again the role of the operator ${\cal S}$ in Rens et al. (2018) is played by the operator ${\cal 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.

Acknowledgements.

We warmly thank G. Düring and C. Coulais for fruitful discussions. E. L. acknowledges support from the Netherlands Organisation for Scientific Research (NWO) (Vidi grant no. 680-47-554/3259). R.R. and E.L. acknowledge support from the Delta Institute for Theoretical Physics (D-ITP consortium), a program of NWO that is funded by the Dutch Ministry of Education, Culture and Science (OCW). Author contribution statement: both R.R. and E.L. designed and carried out the research, and wrote the article.

Appendix A Network generation protocol

The protocol for generating homogeneous networks consists of a two-step process. In the first step, a packing of soft disks is compressed up to a certain pressure. Using a bidisperse distribution of radii we prevent crystallization and ensures a disordered realization. A contact is assigned to each overlapping disc. We then distillate the contact network to obtain a highly coordinated network of edges and nodes, with an average of between 5 and 6 edges per node.

The second step is a biased bond dilution protocol that aims at maintaining homogeneity of the network, which is typically lost in a truely random dilution. The network with $N$ nodes and a set of edges $\cal E$. An edge $e_{ij}\!\in\!{\cal E}$ connects between a pair of neighboring nodes $i,j$, with ranks $z_{i}$ and $z_{j}$ respectively. For every edge $e_{ij}$ we define the sum of ranks $s_{ij}\!=\!z_{i}\!+\!z_{j}$ and the absolute difference $d_{ij}\!=\!\left|z_{i}\!-\!z_{j}\right|$. The edges with largest value for $s_{ij}$ and than smallest value of $d_{ij}$, will remove fluctuations in coordination of the network when removed. This min-max protocol typically produces a large set of edges that are eligible to be removed. An additional selection condition is checked to suppress the creation of sharp angles between edges. For each edge in the selected set the remaining angle after bond removal is calculated. The bond with the minimal remaining angle is now removed from network. This process is repeated until the desired connectivity is achieved.

As an example, in Fig. 1 (a), one finds the network diluted to a connectivity of $3.95$. To reiterate, this methods serves particularly well in generating a network with low connectivity fluctuations.

Appendix B Equations of motion for nonaffine displacements and clamping torques variations

In this Appendix we show how Eq. 27 is obtained from Eqs. 23, 24 and 25. The requirement of maintaining force balance under an imposed deformation is given by Eq. 24; it reads

[TABLE]

Where $\mathitbf f_{m}$ is the sum of all forces acting on node $m$, which is given by equation 25,

[TABLE]

Each term in Eq. 47 contains a partial derivative of the net force $\mathitbf f$, worked out below using the definition given by Eq. 48, namely

[TABLE]

where terms containing $\tau$ are dropped because in the undeformed state the torque forces are assumed to be identically zero.

To simply further consider the specified potential energy given by

[TABLE]

In undeformed state all pairwise distances in the network are at their rest length, so $\Delta r_{ij}\!=\!0$ for each $\left<i,j\right>$. Setting $k_{r}\!=\!1$, the second order derivative $\partial^{2}_{\mathitbf x,\mathitbf x}U$ reduces to

[TABLE]

and the mixed derivative $\partial^{2}_{\gamma,\mathitbf x}U$ to

[TABLE]

Finally we employ operator and braket notations, where

[TABLE]

Combining all of the above, Eq. 47 can now be compactly written as

[TABLE]

and similarly the invariance of the angles, given by Eq. 25, as

[TABLE]

with $\partial_{\gamma}\!\equiv\!\partial/\partial\gamma$. Combining equations 53 and 53 into one set of equations, we obtain

[TABLE]

as given by Eq. (27) in the main text.

Appendix C Leading order term in Eq. 35

In Eq. 35 the second term on the RHS contains the sum of two terms and reads ${\cal Q}^{T}({\cal Q}{\cal Q}^{T})^{-1}({\cal Q}{\cal S}^{T}|\partial_{\gamma}r\rangle-|\partial_{\gamma}\theta\rangle)$. Here we argue that the vector $|b\rangle=\left({\cal Q}{\cal S}^{T}|\partial_{\gamma}r\rangle-|\partial_{\gamma}\theta\rangle\right)$ can be approximated as $|b\rangle\sim-|\partial_{\gamma}\theta\rangle$ 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 ${\cal Q}^{T}({\cal Q}{\cal Q}^{T})^{-1}$ (and see similar discussion in Rens et al. (2018)), it can be neglected in our scaling analysis.