Force networks and jamming in shear deformed sphere packings

TL;DR

This paper investigates the conditions and structural features leading to shear jamming in sphere packings, demonstrating that isostatic contact networks and rigidity percolation are key factors in both frictionless and frictional systems.

Contribution

It provides a theoretical framework for understanding shear jamming, validating the isostatic condition at jamming, and linking rigidity percolation to the jamming transition in sphere packings.

Findings

Isostatic contact number Z = D + 1 holds at jamming

Shear jamming coincides with percolation of over-constrained regions

Rigidity percolation occurs before shear jamming in 2D systems

Abstract

The formation of self-organised structures that resist shear deformation have been discussed in the context of shear jamming and thickening[1-3], with frictional forces playing a key role. However, shear induces geometric features necessary for jamming even in frictionless packings[4]. We analyse conditions for jamming in such assemblies by solving force and torque balance conditions for their contact geometry. We demonstrate, and validate with frictional simulations, that the isostatic condition for mean contact number Z = D + 1 (for spatial dimension D = 2, 3) holds at jamming for both finite and infinite friction, above the random loose packing density. We show that the shear jamming threshold satisfies the marginal stability condition recently proposed for jamming in frictionless systems[5]. We perform rigidity percolation analysis[6,7] for D = 2 and find that rigidity percolation precedes shear jamming, which however coincides with the percolation of over-constrained regions, leading to the identification of an intermediate phase analogous to that observed in covalent glasses[8].

Tables

Table 1. Table 1: Table of small force exponents ( θ 𝜃 \theta , θ b subscript 𝜃 𝑏 \theta_{b} ) and g ​ ( r ) 𝑔 𝑟 g(r) power law exponents ( γ g subscript 𝛾 𝑔 \gamma_{g} ). The stability criteria require θ > 1 γ g − 2 𝜃 1 subscript 𝛾 𝑔 2 \theta>\frac{1}{\gamma_{g}}-2 , θ b > 1 − 2 ​ γ g subscript 𝜃 𝑏 1 2 subscript 𝛾 𝑔 \theta_{b}>1-2\gamma_{g} . The right hand sides of the inequalities are also shown.

$\varphi$	${\gamma }_g$	$\frac{1}{{\gamma }_g}-2$	$\theta$	$1-2{\gamma }_g$	${\theta }_b$
$0.55$	$0.4212$	$0.374$	–	$0.1576$	–
$0.56$	$0.4232$	$0.363$	$0.404$	$0.1536$	–
$0.57$	$0.4288$	$0.332$	$0.925$	$0.1424$	$0.366$
$0.58$	$0.4341$	$0.304$	$1.199$	$0.1328$	$0.435$
$0.59$	$0.443$	$0.2573$	$1.3654$	$0.114$	$0.571$
$0.61$	$0.4574$	$0.1863$	$1.7242$	$0.0852$	$0.773$
$0.627$	$0.46335$	$0.158$	$1.9782$	$0.0733$	$1.164$

Equations

$$\overrightarrow{f_i}={\Sigma }_j({\hat{n}}_{ij}{f}_{ij}^n+{\hat{\theta }}_{ij}{f}_{ij}^{\theta }+{\hat{\varphi }}_{ij}{f}_{ij}^{\varphi })=0$$

$$\overrightarrow{{\Gamma }_i}={\Sigma }_j\overrightarrow{R_i}\times \overrightarrow{f_{ij}}=0$$

$$\overrightarrow{{\Gamma }_i}={\Sigma }_j{R}_i(f_{ij}^{\theta }{\hat{\varphi }}_{ij}-f_{ij}^{\varphi }{\hat{\theta }}_{ij})$$

$$M=\left(\begin{array}[]{cccc}n_{12}{}^{x}&\theta_{12}{}^{x}&\phi_{12}{}^{x}&\cdots\\ -n_{12}{}^{x}&-\theta_{12}{}^{x}&-\phi_{12}{}^{x}&\cdots\\ n_{12}{}^{y}&\theta_{12}{}^{y}&\phi_{12}{}^{y}&\cdots\\ -n_{12}{}^{y}&-\theta_{12}{}^{y}&-\phi_{12}{}^{y}&\cdots\\ n_{12}{}^{z}&\theta_{12}{}^{z}&\phi_{12}{}^{z}&\cdots\\ -n_{12}{}^{z}&-\theta_{12}{}^{z}&-\phi_{12}{}^{z}&\cdots\\ \vdots&\vdots&\vdots&\vdots\\ 0&R_{1}\phi_{12}{}^{x}&-R_{1}\theta_{12}{}^{x}&\cdots\\ 0&R_{2}\phi_{12}{}^{x}&-R_{2}\theta_{12}{}^{x}&\cdots\\ 0&R_{1}\phi_{12}{}^{y}&-R_{1}\theta_{12}{}^{y}&\cdots\\ 0&R_{2}\phi_{12}{}^{y}&-R_{2}\theta_{12}{}^{y}&\cdots\\ 0&R_{1}\phi_{12}{}^{z}&-R_{1}\theta_{12}{}^{z}&\cdots\\ 0&R_{2}\phi_{12}{}^{z}&-R_{2}\theta_{12}{}^{z}&\cdots\\ \vdots&\vdots&\vdots&\vdots\\ \end{array}\right)$$

$$M\mid F⟩=0$$

$$f_i:{\Sigma }_{j=1}^{D_{ns}}{X}_{ij}{x}_j-y_i=0,y_i\ge 0$$

$$M^{{}^{\prime }}\mid xy⟩=0,$$

$$M^{{}^{\prime}}=\left(\begin{array}[]{ccccccc}X_{11}&X_{12}&\cdots&X_{1k}&-1&0&\cdots\\ X_{21}&X_{22}&\cdots&X_{2k}&0&-1&\cdots\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ X_{C1}&X_{C2}&\cdots&X_{Ck}&\cdots&0&-1\\ \end{array}\right)$$

$$\mid xy\rangle=\left(\begin{array}[]{c}x_{1}\\ :\\ :\\ x_{k}\\ y_{1}\\ :\\ :\\ y_{C}\\ \end{array}\right)$$

$$f_t-\mu f_n\le 0$$

$$-f_t-\mu f_n\le 0$$

$$\Vert f_t\Vert \le \mu f_n$$

$$(f_{\theta }^2+f_{\varphi }^2)\le {\mu }^2{f}_n^2$$

$$f_{\varphi }^2\le (\mu f_n-f_{\theta })(\mu f_n+f_{\theta })$$

$$\mid f_{\theta }\mid \le \mu f_n$$

$$\mid f_{\varphi }\mid \le (\mu f_n-f_{\theta })$$

$$\mid f_{\varphi }\mid \le (\mu f_n+f_{\theta })$$

$$\hat{\sigma }=\frac{1}{V}{\Sigma }_{i\ne j}{\vec{r}}_{ij}\otimes {\vec{f}}_{ij}$$

$$SA=(P_1-P_3)/(P_1+P_2+P_3)$$

$$C_{\alpha \beta }^{(\sigma )}=\text{textlangle}\sum _{k\ne 0}{\sigma }_{\alpha \beta }^{(0)}{\sigma }_{\alpha \beta }^{(k)}\delta ([X^{(0)}-X^{(k)}]-X)\text{textrangle}$$

$$\hat{R}=\frac{1}{N}{\Sigma }_{i\ne j}\frac{\overrightarrow{r_{ij}}}{\mid r_{ij}\mid }\otimes \frac{\overrightarrow{r_{ij}}}{\mid r_{ij}\mid }$$

$$\hat{\sigma }=\frac{1}{V}{\Sigma }_{i\ne j}\overrightarrow{r_{ij}}\otimes \overrightarrow{f_{ij}},$$

$$FA=(C_2-C_1)/(C_1+C_2)$$

$$SA=(P_2-P_1)/(P_1+P_2)$$

Full text

The formation of self-organised structures that resist shear deformation have been discussed in the context of shear jamming and thickeningbi2011 ; mari2015 ; peters2016 , with frictional forces playing a key role. However, shear induces geometric features necessary for jamming even in frictionless packingsvinu2016 . We analyse conditions for jamming in such assemblies by solving force and torque balance conditions for their contact geometry. We demonstrate, and validate with frictional simulations, that the isostatic condition for mean contact number $Z=D+1$ (for spatial dimension $D=2$, $3$) holds at jamming for both finite and infinite friction, above the random loose packing density. We show that the shear jamming threshold satisfies the marginal stability condition recently proposed for jamming in frictionless systemswyart2012 . We perform rigidity percolation analysisthorpe1983 ; henkes2016 for $D=2$ and find that rigidity percolation precedes shear jamming, which however coincides with the percolation of over-constrained regions, leading to the identification of an intermediate phase analogous to that observed in covalent glassesintermediate .

Force networks and jamming in shear deformed sphere packings

H. A. Vinutha

Jawaharlal Nehru Center for Advanced Scientific Research, Jakkur Campus, Bengaluru 560064, India.

TIFR Center for Interdisciplinary Sciences, 21 Brundavan Colony, Narsingi, Hyderabad 500075, India.

  

Srikanth Sastry

Jawaharlal Nehru Center for Advanced Scientific Research, Jakkur Campus, Bengaluru 560064, India.

Jamming is the process by which disordered assemblies of particles become rigid and resist externally imposed stresses, for instance when their density becomes large enough. It has been widely investigated, both as a phenomenon that occurs in granular matter, and as a particular aspect of the emergence of rigidity in disordered matter, e. g. colloidal suspensions, foams, glass formers and gels and to understand the rheological properties of thermal and athermal driven systems liu-2010 ; zhang2010 ; bi2011 ; wyart2012 ; mari2015 ; peters2016 . The jamming of frictionless sphere assemblies is particularly well studied and occurs at a packing fraction of $\phi_{J}\approx 0.64$, referred to as random close packing (RCP) or the jamming point. In the presence of friction, jamming is expected to occur down to a significantly lower density, which is $\sim 0.54$ (in $3D$) silbert2010 ; makse2008 in the isotropic case, also known as the random loose packing density (RLP), but strong dependences on friction and protocol lead to a wide range of estimates of this density, [$0.54-0.61$]. A rather different scenario was envisaged by Cates et al cates1998 for jamming in systems subjected to external stress, in which the application of external stress itself leads to a self organisation of particles that could resist stress, and thus lead to jamming. Such a scenario of shear jamming has been studied recently experimentally and theoretically zhang2010 ; bi2011 ; sarkar2013 ; sarkar2016 for sheared granular packings in the presence of friction, but also in the case of frictionless spheres bertrand2016 ; kumar2016 ; baity2016 ; urbani-2017 . However, our understanding is yet incomplete concerning various central issues, such as: (i) the range of densities over which shear jamming may occur, and the corresponding conditions, (ii) the differences and similarities between shear jamming and the isotropic frictionless as well as frictional jamming, (iii) a geometric description of the self organization of particles that lead to jamming behaviour, and (iv) the origins of the geometric organisation observed.

We address these issues in the present work, by exploiting the observation vinu2016 that athermally sheared frictionless sphere assemblies develop structural features that correspond to shear jamming, when frictional forces are present to provide mechanical stability. Thus, the geometry of such assemblies, encapsulated in their contact network can be interrogated to understand geometric and mechanical conditions necessary for shear jamming, avoiding some of the ambiguities attendant simulation and experimental investigations of frictional packings. We do so by solving force and torque balance conditions that must be satisfied, with finite contact forces, by jammed states, for both finite and infinite friction. We employ a new method which improves the accuracy of such computations, and validate our results through simulations of frictional particles. Details of our computations are provided in the Methods section and in the Suppmentary Information (SI). As described below, we obtain a precise, detailed geometric characterisation of the shear jamming transition that is at variance in key aspects with results for isotropic frictional jamming makse2008 , and makes contact with analyses of rigidity in the rather different context of covalent network glasses intermediate . We also obtain a mechanical characterisation, that is consistent with criteria for marginal stability analysed for frictionless packings but not hitherto applied to shear or frictional jamming.

We first consider three dimensional sphere assemblies that are athermally sheared, and consider force and torque balance conditions as a function of strain to estimate the jamming strain, in the limit of infinite friction (friction coefficient $\mu\rightarrow\infty$). Fig. 1(a) shows the mean contact number $Z$ for a range of densities from $0.56$ to $0.627$. Fig. 1(b) and (c) show respectively the shear stress $\sigma_{xz}$, where $xz$ is the shear plane and the average contact force $\langle f_{c}\rangle$ obtained from the solutions to the balance conditions, where averages are performed over contacts and over independent solutions. The shear stress and the average contact force exhibit sharp increases at density dependent strain values. We identify the shear jamming transition from the discontinuous change in $\langle f_{c}\rangle$, which closely corresponds to the strain values where stress shows a sharp increase. The corresponding strain value (vertical lines) is termed the jamming strain. The jamming strain at each density closely corresponds to a strain value at which $Z=D+1$ (as shown in Fig. 1(a)) which is the isostatic value in the infinite friction limit, based on the constraining counting argument due to Maxwell maxwell1864 . In Fig. 1(d), we show the jamming strain values along with the strain values at which $Z=D+1$. We next consider the finite friction case, with $\mu=1$. In Fig. 2(a)(b), we show $Z$ and $<f_{c}>$ as a function of strain for each density, obtained from solutions to the balance conditions (marked BCM $\mu=1$) and from frictional simulations (marked DEM $\mu=1$). Note that in the case of BCM, only contacts with finite contact forces (with threshold $f_{c}<10^{-10}$) are counted, whereas with DEM, contacts are lost during the frictional simulations. As shown in Fig. 1(d), It is clear that for finite friction too, jamming strain values correspond to $Z=4(=D+1)$, at all densities (however, the jamming strain values differ for BCM and DEM, a feature that may depend on protocol details of DEM and need further investigation). This result for shear jamming is at variance with arguments and results for isotropic frictional jamming, wherein $Z$ varies continuously from $6(2D)$ to $4(D+1)$ as $\mu$ is varied from [math] to $\infty$ silbert2010 ; makse2008 .

We characterize the ensemble of independent force solutions by measuring the mean angle $\alpha_{mean}$ between pairs of independent solutions, In Fig. 3(a), we show for two different system sizes the mean angle as a function of density, for $\mu\rightarrow\infty$. $\alpha_{mean}$ decreases with a decrease in density, approaching zero as the lower limit density limit $\phi=0.55$ is approached, indicating that decrease of the solution space volume over and above the decrease in its dimensionality (see Methods). The rapid decrease below $\phi=0.58$ also helps explain the difficulty of finding force balanced configurations at lower densities, as noted in previous work vinu2016 ; vinu-sjperc . Another characteristic feature of forces that signals shear jamming is the saturation of the spatial anisotropy of stresses. In Fig. 3(b), we show the stress anisotropy (defined in Methods) as a function of strain for different densities. Initially we observe a linear increase in the stress anisotropy with strain, which flattens out above a strain value which closely corresponds to $Z=D+1$, indicated by dashed lines for each density.

In Fig. 4(a), we show distribution of contact forces (whose magnitudes include normal and tangential components) for different densities in the steady state. These distributions display peaks at finite forces, indicative of jamming ohern-2002 , above the density $\phi=0.55$ which thus marks the lower density limit to jamming. Recently, attention has been focussed on the distribution of forces at the small force limit, observed to obey a power law distribution $P(f)\sim f^{\theta}$, which we consider next. The exponent $\theta$, together with the power law exponent $\gamma_{g}$ governing the near contact singularity of the pair correlation function $g(r)$ have been related through a stability criterion by Wyart and co-workers wyart2012 ; lern2013 . Considering only extended mode instabilities, the inequality expressing the criterion for stability is $\gamma_{g}\geq\frac{1}{2+\theta}$, whereas considering local buckling modes, the stability criterion is $\gamma_{g}\geq\frac{1-\theta_{b}}{2}$. With a view of studying the extent to which these stability criteria may correspond to the jamming thresholds we observe, we compute the small force distributions, for steady state strains, which are shown in Fig. 4(c). We observe a regime in small forces which is described by a power law, whose slope decreases as density decreases. The behaviour and significance of the distributions at smaller forces than the power law regime is difficult to analyse with confidence, which we discuss briefly in the SI. In Fig. 4(d), we show the small force exponent values as a function of density. We observe that the inequality $\gamma_{g}\geq\frac{1}{2+\theta}$ is observed at higher densities, and approaches an equality as $\phi=0.55$ is approached, marking it as the lower density limit to shear jamming. We monitor the evolution of force distributions as a function of strain (approaching the jamming strain from above) to investigate whether the jamming strain corresponds to marginal stability. In Fig. 4(e)(f), we show small force distributions (evaluated over bins of strain value for statistics) and the small force exponent as a function of strain, for $\phi=0.61$. The strain value above which the stability criterion is met is the jamming strain identified earlier, as we show for three densities Fig. 4(b). Thus, the marginal stability criterion, analysed for frictionless systems wyart2012 also describes the shear jamming threshold. This remarkable agreement is not a priori obvious, and prompts theoretical analysis of the stability criterion for frictional systems under the application of shear deformation.

The results above characterise the force networks we obtain and the shear jamming limit density of $\phi=0.55$. In order to elucidate further the nature of the force networks we generate, we perform a rigidity analysis, as we describe next. Rigidity percolation analysis thorpe1983 ; jacobs1995 ; mouka-1995 has previously been used to study the rigidity of covalent glass networks and jammed packings with and without friction henkes2016 . We perform this analysis, along with force balance analysis, for a two dimensional soft disc system, to take advantage of reliable methods for rigidity analysis for two dimensional constraint networks. The shear jamming transition for AQS configurations, for $\mu\rightarrow\infty$ (BCM), occurs at strain values very close to $Z=D+1$, which is consistent with the $3D$ results (see also zhang2010 ), as shown in Fig. 5(a). Rigidity analysis is performed using the pebble game jacobs1997 ; henkes2016 , as described in the Methods section. In this analysis, jamming corresponds to the emergence of a system spanning rigid cluster. In Fig. 5(a), we show that percolation of rigid clusters, along all directions, (marked rigid) occurs before the isostaticity condition, $Z=(D+1)$, is reached, as has also been observed for sheared frictional packings in henkes2016 . To understand this discrepancy, we also consider percolation of over-constrained (marked rigid+stress) regions in our configurations. It has been noted by Moukarzel et al. mouka-1995 that the onset of stress transmission through a lattice of springs occurs when stressed clusters of macroscopic size are present. In Fig. 5(a), we observe that strain values of rigid+stress percolation closely correspond to the shear jamming transition. We perform a system size analysis of percolation probabilities, as a function of contact number $Z$, and find threshold values of $Z\approx 2.89$ for rigid, and $Z\approx 3.0$ for rigid+stress percolation, see Fig. 5(b). Thus, shear jamming corresponds to rigid+stress percolation, preceded by an intermediate regime of rigidity percolation without stress propagation or shear jamming. This observation has not previously been reported, although results that suggest such an intermediate regime have been reported for compressed granular packings bandi . The presence of an intermediate phase was previously observed in chalcogenides and oxide glasses intermediate ; thorpe2000 ; chuby2006 ; briere2007 . In the case of shear jamming, the role and implications of the intermediate phase is not clear and merit further investigation. An appealing possibility is that fragile force networks form in the intermediate phase, see SI Fig. $S13$ and $S14$. To see the relation between force balance conditions and rigidity percolation analysis, we overlay the strong force network $f_{c}>\langle f_{c}\rangle$ onto the network of rigid and floppy regions, and observe that contacts with strong forces are mostly concentrated on discs that belong to the over-constrained regions, see Fig. 5(c) and SI Fig. $S15$. We further characterize the nature of the shear jamming transition by the cluster size distribution of rigid+stress percolation, see Fig. 5(d), which show features characteristic of a continuous percolation transition.

In summary, we develop a new approach to solving force balance conditions which improves the accuracy of such calculations, and show that sheared frictionless spheres evolve self-organized structures that can support external stress, and be jammed, should frictional forces also be present. We show convincingly that the random loose packing density of $0.55$ is the low density limit of shear jamming. We find that the mean contact number required for shear jamming is $Z=D+1$, the isostaticity condition for frictional particles, independent of friction coefficient and shear jamming protocols. This result appears to be valid for shear jamming, but not in general for frictional jamming makse2008 . We also show that the stability criterion proposed by Wyart wyart2012 is valid for the shear jammed states we investigate. This is particularly interesting as the role of near contact neighbours in the analysis of Wyart et al. wyart2012 ; lern2013 also appears to hold for shear jammed states. We compare our analysis of force balance conditions with rigidity percolation analysis and find that rigidity percolation precedes shear jamming in strain at any given density. More interestingly, we show that shear jamming corresponds to the percolation of over-constrained regions, implying also the presence of an intermediate phase in shear jamming systems, analogously to the case of covalent glass formers. Our results thus reveal many interesting geometric aspects of shear jamming in sharp detail, although the implications of some of these features require further investigation to elucidate.

Methods:

The starting point for our analysis is the specification of the contact network for athermally sheared sphere and disc configurations, which we generate using shear deformation of frictionless mono-disperse soft spheres (and binary soft discs, specified below) using the athermal quastistatic shearing (AQS) protocol, as detailed in previous work vinu2016 ; vinu-sjperc , which we refer to for further details. The contact networks obtained from sheared frictionless soft spheres/discs are used to obtain force balance solutions. We study sphere packings for a wide range of densities [$0.55-0.63$] and system sizes $N=256,2000$. The number of samples used to obtain Fig. 1 and Fig. 2 is between $5$ and $10$. For Fig. 3(a), the data for $N=256$ are averaged over all independent solutions of $5$ steady state (SS) configurations (by which we mean the large strain regime in which the properties of the packings, e. g. the mean contact number, do not statistically change with changes in strain) and for $N=2000$ independent solutions are obtained from $6$ configurations. The number of samples used to obtain Fig. 3(b) is above $5$. The number of samples used to obtain Fig. 4(a),(c) is $4000$ force configurations and $20$ contact geometries. Force distributions for each $Z$ window in Fig. 4(e) is obtained by averaging at least $20$ contact geometries and more than $50$ force configurations. In order to perform rigidity analysis, we analyze sheared frictionless packings in 2D of a $50:50$ binary mixture of soft discs, the diameter ratio being $1.4$, for densities $0.79-0.835$ and system sizes ($N=2000,20000,10^{5},10^{6}$). The number of configurations used to obtain percolation data in Fig. 5 are $300,75,10,10$ for $N=2000,20000,10^{5},10^{6}$ respectively.

Force Balance Solutions – The null space method: We develop a new method based on projecting the problem onto the null space of the contact matrix, described below. Various simulation techniques and numerical methods have been employed in previous work to obtain force balanced contact geometries, and force networks for a given contact geometry lerner2012 ; lerner2013 ; gendelman2016 ; jean1999 ; unger2003 ; hurley2016 ; snoeijer2004 . Our method vastly improves the accuracy of the solutions, see SI Fig. $S1$(b). To obtain contact forces for a given contact network, first we express the vector $\vec{r}_{ij}$ joining the center of two spheres $i$ and $j$ in spherical polar coordinates ($\hat{n}_{ij},\hat{\theta}_{ij}$,$\hat{\phi}_{ij}$). Let $\vec{f_{i}}$ and $\vec{\Gamma_{i}}$ denote the total force and the total torque on the $i^{th}$ particle and $\vec{f_{ij}}$ is the force exerted on particle $i$ from particle $j$. Then the force and torque balance conditions are written as follows.

[TABLE]

Where $R_{i}$ is the radius of particle $i$. The matrix $M$ is constructed from the unit vectors between particles in contact. For a single contact ( SI Fig. $S1$(a)) the matrix $M$ is shown below.

[TABLE]

Then, we write down the force balance and torque balance conditions in the matrix form.

[TABLE]

where $M$ is $((\frac{D(D+1)}{2})N\times DC)$ matrix, $D$ is the dimensionality of space, $C$ is the number of contacts and $\mid F\rangle$ is a vector of size $DC\times 1$, with $3$ (for $D=3$, 2 for $D=2$) force components ($f^{n},f^{\theta},f^{\phi}$) for each contact. The normal forces $f^{n}$, which form the first $C$ elements of the matrix, need to be positive, which is a constraint to be imposed on all solutions. Now, we construct a matrix $H=M^{T}M$, which is of dimension $DC\times DC$. Using the matrix $M$, we construct an energy function $E=\langle F\mid M^{T}M\mid F\rangle$ to directly obtain contact forces (Direct Minimization method). But the null space method we describe here offers a more accurate method. To implement this method, we diagonalize the matrix $H$ and obtain eigenvectors with zero eigenvalues, as they satisfy conditions of mechanical equilibrium, which form the basis of null space. Hence any vector (force solution) obtained from a linear combination of these eigenvectors in the null space also satisfies the force balance conditions. Let ${\bf X_{j}}$ represent eigenvector $j$ and has $DC$ elements. The first $C$ elements do not in general satisfy positivity constraints and hence are not physical forces. We obtain physical force solutions in the null space by finding the coefficients $x_{j}$ of the eigenvectors that will satisfy the positivity constraint on $f_{n}$. In other words:

[TABLE]

where $i$ represents the contact number, $D_{ns}$ is the number of eigenvectors in the null space, and $y_{i}$ are auxiliary variables introduced to impose positivity. The above set of equations can be written in a matrix form.

[TABLE]

where the matrix $M^{{}^{\prime}}$ has the dimension $C\times(D_{ns}+C)$, the vector $\mid xy\rangle$ has dimension $(D_{ns}+C)\times 1$ and are given by

[TABLE]

and

[TABLE]

By defining an energy function $E^{{}^{\prime}}=\langle xy\mid M^{{}^{\prime}T}M^{{}^{\prime}}\mid xy\rangle$ which is a quadratic function with positivity constraints on $y_{i}$. The energy functions $E^{{}^{\prime}}$ and $E$ are minimized using reflective Newton method for bound constraint minimization (BCM) coleman1996 , and we use the label BCM to refer to either one of the methods above. We find that independent initial guesses generate independent (but not orthogonal) solutions. We use random forces as well as forces from DEM simulations as initial guesses. The force scale is set by the magnitude of the initial guess. For SS packings, using the null space method, we obtain all the independent solutions for different densities.

We also implement the Coulomb criterion that restricts the magnitude of the tangential forces, in two and three dimensions. We show data for $\mu=1$, a physical value of the friction coefficient, in addition to solutions for the infinite friction case.

In $2D$, the Coulomb criterion, $\mid f_{t}\mid\leq\mu f_{n}$, can be written as a set of linear constraints as follows:

[TABLE]

We include these constraints in the matrix $M$ or $M^{{}^{\prime}}$, by introducing auxiliary variables as above.

In 3D, the Coulomb criterion is a quadratic constraint, but can be expressed as linear constraints as follows. We have

[TABLE]

which can be written as

[TABLE]

The last three equations are implemented in a similar way as the $2D$ case. This increases the dimension of the matrix $M$ or $M^{{}^{\prime}}$ by $2\times(3C)$ and hence it is computationally expensive. The quadratic constraint can be directly imposed and solved using interior point methods developed for solving quadratically constrained quadratic program (QCQP). Solutions for these constrained optimization problems are obtained using the optimization library in MATLAB.

Using the above method, we first obtain force balance solutions using the forces obtained from DEM simulations cundall1979 for the same initial SS configurations. In SI Fig. $S2$, we show that force solutions obtained from the two methods match well.

Stress anisotropy: Stress anisotropy is computed from the eigenvalues of the stress tensor. Stress tensor is defined as follows:

[TABLE]

Where $\vec{r}_{ij}$ is the distance between the centers of spheres $i,j$. We diagonalize the stress tensor and obtain principal components of the stress tensor and the corresponding eigenvalues. Eigenvalues of the stress tensor are $P_{1}>P_{2}>P_{3}$, where $P_{1}$ is along the compressive direction, $P_{2}$ is along the transverse or vorticity direction, and $P_{3}$ is along the dilative direction. The stress anisotropy is defined as follows:

[TABLE]

The pebble game algorithm: We implement a ($k=3,l=3$) pebble game jacobs1997 ; henkes2016 , which we describe briefly. The algorithm is based on Laman’s theorem, which states that the network of $N$ vertices is generically, minimally rigid in two dimensions if and only if it has $3N-3$ bonds and no subgraph of $n$ vertices has more than $3n-3$ bonds laman1970 . The pebble game is implemented as follows: Each disc (from now on called as a site) is assigned $k=3$ pebbles, each pebble corresponding to one degree of freedom, $2$ translational and $1$ rotational. The quantity $l=3$ corresponds to the total number of global translational and rotational degrees of freedom that are always present for a rigid body and must be accounted. Each contact represents one constraint for translational motion and one constraint for rotational motion. Hence, each contact in the contact network is replaced by two bonds, representing the two constraints. A series of steps are used to assign pebbles to bonds, such that each bond that is covered by a pebble is an independent constraint restricting one degree of freedom. When a pebble is assigned to a bond, it is assigned a direction, If the bond $E_{ab}$ is covered by a pebble from site $a$, then the bond is directed from $a$ to $b$. Pebbles that remain on the sites are free and can be used to cover other bonds, but once a bond is covered by a pebble, it continues to be covered, and a pebble covering a bond can be moved only with another pebble taking its place. A bond can be covered by a free pebble from one of the two sites at either end of it, if the total number of pebbles at its sites is $(l+1)$. If for a bond the total number of pebbles at its sites is less than $l+1$, we search for a free pebble to cover this bond, but the search is conducted only along bonds that are directed away from the sites adjacent to the bond (ı. e., we attempt to retract a pebble previously assigned to a bond, and replace it with another from elsewhere). From each site $v$, we can search along ($k=3$) bonds directed away from $v$. The search for a free pebble continues until a pebble is found and a sequence of swaps allows the bond under consideration to be covered and marked as an independent bond. If a free pebble is not found, due to the search process encountering a set of closed loops that takes the search back to the initial sites, then the bond is marked as redundant. The algorithm terminates when each bond is marked as an independent bond or a redundant bond. The presence of redundant bonds leads to over-constrained regions. The sites visited during the failed pebble search belongs to the over-constrained (rigid +stress) regions.

We map out rigid clusters from the network of redundant and independent bonds. Only the independent bonds are used to identify rigid clusters, and we seek to label them so that all the independent bonds belonging to a rigid cluster has the same cluster label. We start with an unlabelled bond $E_{ab}$ and we perform a pebble search to obtain $k=3$ pebbles (which can always be found) and pin the three pebbles at the sites $a$ and $b$ (i.e., these three pebbles are not free). Sites $a$ and $b$ are marked rigid and bond $E_{ab}$ is assigned a cluster label. Next, we inspect neighbours of sites $a$ and $b$ to mark them floppy or rigid with respect to $a$ and $b$. For $a_{1}$, a neighbour of site $a$, we perform a pebble search and attempt to free a pebble for the bond $E_{a,a_{1}}$. If a free pebble is found then site $a_{1}$ is marked floppy. If a free pebble is not found then the site $a_{1}$ is marked rigid and also all the sites that are visited during the failed search are marked as rigid with respect to the initial sites. The above procedure is repeated until all the neighbours of the rigid sites are marked floppy. All the bonds between pairs of sites marked rigid are given the same cluster label as $E_{ab}$. The procedure is repeated until all the bonds are assigned a cluster label, with the labelling as rigid or floppy being removed when a new cluster search is initiated each time. Clusters containing only two sites (one bond) belong to floppy regions. More details of the algorithm are in the reference jacobs1997 .

After the pebble game, we have sites (or discs) belonging to the biggest rigid cluster and also sites that are identified as over-constrained, which we then test for percolation. The percolation probability is computed by first identifying the biggest cluster by connecting particles that are in contact using Lees-Edwards periodic boundary conditions. Then we check if the biggest cluster percolates not only in the simulation box (by inquiring whether it extends from one edge of the box to the other, as done in previous studied vinu-sjperc ), but also considering an extended system composed of the simulation box, surrounded by periodic copies of the system, and testing whether the largest cluster percolates across the extended area so defined. These two procedures do not produce significantly different results in cases studied before vinu2016 ; vinu-sjperc , but we find that for rigidity percolation analysis, the thresholds depend more significantly on the definition used. We thus use the more restrictive definition of percolation across the extended area. Since the force network in the sheared systems we study are anisotropic, we study percolation independently along one direction ($x$ or $y$ or compressive (diagonal)) and all directions ($x$ and $y$ or compressive and dilative directions).

To obtain cluster size distributions, shown in Fig. 5, we first compute the percolation threshold value of the contact number $Z_{c}$ seperatly for each configuration as the value where where the rigid+stress particles percolate i.e., cluster spans at least in one direction. The average value of $\langle Z_{c}\rangle=2.858$, is used in the plot shown.

Force networks and jamming in shear deformed sphere packings(Supplementary Information)

H. A. Vinutha and Srikanth Sastry

Here we present additional information regarding various aspects of our analysis of force balance solutions, namely: (i) Accuracy of the null space method, (ii) Comparison with DEM forces, (iii) Dimensionality of the null space, (iv) Small force exponents and stability criteria, (v) Small force distributions as a function of strain, (vi) Spatial correlations of stress, (vii) Density-Strain phase diagram for the two dimensional system, (viii) Rigid and rigid+stress percolation in two dimensions, (ix) Intermediate phase in shear jamming and (x) Evolution of rigid, rigid+stress, and force networks.

I Accuracy of the Null Space Method

In Fig. S1, we show the average contact force $\langle f_{c}\rangle$ and the average of total force on a particle $\langle F_{t}\rangle$, for SS configurations at different densities, for force solutions that are obtained using the null space method and from direct minimization. The magnitude of the average contact force is set by the initial guess used by the minimization protocol. Figure S1 shows that force solutions obtained using the null space method are more accurate and are force balanced up to an accuracy of $\approx 10^{-14}$, whereas the force solutions obtained from direct minimization are force balanced up to an accuracy of $\approx 10^{-11}$.

II Comparison with DEM forces

Using the BCM method, we first obtain force balance solutions and compare the solutions to the forces obtained from DEM simulations for the same initial SS configurations, obtained quasistatically with a small strain step to ensure there is a minimal change in the contact network during DEM vinu2016 .The initial guess for forces is taken from DEM simulations, and the simulation details and DEM parameters are as mentioned in vinu-sjperc . The normal and tangential components of the contact forces, $f_{n}$ and $f_{t}=\sqrt{(}(f^{\theta})^{2}+(f^{\phi})^{2})$, are computed for densities $\phi=0.627,0.61$, and $0.58$, including in the matrix $M$ only those contacts that remain at the end of the DEM simulations. The forces at each contact estimated using BCM are compared with those obtained from DEM simulations in Fig. S3, demonstrating that they agree very well with each other. This analysis supports two conclusions: (i) Forces obtained by solving force balance conditions for SS configurations agree well with DEM forces, implying that the SS contact network can support the forces observed in the DEM simulations, which in turn correspond to shear jamming. (ii) The contact networks at the end of the DEM simulations, which support finite contact forces and stresses, are those that are present in the SS configurations. The only changes in the contact network is that some of the SS contacts are lost during the DEM simulations, and we can account for all of them. They are either contacts of rattlers, or contacts that require higher friction coefficients to be retained than what we employ during the DEM simulations.

III Dimensionality of the Null Space

In Fig. S3, we show the null space dimension as a function of density for $N=256$ and $2000$, averaged over $5$ initial SS configurations. Observe that the limiting density when the number of force solution is one is close to $\phi=0.55$, indicated by the green fit curve for the $D_{n}s/6N$ ($N=2000$) data. This is consistent with other indicators discussed in the paper.

IV Small Force Exponents and Stability Criteria

In Fig. $1$(d) of the paper, we showed force distributions obeying a power law distribution $P(f)\sim f^{\theta}$. The exponent $\theta$ and the power law exponent $\gamma_{g}$ of the near contact singularity of the pair correlation function $g(r)$ have been related through the stability criterion $\theta>\frac{1}{\gamma_{g}}-2$. In Fig. $3$(b) of the paper only the variation of small force exponent $\theta$ is shown, which corresponds to extended mode instabilities. In addition to extended modes, a small force distribution arising from buckling modes has also been discussed. We find in our force distributions that the smallest force regime is characterised by a smaller exponent, which we tentatively identify as the buckling mode exponent $\theta_{b}$, which is obtained by a fit to small contact forces (the fit curves in maroon, see Fig. S4 (a)). In Fig. S4, we show the variation of the small force exponent $\theta_{b}$, as the lower density limit is approached. Similar to the $\theta$ exponent, $\theta_{b}$ also shows that as the lower density limit of $\phi=0.55$ is reached the system ceases to obey the stability criterion. In table $1$, we tabulate different exponent values corresponding to extended and buckling modes and the exponent values obtained from relations of marginal stability. Note that we do not have data to directly substantiate the association of forces as due to extended or buckling modes.

V Small force distributions as a function of strain

In Fig. $4$(c) of the paper, we showed force distributions obeying a power law distribution $P(f)\sim f^{\theta}$, as the jamming strain is approached. In Fig. S5, we show $g(r)$ having a power law form near contact for different $Z$ (or $\gamma$) windows. Using the power law exponents of $g(r)$ ($\gamma_{g}$), we compute the force distribution exponent required for stability against extended modes. In Fig. S6 and S7, we show small force distributions and its exponents for $\phi=0.57$ and $\phi=0.58$. As we approach the shear jamming strain, the system loses its stability. Which also occurs when the low density limit is approached.

VI Spatial Correlations of Stress

The presence of non-zero shear stress and non-zero spatial stress correlation of bonds in the shear plane, even in the thermodynamic limit, distinguishes shear jammed packings from isotropically jammed packings, baity2016 . We compute spatial correlation of stresses in the shear plane ($xz$) for the SS configurations at different densities. The correlation is defined as follows:

[TABLE]

where $\sigma^{k}$ is the stress on the bond $k$ and $X^{k}$ is the position of bond $k$. Along the shear plane $xz$, we observe that the correlation are finite and is zero in the other planes. Along the shear plane $xz$, we observe that at large $X$ the correlations are finite for densities above $\phi=0.55$, see Fig. S8(a). In Fig. S8(b), we show spatial stress correlation of bonds along different planes for $\phi=0.61$, which are finite at large $X$ for the shear plane $xz$ and zero in the other planes.

VII Density-Strain Phase Diagram for the Two Dimensional System

In this section, supporting plots for the 2D phase diagram discussed in the main text. In Fig. S9, we show supporting data for the BCM results shown in Fig. $5(a)$ of the paper. The average contact force and stress obtained from the BCM method are shown. The data is averaged over $5$ configurations. For $\phi=0.8$ average is done over $15$ configurations. The BCM jamming strain is obtained from the strain values marked by the vertical lines in plot of $<f_{c}>$, which is very close to $Z=D+1$. In Fig. S10, contact force data for the finite friction coefficient $\mu=1$ from BCM is shown along with DEM data for $\mu=1$, showing that the jamming strain values for BCM and DEM are close to each other and to $Z=3(=D+1)$.

VIII Rigid and Rigid+Stress Percolation in Two Dimensions

In Fig. S11, we show rigid and over-constrained (rigid+stress) percolation probabilities as a function of strain, for $N=2000$. The rigid and rigid+stress strain values in Fig. $5$(a) of the paper corresponds to percolation probability of $0.5$ of case all (where we require percolation to occur along both axes). The percolation along one and all directions occurs at different strain values. This separation is more apparent as a function of $Z$ rather than $\gamma$, as shown in Fig. S13. To show the presence of an intermediate phase for different densities, we show the percolation of rigid clusters and over-constrained regions (rigid+stress) as a function of $Z$ (for the all case) in Fig. S12. Since the configurations are anisotropic, the rigid and rigid+stress percolations along one and all directions are well separated, as shown in Fig. S13. The percolation of rigid+stress clusters that percolate along one direction can support stress transmission along the direction of percolation. Propagation of stress along one direction is characteristic of the fragile force networks.

IX Intermediate phase in shear jamming

The possibility of fragile force networks to form in the intermediate phase is supported by two results below. In Fig. S13, we show rigid+stress percolation along one and all directions to occur at different values of $Z$ and $\gamma$ and the rigid+stress percolation along one direction occurs after the rigid percolation. In the intermediate window, the percolating clusters of over-constrained regions can support strong forces along one direction, which is the definition of fragile force network cates1998 ; bi2011 . We characterize the anisotropy in the sheared two dimensional packings using fabric and stress anisotropy, defined as follows:

[TABLE]

where $\vec{r_{ij}}$ is the distance between the centers of pair $(ij)$. We diagonalize the fabric (stress) tensor and obtain principal components of fabric (stress) tensor and the corresponding eigenvalues. Eigenvalues of the fabric tensor are $C_{1}>C_{2}$ and the stress tensor are $P_{1}>P_{2}$. $C_{1}$ ($P_{1}$) is along the compressive direction and $C_{2}$($P_{2}$)is along dilative direction. The stress and fabric anisotropy are defined as follows:

[TABLE]

In Fig. S14, we show the stress anisotropy and fabric anisotropy as a function of strain for different densities. In the stress anisotropy plot, we mark the strain values corresponding to $Z=3$ and rigid percolation. We see that in this window of strain values the anisotropy of the force network starts to deviate from the linear behavior and to saturate, which supports the above idea. These ideas are speculative at the moment and require further investigation to validate them.

X Evolution of Rigid, Rigid+Stress, and Force Networks

In Fig. S15, we show the evolution of rigid clusters and the force network as a function of strain, for one initial configuration strained using the AQS protocol at $\phi=0.82$. The snapshots clearly show that there is a strong spatial correlation between the rigid cluster (red network, left panel) and the force network (right panel). Particularly, contacts with forces $f_{c}>\langle f_{c}\rangle$ are concentrated on the discs that belong to the over-constrained (or stressed) regions, the sites of which are marked with blue dots. Observe that the average contact force increases with increase in the number of over-constrained sites. Above the shear jamming transition, almost all the particles, except a few rattlers, are over-constrained and the system can have multiple force balance solutions for the same contact network. Hence the comparison becomes less revealing at larger strains.