Contact forces in nearly jammed hard disc glasses

(1)

Verslag van Bachelorproject Natuur- en Sterrenkunde

Omvang 15 EC, uitgevoerd tussen 28-04-2015 en 02-07-2015

Contact forces in nearly jammed hard

disc glasses

Author:

Mirte van der Eyden

Student number:

10358307

Institute: Institute for

Theoretical Physics Amsterdam,

FNWI

Supervisor:

Dr. Edan Lerner

Second corrector:

Dr. Tom Kodger

July 2, 2015

(2)

Abstract

Despite the apparent simplicity of hard disc glasses, they form a useful model to investigate information about the contact forces between particles in hard disc glasses. We numerically compute these forces in isostatic soft sphere packings and compare them to those in hard disc glasses, to see that they agree. This was expected because there’s only one solution for the contact forces, unique up to a multiplicative constant, that solves the system of equations that originate from force balance and constant pressure. Furthermore, the average distances hhi between the particles in the glass and the contact forces are proven to relate as hhi ∼ 1

f, which

agrees with theory. The rate at which interacting pairs of particles depart from each other under quasi-static decompression is found to be inversely proportional to the contact force. This means that contacts with a weak force separate in a singular way in the thermodynamic limit. Finally the separation of weak forces into localized and extended forces is not the same as the separation between forces with a low collision frequency and a low average momentum exchange. Further research is needed to see if the separation can be made another way in hard disc glasses. The possibilities for further research are numerous. It would be interesting to study the emergence of new interactions that appear upon quasi-static decompression of isostatic packings.

(3)

1 Balletjes als simpel model voor glas

In heel veel vaste stoffen zitten de atomen vast in een strak geordend rooster. Dit wordt ook wel een kristalstructuur genoemd. Stoffen zoals glas hebben echter helemaal geen structuur, het lijken wel vloeistoffen die op een bepaald moment stil zijn gaan staan. Dit is ook precies wat er gebeurt. Glas, maar ook gel en plastic, is een als vloeistof die zo stroperig is dat de tijd die het duurt om te vervormen langer is dan we ooit zullen meten. Als deeltjes een glas zijn, kan je ze nog verder in elkaar person door deeltjes onder druk te zetten tot ze geen kant meer op kunnen. Dit wordt ook wel de ’jammed state’ genoemd. Hoewel gewone faseovergangen zoals bevriezen en ontdooien helemaal natuurkundig doorgrond zijn, is over deze faseovergangen in amorfe materialen nog erg weinig bekend.

Om meer over deze toestand te weten te komen wordt veel gebruik gemaakt van simulaties op de computer. We hebben een model dat op glass lijkt en daarmee kan je de ’jammed state’ nabootsen. Het model dat wij gebruikt hebben bestaat uit allemaal harde balletjes die als biljartballen tegen elkaar aan kunnen botsen. Je zou zeggen dat dit geen erg goed model is omdat atomen wel wat ingewikkelder zijn dan harde balletjes, maar verrassend genoeg zijn er heel veel dingen die uit zulke systemen kan afleiden. Het mooiste is, je hebt er helemaal geen lastige natuurkunde voor nodig!

Het eerste waar we naar gekeken hebben is wat de krachten zijn tussen alle balletjes. Als de balletjes in de ’jammed state’ zijn en je geeft ze net iets meer ruimte, gaan ze trillen waarbij ze de hele tijd tegen hun buren aan botsen. Bij al deze botsingen geven ze de buren een zetje en krijgen ze een zetje terug. Dit ’zetje’ wordt een stoot genoemd en dit is iets dat de balletjes uitwisselen bij iedere botsing. De kracht die de balletjes op elkaar uitoefenen kan je uitrekenen door de grootte van deze stoot te delen door de totale tijd. Op deze manier zijn alle krachten berekend en dit kan je zien in figuur 1. Hier zijn alle krachten weergegeven als lijntjes en des te dikker het lijntje is, des te sterker de kracht.

Figure 1: Het contact netwerk tussen de balletjes. De dikte van de lijn is evenredig aan de sterkte van de kracht.

Een van de dingen die we ontdekt hebben is dat als je het volume iets vergroot en de balletjes dus meer ruimte geeft, dat de balletjes met de zwakste krachten het snelste van elkaar weg bewegen. Als je hier over nadenkt is dit misschien best vreemd, want intuitief klinkt het logischer dat als twee balletjes heel hard tegen elkaar duwen ze ook snel weg zouden bewegen, net zoals jij snel weg beweegt als je je hard afzet tegen de muur. Dit is dus niet het geval.

(5)

2 Introduction

Amorphous materials are without crystalline structure or long-range order. A few examples are windows, plastic and gel. They have the property that there is no distinct phase transition to the solid phase. The transition that occurs is called the glass transition and there are a still a lot of unsolved questions on this topic. To try and get an understanding in the behaviour of these materials one has to simplify. Therefore we use a model of hard frictionless spheres that collide elastically. Despite the simplicity of the hard sphere model, it allows for the derivation of interesting scaling properties of the solid such as elasticity (the capability of the solid to return to its original shape after deformation), vibrational properties or rearrangements under shear deformations (4). The hard spheres that have a certain free energy can be compressed and make the transition from a liquid into a solid. The rate of this compression determines the structure that will form. A slow quasi-static compression will lead to a crystalline struc-ture, since this structure has less possible configurations and therefore a lower entropy and therefore higher free energy. A fast compression will lead to a glassy state that state can be compressed further, deep into the glassy state, up to the jamming point. At that point the spheres cannot be compressed any further, they are jammed. Near the jamming point, there are a lot op open questions about the exact behaviour of these hard sphere systems.

This report will focus on hard disc glasses near the jamming point. We will take hard frictionless particles as a model for the amorphous material to obtain several physical prop-erties of these systems. Our main focus will lie on the collisions between particles and the information that we can extract about the forces from those collisions. The contact forces are numerically evaluated in 2 dimensional jammed hard disc glasses of 64 particles. To obtain those packings, isostatic soft particle packings are decompressed such that there is a small gap h between the particles to move and collide with their neighbours. The forces in this isostatic packing are proportional to their overlaps. The first task is to show that these isostatic forces are the same as in the hard disc glasses. It is expected that this correspondence holds best when the gap is as small as the numerical calculations allow to make sure that there are as little new contacts as possible. When we deflate the system further, there will be a point where a lot of new interactions occur and the correspondence between the isostatic packing and hard disc glass will break down. The second goal is to find the packing fraction (the fraction of the packing that is occupied) for which this happens.

When the above questions are answered we are sure the system works properly. We will then focus on the relation between the average distance between two particles and their contact force. Intuitively this relation is expected to be inversely proportional. The closer the particles stay together, the more frequently they collide and therefore the bigger the contact force is. In section 4 there will also be a theoretical prediction of this relation. The relation between the average distance and the packing fraction φ, or more precisely δφ = φc− φ

with φc the critical packing fraction, is less intuitive. Because of scaling arguments that are

described in section 7.3, we expect the following relation: dhhi_dφ ∼1

f. The hard-sphere system

will be used to verify this.

(6)

because they are most likely to break and interact with new particles. When an external stress is implied, the small forces have to be understood to be able to predict the rearrangements and describe the new configuration that arises. The distribution and density of these weak forces effects the ability of the solid to flow. It is because of these reasons that it is important to understand the origin of weak contact forces. In order to do that we can just take the packings that have a lot of weak forces in them and use them for our analyses. Edan Lerner, Gustavo D¨uring and Matthieu Wyart in Ref (5) made a distinction between two kinds of low forces: extended and localized forces. When breaking extended forces there will be a spatially extended response in the whole system, while perturbing localized forces only has significant influence on the nearby particles, see figure 2.

Figure 2: The displacement field of the system caused by creation of a floppy mode by breaking a contact while leaving all other contact intact. a) extended forces have a great influence on the whole packing when breaking them, b) localized forces have only influence in the direct neighbourhood. (5)

At first the assumption was made that when external stress is imposed, the contacts carry-ing the small forces will open (5). This leads to the possibility to calculate several microscopic quantities of the packing in terms of an exponent θ. Later it was however found that only the extended forces play that role, since they are strongly coupled to the system. (4) This is the reason that the distinction between the two kinds of forces is of high interest.

A final goal of this research is to also make a distinction between two kinds of forces in packings with frictionless hard spheres. We expect that the information about the response of the system, that is used to separate between localized and external forces, is somehow obtained in the local geometry of the contact. The forces in this system are calculated by dividing the total momentum exchange (impulse) over time. Therefore we can find low

(7)

con-tact forces either where the average impulse is low or where the collision frequency is low. The aim is to show that the contacts with low frequency correspond directly to the extended forces and the forces with low impulse to the localized forces. If that is the case, all that is known about the distribution of localized and extended forces is then immediately also appli-cable to the contact forces in the packings in this research, which has never been shown before.

The report is organized as follows. First the two models we use are described in detail. The next section contains theory about the partition function of hard sphere systems, jammed systems and the contact forces in those systems. Then the research goals and hyptheses will be specified in section 5. Our numerical method to answer our questions will be presented in section 6 and our results are shown in section 7. Finally there will be suggestions for further research and a conclusion.

3 Model

We use two different kinds of models near the jamming point. The first one, the isostatic packing is an athermal system of soft spheres that can overlap and the second one, the hard disc glass consists of hard discs with a certain temperature that collide with each other. In this section both models are described and the quantities of interest and their notations are introduced.

3.1 The isostatic packing

The two-dimensional isostatic packing is a box that consists of 64 discs/particles, made up of two subsets of different sizes that are of ratio 1:1.4. On top of those sizes they have a small random fluctuation on their radii, to make sure no hexagonal configurations of particles from one subset will form (7). The masses m of all particles are 1. The particles are ’soft’, which means they are elastic and the distance between the centers of two particles rijcan be

smaller than the sum of their radii. This means there can be an overlap between the 2D-discs. The forces between those particles are repulsive. They are calculated like the particles are connected by springs and the overlap is the amount that the spring is displaced from its relaxed position. This relaxed position is when there is no overlap. We take the spring-constant k = 1 for all particles, such that the force from particle i acting on j can be calculated by:

~

Fij = ∆~xij, (1)

with ∆~xij≡ (σi+ σj) − rij the overlap, with σithe radius of particle i.

The packing fraction φ is defined as the fraction of the total volume (area in 2D) that is occupied by particles, φ = Vo

V. At the jamming point, this packing fraction has a critical value

φcthat is around 0.83 for 2D systems. The isostatic packing we use is just at the jamming

point. The pressure p in the system is smaller than 10−12so the particles are still overlapping but the contact forces are very low (of the order of p). The coordination number z is the average number of contacts for a particle in the packing. The packing is called isostatic when

(8)

z = 2d, with d the spacial dimension so this case, z ≥ 4. The other requirement that comes from the solidity of the glass is that the system is in mechanical equilibrium. This means that the net force on each of the particles is zero:

∀i Xfijn~ij = 0. (2)

The force vector is here separated in its magnitude fij and its direction ~nij, which is a

unit vector pointing from the center of particle i to the center of j. The boundaries of the box are periodic. The system has no thermal energy in it, it has a temperature of zero. All the energy there is is stored in the potential energy from the interactions.

3.2 The hard disc glass

The two-dimensional hard disc glass has a lot of the same properties as the isostatic packing. It is created by inflating the isostatic packing (section 6) and therefore it’s built of the same 64 particles with the same sizes and masses, but in this case the particles are hard. This means they cannot overlap and the φc is the maximum of the packing fraction. We define

δφ = φc− φ to be the relative packing fraction, compared to the critical one. In this system

the particles have a certain kinetic energy that they use to move around. When the particles move and the only potential there is is a discontinuous one that has a repulsive influence when two particles touch. If the distance between two particles, rij, is exactly equal to the

sum of their radii, they collide elastically. The two particles exchange momentum that can be calculated using the conservation of momentum and energy. A full description of the 2D elastic collision can be found in appendix A.

This system has periodic boundary conditions too. When δφ is small, the glass is close to jamming. Since the glass is in its solid state, the particles still move but their time-averaged position is constant. They will have well-defined neighbours that are defined as the particles with which the particle collides.

Using classical statistical mechanics, a lot of properties of these hard disc systems can be derived. This will be covered in section 4.2.

4 Theory

4.1 Contact forces in the isostatic packing

In the isostatic packing, the forces are in balance. Therefore, knowing all the positions and sizes of the particles, the contact forces can be calculated by solving the set of equations, consisting of N × d linear equations which can be written in bra-ket notation as

ST|f i = 0, (3) defining:

S ≡ ∂rij ∂ ~Rk

(9)

To solve this linear set we need the constraint that the pressure is fixed:

X

<ij>

fijrij = pΩd, (5)

where Ω is the volume and d the dimension. We want to solve the combined equation:

H−1|f i = |0, pΩdi (6) for the forces |f i. To do that it is more convenient to use the symmetric form

HHT|f i = H |0, pΩdi = pΩd |ri , (7) where |ri is the pairwise distance between the centers of two neighbours. The solution for the contact forces is unique upto a multiplicative constant, because if |f i is a solution of 3, than also c |f i.

4.2 Hard sphere physics

4.2.1 Analogy between hard sphere glasses and elastic networks

In this section we describe how we cam make an analogy between a hard sphere glass and a system of elastic springs (2). This analogy will enable us to derive a relation between the contact force f and the average gap hhi between two particles in the hard sphere system. We consider a hard sphere glass where the collision time among neighbour particles is much smaller than the time at which the structure rearranges. Therefore we can define a contact network, since every particle has a constant set of neighbours. The coordination number z is therefore defined in this meta stable state. In this state with the contact network, all possible configurations for which particles in contact do not interpenetrate, which is demanded in a hard sphere system, are equally likely. Those configurations satisfy:

Y

hiji

Θ(|| ~Ri− ~Rj|| − σ) = 1, (8)

where Θ is the Heaviside step function, σ is the sum of the radii of the two particles and the product is made of all contacts hiji. This equation says that the distance between two particles must be equal to the sum of their radii, since they’re exactly touching. Since this equation is 1, we can put in in the partition function without consequences. We can now write the partition function of the system:

Z(p, β) = Z dV Y i Z D ~Ri Y hiji Θ(|| ~Ri− ~Rj|| − σ)e−βpV. (9) Here β ≡ 1

kbT with T the temperature and kbthe Bolzmann constant. In one dimension,

we can make a coordinate transformation using hij= Ri− Rj. Using this, we have a

(10)

Y i D ~Ri ∝ Y ij dhijδ( X ij hij− (V − V0)). (10)

V0 is the volume of the packing when p → inf. Substituting this in equation 9 gives us:

Z = Y

ij

Z

hij≥0

dhije−βhijp. (11)

This tells us that p = 1 βhij =

kbT

hij, in one dimension. However, in higher dimensions the

mapping we used in equation 10 is not one-to-one and linear. Because we look at the hard sphere system in the jamming state, we can make an exception. The system is isostatic and therefore the number of contacts is exactly equal to the number of degrees of freedom in the system. We can thus map the positions of the particles one-to-one to the gaps between the particles in contact. When we are nearly at maximal packing, this is a linear map. Furthermore, since we have an isostatic state, there is mechanical equilibrium. This means the net force on every particle is zero and the work that has to be done to move a particle has to be zero as well. The virtual work theorem dW =P

ijfijdhij− pdV = 0 leads to the

volume constraint δ(P

ijfijhij− p(V − V0)). This volume constraint is the generalized form

of δ(P

ijhij− (V − V0)) from equation ?? for higher dimensions, with fij the contact force

between particles i and j. Using this, we can define a generalized one-to-one linear mapping for an isostatic packing:

Y i D ~Ri ∝ Y ij dhijδ( X ij fijhij− p(V − V0)), (12)

with fij the contact force between particles i and j. Substituting this in the partition

function 9 leads to:

Z = Y

ij

Z

hij≥0

dhije−βfijhij (13)

and therefore we have:

fij =

kbT

hij

. (14)

The force is therefore predicted to be inversely proportional to the average gap between the particles. It is one of the goals of this research to check this relation numerically. The analogy between the hard sphere system with these properties and an elastic system, is that since we know the relation between f and h is inversely proportional, we can compare the system to a network of logarithmic springs. The springs are logarithmic because instead of the harmonic spring, where the potential U ∼ x2, we have a spring were U ∼ log(x).

(11)

4.3 Extended and localized forces in isostatic packings

Most of the contact forces in the isostatic packing are of the order of the pressure, so f_p is of order 1. The distribution of these contact forces varies more than two orders of magni-tude. The small forces play a special role in the dynamics of isostatic packings. As you apply external stress to the system, these forces are less stable and therefore more likely to open. Based on the geometry, a distinction can be made between two kinds of weak forces (5). The distribution of both kind of low forces contains a lot of information about the geometry and structure of the packings and are therefore interesting observables.

The system used in Ref. (5) consists of N frictionless hard particles in a box of volume Ω with rigid walls. Our system has periodic boundary conditions, but the for the theory presented here, this does not matter. The d-dimensional packing is created by reducing the box size of an initial low dense packing. The boundaries apply forces on the particles in contact with it. The position of the ith _{particle is denoted by ~}_R

i and ~Rij = ~Ri− ~Rj is the

distance between a pair of particles.

In a mechanical stable jammed packing all forces are in balance:

∀i, _F~_i₊X

j(i)

fij~nij = 0. (15)

Here the sum is over all particles j(i) that are in contact with i, ~Fiis the force applied by

the wall (if present) and fij is the magnitude of the contact force between particle i and j.

This force is purely repulsive and is directed along ~nij, see figure 3.

Figure 3: The unit vector ˆnij is the direction of the momentum exchange and therefore also the

force that particle i exerts on j during the collision.

Another way of defining mechanical stability is to require that there are no floppy modes. Floppy modes are collective motions of the degrees of freedom of the system for which the distances between objects in contact are fixed. This means all particles that are in contact stay in contact, but they can move collectively. If there would be such a floppy mode, the

(12)

system would rearrange until the floppy mode is gone. Packings of hard particles that are rapidly compressed or poly-disperse are isostatic, which means that the average number of contacts between particles is just sufficient for mechanical stability. This number is called the coordination number z and in an isostatic packing z = zc= 2d. Removing one contact in the

packing leads to the existence of exactly one floppy mode. The impact of this floppy mode is dependent on the force that was in the contact before it was broken.

Figure 4: The breaking of a contact, creating a gap h between the two particles.

The way to generate a floppy mode is by pushing two particles 1 and 2 apart thus breaking the contact (figure 4) α between them, while all the other contacts in the system remain intact. δ ~R_i(α)(s) is the displacement of particle i after opening the contact α by a distance s. The displacements for all particles in the packing form a displacement field, which is uniquely defined and satisfies the following equation:

δ ~Rij(α)(s) · ~nij+ O(s2) = sδα,ij. (16)

In the above equation δα,ij= 1 if the pair hiji is the opened pair (so equal to the pair α)

and is zero for all other pairs. ~nij≡ ~Rij/rij is the unit vector directed from particle i to j.

The equation tells us that for every pair of particles, the change in their distance projected on their contact vector ~nij is zero for all pairs except pair α, for which it is equal to s. The

rest term O(s2) comes from the Taylor expansion for small s.

Using the above, the goal is to find two kinds of floppy modes that correspond to two kind of low forces in the packing. We start with multiplying equation 15 with a displacement field δ ~Ri. This leads to:

∀i, _F~_i_{· δ ~}_R_i_{+ (}X

j(i)

fij~nij)δ ~Ri = 0. (17)

Summing over all particles gives us the virtual work theorem:

X i ~ Fi· δ ~Ri+ X hiji δ ~Rij· ~nijfij = 0, (18)

(13)

where P

<ij> means summing over all contacts. The external forces are related to the

applied pressure p from the rigid boundary by:

X

i

~

Fi· δ ~Ri = −pδΩ. (19)

Substituting equation 16 into equation 18 and using the above gives:

pδΩ(α) = −X

i

~

Fi· δ ~R(α)i = sfα+ O(s2). (20)

To see the relevance of the previous steps, it is convenient to introduce the linear floppy mode: ~ V_i(α) ≡ d ~R (α) i ds |s=0. (21) This is a measure of how much particle i will displace if a contact α is broken with a distance s, evaluated at the very beginning of the breaking (s=0). If we take the derivative of equation 20 with respect to s, we obtain:

fα = −

X

i

~

Fi· ~Vi(α). (22)

Furthermore, it is convenient to introduce a vector field ~Vi∗(α)for each particle i that restricts

the floppy modes to be at the boundary. Because we already know that for particles in the bulk (so with no contact with the boundary) ~Fi is zero, these terms will not participate in

the sum in equation 22. Therefore, with ~Vi∗(α)= ~V (α)

i for all particles i on the boundary and

~

V_i∗(α)= 0 for all particles i in the bulk, we can write:

fα = − X i ~ Fi· ~V ∗(α) i = −||F || · ||V ∗(α) || cos θα, (23) whereP iX~ 2 i ≡ ||X||2 so P iF~i· ~V ∗(α)

= ||F || · ||V∗(α)|| and θα is the angle between the

two vector fields.

What can be concluded from this is that like before, a contact force fα can be small for

two reasons.

First of all the force is low when ||V∗(α)|| is small compared to the average of the force between other contacts β:

||V∗(α)|| h||V∗(β)||iβ, (24)

where hi >β means averaging over all contacts in the packing. The typical value of ~V (α) i

is of order 1 in an isostatic system. This means that if we remove one contact, all the parti-cles in the packing will displace by roughly the same amount. In special cases, the partiparti-cles around the contact are distributed such that the coupling of the local displacements to the rest of the packing is weak. This means that the displacement is only significant in the direct

(14)

neighbourhood of the broken contact. Both modes are shown in figure 2, where in figure a an extended force is broken and in figure b a localized force.

To make a clear distinction it is convenient to quantify these behaviours. Therefore, we introduce the dimensionless quantity:

bα =

||~V∗(α)|| h||V∗(β)_||i β

. (25)

This way, contacts with bα 1 are mechanically isolated and forces applied on the

boundaries of the system have little impact on the contact force fα.

The second way a contact force fα can be small compared to other forces is (see eq. 23)

when cos θα hcos θβiβ. (26) Defining: Wα ≡ − cos θα hcos θβiβ and (27)

ftyp ≡ ||F ||h||V∗(β)||iβhcos θβiβ (28)

then equation 23 can be written as:

fα

ftyp

= Wαbα. (29)

In this above formula, the two types of low forces are easily distinguishable. We expect the distributions of b and W to decay exponentially when their values are larger than typical (b > 1, W > 1).For b 1 and W 1 we assume:

P (b) = bθb_, ₍₃₀₎

P (W ) = WθW_. ₍₃₁₎

This θband θW relate to the geometry and structure of the packings. They are calculated

numerically in ref (5) with the results θb = 0.17 and θW = 0.44, see figure 5 and 6. These

distributions will be used in section 7 to compare with distributions of physical quantities in the hard disc model.

(15)

Figure 5: The distribution of b, defined in equation 25. When b 1 for a certain con-tact, breaking it would have local influence only.

Figure 6: The distribution of W , defined in equation 27. When W 1 for a certain contact, breaking it would have impact on the entire packing.

4.3.1 Numerical methods

In this section we describe the way of generating the jammed packings and calculating the quantities fα, bαand Wαthat are used in Ref (5) for the analysis above. The jammed packing

is created by isotropically compressing a system of hard frictionless spheres, until the jamming state is reached. The forces are calculated using the theory described in section 4.1.

To calculate b, the response to an opening of contact α with gap h in the far field, the displacement ||~Vα

ij|| where V ≡ dRdh is evaluated for all pairs < ij > within a ring with a

radius of L/2 around the perturbed pair and a width of the size of a particle. b is defined as the mean of these ||~Vijα||. Another option is to take into account all pairs in the packing and

calculating the median instead of the mean, because in that way the few particles near the contact α with V ∼ 1 do not affect b.

Finally the last quantity of interest W is calculated using the above:

Wα =

fα

bα

. (32)

4.4 Jammed packings

We’re interested in the behaviour of the system at the jamming point. At this point the particles are locked at a certain position, with only a little space to move around and collide repeatedly with a small number of neighbours. The forces between the colliding neighbours form the contact network of the packing (see figure 10 and 11).

There are various ways to create a jammed packing (3), that result to packings with different density and structure. The options that are used the most will be described below and then our approach will be thoroughly explained in the next paragraphs.

(16)

• Isostatic Compression. A packing with a random packing fraction below the critical packing fraction φc is compressed until the jammed state is reached and the particles

are in mechanical stability. This method is used in the analysis of Ref (5).

• Energy Minimization. Here instead of hard spheres, one uses soft spheres and then looks for the potential minima. The interaction potential can be chosen to be very stiff (like the hard sphere potential).

Our approach is one of energy minimization, we start with a soft sphere model where the energy is minimized every time-step. The three steps that are taken to create the jammed packing are described in section 6.

5 Research questions and predictions

In this section the goals of our research are presented in more detail and the predicted answers of the questions are argued.

5.1 Comparison of forces in the isostatic packing and hard disc

glass

The contact force ~Fijbetween two particles i and j is calculated in two ways:

• In the isostatic packing: ~Fij= −k∆~xij. with ∆~xijthe overlap between particle i and j.

• In the hard disc glass: ~Fij= P

c∆~pc

T , where the sum runs over all the collisions between

particle i and j and ∆~pcis the momentum exchange.

These two ways of calculating the forces should give the same contact network because both systems are in mechanical equilibrium. In the isostatic situation, the particles do not move so all the energy is stored in potential spring-energy. In the hard disc model, the energy is kinetic and the particles move and collide elastically with their neighbours. Both systems are (nearly) jammed. The set of equations given in section 4.1 have a unique solution for the contact forces (up to a multiplicative constant). Therefore the contact networks of the two systems should agree and have a linear relation.

5.2 The variance of δφ

Before launching the hard disc dynamics, the packing fraction φ to which the system is expanded compared to the critical packing fraction φc has to be chosen. In order for the

contact-forces to be the same, the system has to get as few new contacts as possible during the hard disc dynamics. Therefore we expect that δφ has to be very small, so that the particles are vibrating on a very small area and have no room to move around and interact with new neighbours. To test for the optimal gap, we did the dynamics for δφ varying logarithmically from 10−1 down to 10−10 and compared the isostatic forces to the forces in the hard disc dynamics. The sum of squares due to error (SSE) is defined as:

(17)

SSE =

N

X

i=1

wi(yi− ˆyi)2, (33)

with N the number of data points, yi the data and ˆyi the response value from the fit. We

expect that for a specific δφ the comparison will break down, because from that moment in the quasi-static deflation new contacts will appear.

5.3 The relation between gaps and forces

A very intuitive prediction is that the bigger the average distance between two particles is, the lower the contact force. This is also the result that is calculated in section ??. The intuitive argument for this is that a ball colliding between walls exerts a weaker force on that wall when we move the walls apart and keep the internal energy of the ball fixed. The average momentum exchange between the wall and the ball will be the same, but the collision frequency will be lower. We will verify this relation numerically.

5.4 The relation between gaps and packing fraction

We consider an isostatic packing, which volume is increased by a very small δφ. We then launch the hard sphere dynamics. In the following we will consider the coordinates |Ri to be the time-averaged mean positions of the particles. The corresponding gaps between the particles are the mean difference between the pairwise distances rij and the sum of the radii

of the pair of particles between which there is a gap dα≡ di+dj

2 (with dithe diameter of the

ithparticle), namely

hα = rα− dα. (34)

The gaps h must grow linearly with δφ = φc− φ. The argument is that they cannot grow

as h ∼ δφx _{with x < 1, because then for small increments δφ there will be contacts that}

open by a huge amount. If x were larger than 1, you would find that for small increments δφ some contacts did not open at all, which is impossible when we think about entropy. Also, when we create the hard disc packing, we increase the length of the box from L to (1 + γ)L. This means the volume changes to V0(1 + γ)3, with V0 the previous volume. When we Taylor

expand this, we get: V ≈ V0(1 + 3γ) so the volume and therefore the packing fraction φ

depends approximately linearly on h (∼ γ).

The fact that h ∼ δφ does not mean that the proportionality constant of all different gaps is the same, and in particular

hα = cαδφ , (35)

for every contact α. We’ve seen in section ?? and will show in the results that the effective force in the αth pair is just the inverse of the gap, fα= h−1α , up to a temperature dependent

factor that we take it be unity since it only sets the time scale. The forces due to collisions in the hard disc glass should converge to the contact forces in the jammed isostatic packing,

(18)

and therefore the equations hα = cαδφ (36) and hα∼ 1 fα (37) lead to hα ∼ 1 fα φ. (38)

Here we identify the pre-factors cα of the linear dependence of the αthgap with δφ as the

inverse of the force in the αth contact in the isostatic state.

Imagine now that we start with an isostatic packing and expand the system such that the packing fraction changes by a small increment δφ. We just concluded that not all gaps grow with the same rate, but instead with rates _f1

α; this means that the contacts carrying

the weakest forces will grow gaps with the largest rate.

The relations between h, f and δφ are checked numerically. The results are presented in chapter 7.

5.5 The agreement between the separations of weak forces

The information about influence on the displacement field when perturbing a contact should in our hypothesis somehow also be contained directly from the way the two particles that form the contact collide. In the hard disc dynamics the contact-forces are calculated with use of ~Fij =

P

c∆~pc

T and depend on two variables: the frequency of collision ν ≡ X T with

X the number of collisions and average momentum that is exchanged during a collision ¯p ≡

1 X

P

c∆~pc. Because of these two variables, we assume a separation of two kinds of weak

forces:

• Forces with small ν. Particles that rarely collide have a weak contact force between them.

• Forces with small ¯p. The configuration in figure 7 is expected to give rise to a weak force. The particles will move mostly in the horizontal direction and therefore the momentum in the direction of ~n is small, resulting in a small impulse.

In (5) there has already been made a separation between extended and localized forces, which is thoroughly described in section 4.3. To evaluate whether the assumption that these separations agree is correct, we calculate the frequency and the average momentum for all contacts in the hard disc glasses. We look at the distributions of the number of collisions and average momentum exchange in the system and compare those to the distributions shown in figure 5 and 6. Then we select out the contacts that are extended or localized and see if the distribution of their momentum exchange, for a specific contact, is different. It could for example be that the distribution of momentum exchange for the localized force is very broad and for the extended force is very narrow.

(19)

Figure 7: Even if the force f0 in the surrounding contacts is on the order of the mean force, the

force in the horizontal contact can be small if the angle φ is small, and displacements resulting from opening that contact will be of order δR ∼ sin φ.

The analysis of the extended and localized forces is done on systems that are larger than 64 particles. In order to be able to make a legitimate conclusion, we make sure that the distributions of these variables do not depend on system size.

6 The numerical approach

To test our hypotheses, we used numerical simulations of isostatic packings and hard disc glasses. In this section we first describe how these are created. Then we will cover the methods that are used to calculate the contact forces, average gaps, collision frequencies and impulses that we need.

6.1 Initial packing

The first stage in creating isostatic and later on hard disc glasses is the initial packing. In this experiment, the two dimensional packing consists of 64 discs in a box of size L × L. There are two types of discs with different size. At first all discs are given a place in the 2d-grid and the sizes are determined with small random fluctuations on top of the two given sizes to avoid hexagonal unions of similar particles(7). The discs are ’soft’, which means there can be overlap between them. We used the fast inertial relaxation engine (FIRE) to minimize the energy in the packing for a large number of time-steps (6). This creates a packing where the discs still have overlaps, shown in figure 8. The packing fraction here is larger than φc. The

size and position as well as the packing fraction of all the discs are saved in a file and are used to make an isostatic packing out of the initial packing.

(20)

Figure 8: Initial packing of 64 2D soft discs with overlaps, so φ > φc.

6.2 Isostatic packing

The initial packing is turned into an isostatic packing by letting the system evolve in time as long as the overlaps are still big. As soon as the pressure is below a small target value of 10−12, the isostatic packing is reached. The target-pressure should be small, but not so small that it is in the order of the error of the computer. In the end, all the contact forces between the particles are calculated. Together the form the contact network of the packing, shown in figure 10. For every packing, again a file with size and position of all particles is saved, but also a file with all the contact forces. The latter will be used to compare these to the forces in the hard disc glass, which is the first goal.

The generation of this packing is done in 128 bit precision, which supports data with 34 digits. This precision is very important in this specific step, because of the following reason. Every particle i has z neighbours j(i) that have a certain overlap with i. These overlaps ∆~xij= (~ri− ~rj) − (sizei+ sizej) are of order 10−10. This means there are only four digits

left of the fourteen we have in the regular double precision. In the next step, the overlaps are used to calculate the forces. The total force on a particle Fiis ideally zero, but on a computer

that cannot happen. Compared to the contact forces however, the total force is very weak (Fi

fij ≈ 10

−8

). We need more than four digits to calculate these values.

6.3 Hard disc glass

The next step is to see what happens if you give the particles a little more space and let them move around to collide as hard discs. As already mentioned, the hypothesis is that the contact forces calculated in this system will be exactly the same as in the isostatic packing, because the systems are both nearly at the jamming point and there is only one solution for the contact forces with a constant pressure, as described in section 4.1.

At first the isostatic packing has to be deflated, because the particles cannot have any overlap anymore. This is done just by increasing the linear size of the system according to:

(21)

where λ(h) is a function of the minimum gap h we want between the discs, given by:

λ(h) = h

2m, (40)

where m is the size of the biggest particle. Because h(order δφ) will be several orders of magnitude bigger than the overlap between the particles (which is of order 10−12), this overlap can be totally neglected. The minimal gap h just has to be large enough that the distance between the two larges particles (with radius m) is equal to h.

The variable h is directly related to the packing fraction φ. Since the overlaps are so small, the packing fraction in the isostatic packing is by very good approximation equal to φc. Expanding the system, phi will decrease, or δφ = φc− φ will increase.

After the deflation the packing is put into the hard disc evolution. This code takes the particles and gives them a certain random velocity, in such a way that the temperature is 1. Then, every particle moves in straight lines, until it has a perfect elastic collision with one of its neighbours. The loop runs for a certain amount of collisions (2 · 106 _{in our analysis).}

So after every event, the time till next collision is calculated (with use of the velocities and distances between particles) and the first particle that hits a neighbour is the next event. The interesting part in this is the exchange of momentum (impulse). The impulse is exactly the difference between the momentum of both particles, projected on the direction of the collision ˆ

n, which leads to:

∆~p2 =

∆pxdx2

dx2_{+ dy}2x +ˆ

∆pydy2

dx2_{+ dy}2y.ˆ (41)

For the complete calculation, see appendix A. The forces are calculated by dividing the total impulse by the time T the simulation runs, according to:

~ Fij =

P

c∆~pc

T . (42)

The simulation is running for a given number of collisions (2 ∗ 105). Every time two particles collide, the time until the next collision is calculated. Dividing the impulse by the total time gives us the right units of force.

7 Results

7.1 Comparison of contact-forces in isostatic and hard disc

glasses

In the isostatic- as well as the hard disc model, all contact forces are calculated during the simulation. We only look at the forces between the particles that were in contact in the isostatic case and compare them to the hard disc case. Because the hard discs are free to move, they will eventually collide with other particles than their original neighbours. We want this to happen as little as possible since it will mean the particles collide less with the neighbours that correspond to those in the isostatic packing, and the comparison will be less accurate.

(22)

When the forces are plotted against each other in figure 9 there is a linear correlation between them. This means the two ways of calculating the forces are actually equivalent up to a multiplicative constant. The equation of the linear line on a log scale is of the form fhard = a · fis0c . The slope of the linear line corresponds to c, so c = 1. The offset in the

logarithmic plot however corresponds to a. In the simulations we put the spring constant k to 1, which is the reason that a is not equal to 1. Because the overlaps ∆~x ∼ δφ, the proportionality constant a is proportional to the inverse of δφ. In figure 9, δφ = 10−9 and a = 1.044 ∗ 109_. _{This agreement becomes more intuitive when you think about a hard}

sphere bouncing between two walls. As you make the distance between the walls smaller, the frequency of collisions increases, and therefore the force.

Figure 9: The contact forces in the hard disc glass and the isostatic packing have a linear relation.

The contact network in the isostatic packing and hard disc glass are shown in figure 11. The thickness of the lines is proportional to the magnitude of the force. It is clear that the networks are the same, which they should be because the contact forces correspond perfectly.

(23)

Figure 10: The contact network in the isostatic packing, where the with of the connecting lines is proportional to the magnitude of the force.

Figure 11: The contact network of the same packing as figure 10 in the hard disc evolution. The networks are the same, up to a proportionality constant, except for the boundaries. The magni-tude of the force depends on frequency of collision and average momentum ex-change. Blue contacts have low fre-quency, green contacts have low aver-age impulse and purple contacts have both low frequency and low impulse.

This result gives us the information that our models are working in the right way and enables us to investigate the other physical properties described in section 5 and compare them between the both systems.

7.2 Variation of the packing fraction

For δφ increasing from 10−10tot 10−1with steps of one power of ten, the comparison between the isostatic and hard disc forces is made. The linear correlation is found in all of them, but the goodness of the fit decreasing drastically for larger δφ, as shown in figure 12. It is clear that the forces converge to the same value once the gap gets smaller. Because of the numerical approach it is not true that we can go arbitrarily small. For small gaps the overlaps in the isostatic packing become in the order of the gap (δφ ≈ 10−10). This means that due to numerical errors, there might be particles that are still overlapping when the hard disc dynamics is launched. This gives rise to errors in the hard disc dynamics.

The dynamics work very well for 4 orders of magnitude in δφ and completely breaks down once δφ ≥ 10−4. This invokes the idea that there is a certain length scale that is depending on δφ and sets a scale of the δφsat which new interactions will take place. We can define this

length scale as the size of the system when the packing fraction is δφs. In this system with

(24)

this for other system sizes as well to see how this length scale varies with N .

Figure 12: The goodness of the linear fit between contact forces in both systems, measured by the sum of squared due to error (SSE) is dependent of the increment δφ.

7.3 Inverse proportionality between gaps and forces

The time averaged gap between all particles (hhi) is calculated for 2 million collisions and in figure 13 the relation is shown between hhi and the contact force f . The slope of the graph is −1 which implies the relation f ∼ 1/hhi. This is the same result we derived using the partition function in section ?? and is therefore what we expected. The fit is very accurate for over 3 orders of magnitudes. For the smaller forces it fails, due to the low statistics because they collide less frequently.

Recall that this result is also very intuitive when thinking about a ball colliding between two walls. When the walls are closer to each other, the collision frequency and therefore the force will be higher.

(25)

Figure 13: The relation between the contact force f and the average gap hhi is inversely propor-tional, as indicated by the slope of -1 in the above figure.

7.4 Relation between gaps and packing fraction

Another property we looked at is how for a specific contact, the average gap between depends on φ. For 10 random contacts, the typical dependence is shown in figure 14. The relation is linear and the offset in the vertical axis is proportional dhhi_dδφ. A conclusion that is intuitively surprising is that the slope dhhi_dδφ ∼ 1

f, shown in figure 15. This means that when compressing

a hard sphere system quasi-statically towards the jammed state, the contacts carrying the weakest forces close the quickest. One would perhaps expect that if a contact closes quickly, the particles would collide ’harder’ (i.e. with a stronger force), but this is not the case. For the small δφ < 10−4, the results look very clean. When the system is expended more, new lines and therefore new contacts appear. This is exactly the same length as were the forces comparison (section 7.1) broke down.

(26)

Figure 14: The average gap hhi between 10 pairs of particles is plotted against the relative packing fraction δφ. The slope of the lines is proportional to the inverse of the force between the pair of particles. For δφ ≥ 10−4 the dynamics are not accurate anymore because new interactions appear. These new interactions can be seen in the figure in the contacts that start existing at a certain δφ. The graphs with the largest offset, the upper 2 lines, correspond to the steepest slope on a linear scale and therefore correspond to the lowest forces.

Figure 15: c = dh

dφ is plotted against 1

f for every contact force. The dependence is linear, which

(27)

The result that dhhi_dδφ ∼ 1

f gives a reason to take a look at the distribution of 1

f in central

limit theorem. We start with P (f ) = fθ _{and define y =} 1 f. Then: P (y) = df dyP (f ) = − 1 y2P (f ) = y −(θ+2) (43) Since 0 < θ < 1, this distribution decreases more slowly than y−3and therefore we can call this distribution broad (1). We are interested in the mean ’velocity’ of the closing/opening of contacts in the system, where we mean the variation of the gap with respect to φ, defined by: dhhi_dφ . The quantity we are looking for is therefore:

1 N N X i dhhii dφ = 1 N N X i 1 fi = 1 N N X i yi, (44)

where the sum runs over all contacts. We define Xn=PNi yi. Since P (y) is of the form

y−(1+µ)with µ = 1 + θ, we have the situation that 1 < µ < 2. Using central limit theorem, this implies that hyi is finite and ¯XN = hyiN . Therefore, the expected value for dhhi_dφ is just

hyi. The most interesting is the behaviour of Xn itself, which is the mean ’velocity’ of the

particles during quasi-static decompression. If µ has a certain value between zero and one, the sum would be completely defined by its largest value and therefore the mean velocity would be singular. However, since 1 < µ < 2, we have XN∼ N

1

µ _{so the mean velocity is:}

1 N N X i dhhii dφ = 1 NXN∼ N 1 µ−1_{= N}− θ 1+θ_. ₍₄₅₎

Therefore the velocity does not only depend on the smallest forces, so there are no singular behaviours.

7.5 Comparison of two separations of weak forces

The separation in localized and extended weak forces can be compared to the separation in forces with low frequency and low momentum exchange by looking at their distributions. The distributions of b and W are shown in figures 5 and 6.

7.5.1 Distribution of contact forces

It has been seen empirically but also in simulations that for small f the contact forces between particles in glasses are distributed like a power law: P (f ) ∼ fθ (4). The distribution is independent of the system size, shown in figure 16, which has also been shown in (4). There is a small error in the order of noise but the exponent seems not to change, θ = 0.22 for both systems. For that reason the packings of 64 discs are sufficient to do the analyses on the (weak) forces and to be able to make conclusions that also apply to packings with more particles. This is also the reason that we can compare the distributions we derived from our analysis to the ones already shown in figure 5 and 6.

(28)

Figure 16: The distribution of forces in packings with 64 and 2500 particles. Clearly the distri-bution does not depend on system size and therefore all analysis of quantities of interest can be done using the 64 disc packings.

7.5.2 Distributions of frequency and average impulse

The distributions of the number of collisions in a packing is shown in figure 17 and 18 shows the distribution of average impulse. The distribution of the number of collisions, which is the same as frequency upto a constant, looks the same as the distribution in figure 16. The distribution has a range of more than two orders of magnitude. However, the distribution of momentum exchange, figure 17, is very narrow. This is due to the fact that it is an average momentum exchange and therefore there can be collisions where the momentum exchange is low, but on average it will have the same value for all contacts. The conclusion is therefore that we cannot make a separation in the weak contact forces based on this momentum exchange. All the information about the magnitude of the force is contained in the frequency.

7.5.3 Distributions of momentum exchange for specific contact

Since the separation between weak forces with low average momentum exchange and forces with low frequency is not of interest and does not contain the information about localized and extended forces, we looked at the distribution of the momentum exchange of a specific contact during the entire simulation. We analyzed one localized contact with b 1 and one extended contact with W 1. The results are shown in figure 19. We wanted to check if there was a separation to make in these distributions. It could for example have been the case that localized have a much broader distribution of impulse than extended contacts. It turns out that the extended contact did have a much lower frequency and therefore collided less and the distribution is lower, but the shape of both graphs are the same. Therefore, we can also not make a separation looking at the distribution of momentum exchange.

(29)

Figure 17: Distribution of average mo-mentum exchange for all contacts in a system. The distribution is very nar-row compared to the force distribution in figure 16

Figure 18: Distribution of the number of collisions for all contacts in a system.

Figure 19: Distribution of momentum exchange for a single contact with b 1 and one with W 1. The distributions have the same shape.

8 Further research

There are numerous possibilities for research on these systems.

First of all, because the impulse is a vector, the direction of this vector constantly changes and therefore it can average itself. It would be interesting to investigate this and see whether

(30)

this could help in finding a separation between the two kinds of low forces. One could keep looking for other ways to separate between two kinds of low forces in the hard disc packing. The information about the localized/extended behaviour of the contacts should be contained somewhere in that packing, since the forces in the hard disc packing correspond so perfectly to the ones in the isostatic packing.

Other suggestions for further research would be to investigate the appearance of new interactions when quasi-statically decompressing the hard disc glass. It would be interesting to keep track of the number of new interactions during decompression and see in what rate they appear. Furthermore, the forces of the new interactions seem to behave differently from the original forces. Looking at figure 14 it is interesting to check whether or not the new interactions obey the equation < h >∼ _f1φ and what the magnitudes of the new forces are.

We found in our system that the comparison of the isostatic packing and hard disc glass breaks down for δφ ≥ 10−4 and mentioned that this could give rise to a certain length scale. A suggestion would be to check for this breakdown in systems that vary in number of particles N . That way one could investigate whether such a length scale, that is defined as the size of the system at which the comparison breaks down, is present.

A final suggestion would be to check whether the relations h ∼_f1 and h ∼ 1_fφ also break down for a certain δφ. There’s no comparison with the isostatic packing anymore, so the new interactions cannot be the reason for the relations to fail. If there is a certain length scale in this too, it would be interesting to investigate it.

9 Conclusion

The two models we used are the isostatic packing and the hard disc glass. We investigated a variety of properties of these systems. The first aim was to confirm that the contact forces in the isostatic packing and the hard disc glass are the same up to a multiplicative constant. This was expected because there should be only one solution to the set of equations that are based on mechanical equilibrium and constant pressure. The confirmation was found for hard disc glasses with a packing fraction close to the critical packing fraction (φc− φ ≤ 10−5).

When the packing was more dilute, new interaction appeared and the agreement breaks down. The theoretical predicted relation f ∼ 1

hhi is numerically verified, see figure 13. The relation

between the packing fraction and the gap is indeed linear and has the contact-depending scaling factor 1 f. The distribution of 1 f is flat tailed, P ( 1 f) = ( 1 f) −(2+θ)

, but since the exponent is between 2 and 3, instead of 1 and 2, the average _f1 is not dominated by the lowest force only. The separation between weak contact forces in localized and extended forces does not agree with the separation in low collision frequency and low momentum exchange. We still expect the information of the impact on the whole packing, which is used do distinct between localized and extended, to be maintained in the direct geometrical properties of the contact. It is however not contained in the frequency and average impulse and also not in the distribution of the non-averaged impulse for the specific contacts.

(31)

Appendices

A

Solving 2D collisions with conservation laws

When the initial velocities are known, the final velocities of the particles in a 2D elastic collision can be calculated using the conservation of energy and momentum. When thy collide, the particles will exchange a certain amount of momentum. Because we have particles with no rotation or friction, it is clear that there is only momentum exchanged along the normal ˆn (see figure 3). In the other direction, ˆn⊥, the momentum is perfectly conserved since no force is

acting in that direction. Therefore, we can simplify the 2D-collision to a 1D head-on collision, where we look at the momenta in the ˆn direction only. ˆn is the normal unit vector given by:

ˆ

n = dx

pdx2_{+ dy}2x +ˆ

dy

pdx2_{+ dy}2y.ˆ (46)

And we write ~pi,1= px,i,1x + pˆ y,i,1yˆ (47)

for the momenta. Let dx, dy be the distances between the particles in the respectively x-and y direction, pr,i/f,nbe the initial/final momentum in the r-direction of particle n.

The interesting component of the momentum is:

~ pi,1,k = (~pi,1· ˆn)ˆn (48) = ( px,i,1dx pdx2_{+ dy}2 + py,i,1dy pdx2_{+ dy}2)ˆn (49) = px,i,1dx 2 dx2_{+ dy}2ˆx. (50)

Part of this momentum will be exchanged, wheres the momentum in the perpendicular direction is constant during the collision for both particles. Because we now look at only one direction, we have a simple head-on collision, where we have two unknowns (~pf,1,kand ~pf,2,k)

and two equations to solve:

~ p2i,1,k+ ~p 2 i,2,k = ~p 2 f,1,k+ ~p 2 f,2,k, (51) ~ pi,1,k+ ~pi,2,k = ~pf,1,k+ ~pf,2,k. (52)

This set can easily be solved for ~pf,2,kand gives two possible solutions:

~

pf,2,k = ~pi,1,kor (53)

~

pf,2,k = ~pi,2,k. (54)

Equation 54 is the trivial solution for when the collision has not taken place. So it appears that the two particles interchange their momentum:

~

pf,2,k = ~pi,1,k (55)

and because of conservation of momentum: ~

(32)

Now we only want to know the exact momentum that is transported from one particle to the other. Define:

∆prˆr = ~pr,i,1− ~pr,i,2 (57)

(58) for particle n, and we finally have:

∆~p2 = ~pf,2,k− ~pi,2,k (59) = ∆pxdx 2 dx2_{+ dy}2x +ˆ ∆pydy2 dx2_{+ dy}2y.ˆ (60)

(33)

Acknowledgements

I would like to take this opportunity to thank my supervisor, Edan Lerner, for the support and the enthusiasm with which he speaks and thinks about the subject. Every meeting you made my motivation grow. What I will remember most of all is the spontaneous way he came up with new ideas and thinks them through immediately.

(34)

References

[1] G. Bouchaud. Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications.

[2] C. Brito and M. Wyart. Geometric interpretation of previtrification in hard sphere liquids. The Journal of Chemical Physics, 131(2):–, 2009.

[3] A. Donev. Jammed packings of hard particles. 2006.

[4] C. B. M. W. E. DeGiuli, E. Lerner. The force distribution affects vibrational properties in hard sphere glasses. Proc. Natl. Acad. Sci. USA,111, 17054, 2014.

[5] M. W. Edan Lerner, Gustavo D¨uring. Low-energy non-linear excitations in sphere pack-ings. Soft Matter, 9, 2013.

[6] F. G. M. M. Erik Bitzek, Pekka Koskinen and P. Gumbsch. structural relaxation made simple. 2006.

[7] D. G. W. M. Lerner, E. Simulations of driven overdamped frictionless hard spheres. Computer Physics Communications, 184, 2013.

Contact forces in nearly jammed hard disc glasses

Verslag van Bachelorproject Natuur- en Sterrenkunde

Omvang 15 EC, uitgevoerd tussen 28-04-2015 en 02-07-2015