# A microscopic mean-field theory of the jamming transition

Hugo Jacquin, Ludovic Berthier, Francesco Zamponi

► To cite this version:

Hugo Jacquin, Ludovic Berthier, Francesco Zamponi. A microscopic mean-field theory of the jamming transition. Physical Review Letters, 2011, 106, pp.135702. 10.1103/PhysRevLett.106.135702. hal-00597233

HAL Id: hal-00597233

https://hal.science/hal-00597233v1

Submitted on 7 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
A microscopic mean-field theory of the jamming transition

Hugo Jacquin,¹ Ludovic Berthier,² and Francesco Zamponi³

¹Laboratoire Matière et Systèmes Complexes, UMR CNRS 7057, Université Paris Diderot – Paris 7,
10 rue Alice Domon et Léonie Duquet, 75205 Paris cedex 13, France

²Laboratoire Charles Coulomb, UMR 5221 CNRS and Université Montpellier 2, Montpellier, France

³Laboratoire de Physique Théorique, École Normale Supérieure,
UMR CNRS 8549, 24 Rue Lhomond, 75231 Paris Cedex 05, France

Dense particle packings acquire rigidity through a nonequilibrium jamming transition commonly observed in materials from emulsions to sandpiles. We describe athermal packings and their observed geometric phase transitions using equilibrium statistical mechanics and develop a fully microscopic, mean-field theory of the jamming transition for soft repulsive spherical particles. We derive analytically some of the scaling laws and exponents characterizing the transition and obtain new predictions for microscopic correlation functions of jammed states that are amenable to experimental verifications, and whose accuracy we confirm using computer simulations.

About 50 years ago, Bernal [1] used dense disordered sphere packings as model systems to understand the liquid state, at a time where the statistical mechanics of liquids was still in its infancy. Today, the idea that jammed materials share deep similarities with dense liquids and glasses remains popular [2]. However, while liquid state theory grew as a cornerstone of theoretical physics [3], no equivalent theory is available for jammed matter, because this is a nonequilibrium, amorphous, athermal state of matter—a theoretical challenge overlooked by Bernal. Thus, despite intense research activity [4–6] with a large body of numerical and experimental observations [7, 8], it is not yet clear what is the appropriate theoretical framework to understand dense athermal packings and the intriguing phase transitions they undergo, although foams, pastes and emulsions are familiar materials.

We address the purely geometrical packing problem of soft spheres, and suggest to study first their statistical mechanics at finite temperatures, $T$, before taking the $T \to 0$ limit where jamming occurs. A similar approach is frequently used in combinatorial optimization problems [9], because powerful statistical mechanics tools can then be used in a context where they are not a priori relevant [10]. We investigate the statistical mechanics of the system at $T \ge 0$ using mean-field theory [6, 11], and develop a fully microscopic theoretical scheme to predict the structure of non-equilibrium configurations of soft repulsive spheres at zero [12] and finite [13] temperatures.

Our microscopic approach is thus markedly different from recent theoretical works [4, 5], which are based on phenomenological and scaling considerations. Similarly to Landau-Ginzburg theory of phase transitions, our aim is to derive, from first principles, the correct qualitative description of the transition, and accurate quantitative predictions for several observables. However, as a mean-field theory, our approach does not describe well all fluctuations near the transition, and the associated scaling laws [5, 7]: these should be included by developing a renormalization group treatment around mean-field theory, which is currently under construction [14].

To make our approach concrete, we study an assembly of $N$ spherical particles of diameter $\sigma$ enclosed in a volume $V$ in three spatial dimensions, interacting with a soft repulsion of finite range. To fix ideas we choose

$$V(r \le \sigma) = \epsilon(1 - r/\sigma)^{\alpha}, \quad V(r > \sigma) = 0, \quad (1)$$

with $r$ the interparticle distance, $\epsilon$ the strength of the repulsion, and $\alpha = 2$ (harmonic repulsion). Although several systems are described by a Hertzian repulsion ($\alpha = \frac{3}{2}$), the harmonic model originally proposed to describe wet foams [15] has become a paradigm in numerical studies of the $T = 0$ jamming transition [7, 8]. It was also studied at finite temperatures [13, 16], and finds experimental realizations in emulsions and soft colloids. The choice $\alpha = 2$ is also technically more convenient, but we emphasize that our approach is easily generalized to any repulsive potential. The model has two control parameters: the temperature $T$, and the fraction of the volume occupied by the particles in the absence of overlap: $\varphi = \pi N \sigma^3/(6V)$. We set $\sigma$ and $\epsilon$ to unity.

Over the last decade, a number of numerical observations were reported for this model [8]. A jamming transition is observed at $T = 0$ at some critical volume fraction, $\varphi_j$, the density above which packings carry a finite density of particle overlaps. Numerically, energy density, $e_{gs}$, and pressure, $P$, are found to increase continuously from zero above $\varphi_j$ as power laws [7]. The pair correlation function of density fluctuations [3], $g(r)$, develops singularities near $\varphi_j$ [12], which are smoothed by thermal excitations [13]. In particular, $g(1) = \infty$ at $\varphi_j$ and $T = 0$, which implies that the density of contacts between particles, $z$, jumps discontinuously from 0 to a finite value, $z_c$, at $\varphi_j$. Above $\varphi_j$, $z$ increases algebraically with $\varphi$ [7, 8]. Thus, jamming appears as a phase transition taking place in the absence of thermal motion, with a peculiar critical behaviour and observable physical consequences [8].

The success of our approach relies on its ability to accurately describe dense systems of harmonic spheres at very low $T$, which is theoretically challenging [17]. Simple liquid state theories, such as integral equations [3], work
<figure><img src="image_2.png" /><figcaption>FIG. 1: Sketch of the derivation of the replicated free energy and effective potential in Eq. (4). Each particle in the original liquid is replicated m times (dashed spheres). Assuming that the replicated particles form a “molecule” of average cage size A, we trace out in the partition sum the degrees of freedom of (m − 1) copies of the liquid to obtain an effective one-component liquid (black spheres) interacting with an effective pair potential V$_{eff}$ ( r ) (green lines).</figcaption></figure>

well for dense systems only when T is not too low [18]. At lower T , numerical simulations [16] indicate the ap- pearance of a complex free energy landscape associated with slow dynamics, as found in glass-forming liquids. Theoretically, the mode-coupling theory of glasses can be applied [19] but only gives limited dynamical insights, in particular failing to identify the jamming transition [17]. For structure, a possible path is a recently developed an- alytical approach based on replica calculations [11, 20]. While generic in principle, the method requires in prac- tice specific approximations. We have found that prior works on Lennard-Jones [11] and hard sphere [6] models fail when directly applied to models of harmonic spheres near jamming [17]. The main technical achievement of the present paper is the derivation of a new analytical scheme to compute the structure and thermodynamics of the model (1) over a range of parameters broad enough to allow the study of the T = 0 jamming transition.

We introduce m copies of the system of harmonic spheres as a mathematical tool to probe its complex free energy landscape [20], and develop approximations to study the statistical mechanics of the replicated liquid using an effective potential approach [6], as sketched in Fig. 1. We make a Gaussian ansatz for the probability distribution of replicated particles around the center of mass of the “molecules” shown in Fig. 1, ρ ( x$_{1}$ · · · x$_{m}$ ) = ∫ d $^{3}$X ∏ m a =1 (2 πA ) − 3 / $^{2}$e − ( x$_{a}$ − x ) $^{2}$/ (2 A $^{)}$, which defines the cage size A . Our central approximation is now per- formed, in which only two-body interactions between par- ticles in copy 1 induced by the coupling to the other ( m − 1) copies are retained, see Fig. 1. Consider two “molecules”, each composed of m particles with positions ( x$_{1}$ · · · x$_{m}$ ) and ( y$_{1}$ · · · y$_{m}$ ): the effective potential be-

tween the particles of replica 1 is obtained by averaging the total interaction ∑ m a =1 V ( x$_{a}$ − y$_{a}$ ) over the positions of particles within the ( m − 1) remaining replicas:

$$
\begin{equation}
\begin{split}
e^{-\beta V_{\text{eff}}(\mathbf{x}_1 - \mathbf{y}_1)} &\equiv \int d^3\mathbf{x}_2 d^3\mathbf{y}_2 \cdots d^3\mathbf{x}_m d^3\mathbf{y}_m \\
&\quad \times \left\{ \rho(\mathbf{x}_1 \cdots \mathbf{x}_m) \rho(\mathbf{y}_1 \cdots \mathbf{y}_m) \prod_{a=1}^m e^{-\beta V(\mathbf{x}_a - \mathbf{y}_a)} \right\}.
\end{split}
\tag{2}
\end{equation}
$$

Thanks to the Gaussian form of the integral the latter expression can be rewritten as follows:

$$
e^{-\beta V_{\text{eff}}(r)} = \frac{e^{-\beta V(r)}}{r\sqrt{4\pi A}} \int_0^{\infty} du \left[ e^{-\frac{(r-u)^2}{4A}} - e^{-\frac{(r+u)^2}{4A}} \right] u q^{m-1}(u), \quad (3)
$$

where q ( u ) = ∫ d $^{3}$t e − $^{βV (}$u − t ) e − t $^{2}$/ (4 A ) / (4 πA ) 3 / 2 has an explicit expression in terms of error functions, and β = 1 /T . Finally, the free energy F ( m, A ; ϕ, T ) is obtained by considering V$_{eff}$ ( r ) − mV ( r ) as small, which becomes ex- act when A → 0, see Eq. (3), and doing standard pertur- bation theory [3] around the liquid with potential mV ( r ), which is equivalent to a liquid with potential V ( r ) at tem- perature T/m . We obtain an effective one-component system with a free energy parametrized by A and m :

$$
\begin{equation}
\begin{split}
F(m, A; \varphi, T) = F_{\text{harm}}(m, A) + F_{\text{liq}}(\varphi, \frac{T}{m}) & \\
& - \frac{3\varphi T}{\pi} \int dr g_{\text{liq}}(r, \varphi, \frac{T}{m}) [e^{-\beta[V_{\text{eff}}(r) - mV(r)]} - 1],
\end{split}
\tag{4}
\end{equation}
$$

where F$_{liq}$ ( ϕ, T ) and g$_{liq}$ ( r, ϕ, T ) are respectively the free energy and pair correlation function of the original (non- replicated) fluid, and F$_{harm}$ = − 3 T 2 [( m − 1) ln(2 πA ) + m − 1 + ln m ] is the ideal gas contribution for the repli- cated system [11]. Thus the core of the approximation is embodied by the effective potential V$_{eff}$ ( r ). Physically, the presence of ( m − 1) replicas induces near jamming a strong short-range effective attraction, similar in spirit to depletion forces in colloid-polymer mixtures [6].

Our task becomes the study of a complicated effective fluid described by Eq. (4). To simplify calculations we perform a standard approximation,

$$
g_{\text{liq}}(r, \varphi, T) \equiv e^{-\beta V(r)} y(r, \varphi, T) \approx e^{-\beta V(r)} y(1, \varphi, 0), \quad (5)
$$

well-suited to study the T → 0 limit [3]. To compute F$_{liq}$ and y (1 , ϕ, 0) analytically we choose the hypernetted chain approximation [3], although more elaborate closure relations [3] could be used. This could change slightly the location of the transition, but not its nature or the scaling predictions we derive. Finally, to obtain concrete results for a given state point ( ϕ, T ), we minimize the free energy with respect to the cage size A and, in the glass phase, to the replica number m [20].

We first determine the location of the transition be- tween the fluid and glass phases, signaled by the appear- ance of a free energy minimum with m < 1 [11, 20],
<figure><img src="image_3.png" /><figcaption>FIG. 2: Theoretical phase diagram of soft repulsive spheres. The glass transition temperature $T_K$ separates the liquid and glass phases with $T_K \sim (\varphi - \varphi_K)^2$ near $\varphi_K \approx 0.577$. At $T=0$, the glass jams under compression across $\varphi_{gcp} \approx 0.633$, above which no glass state with no particle overlap exists at $T=0$. Thus, the ground state energy $e_{gs}$ increases continuously from 0 as $e_{gs} \sim (\varphi - \varphi_{gcp})^2$.</figcaption></figure>

see Fig. 2. A finite temperature glass transition, $T_K$, emerges continuously from zero above $\varphi_K \approx 0.577$, as $T_K \sim (\varphi - \varphi_K)^2$. We obtain the full thermodynamic behaviour in the glass phase (notably energy, pressure, specific heat, glass fragility) which compares qualitatively well with numerical results [16]. In particular, the ground state energy and pressure remain zero across $\varphi_K$, showing that just above $\varphi_K$, $T=0$ glasses are not jammed. In these glassy states, like in a hard sphere crystal, particles can vibrate near well-defined (but random) positions, and the system is not jammed [6].

We now concentrate on the $T \to 0$ limit at large volume fraction in the glass phase. We obtain the ground state energy shown in Fig. 2. As found in simulations [7], it grows continuously from zero above a critical packing fraction, $e_{gs} \sim (\varphi - \varphi_{gcp})^2$, so that the pressure increases linearly, $P \sim (\varphi - \varphi_{gcp})$. The ‘glass close packing’ [6], $\varphi_{gcp}$, represents in our calculations the largest density where $T=0$ glasses with no particle overlap exist. Within the present approximation we obtain $\varphi_{gcp} = 0.633 > \varphi_K$. Thus, from the sole knowledge of $V(r)$ in Eq. (1), our theory predicts the existence and location of a jamming transition deep in the glass phase, and accounts for its critical nature. Since a large number of metastable states exist in the glass phase, our approach also directly explains the strong protocol dependence of the critical jamming density $\varphi_j$ observed in simulations [7, 21], which get arrested in nonequilibrium amorphous states, and thus jam at a critical packing fraction $\varphi_j < \varphi_{gcp}$. However, the results we obtain near $\varphi_{gcp}$ below are found to hold for any metastable glass, and therefore also hold near any protocol-dependent $\varphi_j$.

<figure><img src="image_4.png" /><figcaption>FIG. 3: Evolution of the maximum of the pair correlation function near contact with $T$ and $\varphi$. While $g_{max}$ diverges on both sides of the transition at $T=0$ as $g_{max} \sim |\varphi - \varphi_{gcp}|^{-1}$, this divergence becomes a smooth maximum at finite $T$ near the transition whose position shifts as $\sqrt{T}$ (dashed line), as observed numerically [12, 13] and experimentally [13, 22].</figcaption></figure>

We now turn to the calculation of the pair correlation function near $\varphi_{gcp}$. Within our approximation, $g(r)$ is directly related to the effective potential:

$$g(r) = e^{-\beta V_{\text{eff}}(r)} y(1, \varphi, 0), \quad (6)$$

and comes as a direct result of the free energy minimization. We concentrate on the physics of interparticle contacts and thus focus on distances $r \approx 1$, in the vicinity of the jamming transition ($T \ll 1$, $\varphi \approx \varphi_{gcp}$).

At $T=0$, we find that $g(r)$ develops a diverging peak near contact, which obeys the following scaling law:

$$g(r) \approx |\delta\varphi|^{-1} \mathcal{F}_{\pm} \left[ \frac{r-1}{|\delta\varphi|} \right], \quad (7)$$

where $\mathcal{F}_{\pm}(x)$ are asymmetric scaling functions which depend on the sign of $\delta\varphi \equiv \varphi - \varphi_{gcp}$ and can be computed analytically. In particular, $\log \mathcal{F}_{+}(x) \propto -x^2$ and $\mathcal{F}_{-}(x) \propto x^{-2}$ when $x \gg 1$. Note that $\mathcal{F}_{-}(x)$ can be derived using the hard sphere potential [6]. The scaling (7) means that the peak height, $g_{max}$, diverges as $|\delta\varphi|^{-1}$ on both sides of the transition at $T=0$, see Fig. 3, while its width vanishes as $|\delta\varphi|$. This behaviour was found in simulations [12]. In Fig. 4 we show that numerical results not only obey the scaling form in Eq. (7), but also that the asymmetric shape of the scaling functions compares extremely well with our theoretical predictions.

At very low, but finite temperature, the peak divergence is smoothed by thermal fluctuations, which was the focus of a recent study [13]. In Fig. 3 we show the predicted smooth evolution of $g_{max}$ when $\varphi_{gcp}$ is crossed at finite $T$. A nonmonotonic evolution with density is obtained, as in experiments [13, 22]. The position of the maximum evolves with $T$ with a scaling in perfect agreement with numerical work [13]. We go further and predict the thermal broadening of $g(r)$, see Fig. 4. We obtain nearly perfect agreement of theory with simulations
<figure><img src="image_5.png" /><figcaption>FIG. 4: The pair correlation near jamming predicted by the theory (full lines) and measured in numerical simulations (symbols). Left panels: scaling behaviour at $T = 0$ above (top) and below (bottom) the jamming transition showing the convergence of the first peak near $r = 1$ to a delta function with asymmetric scaling functions on both sides of the transition. Right panels: the first peak of the pair correlation broadens when $T$ increases at constant $\varphi$ above (top: $\delta\varphi = 2.8 \cdot 10^{-4}$) and below the transition (bottom $\delta\varphi = -3.5 \cdot 10^{-4}$). To ease visualization, we show the evolution of $g(r)/g_{max}$, where $g_{max}$ can be read from Fig. 3.</figcaption></figure>

over several decades of temperatures with no adjustable parameter. The full scaling of $g(r)$ near contact, as a function of $T$ and $\delta\varphi$ near the jamming transition, is the main new achievement of the present work.

Finally we obtain the number of contacts per particle by integration: $z = 24\varphi \int_0^\infty dr r^2 g(r)$. The diverging peak described by Eq. (7) gives a discontinuous jump of $z$ from 0 to $z_c = 6$, the celebrated isostatic value, at $\varphi_{gcp}$, as observed [1, 7, 8] and already derived in [6]. Above the transition we find $z-z_c \propto (\varphi-\varphi_{gcp})^\gamma$ with $\gamma = 1$. The exponent is in quantitative disagreement with the observed $\gamma = \frac{1}{2}$ [7]. Indeed, this exponent has been related to the presence of fluctuations [5] that are presumably not well captured by our mean-field theory, indicating that more detailed calculations (possibly based on renormalization group [14]) should be developed to predict $g(r)$ over a broader range of interparticle distances.

Our results show that the fully nonequilibrium problem of soft particle packings, relevant to understanding the mechanical properties of many soft materials, can be successfully address using equilibrium statistical mechanics tools. As an application we have developed a many-body theory of the jamming transition of soft repulsive spheres which satisfactorily derives, from first principles, the existence and location of a jamming transition, and some of its peculiar critical behaviour, and makes new predictions for correlation functions of jammed states. Our approach

is general enough that it can be systematically improved and generalized to various models, such that new, or more precise predictions could be made, hopefully fostering more numerical or experimental work. While Bernal saw packings as simplified models for atomic liquids, it is equally useful to consider packings as a special class of disordered ground states.

We thank P. Chaudhuri for providing the configurations used in numerical verifications of theoretical predictions in Fig. 4. H. Jacquin acknowledges financial support from Capital Fund Management (CFM) Foundation. L. Berthier is partially funded by ANR Dynhet and Région Languedoc-Roussillon.

[1] J. D. Bernal, Nature **183**, 141 (1959); J. D. Bernal and Mason, Nature **188**, 910 (1960).

[2] A. J. Liu and S. R. Nagel, Nature **396**, 21 (1998).

[3] J. P. Hansen and I. R. McDonald, *Theory of Simple Liquids* (Elsevier, Amsterdam, 1986).

[4] M. Clusel, E. I. Corwin, A. O. N. Siemens, and J. Brujic, Nature **460**, 611 (2009); C. Song, P. Wang, and H. A. Makse, Nature **453**, 629 (2008).

[5] M. Wyart, L. Silbert, S. R. Nagel, and T. Witten, Phys. Rev. E **72**, 051306 (2005).

[6] G. Parisi and F. Zamponi, Rev. Mod. Phys. **82** 789 (2010).

[7] C. S. O'Hern, S. A. Langer, A. J. Liu, and S. R. Nagel, Phys. Rev. Lett. **88**, 075507 (2002).

[8] M. van Hecke, J. Phys. Condens. Matter **22**, 033101 (2010); A. J. Liu, S. R. Nagel, W. van Saarloos, and M. Wyart, arXiv:1006.2365.

[9] M. Mézard and A. Montanari, *Information, Physics, and Computation* (Oxford University Press, Oxford, 2009).

[10] F. Krzakala and J. Kurchan, Phys. Rev. E **76**, 021122 (2007).

[11] M. Mézard and G. Parisi, J. Chem. Phys. **111**, 1076 (1999).

[12] A. Donev, S. Torquato, and F. H. Stillinger, Phys. Rev. E **71**, 011105 (2005); L. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E **73**, 041304 (2006).

[13] Z. Zhang et al., Nature **459**, 230 (2009).

[14] M. Castellana et al., Phys. Rev. Lett. **104**, 127206 (2010); C. Cammarota, G. Biroli, M. Tarzia, G. Tarjus, Phys. Rev. Lett. **106**, 115705 (2011).

[15] D. J. Durian, Phys. Rev. Lett. **75**, 4780 (1995).

[16] L. Berthier and T. A. Witten, EPL **86**, 100001 (2009); Phys. Rev. E **80**, 021502 (2009).

[17] L. Berthier, H. Jacquin and F. Zamponi, J. Stat. Mech. (2011) P01004.

[18] H. Jacquin and L. Berthier, Soft Matter **6**, 2970 (2010).

[19] L. Berthier, E. Flenner, H. Jacquin, and G. Szamel, Phys. Rev. E **81**, 031505 (2010); W. T. Kranz, M. Sperl, and A. Zippelius, Phys. Rev. Lett. **104**, 225701 (2010).

[20] R. Monasson, Phys. Rev. Lett. **75**, 2847 (1995).

[21] P. Chaudhuri, L. Berthier, and S. Sastry, Phys. Rev. Lett. **104**, 165701 (2010).

[22] X. Cheng, Phys. Rev. E **81**, 031301 (2010).