Nonaffinity in amorphous solids close to the jamming transition

Nonaffinity is known to be an integral part of the response of amorphous solids. Its role is particularly relevant in particulate systems close to their jamming transition, where it dominates the elastic response. Thus, to determine the elastic properties of amorphous solids it is essential to rationalize the features of their nonaffine response. Via numerical simulations we investigate the relation between the non affine response and the vibrational properties of model amorphous materials. We show that, contrary to previous speculations, modes below the Boson peak are those mostly responsible for the nonaffine response.


Introduction
Crystalline solids are currently very well understood in terms of an ordered array of atoms that repeats itself in space. Linear elasticity and dislocations account for the deformations of the solid when submitted to external stresses. The situation is rather different in amorphous solids [1]. Lacking an ordered structure, the constituents of an amorphous solid perform nonaffine movements when the boundaries are deformed. These nonaffine displacements are not contemplated in linear elasticity theory and lead to the softening of the system.
In a system of particles interacting through repulsive forces, such as granular materials [2,3], nonaffinity might dominate the elastic response. Specifically, when particles are packed at high density they form a solid able to resist shear and compression. As the density or packing fraction is decreased, the nonaffinity increases, which, in turn, affects the elastic moduli. At a critical packing fraction φ c the system ceases to be rigid, as the elastic moduli go to zero. Measuring the distance to the transition in terms of Δφ = φ − φ c , in harmonic solids the bulk modulus is constant but drops to zero discontinuously at the critical point. While the shear modulus vanishes continuously as Δφ 1/2 . The ratio between both moduli defines a length scale that diverges as Δφ vanishes.
Here we investigate the relation between the nonaffine response and the normal modes of the dynamical matrix [4], i.e., the Hessian or matrix of second derivatives of the potential respect to displacements of the particles. In first place, we analyze the features of the density of vibrational states D(ω) which is known to differ from that of crystals. The latter presents a decay D(ω) ∝ ω d−1 where d is the dimensionality of the system. In amorphous solids, instead, D(ω) shows an increased number of low frequency modes e-mail: massimo@ntu.edu.sg forming a plateau known as the boson peak. The powerlaw decay, that we will refer to as Debye scaling, is recovered only below a frequency ω * ∝ Δφ 1/2 [2,3]. Modes below this frequency are supposed to be plane waves, while above, they are heterogeneous. It has been suggested [5] that the modes controlling the nonaffine response of granular systems are those lying in the plateau and with frequency close to ω * .
Here we investigate this speculation in light of recent results showing that the density of states does not follow Debye scaling up to zero frequency. Indeed, it has been recently [6] observed that a very low frequencies the density of states crossovers towards a D(ω) ∝ ω 4 regime below a frequency ω MGM . The corresponding modes are described as quasilocalized in 3D, with localization being independent of frequency for the lower modes. This regime had been predicted for general glasses at high density [7][8][9][10].
After showing that this crossover is also present in a 2D dimensional model of frictionless granular system, we show that the heterogeneity of the elastic regime is due to the modes in the frequency range [ω MGM , ω * ], and not to the modes with ω ≥ ω * as previously speculated.

Methods
We simulate [11] a 50 : 50 mixture of bidisperse soft disks in 2D with size ratio R 1 /R 2 = 1.4. We begin by placing the disks at random positions and minimizing the energy using conjugate gradient method to a target packing fraction. We use this configuration to compute and subsequently analyze the dynamic matrix. Next, we apply a small deformation to the box such that we remain in the linear regime. For simple shear, we apply a strain of magnitude δγ = 10 −5 . For compression, we increase the packing fraction of the system by an amount δφ = 2 · 10 −6 . After the deformation we record the intermediate state of the system, including the unbalanced forces on each particle. Finally, we minimize the energy again using the conjugate gradient method, thus obtaining the final configuration. The particle displacements occurred during this last deformation are their nonaffine displacements.
Our system consists of N = 1000 particles interacting by means a harmonic potential V(r i j ) where r i j is the distance between contacting particles and R i is the radius. We use periodic boundary conditions and, in the case of shear, the Lees-Edwards boundary condition is applied in the transversal direction. The dynamical matrix is diagonalized using the SLEPc library. We run around 4000 simulations for each density, and eigenvectors are computed for around 2000.
The nonaffine displacements can be computed from the normal modes using [5]: Here R i j is the nonaffine part of the relative displacement of neighbouring particles. |δF are the unbalanced forces, |δR(ω) is the normal mode with frequency ω. This expression can be simplified further since modes are very heterogeneous and uncorrelated, then δ R i j (ω) 2 ∼ δ R i (ω) 2 = 1/N, where the last equality is due to normalization. Using again low spatial correlation of the modes, one derives δR(ω)|δF 2 ∼ ω 2 γ 2 (see [5] for more details). With these approximations, one finally has The key assumption of this estimation is that the dominant contribution comes from modes around ω * where the density of states is a constant, while modes at lower frequencies give subdominant contributions [5]. Thus, the above integral is only evaluated for ω > ω * . While this result is in agreement with measurements of the displacement field under shear, and allows to rationalise the scaling of the shear modulus, its underlying assumption has never been explicitly checked. In the remainder of the paper we perform this check, and we will use as a measure of the non-affine response normalized by the deformation applied, i.e., χ = δ R i j 2 /δφ 2 and χ = δ R i j 2 /δγ 2 , for compression and shear respectively.

Results
We compute the theoretical values of the nonaffine displacements using Eq.1 after diagonalizing the dynamical matrix. We compare these values with those obtained the from the actual relative displacement of the particles occurring when minimizing the energy after imposing the shear or the compression deformation. Results are shown in Fig.1. The data show the expected the scaling Δφ −1/2 for shear and we obtain a scaling Δφ −1 for compression. The assumptions used to obtain Eq.1 work very well and turn systematically better at lower densities.
Similar critical-like scalings have been found for related systems. Granular packings under quasi-static shear have been shown to display critical behavior in the dynamical quantities, see e.g. [12][13][14][15][16]. In [17] additional exponents and hyper-scaling relations between them were reported. Related studies of the incremental response of amorphous solids can be found, e.g., in [18].
Length scales have also been discussed recently. The length l ∝ Δφ −1/2 relating bulk and shear moduli was first rationalized in [19] using a variational argument. Subsequent simulations [20] found that the same length scale described the response of granular systems to point forces. However, in [21] a length scale l ∝ Δφ −1/4 was identified studying the response to dipolar forces. A relation between length scales and normal modes was established in [22], where two lengths scaling as Δφ −1/2 and Δφ −1/4 were related to longitudinal and transversal vibrations, respectively. Finally, similar scalings have been related to stability under compression [23] and shear [24].
Before considering how modes with different frequency contribute to the nonaffine response, we analyse the density of states. In Fig.2 we show the low frequency regime D(ω) for all the densities studied. A plateau is apparent for densities smaller than ∼ 0.87. This is the well-known departure of amorphous solids from the Debye scaling. At some lower frequency ω * the linear regime reappears. Finally, below a frequency ω MGM = 2π μ/ρ/L [6], the density of states decreases as ω 4 for all densities. The latter limiting frequency is the minimum Goldstone mode, i.e., the frequency of the largest plane wave that can be fit in the system. Hence, it depends on the system size N, the density ρ and the shear modulus μ. For harmonic solids, given the dependence of the shear modulus on the density, one has ω MGM ∝ Δφ 1/4 . In Fig.2b we rescale the horizontal axis by the minimum Goldstone mode, and the vertical axis to collapse the data below ω MGM . At higher frequencies the data do not collapse because the frequency ω * above which the density of states reaches a plateau scales as Δφ 1/2 .
Let us take a look at the structure of the eigenmodes. In Fig.3 we represent four different eigenmodes for φ = 0.847. In Fig.3a the mode has frequency ω < ω MGM . In Fig.3b,c modes have a frequency in the range ω MGM < ω < ω * , such that Fig.3b is close to the lower frequency, while Fig.3c is close to ω * . Finally, Fig.3d shows a mode in the plateau. The structures in Fig.3a,b look rather similar, with vortices and streams of relatively large displacements mixed with areas of very small displacements. Neither of them seems a neat plane wave and can be more aptly described as quasi-localized or hybrid [25][26][27]. In Fig.3c,d we see that at higher frequencies vortices smear out and large displacements shrink, making the field look more homogeneous and more similar to a plane wave.  In Fig.4 we show the cumulative contribution to the nonaffine parameter, χ = δ R i j 2 /δφ 2 for shear, and χ = δ R i j 2 /δγ 2 , where δ R i j 2 is computed using Eq.1. The modes up to ω MGM give a small contribution that, in addition, is expected to vanish in the thermodynamic limit. Modes in the Debye regime (up to ω * ) give the greatest contribution. This is clearly the steepest part of the sum function, and is followed shortly after by a gentle transition to a plateau. This result is different from the conclusion of the estimation given in [5], where the main contribution to nonaffinity is attributed to modes above ω * .
Converting the sum in Eq.1 into an integral, as it is done in [5], we can gain insight into which modes are contributing more to the nonaffinity. Simplifying the notation, let us write stands for the product of the two averages inside summation. We plot this function in Fig.5 scaled in the horizontal axis by ω MGM (marked by the vertical lines) and in the vertical axis to achieve a collapse of the data. The complicated shape below the minimum Goldstone mode should not concern us since it will not contribute in the thermodynamic limit. Above ω MGM there is a small plateau and afterwards the function grows quadratically with the frequency.
In order to perform the integral we consider the intervals [ω MGM , ω * ] and [ω * , ∞]. Note that to take into account the transition from a plateau to a power law of Π(ω) we need to introduce a frequency ω 2 in the first interval. Since the data appear collapsed in Fig.5 we know that ω 2 ∼ ω MGM ∝ Δφ 1/4 /L. For D (ω), as discussed in reference to Fig. 2, we use the scaling shown in Fig. 2a close to ω MGM . While close to ω * we use the scaling imposed by this frequency, i.e., D (ω) ∝ ω/Δφ 1/2 . Finally, note that due to the normalization of the modes Π(ω) scales like 1/N. Summarising, we split the integral in three parts: Plugging all the scalings in the integral, we have (for shear): This is effectively the correct scaling for shear. For compression Π (ω) is scaled by Δφ 0 instead of Δφ 1/2 , so we obtain immediately R 2 ∼ Δφ −1 .
Summarising, in this work we have carried out a detailed investigation of the contribution of normal modes to the nonaffine displacement. We have shown that modes below ω * cannot be disregarded since they give an important contribution, and set the correct scaling of nonaffinty with density. Modes above ω * increase the value of the nonaffinity but the scaling is unaffected. There remains to investigate the thermodynamic limit. We expect our conclusions to hold since in that case the Debye regime is larger, because ω 2 ∼ ω MGM decreases faster that ω * .
We have seen that the structure of the modes changes slowly a one moves from low to high frequencies. As a consequence, modes below the minimum Goldstone mode and even those in the Debye regime cannot be neatly described as plane waves. It would be interesting to investigate the relationship, if any, between the structure of the modes and their contribution to nonaffinity.