Granular matter like sand is composed of a large number of interacting grains and is, thus, expected to be amenable to a statistical physics treatment. Yet, the frictional properties of grains make the statistical physics of granular matter significantly different from the equilibrium statistical physics of atomic or molecular systems. We use three simple models to illustrate some of the key concepts of the statistical physics introduced by Edwards and co-workers more than 30 years ago to describe shaken granular piles: non-interacting frictional grains attached to a wall by a spring, a chain of frictional grains connected by springs, and a simplified mean-field model of a granular packing. We observe that a chain of frictional grains connected by springs exhibits a critical point at an infinite effective temperature (i.e., infinitely strong shaking) at odds with the zero-temperature critical point generically found in one-dimensional systems at thermal equilibrium in the presence of local interactions.

Whoever has played with sand on a beach may have been amazed by its mechanical properties. For instance, dry sand in an open hand flows almost like water. Yet, the beach forms an essentially solid surface on which you can walk, leaving some footprints that are only partially erased by the spontaneous relaxation dynamics of dry sand. One may, thus, wonder how materials like sand that are formed of grains of macroscopic size (and thus, called “granular materials”) can exhibit mechanical properties that are intermediate between solids and liquids.1 The reason lies in the frictional properties of grains. Two grains in contact experience dry friction, which means that they can support tangential forces at contact without sliding.2 

The fact that granular materials are made of a large number of grains naturally calls for a statistical description, which is, however, expected to differ from the equilibrium statistical physics formalism that describes materials made of atoms or molecules. While the latter has a conservative dynamics, macroscopic particles like grains have a dissipative dynamics, which deeply modifies their large-scale statistical properties.3 In more abstract terms, the time reversal symmetry is broken in granular materials.4 This important difference has been taken into account through a minimal generalization of the equilibrium statistical physics framework, as proposed by Edwards and co-workers in the late 1980s.5,6 The goal of this paper is to illustrate the Edwards theory of granular matter with simple and pedagogical examples.

The Edwards theory introduces a thermodynamic parameter called compactivity that controls the statistics of the volume of the granular pile,5,6 just like temperature controls the energy statistics in equilibrium systems. Some generalizations of Edwards's original proposal have also introduced an effective temperature to characterize the statistics of shaken granular packings.7,8 This effective temperature is many orders of magnitude larger than the thermodynamic temperature of the grains. Thermal fluctuations associated with the thermodynamic temperature are, thus, negligible. Indeed, they are completely unable to lift a grain over a height of the order of its diameter.2 

The article is organized as follows. Section II briefly introduces the general framework of the Edwards statistical mechanics for systems with dry friction and puts it in the general perspective of the description of large assemblies of interacting particles, whose equilibrium statistics is given by the well-known Boltzmann–Gibbs distribution. A simplified version of the Edwards statistics is then introduced for pedagogical purposes. Section III describes a simple model of non-interacting frictional particles attached to a spring, which generalizes the harmonic oscillator model of statistical mechanics. Then, Sec. IV discusses a more complicated model of frictional particles connected by springs. The presence of interactions between particles generates strong correlations in a high temperature regime, at variance with usual equilibrium systems where correlations appear at low temperature. Finally, Sec. V presents a simplified model of a granular packing, which is treated using mean-field approximations.

We start here, for completeness, by recalling the basic statistical description of systems of interacting particles at equilibrium, called the Boltzmann–Gibbs statistics. This will then allow us to put the Edwards statistics into context and to emphasize its similarities and differences with the Boltzmann–Gibbs statistics.

Consider N identical particles of mass m with positions r i and velocities v i ( i = 1 , , N). The microscopic configuration of the system is defined by the list C = ( r 1 , , r N , v 1 , , v N ) of particle positions and velocities. Each particle is subjected to an external potential U ext ( r i ). Pairs of particles (i, j) interact via an interaction potential U int ( r i r j ). The total potential energy of the system is, thus,
U ( r 1 , , r N ) = i = 1 N U ext ( r i ) + i < j U int ( r i r j ) .
(1)
The total energy of the system in a configuration C is
E ( r 1 , , r N , v 1 , , v N ) = i = 1 N 1 2 m v i 2 + U ( r 1 , , r N ) .
(2)
If the system is at equilibrium with a thermal reservoir at temperature T, the probability distribution of a configuration C is given by the following equation:9–11,
P eq ( r 1 , , r N , v 1 , , v N ) = 1 Z exp [ β E ( r 1 , , r N , v 1 , , v N ) ] ,
(3)
where β = ( k B T ) 1 is the inverse temperature (with kB the Boltzmann constant), and
Z = d r 1 d r N d v 1 d v N exp [ β E ( r 1 , , r N , v 1 , , v N ) ]
(4)
is a normalization constant called the partition function ( d r i and d v i are multidimensional integration elements). Equation (3) is called the Boltzmann–Gibbs (or canonical) equilibrium distribution. Average observables can be evaluated from Eq. (3) and can be formally written as derivatives of ln Z. For instance, the average energy E is obtained as
E = ln Z β .
(5)
A particularly simple case consists of considering non-interacting particles ( U int = 0), because the partition function then factorizes as a product of one-particle partition functions. A typical example is the case of non-interacting particles in a potential U ext ( r ) = 1 2 k r 2, a situation referred to as non-interacting harmonic oscillators (the parameter k is called the stiffness of the harmonic potential). In this case, Z = Z 1 N, where Z1 is the one-particle partition function. Considering for the sake of simplicity a one-dimensional geometry, one has
Z 1 = d x d v exp [ β ( 1 2 m v 2 + 1 2 k x 2 ) ] ,
(6)
yielding Z 1 = 2 π / ( β m k ). It follows from Eq. (5) that
E = N k B T .
(7)
This is a particular instance of the equipartition theorem, according to which each degree of freedom associated with a quadratic energy (here, the 2N positions or velocities) carries an average energy 1 2 k B T. As a result, the specific heat C V = d E / d T = N k B is constant, a result called the Dulong and Petit law in condensed matter physics.12 Note that a constant specific heat is obtained only in the framework of classical mechanics. When considering quantum mechanical effects, the discreteness of energy states leads to a drop of specific heat at low temperature.12 

For interacting particles, the evaluation of the partition function is in general much more involved. However, for particles placed on a lattice and interacting through a harmonic potential, U int ( r ) = 1 2 k r 2, one can define non-local degrees of freedom called phonons that do not interact with each other, which makes the determination of thermodynamic properties simpler.13 

Below, we turn to the statistical description of granular piles and discuss how the Edwards description generalizes the equilibrium approach to this non-equilibrium situation.

A statistical description of a granular pile is meaningful if the pile is able to explore many different configurations. Experimentally, it is convenient to inject energy through a “tapping” protocol, by which the pile is repeatedly shaken and then let relax (after switching off the shaking mechanism) to a mechanically stable configuration, also called a blocked configuration. A blocked configuration is such that the sum of all forces acting on any given grain is equal to zero. The tapping protocol allows for many different blocked configurations to be explored. The pile can then be described by a statistics of blocked configurations, which allows for the prediction of average values of observables like the height of the pile or the force exerted by the grains on the container.

A blocked configuration of a granular pile made of N grains is defined by the list ( r 1 , , r N ) of all grain positions (all velocities are then null). The Edwards approach first postulates that only blocked configurations have a non-zero probability. This is justified by the tapping protocol, which samples only blocked configurations. Then, taking inspiration from equilibrium statistical mechanics, the idea of Edwards and co-workers5,6 is to assume that the statistics of blocked configurations takes the simplest possible form, given the macroscopic constraints. In the case of hard grains, one may assume that the most relevant macroscopic observable is the volume. Edwards and co-workers postulated that the granular pile should be described by the most likely probability distribution of blocked configurations with a given average value of the total volume, namely5,6,14 (see Ref. 4 for a review),
P ( r 1 , , r N ) = 1 Z V exp [ V ( r 1 , , r N ) / X ] F ( r 1 , , r N ) ,
(8)
where V ( r 1 , , r N ) is the volume occupied by the grains in configuration ( r 1 , , r N ), and ZV is the normalization factor
Z V = d r 1 d r N exp [ V ( r 1 , , r N ) / X ] F ( r 1 , , r N ) .
(9)
The parameter X is called compactivity. The indicator function F ( r 1 , , r N ) has been included to select blocked configurations among all possible configurations: F ( r 1 , , r N ) = 1 for a blocked configuration ( r 1 , , r N ), whereas F ( r 1 , , r N ) = 0 otherwise. Although this general definition is simple, the explicit form of F ( r 1 , , r N ) may be quite involved in practice.15–19 

Both experimental20–23 and numerical7,24–27 tests of the Edwards probability distribution have been performed. At a qualitative level, these tests confirm that the Edwards distribution captures most of the phenomenology of granular matter, although deviations from the predictions of the Edwards theory are also observed. However, a quantitative assessment of the predictions of the Edwards theory is made difficult by the complexity of the function F ( r 1 , , r N ),15–19 which can often be evaluated only through rather strong approximations.15,28 It may then be hard to disentangle the discrepancies resulting from the Edwards prescription and those resulting from the approximations made in the evaluation of F ( r 1 , , r N ), when comparing predictions with experimental or numerical results. The consistency of the Edwards framework has been validated experimentally: the compactivity X has been measured, based, e.g., on a measurement of the variance of volume fluctuations, which can be related to the derivative of the average volume with respect to compactivity.26,29–31 However, it has also been shown experimentally that the measured compactivity fails to equalize between a subsystem of grains and a surrounding “bath,” when grains in the subsystem have different frictional properties compared to grains in the bath.32 This lack of equilibrium points to a weakness of the Edwards theory if the latter is conceived as a thermodynamic-like theory, where compactivity is understood as a temperature-like parameter, as suggested by the form of Eq. (8), which should equalize between a system and a bath. The experiment performed in Ref. 32 also shows that a parameter similar to the compactivity but conjugated to the mechanical stress, called angoricity, does equalize between subsystem and bath. This result suggests that an Edwards-like theory based on the mechanical stress33–37 rather than on the volume may be more successful from a thermodynamic perspective—see Ref. 4 for a review.

The complexity of the indicator function F ( r 1 , , r N ) is not only a difficulty for researchers working on granular matter but it also constitutes a hurdle in order to teach the Edwards statistical approach in university courses, because simple pedagogical examples illustrating the approach are lacking. The goal of the present paper is to propose such elementary models that may be suitable for teaching purposes.

The indicator function F ( r 1 , , r N ) actually encodes both the complex contact network in a grain packing and the effect of dry friction. In the absence of dry friction, the balance of contact forces characterizing blocked states boils down to the local minimization of the potential energy describing repulsive forces between particles. In the presence of dry friction, blocked configurations are not obtained by such a minimization, but they rather form a continuum of configurations whose extent is directly related to the dry friction coefficient.

Aiming at simple and pedagogical models illustrating some basics of the Edwards statistics, we choose here to maximally simplify the contact network geometry and to focus on dry friction as the key ingredient to generate statistical properties that differ from the Boltzmann–Gibbs ones. Alternatively, one may also treat the contact network using a mean-field approximation, as done, for instance, in the mean-field description of the Ising model.9 We give a flavor of such a mean-field approach in Sec. V.

As a further simplification, we choose in a large part of this paper (except in Sec. V) to work with a formulation of the Edwards distribution in terms of energy E ( r 1 , , r N ) rather than volume V ( r 1 , , r N ). The modified Edwards distribution reads
P ( r 1 , , r N ) = 1 Z E exp [ E ( r 1 , , r N ) / T eff ] F ( r 1 , , r N ) ,
(10)
where ZE is a normalization factor (or effective partition function). The energy E ( r 1 , , r N ) corresponds to the potential energy associated with gravity and, for soft grains, with repulsive elastic forces between grains. The parameter T eff plays a role similar to the thermodynamic temperature T in equilibrium systems or more precisely to k B T. For these reasons, T eff is called effective temperature, although it has the dimension of an energy. As shown below, T eff can in some regimes be determined from the equipartition of energy between grains.

Working with an observable different from the volume does not actually rule out the Edwards approach. While the original approach was formulated in terms of the volume, further developments included the energy,7,8,14 and modern formulations of the Edwards statistics also take into account the stress tensor as a key observable, as mentioned in Sec. II B. In the pedagogical perspective of the present paper, working with energy also has further advantages. First, having in mind the simplest case of non-interacting particles, it is conceptually easy to define a potential energy associated with a single particle: one may simply consider the gravitational potential energy of the particle or attach the particle to a spring and consider the elastic potential energy of the spring. The corresponding physical situation is clear. In contrast, it is much more difficult to define the volume associated with a single particle without considering it part of a dense packing of particles. Hence, the volume of non-interacting particles is ill-defined. Second, the energy formulation of the Edwards distribution allows for a more straightforward comparison of the phenomenology resulting from the Edwards statistics and the one obtained from the equilibrium Boltzmann–Gibbs statistics. As a simple and concrete example, one may compare the Edwards and Boltzmann–Gibbs statistics for non-interacting harmonic oscillators. This is the topic of Sec. III.

Equilibrium statistical physics lectures usually start by describing the simplest possible examples, namely, systems consisting of a large number of non-interacting particles.9–11 These include the ideal gas as well as assemblies of non-interacting harmonic oscillators that may be interpreted for instance as the Einstein model for crystalline solids.12 

To illustrate the statistical physics of frictional particles proposed by Edwards and co-workers, we now describe a model of non-interacting harmonic oscillators with dry friction (Fig. 1). Let us consider an ensemble of N identical masses m moving on a horizontal substrate (e.g., a table) with dry friction coefficient μ, and attached to a spring of stiffness k. Masses are constrained to move parallel to the same axis x. Blocked configurations are sampled using a driving protocol that periodically injects energy into the system. Each period is decomposed into a driving phase of duration τ, during which a strong external force is applied to each block, and a relaxation phase during which the system relaxes to a blocked configuration. The position xi of mass i evolves according to
m d 2 x i d t 2 = k x i μ m g sign ( d x i d t ) + f i ( t ) .
(11)
Fig. 1.

Sketch of the frictional harmonic oscillator model (top view). Each mass i, at position xi, is moving along a horizontal direction parallel to the axis x. The driving force f i ( t ) injects energy to the mass i during the tapping phase.

Fig. 1.

Sketch of the frictional harmonic oscillator model (top view). Each mass i, at position xi, is moving along a horizontal direction parallel to the axis x. The driving force f i ( t ) injects energy to the mass i during the tapping phase.

Close modal

The first term on the rhs of Eq. (11) corresponds to the force exerted by the spring, assuming that xi = 0 corresponds to its rest position. The second term corresponds to the dynamic dry friction force, where mg is the mass weight (with g being the gravitational acceleration). Note that the absolute value of the force is independent of the speed, at variance with viscous friction. Finally, the last term f i ( t ) corresponds to the driving force used to inject energy during the tapping protocol [ f i ( t ) = 0 during the relaxation phase leading the system to a blocked state]. It could model, for instance, a horizontal vibration of the table on which the masses are placed. Note that fi plays no role in the energy E appearing in Eq. (10), since fi = 0 in blocked states.

When at rest, the mass is subjected to a static dry friction force, which exactly compensates the other horizontal forces as long as these forces do not overcome in the absolute value a threshold value equal to μmg. Note that in realistic granular systems, the static friction coefficient defining the threshold force slightly differs from the dynamic one, but we neglect this slight difference here for the sake of simplicity. This is a common assumption when modeling granular systems. When the threshold is exceeded, the mass starts to move. In the absence of driving force ( fi = 0), the mass starts to move if k | x i | > μ m g. Conversely, blocked configurations correspond to static configurations that verify the condition k | x i | < μ m g. In other words, blocked configurations satisfy | x i | < a, with a characteristic length scale
a = μ m g k .
(12)
Note that the length scale a is proportional to the friction coefficient μ, so that a is nonzero for frictional systems only.
We now turn to the determination of the stationary probability distribution P ( x 1 , , x N ) associated with the set ( x 1 , , x N ) of positions of the N masses, using the Edwards prescription [Eq. (10)] and replacing the position vector r i by the one-dimensional position variable xi. The total energy of a configuration is given by
E ( x 1 , , x N ) = i = 1 N 1 2 k x i 2 .
(13)
The indicator function F ( x 1 , , x N ) appearing in Eq. (10) can be formally written as
F ( x 1 , , x N ) = i = 1 N Θ ( a | x i | ) ,
(14)
where Θ ( x ) is the Heaviside function, Θ ( x ) = 1 for x 0 and Θ ( x ) = 0 for x < 0. Equation (14) simply means that a configuration ( x 1 , , x N ) is a blocked configuration if all xi satisfy | x i | < a. It follows that the Edwards distribution Eq. (10) factorizes as
P ( x 1 , , x N ) = i = 1 N p ( x i ) ,
(15)
with a one-body distribution p(x) given by
p ( x ) = 1 Z 1 exp [ β k x 2 / 2 ] if | x | < a ,
(16)
and p(x) = 0 otherwise, with the notation β = 1 / T eff. The normalization factor Z1 defined as
Z 1 = a a d x exp [ β k x 2 / 2 ]
(17)
plays the role of a one-body partition function.
The average energy of the N harmonic oscillators is then obtained as
E = N Z 1 a a d x k 2 x 2 exp [ β k x 2 / 2 ] .
(18)
At this stage, it is convenient to introduce the characteristic effective temperature
T * = 1 3 k a 2
(19)
(the reason for introducing the 1/3 factor is discussed below). Note that the temperature scale T * originates from the presence of static friction, as can be seen from the expression of a given in Eq. (12).
Thanks to the change of variable y = β k x, the average energy E can be rewritten in the form (see Fig. 2),
E = 1 2 N T eff f h ( T eff T * ) ,
(20)
having defined the auxiliary function f h ( u ) as
f h ( u ) = 0 3 / u d y y 2 exp [ y 2 / 2 ] 0 3 / u d y exp [ y 2 / 2 ] .
(21)
From this integral expression of the function f h ( u ), its asymptotic behaviors can be determined. One finds f h ( u ) 1 for u 0 as well as f h ( u ) 1 / u when u + .
Fig. 2.

Dimensionless energy density E / N T * in the frictional harmonic oscillator model, as a function of the reduced temperature T eff / T *, showing the energy saturation at high temperature to the value E / N T * = 1 / 2 (dotted line). Inset: specific heat C V / N vs T eff / T *, showing a drop from the classical value 1/2 (equivalent to the Dulong and Petit law) at low temperature to zero at high temperature.

Fig. 2.

Dimensionless energy density E / N T * in the frictional harmonic oscillator model, as a function of the reduced temperature T eff / T *, showing the energy saturation at high temperature to the value E / N T * = 1 / 2 (dotted line). Inset: specific heat C V / N vs T eff / T *, showing a drop from the classical value 1/2 (equivalent to the Dulong and Petit law) at low temperature to zero at high temperature.

Close modal
The average energy E then satisfies the generalized equipartition relations in the two asymptotic regimes T eff T * and T eff T *,
E 1 2 N T eff for T eff T * ,
(22)
E 1 2 N T * for T eff T * .
(23)
It follows that for low effective temperatures, a form analogous to the equilibrium equipartition theorem E = 1 2 N k B T for harmonic oscillators is obtained: The average energy is proportional to the effective temperature. In contrast, for a high effective temperature, the average energy reaches a maximum value, set by the dry friction coefficient—see Eqs. (19) and (12). Note that the 1/3 factor in the definition (19) of T * has been specifically introduced so that Eqs. (22) and (23) bear the same form.

The breaking of the energy equipartition relation (22) can be evidenced more clearly by considering the specific heat C V = d E / d T eff. As recalled in Sec. II A, energy equipartition corresponds to a constant specific heat C V = 1 2 N (the Boltzmann constant kB has been absorbed into the definition of the effective temperature T eff). In the present frictional harmonic oscillator model, a constant specific heat C V = 1 2 N is observed at low effective temperature ( T eff T *). Conversely, in the high temperature regime, the specific heat decreases and eventually goes to zero. Using the small-u expansion f h ( u ) u 2 5 u 2, we obtain C V 1 5 ( T * / T eff ) 2 for T eff T * (see inset of Fig. 2). Thus, the presence of an upper energy bound leads to a departure from a constant specific heat at high effective temperature, contrary to quantum harmonic oscillators at thermal equilibrium, for which departure from a constant specific heat occurs at low temperature due to the discreteness of energy states.12 

We now go beyond the above non-interacting case and turn to a second example of frictional models, to explore further the interesting phenomenology emerging from dry friction. We consider a one-dimensional chain of N + 1 frictional blocks of mass m on a substrate. Each block i = 0 , , N, located at position x i ( t ), experiences dry friction from the substrate. Blocks are connected by springs of stiffness k and rest length l0, as illustrated in Fig. 3. Springs, thus, induce interactions between blocks. This model has been initially introduced in the context of earthquake modeling38,39 and further studied in Refs. 40 and 41 as a generic model of frictional particles with elastic interactions. To avoid dealing with the no-crossing condition between blocks, we assume the rest length l0 to be significantly larger than the typical value of spring extensions. Similarly to the frictional harmonic oscillator model, a block at rest starts moving if the total force exerted on it by neighboring springs overcomes a threshold force equal to the weight mg times the static friction coefficient μ. We can, thus, define a blocked configuration as a configuration ( x 0 , , x N ), such that resulting spring forces on each block do not exceed μmg.

Fig. 3.

Sketch of the frictional spring–block model. Masses at positions xi are connected by springs of extension ξ i = x i x i 1 l 0, with l0 being the rest length of the springs. External forces f i ( t ) are applied during the driving phase to inject energy into the system. These forces are null during the relaxation phase.

Fig. 3.

Sketch of the frictional spring–block model. Masses at positions xi are connected by springs of extension ξ i = x i x i 1 l 0, with l0 being the rest length of the springs. External forces f i ( t ) are applied during the driving phase to inject energy into the system. These forces are null during the relaxation phase.

Close modal
The driving and relaxation phases of the tapping protocol are described by the following dynamics:
m d 2 x i d t 2 = μ m g sign ( d x i d t ) + k ( x i + 1 + x i 1 2 x i ) + f i ( t ) ,
(24)
where f i ( t ) is the external force applied during the driving phase [ f i ( t ) = 0 during the relaxation phase]. One may define x 1 x 0 l 0 and x N + 1 x N + l 0 so that Eq. (24) formally holds also for i = 0 and i = N. In principle, a configuration of the spring–block model is defined by the list of positions ( x 0 , x 1 , , x N ). However, two configurations ( x 0 , x 1 , , x N ) and ( x 0 + b , x 1 + b , , x N + b ) differing by a global translation can be considered as equivalent, because the dynamics given in Eq. (24) are invariant under such a translation. We, thus, rather characterize a configuration by the list of spring elongations ( ξ 1 , ξ 2 , , ξ N ), with ξ i = x i x i 1 l 0.
As mentioned above, a blocked configuration is defined by the set of conditions
k | ξ i + 1 ξ i | < μ mg ( i = 1 , , N 1 ) .
(25)
This condition can be reformulated by introducing the same length scale a = μ m g / k as in Eq. (12), leading to
| ξ i + 1 ξ i | < a ( i = 1 , , N 1 ) .
(26)
According to the Edwards postulate, the distribution P ( ξ 1 , ξ 2 , , ξ N ) of microscopic configurations is given by Eq. (10), where the energy E ( ξ 1 , ξ 2 , , ξ N ) is the total elastic energy
E ( ξ 1 , , ξ N ) = i = 1 N 1 2 k ξ i 2 ,
(27)
and the indicator function F ( ξ 1 , ξ 2 , , ξ N ) can be written as
F ( ξ 1 , , ξ N ) = i = 1 N 1 Θ ( a | ξ i + 1 ξ i | ) .
(28)
Equation (28) is a formal and compact way to express the list of conditions given in Eq. (26).
The Edwards probability distribution P ( ξ 1 , ξ 2 , , ξ N ), thus, reads
P ( ξ 1 , , ξ N ) = 1 Z E exp [ i = 1 N k ξ i 2 2 T eff ] i = 1 N Θ ( a | ξ i + 1 ξ i | ) ,
(29)
with ZE being a normalization factor defined as
Z E = d ξ 1 d ξ N exp [ i = 1 N k ξ i 2 2 T eff ] i = 1 N Θ ( a | ξ i + 1 ξ i | ) ,
(30)
playing the same role as a partition function in equilibrium systems. As before, the effective temperature T eff is at this stage an auxiliary parameter that cannot be measured directly in a numerical simulation, and one needs to find a relation connecting T eff to the energy density ε = E / N, which is a measurable observable.
Due to the formal analogy between the Edwards distribution (10) and the usual equilibrium canonical distribution, similar relations hold between energy and temperature. With the notation β = T eff 1, the average energy is obtained as
E = ln Z E β ,
(31)
as can be checked by a direct calculation. The effective partition function ZE can be evaluated using a transfer operator technique. This method generalizes the transfer matrix technique classically used to solve for instance the one-dimensional Ising model,13 replacing the transfer matrix by an infinite dimensional operator that can be handled using numerical methods. The interested reader is referred to Ref. 41 for details. An alternative method, which we now describe, uses an approximation that leads to analytically tractable calculations. The basic idea is to replace the “door” function Θ ( a | Δ ξ | ) appearing in Eq. (28) by a Gaussian function exp [ ( Δ ξ ) 2 / 2 a 2 ] with a width a.
Under this approximation, the partition function ZE defined in Eq. (30) can be rewritten as a Gaussian multidimensional integral
Z E = d ξ 1 d ξ N exp [ β E ̃ ( ξ 1 , , ξ N ) ] ,
(32)
where the quantity E ̃ ( ξ 1 , , ξ N ) plays the role of an effective energy, which is quadratic in the degrees of freedom ξi,
E ̃ ( ξ 1 , , ξ N ) = i = 1 N [ 1 2 k ξ i 2 + 1 2 a 2 β ( ξ i + 1 ξ i ) 2 ] .
(33)
Note that we have used here for the sake of simplicity periodic boundary conditions ξ N + r ξ r, but this choice has no consequence on the evaluation of global observables like the energy, for N + . Thanks to this reformulation, the indicator function F ( ξ 1 , , ξ N ) has been reabsorbed into the quadratic effective energy E ̃ ( ξ 1 , , ξ N ), which makes the mathematical treatment of the problem much easier. In more formal terms,
Z E = d ξ 1 d ξ N exp [ 1 2 ξ · A ξ ] ,
(34)
having introduced the vector ξ = ( ξ 1 , , ξ N ) and the symmetric matrix A defined as
( A ξ ) j = β k ξ j + 1 a 2 ( 2 ξ j ξ j + 1 ξ j 1 ) .
(35)
The general formula for multidimensional Gaussian integrals like the one of Eq. (34) reads13 
Z E = ( 2 π ) N / 2 ( det A ) 1 / 2 ,
(36)
where det A is the determinant of the matrix A, which may be obtained as the product of its eigenvalues. Eigenvectors here correspond to Fourier modes ξ j ( q ) = e i q j ( i 2 = 1), and the corresponding eigenvalue reads
λ q = β k + 2 a 2 ( 1 cos q ) ,
(37)
where q = 2 π n / N ( n = 0 , , N 1). One then finds,
ln det A = n = 0 N 1 ln [ β k + 2 a 2 ( 1 cos 2 π n N ) ] ,
(38)
eventually yielding in the large N limit, using Eq. (36) and the increment Δ q = 2 π / N,
ln Z E = N 2 ln ( 2 π ) N 4 π 0 2 π d q ln [ β k + 2 a 2 ( 1 cos q ) ] .
(39)
The average energy E then follows by differentiating ln Z E with respect to β, according to Eq. (31),
E = 1 2 N T eff f s ( T eff T * )
(40)
with
T * = k a 2
(41)
[note the slightly different definition of T * with respect to Eq. (12)], and a function f s ( u ) defined as
f s ( u ) = 1 π 0 π d q 1 + 2 u ( 1 cos q ) = 1 1 + 4 u ,
(42)
where the second equality uses an integration formula from Ref. 42. The form Eq. (40) of the average energy E is plotted in Fig. 4. We, thus, obtain the two asymptotic regimes:
E 1 2 N T eff for T eff T * ,
(43)
E 1 4 N T * T eff for T eff T * .
(44)
Fig. 4.

Dimensionless energy density E / N T * in the spring–block model, as a function of the reduced temperature T eff / T *. The low- and high-temperature regimes are emphasized. Inset: specific heat C V / N vs T eff / T *, showing a slow decrease from the classical value 1/2 at low temperature to zero at high temperature.

Fig. 4.

Dimensionless energy density E / N T * in the spring–block model, as a function of the reduced temperature T eff / T *. The low- and high-temperature regimes are emphasized. Inset: specific heat C V / N vs T eff / T *, showing a slow decrease from the classical value 1/2 at low temperature to zero at high temperature.

Close modal

Hence, one finds a low energy regime where the energy density ε = E / N is proportional to the effective temperature T eff similarly to the equipartition relation valid at equilibrium, and a high energy regime where ε is proportional to T eff. At odds with the frictional oscillator model, the energy density ε does not saturate to a maximal value when increasing the effective temperature T eff, because dry friction does not impose an upper bound on spring extension [strong forces applied on both sides of a mass may almost exactly compensate and satisfy the dry friction stability criterion Eq. (26)]. Despite the lack of an upper bound on the energy, the specific heat C V = d E / d T eff goes to zero at high effective temperature, as C V 1 / T eff (see the inset of Fig. 4). This slow decay of CV, resulting from the absence of upper bound for the energy, is to be compared with the faster decay C V 1 / T eff 2 for the frictional harmonic oscillator model of Sec. III, for which the energy has an upper bound.

In addition to calculating effective thermodynamic properties like the average energy, it is also of interest to determine the extent of spatial correlations in the system as a function of the effective temperature. In one-dimensional equilibrium systems with local interactions, like the one-dimensional Ising model for instance,13 the correlation length diverges when the temperature goes to zero, which corresponds to a zero-temperature critical point. Given the formal similarity between the Edwards probability distribution (10) and the Boltzmann–Gibbs one, one might expect the spring–block model to also have a zero-temperature critical point. We will see below that, quite unexpectedly, the spring–block model has a critical point at infinite effective temperature, due to the constraints imposed by frictional properties.

Let us define the correlation function Cr of the elongation of springs separated by a distance r,
C r = 1 N j = 1 N ξ j ξ j + r ,
(45)
again using periodic boundary conditions. Note that the distance r is measured in numbers of springs, rather than as a geometric length. To proceed further, it is convenient to introduce the discrete Fourier transform ξ ̂ q as
ξ ̂ q = 1 N j = 0 N 1 exp [ i q j ] ξ j ,
(46)
with q = 2 π n / N ( n = 0 , , N 1). Similarly, one defines the discrete Fourier transform C ̂ q of the correlation function Cr,
C ̂ q = 1 N r = 0 N 1 exp [ i q r ] C r .
(47)
Using Eqs. (45) and (46), one finds C ̂ q = | ξ q | 2 . Since Fourier modes diagonalize the matrix A, standard properties of Gaussian multidimensional integrals13 lead to | ξ q | 2 = 2 / λ q, where λq is the eigenvalue defined in Eq. (37). One, thus, finds for small | q | (assuming N to be large),
C ̂ q = 2 ( β k + q 2 2 a 2 ) 1 .
(48)
In the following, we approximate at large N the discrete Fourier transform by a continuous Fourier transform. Noticing that the continuous Fourier transform of an exponential correlation function
C ( r ) = C 0 exp [ | r | / ]
(49)
reads
C ̂ ( q ) = 2 C 0 1 + ( q ) 2 ,
(50)
we identify from Eq. (48) the expression = T eff / ( 2 T * ) of the correlation length, using the definition (41) of T *. A similar scaling has also been found in the presence of viscous damping, although the Edwards statistics has to be generalized.43 

Note that the correlation length is dimensionless because it is measured in numbers of springs rather than in geometric length. To approximately convert it to a geometric length, one may simply multiply it by the rest length l0.

The continuous approximation of the Fourier transform used to obtain Eq. (50) is meaningful only when the correlation length is large, 1 (so that C(r) decays smoothly, almost like a continuous function), which implies T eff T *. This regime is achieved through a strong energy injection rate, i.e., large external forces f i ( t ) during the driving phase. Given that, in this regime, the energy density also scales as ε T eff according to Eq. (44), we end up with the simple scaling relation ε. One, thus, concludes that the model exhibits a critical point at infinite effective temperature, or infinite energy density. This unexpected property may be interpreted as follows. In the high energy regime, spring extensions ξi are typically much larger than the length a characterizing blocked configurations. Since the probability distribution is restricted to blocked configurations, condition (26) is satisfied, with both ξi and ξ i + 1 much larger than a. This implies that ξi and ξ i + 1 are nearly equal, and are, thus, strongly correlated. The decorrelation of the spring elongations ξi and ξ i + 1 occurs only over distances r .

As a last example of models of frictional granular matter, we consider a simplified model of granular packing, as schematically illustrated in Fig. 5(a). In this model, the Edwards distribution is naturally formulated in terms of the volume, as given in Eq. (8). As already mentioned in Sec. II C, a mathematical advantage of the energy as a macroscopic observable is that the total energy can be expressed as a sum of single-particle and particle-pair contributions. This has been exploited in Secs. III and IV. In contrast, a decomposition of the total volume into contributions associated with each grain leads to multi-particle interactions that cannot be dealt with without resorting to rather strong approximations. We briefly describe below such a simple model, following Ref. 28. Given a configuration ( r 1 , , r N ) of the N grains, the volume V ( r 1 , , r N ) can be expressed as a sum of individual volumes attached to each grain using the Voronoi tessellation22,28
V ( r 1 , , r N ) = i = 1 N v i ( r 1 , , r N ) ,
(51)
where v i ( r 1 , , r N ) is the volume of the Voronoi cell surrounding grain i. The construction of a Voronoi cell around a given grain is illustrated in Fig. 5(b). Quite importantly, the individual volume v i ( r 1 , , r N ) depends not only on the position of grain i but also on the positions of all the nearest neighbor grains. A natural question in this context may be to determine the distribution P 1 ( v ) of the individual volume vi (the distribution P 1 ( v ) is assumed to be the same for all grains). Starting from the N-particle volume distribution given in Eq. (8) according to the Edwards postulate, one can re-express it using Eq. (51). One should in principle integrate the distribution Eq. (8) over all positions ( r 1 , , r N ) compatible with a given value v i ( r 1 , , r N ) = v of the individual volume of a fixed grain i. Unfortunately, under this form, the calculation of the distribution P 1 ( v ) is essentially intractable. One, thus, needs to resort to some approximations. Following again Ref. 28, we use a rough mean-field approximation in order to get a simple result that hopefully still captures the essence of the physical phenomenon. First of all, we note that the determination of the distribution P 1 ( v ) would be much simpler if we could interpret the N-particle distribution Eq. (8) as a joint distribution P ( v 1 , , v N ) of the individual volumes vi rather than as a joint distribution of positions. For this approximation to be meaningful, one would need to be able to express the indicator function F ( r 1 , , r N ) in terms of the local volumes only. This is not possible for the exact indicator function, but one may try to approximate it in the following way. We first note that the individual volume vi cannot be smaller than a minimal value v min corresponding to the densest local packing, since grains are assumed to be hard particles that cannot be deformed. For two-dimensional disks of radius r0, the densest local packing is the hexagonal one, corresponding to a minimal volume v min 2 D = 12 r 0 2. In three dimensions, the most compact packing of hard spheres of radius r0 is the face centered cubic packing, corresponding to v min 3 D = 32 r 0 3.28 One may also consider non-circular or non-spherical grains, but in this case the determination of the densest packing is much more challenging. Based on physical intuition, the accessible volume also needs to have a maximal value, since if the volume of the Voronoi cell around a given grain becomes large, the grain would no longer be in contact with neighboring grains and stability of the packing would be lost. The maximum value of the local volume is expected to depend on the dry friction coefficient between grains, as supporting tangential forces at contacts enhances the packing stability. Accordingly, the maximal local volume v max should be an increasing function of the dry friction coefficient μ. An approximate expression of v max has been obtained in Ref. 28 using a simplified geometrical reasoning, by piling up one grain on top of two others over a plane. This is by no means a realistic evaluation of v max, but it has the advantage of providing an explicit expression of the latter, and to highlight the role of the friction coefficient μ. For two-dimensional circular grains of radius r0, one finds28 
v max = 2 π r 0 2 μ 1 μ 2 [ cos 1 ( 1 μ 2 1 + μ 2 ) ] 1 ,
(52)
an approximate expression valid in the realistic range of values of the dry friction coefficient, 0.3 μ 0.5. It can be checked that v max is an increasing function of μ, as expected. Note that in this two-dimensional case, the volume v max is actually an area.
Fig. 5.

(a) Sketch of two-dimensional granular packing in a container. (b) Illustration of the Voronoi cell surrounding a grain. The boundary of the Voronoi cell is the set of points that are at equal distance from the center of the focus grain and from the centers of the nearest neighbor grains.

Fig. 5.

(a) Sketch of two-dimensional granular packing in a container. (b) Illustration of the Voronoi cell surrounding a grain. The boundary of the Voronoi cell is the set of points that are at equal distance from the center of the focus grain and from the centers of the nearest neighbor grains.

Close modal
Under the above-mentioned assumptions, the joint distribution P ( v 1 , , v N ) of the N individual volumes vi is then approximated as
P ( v 1 , , v N ) 1 Z V i = 1 N exp [ v i / X ] ,
(53)
with v min < v i < v max for all i = 1 , , N to take into account the indicator function, and where X is the compactivity (ZV is a normalization constant). The individual volumes vi are, thus, approximated as independent and identically distributed random variables. The single-volume distribution P 1 ( v ) then reads
P 1 ( v ) = 1 Z 1 e v / X ( v min < v < v max ) ,
(54)
with Z1 being a normalization constant satisfying Z V = Z 1 N. One can then deduce from Eq. (54) the average value v of the individual volume, as well as possibly higher moments v n . One finds, in particular,5,28
v = 1 2 ( v min + v max ) + X Δ v coth ( Δ v X ) ,
(55)
where we have defined Δ v = ( v max v min ) / 2. As a result, one may determine the volume fraction ϕ defined as ϕ = v 0 / v (with v0 being the grain volume) as a function of the compactivity X and friction coefficient μ.

The above-mentioned model can be generalized by including two different types of grains to describe the phenomenon of granular segregation, namely, the fact that grains with different friction coefficients may gather into large clusters. The interested reader is referred to Ref. 28 for a presentation of this generalized model, which shares formal similarities with the Ising model of ferromagnetism.

We have seen that in spite of its formal similarity with the equilibrium Boltzmann–Gibbs statistics, the Edwards statistics has a strong impact on the phenomenology of frictional systems, due to its restriction to blocked configurations. Examples include the emergence of an infinite-temperature critical point in the spring–block model of Sec. IV. In practice, one of the main difficulties of the Edwards theory precisely lies in the evaluation of the indicator function F ( r 1 , , r N ) characterizing blocked configurations. The models considered in Secs. III and IV, being either without interactions or with a one-dimensional geometry, manage to keep this difficulty at a reasonable level. The study of more realistic models of granular piles requires to face this difficulty and to find appropriate approximations (see, e.g., Ref. 4 for a review). We gave an elementary example of this approach in Sec. V. Finally, we note that the Edwards approach is not the only way to tackle the statistical description of granular piles. For instance, an interesting approach based on the statistics of contacts of a given grain with neighboring grains, modeled as a stochastic process, has been shown to successfully reproduce measured density fluctuations in polydisperse emulsion packings.44 

The Edwards statistics for granular materials may be the topic of an introductory course, at a graduate level, to the statistical physics of driven athermal systems, which also deals with soft materials like foams or emulsions,45 as well as a vast variety of active matter systems ranging from active colloids to bacterial colonies and cell assemblies like biological tissues.46 The Edwards approach, whose mathematical formalism remains close to the one of equilibrium statistical physics, bridges the gap with those more advanced topics which require the development of specific and more involved methods.

The author has no conflict of interest to disclose.

1.
S.
Herminghaus
,
Wet Granular Matter: A Truly Complex Fluid
(
World Scientific
,
Singapore
,
2013
).
2.
J.
Duran
,
Sands, Powders, and Grains: An Introduction to the Physics of Granular Materials
, (
Springer
,
Heidelberg
,
2012
).
3.
Granular Matter: An Interdisciplinary Approach
, edited by
A.
Mehta
(
Springer-Verlag
,
2012
).
4.
D. P.
Bi
,
S.
Henkes
,
K. E.
Daniels
, and
B.
Chakraborty
, “
The statistical physics of a thermal materials
,”
Annu. Rev. Condens. Matter Phys.
6
(
1
),
63
83
(
2015
).
5.
S. F.
Edwards
and
R. B. S.
Oakeshott
, “
Theory of powders
,”
Physica A
157
(
3
),
1080
1090
(
1989
).
6.
A.
Mehta
and
S. F.
Edwards
, “
Statistical mechanics of powder mixtures
,”
Physica A
157
(
3
),
1091
1100
(
1989
).
7.
J.
Kurchan
and
H.
Makse
, “
Testing the thermodynamic approach to granular matter with a numerical model of a decisive experiment
,”
Nature
415
(
6872
),
614
617
(
2002
).
8.
C.
Song
,
P.
Wang
, and
H. A.
Makse
, “
Experimental measurement of an effective temperature for jammed granular materials
,”
Proc. Natl. Acad. Sci. U. S. A.
102
(
7
),
2299
2304
(
2005
).
9.
D.
Chandler
,
Introduction to Modern Statistical Mechanics
(Oxford University,
Oxford
,
1987
).
10.
D. A.
McQuarrie
,
Statistical Mechanics
(
Harper and Row
,
New York
,
1976
).
11.
T. L.
Hill
,
An Introduction to Statistical Thermodynamics
(
Dover
,
Mineola, NY
,
1986
).
12.
N. W.
Ashcroft
and
N. D.
Mermin
,
Solid State Physics
(
Saunders College Publishing
,
Philadelphia
,
1976
).
13.
P. M.
Chaikin
and
T. C.
Lubensky
,
Principles of Condensed Matter Physics
(
Cambrigde U. P
.,
Cambridge
,
1995
).
14.
A.
Barrat
,
J.
Kurchan
,
V.
Loreto
, and
M.
Sellitto
, “
Edwards' measures for powders and glasses
,”
Phys. Rev. Lett.
85
(
24
),
5034
5037
(
2000
).
15.
R.
Blumenfeld
and
S. F.
Edwards
, “
Granular entropy: Explicit calculations for planar assemblies
,”
Phys. Rev. Lett.
90
(
11
),
114303
(
2003
).
16.
R.
Blumenfeld
and
S. F.
Edwards
, “
Geometric partition functions of cellular systems: Explicit calculation of the entropy in two and three dimension
,”
Eur. Phys. J. E
19
(
1
),
23
30
(
2005
).
17.
C.
Briscoe
,
C. M.
Song
,
P.
Wang
, and
H. A.
Makse
, “
Entropy of Jammed matter
,”
Phys. Rev. Lett.
101
(
18
),
188001
(
2008
).
18.
P.
Wang
,
C. M.
Song
,
Y. L.
Jin
, and
H. A.
Makse
, “
Jamming II: Edwards' statistical mechanics of random packings of hard spheres
,”
Physica A
390
(
3
),
427
455
(
2011
).
19.
D.
Asenjo
,
F.
Paillusson
, and
D.
Frenkel
, “
Numerical calculation of granular entropy
,”
Phys. Rev. Lett.
112
(
9
),
098002
(
2014
).
20.
E. R.
Nowak
,
J. B.
Knight
,
E.
Ben-Naim
,
H. M.
Jaeger
, and
S. R.
Nagel
, “
Density fluctuations in vibrated granular materials
,”
Phys. Rev. E
57
(
2
),
1971
1982
(
1998
).
21.
M.
Schröter
,
D. I.
Goldman
, and
H. L.
Swinney
, “
Stationary state volume fluctuations in a granular medium
,”
Phys. Rev. E
71
(
3
),
030301(R)
(
2005
).
22.
F.
Lechenault
,
F.
da Cruz
,
O.
Dauchot
, and
E.
Bertin
, “
Free volume distribution and compactivity measurement in a bidimensional granular packing
,”
J. Stat. Mech.
2006
,
P07009
.
23.
S.
McNamara
,
P.
Richard
,
S.
de Richter
,
G.
Le Caër
, and
R.
Delannay
, “
Measurement of granular entropy
,”
Phys. Rev. E
80
(
3
),
031301
(
2009
).
24.
P. T.
Metzger
, “
Granular contact force density of states and entropy in a modified Edwards ensemble
,”
Phys. Rev. E
70
(
5
),
051303
(
2004
).
25.
P. T.
Metzger
and
C. M.
Donahue
, “
Elegance of disordered granular packings: A validation of Edwards' hypothesis
,”
Phys. Rev. Lett.
94
(
14
),
148001
(
2005
).
26.
M.
Pica Ciamarra
,
A.
Coniglio
, and
M.
Nicodemi
, “
Thermodynamics and statistical mechanics of dense granular media
,”
Phys. Rev. Lett.
97
(
15
),
158001
(
2006
).
27.
V.
Becker
and
K.
Kassner
, “
Protocol-independent granular temperature supported by numerical simulations
,”
Phys. Rev. E
92
(
5
),
052201
(
2015
).
28.
Y.
Srebro
and
D.
Levine
, “
The role of friction in compaction and segregation of granular materials
,”
Phys. Rev. E
68
(
6
),
061301
(
2003
).
29.
K. E.
Daniels
and
R. P.
Behringer
, “
Characterization of a freezing/melting transition in a vibrated and sheared granular medium
,”
J. Stat. Mech.
2006
,
P07018
.
30.
L. A.
Pugnaloni
,
D.
Maza
,
I.
Sánchez
,
P. A.
Gago
,
J.
Damas
, and
I.
Zuriguel
, “
Static granular packings: The relevant state variables
,”
Phys. Rev. E
82
(
5
),
050301(R)
(
2010
).
31.
J. G.
Puckett
,
F.
Lechenault
, and
K. E.
Daniels
, “
Local origins of volume fraction fluctuations in dense granular materials
,”
Phys. Rev. E
83
(
4
),
041301
(
2011
).
32.
J. G.
Puckett
and
K. E.
Daniels
, “
Equilibrating temperature-like variables in jammed granular subsystems
,”
Phys. Rev. Lett.
110
(
5
),
058001
(
2013
).
33.
S.
Henkes
,
C. S.
O'Hern
, and
B.
Chakraborty
, “
Entropy and temperature of a static granular assembly: An ab initio approach
,”
Phys. Rev. Lett.
99
(
3
),
038002
(
2007
).
34.
S.
Henkes
and
B.
Chakraborty
, “
Statistical mechanics framework for static granular matter
,”
Phys. Rev. E
79
(
6
),
061301
(
2009
).
35.
R.
Blumenfeld
and
S. F.
Edwards
, “
On granular stress statistics: Compactivity, angoricity, and some open issues
,”
J. Phys. Chem. B
113
(
12
),
3981
3987
(
2009
).
36.
R.
Blumenfeld
,
J. F.
Jordan
, and
S. F.
Edwards
, “
Interdependence of the volume and stress ensembles and equipartition in statistical mechanics of granular systems
,”
Phys. Rev. Lett.
109
(
23
),
238001
(
2012
).
37.
D. P.
Bi
,
J.
Zhang
,
R. P.
Behringer
, and
B.
Chakraborty
, “
Fluctuations in shear-jammed states: A statistical ensemble approach
,”
Europhys. Lett.
102
(
3
),
34002
(
2013
).
38.
R.
Burridge
and
L.
Knopoff
, “
Model and theoretical seismicity
,”
Bull. Seismol. Soc. Am.
57
(
3
),
341
371
(
1967
).
39.
J. M.
Carlson
and
J. S.
Langer
, “
Mechanical model of an earthquake fault
,”
Phys. Rev. A
40
(
11
),
6470
6484
(
1989
).
40.
B.
Blanc
,
L.-A.
Pugnaloni
, and
J.-C.
Géminard
, “
Creep motion of a model frictional system
,”
Phys. Rev. E
84
(
6
),
061303
(
2011
).
41.
G.
Gradenigo
,
E. E.
Ferrero
,
E.
Bertin
, and
J.-L.
Barrat
, “
Edwards thermodynamics for a driven athermal system with dry friction
,”
Phys. Rev. Lett.
115
(
14
),
140601
(
2015
).
42.
I. S.
Gradshteyn
and
I. M.
Ryzhik
,
Tables of Integrals, Series, and Products
, 5th ed. (
Academic Press
,
London
,
1994
).
43.
G.
Gradenigo
and
E.
Bertin
, “
Generalized Edwards thermodynamics and marginal stability in a driven system with dry and viscous friction
,”
Phys. Rev. E
95
(
3
),
030106(R)
(
2017
).
44.
M.
Clusel
,
E. I.
Corwin
,
A. O. N.
Siemens
, and
J.
Brujić
, “
A “granocentric” model for random packing of jammed emulsions
,”
Nature
460
(
7255
),
611
615
(
2009
).
45.
A.
Nicolas
,
E. E.
Ferrero
,
K.
Martens
, and
J.-L.
Barrat
, “
Deformation and flow of amorphous solids: Insights from elastoplastic models
,”
Rev. Mod. Phys.
90
(
4
),
45006
(
2018
).
46.
M. C.
Marchetti
,
J. F.
Joanny
,
S.
Ramaswamy
,
T. B.
Liverpool
,
J.
Prost
,
Madan
Rao
, and
R. Aditi
Simha
, “
Hydrodynamics of soft active matter
,”
Rev. Mod. Phys.
85
(
3
),
1143
1189
(
2013
).