Phonon gap and localization lengths in floppy materials

Abstract

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). l_c 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 l_c and l^* have also been identified but where their roles are not fully understood.

---

1 Introduction

In 1864 Maxwell¹ showed that, in order to be mechanically stable, the average coordination value z of structures made of points connected by rigid bars must be larger than a threshold value z_c. Most common solids satisfy a microscopic version of this constraint.^2,3 If the constraint is violated, collective modes with no restoring forces (floppy modes) exist. These materials will be referred to as strictly floppy systems. However, in some cases these modes are stabilized by weak interactions: the bending energy of semi-flexible polymers confers a finite elasticity to gels,⁴ as do van der Waals interactions in weakly coordinated covalent glasses.⁵ Floppy modes can also be stabilized by the pre-stress applied on the system, as is the case for gently compressed packings of ellipsoid particles.^6–9 Fluids, on the other hand, can display modes that are strictly floppy. For example, in granular systems or suspensions of hard particles, large and sometimes percolating clusters of connected particles can be formed,^10,11 and motion within these clusters can only occur along floppy modes where no particles overlap. As the density increases toward jamming, the viscosity^12–14 and the length scale^15–17 characterizing the correlation of the dynamics diverge, up to the point where floppy modes disappear and the dynamics stops.

In simplified numerical models of ellipsoid particles,^7,8 covalent networks,¹⁸ gels of semi-flexible polymers¹⁹ and suspension flows,²⁰ 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.²⁰ An early work by Garboczi and Thorpe²¹ 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 methods^21–25 and numerical simulations, that strictly floppy elastic networks present a vibrational gap between floppy modes (zero frequency modes) and a frequency ω^* ∼ z_c − 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 ∼l_c. Even though l_c 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, l_c 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 l_c and l^* both in floppy materials and in jammed packings, which are discussed in the last section.

2 Model description and simulation

We consider a network of N point particles of mass m, connected by N_c harmonic un-stretched springs in spatial dimension d. The elastic energy following a deformation vector field {δRⁱ}_i=1…N is:where the index 〈ij〉 runs over the N_c springs, whose stiffnesses are k_ij and whose directions are along the unit vectors n_ij. Floppy modes correspond to displacements for which δE = 0. Since the energy is the sum of N_c positive definite terms, it vanishes only if all these terms are zero. Such modes always exist if the number of degrees of freedom N_d is larger than the number of constraints N_c, or equivalently if z ≡ 2N_c/N < z_c = 2d.

The force field ∣F〉 generated by displacements can be obtained by taking the derivative of eqn (1). One obtains ∣F〉 = [scr M, script letter M]∣δR〉, where:

is the stiffness matrix. We have used the bra-ket notation for which 〈i∣δR〉 = δR_i. The normal modes are given by the eigenvectors of [scr M, script letter M], and their frequencies are given by the square root of their associated eigenvalues. Herein, we will investigate the spectrum of [scr M, script letter M] and the spatial properties of its eigenvectors when z < z_c.

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 z_c. 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 [scr M, script letter M] numerically to compute the density of states D(ω) for various z < z_c. 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 → z_c. 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 → z_c. 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 000 nodes in two dimensions. The continuous line corresponds to the theoretical prediction of eqn (9). Inset: non-rescaled density of states D(ω) vs. ω. Floppy modes lead to a delta function at ω = 0 and are not presented.

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.

3 Effective medium theory (EMT)

3.1 General formalism

We use EMT, also known as coherent potential approximation,^22,23,27 to investigate the behavior of floppy networks as the one observed in the previous section. EMT is a mean field approximation and neglects spatial fluctuations in the coordination. Therefore, it does not properly describe rigidity percolation.²³ However, it was shown to describe materials for which spatial fluctuations of coordination are small.²⁴ EMT attempts to describe disordered materials as ordered materials with a frequency-dependent effective stiffness k̃_eff. We apply this technique to isotropic lattices of coordination z_in ≫ z_c, where bonds are then randomly removed with a probability (1 − p), so that the stiffness coefficients k_ij take the values 0 or k with a probability (1 − p) or p respectively, and the final coordination is z = pz_in. The Green's function G(ω) of the disordered floppy system is defined as [−mω² + [scr M, script letter M]]G(ω) = −1. G₀ is the Green's function of the effective medium, corresponding to the initial ordered lattice with an undetermined effective stiffness coefficient k̃_eff. Standard calculations lead to the Dyson relationship G = G₀ + G₀[scr T, script letter T]G₀,²⁸ where the operator [scr T, script letter T] can be expressed as an infinite series in terms of increasing numbers of interacting contacts:

The transfer matrix is found to be:

We seek an effective stiffness k̃_eff that captures the average behavior of the system, i.e. 〈G〉 = G₀, where the average is taken over the disorder on the stiffness coefficients k_ij. This condition leads to 〈[scr T, script letter T]〉 = 0. In the EMT this constraint is approximated by 〈T_〈ij〉〉 = 0. Using standard identities for the Green's function on isotropic lattices,²⁷ one can express this condition as:²³

where tr[·] stands for the trace, and k_eff ≡ k̃_eff/k.

Since we are interested in low frequencies , we approximate G₀ by its continuum limit and use a Debye cut-off q_D. We introduce the bulk modulus K and the shear modulus μ of the ordered lattice without an effective stiffness, and use the approximation:

The sum is taken over all polarization vectors q_α, where for pressure waves and for shear waves.

3.2 Scaling analysis near z = z_c

We now perform a scaling analysis of eqn (4) and (5) as z → z_c from below. Eqn (5) implies that , where the function f is independent of z and ω. As ω → 0, the elastic moduli of floppy systems vanish and k_eff → 0. Eqn (4) then leads to , where ε satisfies the equation

ε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ω²/k_eff∣ ≪ 1 and ∣k_eff∣ ≪ 1 (which can be shown to be true a posteriori in the limit of δz ≪ 1 and ) one finds:

For d = 2, eqn (5) leads to a logarithmic divergence in the small ε limit:

where Γ is a constant that depends on the elastic moduli and the Debye cut-off. The weak divergence of f does not modify the perturbation analysis performed for d ≥ 3, except that f(0) must now be replaced by f(ε) in eqn (7). Only f(ε) remains undetermined.

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 Λ(c_p,c_s,Γ).† To compare the theory and numerical simulations, we use Λ(c_p,c_s,Γ) 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:

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 Λ(c_p,c_s,Γ) = 1.3, which is fixed by this measurement.

3.3 Linear response to a local strain

The response to a local force G₀(ω,r_ij) diverges as ω → 0 due to the presence of floppy modes. One must rather consider the response field due to an imposed displacement. Changing the rest length of a spring placed between nodes i and k at a fixed frequency ω₀ corresponds to imposing a displacement (δR_i − δR_k)·n = e^iω₀t, with n being the unit vector along the connecting bond. The force required to impose such a displacement is

In the small frequency regime, the last expression is found to be |F〉 ∼ k_effe^iω₀tn(|i〉 − |k〉). The magnitude of the force vanishes as ω₀ → 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. rq_D ≫ 1), can be calculated using the continuum limit

〈δR(r)〉 = g(r − r_i) − g(r − r_k), (10)

where
and K₀ is the modified Bessel function of the second kind, that behaves exponentially at long distances. Thus, the mean response to a local strain decays exponentially in floppy systems with a characteristic length . The mean perturbation induced by changing the rest length of a spring has a quadrupolar symmetry shown in Fig. 3a. After averaging over 6000 realizations, we obtain a good agreement with our theoretical prediction, shown in Fig. 3b.

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 c_s and c_p 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 000 nodes and δz = 0.05.

The asymptotic solution in any dimension has an exponential decay. Indeed, taking the angular average one gets

where the over-line stands for angular average. In two dimensional networks, rescaling the distance by l_c and the amplitude by collapses the response to a local strain at different coordination values into a single curve, as shown in Fig. 4a and b.

Fig. 4 (a) Mean response δR_mvs. the distance r from the imposed strain with N = 40 000 nodes. (b) Rescaled average displacement of the mean response δR_mvs. the rescaled distance r/l_c, using Λ(c_p,c_s,Γ) = 1.3 as extracted from Fig. (1). (c) The fluctuations around the mean response are characterized by δR_t, defined in eqn (16). δR_t is plotted vs. the distance from a local strain r. (d) δR_tvs. the rescaled distance r/l_c, using Λ(c_p,c_s,Γ) = 1.3.

The zero frequency limit can be extended to finite frequencies by replacing ε → −mω²/k_eff in eqn (10). The asymptotic behavior is then given by log(〈δR(r)〉) ∼ −r/l_c(ω) + 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 l_c. (ii) For ω ≫ ω^*, the decay length decreases as l_c(ω) ∼ 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 ≥ z_c.²⁴

We can also consider the case of our networks immersed in a viscous fluid with viscosity η₀ and without hydrodynamic interactions. The response in the over-damped regime is calculated by replacing the inertial term mω² by the viscous force −iη₀ω in eqn (5) and then in the effective stiffness in eqn (7). In this limit the elastic modulus is proportional to k_eff, whose real and imaginary parts give the storage and loss moduli G′ and G′′ respectively. Two different regimes are observed. (i) For ω ≪ ω^*2m/η₀, the storage moduli follow G′ ∼ l_c⁶ω²/f(ε) and the loss moduli G′′ ∼ l_c²ω. The response is given by l_c(ω) ∼ l_c and v(ω) ∼ f(ε)l_c⁻³, displaying a similar exponential decay to that observed in an undamped system, but with propagating waves. (ii) For ω ≫ ω^*2m/η₀, both the storage and loss moduli are proportional to , as also observed in ref. 29. The decaying length l_c(ω) ∼ f(ε)^1/4ω^1/4, while the velocity v(ω) ∼ f(ε)^1/4ω^3/4. Our scaling for the characteristic length scale l_c(ω) differs from the recent result of Tighe.²⁹ In our opinion, the difference stems from an incorrect definition of the length scale characterizing the elastic response in floppy materials (see footnote‡).

4 Beyond the mean field

4.1 Fluctuations

The obvious difference between the mean response to a local strain (Fig. 3) and a typical one (Fig. 2) indicates large fluctuations. A priori mean field models as the effective medium do not enable to capture those. However, it is possible to combine the EM results with additional considerations to estimate the amplitude of these fluctuations.

We denote the dipole of forces ∣F_ij〉 generated by changing the rest length of the spring ij by a distance one:

|F_ij〉 = n_ij(|i〉 − |j〉), (12)

where the stiffness coefficient k is set to unity. From the definition of the stiffness matrix, eqn (2) and (12), we can write:where the sum is taken over all the bonds. Note that any floppy mode |δR₀〉 has by definition no restoring force [scr M, script letter M]|δR₀〉 = 0, implying that 〈δR₀|F_ij〉 = 0 for all contacts ij.

Thus, the response |δR_ij〉 to the elongation of a spring ij, which is equivalent to the response to a force dipole |F_ij〉, has no components along floppy modes. Therefore the equation [scr M, script letter M]|δR_ij〉 =  |F_ij〉 can be inverted. Using the spectral decomposition of [scr M, script letter M], one gets:

where ω² and ∣δR_ω〉 are the non-zero eigenvalues and the corresponding eigenvectors of the stiffness matrix [scr M, script letter M].^20,26Eqn (14) implies that the norm of the response follows 〈δR_ij|δR_ij〉 = ∑_ω>0〈δR_ω|F_ij〉²/ω⁴. Introducing the average amplitude of the response , one gets using eqn (13):

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)²〉 decays with the same characteristic length ∝ l_c 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:

where log(h(x)) ∼ −x for x ≫ 1, and the ε dependence is determined by eqn (15). Rescaling the distance by l_c leads to a very good collapse of the response at different coordination values, as shown in Fig. 4c and d. Note that in two dimensions the amplitude does not need to be rescaled.

4.2 Pinning boundaries

We now turn to the spatial properties of floppy modes. The response to the stretch of a spring, exemplified in Fig. (2), is a floppy mode of the network where this spring is removed. We have thus shown that floppy modes can be localized (in the sense of presenting an exponential decay) on a length scale (with a log correction for d = 2). We now extend a previous counting argument,^3,30 tested for packings with z > z_c in ref. 31, to floppy networks (see also ref. 29) and show that floppy modes are also characterized by another length scale l^* ∼ 1/δz ≫ l_c. Below l^*strict localization is impossible: floppy modes cannot have a smaller compact support. Essentially, l^* is the length at which the number of constraints r^d−1_f that result from freezing the boundary of a system of size r_f is equal to the number of floppy modes in the bulk δzr^d_f. This counting argument suggests that if the boundaries are frozen on a scale l^* or smaller, floppy modes must vanish.

To test this prediction, we fix all the nodes outside a circle of radius r_f. A frozen boundary induces a rigid region where floppy modes are forbidden. We determine this region using the pebble algorithm,³² as shown in Fig. 5a. As r_f 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 > l_c 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^*/l_c) ∼ .

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.

5 Weak interactions

As discussed in the Introduction, in some materials, floppy modes are stabilized by weak interaction, as is the case in covalent glasses. In order to model these weak interactions, we consider (see for example ref. 26) floppy networks with k = 1 and add a number of weak springs of stiffness k_weak ≪ mω^*2, which gives a finite elasticity to the network. Floppy modes then gain finite frequencies, of order the characteristic frequency scale associated with the weak interaction . A gap in the density of states remain present in the frequency range [ω_c,ω^*], if the weak interaction is weak enough, i.e.. In this limit the elastic moduli scale as μ_weak ∼ k_weak/δz,²⁶ where δz is the excess coordination of the network of strong interaction. At long enough wavelengths the material, to good approximation, must behave as a continuous elastic medium, and for ω ≪ ω_c the density of states must follow a Debye behavior D(ω) ∼ ω^d−1. Extracting a velocity of sound from the elastic moduli and computing the wavelength at ω ∼ ω_c, one obtains a wavelength of order l_c, which thus characterizes as well the length scale above which a continuum description becomes a good approximation.

In Fig. 6a the response to a local strain for a system with weak springs with k_weak = 10⁻⁷ is shown. The deformation field clearly differs from the response of the same network with k_weak = 0 (Fig. 6b). We now argue that the limit k_weak → 0 is singular: for any k_weak > 0, the response decays as a power-law at large distances r ≫ l_c, even in the limit k_weak → 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 k_weak = 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 F_weak ∼ k_weakδR have appeared on the nodes. In the limit k_weak → 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 k_weak = 0) of characteristic frequency ω_c. These modes relax due to the unbalanced weak forces with amplitude x that must satisfy mω²_cx ≡ k_weakx ∼ F_weak ∼ k_weakδR, implying that x is independent of k_weak. The convergence of the response δR_t(r,k_weak) in the small k_weak limit is shown in Fig. 6d. We can define

to quantify the difference between the response of the networks with weak springs and the ones strictly floppy. In the inset of Fig. 6d one can clearly see that C_weak converges to a non-zero value showing the singularity of this limit.

Fig. 6 (a) Example of the response to a local strain for δz = 0.2 and k_weak = 10⁻⁷. (b) Response to the same strain in the absence of weak spring (k_weak = 0), as studied in the last sections. (c) Difference between the displacement field of (a) and (b). (d) Average displacement δR_tvs. the distance R from a local strain for floppy networks (δz = 0.2) with weak springs, as indicated in legend, with N = 90 000 nodes. Inset: C_weakvs. k_weak (see text for definition), where L₀ has been taken to be 80 spring lengths.

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 ≪ l_c, 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 δR_tvs. the distance from a local strain r for floppy networks plus weak springs (see text) with N = 90 000 nodes and k_weak/k = 10⁻⁷. (b) Average displacement δR_tvs. the rescaled distance r/l_c, using Λ(c_p,c_s,Γ) = 1.3. The response of strictly floppy networks (filled symbols) has been included for comparison. Inset: the initial decay magnified.

6 Rheology of dense suspensions

6.1 Random networks under a global shear

It is interesting to consider the behavior of the present networks under shear, if they were placed in a viscous solvent. This problem is formally related²⁰ to the affine solvent model of suspension flow of hard particles, where hydrodynamic interactions between particles are neglected. We consider that our networks are made of point particles, subjected to a viscous drag F proportional to the particles velocities V in the reference frame of the solvent, which is assumed to follow an affine shear. One can show, from the definition of this model, that the viscosity η is proportional to the non-affine velocity squared of the particles.^20,29,33 The viscosity is proportional to the ratio P/[small gamma, Greek, dot above]², where P is the total power dissipated and [small gamma, Greek, dot above] is the strain rate. The power dissipated follows P ∼ 〈F·V〉 ∼ ‖V‖². It is convenient to write the velocity field in terms of the non-affine displacement following an infinitesimal strain δR/δγ, i.e.V = [small gamma, Greek, dot above]δR/δγ. One can find that η ∼ 〈(δR/δγ)²〉, which relates the viscosity to the non-affine response to shear.

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 > z_c,²⁶ 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 l_c. We consider an affine shear strain δγ applied on the network. After such a strain, unbalanced forces appear on the nodes: ∣F_δγ〉 = ∑_〈ij〉γ_ij∣F_ij〉 where the sum is taken over all the bonds and ∣F_ij〉 correspond to a dipole of force as defined in eqn (12). The coefficients γ_ij are equal to n_ij·Γ·n_ij 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∣δR_ij〉, where ∣δR_ij〉 is defined in eqn (14). Two-point displacement correlations in space obey

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:

where we used γ_ij² ∼ δγ². Combining with eqn (15) one finds
implying that the viscosity diverges as η ∼ C_δγ(0)/δγ² ∼ ε⁻¹.

Moreover, eqn (18) indicates that the correlation length in eqn (18) is essentially the length scale appearing in the response to a dipole l_c. 〈δR_ij(x + r)·δR_ij(x)〉 must vanish when r is larger than the length l_c where the response to a dipole is localized. On the other hand, correlations are not expected to vanish on a scale much smaller than l_c, since the mean dipolar response will already give correlations on that scale:

6.2 ASM model of suspension flows

We predict a plateau of modes above a frequency ω^* ∼ δz. This plateau also appears in models of suspension flows.²⁰ This observation supports that the plateau is generically present in floppy materials, and that the frequency response above ω^* will be well described by the effective medium predictions derived here. There is, however, a crucial difference between the spectra of our networks and those of shear flows: in the latter one mode appears at a frequency ω_min ≪ ω^* that strongly couples to a global shear and leads to a divergence of the viscosity η ∼ 1/δz^2.85,²⁰ which is sharper than in isotropic networks where η ∼ 1/δz. Therefore, typical configurations visited in dense flows differ qualitatively from our random networks, for reasons explained in ref. 34. Nevertheless, we expect some of our present predictions to hold in flow.

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 l_c. This can be shown using spectral decomposition for the response under a dipole of force ∣F_ij〉. Decomposing the spectrum in two parts (above and below ω^*) eqn (14) reads

Numerical results show the appearance of a single mode below ω^*,²⁰ however a theoretical estimation in the large N limit supports that the distribution of modes between ω_min and ω^* is given by D_min(ω) ∝ ω(ω² − ω_min²)^(d−2)/2.³⁴ Following the argument of Section 4.1, the average amplitude of the response is given by

The relative contribution follows I_min/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 l_c, 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.

7 Discussion and open questions

We have argued that the elasticity of floppy networks is characterized by a gap in the vibrational spectrum, and by two length scales l_c and l^* that diverge near jamming (i.e. δz → 0). The existence of two lengths raises the question of which properties are governed by which scale. Our work supports that the length scale characterizing the response to imposed forces, and to most standard observations, is l_c. On the other hand l^* characterizes the effect of pinning boundaries. Such effects are subtle, and can depend on surprising ways on the type of elastic networks considered.³⁵

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 ≪ l_c, 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 > z_c it is known that G ∼ z − z_c,³⁶ and thus we expect G ∼ 1/L. The behavior of G at intermediate length l^* ≫ L ≫ l_c remains to be explored.

Finally we compare our results with previous works in amorphous solids made of soft repulsive particles, for which z > z_c. Near the unjamming transition where pressure vanishes the coordination approaches the Maxwell threshold from above (z → z_c). A vanishing frequency scale ω^* ∼ z − z_c was predicted to characterize the low-frequency part of the spectrum,^37,38 as confirmed numerically.³⁹ The same theoretical argument indicated that boundaries affect elasticity on a length scale l^* ∼ 1/(z − z_c). It was later argued,⁴⁰ based on numerical observations, that l^* characterizes the response to a point force in packings of particles. Another length scale l_c was observed numerically to characterize the response at a frequency ω^*39 and to affect transport,⁴¹ behaviors that were well-captured by effective medium.²⁴ Our work extends these results to floppy materials with z < z_c, 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 l_c characterize the response to a local force, whereas l^* appears to characterize the response to a point force in packings.⁴⁰ More numerical and theoretical investigations are needed to understand this difference.

Acknowledgements

We thank Y. Elmatad, C. Falcón, P. Rowghanian, A. Grosberg and P. Chaikin for comments on the manuscript. This work has been supported by the Sloan Fellowship, NSF DMR-1105387, and Petroleum Research Fund #52031-DNI9. This work was also supported partially by the MRSEC Program of the National Science Foundation under Award Number DMR-0820341.

References

1. J. Maxwell, Philos. Mag., 1864, 27, 294–299.
2. S. Alexander, Phys. Rep., 1998, 296, 65–236.
3. M. Wyart, Ann. Phys., 2005, 30(3), 1.
4. C. Heussinger and E. Frey, Phys. Rev. Lett., 2006, 97, 105501.
5. J. C. Phillips and M. F. Thorpe, Solid State Commun., 1985, 53, 699–702.
6. A. Donev, I. Cisse, D. Sachs, E. A. Variano, F. H. Stillinger, R. Connelly, S. Torquato and P. M. Chaikin, Science, 2004, 303, 990–993.
7. M. Mailman, C. F. Schreck, C. S. O'Hern and B. Chakraborty, Phys. Rev. Lett., 2009, 102, 255501.
8. Z. Zeravcic, N. Xu, A. J. Liu, S. R. Nagel and W. van Saarloos, Europhys. Lett., 2009, 87, 26001.
9. C. F. Schreck, M. Mailman, B. Chakraborty and C. S. O'Hern, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 061305.
10. X. Cheng, J. H. McCoy, J. N. Israelachvili and I. Cohen, Science, 2011, 333, 1276–1279.
11. D. Bonamy, F. Daviaud, L. Laurent, M. Bonetti and J. P. Bouchaud, Phys. Rev. Lett., 2002, 89, 034301.
12. E. Brown and H. Jaeger, Phys. Rev. Lett., 2009, 103, 086001.
13. K. N. Nordstrom, E. Verneuil, P. E. Arratia, A. Basu, Z. Zhang, A. G. Yodh, J. P. Gollub and D. J. Durian, Phys. Rev. Lett., 2010, 105, 175701.
14. F. Boyer, E. Guazzelli and O. Pouliquen, Phys. Rev. Lett., 2011, 107, 188301.
15. O. Pouliquen, Phys. Rev. Lett., 2004, 93, 248001.
16. P. Olsson and S. Teitel, Phys. Rev. Lett., 2007, 99, 178001.
17. R. Lespiat, S. Cohen-Addad and R. Höhler, Phys. Rev. Lett., 2011, 106, 148302.
18. K. O. Trachenko, M. T. Dove, M. J. Harris and V. Heine, J. Phys.: Condens. Matter, 2000, 12, 8041–8064.
19. E. M. Huisman and T. C. Lubensky, Phys. Rev. Lett., 2011, 106, 088301.
20. E. Lerner, G. Düring and M. Wyart, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 4798–4803.
21. E. J. Garboczi and M. F. Thorpe, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 32, 4513.
22. X. Mao, N. Xu and T. C. Lubensky, Phys. Rev. Lett., 2010, 104, 085504.
23. S. Feng, M. F. Thorpe and E. Garboczi, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 31, 276–280.
24. M. Wyart, Europhys. Lett., 2010, 89, 64001.
25. C. P. Broedersz, X. Mao, T. C. Lubensky and F. C. MacKintosh, Nat. Phys., 2011, 7, 983–988.
26. M. Wyart, H. Liang, A. Kabla and L. Mahadevan, Phys. Rev. Lett., 2008, 101, 215501.
27. S. Kirkpatrick, Rev. Mod. Phys., 1973, 45, 574–588.
28. I. Webman, Phys. Rev. Lett., 1981, 47, 1496–1499.
29. B. P. Tighe, ArXiv e-prints, 2012, 1203.3411.
30. M. Wyart, S. Nagel and T. Witten, Europhys. Lett., 2005, 72, 486–492.
31. M. Mailman and B. Chakraborty, J. Stat. Mech.: Theory Exp., 2012, P05001.
32. D. J. Jacobs and B. Hendrickson, J. Comput. Phys., 1997, 137, 346–365.
33. B. Andreotti, J.-L. Barrat and C. Heussinger, ArXiv e-prints, 2011, 1112.1194.
34. E. Lerner, G. Düring and M. Wyart, ArXiv e-prints, 2011, 1201.3650.
35. C. F. Moukarzel, Europhys. Lett., 2012, 97, 36008.
36. A. J. Liu, S. R. Nagel, W. van Saarloos and M. Wyart, Dynamical Heterogeneities in Glasses, Colloids, and Granular media, Oxford, 2010.
37. M. Wyart, S. R. Nagel and T. A. Witten, Europhys. Lett., 2005, 72, 486.
38. M. Wyart, L. E. Silbert, S. R. Nagel and T. A. Witten, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 051306.
39. L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2005, 95, 098301.
40. W. G. Ellenbroek, E. Somfai, M. van Hecke and W. van Saarloos, Phys. Rev. Lett., 2006, 97, 258001.
41. V. Vitelli, N. Xu, M. Wyart, A. J. Liu and S. R. Nagel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 021301.

---

Footnotes

† The parameter

‡ Tighe²⁹ 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 l_c(ω) 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 , as done by Tighe. These two quantities differ because the mean response is much smaller than the fluctuations around it (see Section 4.1). Therefore, associating the quantity with a length appears unjustified. In the zero frequency limit we can actually calculate this last expression using the formalism developed in ref. 20. For isotropic random floppy networks one gets 〈f²〉 ∼ ω^*−3 and 〈F²〉 ∼ ω^*−1, leading to a ratio consistent with the numerical results of Tighe.

§ In the case of ellipses the pre-stress² stabilizes the system, shifting the zero frequencies of floppy modes to finite values.⁹ 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 α.

---