Emergent Elasticity in Amorphous Solids
Jishnu N. Nampoothiri, Yinqiao Wang, Kabir Ramola, Jie Zhang, Subhro Bhattacharjee, Bulbul Chakraborty
EEmergent Elasticity in Amorphous Solids
Jishnu N. Nampoothiri , , Yinqiao Wang , Kabir Ramola , Jie Zhang , Subhro Bhattacharjee ,Bulbul Chakraborty Martin Fisher School of Physics, Brandeis University, Waltham, MA 02454 USA Centre for Interdisciplinary Sciences, Tata Institute of Fundamental Research, Hyderabad 500107,India Institute of Natural Sciences and School of Physics and Astronomy, Shanghai Jiao Tong Univer-sity, Shanghai, 200240 China International Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bengaluru560089, India
The mechanical response of naturally abundant amorphous solids such as gels, jammedgrains, and biological tissues are not described by the conventional paradigm of broken sym-metry that defines crystalline elasticity
1, 2 . In contrast, the response of such athermal solidsare governed by local conditions of mechanical equilibrium, i. e., force and torque balance ofits constituents. Here we show that these constraints have the mathematical structure of elec-tromagnetism, where the electrostatic limit successfully captures the anisotropic elasticity ofamorphous solids . The emergence of elasticity from constraints offers a new paradigmfor systems with no broken symmetry, analogous to emergent gauge theories of quantumspin liquids . Specifically, our U (1) rank-2 symmetric tensor gauge theory of elasticitytranslates to the electromagnetism of fractonic phases of matter with stress mapped to elec- a r X i v : . [ c ond - m a t . d i s - nn ] A p r ric displacement and forces to vector charges. We present experimental evidence indicatingthat force chains in granular media are sub-dimensional excitations of amorphous elasticitysimilar to fractons in quantum spin liquids . Solids that emerge in strongly nonequilibrium processes such as jamming or gela-tion
17, 18 , are characterized by strong stress heterogeneities, often referred to as force chains. Theyare rigid in that they can sustain external shear, yet they are often fragile
14, 15 . Analogous to clas-sical elasticity theory, it is plausible to ask whether a long wavelength field theoretic descriptionexists for the mechanical response of such athermal solids and if so, what are its characteristicsand universal features, and what are the appropriate variables that can account for the underlyingkinetic constraints in the emergent field theory? Any attempt to construct such a field theory mustanswer: (a) how to obtain the stress state, and (b) how to incorporate microscopic informationabout the structural disorder, accounting for the mechanical constraints into a continuum formu-lation. This second problem, in particular, has a close resemblance with kinetically constrainedmodels such as hard-core dimer models on lattices where the hard-core constraint naturally allowsfor an emergent gauge theory description at low energy and long wavelengths
9, 10, 12 .In this work, we develop a theory of stress transmission in disordered granular solids, bothwith and without friction, where the local constraints of mechanical equilibrium are paramount–every grain satisfies the constraints of force and torque balance. These local constraints imply thatthe grain-level stress tensor ˆ σ g , is symmetric
19, 20 and satisfies ( ∇ · ˆ σ ) g = (cid:88) c ∈ g f g,c = f ext g , (1)2here f g,c are the contact forces acting on grain g and f ext g is the total external force acting on thegrain. The divergence operator above is a discrete version of the differential operator, which isdefined over the underlying contact network, as detailed in the SI. The condition encapsulatingthe local kinetic constraints is the exact analog of the celebrated Gauss’s law in the theory ofelectromagnetism. In granular solids, the external boundary forces give rise to the internal stressfield; analogous to charges generating a non-zero electric field in electromagnetism.The kinetic constraints in Eq. 1 on ˆ σ g is the grain-level realization of the continuum equationof mechanical equilibrium: ∂ i σ ij ( r ) = f j . Following the electromagnetic analogy further, it istempting to identify this condition with the Gauss’s law for a U (1) , symmetric rank-2 tensor gauge-theory with vector charges (VCT) , with the mapping E ij = E ji ? ↔ σ ij and vector charges tounbalanced forces i.e. ρ ↔ f , satisfying ∂ i E ij = ρ j , (2)Eq. 2 correctly incorporates, by construction, both the conservation of charge/force ( (cid:82) d r ρ = 0 ),and charge angular momentum/torque ( (cid:82) d r ( r × ρ ) = 0 )
21, 22 . Incidentally, as pointed out recently,the two conservation laws lead to sub-dimensional propagation of charges– a feature of recentlydiscussed fractonic phases of matter as well as topological defects in elastic solids . In thepresent context, it also provides a natural explanation for the visually striking “force chains” (seeFig. 3) observed in photoelastic images of granular solids , as our analysis will demonstrate.This reformulation of the stress equation in terms of the Gauss’s law for symmetric tensorelectromagnetism gives us a natural starting point for understanding the mechanical response of3ranular solids and a derivation of the correct continuum theory. The above formulation is similarto the problem of frustrated magnets and/or dimer models where due to non-trivial local ener-getic/kinetic constraints, the individual spins/dimers cease to be the right degrees of freedom andhence fail to describe the low energy theory, which in turn is often described by emergent gaugefields that naturally capture the constraints . Similarly, the displacement of the individual grains–the mainstay of the theory of elasticity of crystalline solids– cease to be the right variables in theabsence of broken translation symmetry. However, the long-range stress correlations generated byNewton’s laws of force and torque balance are captured by the emergent electromagnetism.Eq. 2 does not provide enough equations to solve for the field E ij , since there are only d equations in d dimensions, for the d ( d + 1) / components of a symmetric tensor . This lack ofequations is an issue that has plagued the solution of the granular stress problem . In VCT, theGauss’s law constraint generates a gauge transformation for the symmetric tensor gauge potential, A ij → A ij + ( ∂ i ψ j + ∂ j ψ i ) . Since E ij is the momentum conjugate to A ij , we have the analog ofthe Maxwell-Faraday condition: ∂B ij ∂t = − (cid:15) iab (cid:15) jcd ∂ a ∂ c E bd , where B ij is the magnetic field: B ij = (cid:15) iab (cid:15) jcd ∂ a ∂ c A bd . Invoking the magnetostatic condition (cid:15) iab (cid:15) jcd ∂ a ∂ c E bd = 0 , leads to the potentialformulation of electrostatics: E ij = 12 ( ∂ i φ j + ∂ j φ i ) , (3)which can be used to obtain E ij for any charge configuration .The theory described above, although self-consistent, does not capture the full complexityof granular mechanics. Crucial to their properties are the twin aspects of a granular solid: (1) it4s only defined under external pressure (as a packing of grains with purely repulsive interactionswill fall apart in the absence of such boundary forces), and (2) it can support internal stresses. Thistranslates, in terms of the VCT, to an assembly being subject to well defined boundary charges de-veloping internal charge dipoles , akin to the response of a polarizable medium (dielectric). Stateddifferently, once under external pressure, although the granular solid is free of “charges” sinceevery grain satisfies the constraints of force and torque balance, “bound charges” exist as pairs ofequal and opposite forces at every contact of the disordered granular network. We need to includethese force or charge dipoles in our theory through a tensorial dipole moment P ij , whose diver-gence defines the bound charges. The structure of P ij is influenced by various microscopic detailsof the system such as the features of the underlying contact network and the nature of contact forceswhich, for example can be purely repulsive or attractive and frictionless or frictional. To constructa continuum theory, we assert that P ij is related to E ij , through a fourth-rank polarizability tensor,as in linear dielectrics. For such a polarizable medium, Gauss’s law becomes: ∂ i E ij = ρ freej + ρ boundj , (4) ∂ i P ij = − ρ boundj . A complete derivation of these relations, and a detailed discussion of the structure of the theorywill be presented in a future paper. The polarizability tensor, ˆ χ , defined via P ij = χ ijkl E kl , leadsto the definition of a “Displacement” tensor, D ij = ( δ ik δ jl + χ ijkl ) E kl , satisfying ∂ i D ij = ρ freej . ∂ i D ij = 0 ; (cid:15) iab (cid:15) jcd ∂ a ∂ c (Λ D ) bd = 0 ,E ij = ( δ ik δ jl + χ ijkl ) − D kl ≡ Λ ijkl D kl . (5)Since both ˆ D and ˆ E are symmetric tensors, ˆΛ has to satisfy Λ ijkl = Λ jikl = Λ ijlk = Λ jilk . Since the5nherent stresses in a granular solid satisfy the first of the equations above as a direct consequenceof force balance, we need to interpret D ij as the Cauchy stress tensor measured from contact forcesand contact vectors inside the material: σ ij ↔ D ij in Eq. 2. In the presence of body forces such asgravity, the divergence of ˆ D would be equated with the body force. We can compare this theory toanisotropic elasticity : ∂ i σ ij = 0 ; (cid:15) iab (cid:15) jcd ∂ a ∂ c U bd = 0 ,σ ij = (Λ) − ijkl U kl . (6)where U ij is the macroscopic strain tensor. The tensor ˆΛ in Eq. 6 is the inverse of the elasticmodulus tensor. Identifying E ij ↔ U ij , and D ij ↔ σ ij , the VCT-based framework is identicalto the anisotropic elasticity theory with φ , a gauge potential, replacing the displacement field inelasticity theory, and providing the “missing” compatibility equations that allow us to solve thegranular stress response problem without invoking a displacement field.Eq. 5, our main theoretical result, demonstrates that an elasticity theory capturing the fullrange of stress responses in granular solids emerges from VCT. Although ˆ D and ˆ E are identifiedwith ˆ σ and ˆ U , the elasticity emerges from local constraints and not from broken symmetry: it isa stress-only description that needs no reference to a stress-free state or displacement fields. Al-though ˆΛ is an effective elastic modulus, it depends on protocol, and there are no symmetry require-ments imposed by a free-energy. Within the theory, protocols translate to stress-ensembles
3, 4, 19 characterized by externally imposed stresses such as isotropic compression or pure shear.Below, we compare the predictions for stress-stress correlations obtained from Eq. 5 to6xperimental and numerical data. Through this analysis, we extract ˆΛ for frictionless and frictionalgranular solids prepared under different protocols. A hallmark of the VCT in free space, both in 2Dand 3D, is the appearance of pinch point singularities in the Fourier transforms of E ij correlators : C freeijkl ( q ) ≡ (cid:104) E ij ( q ) E kl ( − q ) (cid:105)∝ ( δ ik δ jl + δ il δ jk )2 + q i q j q k q l q − (cid:18) δ ik q j q l q + δ jk q i q l q + δ il q j q k q + δ jl q i q k q (cid:19) . (7)Eq. 7 is obtained by imposing the Gauss’s law constraint, ∂ i E ij = 0 , on Eq. 3, and assumingthat all states are equiprobable , i.e. the Edwards measure . Earlier granular field theories
3, 4, 26 based on this measure used a dual formulation of VCT where the emergence of elasticity is notevident. Since C ijkl ( q ) is independent of | q | , it is straightforward to show that the correlations inreal-space decay as /r d . A stringent test of the theory, therefore is the pinch-point structure of thecorrelation functions.For granular solids, we computed the correlators C ijkl ( q ) = (cid:104) D ij ( q ) D kl ( − q ) (cid:105) using Eq.5, and tested the predictions in ensembles of 2D and 3D isotropically compressed soft particles(numerically), and in ensembles of 2D packings of frictional grains (experimentally). Pinch pointsingularities are clearly exhibited in both 2D (Figs. 1 and 3) and 3D (Fig. 2). We have determined ˆΛ through detailed comparisons between theory and measurements of C ijkl (Figs. 1- 3). Fig. 1demonstrates that for packings created under isotropic compression, ˆΛ = λ , with being theidentity tensor. Additional tests of the theory are presented in the SI.7 ab Figure 1: Comparisons in Fourier space between the theoretical predictions (solid black line) andthe disorder-averaged angular dependent stress-stress correlations C xxxy ( θ ) and C xxyy ( θ ) in thenumerical: a , and the experimental results (red symbols): b , for isotropically jammed systems.The range of pressure for the numerical data is P ∈ [0 . , . and the results are displayedfor five different system sizes N = 512 , , , , . The experimental data is fromfrictional packings with a range of pressure P ∈ [1 . × − , . × − ] . All correlation functionsare normalized by the peak value of C xxxx ( θ ) .To illustrate the sensitivity of the ˆΛ tensor to protocol (stress ensemble) we generated shearedpackings of the same grains used in the isotropic compression. Under the experimental conditionsof pure shear, with principle stress along x and y , a diagonal form with different values of λ ii provides an excellent description of the experimental observations (Fig. 3). We find that λ ii , satisfy8 b c d Figure 2: Comparisons in Fourier space between theoretical predictions (top panels) and numer-ical results (bottom panels) from jammed packings of frictionless spheres in three dimensions.The figures display the radially averaged correlation functions a : C xxxx ( θ, φ ) , b : C xyxy ( θ, φ ) , c : C xxxz ( θ, φ ) and d : C xyyz ( θ, φ ) respectively. The coordinates ( H x , H y ) represent a Hammer projec-tion of the ( θ, φ ) shell onto the plane. The results are presented for system size N = 27000 , andhave been averaged over configurations. The range of packing fractions for these configura-tions is φ ∈ [0 . , . and the range of pressure per grain is P ∈ [0 . , . . Results forthe ˆΛ tensor have not been presented due to the small effective system size: × × . Theblank regions at the poles for the numerical results is due to the difficulty in sampling for distinct φ at θ = 0 and θ = π sheared frictional packings and the theoretical predictions (black line). a : Photoelastic images producedfrom a sheared packing. b : Contour plots of the stress-stress correlation functions: C xxxx ( q, θ ) ( Top ), C yyyy ( q, θ ) ( Middle ) and C xyxy ( q, θ ) ( Bottom ). c : Corresponding angular plots. d : Angularplots of the correlation functions C xxxy ( θ ) ( Top ), C xxyy ( θ ) ( Middle ) and C xyyy ( θ ) ( Bottom ). Allthese correlations are normalized by the peak value of C xxxx ( θ ) . The ˆΛ tensor obtained from thefit is diagonal with λ = 1 , λ = 4 . , and λ = 1 . . The ratio of the boundary stresses is Σ xx / Σ yy = 1 . , which satisfies the bound of positivity
3, 26 λ λ ≥ ( Σ xx Σ yy ) .10 set of bounds imposed by the constraint of positivity of normal forces in granular media
3, 26 .A consequence of the pinch point singularities is that in real-space, (cid:104) D ij ( q ) D ij ( − q ) (cid:105) isnegative in transverse directions and positive along longitudinal directions, as shown in the SI. Itis this property that is strikingly demonstrated in photoelastic images of force chains in Fig. 3 andFigs. S2, S3 and S4 in the SI.To summarize, we demonstrate that elasticity of athermal amorphous solids emerges fromthe local constraints of force and torque balance, and not from broken symmetry. The theoryis an exact analog of the electrostatics of a fractonic U (1) gauge theory in polarizable media.Although our analysis focused on granular solids, it marks a paradigm shift in our understandingof “amorphous elasticity” for a much broader class of solids such as jammed suspensions , andcolloidal gels where thermal fluctuations are completely irrelevant.The current theoretical framework can be extended to thermal amorphous solids such aslow-temperature glass formers . As in frustrated magnets , thermal fluctuations lead to alength-scale characterizing the distance between particles at which force and torque balance areviolated, and wash out the pinch-point singularity, as shown in the SI. A fully dynamical theoryof amorphous materials can be constructed by extending the “electrostatics” to “electrodynamics”through the identification of the analog of a magnetic field, and including unbalanced forces ascharged excitations . 11 ethods The main quantity of interest in this study, for a given packing is the stress tensor field in Fourierspace given by ˆ σ p ( q ) = N pG (cid:88) g =1 ˆ σ pg exp (cid:0) i q · r pg (cid:1) . (8)Here, ‘ p ’ denotes a particular realization or packing of N pG grains, while g denotes a particulargrain in the packing located at r pg . ˆ σ pg represents the force moment tensor for the grain g , given by ˆ σ pg = n gc (cid:88) c =1 r gc ⊗ f gc . (9)Here r gc denotes the position of the contact c from the center of the grain g and f gc denotes theinter-particle force at the contact. Numerical Methods:
We generate jammed packings of frictionless spheres interacting throughone-sided spring potentials in two and three dimensions. Our implementation follows the standardO’Hern protocol
13, 31, 32 , with energy minimization performed using two procedures (i) conjugategradient minimization, and (ii) a FIRE
33, 34 minimization implementation in LAMMPS . Wehave verified that these differences in protocol do not modify our results.We simulate a 50:50 mixture of grains with diameter ratio 1:1.4. In our simulations, thesystem lengths are held fixed at L x = L y = 1 in 2D and L x = L y = L z = 1 in 3D. We imposeperiodic boundary conditions in each direction, setting a lower cutoff between points in Fourierspace q min = 2 π . We choose an upper cutoff q max = π/d min so as to not consider stress fluctua-tions occurring at length scales shorter than d min , the diameter of the smallest grain in the packing.12e have presented data for system sizes N = 512 , , , , in 2D, averaged overatleast configurations for each system size. The results obtained for different system sizes havebeen collapsed (see Fig. 1) using the system size N and q max as scaling parameters. This showsthat the data presented is not significantly affected by finite size effects. All the 2D packings havea pressure per grain P ∈ [0 . , . and packing fraction φ ∈ [0 . , . . In 3D, the datafor N = 27000 is presented in Fig. 2, the data have been averaged over configurations. Therange of packing fractions for these configurations is φ ∈ [0 . , . , with a pressure per grain P ∈ [0 . , . . Experimental Methods:
The experimental results were produced from the analyses of isotrop-ically jammed packings and pure-sheared packings, which were both prepared using a biaxialapparatus whose details can be found in Wang Et al. 2018 . This apparatus mainly consists of arectangular frame mounted on top of a powder-lubricated horizontal glass plate. Each pair of par-allel walls of the rectangular frame can move symmetrically with a motion precision of 0.1 mm sothat the center of mass of the frame remains fixed. To apply isotropic compression, the two pairs ofwalls are programed to move inwards symmetrically. To apply pure shear, one pair of walls movesinwards, and the other pair of walls moves outwards, such that the area of the rectangle is keptfixed. The motion of walls is sufficiently slow to guarantee that the deformation is quasi-static.About . m above the apparatus, there is an array of 2 × φ , which is the ratio between the areaof disks and that of the rectangle. To minimize the potential inhomogeneity of force chains in thejammed packing, we constantly applied mechanical vibrations before the φ exceeded the jammingpoint φ J ≈ . of frictionless particles. The final isotropically jammed packing is confined in asquare domain of 67.2 cm × φ ≈ . , the mean coordination number is around4.1, the pressure is around 12 N/m, and the corresponding dimensionless pressure is × − .Once the isotropically jammed packing was prepared, we then applied pure shear of strain . to the packing to produce the pure-sheared packing. For both types of packings, two differentimages were recorded. Disk positions were obtained using the normal image , recorded withoutpolarizers. Contact forces were analyzed from the force-chain image , recorded with polarizers,using the force-inverse algorithm . 14 upplementary Information In this supplementary document, we describe in detail several key aspects of the theoretical frame-work and analysis of numerical and experimental data. In Section 1, we outline the derivation ofthe Gauss’s law constraint on the Cauchy Stress tensor starting from the constraints of force andtorque balance on every grain and discuss the mapping of grain-level properties to the continuumtheory. In Section 2, we present results for the correlations of the electric displacement tensor ˆ D ,in a polarizable medium characterized by ˆΛ . Further in Section 3, we present experimental datafor stress correlations from individual configurations. Finally, Section 4 describes the numericalresults for the 2D system at finite temperature. Figure S1:
A section of a jammed configuration of soft frictionless disks in 2D. The centers of the grains arelocated at positions r g . The contact points between grains are located at positions r c . The triangle formedby the points r g , r g (cid:48) , r c (shaded area) is uniquely assigned to the contact c and has an associated area a g,c .
5, 7 . Here, we demonstrate theemergence of this Gauss’s law from local constraints of mechanical equilibrium, for the specific ex-ample of disordered granular solids. The arguments can be easily generalized to other amorphouspackings at zero temperature. Granular materials consist of an assembly of grains that interact witheach other via contact forces, as shown in Fig. S1. In a granular solid, each grain is in mechanicalequilibrium and thus, satisfy the constraints of force and torque balance. The constraints of forceand torque balance on a grain g , with no “body forces” can be written as: (cid:88) c ∈ g f g,c = 0 , (cid:88) c ∈ g r g,c × f g,c = 0 , (S1)respectively. Here, f g,c is the contact force, and r g,c the vector joining the center of grain g to thecontact c (Fig. S1). This places dN + d ( d − N nontrivial constraints on the N grains that arepart of the contact network. A grain is said to be a part of the contact network if it has more thanone contact and grains which are not part of the contact network are defined to be “rattlers”. In ourrepresentation, the rattlers become part of voids. Given a set of f g,c and r g,c , one can define a stresstensor for a grain with area A g : ˆ σ g = (1 /A g ) (cid:88) c ∈ g r g,c ⊗ f g,c . (S2)The coarse-grained stress tensor field, ˆ D ( r ) is obtained by summing ˆ σ g over all grains included ina coarse-graining volume, Ω r , centered at r : ˆ D ( r ) = 1Ω r (cid:88) g ∈ Ω r A g ˆ σ g . (S3)16he symmetry of ˆ σ g is easy to establish by writing every contact force as the sum of a normalforce, which is along the contact vector r g,c , and a tangential force perpendicular to it. The normalpart leads to a symmetric contribution to ˆ σ g . Using the torque-balance equation, Eq. S1, thecontribution from the tangential forces sum up to zero. To establish the divergence free condition,we follow the approach outlined in Degiuli, E. and McElwaine, J. 2011 by first subdividing ˆ σ g into contributions from each contact. As seen from Fig. S1, we can associate a triangle of area a g,c with each contact, and A g = (cid:80) c ∈ g a g,c . Adopting a convention that we traverse around a grain ina counterclockwise direction, we associate with contact c , the triangle that is defined by c and thecontact c (cid:48) that follows it. We can then write: A g ˆ σ g = (cid:80) c ∈ g a g,c ˆ σ c , where ˆ σ c is yet to be defined.Comparing to Eq. S2, we see that a g,c ˆ σ c = r g,c ⊗ f g,c , therefore ˆ σ c = r g,c ⊗ f g,c /a g,c . The signedarea a g,c is given by a g,c = (1 / r g,c × ( r c (cid:48) − r c ) . The divergence theorem is: (cid:82) V ∂ i σ ij = (cid:82) ∂V n i σ ij ,where ˆ n is the unit normal to ∂V , which can be written as as (cid:82) V ∇ · ˆ σ = (cid:82) ∂V ( d r × ˆ σ ) j . We canapply the discrete version of this theorem to ˆ σ g to get: A g ( ∇ · ˆ σ ) g = (cid:88) c ∈ g ˆ σ c × ( r c (cid:48) − r c ) = (cid:88) c ∈ g f g,c = f ext . (S4)In the absence of external forces, ˆ σ g is divergence free. This grain-level condition leads to asimilar condition on ˆ D ( r ) : Ω r ∇ · ˆ D ( r ) = (cid:80) c ∈ ∂ Ω f c , where the sum is over the contact forces onthe boundary of Ω , which is still discrete.To map to the continuum theory, we posit that disorder averaging over all discrete networksthat occur under given external conditions leads to17 i ( ˆ D ( r )) ij = ( f ext ) j . We expect this mapping to be accurate if the coarse-graining volume Ω is much larger thana typical grain volume. The excellent correspondence between disorder-averaged ˆ D correlationsmeasured in granular packings and theoretical predictions, shown in the main text justifies theabove mapping. In Section 3 of this Supplementary Information, we present experimental mea-surements of ˆ D correlations in individual configurations to show that self-averaging is a very goodapproximation for internal stresses in granular media. In this section we present expressions for the correlations of the ˆ D tensor, analogous to the ex-pressions for the ˆ E correlations in vacuum (Eq. 7 in main text). The starting point is Eq. 5 in themain text: Gauss’s law and the magnetostatic condition for a polarizable medium characterized bythe rank-4 tensor, ˆΛ . In the vacuum theory , the strategy is to project out the divergence modefrom the completely isotropic rank-4 tensor, using the magnetostatic condition. This condition in q - space, for a polarizable medium is given by D ij ( q ) = ( ˆΛ − ˆ A ) ij ( q ); A ij ( q ) ≡ q ⊗ φ , (S5)where φ is the electrostatic gauge potential, as discussed in the main text. Since ˆΛ has to obeythe symmetry ij → ji , it is simpler to write the components of ˆ D as a vector of length in 2D: ( D xx , D yy , D xy ), and a vector of length in 3D. The rank-4 tensor can be then expressed as a18 × (2D) and a × (3D) matrix . Furthermore, if ˆΛ is a symmetric matrix in this representation,then the ˆ D − ˆ D correlations can be obtained from the ˆ E − ˆ E correlations by a transformation ofthe metric: q → ¯q ( ˆΛ) . Such a transformation is reminiscent of the emergence of birefringence inquantum spin ice in the presence of an applied electric field . For the more general situation thatcan occur in granular media the matrix is not symmetric, and a cleaner approach is to use the dualformalism in which the potential is obtained by solving Gauss’s law . In this dual formalism,potentials in 2D and 3D appear differently: a scalar in 2D and a second-rank symmetric tensor in3D. The expression for the correlations of the potentials can be worked out explicitly, and from thatthe ˆ D − ˆ D correlations can be obtained in a straightforward manner. In 2D, ∂ i D ij = 0 is solved byintroducing a potential
3, 26, 37, 38 , ψ : D ij = (cid:15) ia (cid:15) jb ∂ a ∂ b ψ . The potential in 3D is a symmetric tensor, ψ ij : D ij = (cid:15) iab (cid:15) jcd ∂ a ∂ c ψ bd Here, we present the explicit construction of the correlations of D ij in 2D
3, 26, 37 . The magne-tostatic condition implies that ˆΛ acts as a stiffness tensor in a Gaussian theory. Using the q -spacerepresentation: D ij ( q ) = (cid:15) ia (cid:15) jb q a q b ψ ( q ) , The correlations (cid:104) ψ ( q ) ψ ( − q ) (cid:105) can be computed, andgive: (cid:104) ψ ( q ) ψ ( − q ) (cid:105) = [ A ij ( q )Λ ijkl A kl ( − q )] − ,A ij = q δ ij − q i q j . (S6)The correlations of D ij then follow as: (cid:104) D ij ( q ) D kl ( − q ) (cid:105) = (cid:15) ia (cid:15) jb (cid:15) kc (cid:15) ld q a q b q c q d (cid:104) ψ ( q ) ψ ( − q ) (cid:105) . ˆΛ being a diagonal tensor with components λ i , i = xx , yy , xy , thecorrelations simplify to: C xxxx ( q ) = (cid:104) D xx ( q ) D xx ( − q ) (cid:105) = q y λ xx q y + λ yy q x + 2 λ xy q x q y ,C xyxy ( q ) = (cid:104) D xy ( q ) D xy ( − q ) (cid:105) = q x q y λ xx q y + λ yy q x + 2 λ xy q x q y ,C yyyy ( q ) = (cid:104) D yy ( q ) D yy ( − q ) (cid:105) = q x λ xx q y + λ yy q x + 2 λ xy q x q y , (S7) C xxxy ( q ) = (cid:104) D xx ( q ) D xy ( − q ) (cid:105) = − q x q y λ xx q y + λ yy q x + 2 λ xy q x q y ,C xxyy ( q ) = (cid:104) D xx ( q ) D yy ( − q ) (cid:105) = q x q y λ xx q y + λ yy q x + 2 λ xy q x q y ,C xyyy ( q ) = (cid:104) D xy ( q ) D yy ( − q ) (cid:105) = − q y q x λ xx q y + λ yy q x + 2 λ xy q x q y . The experimental and numerical measurements of correlations in 2D, shown in Fig. 1 ofthe main text and in Fig. S2 of the supplementary, have been fit to the above forms. To analyzethe correlations in isotropically compressed 3D packings, we assume that ˆΛ is the identity tensorand use Eq. 7 of the main text, which gives the correlations in vacuum with an overall stiffnessconstant, λ . A striking consequence of the anisotropic correlations in q -space is evident if we analyze thecorrelations of the normal stresses, D xx and D yy in real space. The Fourier Transform of C xxxx in20igure S2: Comparisons in Fourier space between the theoretical predictions (black line) with
Λ = , and the numericaland the experimental results (symbols) of the stress-stress correlations in 2D, isotropically jammed systems. a , Photo-elasticimages, in which each grain is shaded according to the magnitude of its normal stress, exhibit clear filamentary structures that arenormally referred to as force chains. b , Theoretical predictions of C xxxx ( q, θ ) and C xyxy ( q, θ ) , which are independent of q , andthe corresponding angular functions C xxxx ( θ ) and C xyxy ( θ ) . c , Numerical data of the frictionless jammed packings within therange of pressure P ∈ [0 . , . . The results of the five different system sizes N = 512 , , , , are shownin the angular plots. d , Experimental data from frictional packings within the range of pressure P ∈ [1 . × − , . × − ] .All correlation functions are normalized by their peak values of C xxxx ( θ ) . The units of q are π/L , where L is the system size: L ≈ d min in simulations, L = 40 d min in experiments. Here d min is the diameter of the small particle. Both the numericaland experimental data start to deviate from the theoretical predictions around q ≥ π/ d min , indicating the breakdown of thecontinuum limit. ˆΛ = λ illustrates the point: C xxxx ( r x , r y ) = 32 λr x for r x (cid:29) r y ,C xxxx ( r x , r y ) = − λr y for r y (cid:29) r x . (S8)The reverse is true for C yyyy . The consequence of this feature is that the transverse correlationsbecome negatively correlated. The photo-elastic images from 2D experiments, shown in the maintext and in Figs. S3 and S4, are a striking visual representation of this stark difference betweenlongitudinal and transverse correlations, which in turn is a manifestation of the conservation of“charge-angular-momentum”, and the resulting sub-dimensional propagation . The U (1) gaugetheory with vector charges, therefore, clarifies the meaning of force-chains within a continuum,disorder-averaged theory. Additional Analysis of Experiments
In this subsection, we present results of stress correlationsfrom individual configurations in the sheared experimental packings to illustrate how well self-averaging works in these jammed packings. We note that our systems are deep in the jammedregion: we do not address the possible breakdown of self-averaging close to the unjamming tran-sition. 22igure S3: Experimental measurements of correlations in Fourier space, for a single packing inthe ensemble of packings, used to generate the averaged correlations shown in the main text (Fig.3). The features observed in these averaged correlations, are seen to emerge in a single pack-ing, demonstrating the self-averaging property of the stress in these packings that are deep in thejammed regime. 23igure S4: Experimental measurements of correlations in Fourier space, for a second packingcreated under the same external conditions as in Fig. S324
Finite Temperature Results
Pinch point singularities are one of the salient features of the VCT correlation functions. Thesesingularities originate from the strict constraints of mechanical equilibrium imposed on athermalsystems. For a system at finite temperature however, these constraints can be violated and hencewe expect the pinch point singularities to disappear at finite temperatures. Thus, the presence ofa pinch point singularity is a hallmark of an athermal system. The numerically generated stresscorrelations from a 2D system at finite temperature is shown in Fig. S5 and it can be clearly seenthat the pinch point singularity has vanished at this temperature ( E thermal /E compression = 3 . .The numerical simulations were carried out in LAMMPS and the finite temperature wasimposed through a Nos´e -Hoover thermostat. The protocol is to start with a valid athermal T = 0 configuration, generated following the procedure described in the Numerical Methods Section ofthe main text and then perform finite temperature dynamics to compute the stress correlations at anon-zero temperature. This procedure is then repeated over multiple initial athermal configurationsand ensemble averaged to obtain the finite temperature stress correlations. The results displayedare obtained for packings of disks with an average pressure per grain P ∈ [0 . , . .The results shown have been averaged over 50 starting athermal configurations in 2D with 50finite temperature configurations sampled during the finite temperature molecular dynamics run,for each of the 50 starting configurations. 25 .10.20.30.40.50.60.70.80.91 -0.3-0.2-0.100.10.20.3 -0.3-0.2-0.100.10.20.3 -0.06-0.04-0.0200.020.040.060.080.1 a b c d Figure S5: Comparisons in Fourier space between stress correlations at zero (
Top ) and finite (
Bot-tom ) temperatures . The columns a , b , c , and d show the results for correlation functions C xxxx , C xyxy , C xxxy and C xxyy respectively. The packings used have an average compression energy pergrain E compression ≈ − and the finite temperature results have an average thermal energy pergrain E thermal ≈ . × − . References
1. Landau, L. D., Pitaevskii, L. P., Kosevich, A. M. & Lifshitz, E.
Theory of Elasticity (Elsevier,2012).2. Ashcroft, N. W. & Mermin, N. D.
Solid State Physics (Saunders College Publishing, 1976).3. Henkes, S. & Chakraborty, B. Statistical mechanics framework for static granular matter.
Phys. Rev. E , 061301 (2009). 26. DeGiuli, E. Field Theory for Amorphous Solids. Phys. Rev. Lett. , 118001 (2018).5. Lemaˆıtre, A. Stress correlations in glasses.
Journal of Chemical Physics , 104107 (2018).6. Geng, J. et al.
Footprints in sand: The response of a granular material to local perturbations.
Phys. Rev. Lett. , 035506 (2001).7. Bouchaud, J.-P. Granular media: some ideas from statistical physics. In Bouchaud, J., Barrat,J. L., Feigelman, M., Kurchan, J. & Dalibard, J. (eds.) Slow Relaxations and NonequilibriumDynamics in Condensed Matter , vol. 77, 185–202 (Les Ulis: EDP Sciences, 2003).8. Otto, M., Bouchaud, J. P., Claudin, P. & Socolar, J. E. Anisotropy in granular media: Classicalelasticity and directed-force chain network.
Phys. Rev. E , 24 (2003).9. Diep, H. et al. Frustrated spin systems (World Scientific, 2013).10. Huse, D. A., Krauth, W., Moessner, R. & Sondhi, S. L. Coulomb and liquid dimer models inthree dimensions. Phys. Rev. Lett. , 167004 (2003).11. Domb, C., Green, M.S. & Lebowitz, J.L. Phase Transitions and Critical Phenomena Vol. 13edited by by C. Domb and J. Lebowitz (Academic Press, London, 1989).12. Castelnovo, C., Moessner, R. & Sondhi, S. L. Magnetic monopoles in spin ice.
Nature ,42 (2008).13. O’Hern, C. S., Silbert, L. E., Liu, A. J. & Nagel, S. R. Jamming at zero temperature and zeroapplied stress: The epitome of disorder.
Phys. Rev. E , 011306 (2003).274. Bi, D., Zhang, J., Chakraborty, B. & Behringer, R. P. Jamming by shear. Nature , 355(2011).15. Cates, M. E., Wittmer, J. P., Bouchaud, J.-P. & Claudin, P. Jamming, force chains, and fragilematter.
Phys. Rev. Lett. , 1841–1844 (1998).16. Peters, I. R., Majumdar, S. & Jaeger, H. M. Direct observation of dynamic shear jamming indense suspensions. Nature , 214 (2016).17. Cates, M. E., Fuchs, M., Kroy, K., Poon, W. C. K. & Puertas, A. M. Theory and simulation ofgelation, arrest and yielding in attracting colloids.
Journal of Physics: Condensed Matter ,S4861 (2004).18. Zhang, S. et al. Correlated Rigidity Percolation and Colloidal Gels.
Physical Review Letters , 33–35 (2019).19. Bi, D., Henkes, S., Daniels, K. E. & Chakraborty, B. The Statistical Physics of AthermalMaterials.
Annual Review of Condensed Matter Physics , 63–83 (2015).20. DeGiuli, E. Continuum limits of granular systems (Thesis, U. British Columbia, (2012).21. Pretko, M. Generalized electromagnetism of subdimensional particles: A spin liquid story. Phys. Rev. B , 1–26 (2017).22. Pretko, M. The fracton gauge principle. Phys. Rev. B , 115134 (2018).23. Pretko, M. & Radzihovsky, L. Fracton-Elasticity Duality. Phys. Rev. Lett. , 195301(2018). 284. Majmudar, T. S. & Behringer, R. P. Contact force measurements and stress-induced anisotropyin granular materials.
Nature , 1079–1082 (2005).25. Prem, A., Vijay, S., Chou, Y. Z., Pretko, M. & Nandkishore, R. M. Pinch point singularitiesof tensor spin liquids.
Phys.Rev. B , 1 (2018).26. Lois, G. et al. Stress correlations in granular materials: An entropic formulation.
Phys. Rev. E , 060303(R) (2009).27. Xu, C. Novel Algebraic Boson Liquid phase with soft Graviton excitations arXiv: 0602443v4 (2006).28. Gelin, S., Tanaka, H. & Lemaˆıtre, A. Anomalous phonon scattering and elastic correlations inamorphous solids. Nature Materials , 1177 (2016).29. Maier, M., Zippelius, A. & Fuchs, M. Emergence of Long-Ranged Stress Correlations at theLiquid to Glass Transition Phys. Rev. Lett. , 265701 (2017).30. Wu, B., Iwashita, T. & Egami, T. Anisotropic stress correlations in two-dimensional liquids.
Phys. Rev. E , 032301 (2015).31. O’Hern, C. S., Langer, S. A., Liu, A. J. & Nagel, S. R. Random packings of frictionlessparticles. Phys. Rev. Lett. , 075507 (2002).32. Ramola, K. & Chakraborty, B. Stress Response of Granular Systems. Journal of StatisticalPhysics (2017). 293. Gu´enol´e, J. et al.
Assessment and optimization of the fast inertial relaxation engine (fire)for energy minimization in atomistic simulations and its implementation in lammps arXiv:1908.02038 (2019).34. Bitzek, E., Koskinen, P., G¨ahler, F., Moseler, M. & Gumbsch, P. Structural relaxation madesimple.
Phys. Rev. Lett. , 170201 (2006).35. Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. Journal of Compu-tational Physics , 1 (1995).36. Wang, Y., Hong, L., Wang, Y., Schirmacher, W. & Zhang, J. Disentangling boson peaks andvan hove singularities in a model glass.
Phys. Rev. B , 174207 (2018).37. DeGiuli, E. Edwards field theory for glasses and granular matter. Phys. Rev. E , 33001(2018).38. Xu, C. Gapless bosonic excitation without symmetry breaking: An algebraic spin liquid withsoft gravitons. Physical Review B - Condensed Matter and Materials Physics , 1–11 (2006).39. DeGiuli, E. & McElwaine, J. Laws of granular solids: Geometry and topology. PhysicalReview E - Statistical, Nonlinear, and Soft Matter Physics (2011).40. Lantagne-Hurtubise, ´E., Bhattacharjee, S. & Moessner, R. Electric field control of emergentelectrodynamics in quantum spin ice. Physical Review B , 1–20 (2017). Acknowledgements
The authors acknowledge Michael D’Eon, Albion Lawrence, Debarghya Banerjee,Chandan Dasgupta, Kedar Damle, Roderich Moessner, Nandagopal M, Prasad Perlekar, Samriddhi Sankar ay and Vijay Shenoy for fruitful discussions. JN and BC acknowledge support from NSF-CBET-1916877and BSF-2016188. BC was supported by a Simons Fellowship in Theoretical Physics. SB acknowledgesfinancial support through Max Planck partner group on strongly correlated systems at ICTS; SERB-DST(Govt. of India) Early Career Research grant (No. ECR/2017/000504) and the Department of AtomicEnergy, Government of India, under project no.12-R&D- TFR-5.10-1100. YW and JZ acknowledge thesupport from NSFC under (No. 11774221 and No. 11974238). BC, JNN, KR and SB would like to thankthe ICTS program: Entropy, Information and Order in Soft Matter (August-November 2018) during whicha part of the work was done. Competing Interests
The authors declare that they have no competing financial interests.
Author Contributions
SB and BC conceived and constructed the theory. JN and KR designed and per-formed the numerical simulations and analyzed the numerical data. YW and JZ conceived and executed theexperimental project. All authors contributed equally to the manuscript.
Correspondence