Contact changes of sheared systems: Scaling, correlations, and mechanisms

We probe the onset and effect of contact changes in two-dimensional soft harmonic particle packings which are sheared quasistatically under controlled strain. First, we show that, in the majority of cases, the first contact changes correspond to the creation or breaking of contacts on a single particle, with contact breaking overwhelmingly likely for low pressures and/or small systems, and contact making and breaking equally likely for large pressures and in the thermodynamic limit. The statistics of the corresponding strains are near-Poissonian, in particular for large-enough systems. The mean characteristic strains exhibit scaling with the number of particles $N$ and pressure $P$ and reveal the existence of finite-size effects akin to those seen for linear response quantities [C. P. Goodrich , Phys. Rev. Lett. 109, 095704 (2012); C. P. Goodrich et al., Phys. Rev. E 90, 022138 (2014)]. Second, we show that linear response accurately predicts the strains of the first contact changes, which allows us to accurately study the scaling of the characteristic strains of making and breaking contacts separately. Both of these show finite-size scaling, and we formulate scaling arguments that are consistent with the observed behavior. Third, we probe the effect of the first contact change on the shear modulus $G$ and show in detail how the variation of $G$ remains smooth and bounded in the large-system-size limit: Even though contact changes occur then at vanishingly small strains, their cumulative effect, even at a fixed value of the strain, are limited, so, effectively, linear response remains well defined. Fourth, we explore multiple contact changes under shear and find strong and surprising correlations between alternating making and breaking events. Fifth, we show that by making a link with extremal statistics, our data are consistent with a very slow crossover to self-averaging with system size, so the thermodynamic limit is reached much more slowly than expected based on finite-size scaling of elastic quantities or contact breaking strains.

Figure 1

(a) The first contact change in a sheared packing ($N=64, P={10}^{-6}$) occurs at a strain ${\gamma }_{*}=9.003851\left(2\right)\times {10}^{-7}$, when the two marked particles lose their contact. (b) The corresponding stress-strain curve remains continuous but exhibits a sharp kink; we define $G_0$ as the shear modulus of the undeformed packing and $G_1$ as the shear modulus of the packing just above ${\gamma }_{*}$.

Figure 2

(a) Zoom-in of a packing where particle $r$ becomes a rattler after the first contact change ($N=22, P=1.5\times {10}^{-5}$). Neighboring particles $A, B$, and $C$ are indicated. (b) Overlap of $r$ with the neighboring particles $A, B$, and $C$ as a function of strain $\gamma$. Markers are direct numerical simulation (DNS) data points, lines indicate the linear response prediction. The DNS and LR predictions for rattler creation are ${\gamma }_{*}^{\mathrm{DNS}}=2.45\times {10}^{-5}$ and ${\gamma }_{*}^{\mathrm{LR}}=2.41\times {10}^{-5}$.

Figure 3

(a) Cumulative distribution functions $Pr({\gamma }_{*}<\gamma )$ of the contact change strain ${\gamma }_{*}$ for $N=256, P={10}^{-6}$ (left) $...{10}^{-2}$ (right). (b) Stacked probabilities for the first contact change being a break event (blue striped), a make event (red striped), a make event involving a rattler (red), and a mixed event, where contacts are both broken and created (black), for $N=256$ ensembles.

Figure 4

(a) Rescaled complementary cumulative distribution functions (CCDF), $N=256, P={10}^{-6}...{10}^{-2}$ (highest pressures have lowest values for $k=0.1$). The dotted line gives the CCDF for an exponential distribution. (b) Same as in (a) for $N=16$ systems at various pressures. Inset: Result of the Anderson-Darling test. Ensembles that fail the test are indicated with a red cross. Other ensembles are indicated with a dot ($\le 100$ samples), open circle (100–1000 samples), or filled circle ($\approx 1000$ samples). The blue line indicates the finite-size threshold $N^{2}P{log}_{10}{\left(N\right)}^{-0.7}=1$ (see Sec. 3c).

Figure 5

(a) Scaling of the strain at first contact change ${\gamma }_{\mathrm{cc}}$ as function of $N$ and $P$. Symbols and shades (colors) indicate packing sizes. Lines indicate power-law functions with exponent 1 (lower branch) and 0.5 (upper branch). (b) Log corrections improve the collapse. Inset: Probability of the first contact change creating a new contact. At high $N^{2}P{log}_{10}{\left(N\right)}^{-0.7}, Pr\left(\mathrm{mk}\right)\approx Pr\left(\mathrm{bk}\right)\approx 0.5$, but at low $N^{2}P{log}_{10}{\left(N\right)}^{-0.7}$, breaking strongly dominates.

Figure 6

Excess number of contacts $N\mathrm{\Delta }z/2$ as a function of $N^{2}P$ (blue curve, based on Refs. [31, 32, 33]). Arrows indicate volumetric strains corresponding to a single contact change.

Figure 7

${\sigma }_{\mathrm{cc}}$ as a function of $P/N$. Symbols and shades (colors) indicate system size. The data support an overall scaling $\sigma \sim P/N$, but the lack of a good collapse suggests this does not describe the entire behavior—larger $N$ tend to have lower ${\sigma }_{cc}$.

Figure 8

Stress response for (a) a packing with typical $Q=0.014$ and (b) a packing with a very strong nonlinearity ($Q=0.267$; both $N=16,P={10}^{-2}$). The simulation data (red $\times$) are fitted with the second-order polynomial (blue solid curves) $\sigma =G_2\gamma +\lambda {\gamma }^2$. The black dotted curves are the linear contribution $\sigma =G_2\gamma$; the green dash-dotted curves are the linear response predictions $\sigma =G_{\mathrm{LR}}\gamma$. The gray vertical lines indicate the strain at the first contact change ${\gamma }_{*}$.

Figure 9

(a) Width of the distribution of $Q$, the relative deviation from linear response at the first contact change for different ensembles. Clearly, the stresses are very well described by linear response for small systems at low pressures. More significant deviations occur for small systems at high pressures and large systems at low pressures. (b) Width of the distribution of $\lambda$. The quadratic component is on the order of ${10}^{-1}$ in most cases, but increases for large systems at low pressures (i.e., close to jamming).

Figure 10

PDFs of ${\gamma }_{*}^{\mathrm{LR}}/{\gamma }_{*}^{\mathrm{DNS}}$ for various system sizes as indicated at (a) $P={10}^{-2}$ and (b) $P={10}^{-6}$. For each PDF, the standard deviation $\sigma$ is indicated.

Figure 11

(a) Scaling of ensemble averaged breaking ($▽$) and making ($△$) strains ${\gamma }_{\mathrm{bk}}^{\mathrm{LR}}=\langle {\gamma }_{*,\mathrm{mk}}^{\mathrm{LR}}\rangle$ and ${\gamma }_{\mathrm{mk}}^{\mathrm{LR}}=\langle {\gamma }_{*,\mathrm{mk}}^{\mathrm{LR}}\rangle$ from linear response. (b) As in Fig. 5, log corrections significantly improve the quality of the collapse. Inset: $Pr{\left(\mathrm{mk}\right)}^{\mathrm{LR}}=1/(1+{\gamma }_{\mathrm{mk}}/{\gamma }_{\mathrm{bk}})$ is approximately 0.5 for high $N^{2}P{log}_{10}{\left(N\right)}^{-0.7}$ and scales as ${\gamma }_{\mathrm{bk}}$ for small $N^{2}P{log}_{10}{\left(N\right)}^{-0.7}$.

Figure 12

Alternate scaling plots for $N^q\gamma \sim F\left(N^{r}P\right)$, with $q=1.8...2.2$ (vertical) and $r=1.6...2.0$ (horizontal). $\aleph$ is a measure for collapse quality (see text), and the best collapse is found for $N^2\gamma \sim F\left(N^{1.8}P\right)$ (blue border), and all collapses with $0\le r-q\le 0.4$ are reasonable (green border).

Figure 13

(a) Asymptotical behavior of $F\left(x\right)$. Black lines show the result from the power-law fit: $F\left(x\right)\sim x^{1.0}$ for low $x$ and $F\left(x\right)\sim x^{0.5}$ for high $x$. The crossover between the two regimes is at $x=0.4$. (b) Residual plot $F\left(x\right)/x^{1.0}$ [(dark)blue] and $F\left(x\right)/x^{0.5}$ [(light)red] show the fitted power laws match the behavior very well in their respective regimes, as they scatter around a constant value. (c) Log-correction $c\left(N\right)={log}_{10}{\left(N\right)}^{-0.7}$ and power-law correction $c\left(N\right)=N^2/N^{1.8}=N^{-0.2}$ as function of system size $N$. Both vary roughly by a factor of two in the range of $N$ we probe. (d) The ratio of the two varies by less than 35%.

Figure 14

(a) Probability distribution functions for $G_1/G_0$, the relative shear modulus after the first contact change. For small systems at low pressures (bottom), we find $0\le G_1/G_0\le 1$; for intermediate systems we find $G_1/G_0$ is typically smaller than 1 but can become negative (indicating an unstable system). For large systems at high pressures (top), we find $G_1/G_0\approx 1$. The creation of contacts [(dark)blue] correlates with an increase in $G$, while the breaking of contacts [(light)red] correlates with a decrease in $G$. (b) The fraction of events where $G_1<0$ peaks around $N^{2}P{log}_{10}{\left(N\right)}^{-0.7}\approx 1$. (c) The standard deviation of $G_1/G_0$. For small systems at low pressures, $\sigma \approx 0.3$, whereas for large systems $\sigma \sim {\left(N^{2}P\right)}^{-\beta }$ with $\beta =0.35\pm 0.01$.

Figure 15

The number of contacts, $N_c$, for systems with $N=16$ particles at $P={10}^{-6}$ as function of the cumulative number of contact changes. The circle area represents the fraction of systems with a given number of contacts; the thickness of the lines represent transition probabilities. Initially, the systems start off with the minimum number of contacts $2N+1=33$ (31 or 29 when there are one or two rattlers, respectively). In the first and second contact change, the system loses one contact (three when a rattler is created). In the following events, the system alternately gains and loses a contact.

Figure 16

(a) CDF of ${\gamma }_{ij}^{†}$ for every contact for every packing in the $N=1024, P={10}^{-2}$ ensemble. The strain at which $Pr({\gamma }_{ij}^{†}<{\gamma }_{\mathrm{bk}}^{\mathrm{dist}})=1/\langle N_c\rangle$ is the expected contact breaking strain for this ensemble: ${\gamma }_{\mathrm{bk}}^{\mathrm{dist}}=1.4\times {10}^{-4}$. The mean breaking strain from linear response is ${\gamma }_{\mathrm{bk}}^{\mathrm{LR}}=1.5\times {10}^{-4}$ and is indicated with the dashed line. (b) Colored symbols: Resulting scaling of ${\gamma }_{\mathrm{bk}}^{\mathrm{dist}}$. Gray background: Scaling of ${\gamma }_{\mathrm{bk}}^{\mathrm{LR}}$, as in Fig. 11. (c) The ratio ${\gamma }_{\mathrm{bk}}^{\mathrm{dist}}/{\gamma }_{\mathrm{bk}}^{\mathrm{LR}}$ varies slowly with $N^{2}P{log}_{10}{\left(N\right)}^{-0.7}$, from ${\gamma }_{\mathrm{bk}}^{\mathrm{dist}}/{\gamma }_{\mathrm{bk}}^{\mathrm{LR}}\approx 0.5$ to ${\gamma }_{\mathrm{bk}}^{\mathrm{dist}}/{\gamma }_{\mathrm{bk}}^{\mathrm{LR}}\approx 1.0$.

Figure 17

(a) Scatter plot of each positive contact breaking strain ${\gamma }_{ij}^{†}$ for 100 synthetic systems drawn (bootstrapped) from the distribution $\rho \left({\gamma }_{ij}^{†}\right)$ for $N=1024, P={10}^{-2}$ (black dots). For each system, ${\gamma }_{*}\equiv min{\gamma }_{ij}^{†}$ is indicated with a red $+$. All values below the $1/N_c$ percentile are indicated with a blue $\circ$. (b) The PDF $\rho \left(x_{ij}\right)$ (black). The distribution of per system minima ($\rho \left({\gamma }_{*}\right)/M$, red dashed) and values below the $1/N_c$ percentile ($\rho \left({\gamma }_{<}\right)/M$, blue dash-dotted) as part of the whole are indicated. (c) Same as (b) but with a linear PDF axis. [(d)–(f)] Same as (a)–(c), with numerical data from the $N=1024, P={10}^{-2}$ ensemble. [(g)–(i)] Same as (a)–(c) with numerical data from the $N=16, P={10}^{-6}$ ensemble.

Figure 18

$Pr({\gamma }_{*}\ge \gamma )$ as a function of $\langle N_{\mathrm{bk}}\rangle Pr({\gamma }_{ij}<\gamma )$ (see text). (a) Solid black: Synthetic data, drawn from $\rho \left({\gamma }_{ij}^{†}\right)$ in the $N=1024, P={10}^{-2}$ ensemble ($\langle N_{\mathrm{bk}}\rangle =1147$). For the same ensemble, data with a single value from a distribution with lower mean (dot-dashed blue) and for systems with an overall per-system scale (dashed purple) are also shown. Dotted red: Synthetic data, from $\rho \left({\gamma }_{ij}^{†}\right)$ in the $N=16, P={10}^{-6}$ ensemble ($\langle N_{\mathrm{bk}}\rangle =16$). The gray line indicates $Pr({\gamma }_{*}\ge \gamma )=exp(-\langle N_{\mathrm{bk}}\rangle Pr({\gamma }_{ij}<\gamma ))$. (b) Data from our simulations. We observe the curves decay slower than exponential, indicating correlations between contacts. Curves from top to bottom: $N=1024, P={10}^{-6}$; $N=16, P={10}^{-6}$; $N=16, P={10}^{-2}$; $N=1024, P={10}^{-2}$; (c) Data from (b) but with all strains rescaled to the mean of strains within one system. This reduces the effect of a per-system scale (dot-dashed red) but does not completely negate it. The behavior for the packing-derived data is unchanged as compared to (b).

Figure 19

Top: Distributions of $u_{\parallel ,ij}$, rescaled by their standard deviation ${\sigma }_{\parallel }$, for ensembles with $N=16$, 256, or 1024 particles at $P={10}^{-6}...{10}^{-1}$. ${\sigma }_{\parallel }$ is indicated in each figure. The distributions develop a sharp kink around 0 for low pressures and become smooth for $P⪆{10}^{-2}$. There is a weak dependence on $N$, with the distribution becoming more peaked for high $N$. Bottom: Same as in the top panel but for $u_{\perp ,ij}$. Here the distributions depend less on $N$ and $P$, although also here the distribution gains weight near 0 for decreasing $P$.

Figure 20

(a) Scaling of the standard deviation ${\sigma }_{\parallel }$ as a function of $N$ and $P$. At low $N^{2}P, {\sigma }_{\parallel }$ is independent of pressure and at high $N^{2}P$ we recover a scaling to ${\sigma }_{\parallel }\sim {\delta }^{0.25}$, consistent with [30]. (b) Same as in (a) but for ${\sigma }_{\perp }$. At low $N^{2}P, {\sigma }_{\parallel }$ is independent of pressure. At high pressure we find a scaling ${\sigma }_{\parallel }\sim {\delta }^{-0.2...-0.15}$, somewhat slower than the ${\delta }^{-0.25}$ found in Ref. [30].

Figure 21

(a) Data rescaled as in Schreck et al. [12] [Eq. (D7)]. Black lines indicate power laws with exponent 1 and 0.75. (b) The residuals $F\left(x\right)/x^{1.0}$ [(dark) blue] and $F\left(x\right)/x^{0.75}$ [(light) red] do not have a plateau, indicating these power laws do not well describe the data. [(c) and (d)] Same as in (b) but with data rescaled as in Wyart [58] [Eq. (D10)], i.e., with $q-r=0.15$. We have chosen $r=1.8$, as in our best collapse.