Gustavo
Düring
,
Edan
Lerner
and
Matthieu
Wyart
*
New York University, Center of Soft Matter Research, 4 Washington Place, New York, NY 10003, USA. E-mail: mw135@nyu.edu
First published on 15th October 2012
Gels of semi-flexible polymers, network glasses made of low valence elements, softly compressed ellipsoid particles and dense suspensions under flow are examples of floppy materials. These systems present collective motions with almost no restoring force. To study theoretically and numerically the frequency-dependence of the response of these materials and the length scales that characterize their elasticity, we use a model of isotropic floppy elastic networks. We show that such networks present a phonon gap for frequencies smaller than a frequency ω* governed by coordination, and that the elastic response is characterized, and in some cases localized, on a length scale that diverges as the phonon gap vanishes (with a logarithmic correction in the two dimensional case). lc also characterizes velocity correlations under shear, whereas another length scale l* ∼ 1/ω* controls the effect of pinning boundaries on elasticity. We discuss the implications of our findings for suspension flows, and the correspondence between floppy materials and amorphous solids near unjamming, where lc and l* have also been identified but where their roles are not fully understood.
In simplified numerical models of ellipsoid particles,7,8 covalent networks,18 gels of semi-flexible polymers19 and suspension flows,20 it has been observed that the density of vibrational modes D(ω) displays a gap at low frequency (in addition to the floppy or nearly floppy modes present at zero-frequency). In suspension flows the amplitude of the gap was shown to affect the divergence of the viscosity near jamming.20 An early work by Garboczi and Thorpe21 supported the idea that a gap of modes exists in floppy materials. However, the dependence of the gap on the microscopic structure and its consequences on the material properties and the different length scales characterizing the elastic response are not understood.
In this manuscript, we study the elastic properties of isotropic harmonic spring networks that are strongly disordered, but where spatial fluctuations of coordination are weak. The model networks we use are presented and analyzed numerically in Section 2. In Section 3, we show, using both mean-field methods21–25 and numerical simulations, that strictly floppy elastic networks present a vibrational gap between floppy modes (zero frequency modes) and a frequency ω* ∼ zc − z ≡ δz. The absence of low frequency phonons suppresses the elastic propagation of the response at frequencies smaller than ω*. As a consequence, we show that the mean response to a local perturbation (local strain) is localized, displaying an exponential decay with a characteristic length . We also predict a logarithmic correction to the scaling for two dimensional systems. In Section 4 we go beyond the mean field and estimate the fluctuations around the mean response. We show that fluctuations dominate the amplitude over the mean value, however, the fluctuation around the mean response also decays with the same length scale ∼lc. Even though lc can be considered to be the localization length of floppy modes, we show that these modes cannot exist in a region of typical radius smaller than l* ∼ 1/δz. Surprisingly, floppy modes can be strictly localized (have a compact support) only on the much larger scale l*. In Section 5 we study the floppy networks stabilized by weak interactions, which confer a finite but small restoring force to the material. We model the weak interaction by adding springs whose stiffnesses are very small. We find that localization is lost even when weak interactions have a vanishingly small amplitude, showing that the response of strictly floppy systems is a singular limit. However, lc still characterizes the response near the applied strain. In Section 6 we compare our results with recent observations in the affine solvent model, a simplified model for suspension flows of hard spheres. For this strictly floppy system, the spectrum associated with the evolution operator shares both commonalities and striking differences with the density of states of the isotropic floppy networks considered here, leading to a prediction on some dynamical length scale in flow. Finally, our results raise questions associated with the respective role of lc and l* both in floppy materials and in jammed packings, which are discussed in the last section.
![]() | (1) |
The force field ∣F〉 generated by displacements can be obtained by taking the derivative of eqn (1). One obtains ∣F〉 = ∣δR〉, where:
![]() | (2) |
As a model system we consider isotropic disordered elastic networks, generated by following the method of ref. 26: we prepare amorphous packings of compressed soft elastic particles with a coordination value significantly larger than zc. The centers of the particles are the nodes of our networks, and the contacts between particles are replaced by un-stretched springs of stiffness k. Springs are then removed until the desired coordination z is reached. Removal takes place randomly from the set of most connected pairs of nodes, leading to isotropic networks with low heterogeneity in density and coordination. Such networks are appropriate models of amorphous solids for which large spatial heterogeneities are not energetically favorable. In rigidity percolation,21,23 spring removal is completely random, resulting in networks with large fluctuations that affect the elastic properties.
For these networks, we diagonalize numerically to compute the density of states D(ω) for various z < zc. As shown in the inset of Fig. 1, we find that a gap in the vibrational spectrum appears below some frequency ω* that decreases as z → zc. If spatial fluctuations of coordination were large, one would expect this gap to be filled-up in the thermodynamic limit due to Griffiths-like singularities, but we do not observe this effect in the network considered here. The role of this gap in elasticity can be analyzed by considering the response to a local strain. We change the rest length of one spring and let the system relax to zero energy under over-damped dynamics. Such a response can also be obtained in the absence of damping, by imposing a local oscillatory strain at a vanishing small frequency. As shown in Fig. 2, the elastic information does not propagate: the response is localized on some length scale that appears to diverge as the gap vanishes, i.e. as z → zc. The response appears to be very heterogeneous. In what follows we will investigate its mean and its fluctuations.
![]() | ||
Fig. 1 Rescaled density of states vs. rescaled frequency for z ∈ [3.2, 3.95] using N = 10![]() |
![]() | ||
Fig. 2 Displacement field caused by a local strain in a floppy network. An over-damped relaxation was used following the elongation of one spring for (a) δz = 0.2 and (b) δz = 0.05. |
The transfer matrix is found to be:
![]() | (3) |
We seek an effective stiffness eff that captures the average behavior of the system, i.e. 〈G〉 = G0, where the average is taken over the disorder on the stiffness coefficients kij. This condition leads to 〈
〉 = 0. In the EMT this constraint is approximated by 〈T〈ij〉〉 = 0. Using standard identities for the Green's function on isotropic lattices,27 one can express this condition as:23
![]() | (4) |
Since we are interested in low frequencies , we approximate G0 by its continuum limit and use a Debye cut-off qD. We introduce the bulk modulus K and the shear modulus μ of the ordered lattice without an effective stiffness, and use the approximation:
![]() | (5) |
The sum is taken over all polarization vectors qα, where for pressure waves and
for shear waves.
εf(ε) = δz/2. | (6) |
For spatial dimensions d ≥ 3, eqn (5) implies that f(0) is a positive constant, therefore in the regime δz ≪ 1 one finds
ε ≈ δz/(2f(0)). |
Using eqn (4) together with the assumption that ∣mω2/keff∣ ≪ 1 and ∣keff∣ ≪ 1 (which can be shown to be true a posteriori in the limit of δz ≪ 1 and ) one finds:
![]() | (7) |
For d = 2, eqn (5) leads to a logarithmic divergence in the small ε limit:
![]() | (8) |
The asymptotic value for ε in the limit δz → 0 can be obtained from eqn (6), which leads to and therefore f(ε) ∼ −log(δz). This behavior is only valid for values of |log(δz)| ≫ 1, a difficult limit to observe empirically. Therefore, to compare our theoretical predictions with two dimensional numerical observations we compute ε(δz) by solving eqn (6) numerically using eqn (8). We find that the relationship ε(δz) depends on the initial network via only one parameter Λ(cp,cs,Γ).† To compare the theory and numerical simulations, we use Λ(cp,cs,Γ) as a fitting parameter.
The mechanical response of floppy systems is therefore fully determined by the Green's function in eqn (5) and the effective stiffness given by eqn (7). We start by computing the density of states D(ω) using the relationship . For 0 < ω < ω*, where
, we find that D(ω) = 0. For ω > ω* we obtain:
![]() | (9) |
Thus the size of the gap scales linearly with δz for d ≥ 3, where f(ε) ≈ f(0), whereas for d = 2, a logarithmic correction exists. According to eqn (9), rescaling the frequency by ω* and D(ω) by f(ε)1/2 should collapse the low-frequency part of the spectrum. Fig. 1 shows that the quality of the collapse is very good. For all considered coordinations, we used the same fitting parameter Λ(cp,cs,Γ) = 1.3, which is fixed by this measurement.
In the small frequency regime, the last expression is found to be |F〉 ∼ keffeiω0tn(|i〉 − |k〉). The magnitude of the force vanishes as ω0 → 0, which is consistent with the existence of floppy modes.
The response at the zero frequency limit in two dimensions, at distances r = |r| larger than the typical springs length (i.e. rqD ≫ 1), can be calculated using the continuum limit
〈δR(r)〉 = g(r − ri) − g(r − rk), | (10) |
![]() | ||
Fig. 3 Mean response following the elongation of a spring placed at the center, along the vertical axis. (a) Theoretical prediction given by eqn (10) using cs and cp of a hexagonal lattice. Arrows show the displacement field and continuous lines are stream lines of the vector field. (b) Average response over 6000 independent realizations of networks with N = 40![]() |
The asymptotic solution in any dimension has an exponential decay. Indeed, taking the angular average one gets
![]() | (11) |
![]() | ||
Fig. 4 (a) Mean response δRmvs. the distance r from the imposed strain with N = 40![]() |
The zero frequency limit can be extended to finite frequencies by replacing ε → −mω2/keff in eqn (10). The asymptotic behavior is then given by log(〈δR(r)〉) ∼ −r/lc(ω) + iωr/v(ω), where
From eqn (7) one can prove the existence of two regimes. (i) For ω < ω*, the imaginary part 1/v(ω) = 0: one finds a pure exponential decay, with a characteristic length of order lc. (ii) For ω ≫ ω*, the decay length decreases as lc(ω) ∼ f(ε)1/4ω−1/2, whereas the velocity of the corresponding vibrations grows as v(ω) ∼ f(ε)1/4ω1/2. A similar behavior above ω* has been predicted by one of us for z ≥ zc.24
We can also consider the case of our networks immersed in a viscous fluid with viscosity η0 and without hydrodynamic interactions. The response in the over-damped regime is calculated by replacing the inertial term mω2 by the viscous force −iη0ω in eqn (5) and then in the effective stiffness in eqn (7). In this limit the elastic modulus is proportional to keff, whose real and imaginary parts give the storage and loss moduli G′ and G′′ respectively. Two different regimes are observed. (i) For ω ≪ ω*2m/η0, the storage moduli follow G′ ∼ lc6ω2/f(ε) and the loss moduli G′′ ∼ lc2ω. The response is given by lc(ω) ∼ lc and v(ω) ∼ f(ε)lc−3, displaying a similar exponential decay to that observed in an undamped system, but with propagating waves. (ii) For ω ≫ ω*2m/η0, both the storage and loss moduli are proportional to , as also observed in ref. 29. The decaying length lc(ω) ∼ f(ε)1/4ω1/4, while the velocity v(ω) ∼ f(ε)1/4ω3/4. Our scaling for the characteristic length scale lc(ω) differs from the recent result of Tighe.29 In our opinion, the difference stems from an incorrect definition of the length scale characterizing the elastic response in floppy materials (see footnote‡).
We denote the dipole of forces ∣Fij〉 generated by changing the rest length of the spring ij by a distance one:
|Fij〉 = nij(|i〉 − |j〉), | (12) |
![]() | (13) |
Thus, the response |δRij〉 to the elongation of a spring ij, which is equivalent to the response to a force dipole |Fij〉, has no components along floppy modes. Therefore the equation |δRij〉 =
|Fij〉 can be inverted. Using the spectral decomposition of
, one gets:
![]() | (14) |
![]() | (15) |
For comparison, the total amplitude of the mean response can be calculated from eqn (11), and one obtains for the norm square . Thus for d ≥ 2, the norm of the mean response vanishes as δz → 0, whereas the norm of the fluctuations diverge. Thus relative fluctuations must diverge at small ε, as observed in our data.
The most simple scenario is that the fluctuations of the response 〈δR(r)2〉 decays with the same characteristic length ∝ lc characterizing the mean response. Making this assumption, which is numerically verified (see below), and using the result of eqn (15), the angular average must read:
![]() | (16) |
To test this prediction, we fix all the nodes outside a circle of radius rf. A frozen boundary induces a rigid region where floppy modes are forbidden. We determine this region using the pebble algorithm,32 as shown in Fig. 5a. As rf decreases, all floppy modes eventually vanish. For any z we can define n*, the number of nodes involved in the last floppy mode, and for d = 2. Our measurements are shown in Fig. 5b and follow the prediction l* ∼ 1/δz. This result implies that exponentially small displacement at distances r > lc cannot be neglected when the rigid–floppy transition induced by freezing boundary conditions is considered. On the other hand, our results imply that floppy networks that are stabilized by pinning boundaries on the scale l* present soft modes with exponentially small frequencies near the center of the sub-system, since there are floppy modes with displacements of very tiny amplitude near the boundaries, of order exp(−l*/lc) ∼
.
![]() | ||
Fig. 5 Rigid regions, shown as red diamonds and magenta square, are induced by freezing the nodes outside an external radius (red) and an internal radius (magenta) respectively. Circular nodes (blue) show the minimal floppy region. Freezing any extra particle rigidifies the entire system. (b) l*, as defined in the text, vs. δz. The red line corresponds to l* ∼ 1/δz. |
In Fig. 6a the response to a local strain for a system with weak springs with kweak = 10−7 is shown. The deformation field clearly differs from the response of the same network with kweak = 0 (Fig. 6b). We now argue that the limit kweak → 0 is singular: for any kweak > 0, the response decays as a power-law at large distances r ≫ lc, even in the limit kweak → 0. The singularity of this limit can be seen by decomposing the response to a local strain in two parts. First, we consider the displacement of characteristic amplitude δR that would follow such a strain in the absence of weak interactions, i.e. with kweak = 0 (Fig. 6b). Second, we consider the additional displacement induced by the presence of weak springs (Fig. 6c). Indeed due to the weak springs forces are not balanced after step one, and forces of order Fweak ∼ kweakδR have appeared on the nodes. In the limit kweak → 0, these forces are vanishingly small and lead to no relaxation on the modes of non-vanishing frequency ω > ω*. However the spectrum now presents modes (stemming from the floppy modes that exist when kweak = 0) of characteristic frequency ωc. These modes relax due to the unbalanced weak forces with amplitude x that must satisfy mω2cx ≡ kweakx ∼ Fweak ∼ kweakδR, implying that x is independent of kweak. The convergence of the response δRt(r,kweak) in the small kweak limit is shown in Fig. 6d. We can define
![]() | ||
Fig. 6 (a) Example of the response to a local strain for δz = 0.2 and kweak = 10−7. (b) Response to the same strain in the absence of weak spring (kweak = 0), as studied in the last sections. (c) Difference between the displacement field of (a) and (b). (d) Average displacement δRtvs. the distance R from a local strain for floppy networks (δz = 0.2) with weak springs, as indicated in legend, with N = 90![]() |
Although localization is lost, the coordination continues to influence the response to a local strain (see Fig. 7a). In particular, near the imposed strain for r ≪ lc, the response is weakly influenced by the weak interaction, and decays with the same characteristic length, as shown in Fig. 7b. Note that we expect this scenario to hold also for gently compressed, hypostatic soft ellipses.§
![]() | ||
Fig. 7 (a) Average displacement δRtvs. the distance from a local strain r for floppy networks plus weak springs (see text) with N = 90![]() |
For isotropic floppy networks it was found numerically that , implying η ∼ 1/δz. Using our previous result on the response to a dipole in floppy materials and a simple hypothesis it is straightforward to derive this result, thus extending a previous derivation valid for z > zc,26 see also ref. 29. Here we perform this calculation for completeness, and because the present argument can be extended to predict that the correlation length under shear is lc. We consider an affine shear strain δγ applied on the network. After such a strain, unbalanced forces appear on the nodes: ∣Fδγ〉 = ∑〈ij〉γij∣Fij〉 where the sum is taken over all the bonds and ∣Fij〉 correspond to a dipole of force as defined in eqn (12). The coefficients γij are equal to nij·Γ·nij where Γ is the strain tensor, which is linear in δγ at first order. The displacement field can be written as a linear combination of responses to local dipoles ∣δRδγ〉 = ∑〈ij〉γij∣δRij〉, where ∣δRij〉 is defined in eqn (14). Two-point displacement correlations in space obey
![]() | (17) |
We now make the assumption that in random isotropic networks (as those considered here), the response of different dipoles is weakly correlated (this assumption turns out to be incorrect for flow of particles where subtle correlations are present in the structures visited by the dynamics). Using this assumption:
![]() | (18) |
Moreover, eqn (18) indicates that the correlation length in eqn (18) is essentially the length scale appearing in the response to a dipole lc. 〈δRij(x + r)·δRij(x)〉 must vanish when r is larger than the length lc where the response to a dipole is localized. On the other hand, correlations are not expected to vanish on a scale much smaller than lc, since the mean dipolar response will already give correlations on that scale:
For instance, the counting argument comparing surface effects and bulk degrees of freedom should hold, and l* is expected to characterize the effect of freezing boundaries.
We do not know what the spatial correlations of the velocities should be in flows, because in this case the velocity is dominated by the lowest frequency mode that is not captured by the present analysis.20,34 However the response to a local disturbance (such as the formation of a new contact) must couple predominantly to the modes above ω*, and is thus expected to affect the flow on a length lc. This can be shown using spectral decomposition for the response under a dipole of force ∣Fij〉. Decomposing the spectrum in two parts (above and below ω*) eqn (14) reads
![]() | (19) |
Numerical results show the appearance of a single mode below ω*,20 however a theoretical estimation in the large N limit supports that the distribution of modes between ωmin and ω* is given by Dmin(ω) ∝ ω(ω2 − ωmin2)(d−2)/2.34 Following the argument of Section 4.1, the average amplitude of the response is given by
The relative contribution follows Imin/I* ∼ δz(d−1) (with a logarithmic correction in two dimensions), indicating that in the limit δz ≪ 1 the contribution of the modes below ω* becomes vanishingly small and can be neglected. The modes above ω* in flow have the same density of states and are expected to have the same properties than those of isotropic networks. Thus the response to a local disturbance in flow must decay exponentially with a typical length lc, as could be tested empirically in two dimensional granular flows where imaging is possible, or perhaps using confocal imaging in slow emulsion flows. This simple argument does not hold for the velocity correlation in flow: the modes in the plateau do not contribute significantly to the response to shear, which is rather dominated by the lowest frequency modes.
For the strongly disordered floppy isotropic networks that we consider, one example of a question that remains to be explored is the evolution of elasticity (e.g. the shear modulus G) when boundaries are pinned at a distance L. For L > l*, the system remains floppy and G = 0. For L ≪ lc, one expects a mean-field argument to apply: pinning boundaries is equivalent to adding springs, leading to an increase of coordination Δz ∼ 1/L. For z > zc it is known that G ∼ z − zc,36 and thus we expect G ∼ 1/L. The behavior of G at intermediate length l* ≫ L ≫ lc remains to be explored.
Finally we compare our results with previous works in amorphous solids made of soft repulsive particles, for which z > zc. Near the unjamming transition where pressure vanishes the coordination approaches the Maxwell threshold from above (z → zc). A vanishing frequency scale ω* ∼ z − zc was predicted to characterize the low-frequency part of the spectrum,37,38 as confirmed numerically.39 The same theoretical argument indicated that boundaries affect elasticity on a length scale l* ∼ 1/(z − zc). It was later argued,40 based on numerical observations, that l* characterizes the response to a point force in packings of particles. Another length scale lc was observed numerically to characterize the response at a frequency ω*39 and to affect transport,41 behaviors that were well-captured by effective medium.24 Our work extends these results to floppy materials with z < zc, where these lengths and frequency scales characterize the phonon gap and the localization of floppy modes. However the comparison underlines an important disparity: for floppy networks we find both numerically and theoretically that lc characterize the response to a local force, whereas l* appears to characterize the response to a point force in packings.40 More numerical and theoretical investigations are needed to understand this difference.
Footnotes |
† The parameter ![]() |
‡ Tighe29 defines a length scale λf = f/F as the ratio between the typical contact force f carried by springs to the viscous forces F exerted by the fluid. F is proportional to the velocity of the particles, and therefore to their displacements time the frequency ω. It is thus readily extractable from our results. The contact forces f are such that the forces are balanced on each node, which implies on average that the spatial derivatives of the contact forces f are the external forces F. If one considers for example the response to a local perturbation, our result implies that 〈f〉/|〈F〉| is of the order of the length scale lc(ω) on which the mean elastic response decays. Thus this definition is consistent with our results. However a different result is obtained if one considers the ratio ![]() ![]() ![]() |
§ In the case of ellipses the pre-stress2 stabilizes the system, shifting the zero frequencies of floppy modes to finite values.9 We believe that the pre-stress plays a similar role to the weak springs of our network model. Accordingly, we expect the response to a local strain to be extended, even when the pre-stress is vanishingly small. Finally, note that our theoretical description of the vibrational spectrum does not capture the singular coupling that occurs between translational and rotational modes when the ellipticity α is small. We thus expect our approach to apply only for rather large values of α. |
This journal is © The Royal Society of Chemistry 2013 |