Unjamming of the Bump model

(1)

MSc Physics

Theory Track

Unjamming of the Bump model

by

Jelle Conijn

11052821

August 29, 2020

60 ec

Supervisor/Examiner:

Examiner:

(2)

Abstract

Jammed soft sphere packings are known to exhibit universal scaling laws of their elastic moduli in the athermal limit when being unjammed. Canonical models of unjamming pairwise poten-tials have a cut-off point to simulate the size of the particles. Therefore derivatives of the pairwise potentials with respect to interparticle distance are discontinuous at the cut-off point. The bump model potential and all its derivatives vanish smoothly at the cut-off point. The form of the equa-tions of the elastic moduli suggest a dependency on the derivatives and perhaps different scaling laws compared to those seen in canonical models. In this research it is found that the bulk and shear moduli scale similar to previously found laws. However, the coordination number and packing fraction do not, as they scale as δz ∼ (p/B)0.3_{and δφ ∼}_{pp/B respectively. A concrete}

(3)

1. Introduction

The jamming transition is a phase transition that governs the change of a liquid to a solid state. This is different from other phase transitions as this transition can also be induced by (local) changes in pressure and shear stress, where a solid material can temporarily unjam and flow as a liquid. Accompanied by unjamming is the abrupt loss of shear rigidity at the jamming transition. Materials known to exhibit this behaviour are foams, emulsions, suspensions and granular materials.

The topic of jamming in disordered media is a relatively new area in solid state physics. The inherent disorder found in jammed media does not lend itself well for analytic solutions that rely on periodicity. Mean-field approaches also fall short when looking at the scaling of shear-to-bulk moduli ratio G/B, as they don’t capture the self-organization processes that determine their scal-ing. The introduction of the computer and its ever increasing capacity [1] allows for accurate simulation of the jamming transition and its emergence of various phenomena such as low fre-quency glassy modes [2], diverging length scales [3–5] and scaling laws of elastic moduli [2,6]. Perhaps, the simplest model where these phenomena have been found is a repulsive soft sphere model, where each particle has a potential that determines its repulsiveness to one another. Vari-ous potentials have been studied in previVari-ous works [2,7,8]. However these models have a cut-off introduced in the potential at some length of interparticle distance, giving rise to discontinuities in derivatives to interparticle distance. The form of the equations for the elastic moduli suggest that a model that lacks these discontinuities might show interesting behaviour in the scaling laws of these elastic moduli.

The bump model is such a model, previously studied in two dimensions by S. Kooij and E. Lerner [8]. In this work we aim to study the various behaviours of the 3D bump model as it gets unjammed. While some scaling laws stay universal, other established laws are broken.

(5)

2. Amorphous solids

2.1 What is disorder?

Most of the theoretical understanding in the field of solid state physics rely on systems’ regular structure that infinitely repeats itself. However, not all solids in nature are orderly. Disordered media can take many forms, for example: Granular media, toothpaste, mayonnaise, shaving foam, plastics, glass, etc. These media are jammed. This means that their underlying structure is mechanically stable. Instead of crystallizing, however, they undergo a jamming transition. Which leaves disorder in the particles that make up their structure. The question that arises is how to handle the disorder in such materials. Effective medium theory, a mean-field approach which assumes that local deformations and forces scale similarly as global deformations and stresses, fall short as disordered media have a non-affine response to applied stresses that make them diverge from long range order media [9–11].

Figure 2.1: Taken from [12]. The left figure shows a crystalline structure. The right figure shows a disordered structure. The disordered structure does not allow long range solutions.

2.2 The jamming transition

The jamming transition is a second order phase transition that occurs when a liquid-like system of particles stops moving because its jammed. The particles in these system can be microscopic particles like atoms and molecules, but macroscopic particles like bubbles or sand grains as well. This means the phase transition of glass doesn’t only apply to glass-like materials, but also ap-plies to foams, emulsions, colloidal suspensions, pastes and granular media [13]. Jamming occurs when the temperature and density of a system becomes such that particles are stuck in their po-sition and can only vibrate. Figure2.2ashows a possible phase diagram for a system that can jam. Unique to the jamming phase diagram there is a third dimension signifying structural load (Σ) that influences jamming. For example applying a shear stress greater than the yield stress (γ > γc) on a gel-like substance, will allow it to flow for as long as that stress is applied. But

as soon as the stress is lower than the yield stress (γ < γc), the flowing stops and the system is

(6)

2.2. THE JAMMING TRANSITION AMORPHOUS SOLIDS

characteristics of a real temperature in amorphous solids [2,14–17].

Computer models are commonly used to study the jamming transition. These models can have varying types of particles, such as hard spheres, colloidal shaped particles, soft spheroids, soft ellipsoids and soft spheres. These models can exist in varying dimensionality ( ¯d) and can also have friction or be frictionless. For a model with soft frictionless spherical particles, which the bump model is, jamming occurs when the coordination number (z) of the particles reaches a point larger than zcrit ≡ 2 ¯d: For every dimension a sphere needs to be held in place by at least two

other spheres. The point J , where the jamming transition occurs is easy to define for athermal soft frictionless particles with a repulsive potential. The repulsive nature of the particles assures the particles are only in contact if the pressure is greater than zero (p > 0). For soft particles the amount of deformation is a direct measure of the rigidity in the system, where more deformation means the system is more rigid. In figure2.2b, a model with soft spherical particles shows how these particles will have their original shape when they are closer to the jamming transition.

(a) Taken from [2].

(b)

Figure 2.2:a: This figure shows a possible phase transition diagram for jamming. The diagram has a third dimension denoting stress (Σ). This is due to the effect where a stress can temporarily unjam a system. Low temperature and high density are both inducers of jamming.b: As density increases, the coordination number (z) also increases. At the jamming point J , φ = φc and z = zc. The amount of deformation increases as the system becomes more jammed.

The density or packing fraction (φ) has, similar to the coordination, also a critical value at the jamming point J , where φ = φc. For soft particles without deformation, the random ordering of

the particles gives a value of φc∼ 0.83 in 2D [18] and φc∼ 0.64 in 3D [13].

Various scaling laws have been found between pressure, coordination, density and elastic moduli for which O’Hern et al. [2,19] laid most of the groundwork. These scaling laws grow with a power-law from zero and have the form

f (x) = Axk+ C. (2.1)

As the equation has the form of a power-law, log-log plots are used to transform the exponent k into the slope of a straight line on such a plot. The constant C is in most cases zero except when the jamming phenomenon has a critical value associated with an observable, for example φcand

zcrit. For these variables δφ ≡ φ − φcand δz ≡ z − 2 ¯d are usually plotted. The scaling of these

observables are consistently δz ∼ (δφ)0.50±0.03_{as can be seen in figure}_2.3a_{. It is noteworthy that}

the value of δz fluctuates quite substantially between individual packings [20]. For this reason every point denotes the average of a several number of packings.

(7)

2.3. CALCULATION OF OBSERVABLES AMORPHOUS SOLIDS

(a) Taken from [2]. (b) Taken from [8].

Figure 2.3: a: Shows the scaling of δz with δφ. The scaling follows a universal δz ∼ √φbehaviour. The points are generated using monodisperse and bidisperse packings of 512 soft spheres in three dimensions with various interaction potentials. The bottom points are generated using bidisperse packings of 1024 soft discs in two dimensions. b: Shows the scaling of G/B with p/B. The curve follows a universal scaling of G/B ∼pp/B. The points are generated using packings of 1600 soft disks with various potentials in 2D.

Other observables such as the bulk (B) and shear (G) moduli start from zero and increase as the system becomes more jammed. The scaling of the elastic moduli are closely related to coordina-tion and pressure and follow a GB ∼ δz ∼ pp/B scaling [2,7,8,13,21,22]. This can be seen in

figure2.3b. Substituting p/B for δz gives a linear scaling of φc∼ p/B.

2.3 Calculation of observables

2.3.1 Elastic moduli

(a) Bulk modulus (b) Shear modulus

Figure 2.4: Taken from [23]. Example of macroscopic bulk and shear modulus. The bulk modulus denotes the resistance of a material to a volume change. The shear modulus denotes the resistance of a material to shear force.

Amorphous solids seemingly behave in many ways as a regular crystalline solid. For instance they will also exert a restoring force when applying compression or a shear force. However, un-like crystals which only have an affine response to deformation, amorphous solids have an addi-tional non-affine response. The bulk modulus (B) represents the reaction of a system responding to compressive strain and is defined as

(8)

2.3. CALCULATION OF OBSERVABLES AMORPHOUS SOLIDS B ≡ −Ωdp dΩ ≡ 1 Vd¯2 d2_U dη2. (2.2)

similarly the shear modulus (G) is defined as

G ≡ F l A∆x ≡ 1 V d2_U dγ2. (2.3)

Where p is pressure, Ω is volume, F is force, A is the surface of the object which is getting sheared and ∆x is the amount of shearing (see also fig:2.4aand2.4b). η and γ represent the infinitesimal changes in compression and shear strain respectively. The middle part of the equation denotes the more generally used definition. While the latter part of the equation is more applicable to microscopic systems using the potential energy (U ). Although the more general definition is easy to understand, the microscopic definition is less straightforward. The transformations η and γto be used at a particle at position r are defined in matrix form as

H(η) =   eη 0 0 0 eη 0 0 0 eη   and H(γ) =   1 γ 0 0 1 0 0 0 1  . (2.4)

They transform the coordinates of a particle such that r −→ H · r. The strain tensor defined as

≡ H T_{H − I} 2 (2.5) then gives (η) = 1 2   e2η− 1 0 0 0 e2η− 1 0 0 0 e2η− 1   and (γ) = 1 2   0 γ 0 γ γ2 0 0 0 0  . (2.6)

The next step is to write the potential energy U in terms of the strain tensor. The Taylor expansion with respect to η and γ gives

δU/Ω ' Uηη + 1 2Uηηη 2 and δU/Ω ' Uγγ + 1 2Uγγγ 2_. (2.7) Where Uα≡ dU_dα and Uαα≡ d 2_U

dα2.The free energy expansion in terms of the strain is then written

as δU/Ω 'X ij Cijij+ 1 2 X ijkl Cijklijkl. (2.8)

Note that the free energy F = K + U is equal to U here since the systems studied are in the athermal limit of T → 0.

(9)

This gives the coefficients:

Uη/Ω = Cxx+ Cyy+ Czz (2.9)

Uηη/Ω = Cxx+ Cyy+ Czz+ Cxxxx+ Cyyyy+ Czzzz+ 2(Cxxyy+ Cxxzz+ Cyyzz) (2.10)

Uγ/Ω = Cxy (2.11)

Uγγ/Ω = Cyy+ Cxyxy. (2.12)

From these expressions it is seen that equation (2.10) corresponds to the bulk modulus and equa-tion (2.12) to the shear modulus. The expression of (2.9) corresponds to pressure as seen in equa-tion (2.13). From this it is concluded that the shear modulus is a combination of a one dimensional pressure term and a combined term.

Beginning with the bulk modulus, the bulk modulus needs to be expressed in terms of a deriva-tive to the strain. Starting from the middle part of equation (2.2), the pressure needs to be written in terms of the energy density. Then the chain rule is used to change the derivative variable to η

p ≡ −dU dΩ = − ∂U ∂η dη dΩ. (2.13)

The volume Ω ≡ Ld¯_{can be expressed in terms of the dilatant strain η as}

Ω(η) = L(η)d¯ = (eηL0)d¯ = eηd¯Ω0. (2.14)

Where Ω0 ≡ Ld0¯ is the volume before any deformation was imposed. Inverting this equation

gives

η = 1 ¯

dlog(Ω/Ω0). (2.15)

Taking the derivatives to η then gives

∂η ∂Ω = 1 Ωd¯ and ∂2η ∂Ω2 = − 1 Ω2_d_¯. (2.16)

(10)

2.3. CALCULATION OF OBSERVABLES AMORPHOUS SOLIDS B = −Ωdp dΩ (2.17) = Ω d dΩ dU dη ∂η ∂Ω ! (2.18) = Ω d 2_U dη2 ∂η ∂Ω 2 +dU dη ∂2_η ∂Ω2 ! (2.19) = Ω d 2_U dη2 1 (Ωd)¯2 − dU dη 1 Ω2_d_¯ ! (2.20) = 1 Ωd¯2 d2_U dη2 + p. (2.21)

As the jamming transition occurs at p → 0 the pressure is omitted from here on. For perfectly spherical particles the potential energy will depend on only the interatomic distance (rij) and

will be a summation of all the potentials from all the particles

U = N X i=0 N X j>i ϕij(rij). (2.22)

Where N is the total amount of particles and ϕijis the potential energy between a particle i and

j, with distance rijbetween them. The derivatives to η are then

∂U ∂η = ∂U ∂r ∂r ∂η = N X i=0 N X j>i ϕ0_ijrij. (2.23) ∂2_U ∂η2 = N X i=0 N X j>i ϕ00_ijr2_ij+ ϕ0_ijrij . (2.24)

Note that ∂r∂η = ras r gets transformed by the compressive strain as r → e η_r_.

The total net force in the system F is defined as

F = − N X i=0 N X j>i ∂ϕij ∂ri = 0. (2.25)

Inserting this force in the total derivative to η and demanding it to be zero because of mechanical equilibrium then gives

dF dη = ∂F ∂η + dxna dη · (2.26) ∂F ∂η = − " ∂F ∂xna #−1 · dxna dη = −M −1_·dxna dη . (2.27)

(11)

Where the second part of (2.26) is the non-affine displacement that occurs to keep the system in equilibrium from the imposed strain. The matrix M = ∂2U

∂xi∂xj is known as the Hessian or

dynam-ical matrix, which will be explained more in depth in section2.4. Finally the second derivative of U to η in equation (2.24) provides the final equation to solve the bulk modulus

B = 1 Ωd¯2 d2_U dη2 = 1 Ωd¯2 ∂2_U ∂η2 + ∂2_U ∂x_na∂η · ∂xna ∂η = 1 Ωd¯2 ∂2_U ∂η2 − ∂F ∂η · M −1_·∂F ∂η . (2.28) The expression ∂F ∂η is solved in section6.1.

A similar equation can be made for the shear modulus [24] G = 1 Ω ∂2_U ∂γ2 − ∂F ∂γ · M −1_·∂F ∂γ . (2.29)

Which does not have a ¯d2in the denominator. The form of both equations suggests a form where one part represents the displacement imposed and another that corresponds to the equilibrium restoring behaviour of the system. Therefore the bulk and shear modulus can be written as

B = Ba− Bna (2.30)

G = Ga− Gna. (2.31)

Where Ba and Gacorrespond to _Ω1∂

2_U

∂η2 and

1 Ω

∂2_U

∂γ2 respectively. The non-affine parts are defined

by Bna≡ ∂F_∂η · M−1· ∂F_∂η and Gna≡ ∂F_∂γ · M−1·∂F_∂γ.

Figure 2.5: Taken from [8]. Scaling ofBna B with

p B. It is interesting to look at what influence the

in-teratomic potential has on the non-affine part of the bulk modulus. A recent result from Kooij and Lerner [8] proved that B/Bna ∼ (r dc/dr)2

re-gardless of dimension, where c(r) ≡ ϕ0/ϕ00r. For an inverse power law (IPL) model, there is no non-affine term due to it being the same as the force in that system up to a factor, which under equilibrium is zero. The factor c is therefore a constant. A more detailed look at different po-tential models can be found in section2.5. Next to the bulk and shear moduli, there is also pressure, which in the athermal limit can conve-niently be written in terms of a sum over the inter particle force and distance [25,26]

p = 1 Ωd¯ N X i=0 N X j>i Fij· rij. (2.32)

Where ¯d is the dimensionality, Ω is the volume of the total system (including empty space) and Fijthe force between a particle pair.

(12)

2.4. THE HESSIAN AMORPHOUS SOLIDS

2.3.2 Packing fraction and coordination number

The packing fraction (φ) is a definition of how much space all the particles take up compared to the total volume. It is a sum over the volume of the particles divided by the total volume. For three dimensions the equation reads

φ = 1 Ω N X i 4π 3 R 3 i. (2.33)

Where Riis the size of a particle i. Finally the coordination number (z) is defined as the mean

amount of particles each particle overlaps with. Overlapping is defined as rij ≤ Ri+ Rj.

2.4 The Hessian

The Hessian is defined as

M = ∂ 2_U ∂xi∂xj ¯ d=3 =               ∂2_U ∂x1∂x1 ∂2_U ∂x1∂y1 ∂2_U ∂x1∂z1 . . . ∂2_U ∂x1∂xN ∂2_U ∂x1∂yN ∂2_U ∂x1∂zN ∂2U ∂y1∂x1 ∂2U ∂y1∂y1 ∂2U ∂y1∂z1 . . . ∂2U ∂y1∂xN ∂2U ∂y1∂yN ∂2U ∂y1∂zN ∂2_U ∂z1∂x1 ∂2_U ∂z1∂y1 ∂2_U ∂z1∂z1 . . . ∂2_U ∂z1∂xN ∂2_U ∂z1∂yN ∂2_U ∂z1∂zN .. . ... ... . .. ... ... ... ∂2U ∂xN∂x1 ∂2U ∂xN∂y1 ∂2U ∂xN∂z1 . . . ∂2U ∂xN∂xN ∂2U ∂xN∂yN ∂2U ∂xN∂zN ∂2_U ∂yN∂x1 ∂2_U ∂yN∂y1 ∂2_U ∂yN∂z1 . . . ∂2_U ∂yN∂xN ∂2_U ∂yN∂yN ∂2_U ∂yN∂zN ∂2_U ∂zN∂x1 ∂2_U ∂zN∂y1 ∂2_U ∂zN∂z1 . . . ∂2_U ∂zN∂xN ∂2_U ∂zN∂yN ∂2_U ∂zN∂zN               . (2.34)

Where xiand xjare derivatives to a dimension of particle i and j. This makes the Hessian a rank

2 tensor with a size of (N ∗ ¯d) × (N ∗ ¯d). The eigenvalues of this matrix (ω2_{) are the square of}

vibrational frequencies that play an important role in many mechanical [27], transport [28], and thermodynamic [29] properties of amorphous solids [30]

M · Ψ = ω2Ψ. (2.35)

Using radially symmetric pairwise potentials as introduced in (2.22) it is then possible to write the Hessian as Mkl= N X i=0 N X j>i (δjk− δik)(δjl− δil) " ϕ00_ij r2 ij −ϕ 0 ij r3 ij xijxij ! +ϕ 0 ij rij I # . (2.36)

Where xij is an arbitrary vector and I the identity matrix. How this is derived can be found in

(13)

2.4. THE HESSIAN AMORPHOUS SOLIDS

The distribution of the vibrational frequencies (modes), the density of states (DOS), describes how many states are available for occupancy at a certain energy level. Higher frequencies cor-respond to higher energies. The Debye model for elastic solids predicts a scaling of the DOS of D(ω) ∼ 1

ωd¯ D

ωd−1¯ _{, where ω}

D is the Debye frequency [31]. The Debye frequency is defined via

RωD 0 D(ω)dω = 1, giving ωD= 18π2(N/V ) 2c−3s + c−3_l !1/3 (2.37)

in three dimensions, where V is the volume. The terms in the denominator are the shear (cs) and

longitudinal (cl) sound speeds and are given by

cs≡ s G ρ and cl≡ s B + (2 −2_d_¯)G ρ . (2.38)

Where ρ ≡ mN/Ld¯_{is the mass density, with m the mass of a particle.}

The modes predicted by the Debye model are called Goldstone modes and are generally well understood. For amorphous solids there exists an additional scaling at the lower frequency part of the spectrum of D(ω) ∼ ω4_{. This additional scaling is independent of dimensionality or}

prepa-ration protocol [8,32–35] and was predicted as far back as 1991 [36–38]. These so called soft glassy modes are modes that stem from the disordered nature of glass. The lowest Goldstone mode is estimated to be

ωc≡ 2π

pG/ρ

L , (2.39)

where µ is shear modulus.

By making the system size L smaller, the lowest Goldstone mode possible will then be pushed up, the glassy modes will get separated from the other modes and the w4_{scaling will be exposed [39].}

The glassy modes are also different than regular phonons in that they are localized instead of spanning the entire system as can be seen in fig2.6d. Each of the dimension also has a translation mode that come from breaking the translational symmetry in any condensed matter system as shown in figures2.6a,2.6band2.6c.

(14)

2.5. POTENTIAL MODELS AMORPHOUS SOLIDS

(a) (b) (c) (d)

Figure 2.6: a,b&c: Three translational modes corresponding to the three dimensions of the system. The vector fields are orthogonal to each other. Note that forc, the vector field points outwards. d: Localized mode, low frequency vibrations do not span across the entire structure, but are instead localized within the system.

To quantify the degree of localization in a single number there is the participation ratio

e ≡ 1

NP

i(zi· zi)2

. (2.40)

Where ziare the ¯d-dimensional Cartesian components of z, the eigenvectors of the Hessian,

per-taining to the i-th particle. Delocalized vector fields will have e of order unity, where localized, glassy modes feature fields with e ∼ N−1in 3 and 4 dimensions and e ∼ (logN )2_/N_{in 2 [33]. This}

can be seen in figure2.7, where the lower frequency modes feature lower participation ratios.

Figure 2.7: Taken from [33]. Participation ratio’s per frequency for 2D (a), 3D (b) and 4D (c). ω0= 1.0, 1.6, 2.8 for 2D, 3D and 4D, respectively, for visualization purposes. The peaks of the participation ratio correspond to the lowest Goldstone modes.

2.5 Potential models

In molecular simulations the interaction potential of the particles play a large role in their be-haviour. This section aims to explain some of the different potential functions used in amor-phous solid and unjamming literature. The function of the potential is graphed against pairwise distance in figures2.8and2.9and is denoted by a blue line. Present also in these figures there is a red line representing the force profile associated with that potential, which is calculated

(15)

us-2.5. POTENTIAL MODELS AMORPHOUS SOLIDS

ing −Fij = ∂ϕij

∂rij. The green line denotes the second derivative of the potential with respect to

pairwise distance.

Perhaps the most famous potential used in molecular simulations is the Lennard–Jones potential

ϕLJ(r) = 4 " σ r 12 − σ r 6# . (2.41)

Where ε denotes the microscopic units of energy and σ the van der Waals radius, which deter-mines the point where the attractive potential changes in a repulsive potential and r the distance between two particles. While the potential is analytic, meaning all derivatives to pairwise dis-tance are smooth, cutting the potential off at long range will introduce discontinuities in them. For unjamming a challenge with this potential is the attractive force at r = 1.122σ. The transition point J is more easily defined in purely repulsive potentials, where a pressure p > 0 jams the system at T = 0.

For repulsive potentials typically the force exerted by the particles increases as they move closer together. Examples of these potentials are the harmonic (n = 2) and cubic (n = 3) potential

ϕn(rij) = n ε nλn (Ri+ Rj) − rij n , rij≤ Ri+ Rj 0 , rij> Ri+ Rj . (2.42) (a) (b) (c)

Figure 2.8: a: Shows the Lennard–Jones potential. The black dotted line represents the zero line. Force is slightly attractive until the r = 1.122σ point. b: Shows the harmonic model. The force profile is a straight line.c: Shows the cubic model with quadratically increasing force.

With λ denoting microscopic units of length, rijthe distance between the centers of two particles

and Ri and Rj their respective size. The exp model is model with an exponential function as

potential and can be found in [7]

ϕEXP(rij) =

n _ij_(r_ij_{) exp (−r}_ij_/l_ij_{) + a}_ij_r_ij_{+ b}_ij _{, r}_ij≤ rc

0 , rij> rc (2.43)

(16)

2.5. POTENTIAL MODELS AMORPHOUS SOLIDS ϕij =    ε0 _r ij lij −β +P3 k=0c2k _r ij lij 2k ,rjj ljj ≤ xc 0 ,rij lij > xc . (2.44)

Where the size of the particles Ri and Rj is encoded in an interaction length variable lij. The

exp model features an extra interaction strength parameter εijwhich is different for the types of

particles used. For information how they are used see [7]. Both models feature next to the exp term and IPL term additional terms. These terms are chosen in such a way that at rij → rcor

rij/lij → xcthe potential ϕij → 0. The coefficients c2k are such that all derivatives to the pairwise

distance are continuous at rij/lij = xc.

Equations (2.28) and (2.29) suggest the bulk and shear moduli are dependent on second deriva-tive behaviour. Ideally a potential that does not need terms of increasing complexity to keep all the derivatives continuous would be used. The bump model is a potential with an exponential that vanishes at rij= Ri+ Rj, named for its bump like feature in the potential

ϕBUMP(rij) = n _{exp (− 1 −} rij Ri+Rj 2₋₁ ) , rij ≤ Ri+ Rj 0 , rij > Ri+ Rj . (2.45)

An advantage of the model is that its derivatives are continuous for infinite order due to nature of the exponent. The disadvantage is however that the force goes to 0 as r → 0. The potential gradient or repulsive force starts to decrease at r = 0.76, this puts a constraint on the maximum pressures that can be generated by this model. Using the equations (2.32) and (2.33) to calculate the volume at high pressures. The maximum possible pressure with the bump model is p ≈ 0.89. Take note this is in ideal circumstances where each particle will be the perfect distance apart to create the highest amount of pressure (for the calculation see6.2). Which disordered systems inherently do not. Care must also be taken when generating a packing with a prior thermalization phase as a too high temperature runs the risk that particles will be too close together as the system gets frozen. This will cause them to get pushed together on one spot. Another way to subvert this issue is by adding a repulsive term to the potential such that the repulsive force only increases. A potential downside to this approach is ruining the smoothness the model gives.

(a) (b) (c)

Figure 2.9:a: Shows the exponential model, the profile is similar to2.8c.b: Shows the IPL model, different from other models since the potential and force increase to infinity as r → 0. c: Shows the bump model. Since the potential gradient decreases at r = 0.76, the force decreases in magnitude until r = 0. This puts an upper limit to the pressures that can be generated by this model

(17)

3. Soft sphere packing generation

and algorithms

The generation and analysis of soft sphere packings comes arises from classical molecular simu-lations. Similarly as in hydrodynamics, the chaotic nature of amorphous solids ensures that little or no analytic solutions exists for many systems. Because of this limitation, the system sizes that can be analyzed scale directly with computer power available. To ensure that computer power is utilized to its full potential smart algorithms are used to study amorphous solids.

3.1 Leapfrog algorithm

A common algorithm that is well known by physicists that simulate molecular systems is the leapfrog algorithm. This algorithm discretizes the equation of motion with small timesteps (∆t)

t → t + ∆t. (3.1)

While the position (ri) and velocity (vi) of a particle i are governed by

ri(t) = xi(t − ∆t) + ∆t · vi (3.2) vi(t + ∆t/2) = vi(t − ∆t/2) + ∆t · Fi mi . (3.3) Where Fiis calculated by Fi= − ∂U ∂ri . (3.4)

Note that the velocity is calculated at half time steps so the average of two timesteps is taken. This is known as the midpoint method or modified Euler method [40]. The miis the mass of the

particle, which is equal to 1 in our simulations, so henceforth it is omitted.

The method has a global error associated with it of order O(∆t2₎_{. This means that for smaller}

timesteps the error decreases quadratically. Which can be seen in figure3.1. The figure is gener-ated by starting a system with particles with a bump potential and running it with a non zero tem-perature without thermostating. After a thermalization phase, energy measurements are taken. As it is a closed system no heat or energy should go in or out. However due to the discretization of time some fluctuations occur in energy. This is shown by taking the standard deviation of the measurements.

(18)

3.2. OTHER METHODS GENERATION AND ALGORITHMS

Figure 3.1: Energy and standard error at different ∆t. Between each order of magnitude there are 100 points separated linearly from each other. Each point has 63 measurements of the energy that were taken in such a way that fluctuations from initial conditions have no effect. The blue points are generated by taking the standard deviation of the 63 measurements. The blue points follow a quadratic scaling consistent with Euler method simulations.

3.2 Other methods

Leapfrog simulations work well for soft spheres, since the particles can exchange their momen-tum over enough timesteps such that the discretization has a minimal effect. However for some models, such as the hard sphere models the collisions and exchanges are instant. For this type of models other algorithms are used, such as the event driven algorithm [41]. Here instead of leapfrogging forward through time, the velocity and size of the particles is used to determine the next collision at time t + ∆τ . The collision is then executed and the whole system is moved forward with time ∆τ . After this, the next collision is calculated again and so the system moves forward.

3.3 Cell subdivision

Calculating when the next collision is going to happen or calculating the force on each particle with equation (3.4) can take a long time with large number of particles. It is however highly unlikely that particles at long distances will meet or have any force interaction at all. A method for improving the speed of the simulation is to divide the box with the particles into smaller boxes. If the size of the smaller boxes is chosen in the right way. The only interactions that need to be checked are those in the smaller boxes and all the smaller boxes next to that box. For example if 1000 particles are divided into 50 small boxes, each box will on average have 20 particles in them. In three dimensions every box will have 14 neighbour boxes. This means that instead of comparing 999 particle pairs, now only 15 · 20 − 1 = 299 pairs need to be checked. This concept can now be taken further by keeping track of neighbouring particles. Where being a neighbour is defined by

(19)

3.4. SOFT SPHERE GENERATION GENERATION AND ALGORITHMS

neighbour < 2Rmax+ ∆R. (3.5)

Where Rmaxis the maximum size of a particle in the system and ∆R is a constant that tunes the

neighbour inclusion ratio. The neighbours are tracked in a list so that for every particles only close neighbours are checked. In practice this brings the amount of particle pairs that need to be compared each timestep down to less than 64.

However, care must be taken that a particle does not enter the potential before it is included as a neighbour. Before that happens the neighbour list has to be updated. Updating the neighbourlist is a somewhat costly operation so ideally it is done as little as possible. An update is required every time a distance could be covered that is equal to ∆R. This is done with equation

2 · v_max· n · ∆t > ∆R. (3.6)

Where vmaxis the maximum velocity recorded for n timesteps. The constant 2 is included since

two particles could move towards each other with vmax. Careful tuning of the ∆R constant gives

control over the update frequency. However if this parameter is set too large the amount of neighbours will be to large and the benefits of this method are lost.

3.4 Soft sphere generation

Soft sphere packings are generated computationally. The packing size that can be generated has a hard limit because of the available computer hardware. Periodic boundary conditions are used to make the space particles reside in repeat itself in every dimension. To generate a packing that is sufficiently amorphous, a temperature is introduced that liquefies the system. To prevent crystallization, the particles will be bidisperse, meaning two different shapes. In the three dimensional bump model particle sizes of Rsmall _{= 0.5}_{and R}large _{= 0.7}_{were used, which}

have been found to sufficiently prevent crystallization.

3.4.1 Thermalization phase

The first step in the making of a soft sphere packing is to set up the particles in a crystal-like structure, equally spaced from each other as seen in figure3.2a. Because the particles overlap and half the particles are a different size, there will already be some forces present in the system. When the system gets moved forward in time, this will result in a chaotic movement of the particles otherwise known as heat. To manage the magnitude of the movement a Berendsen thermostat [42] is applied according to equation (3.7) and (3.8). A snapshot of the liquid phase can be seen in figure3.2b. For the bump model T was set at 0.01 and was run for 2000 seconds

ξT = s 1 + ∆t ∗ Ttarget Tcurrent − 1 /τT (3.7) v → ξT · v. (3.8)

(20)

3.4. SOFT SPHERE GENERATION GENERATION AND ALGORITHMS

(a) (b)

Figure 3.2:a: Shows the starting state of the system. Particles of size 0.5 (blue) and 0.7 (red) are randomly spread in a grid-like structure. The different sizes frustrate the system preventing crystallization.b: Shows the liquid phase of the two dimensional bump model, particles bump in each other and flow through the system.

3.4.2 Jamming the system

To induce the jamming transition the temperature is abruptly reduced to zero. An easy way to accomplish this is to set all particle velocities to zero. However this state is not yet in mechan-ical equilibrium as some particles will start moving again if they are not kept still. To achieve equilibrium the fast inertial relaxation engine (FIRE) [43] algorithm is used. This algorithm re-lies on inertia to achieve fast minimization and is generally much faster than other methods like conjugate gradient (CG) and the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithms for molecular simulations.

3.4.3 Approaching the jamming point

The initial soft sphere packing is formed at high pressure. The goal now is to approach the jam-ming point from the jammed phase, unjamjam-ming the system. As pressure is a control parameter for repulsive spheres in jamming theory. A target pressure can be set and is then achieved by the use of a Berendsen barostat [42]. Energy is continually minimized until an equilibrium is formed. The pressure is given by

p ≡ 1 V ∗d¯ N X i 2φijr2(Ri+ Rj)2 (r2_{− (R} i+ Rj)2)2 . (3.9)

In effect this means that the size of the space the particles reside in grows slightly larger as p → 0. The first pressure is 1 · 10−1and is lowered discretely to 4 · 10−8in 5 steps per magnitude divided equally on a logarithmic scale. The position of the particles are saved, creating a soft sphere packing that can be analyzed.

(21)

3.5. CALCULATING OBSERVABLES GENERATION AND ALGORITHMS

3.5 Calculating observables

Figure 3.3: Example of a Rattler. Observables are calculated from the soft sphere packings

that were generated at different measures of unjamming. To calculate the packing fraction equation (2.33) is used. For z every pairwise distance smaller than Ri+ Rjcounts

as overlapping. A special case is the occurrence of ’rat-tlers’. Rattlers are particles that are enclosed within other particles, but do not contribute to the elasticity of the sys-tem. To eliminate them a combination of two methods is used. First the coordination is counted for every particle, if a particle has less than the critical coordination z < 2 ¯d, it is removed from the system. Other particles might be-come rattlers because of their removal, so this method is repeated until no particles are left with z < 2 ¯d. The second method is to compare the mean pairwise force of a particle with the mean pairwise force of the system. If the force is 5 magnitudes smaller, it is considered a rattler and is elim-inated. This is only done once as the removal of a particle

with such low forces will only marginally change the mean pairwise force of other particles. The bulk and shear modulus require solving of ∂F∂ · M

−1_· ∂F

∂. However, the full calculation of

M−1_{would be impossible within a reasonable time. Instead a is defined as}

a ≡ M−1·∂F

∂ (3.10)

M · a = ∂F

∂ (3.11)

and solved for using conjugate gradient descent. Then it is a matter of contracting ∂F

∂ · a for the

final scalar value. For the DOS and participation ratio the eigenvalues and eigenvectors of M are needed. The Python3 Scipy library is used to calculate these, which relies on the LAPACK library for numerical linear algebra. For large matrices there is also a sparse variant of the Scipy library, which can calculate a subset of eigenvalues and eigenvectors.

(22)

4. 3D Bump model scaling laws

Features of the three dimensional bump model were measured for systems with 1000 particles and 8000 particles. Of the 1000 particle systems 6400 packings were generated on each pressure, while for the 8000 particle systems 800 packings were generated. The mean of all packings is taken at each point in the graph for each pressure unless otherwise specified. The points are displayed using a triangle which is ankered in the middle of the horizontal side. This is done as points for 1000 and 8000 particles largely overlap (for example figure4.4). For the measurement of the bulk and shear modulus the full equation is used specified in as (2.30) and (2.31). The amount of pressures measured is 33, which are logarithmically spaced from the highest pressure of p = 0.1 until the lowest pressure of p = 4.0 · 10−8_{. Any lines drawn through the points are}

done with Scipy’s curve fit algorithm, which uses a least squares algorithm to determine the line parameters.

4.1 Macroscopic Observables

4.1.1 Bulk and shear modulus

(a) (b)

Figure 4.1:a: Shows in essence G scaled with p but made unit-less with the division by B. The slope is both 0.5for 1000 particle and 8000 particle systems.b: Shows the scaling of non-affine term of the bulk modulus with the pressure. This scaling agrees with the scaling laws found in [8].

(23)

4.1. MACROSCOPIC OBSERVABLES 3D BUMP MODEL SCALING LAWS

The scaling seen in figure4.1aof G/B with p/B seems to be universal of soft sphere packings [2,7,8,13,21,22] and also holds here. Similar to that the scaling seen in figure4.1bis reminiscent of earlier results found in 2019 [8]. Where it was found that this scaling is dependent on the ratio of first and second derivatives of the potential to pairwise distance. As these potentials are the same in varying dimensionality, the expectation is that the dimensionality is not a factor and here it shows that as well.

4.1.2 Critical packing fraction

(a) (b)

Figure 4.2: Points from the packing fraction are plotted on a line and fitted by a power law plus constant. This data is then used to generate figureb. The packing fraction is determined to be about 0.65, which is close to RCP φ ∼ 0.64. The scaling of φ as seen inbdeviates from earlier results for soft sphere packings where the slope was 1.0 [13].

The critical packing fraction is determined by fitting a power law plus a constant to the points. Subtracting the found constant makes the curved line a straight one in the log-log plot. The constant that was found is equal to φc ≈ 0.65. It shows φcconforms to the expectation of random

(24)

4.1. MACROSCOPIC OBSERVABLES 3D BUMP MODEL SCALING LAWS

4.1.3 Coordination number

(a) (b)

Figure 4.3:a: Shows the coordination number z plotted against _Bp. The scaling here deviates from earlier found results. The slope of the power law is ∼ 0.3, while previously found values for soft sphere packings are ∼ 0.5 [13].b: Shows the scaling of δz with δφ. It is shown for easy comparison with [2].

Quite interestingly the scaling of δφ ∼pp/B and δz ∼ (p/B)0.3_{is different from previous}

estab-lished work, where they are δφ ∼ _Bp and δz ∼pp/B. A hypothesis for this behaviour might be the force distribution of the packings [44].

(25)

4.2. MICROSCOPIC OBSERVABLES 3D BUMP MODEL SCALING LAWS

4.2 Microscopic observables

Figure 4.4: Eight P/B are selected at equal logarithmic spacing for calculations that are resource intensive.

To calculate the DOS and participation ratio, the Hessian matrix needs to be solved. Since this is a costly operation, only the selected p/B from figure4.4were used. These selections of packings were made such that p/B can be used as a control parameter. The selections are equally spaced in the log-log plot.

(26)

4.2.1 Density of states

(a) (b)

Figure 4.5: Shows the DOS for 1000 and 8000 particles. The different values for p/B correspond with the red encircled points in figure4.4. The value of ωcis calculated using equation (4.1). There are increased low frequency modes with lower p/B and thus lower jamming. This is also seen in [2,21]. The 1000 particle D(ω)is calculated from 6400 soft sphere packings per p/B with each 3000 modes. The 8000 particle D(ω) is calculated from 800 soft sphere packings per p/B with each 24000 modes.

The density of states shows a ω4_{scaling for the highest pressure before transitioning into normal}

Debye behaviour. The modes ω are divided by the frequency of the first Debye phonon ωc, which

is calculate by

ωc=

cs

a0

(4.1)

Where csis the speed of shear waves (see equation (2.38) and a0≡ (V /N )1/3for three dimensions.

At higher pressures more low frequency modes appear. As the system decreases in stiffness more floppy modes become available. This can also be seen in previously studied models with soft sphere packings [2,21].

(27)

4.2.2 Participation ratio

(a) (b)

Figure 4.6: Participation ratio per frequency for 1000 and 8000 particles. Lower frequency modes are more localized and thus have a lower participation ratio. The 1000 particle lines are calculated from 6400 soft sphere packings per p/B with 3000 modes each. The 8000 particle lines are calculated from 200 soft sphere packings per p/B with 24000 modes each.

The participation ratio decreases for low frequency modes. Thus it can be seen that the lower frequency modes are more localized than the higher ones. This makes sense logically, since the frequency of ωcis the longest standing wave that is available to the system.

(28)

5. Discussion and Outlook

In this thesis the bump model was studied in three dimensions. Packings of 1000 and 8000 par-ticles were made at 33 different pressures logarithmically spaced from one another. Starting at p = 0.1and decreasing to p = 4.0 · 10−8. Lower pressures were attempted, but the computation time increased exponentially for each decreasing step in pressure. Thus making it impractical to generate enough packings for lower pressures. Extrapolation to lower pressure seems possible as the straight lines in the plots of4.1a,4.1b,4.2b,4.3aand4.3bappear to have no curvature in them.

In conclusion, the 3D bump model, when looking at the elastic moduli, behaves similarly as previously studied models. The scaling of shear modulus with pressure follows the universal scaling of G ∼ √p. The non-affine term of the bulk modulus follows a scaling of Bna ∼ p,

identical as the scaling found in the paper from Kooij and Lerner in 2019 [8]. While the vibrational modes were not looked at in detail, the appearance of floppy modes closer to unjamming and the appearance of localised modes for modes with a frequency ω < ωcare clearly shown.

The 3D bump model stands out with its out of the norm scaling of the coordination number δz and packing fraction δφ. While in previous works these observables scaled with δz ∼ pp/B and δφ ∼ p/B, they scale in the bump model as δz ∼ (p/B)0.3 _{and δφ ∼} _{pp/B. A possible}

explanation for this is the force distribution of the bump model, which tapers off in such a smooth way, that particles stay connected and densely packed even if the pressure drops. Future work might explore the origin of these scalings.

5.1 Acknowledgements

I would like to thank my supervisor, Edan Lerner, for always showing great enthusiasm and supervision during my project. Also thanks to Corrado Rainone, David Richard, Geert Kapteijns and Karina Gonzalez Lopez for being a wonderful group. Furthermore thanks to SURFsara for their support in providing the LISA computer cluster. Additionally I would like to thank my theoretical physics friends, Ferran Faura Iglesias, Jasper Kager, Stan de Lange and Hoang Vu. Lastly I would like to thank Froukje Vroom for the support provided during the final stages of my project.

(29)

6. Appendix

6.1 Hessian derivation

Assuming three dimensions (x, y, z) the pairwise distance rijcan be defined as

rij = |ri− rj| =

q

(xj− xi)2+ (yj− yi)2+ (zj− zi)2. (6.1)

A derivative to dimension x and particle k would then take the form of ∂rij ∂xk = (δjk− δik) xk rij (6.2)

with the difference vector

∂xij

∂xl

= (δjk− δik)I. (6.3)

Defining ϕij(rij)to be a radially symmetric function only dependent on rij. Taking a derivative

to xkthen gives ∂ϕij ∂xk = ϕ0ij ∂rij ∂xk = (δjk− δik) ϕ0_ij rij xij. (6.4) Where ϕ0ij ≡ ∂ϕij ∂rij.

Going to second order gives

∂2_ϕ ij ∂xk∂xl = (δjk− δik) xij rij ∂ϕ0_ij ∂xij + ϕ0_ijxij ∂ 1 rij ∂xij +ϕ 0 ij rij ∂xij ∂xij ! (6.5) = (δjl− δil)(δjk− δik) " ϕ00_ij r2 ij −ϕ 0 ij r3 ij # xijxij+ ϕ0_ij rij I ! . (6.6)

The next calculation is the derivative to xkand the displacement field as in equation (2.28).

(30)

6.1. HESSIAN DERIVATION APPENDIX ∂2_ϕ ij ∂η∂xk = (δjk− δik)(ϕ 00 ijxij+ ϕ 0 ijxij) = (δjk− δik)φ 00 ijxij. (6.7)

Note how the second term is the force term from equation (2.25) and is thus 0.

6.1.1 Filling in the bump model potential

The bump model has a radially symmetric function as

ϕij = ( exp " _r ij Ri+Rj 2 − 1−1 # , rij ≤ Ri+ Rj 0 , rij ≥ Ri+ Rj . (6.8)

Defining Rij ≡ Ri+ Rj. The derivatives to rijthen become

∂ϕij ∂rij = −2rij r2 ij− R2ij 2R 2 ijϕij (6.9) and ∂2ϕij ∂r2 ij = −2 r2 ij− R2ij 2 + 8r2 ij r2 ij− R2ij 3 + 4R2 ijrij2 r2 ij− R2ij 4 ! R2ijϕij. (6.10)

The first part of equation (6.6) then reads

" ϕ00_ij r2 ij −ϕ 0 ij r3 ij # = " −2 r2 ij rij2 − R2ij 2+ 8 r2 ij− Rij2 3 + 4R2 ij r2 ij− R2ij 4 + 2 r2 ij r2ij− R2ij 2 # R2_ijϕij (6.11) = " 8 r2 ij− Rij2 3 + 4R2 ij r2 ij− R2ij 4 # R2ijϕij (6.12) = " 8(r2_ij− R2 ij) + 4R 2 ij r2 ij− R2ij 4 # R2ijϕij (6.13) = " 4(2r2_ij− R2 ij) r2 ij− R 2 ij 4 # R2_ijϕij. (6.14)

(31)

6.2. HIGHEST PRESSURE CALCULATION APPENDIX ∂2_ϕ ij ∂xk∂xl = " 4(2r2 ij− R2ij) r2 ij− R 2 ij 4 # xkxl+ 2 r2 ij− R 2 ij 2I ! R2_ijϕij. (6.15)

6.2 Highest pressure calculation

The highest pressure is calculated using equations (2.32) and (2.33). Where Ω is calculated from the packing fraction at a certain pressure

p = 4π 3

R3

φd¯· Fij· rij· z(δφ). (6.16)

Where R3_{= 0.234}_{is the mean of the cubed size of the particles in the system and F}

ij· rij = 1.16

the mean force times distance per particle. For this calculation the values Fij and rij are chosen

in such a way the maximum value is achieved. Figure4.3b was used to calculate z(δφ) and figure 6.1 was used to estimate the packing fraction for higher pressures using the highest 4 pressure points. An initial guess was made for the packing fraction at higher pressures, which was used to calculate the pressure. Then using the straight line the packing fraction and pressure are iteratively changed until convergence. This gave a final maximum pressure of p ≈ 0.89 and packing fraction φ ≈ 0.95.

(32)

Bibliography

[1] G. E. Moore, “Cramming more components onto integrated circuits, Reprinted from Elec-tronics, volume 38, number 8, April 19, 1965, pp.114 ff,” Solid-State Circuits Newsletter, IEEE, vol. 11, pp. 33–35, Oct 1965.

[2] C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, “Jamming at zero temperature and zero applied stress: The epitome of disorder,” Phys. Rev. E, vol. 68, p. 011306, Jul 2003.

[3] W. G. Ellenbroek, E. Somfai, M. van Hecke, and W. van Saarloos, “Critical Scaling in Lin-ear Response of Frictionless Granular Packings nLin-ear Jamming,” Phys. Rev. Lett., vol. 97, p. 258001, Dec 2006.

[4] E. Lerner, E. DeGiuli, G. Düring, and M. Wyart, “Breakdown of continuum elasticity in amorphous solids,” Soft Matter, vol. 10, pp. 5085–5092, Jun 2014.

[5] E. Lerner, “Quasilocalized states of self stress in packing-derived networks,” Eur. Phys. J. E, vol. 41, pp. 93–9, Aug 2018.

[6] W. G. Ellenbroek, Z. Zeravcic, W. van Saarloos, and M. van Hecke, “Non-affine response: Jammed packings vs. spring networks,” EPL, vol. 87, p. 34004, Aug 2009.

[7] S. Kooij and E. Lerner, “Unjamming in models with analytic pairwise potentials,” Phys. Rev. E, vol. 95, p. 062141, Jun 2017.

[8] S. Kooij and E. Lerner, “Characterizing nonaffinity upon decompression of soft-sphere pack-ings,” Physical Review E, vol. 100, Oct 2019.

[9] H. A. Makse, N. Gland, D. L. Johnson, and L. M. Schwartz, “Why Effective Medium Theory Fails in Granular Materials,” Phys. Rev. Lett., vol. 83, pp. 5070–5073, Dec 1999.

[10] D. L. Johnson, H. A. Makse, N. Gland, and L. Schwartz, “Nonlinear elasticity of granular media,” Physica B, vol. 279, pp. 134–138, Apr 2000.

[11] H. A. Makse, N. Gland, D. L. Johnson, and L. Schwartz, “Granular packings: Nonlinear elasticity, sound propagation, and collective relaxation dynamics,” Phys. Rev. E, vol. 70, p. 061302, Dec 2004.

[12] Contributors to Wikimedia projects, “Crystal - Wikipedia,” Aug 2020. [Online; accessed 27. Aug. 2020].

[13] M. van Hecke, “Jamming of soft particles: geometry, mechanics, scaling and isostaticity,” J. Phys.: Condens. Matter, vol. 22, p. 033101, Dec 2009.

[14] L. F. Cugliandolo, J. Kurchan, and L. Peliti, “Energy flow, partial equilibration, and effective temperatures in systems with slow dynamics,” Phys. Rev. E, vol. 55, pp. 3898–3914, Apr 1997.

(33)

BIBLIOGRAPHY BIBLIOGRAPHY

[15] L. Berthier and J.-L. Barrat, “Nonequilibrium dynamics and fluctuation-dissipation relation in a sheared fluid,” J. Chem. Phys., vol. 116, pp. 6228–6242, Apr 2002.

[16] H. A. Makse, J. Kurchan, H. A. Makse, and J. Kurchan, “Testing the Thermodynamic Ap-proach to Granular Matter With a Numerical Model of a Decisive Experiment,” Nature, vol. 415, pp. 614–617, Feb 2002.

[17] I. K. Ono, C. S. O’Hern, D. J. Durian, S. A. Langer, A. J. Liu, and S. R. Nagel, “Effective Temperatures of a Driven System Near Jamming,” Phys. Rev. Lett., vol. 89, p. 095703, Aug 2002.

[18] N. Xu, J. Blawzdziewicz, and C. S. O’Hern, “Random close packing revisited: Ways to pack frictionless disks,” Phys. Rev. E, vol. 71, p. 061306, Jun 2005.

[19] C. S. O’Hern, S. A. Langer, A. J. Liu, and S. R. Nagel, “Random Packings of Frictionless Particles,” Phys. Rev. Lett., vol. 88, p. 075507, Jan 2002.

[20] S. Henkes and B. Chakraborty, “Statistical mechanics framework for static granular matter,” Physical Review E, vol. 79, Jun 2009.

[21] M. Wyart, L. E. Silbert, S. R. Nagel, and T. A. Witten, “Effects of compression on the vibra-tional modes of marginally jammed solids,” Physical Review E, vol. 72, Nov 2005.

[22] A. J. Liu and S. R. Nagel, “The Jamming Transition and the Marginally Jammed Solid,” Annu. Rev. Condens. Matter Phys., vol. 1, pp. 347–369, Jul 2010.

[23] Contributors to Wikimedia projects, “Bulk modulus, Shear modulus - Wikipedia,” May 2020. [Online; accessed 4. Jun. 2020].

[24] S. Karmakar, E. Lerner, and I. Procaccia, “Athermal nonlinear elastic constants of amorphous solids,” Physical Review E, vol. 82, Aug 2010.

[25] G. Hummer, N. Gro/nbech-Jensen, and M. Neumann, “Pressure calculation in polar and charged systems using ewald summation: Results for the extended simple point charge model of water,” The Journal of Chemical Physics, vol. 109, p. 2791–2797, Aug 1998.

[26] M. J. Louwerse and E. J. Baerends, “Calculation of pressure in case of periodic boundary conditions,” Chem. Phys. Lett., vol. 421, pp. 138–141, Apr 2006.

[27] M. L. Manning and A. J. Liu, “Vibrational Modes Identify Soft Spots in a Sheared Disordered Packing,” Phys. Rev. Lett., vol. 107, p. 108302, Aug 2011.

[28] R. C. Zeller and R. O. Pohl, “Thermal Conductivity and Specific Heat of Noncrystalline Solids,” Phys. Rev. B, vol. 4, pp. 2029–2041, Sep 1971.

[29] Yu. M. Galperin, V. G. Karpov, and V. I. Kozub, “Localized states in glasses,” Adv. Phys., vol. 38, pp. 669–737, Jan 1989.

[30] L. Gartner, “Micro-structural characterization of glassy solids,” 2019. [Online; accessed 22. Jun. 2020].

[31] C. Kittel, Introduction to Solid State Physics, 8th Edition. Wiley, Nov 2004.

[32] E. Lerner and E. Bouchbinder, “Effect of instantaneous and continuous quenches on the density of vibrational modes in model glasses,” Physical Review E, vol. 96, Aug 2017. [33] G. Kapteijns, E. Bouchbinder, and E. Lerner, “Universal Nonphononic Density of States in

2D, 3D, and 4D Glasses,” Phys. Rev. Lett., vol. 121, p. 055501, Aug 2018.

[34] E. Lerner and E. Bouchbinder, “Frustration-induced internal stresses are responsible for quasilocalized modes in structural glasses,” Physical Review E, vol. 97, Mar 2018.

(34)

BIBLIOGRAPHY BIBLIOGRAPHY

[35] E. Lerner and E. Bouchbinder, “A characteristic energy scale in glasses,” The Journal of Chem-ical Physics, vol. 148, p. 214502, Jun 2018.

[36] U. Buchenau, Y. Galperin, V. Gurevich, and H. Schober, “Anharmonic Potentials and Vibra-tional Localization in Glasses,” Phys. Rev. B, vol. 43, pp. 5039–5045, Feb 1991.

[37] U. Buchenau, Y. Galperin, V. Gurevich, D. Parshin, M. Ramos, and H. Schober, “Interaction of Soft Modes and Sound Waves in Glasses,” Phys. Rev. B, vol. 46, pp. 2798–2808, Aug 1992. [38] V. L. Gurevich, D. A. Parshin, and H. R. Schober, “Anharmonicity, vibrational instability,

and the Boson peak in glasses,” Phys. Rev. B, vol. 67, p. 094203, Mar 2003.

[39] E. Lerner, G. Düring, and E. Bouchbinder, “Statistics and properties of low-frequency vibra-tional modes in structural glasses,” Physical Review Letters, vol. 117, Jul 2016.

[40] Contributors to Wikimedia projects, “Midpoint method - Wikipedia,” May 2020. [Online; accessed 7. Jul. 2020].

[41] D. C. Rapaport, “The Art of Molecular Dynamics Simulation,” Cambridge Core, Apr 2004. [42] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak,

“Molec-ular dynamics with coupling to an external bath,” J. Chem. Phys., vol. 81, pp. 3684–3690, Oct 1984.

[43] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, “Structural relaxation made simple,” Phys. Rev. Lett., vol. 97, p. 170201, Oct 2006.

[44] E. DeGiuli, E. Lerner, C. Brito, and M. Wyart, “Force distribution affects vibrational prop-erties in hard-sphere glasses,” Proceedings of the National Academy of Sciences, vol. 111, p. 17054–17059, Nov 2014.

Unjamming of the Bump model

MSc Physics

Theory Track