Open Access Article

This Open Access Article is licensed under a

Creative Commons Attribution 3.0 Unported Licence

Karsten
Baumgarten
* and
Brian P.
Tighe

Delft University of Technology, Process & Energy Laboratory, Leeghwaterstraat 39, 2628 CB Delft, The Netherlands. E-mail: k.baumgarten@tudelft.nl

Received
11th August 2017
, Accepted 9th October 2017

First published on 11th October 2017

When weakly jammed packings of soft, viscous, non-Brownian spheres are probed mechanically, they respond with a complex admixture of elastic and viscous effects. While many of these effects are understood for specific, approximate models of the particles' interactions, there are a number of proposed force laws in the literature, especially for viscous interactions. We numerically measure the complex shear modulus G* of jammed packings for various viscous force laws that damp relative velocities between pairs of contacting particles or between a particle and the continuous fluid phase. We find a surprising sensitive dependence of G* on the viscous force law: the system may or may not display dynamic critical scaling, and the exponents describing how G* scales with frequency can change. We show that this sensitivity is closely linked to manner in which viscous damping couples to floppy-like, non-affine motion, which is prominent near jamming.

Linear viscoelasticity is characterized by the frequency dependent complex shear modulus G*(ω) = G′(ω) + ıG′′(ω); its real and imaginary parts are known as the storage and loss modulus, respectively, and quantify the amount of energy stored elastically and dissipated viscously during one cycle of oscillatory driving at angular frequency ω.^{18} The form of the complex shear modulus near jamming was first described by Tighe^{13} for a system of soft spheres interacting via “one-sided” (purely repulsive) springs and linear viscous contact forces; details of the model are presented below. Characteristic features can be seen in Fig. 1, which plots the average G* for states prepared close to jamming. At both low and high frequencies, the storage modulus (filled symbols) and loss modulus (open symbols) resemble a simple Kelvin–Voigt solid (a spring and dashpot in parallel),^{18} with G′ ∼ const and G′′ ∼ ω. There is also a critical regime at intermediate frequencies, in which both G′ and G′′ scale as ω^{1/2}. Similar square root scaling has been observed experimentally in foams, emulsions, and other complex fluids,^{19–23} and has been linked theoretically to strongly non-affine motion.^{19} Plots of the particles' displacements from a static initial condition, evaluated at zero and peak shear stress (Fig. 1, bottom six panels) show that the critical regime represents a broad crossover from highly non-affine motion in the quasistatic limit at vanishing ω, to strongly affine motion at high frequencies.

The square root scaling is anomalous, in the sense that simple linear interactions at the particle scale give rise to nonlinear frequency dependence in the bulk. In contrast, the frequency dependence of G* in a Kelvin–Voigt solid is consistent with a direct extrapolation from the elastic forces (linear in the particle displacements) and viscous forces (linear in the velocities). Moreover, in soft spheres the critical regime broadens on approach to the jamming transition, with its lower bound approaching zero as the confining pressure p goes to zero and the system unjams.^{13} This strongly suggests that critical effects lie at the origin of the square root scaling near jamming.

While spring-like forces are standard in numerical models of foams and emulsions,^{1,3–7,9,13,24} a number of alternate proposals for viscous interactions can be found in the literature.^{1,25–28} This variety is largely due to authors' efforts to strike a balance between physical accuracy and computational complexity. What influence does the viscous force law have on bulk viscoelastic response near jamming? In equilibrium systems near a critical point, growing correlations wash out particle-scale details, so that similar scaling in bulk properties can be found for different interparticle interactions.^{29} Here we show that the nonequilibrium jamming transition is different: the complex shear modulus near jamming is surprisingly sensitive to the form of the viscous force law. Seemingly similar choices can alter the apparent scaling exponents or eliminate dynamic critical scaling entirely. Still others lead to subtler changes in the form of correlation functions.

To probe the role of viscous damping in viscoelasticity near jamming, we implement computer simulations of Durian's bubble model, a widely studied numerical model for foams and emulsions near ϕ_{c}. We investigate linear contact damping for varying ratios of the drag coefficients for normal and transverse motion, Stokes-like drag laws, and finally nonlinear damping of the relative velocities. One of our main conclusions will be to relate floppy-like, non-affine motion in the quasistatic limit to the form of the storage and loss moduli at finite frequency. We further study the role of two-point velocity correlations and effect of pre-stress on the dynamic viscosity.

^{el}_{i} + ^{visc}_{i} = . | (1) |

We consider ensembles of packings of N particles in D = 2 spatial dimensions prepared at a target pressure p. N = 32768 unless indicated otherwise. Initial conditions are generated by minimizing the total elastic potential energy using a nonlinear conjugate gradient algorithm, starting from particle positions placed randomly via a Poisson point process. As is typical in studies of jamming,^{30} the packings are bidisperse to avoid crystallization, with equal numbers of large and small particles and a radius ratio 1.4:1. The systems are bi-periodic, and shear is imposed via Lees–Edwards boundary conditions.

Units are set by the mean particle size d, the particle stiffness k, and a microscopic time scale τ_{1} (the latter two being introduced below). In simulations all three are set to one. However, in some cases we include the microscopic time scale in scaling relations in order to emphasize the dimensionful or dimensionless character of a relation.

All our simulations are performed in D = 2 spatial dimensions, which is the upper critical dimension for the jamming transition.^{31} We therefore expect the critical behavior we describe here, and in particular the values of critical exponents, to remain unchanged for D > 2.

(2) |

For later convenience we note that the elastic energy corresponding to eqn (2) is , where

(3) |

(4) |

Linear contact damping.
We will explore a class of linear viscous contact force laws that damp relative velocities at the point of contact,

See Fig. 2a for an illustration. The quantity is the tangential velocity at the contact and ẑ is the out-of-plane unit vector. θ_{i} is the angular displacement of particle i from its orientation in the initial condition. Dots indicate differentiation with respect to time. The coefficient kτ_{1} controls the damping of relative normal motions. It is defined in terms of a microscopic time scale τ_{1}, which describes the exponential relaxation of two overlapping disks and sets the natural unit of time. The damping coefficient for relative transverse motion βkτ_{1} is defined by its ratio β to the damping coefficient for normal motion.

The Rayleigh dissipation function is used to implement linear damping forces in a Lagrangian formalism. Just as conservative forces are proportional to gradients of the potential energy, dissipative forces are proportional to gradients of the dissipation function.

(5) |

Fig. 2 (a) Relative velocities in the rest frame of particle i. (b) Motion with respect to an affine flow field. |

The case β = 1 describes equal damping of normal and transverse motion. For brevity we refer to this case as “balanced” contact damping. Examples of prior studies employing balanced contact damping include ref. 2, 4, 6–8, 13, 40 and 41. Note that some of these studies apply damping to the relative motion of the particles' centers, neglecting particle rotations. We include rotations, as this seems more physical – however, we have also implemented balanced damping without rotations and find the form of G* qualitatively unchanged from the results presented below.

We also separately consider the case β = 0, in which transverse motion goes undamped. This is not a physically realistic scenario for densely packed foams and emulsions. Nevertheless, this damping law is found in the literature, presumably because it exerts no torque, eliminating the need to keep track of rotational degrees of freedom.^{42–45} In dilute systems with volume fractions outside the range considered here, this same force law is also a means to implement inelastic collisions.

Finally, we also treat the case of arbitrary β. We are not aware of any prior work that has systematically varied this coefficient.

Again for later convenience, we note that the Rayleigh dissipation function corresponding to Eqn (5) is

(6) |

Stokes-like drag forces.
In addition to linear contact drag, we also consider a class of linear viscous force laws in which drag enters as a body force reminiscent of Stokes drag.^{1,46} These can be motivated in two ways.

In this interpretation, the damping coefficient kτ_{1} should be proportional to the fluid viscosity η_{F}, as specified in Stokes' law. The dissipation function is

The second term accounts for dissipation due to shearing of the continuous fluid phase.

In the first interpretation, drag between particles is neglected entirely. Instead drag is assumed to result from the motion of individual particles with respect to the continuous fluid phase, which itself is assumed to flow with an affine velocity profile _{aff}() = y set by the shear rate . A particle at position _{i} then experiences a drag force proportional to the difference between its velocity _{i} and the affine profile (see Fig. 2b),

(7) |

(8) |

An alternative interpretation of eqn (7) known as “mean field drag” was introduced by Durian.^{1} In this view the body force is an approximation to balanced contact damping. One assumes that the velocity of each contacting particle j can be replaced with its average value at that position, which coincides with the affine velocity field. Angular velocities are set to zero. The resulting viscous force law and dissipation function are identical to eqn (7) and (8), with the caveat that kτ_{1} no longer has a fixed proportion to the fluid viscosity. Retaining the fluid viscosity term in the dissipation function is advisable, however, as otherwise the system could deform affinely without dissipating energy.

Regardless of how the Stokes-like drag force is motivated, its advantage is again computational. As the equations of motion in the bubble model are overdamped, they are first order linear differential equations. Generally, these must be solved using matrix inversion (see below). However in the special case of eqn (7), the relevant inversion can be performed by hand. Prior studies using Stokes or mean field drag include ref. 1–3, 5, 17, 28, 41 and 47–49.

Nonlinear contact forces.
The viscous contact force law of eqn (5) is linear in the particle velocities. However, viscous friction laws in real foams are actually nonlinear in the relative velocities. There are two classes of interactions, associated with so-called mobile and immobile surfactants, which give rise to different flow profiles within the thin films of the flow, and therefore dissipate energy differently. The case of immobile surfactants was treated by Bretherton,^{25} whose drag law proportional to the 2/3 power of velocity was subsequently verified experimentally.^{50} More recently, Denkov and co-workers have argued for an exponent 1/2 in the case of mobile surfactants.^{26} Seth et al.^{27} have also suggested a nonlinear force law with exponent 1/2 to account for elastohydrodynamic interactions between deformable particles in soft glassy matter. We therefore consider force laws of the form

where is the relative velocity at the contact. The constant ρ_{0} has units of length and is required for dimensional consistency. We set it to 1.

(9) |

(10) |

(11) |

The Fourier transform of eqn (11) gives

(12) |

The complex shear modulus can be determined by solving eqn (12) for the complex vector |Q*(ω)〉 using standard linear algebra routines. The resulting shear strain is γ*(ω) = 〈|Q*(ω)〉. The complex shear modulus is then

(13) |

Inspired by the above observation, we now restrict our focus to frequencies ω < 1. A more rigorous derivation of the following results is found in ref. 13. Our approach here is more heuristic and begins with the scaling ansatz

(14) |

(15) |

(16) |

(17) |

The scaling ansatz of eqn (14)–(17) is tested in Fig. 3b, which plots G′/p^{μ} and G′′/p^{μ} with μ = 1/2 versus ω/p^{λ} with λ = 1. The resulting collapse is excellent. As expected, the real and imaginary parts of the scaling function are constant and linear, respectively, for low values of the rescaled frequency. There is a crossover around to a power law with exponent Δ ≈ 0.5 (long dashed line). This is the ω^{1/2} scaling discussed above.

The scaling collapse in Fig. 3 empirically determines the values of the critical exponents; they are μ = 1/2, λ = 1, and Δ = 1/2. The value of μ is fixed by the known scaling of G_{0}. The exponent Δ is related to μ and λ. To see this, note that one generally expects the moduli to remain finite except possibly at the critical point, where both p and ω go to zero. In the case where p = 0 and ω > 0, eqn (14)–(17) predict that both moduli scale as p^{μ−λΔ}ω^{Δ}, which remains finite only if Δ = μ/λ = 1/2. It remains to motivate λ = 1, which we do in Section 3.3.

It is useful first to consider response in the absence of the pre-stress term, i.e. by setting f^{0}_{ij} = 0 in eqn (4). The Hessian and damping matrices are then directly proportional, , allowing eqn (12) to be solved exactly in terms of G_{0},

G* = G_{0}(1 + ιτ_{0}ω). | (18) |

We emphasize that a seemingly simple change to the viscous force law, namely setting the damping coefficient for sliding motion to zero, has produced a dramatic and qualitative shift in the viscoelastic response. More precisely, the intermediate regime, identified above when β = 1, has completely vanished. Recall that this regime is a manifestation of dynamic critical scaling and dominates the response for a wide range of frequencies near jamming. In this sense setting β = 0 kills dynamic critical scaling.

What happens for intermediate values of β? In Fig. 4b we plot G* for fixed p and a range of β over seven decades. One sees that the critical regime gradually appears, and for sufficiently large β the moduli resemble their form for β = 1. This suggests that it is reasonable to speak of weakly and strongly damped sliding motion. We quantify this distinction more precisely below.

Packings at the jamming point are isostatic, meaning they have just enough contacts to constrain all particle motions (except for a few individual “rattlers”, which can be removed from the analysis). Consider breaking a contact in a packing at the jamming transition, where all contacting particles are “kissing” and f^{0}_{ij} = 0. The broken contact removes a constraint and therefore introduces a floppy mode, an infinitesimal motion of the particles that can be performed without work. By considering the energy expansion of eqn (4), one sees that all relative normal motions in a floppy mode must be zero – floppy motions are sliding motions, in which all relative motion between particles is transverse to the contact. Jammed packings do not have floppy modes, but the eigenmodes of the Hessian remain “floppy-like,” i.e. transverse/sliding motion dominates.^{12,51} This feature is also found in the response to shear, which is dominated by low frequency modes.^{13} Through careful analysis of the modes, it is possible to show that the shear modulus scales as G_{0} ∼ p^{1/2}.^{13,52} Here we take this scaling relation as a given and, following ref. 51, infer its consequences for the typical relative normal and transverse displacement amplitudes, Δu^{‖} and Δu^{⊥}, as well as the typical amplitude of non-affine displacements u_{na}.

By definition, the change in elastic energy ΔU ≡ U − U_{0} due to an infinitesimal shear strain γ is ΔU = (1/2)G_{0}Vγ^{2}. Momentarily neglecting the pre-stress term in eqn (4), which should be small as p → 0, we anticipate that the typical relative normal motion scales as (Δu^{‖})^{2} ∼ G_{0}γ^{2}, or

(19) |

(20) |

Finally, we consider the typical amplitude of non-affine displacements u_{na}. Bond vectors Δ^{0}_{ij} in the reference packing are randomly oriented, so there is a local competition between energetically favorable sliding at the particle scale, and globally imposed affine motion. Therefore we expect the typical non-affine amplitude to be comparable to the typical relative displacement amplitude, which is dominated by transverse motion, i.e.

(21) |

Eqn (19)–(21) have previously been derived and tested numerically by Ellenbroek et al. and Wyart et al.^{51,54} For completeness we verify them again in Fig. 5, which plots the median of the probability density function of |Δu_{‖}|, |Δu_{⊥}|, and |u_{na}| for varying p while neglecting the pre-stress term. Results including pre-stress show compatible trends, albeit with more noise; we revisit the role of pre-stress below. Plots of the means show the same trend for Δu_{⊥} and u_{na}, but Δu_{‖} develops a plateau at low p due to a long tail of the PDF.

Fig. 5 Scaling of the relative normal, relative transverse, and non-affine motion as a function of pressure. Dashed lines have slopes of ±1/4. |

We now use the quasistatic relations (19)–(21) to determine the critical exponent λ. The ω → 0 limit of the dissipation function is proportional to the dynamic viscosity, R = η_{0}(ωγ_{0})^{2}V/2, where γ_{0} is the maximum strain amplitude. At the same time, from the viscous force law one anticipates R ∼ ^{visc}·Δ ∼ ω^{2}[(Δu^{‖})^{2} + β(Δu^{⊥})^{2}]. Invoking eqn (19) and (20) gives

(22) |

Eqn (22) is compatible with our numerical results for undamped sliding (β = 0), as well. Then only the first term is present and η_{0} ∼ p^{1/2} – it vanishes rather than diverges.

One can also consider the case of arbitrary β. The second term will always dominate for sufficiently low pressure; hence the dynamic viscosity diverges for any finite damping of sliding motion. In this sense the case β = 0 is singular. For arbitrary β > 0 the crossover frequency where the quasistatic regime ends and the linear regime begins scales as ω* ∼ p/β. We have seen above that the critical regime ends at a frequency . Hence the critical regime, with its ω^{1/2} scaling in G′ and G′′, is avoided entirely whenever ω* ≫ 1, or β ≪ β* ∼ p. This crossover is evident in Fig. 4b. The scale β* provides a convenient dividing line between cases of strong and weak damping of transverse motion.

Fig. 6 (a) Critical scaling collapse of the storage and loss moduli for Stokes drag and fluid viscosity η_{F} = 1. (b) Storage and loss moduli for p = 10^{−4} and varying η_{F}. |

As in the previous section, the above result can be rationalized on the basis of quasistatic scaling relations. The key observation is that the typical non-affine motion u_{na} and the relative transverse motion Δu^{⊥} diverge in the same way as the pressure tends to zero; cf.Eqn (20) and (21). For Stokes drag the dissipation function scales as R ∼ (u_{na}ω)^{2}, again giving η_{0} ∼ 1/p^{1/2}.

Recall that if one considers the Stokes drag term to be a mean field approximation for balanced contact damping, then the fluid viscosity η_{F} can vary independently of the damping coefficient kτ_{1}. We probe the dependence of G* on η_{F} in Fig. 6b by varying η_{F} over ten decades. We observe that the fluid viscosity contributes a linear term η_{F}ω to the loss modulus, which is always dominant at sufficiently high frequencies. For large η_{F} and/or low pressures satisfying η_{F} ≫ 1/p^{1/2}, the loss modulus becomes linear for all frequencies. In this event the critical properties of the loss modulus are obscured, but criticality is still apparent in the storage modulus.

For a system undergoing simple shear in the x-direction, correlations of the non-affine displacements between particles separated by a distance δ_{ij} = |x_{i} − x_{j}| can be quantified with the two-point correlation function . Here u_{i,y}′ is the y-component of the real part of the complex displacement vector of particle i with x-coordinate x_{i}. The average 〈·〉 runs over all particle pairs within a narrow “lane”, hence C is a function of |δ_{ij}|. We have verified that C becomes independent of the lane width for sufficiently small values. We have also confirmed that results using the imaginary part are indistinguishable.

Non-affine correlation functions in weakly jammed solids have been studied previously for three cases. DiDonna and Lubensky^{55} and Maloney^{56} showed there is no characteristic length scale in quasistatic linear elastic response; instead C collapses when distances are rescaled by the box size L. Heussinger and Barrat found compatible results for quasistatic shear flow. Olsson and Teitel^{3} found that the same correlation function does select a growing length scale, independent of L, in shear flow at finite rate using Stokes drag. However Tighe et al.^{57} showed that the form of C resembles quasistatic linear response when one uses balanced contact damping instead of Stokes drag. Hence there remain important open questions about correlations at finite driving rate and the role of the viscous force law. Here we fill a gap in the literature, namely linear response at finite rates.

In Fig. 7a we plot C for balanced contact damping at a single pressure, two system sizes, and three values of the frequency ω separated by twelve decades. There is a monotonic decay of the correlations, with little dependence on the frequency. The shape is also independent of the pressure (not shown). The data collapse when plotted as a function of δ/L. Hence two-point displacement correlations provide no evidence of a growing length scale near jamming; snapshots of the velocities display “swirls” with a characteristic radius of approximately one quarter of the box size.

Correlations for Stokes-like drag display a strikingly different shape, as shown in Fig. 7b. C possesses a minimum that shifts to larger distances with decreasing ω. For the lowest plotted frequencies, ω = 10^{−5.5} and 10^{−6.5}, the minimum is no longer clearly identifiable and the shape of C begins to resemble the form for balanced contact damping. One can define a correlation length l from the point where C crosses the x-axis, plotted in Fig. 7c. We find a length scale that grows with decreasing frequency, before reaching a plateau with a height of approximately L/4. Focusing on length scales below this plateau, we find empirically that a reasonable data collapse is achieved by plotting l/(−lnp)^{0.65}versus ω/p^{0.5}, implying that the length scale would diverge at the jamming point (p → 0 and ω → 0) in thermodynamically large systems. We note that log corrections are typical in systems at their upper critical dimension, which is indeed D = 2 for the jamming transition.^{31,58}

The takeaway is that the form of the correlation function at finite rate is strongly sensitive to the viscous force law. For balanced contact damping there is no evidence of a diverging length scale. For Stokes drag there is a growing correlation length that is cut off by the box size as the frequency is sent to zero.

In Fig. 8 the dynamic viscosities for both balanced contact damping and Stokes drag (η_{F} = 1) are plotted for a wide range of pressures and system sizes, both with pre-stress (open symbols) and without pre-stress (filled symbols). In all cases, we find that the data collapse to a master curve when η_{0}/N^{a} is plotted versus p/p* ∼ pN^{2}, implying η_{0} ∼ 1/p^{a/2}. For balanced contact damping, we find the best collapse when a = 1.0, consistent with the scaling η_{0} ∼ 1/p^{1/2} determined above. For Stokes drag we find better collapse for the somewhat higher value a = 1.09. For comparison, we also plot curves with slope −a/2. Provided that η_{0} is an intensive material property, as is typically the case, the master curves must approach this slope for large system sizes. This condition is met for the contact damping data, but for Stokes drag the collapsed data have a slightly shallower slope, particularly for the data with pre-stress. Using a lower value of a brings −a/2 closer to the observed slope, but the data collapse is somewhat worse. Given the small difference in these values and the scatter in our data, we consider it likely that a is in fact equal to 1 for Stokes drag. However, on the basis of present data we cannot exclude the possibility that a > 1 for Stokes drag, or that η_{0} has a weak system size dependence.

For both contact damping and Stokes drag, pre-stress plays a role in the onset of finite size effects. Whereas the data without pre-stress show a sharp crossover around , the crossover in the data with pre-stress is much more gradual. Even for pN^{2} > 10^{3}, a naïve power-law fit to η_{0}versus p would yield a slope that is too shallow. Therefore studying the results of simulations with and without pre-stress, side-by-side, can potentially improve the assessment of critical exponents near jamming at modest system sizes.

Nonlinear equations of motion cannot be written as a matrix equation in terms of and . Molecular dynamics simulations are an option,^{28} but beyond the scope of the present work. Instead, we turn to an approximation known as the method of equivalent damping. The central idea of the approximation is to replace the nonlinear force law with an “equivalent” linear force law with a frequency dependent effective damping coefficient kτ_{α},

^{eff}_{α} = −kτ_{α}Δ^{c}. | (23) |

We now apply the method of equivalent damping to a single degree of freedom system, namely an overdamped oscillator driven by a sinusoidal force with amplitude F_{0} and frequency ω. For the effective damping law of eqn (23), the resulting oscillations have an amplitude

(24) |

(25) |

(26) |

To extend the above insights to soft sphere packings, we make an additional but reasonable assumption that the typical induced force on each contact is proportional to the applied stress, F_{0} ∼ δσ. Under this assumption, the scaling ansatz (14)–(17) remains valid, provided that one takes τ* ∼ τ_{α}/p, instead of τ_{1}/p. Because τ_{α} is a function of frequency and the applied stress, the “bare” storage and loss moduli G′ and G′′ (as opposed to and ) inherit new dependences on ω and δσ. For systems near jamming and the physically relevant case α < 1, 1/τ* is always smaller than ω_{×} and hence τ_{α} ∼ (δσω)^{α−1} in the quasistatic and critical regimes. In the quasistatic regime one finds that the storage modulus G′ ≃ G_{0} is unchanged, while the loss modulus becomes

(27) |

(28) |

(29) |

We have also made predictions for the viscoelastic response in the presence of nonlinear drag laws. Within the context of the method of equivalent damping, we find that dynamic critical scaling survives; however the scaling of the bare storage and loss moduli now depends on the microscopic exponent α. This provides a novel way to infer properties of the dominant dissipative mechanism at the particle scale from the frequency dependence of G*. As the method of equivalent damping is an approximation, these predictions require further testing. As a basic check, we have verified eqn (27)–(29) by directly inserting the effective damping coefficient from eqn (26) in the linear equations of motion. Of course this does not constitute an independent test of the method of equivalent damping, which would require, e.g., molecular dynamics simulations of Durian's bubble model. We leave this as an important task for future work.

Our results suggest that, when performing numerical studies of jammed matter, one must take care to match the form of the viscous force law to the physics of whatever particular material one wishes to model – growing correlations do not wash out this detail. In particular, the linear contact damping law with β = 0 should be avoided, as it significantly alters the viscoelastic response and is difficult to justify on physical grounds, at least in the context of foams, emulsions, and soft colloidal particles.

- D. J. Durian, Phys. Rev. Lett., 1995, 75, 4780–4783 CrossRef CAS PubMed.
- S. Tewari, D. Schiemann, D. J. Durian, C. M. Knobler, S. A. Langer and A. J. Liu, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 60, 4385 CrossRef CAS.
- P. Olsson and S. Teitel, Phys. Rev. Lett., 2007, 99, 178001 CrossRef PubMed.
- V. Langlois, S. Hutzler and D. Weaire, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 021401 CrossRef PubMed.
- T. Hatano, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 050301 CrossRef PubMed.
- B. P. Tighe, E. Woldhuis, J. J. C. Remmers, W. van Saarloos and M. van Hecke, Phys. Rev. Lett., 2010, 105, 088303 CrossRef PubMed.
- J. Boschan, D. Vågberg, E. Somfai and B. P. Tighe, Soft Matter, 2016, 12, 5450–5460 RSC.
- J. Boschan, S. Vasudavan, P. E. Boukany, E. Somfai and B. P. Tighe, Soft Matter, 2017, 13, 6870–6876 RSC.
- C. S. O'Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 011306 CrossRef PubMed.
- F. Bolton and D. Weaire, Phys. Rev. Lett., 1990, 65, 3449–3451 CrossRef CAS PubMed.
- L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2005, 95, 098301 CrossRef PubMed.
- W. G. Ellenbroek, E. Somfai, M. van Hecke and W. van Saarloos, Phys. Rev. Lett., 2006, 97, 258001 CrossRef PubMed.
- B. P. Tighe, Phys. Rev. Lett., 2011, 107, 158303 CrossRef PubMed.
- E. Lerner, E. DeGiuli, G. Düring and M. Wyart, Soft Matter, 2014, 10, 5085–5092 RSC.
- K. Karimi and C. E. Maloney, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 022208 CrossRef PubMed.
- K. Baumgarten, D. Vågberg and B. P. Tighe, Phys. Rev. Lett., 2017, 118, 098001 CrossRef PubMed.
- K. Khakalo, K. Baumgarten, B. P. Tighe and A. Puisto, 2017, arXiv:1706.03932.
- H. A. Barnes and J. F. Hutton, An Introduction to Rheology, Elsevier, 1989 Search PubMed.
- A. J. Liu and S. R. Nagel, Nature, 1998, 396, 21–22 CrossRef CAS.
- S. Cohen-Addad, H. Hoballah and R. Höhler, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57, 6897–6901 CrossRef CAS.
- A. D. Gopal and D. J. Durian, Phys. Rev. Lett., 2003, 91, 188303 CrossRef CAS PubMed.
- J. M. Kropka and M. Celina, J. Chem. Phys., 2010, 133, 024904 CrossRef PubMed.
- J. J. Liétor-Santos, B. Sierra-Martn and A. Fernández-Nieves, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 060402 CrossRef PubMed.
- I. Srivastava and T. S. Fisher, Soft Matter, 2017, 13, 3411–3421 RSC.
- F. P. Bretherton, J. Fluid Mech., 1961, 10, 166 CrossRef.
- N. Denkov, S. Tcholakova, K. Golemanov, K. Ananthapadmanabhan and A. Lips, Phys. Rev. Lett., 2008, 100, 138301 CrossRef CAS PubMed.
- J. R. Seth, L. Mohan, C. Locatelli-Champagne, M. Cloitre and R. T. Bonnecaze, Nat. Mater., 2011, 10, 838–843 CrossRef CAS PubMed.
- B. Andreotti, J.-L. Barrat and C. Heussinger, Phys. Rev. Lett., 2012, 109, 105901 CrossRef PubMed.
- B. I. Halperin and P. C. Hohenberg, Phys. Rev., 1969, 177, 952–971 CrossRef CAS.
- D. J. Koeze, D. Vågberg, B. B. Tjoa and B. P. Tighe, EPL, 2016, 113, 54001 CrossRef.
- C. P. Goodrich, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2012, 109, 095704 CrossRef PubMed.
- M. van Hecke, J. Phys.: Condens. Matter, 2010, 22, 033101 CrossRef CAS PubMed.
- A. J. Liu and S. R. Nagel, Annu. Rev. Condens. Matter Phys., 2010, 1, 347–369 CrossRef.
- D. Morse and T. Witten, Europhys. Lett., 1993, 22, 549 CrossRef CAS.
- M.-D. Lacasse, G. S. Grest and D. Levine, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1996, 54, 5436 CrossRef CAS.
- B. P. Tighe, in Handbook of Granular Materials, ed. S. V. Franklin and M. D. Shattuck, CRC Press, 2016, ch. Wet Foams, Slippery Grains Search PubMed.
- S. Hutzler, R. P. Murtagh, D. Whyte, S. T. Tobin and D. Weaire, Soft Matter, 2014, 10, 7103–7108 RSC.
- R. Höhler and S. Cohen-Addad, Soft Matter, 2017, 13, 1371–1383 RSC.
- J. Winkelmann, F. Dunne, V. Langlois, M. Möbius, D. Weaire and S. Hutzler, Colloids Surf., A, 2017 DOI:10.1016/j.colsurfa.2017.03.058.
- C. E. Maloney and M. O. Robbins, J. Phys.: Condens. Matter, 2008, 20, 244128 CrossRef.
- D. Vagberg and B. P. Tighe, 2017, arXiv:1706.06378.
- P.-E. Peyneau and J.-N. Roux, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 011307 CrossRef PubMed.
- T. Hatano, J. Phys. Soc. Jpn., 2008, 77, 123002 CrossRef.
- M. Otsuki and H. Hayakawa, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 011308 CrossRef PubMed.
- C. Heussinger, P. Chaudhuri and J.-L. Barrat, Soft Matter, 2010, 6, 3050–3058 RSC.
- D. J. Evans and G. Morriss, Statistical Mechanics of Nonequilibrium Liquids, Cambridge University Press, 2008 Search PubMed.
- A. Ikeda, L. Berthier and P. Sollich, Phys. Rev. Lett., 2012, 109, 018301 CrossRef PubMed.
- B. P. Tighe, Phys. Rev. Lett., 2012, 109, 168303 CrossRef PubMed.
- E. Lerner, G. Düring and M. Wyart, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 4798–4803 CrossRef CAS PubMed.
- G. Katgert, M. E. Möbius and M. van Hecke, Phys. Rev. Lett., 2008, 101, 058301 CrossRef PubMed.
- W. G. Ellenbroek, M. van Hecke and W. van Saarloos, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 061307 CrossRef PubMed.
- M. Wyart, Ann. Phys., 2005, 30, 1 Search PubMed.
- S. Dagois-Bohy, B. P. Tighe, J. Simon, S. Henkes and M. van Hecke, Phys. Rev. Lett., 2012, 109, 095703 CrossRef PubMed.
- M. Wyart, H. Liang, A. Kabla and L. Mahadevan, Phys. Rev. Lett., 2008, 101, 215501 CrossRef CAS PubMed.
- B. A. DiDonna and T. C. Lubensky, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 066619 CrossRef CAS PubMed.
- C. E. Maloney and A. Lematre, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 016118 CrossRef PubMed.
- B. P. Tighe, J. H. Snoeijer, T. J. H. Vlugt and M. van Hecke, Soft Matter, 2010, 6, 2908 RSC.
- C. P. Goodrich, S. Dagois-Bohy, B. P. Tighe, M. van Hecke, A. J. Liu and S. R. Nagel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 022138 CrossRef PubMed.
- B. P. Tighe and T. J. H. Vlugt, J. Stat. Mech.: Theory Exp., 2011, P04002 Search PubMed.

This journal is © The Royal Society of Chemistry 2017 |