Characterizing nonaffinity upon decompression of soft-sphere packings

TL;DR

This study investigates how nonaffine deformation contributions to the bulk modulus vary near the unjamming transition in soft sphere packings, revealing model-dependent behaviors and providing a theoretical scaling framework.

Contribution

It demonstrates that the relative nonaffine contribution to the bulk modulus can differ between models, and offers a scaling analysis to predict this behavior based on interaction potentials.

Findings

In canonical models, nonaffine bulk modulus ratio approaches a constant upon decompression.

In some models, the nonaffine contribution vanishes as unjamming is approached.

A theoretical scaling analysis explains the observed behaviors and predicts scaling laws.

Abstract

Athermal elastic moduli of soft sphere packings are known to exhibit universal scaling properties near the unjamming point, most notably the vanishing of the shear-to-bulk moduli ratio $G/B$ upon decompression. Interestingly, the smallness of $G/B$ stems from the large nonaffinity of deformation-induced displacements under shear strains, compared to insignificant nonaffinity of displacements under compressive strains. In this work we show using numerical simulations that the relative weights of the affine and nonaffine contributions to the bulk modulus, and their dependence on the proximity to the unjamming point, can qualitatively differ between different models that feature the same generic unjamming phenomenology. In canonical models of unjamming we observe that the ratio of the nonaffine to total bulk moduli $B_{na}/B$ approaches a constant upon decompression, while in other, less well-studied models, it vanishes. We show that the vanishing of $B_{na}/B$ in non-canonical models stems from the emergence of an invariance of net (zero) forces on the constituent particles to compressive strains at the onset of unjamming. We provide a theoretical scaling analysis that fully explains our numerical observations, and allows to predict the scaling behavior of $B_{na}/B$ upon unjamming, given the functional form of the pairwise interaction potential.

Equations

$$G=G_a-G_{na}B=B_a-B_{na},$$

$$\varphi_{n}(r_{ij},R_{i},R_{j})=\left\{\begin{array}[]{cc}\!\!\frac{\varepsilon}{n\lambda^{n}}\big{(}(R_{i}\!+\!R_{j})\!-\!r_{ij}\big{)}^{n}\!\!,&r_{ij}\!\leq\!R_{i}\!+\!R_{j}\\ 0\quad\quad,\!\!&r_{ij}\!>\!R_{i}\!+\!R_{j}\end{array}\right.\!\!\!,$$

$${\phi }_{\text{EXP}}(r_{ij})={\epsilon }_{ij}{e}^{-r_{ij}/{\ell }_{ij}},$$

$$\varphi_{\text{\tiny BUMP}}(r_{ij},R_{i},R_{j})=\left\{\begin{array}[]{cc}\!\!\varepsilon e^{-\left(1-\left(\frac{r_{ij}}{R_{i}+R_{j}}\right)^{2}\right)^{-1}}\!\!\!\!\!,&r_{ij}\!<\!R_{i}\!+\!R_{j}\\ 0\quad\quad,\!\!&r_{ij}\!>\!R_{i}\!+\!R_{j}\end{array}\right.\!\!\!,$$

$$G\equiv \frac{1}{V}\frac{d^{2}U}{d{\gamma }^2}\text{mbox}andB\equiv \frac{1}{V}\frac{d^{2}U}{d{\eta }^2}$$

$$H=\left(\begin{array}[]{cc}1+\eta&\gamma\\ 0&1+\eta\end{array}\right)\,.$$

$$\epsilon=\frac{1}{2}\left(H^{T}\cdot H-{\cal I}\right)=\frac{1}{2}\left(\begin{array}[]{cc}2\eta+\eta^{2}&\gamma+\gamma\eta\\ \gamma+\gamma\eta&2\eta+\eta^{2}+\gamma^{2}\end{array}\right)\,,$$

$$\delta U/V\simeq \sum _{\kappa \chi }{C}_{\kappa \chi }{ϵ}_{\kappa \chi }+\frac{1}{2}\sum _{\kappa \chi \theta \tau }{C}_{\kappa \chi \theta \tau }{ϵ}_{\kappa \chi }{ϵ}_{\theta \tau }.$$

$$G=C_{yy}+C_{xyxy},$$

$$B=C_{xxxx}+C_{yyyy}+2C_{xxyy}.$$

$$B=\frac{1}{V}\bigg{(}\sum_{i<j}\varphi^{\prime\prime}_{ij}r^{2}_{ij}-\sum_{k,l}{\bf\Xi}_{k}\cdot\mbox{\boldmath${\cal M}$}^{-1}_{kl}\cdot{\bf\Xi}_{l}\bigg{)}$$

$${\mathbf{\Xi }}_k\equiv \frac{{\partial }^{2}U}{\partial \eta \partial \text{mathitbf}x}=\sum _{i<j}{\phi }_{ij}^{\prime \prime }{r}_{ij}\text{mathitbf}{d}_k^{ij}$$

$$\text{mathitbf}{d}_k^{ij}\equiv \frac{\partial r_{ij}}{\partial \text{mathitbf}{x}_k}=({\delta }_{jk}-{\delta }_{ik})\text{mathitbf}{n}_{ij},$$

$$B_{a}\equiv\frac{1}{V}\sum_{i<j}\varphi^{\prime\prime}_{ij}r^{2}_{ij}\,,\ \ \mbox{and}\ \ B_{na}\equiv\frac{1}{V}\sum_{k,l}{\bf\Xi}_{k}\cdot\mbox{\boldmath${\cal M}$}^{-1}_{kl}\cdot{\bf\Xi}_{l}\,.$$

$$\frac{B_{na}}{B}=\frac{\Xi^{2}}{VB}\;\hat{{\bf\Xi}}\cdot\mbox{\boldmath${\cal M}$}^{-1}\cdot\hat{{\bf\Xi}}\,.$$

$$\hat{{\bf\Xi}}\cdot\mbox{\boldmath${\cal M}$}^{-1}\cdot\hat{{\bf\Xi}}=\sum_{\omega}\frac{\big{(}\mbox{\boldmath$\Psi$}_{\omega}\cdot\hat{{\bf\Xi}}\big{)}^{2}}{\omega^{2}}\,.$$

$$\sum_{\omega}\frac{\big{(}\mbox{\boldmath$\Psi$}_{\omega}\cdot\hat{{\bf\Xi}}\big{)}^{2}}{\omega^{2}}\sim\frac{1}{\omega_{0}^{2}}\int\limits_{\omega_{*}}^{\omega_{0}}D(\omega)d\omega\,,$$

$$\hat{{\bf\Xi}}\cdot\mbox{\boldmath${\cal M}$}^{-1}\cdot\hat{{\bf\Xi}}\sim 1/B\,.$$

$$\frac{B_{na}}{B}\sim \frac{{\Xi }^2}{V{B}^2}.$$

$$\mbox{\boldmath$F$}_{k}\equiv-\frac{\partial U}{\partial\mathitbf x_{k}}=-\sum_{i<j}\varphi_{ij}^{\prime}\mathitbf d^{\,ij}_{k}=\mbox{\boldmath$0$}\,,\quad{\bf\Xi}_{k}=\sum_{i<j}\varphi_{ij}^{\prime\prime}r_{ij}\mathitbf d^{\,ij}_{k},$$

$$c(r)\equiv {\phi }^{\prime }/{\phi }^{\prime \prime }r$$

$${\mathbf{\Xi }}_k$$

$$\Xi^{2}=\sum_{k}{\bf\Xi}_{k}\cdot{\bf\Xi}_{k}\sim\frac{1}{c^{4}}\bigg{(}\frac{dc}{dr}\bigg{)}^{2}\sum_{i<j}\left(\Delta r_{ij}\varphi_{ij}^{\prime}\right)^{2}\,,$$

$$\Xi^{2}\sim\frac{1}{c^{4}}\bigg{(}\frac{dc}{dr}\bigg{)}^{2}(\bar{\varphi}^{\prime})^{2}\sum_{i<j}\big{(}\Delta r_{ij}\big{)}^{2}\,.$$

$$\frac{B_{na}}{B}\sim\frac{\Xi^{2}}{VB^{2}}\sim\frac{B^{2}}{p^{2}}\bigg{(}\frac{dc}{dr}\bigg{)}^{2}\langle(\Delta r)^{2}\rangle\,.$$

$$\frac{B_{na}}{B}\sim\bigg{(}r\frac{dc}{dr}\bigg{)}^{2}=\left(1-\frac{\varphi^{\prime}(\varphi^{\prime\prime\prime}r+\varphi^{\prime\prime})}{(\varphi^{\prime\prime})^{2}r}\right)^{2}\,.$$

$$\Lambda \sim \frac{\Delta c}{p/B}\equiv \frac{\sqrt{⟨(c-⟨c⟩{)}^2⟩}}{p/B}\sim \frac{\frac{dc}{dr}\Delta r}{p/B}\sim r\frac{dc}{dr}\sim \sqrt{\frac{B_{na}}{B}},$$

$$\frac{B_{na}}{B}\sim \frac{1}{(n-1{)}^2}\to \text{mbox}constant,$$

$$\frac{B_{na}}{B}\sim \frac{1}{{\bar{r}}^2}\sim \rho \sim {\left(\frac{p}{B}\right)}^2,$$

$$\frac{B_{na}}{B}\sim \frac{p}{B},$$

Full text

I introduction

Many disordered substances such as foams, emulsions, suspensions, and granular materials can jam into a solid-like state, or rather display fluid-like behavior, depending on their confining volume or pressure O’Hern et al. (2003); van Hecke (2010); Liu and Nagel (2010). This diverse class of systems, that exhibit a continuous but abrupt transition between solid and fluid states by decompression, are often modeled theoretically and computationally by assemblies of soft repulsive particles O’Hern et al. (2003). Using these models, it has been shown that the unjamming transition is accompanied by a number of intriguing phenomena, such as the emergence of an excess of low frequency vibrational modes (O’Hern et al., 2003), diverging length scales (Ellenbroek et al., 2006; Lerner et al., 2014; Lerner, 2018) and scaling laws of elastic moduli (O’Hern et al., 2003; Ellenbroek et al., 2009). One of the hallmarks of the unjamming transition of soft sphere packings is the eventual loss of their shear rigidity — reflected by the vanishing shear-to-bulk moduli ratio $G/B$ — upon decompression Ellenbroek et al. (2009). Mean-field approaches, such as effective medium theory Makse et al. (1999); DeGiuli et al. (2014) or rigidity percolation of random networks (Jacobs and Thorpe, 1995), are unable to capture the self-organizational processes that determine the scaling behavior of $G/B$ Ellenbroek et al. (2009).

Previous work Lutsko (1989) has shown that athermal elastic moduli of disordered solids consist of two contributions with different physical origins: an affine term (also referred to often as the “Born” term Born and Huang (1954)), which captures the stiffness of the material with respect to imposed affine deformations, and a nonaffine term that accounts for additional (nonaffine) displacements of the internal degrees of freedom, that are required in order to preserve mechanical equilibrium under the imposed deformation. Thus, we decompose the shear and bulk moduli as

[TABLE]

where the subscripts ‘$a$’ and ‘$na$’ stand for the affine and nonaffine terms, respectively. The aformentioned vanishing of $G/B$ upon unjamming is seen as the near cancellation of $G_{a}$ and $G_{na}$ as the unjamming point is approached, while the nonaffine term of the bulk modulus always remains significantly smaller than the affine term Mizuno et al. (2016); Zaccone and Terentjev (2014).

In this work we show that the relative smallness of the nonaffine term of the bulk modulus $B_{na}/B$ is non-universal across various model systems of purely repulsive soft spheres in two dimensions (2D), and focus on understanding its physical origin. The scaling properties of $B_{na}/B$ have been previously explained in terms of excluded volume correlations (Schlegel et al., 2016; Zaccone and Terentjev, 2014), and in terms of the behavior of states of self stress (Wyart, 2005; Ellenbroek et al., 2009). Here we put forward a different perspective on this problem; below we show that the relative smallness of $B_{na}$ stems from two key ingredients: $(i)$ the approximate proportionality between pairwise forces and pairwise stiffnesses, that can be directly inferred from the pairwise interaction potential, and $(ii)$ that the net force acting on each particles vanishes — the mechanical equilibrium condition, also invoked in the argumentation of (Wyart, 2005; Ellenbroek et al., 2009).

In what follows we study the canonical unjamming model O’Hern et al. (2003): soft spheres that interact via a $\varphi\sim\delta^{n}$ interaction, where $\delta$ denotes the overlap between neighboring spheres, and for most of our study we focus on $n\!=\!2$ and $n\!=\!3$. We find that, in these models, the ratio $B_{na}/B$ approaches a constant in the limit of zero pressure, at the onset of the unjamming transition. We also study two other models that feature the same loss of shear rigidity as observed in the canonical models, while at the same time exhibit qualitatively different scaling behavior of the ratio $B_{na}/B$; upon approaching the unjamming point, $B_{na}/B$ gradually vanishes, in contrast with its behavior in the canonical models.

This paper is organized as follows; in Sect. II we provide descriptions of the different models that we explore and the employed numerical methods. In Sect. III we present the results from our numerical simulations regarding the unjamming phenomenology of the different models, with particular focus devoted to the ratio $B_{na}/B$. In Sec. IV we construct a string of scaling arguments, validated by our numerical simulations, that culminate in a scaling theory that fully explains the observed scaling laws of the ratio $B_{na}/B$ as unjamming is approached. In Sect. V we summarize our findings and discuss future work.

II models, methods and observables

In this Section we describe the numerical models and methods employed in our work, and define the observables of interest.

II.1 Soft sphere models

In this work we make use of four different models of soft, purely repulsive spheres in 2D. Disordered packings were generated for all models by a short high temperature equilibration, following by a minimization of the potential energy by means of the FIRE algorithm Bitzek et al. (2006). In all models except the EXP (introduced below) model, we incorporated the Berendsen barostat Berendsen et al. (1984) into our minimization algorithm, allowing us to generate packings at any desired target pressures, as explained in Brito et al. (2018). For the EXP model we employed 128-bit numerics for obtaining packings at extremely small densities. For all models, at all investigated state points, we employed 1000 independent packings of $N\!=\!1600$ particles. Plots of the employed pairwise potentials with a finite-range cutoff are presented in Fig. 1.

II.1.1 Harmonic and cubic spheres

We employed the canonical model of unjamming O’Hern et al. (2003), in which spheres of radii $R_{i}$ and $R_{j}$ interact via the pairwise potential

[TABLE]

with $\varepsilon$ and $\lambda$ denoting microscopic units of energy and length, respectively, and $r_{ij}$ is the pairwise distance between the centers of particles $i$ and $j$. In this work we mainly focus on models with $n\!=\!2$ and $n\!=\!3$, referred to in what follows as the harmonic and cubic models, respectively. We chose the radii of half the particles to be $7\lambda/5$, and the other half’s radii is $\lambda$. This model undergoes an unjamming transition at a packing fraction $\phi\!\equiv\!V/V_{p}\!\approx\!0.84$ where $V\!=\!L^{2}$ denotes the volume (in 2D), and $V_{p}\!=\!\sum_{i}\pi R_{i}^{2}$ is the volume occupied by the particles.

II.1.2 Exponential spheres

The third model we employed is a binary mixture of particles that interact via an exponentially-decaying pairwise potential

[TABLE]

where the parameters $\ell_{ij}$ and $\epsilon_{ij}$ depend on the species of the pair $i,j$, and can be found in Kooij and Lerner (2017). As shown in Kooij and Lerner (2017), this model undergoes an unjamming transition at vanishing densities $\rho\!\equiv\!N/V\!\to\!0$. We refer to this model in what follows as the EXP model.

II.1.3 The “bump” model

The fourth model we employed is a binary mixture in which pairs of particles interact via the “bump” pairwise potential

[TABLE]

where $\varepsilon$ denotes the microscopic units of energy. As in the case for the harmonic and cubic systems, half of the particles are of size $7\lambda/5$, and the other half of size $\lambda$, where $\lambda$ denotes the microscopic units of length. We refer to this model as the “bump” model. As far as we know, it has not been studied before in the context of the unjamming transition.

II.2 Observables

Athermal elastic moduli were calculated following the formalism put forward in Lutsko (1989); the key principle in deriving expressions for athermal elastic modulus is that the system remains in mechanical equilibrium under imposed deformations. We employed the definitions

[TABLE]

for the shear and bulk modulus, respectively, where $U$ is the potential energy, $V$ is the volume, $\gamma$ is the simple shear strain and $\eta$ is the expansive strain. The latter two parametrize the imposed affine transformation of coordinates $\mathitbf x\!\to\!H\!\cdot\!\mathitbf x$ in 2D as

[TABLE]

Using this transformation, the strain tensor $\epsilon$ is given by

[TABLE]

where ${\cal I}$ represents the identity tensor. In the athermal limit the potential energy density variation \delta U/V\!\equiv\!\big{(}U(\epsilon)\!-\!U(0)\big{)}/V can be approximated by a Taylor expansion in terms of the strain tensor $\epsilon$ and the general elastic coefficients $C_{\kappa\chi}$ and $C_{\kappa\chi\theta\tau}$ of the form

[TABLE]

In terms of the general coefficients $C_{\kappa\chi}$ and $C_{\kappa\chi\theta\tau}$, our definitions of shear and bulk moduli given by Eq. (5) read

[TABLE]

and

[TABLE]

We note importantly that our definition of the bulk modulus $B$ as given by Eq. (5), chosen for the sake of simplicity, slightly differs from the conventional definition of the bulk modulus $K\!\equiv\!-V\,dp/dV$, with $p\!\equiv\!-dU/dV$ denoting the pressure. The two definitions are related via $B\!=\!(K\!-\!p){\,\mathchar 22\relax\mkern-12.0mud}^{2}$, where ${\,\mathchar 22\relax\mkern-12.0mud}$ denotes the spatial dimension. As we shall show below, in the unjamming limit $p/B\!\to\!0$, then the two definitions agree (up to an unimportant factor of ${\,\mathchar 22\relax\mkern-12.0mud}^{2}$).

In our models the potential energy is given by a sum over pairwise radially-symmetric interactions, namely $U\!=\!\sum_{i<j}\varphi_{ij}(r_{ij})$; the atomistic expression for the bulk modulus writes

[TABLE]

where $\varphi^{\prime\prime}_{ij}\!\equiv\!\partial^{2}\varphi_{ij}/\partial r_{ij}^{2}$, \mbox{\boldmath{\cal M}}_{ij}\!\equiv\!\frac{\partial^{2}U}{\partial\mathitbf x_{i}\partial\mathitbf x_{j}}, and

[TABLE]

is the force (linear) response to an imposed expansive strain, written as a weighted sum over dipole vectors

[TABLE]

with $\mathitbf n_{ij}$ the unit vector pointing from particle $i$ towards particle $j$. Following the decomposition of $B$ as spelled out in Eq. (1), we identify the affine and nonaffine contributions to the bulk modulus as

[TABLE]

III Nonuniversality of nonaffinity under compressive strain

One of the key characteristics of the unjamming transition in soft sphere packings is the loss of shear rigidity captured by the vanishing of the shear-to-bulk modulus ratio $G/B$ upon approaching the transition. In Fig. 2a we plot the ratio $G/B$ as a function of $p/B$, the latter will serve in what follows as our central dimensionless control parameter of the unjamming transition, that occurs when $p/B\!\to\!0$. Our data for all four investigated models are consistent with previous results for the canonical models van Hecke (2010); Liu and Nagel (2010), namely that $G/B\!\sim\!\sqrt{p/B}$. In previous work we have also observed this scaling for soft sphere packings with inverse-power-law pairwise interactions Kooij and Lerner (2017), strongly supporting its universality.

In Fig. 2b we plot the ratio of the nonaffine to total bulk modulus $B_{na}/B$, against $p/B$; while the scaling of the shear rigidity $G/B\!\sim\!\sqrt{p/B}$ appears to be universal for soft sphere packings, $B_{na}/B$ displays very different behavior between our different models: the canonical models (harmonic and cubic) show an initial slight increase of $B_{na}/B$ upon decompression, after which the ratio saturates to an $n$-dependent constant upon approached the unjamming point as $p/B\!\to\!0$.

Remarkably, the EXP and bump models show opposite trends; as $p/B\!\to\!0$ the fraction $B_{na}/B$ vanishes instead of converging to a constant, with scalings consistent with $\sim\!(p/B)^{2}$ for the EXP model, and $\sim\!p/B$ for the bump model. Moreover, at the same high $p/B$ the nonaffine term is relatively much smaller for the EXP model compared to the canonical models, by several orders of magnitude.

In the next Section we will build a string of scaling arguments and present supporting data from our numerical tests, that will fully explain the various scalings of $B_{na}/B$ vs. $p/B$ as shown in Fig. 2.

IV Scaling theory for $B_{na}/B$

IV.1 What makes $B_{na}/B$ small?

To explain the scaling of $B_{na}/B$ upon approaching the unjamming transition, we first go back to Eq. (14) for the nonaffine term of the bulk modulus, $B_{na}$; we define $\Xi\!\equiv\!|{\bf\Xi}|$ and $\hat{{\bf\Xi}}\!\equiv\!{\bf\Xi}/\Xi$, then $B_{na}/B$ can be trivially written in the form

[TABLE]

This form makes clear that the smallness of $B_{na}/B$ can stem either from the smallness of $\Xi$ alone, or by weak coupling between the eigenfunctions \mbox{\boldmath\Psi}_{\omega} of ${\cal M}$ and $\hat{{\bf\Xi}}$.

Let us therefore first consider the quadratic form \hat{{\bf\Xi}}\!\cdot\!\mbox{\boldmath{\cal M}}^{-1}\!\cdot\!\hat{{\bf\Xi}}; using the spectral decomposition of the dynamical matrix \mbox{\boldmath{\cal M}}\!=\!\sum_{\omega}\omega^{2}\mbox{\boldmath\Psi}_{\omega}\!\otimes\!\mbox{\boldmath\Psi}_{\omega}, it can be written as

[TABLE]

It has been shown Wyart (2005); Mizuno et al. (2016) that modes near unjamming project onto pairwise directors as $\mathitbf d\cdot\Psi_{\omega}\!\sim\!\omega/\omega_{0}$, where $\mathitbf d$ is a local dipole vector as defined in Eq. (13), and $\omega_{0}\!\sim\!\sqrt{\bar{\varphi}^{\prime\prime}}$ is the square root of a characteristic pairwise stiffness. Since ${\bf\Xi}$ is a weighted sum of dipoles (see Eq. (12)), and is normalized (i.e. $\hat{{\bf\Xi}}\!\cdot\!\hat{{\bf\Xi}}\!=\!1$), we expect that \big{(}\mbox{\boldmath\Psi}_{\omega}\cdot\hat{{\bf\Xi}}\big{)}^{2}\!\sim\!(\omega/\omega_{0})^{2}/V. The sum in Eq. (16) thus follows

[TABLE]

where $D(\omega)$ is the density of states, and $\omega_{*}$ is the well-studied frequency scale above which a plateau appears in $D(\omega)$ in unjamming systems O’Hern et al. (2003); Wyart (2005); Mizuno et al. (2016); van Hecke (2010); Liu and Nagel (2010), known to universally follow $\omega_{*}/\omega_{0}\!\sim\!\sqrt{p/B}$ in decompressed soft sphere packings near unjamming.

The main support of the density of states $D(\omega)$ near unjamming is between $\omega_{*}$ and $\omega_{0}\!\gg\!\omega_{*}$, implying that, to leading order $\int_{\omega_{*}}^{\omega_{0}}D(\omega)d\omega\sim{\cal O}(1)$. Finally, as long as $B\!\gg\!B_{na}$ then in 2D $\omega_{0}\!\sim\!\sqrt{B}$ (ignoring the unimportant units of mass) and therefore Eqs. (16) and (17) imply that

[TABLE]

Combining Eqs. (18) and (15), we conclude that

[TABLE]

In Fig. 3 we test this prediction; we find that not only are $B_{na}/B$ and $\Xi^{2}/(VB^{2})$ proportional to each other, but they also share the same proportionality constant across all investigated systems, which we find to be $5/2$. The conclusion from this analysis is that the smallness of $B_{na}/B$ must stem from the smallness of the rescaled compression-induced forces $\Xi/(B\sqrt{V})$.

IV.2 The roles of mechanical equilibrium and pairwise potential

We continue the discussion by explaining how the condition of mechanical equilibrium, together with the functional form of the pairwise potential, affect the magnitude $\Xi$ of compression-induced forces ${\bf\Xi}$, and therefore control the scaling of $B_{na}/B$ near unjamming. Comparing the mechanical equilibrium equation and the definition of ${\bf\Xi}$ (see Eq. 14):

[TABLE]

we see that they have similar structures; both forces $F$ and ${\bf\Xi}$ are given by a weighted sum of dipoles $\mathitbf d$, where the weight factors are the pairwise forces $\varphi^{\prime}$ in the case of $F$, and are the products of the pairwise stiffnesses and distances $\varphi^{\prime\prime}r$ for the case of ${\bf\Xi}$. Therefore, in the special case $\varphi^{\prime}\!=\!c\,\varphi^{\prime\prime}r$ — with a dimensionless proportionally coefficient $c$ that is independent of the pairwise distance $r$, as holds for widely-employed inverse-power-law (IPL) pairwise potentials, of the form $\varphi_{\mbox{\tiny IPL}}\!\sim\!r^{-\beta}$ Kooij and Lerner (2017) — we expect {\bf\Xi}\!=\!\mbox{\boldmath0} by virtue of mechanical equilibrium.

For other pairwise potentials different from $\varphi_{\mbox{\tiny IPL}}$,

[TABLE]

is generally a function of the pairwise distance $r$, and will thus fluctuate between different pairs of interacting particles (due to fluctuations of pairwise distances), resulting in a nonzero ${\bf\Xi}$. Nevertheless, when approaching the unjamming point, spatial fluctuations in pairwise distances between interacting particles are known to decrease O’Hern et al. (2003); Jacquin et al. (2011), in some cases (depending on the form of the interacting potential, see below) giving rise to an approximate proportionality between $\varphi^{\prime}$ and $\varphi^{\prime\prime}r$, that render $\Xi$ small.

The emergence of an approximate proportionality between $\varphi^{\prime}$ and $\varphi^{\prime\prime}$ as unjamming is approached can be directly observed in our simulations. Before discussing data, we note that characteristic pairwise forces scale as $\bar{\varphi}^{\prime}\!\sim\!p\,\bar{r}$ (in 2D) with $\bar{r}\!\equiv\!1/\sqrt{\rho}$ a characteristic pairwise distance, and, similarly, characteristic pairwise stiffnesses scale as $\bar{\varphi}^{\prime\prime}\bar{r}\!\sim\!B\bar{r}$. We thus expect that characteristic proportionality coefficients $\bar{c}\!\equiv\!\bar{\varphi}^{\prime}/\bar{\varphi}^{\prime\prime}\bar{r}\!\sim\!p/B$. This expectation is validated in Fig. 4, in which the probability distributions $p(c)$ are plotted against the rescaled variable $c/(p/B)$, for the EXP (panel (a)) and cubic (panel (b)) models, as the unjamming transition is approached (i.e. $\rho\!\to\!0$ and $p\!\to\!0$ for the two models, respectively). For the EXP model, there is a clear narrowing of the distribution, showing that correlations between $\varphi^{\prime}$ and $\varphi^{\prime\prime}r$ increase as unjamming is approached, while for the cubic model the distributions approach a limit form as $p\!\to\!0$. The bump and harmonic models show similar behavior (i.e. the distributions narrow for the bump model, but approach a limit form for the harmonic model) and are therefore not shown here.

To quantify the reduction of relative fluctuations of the proportionality coefficients $c$ near unjamming, we determined the full width at half maximum (denoted in what follows by $\Lambda$) of their distributions, for all investigated models and pressures. The modes of the distributions are close to the value two, which simply originates from the factor $1/{\,\mathchar 22\relax\mkern-12.0mud}$ in the expression for the pressure. In Fig. 5 we show that $\Lambda^{2}$ has the same scaling as $B_{na}/B$ in terms of $p/B$, supporting that it is controlled by the reduction of fluctuations of the proportionality coefficients $c$.

IV.3 Scaling theory for $B_{na}/B$

To deduce the scaling of $B_{na}/B$ in terms of $p/B$ directly from the functional form of the pairwise potential, let us first express the compression-induced forces ${\bf\Xi}$ — that we have shown to control $B_{na}/B$ in Sect. IV.1 — in terms of the proportionality coefficients $c\!\equiv\!\varphi^{\prime}/\varphi^{\prime\prime}r$ as

[TABLE]

assuming fluctuations of the proportionality coefficients $c$ around their mean are small, as seen in Fig. 4 and Fig. 5. The procedure followed in Eq. (22) is akin to projecting out the so-called ‘states of self stress’ Tighe and Vlugt (2011) from the weighted sum over dipoles, which bears similarity with the arguments of Wyart (2005) regarding the behavior of the bulk modulus near unjamming. Using that $\frac{dc^{-1}}{dr}\!\propto\!\frac{1}{c^{2}}\frac{dc}{dr}$ and denoting $\Delta r_{ij}\!\equiv\!r_{ij}\!-\!\langle r\rangle$, the squared compression-induced force follows

[TABLE]

where we neglected subleading off-diagonal contributions.

At this point we exploit the observation that pairwise force fluctuations scale with the mean pairwise force, namely $\varphi^{\prime}\!\sim\!\Delta\varphi^{\prime}\!\equiv\!\varphi^{\prime}\!-\!\langle\varphi^{\prime}\rangle$, as demonstrated for all four investigated models in Fig. 6. Since $\Delta\varphi^{\prime}\!\sim\!\varphi^{\prime\prime}\Delta r$, then as long as $\varphi^{\prime\prime}\!\to\!0$ near unjamming, $\Delta\varphi^{\prime}\!\ll\!\Delta r$ (and $\Delta\varphi^{\prime}\!\sim\!\Delta r$ in the harmonic model for which $\varphi^{\prime\prime}$ is constant). We can thus neglect fluctuations of $\varphi^{\prime}_{ij}$ in Eq. (23), and only keep the leading order contribution, that reads

[TABLE]

Since $\bar{c}\!\sim\!p/B$, $\bar{\varphi}^{\prime}\!\sim\!p\,\bar{r}$, and $N/V\!\sim\!1/\bar{r}^{2}$ (in 2D), together with Eq. (19) we obtain

[TABLE]

As discussed above, characteristic fluctuations of pairwise distances follows $\Delta r\!\sim\!\Delta\varphi^{\prime}/\varphi^{\prime\prime}\!\sim\!\varphi^{\prime}/\varphi^{\prime\prime}\!\sim\!\bar{r}\,p/B$, which leads us to our key result for 2D packings

[TABLE]

Our observation that $\Lambda^{2}\!\sim\!B_{na}/B$, as seen in Fig. 5, is now readily explained; $\Lambda$ should scale as a characteristic fluctuation of $c$, rescaled by $p/B$, namely

[TABLE]

in agreement with our observation.

IV.4 Model-dependent scaling laws for $B_{na}/B$

Now we can readily explain the measured scalings laws of $B_{na}/B$ in terms of $p/B$ for the different models, as presented in Fig. 5, using Eq. (26). The harmonic and cubic interactions employ a pairwise potential of the form given by Eq. (2). Using Eq. (26) we see that near unjamming

[TABLE]

in agreement with our data of Figs. 2 and 5. Interestingly, this also shows that the nonaffine contribution to the bulk modulus near unjamming decreases for higher exponents $n$. Fig. 7 shows $B_{na}/B$ near unjamming, for systems with increasing $n$, confirming once more the predicted scaling.

For the EXP model, we make use of the relation $p/B\!\sim\!\sqrt{\rho}$ as derived in Kooij and Lerner (2017) for the EXP model in 2D, to predict that

[TABLE]

in agreement with our data of Figs. 2 and 5.

Finally, it can be shown that for the bump model near the unjamming transition $dc/dr\!\sim\!\sqrt{\varphi^{\prime}/\varphi^{\prime\prime}r}$, leading to

[TABLE]

in agreement with our data of Figs. 2 and 5.

V Discussion and outlook

In this work we investigated the relative contribution of the nonaffine term to the bulk modulus, $B_{na}/B$, across four different models of soft sphere packings in 2D, near their respective unjamming points. We find that as the critical point is approached, the relative contribution of the nonaffine term to the bulk modulus is non-universal: it can either saturate to a constant or vanish, depending on the form of the pairwise interaction potential.

In order to explain the non-universality in the observed scaling laws, we first established that the smallness of $B_{na}/B$ stems from the relative smallness of the compression-induced forces ${\bf\Xi}\!\equiv\!\partial^{2}U/\partial\eta\partial\mathitbf x$, which are a key component of $B_{na}$ (see Eq. (14)). This was done by arguing that $B_{na}/B\!\sim\!\Xi^{2}/VB^{2}$, which was also verified numerically. Interestingly, we found that not only does this scaling relation hold for all investigated systems, but also that the dimensionless proportionality coefficient of this scaling law appears to be universal across all investigated models, suggesting that it might be amenable to an exact calculation in future work.

In two of the four investigated models, namely the EXP and bump models, the compression-induced forces ${\bf\Xi}$ are shown to vanish as the unjamming point is approached. This implies that at the critical point the net (zero) force becomes invariant to compressive or expansive strains. We speculate that a close connection exists between this emergent property of the EXP and bump models, and the approximate ‘isomorph-invariance’ put forward by Dyre and co-workers Dyre (2016); Bacher et al. (2018). It has been argued that high-dimensional soft sphere packings are isomorph-invariant Maimbourg and Kurchan (2016), which might have implications on the dimensionality-dependence of nonaffinity upon compression as studied here. The detailed investigation of these connections is left for future work.

We next proceeded to construct a string of scaling arguments that lead to our key result given by Eq. (26), namely that $B_{na}/B\!\sim\!(r\,dc/dr)^{2}$, where $c(r)\!\equiv\!\varphi^{\prime}/\varphi^{\prime\prime}r$ is the proportionality coefficient between the pairwise force $\varphi^{\prime}$ and the product $\varphi^{\prime\prime}r$, where $\varphi^{\prime\prime}$ and $r$ are the pairwise stiffness and distance between interacting particles, respectively. We showed that the derived scaling relation fully explains the behavior of $B_{na}/B$ near unjamming in all four investigated models.

Our scaling argument is based on the observation, presented in Fig. 6, that the distribution of pairwise (contact) forces, rescaled by $p/\rho^{1/{\,\mathchar 22\relax\mkern-12.0mud}}$, universally assumes a finite width near the unjamming point, independent of model details, as also seen in hard sphere packings Lerner et al. (2013), and in colloidal glasses near jamming Charbonneau et al. (2012). At this point we cannot argue why this must always be the case. An interesting result of this observation, however, is that fluctuations in pairwise distances follow $\Delta r\!\sim\!\bar{r}p/B$, where $\bar{r}$ is a characteristic pairwise distance. In the EXP model $\bar{r}\!\sim\!B/p$ Kooij and Lerner (2017), indicating that, in this particular model, pairwise fluctuations are independent of density.

In two of the four investigated models, namely the EXP and bump models, the finite width of the rescaled pairwise forces translates to a vanishing width of the rescaled dimensionless proportionality coefficients $c/(p/B)\!\sim\!c/\langle c\rangle$. This implies, remarkably, that while pairwise force fluctuations scale as the mean pairwise force, namely $\Delta\varphi^{\prime}\!\sim\!\langle\varphi^{\prime}\rangle$, this is not the case for the proportionality coefficients $c$, whose fluctuations to mean ratio $\Delta c/\langle c\rangle\to 0$ as unjamming is approached. In future work the implications of these reduced fluctuations might be explored.

Acknowledgements.

We would like to thank SURFsara for the support in using the Lisa Computer Cluster. This work is part of the Industrial Partnership Programme Hybrid Soft Materials that is carried out under an agreement between Unilever Research and Development B.V. and the Netherlands Organisation for Scientific Research (NWO). EL acknowledges support from the NWO (Vidi grant no. 680-47-554/3259).