Unjamming dynamics: the micromechanics of a seismic fault model
aa r X i v : . [ c ond - m a t . s o f t ] J un Unjamming dynamics: the micromechanics of a seismic fault model
Massimo Pica Ciamarra, ∗ Eugenio Lippiello, Cataldo Godano, and Lucilla de Arcangelis
3, 4 CNR–SPIN, Department of Physical Sciences, University of Naples Federico II, 80126 Napoli, Italy. † Dep. of Environmental Sciences and CNISM, Second University of Naples, 81100 Caserta, Italy Dep. of Information Engineering and CNISM, Second University of Naples, 81031 Aversa (CE), Italy Institute for Building Materials, Schafmattstr. 6, ETH, 8093 Z¨urich, CH
The unjamming transition of granular systems is investigated in a seismic fault model via threedimensional Molecular Dynamics simulations. A two–time force–force correlation function, and asusceptibility related to the system response to pressure changes, allow to characterize the stick–slipdynamics, consisting in large slips and microslips leading to creep motion. The correlation functionunveils the micromechanical changes occurring both during microslips and slips. The susceptibilityencodes the magnitude of the incoming microslip.
PACS numbers: 45.70.-n; 46.55.+d; 45.70.Ht; 91.30.Px
In a number of industrial processes and natural phe-nomena, such as earthquakes or landslides, disorderedsolid granular systems start to flow. This solid-to-liquidtransition, known as unjamming, occurs either on de-creasing the confining pressure P , or increasing the ap-plied shear stress σ . Understanding the properties ofthis transition is a big challenge due to the absence ofan established theoretical framework for granular mate-rials. A proposed analogy with the glass transition [1] ofthermal systems has recently triggered the study of thejamming transition via numerical investigations of sys-tems of soft frictionless particles at zero applied shearstress [2], where the only control parameter is the pres-sure (or the density). As the unjamming transition isapproached by decreasing the confining pressure, the vi-brational spectrum develops an excess of low frequencymodes, known as soft–modes, leading to the identifica-tion of a length scale which diverges on unjamming [3].This length scale is related to the emergence of an in-creasingly heterogeneous response as the system movestowards the transition [3]. A different approach to thestudy of the unjamming transition has been followed ina two dimensional numerical study [4] and in a numberof experiments [5–8], where the applied shear stress iscontrolled via a spring mechanism, as the one in Fig. 1a.A stick–slip motion characterized by a complex slip sizestatistics [7] is recovered at high confining pressures P and small driving velocities V d . This stick-slip dynamicsis altered by the presence of noise [9]. Analogous resultshave been found at fixed strain rate [10, 11].In this Letter we tackle this problem via three di-mensional Molecular Dynamics simulations of a modelof a seismic fault (Fig. 1a), where grains play the roleof the gouge [5–8, 10, 11]. Numerical details are givenin [13, 14]. The micromecanical mechanisms leading tothe transition are analyzed at a level of spatial and tem-poral resolution not considered before. Stick–slip dynamics –
For the investigated valuesof the parameters [13], the system is characterized bystick–slip dynamics, which we analyse considering that a t X(t) -5 -4 -3 -2 -1 ∆ X / L x -10 -8 -6 -4 -2 n( ∆ X) a bc d → V d ↓ P FIG. 1: (Color online) (a) The system consists of granular par-ticles (light grains) confined between two rigid plates (darkgrains) at constant pressure. The top plate is driven via aspring mechanism, where an extremum of the spring is at-tached to the plate and the other is pulled at constant veloc-ity. See the supplementary materials available online for ananimation of a slip event in the real and in the force space [12].(b) Time evolution of the top plate position X ( t ). The insetindicates that between two large slips there are many mi-croslips. (c) Configuration of normal force network. Strongerforces correspond to thicker and darker lines. (d) Slip sizedistribution for slips and microslips. slip begins and ends when the velocity of the top platebecomes, respectively, larger and smaller than a smallenough threshold. We measure the displacements ∆ X ofthe top plate due to slips, and compute their distribution n (∆ X ) (Fig. 1d). For slips smaller than ≃ . L x , where L x is the system length, the distribution follows a powerlaw, n (∆ X ) ∝ ∆ X − β with β ≃ .
85, in agreement withexperimental values for earthquakes [16]. Larger slips arealmost periodic in time and roughly follow a lognormaldistribution with a characteristic size ∆ X ≃ . L x . Sum-marizing, the dynamics consists in the occurrence of al-most periodic large events, here called slips , and of creep time X ( t ) ∆ t X ( t ) t u t u + δ t t u FIG. 2: (Color online) The top plate position in a slip event X ( t ). The vertical dashed line indicates the unjamming time t u . Inset: Top plate position for two system replicas, whichevolve with V d = 0. The replica made at time t u remains inthe jammed configuration (continuous line), whereas the onemade at time t u + δt slips (dashed line). motion characterized by smaller events, here called mi-croslips , in agreement with previous experiments [7, 11]. Onset of a slip –
To understand the mechanisms actingat the onset of a slip, we need to the identify its precisestarting time t u . Here we describe the analysis performedon the particular slip occurring at time t ≃ t , and follow itstime evolution at zero driving velocity V d = 0. If thereplica made at time t resists to the applied stress, then t ≤ t u . Conversely, if a slip is observed, t > t u . Wedefine t u as the largest time where no slip occurs, andwe identify it (one for each slip) with an accuracy equalto the time step of integration of the equation of motion δt [13] (Fig. 2). This procedure is equivalent to a quasi-static simulation [17] around the un-jamming time t u andgives the value of the shear stress above which the systemstarts to flow.We have performed a number of checks which suggestthat no structural changes occur at t u . For instance, thecomparison of the state of the system at time t u withthe one at shortly earlier and later times, shows that nocontact breaks at t u . We have also considered the distri-bution of the parameter λ = | f t | /µ | f n | where f t and f n are the tangential and the normal forces. When λ > P ( λ ) graduallymoves toward 1 as t u is approached, indicating the weak-ening of the solid [4]. However, neither at t u the numberof contacts with λ ≃ t u supports a scenario in which thesystem is located in an energy minimum which slowlybecomes an inflection point at time t u , and therefore the smallest eigenvalue of the dynamical matrix continuouslydecreases to zero [18, 19]. Evolution in the force space –
When the system sticksand the shear stress increases, no macroscopic motion isobserved. However, the system microscopically changessince it sustains an increasing shear stress. The time evo-lution of the top plate position X ( t ) (Fig. 3a), consists inan elastic deformation, where X ( t ) increases very slowlyin time, interrupted by sudden microslips. To charac-terize the evolution of the system in the force space, weintroduce a two time force-force correlation function forthe normal forces, defined as C n ( t , t ) = P ij | f nij ( t ) || f nij ( t ) | P ij | f nij ( t ) | , (1)where the sum running over all couples of particles ( i, j )corresponds to a spatial average. An equivalent defini-tion holds for the correlation of tangential forces, C tg .Being interested in the unjamming transition of the slipevent shown in Fig. 2, here we fix t = t u , and con-sider the evolution of the correlation function C n ( t u , t )for earlier times, t < t u (Fig. 3b). The force correla-tion function C n (and C tg , not shown) exhibits smalljumps in correspondence of microslips (Fig. 3a), reveal-ing the unusual occurrence of bursts in the reorganiza-tion of the force network. During these bursts, the en-ergy due to the tangential interaction decreases, whereasthe one due to the normal interaction increases. A pos-sible interpretation is in terms of a two force networkscenario, in which the applied stress σ is supported bya stress σ n due to the normal force network, and by astress σ t due to the tangential forces, σ = σ n + σ t . Ina burst, few contacts break, leading to a decrease of σ t , σ t → σ t − δσ . A microscopic slips is observed since thenormal forces quickly adapt and succeed in sustaining theapplied stress, σ n → σ n + δσ . This scenario is supportedby the inset of Fig. 3b, which shows both C n and C tg across a microslip, with t its starting time. C n slightlyincreases and overcomes 1, while C tg exhibits a sharpdrop due to the breaking of several contacts.The force correlation function (Eq. 1) gives also in-sights into the system evolution during a slip. In Fig. 4 weplot C n,tg ( t u , t ) for t > t u , and for comparison the scaledtop plate position ( X ( t ) − X ( t u )) / ( X ( t ∞ ) − X ( t u )), where t ∞ is a time following the slip event, whose precise valuedoes not influence our results. We first notice that theforces evolve on a timescale much shorter than the platemotion. For instance the correlation functions reach thevalue 0 .
1, denoting an almost complete relaxation, whenthe top plate moved only by 10% of its total displace-ment. The presence of different time scales in the relax-ation process is evidenced by the self-scattering correla-tion function F ( q, t ) = N (cid:12)(cid:12)(cid:12)P j exp[ i q · ( r j ( t ) − r j ( t u ))] (cid:12)(cid:12)(cid:12) ,where r j represents the position of the j th particle. Sincethe system is sheared along x and confined along z , we time n (t,t u ) t u ab ∆ X C n C tg FIG. 3: (Color online) Time evolution of the top plate posi-tion (a) and of the correlation function of the normal forces C n ( t, t u ) (b). Dashed lines identify the unjamming time t u .The inset show C n and C tg across a microslip at t = 1579 . have considered wave vectors along y , q = (0 , q, q , F ( q, t ) probes small scale relaxation, and coin-cides with C n ( t u , t ), as shown in Fig. 4. Conversely, atsmall q , F ( q, t ) relaxes on a time scale comparable tothat of the upper plate motion. The relaxation time τ q , F ( q, τ q ) = 1 /e , indeed increases as q decreases. More-over, tangential forces decorrelate before normal ones.This can be explained considering the unjamming tran-sition as a buckling-like instability of the chains of largenormal forces, which are sustained by weaker tangentialcontacts. When the weaker sustaining contacts break,either the normal forces adapt to sustain the extra load,leading to a microslips, or a buckling-like instability oc-curs, giving rise to a slip. The same quantities canbe used to investigate the subsequent jamming transi-tion. To this end, the force network at time t is com-pared with the force network after the slip event study-ing C n,tg ( t, t ∞ ). Tangential forces correlate after normalones, and makes the force network stable. Response to perturbations: slips versus microslips –
The correlation length ξ in equilibrium systems is mea-sured from the response to an external perturbation. Thesusceptibility, for instance, scales like ξ − η near criticalpoints, where η is the correlation function critical expo-nent. Here we measure the response of the system to apressure change at the onset of both slips and microslips.More precisely, at each time t we stop the external drive,setting V d = 0, and introduce a perturbation in the exter-nal pressure P , fixed to P ′ = P (1 − α ) for a time interval δt perb = 0 .
1. The response χ α at time t is defined as χ α ( t ) = 1 αP lim τ →∞ " N X i ( r αi ( t + τ ) − r i ( t + τ )) / (2)where { r αi } and { r i } are the asymptotic states of the per- C ( t , t ) C n C tg time C ( t , t ∞ ) q τ q τ q s ca l e d p l a t e po s iti on t u FIG. 4: (Color online) Correlation functions at the unjam-ming (upper panels) and at the subsequent jamming transi-tion (lower panels). The left panels show the normal and tan-gential force correlation functions C n,tg ( t , t ) (symbols) dur-ing a slip event. The dashed line shows the time evolutionof the plate position scaled between 0 and 1. The plain linesshow the self-scattering function F ( q, t ) for q = 3 , , , τ q are shown in the right panels. t is fixedequal to t u for the unjamming transition, and to t = 1614for the subsequent jamming. turbed of the unperturbed systems. In the unjammedphase, the susceptibility is divergent. In the jammedphase, it measures the size of the region of correlatedparticles that respond to the external perturbation, pro-viding an estimate of the correlation length. It is a staticquantity since the time t only indicates the instant atwhich the perturbation is applied. χ α can be also de-fined without setting V d = 0, provided that the charac-teristic response time of the system is much smaller thanthe timescale over which the applied stress varies. For awide range of α , the response of the system at times farfrom t u is linear in the perturbation, since χ α does not de-pend on α (Fig. 5). In particular, χ α gradually increasesin time and drops in correspondence to microslips. Theinset shows that the microslip size ∆ X depends on χ α evaluated just before the slip, as ∆ X ≃ χ bα . This indi-cates that the size of a microslip is already encoded inthe system state. In fact, considering that χ α increaseson approaching a microslip, the measured χ α provides alower bound for the magnitude of the incoming event.As the unjamming time is approached, the responseis no longer linear. χ α remains roughly constant untilit abruptly increases at a time which depends on α , thesooner the greater α . This increase is consistent with apower-law divergence. Accordingly at each time, namelyat each value of the applied shear stress, there is a min-imum value of the perturbation intensity for slip trig-gering. This behavior is in line with the existence ofa minimum threshold amplitude in the deformation as-sociated with seismic waves for earthquake remote trig- time χ α -3 -2 χ α -3 -2 -1 ∆ X t u FIG. 5: (Color online) Time evolution of the susceptibility χ α for different α . (Inset) Size of a microslip ∆ X versus thecorresponding value χ α . The straight line is χ bα , with b ≃ . gering [20]. A difference between slips and microslips isin how different particles contribute to χ α in the sumof Eq. 2. Indeed, these contributions are very similarfor microslips, and highly heterogeneous for slips. Thepresence of an heterogeneous response is consistent withprevious numerical results found at σ = 0 [3].In conclusion, the absence of precise structural changesat the unjamming time, and the bursts observed inthe prior dynamics, suggest that the increasing exter-nal stress progressively modifies the underlying energylandscape.At each time the system is in an equilibrium positionwhich can be seen as local energy minimum of an effec-tive energy landscape which depends on the applied shearstress. Microslips occur when the local energy minimumflattens down as the applied shear stress increases, lettingthe system fall in a neighbor minimum. If there are noclose minima, a slip occurs, and the system jumps to afar away configuration. The deforming energy landscapepicture also suggests that soft-modes could be found notonly when unjamming is approached at σ = 0 by de-creasing the volume fraction φ [3], but also [21] whenunjamming is approached increasing σ .We thank A. Coniglio for helpful discussions and theUniversity of Naples Scope grid project, CINECA andCASPUR for computer resources. ∗ [email protected] † URL: http://smcs.na.infn.it [1] A. J. Liu, S. R. Nagel, Nature , 21 (1998).[2] C.S. O’Hern et al. , Phys. Rev. Lett , 075707 (2002);C.S. O’Hern et al. , Phys. Rev. E , 011306 (2003).[3] L.E. Silbert, A.J. Liu, and S.R. Nagel, Phys. Rev. Lett. , 098301 (2005); E. Somfai, M. van Hecke, W.G. Ellen- broek, K. Shundyak, and W. van Saarloos Phys. Rev. E , 020301 (2007); W.G. Ellenbroek, E. Somfai, M. vanHecke, and W. van Saarloos, Phys. Rev. Lett. ,B09306 (2004).[5] S. Nasuno, A. Kudrolli, and J.P. Gollub, Phys. Rev. Lett. , 949 (1997).[6] J.P. Gollub, G. Voth, and J.C. Tsai, Phys. Rev. Lett. ,064301 (2003).[7] F. Dalton and D. Corcoran, Phys. Rev. E , 061312(2001); Phys. Rev. E , 031310 (2002).[8] M. Bretz, R. Zaretzki, S.B. Field, N. Mitarai, F. Nori,Europhys. Lett. , 1116 (2006). K. Daniels and N.W.Hayman, J. Geophys. Res., , B11411 (2008).[9] R. Capozza, A. Vanossi, A. Vezzani and S. Zapperi, Phys.Rev. Lett. , 085502 (2009); M.F. Melhus, I.S. Aran-son, D. Volfson and L.S. Tsimring, Phys. Rev. E ,041305 (2009).[10] A.A. Pe˜na, S. McNamara, P.G. Lind and H.J. Herrmann,Granular Matter , 243 (2009).[11] P.A. Johnson and X. Jia, Nature , 871 (2005); P.A.Johnson, H. Savage, M. Knuth, J. Gomberg and C.Marone, Nature , 57 (2008).[12] See EPAPS Document No. E-PRLTAO-XXXXXX for ananimation of a slip in the real and in the force space.[13] We consider N = 10 particles of mass m and diameter D confined in the z -direction by two rigid rough plateswith size L x = 20 D , L y = 5 D . Periodic boundary condi-tions are imposed along x and y . Particles interact in thenormal direction via the standard linear–spring dashpotmodel [14], characterized by an elastic constant k n and bya restitution coefficient e = 0 .
88. In the tangential direc-tion, we use the spring model developed in Ref. [15] withelastic constant k t = 2 / k n and coefficient of friction µ = 0 .
1. Lengths, masses, times and stress are expressedin units of D , m , p k n /m and k n /D , respectively. Theconfining pressure is P = 100 and k d V d = 10, where k d isthe elastic constant of the driving spring. The time stepof integration of the equation of motion is δt = 2 10 − .The slipping velocity threshold is v t = 10 − .[14] M. Pica Ciamarra, L. De Arcangelis, E. Lippiello, C. Go-dano, Int. J. Mod. Phys. B , 5374 (2009).[15] P.A. Cundall and O.D.L. Strack, Geotechnique , 1073 (1975); Y.Y. Kagan, Geophys. J. Int. , 123(1991); C. Godano and F. Pingue, Geophys. J. Int. ,193 (2000).[17] C. Heussinger and J.-L. Barrat, Phys. Rev. Lett. ,218303 (2009).[18] P.R. Welker and S.C. Mcnamara, Phys. Rev. E ,061305 (2009).[19] C.E. Maloney and A. Lemaitre, Phys. Rev. E , 016118(2006).[20] J. Gomberg and P.l Johnson, Nature , 830 (2005).[21] M. Pica Ciamarra and A. Coniglio, Phys. Rev. Lett.23