Master project Theoretical Physics (60EC)
01-7-2015 — 24-6-2016
FNWI
Institute for Theoretical Physics Amsterdam (ITFA)
Jamming for potentials without sharp
cut-offs
Author:
Stefan Kooij
Student ID:
10243593
Supervisor:
Dr. E. Lerner
Second corrector:
Dr. G.J. Stephens
July 11, 2016
Contents
1 Preface 3
2 Introduction 3
3 Theory 4
3.1 The canonical soft and hard sphere models . . . 4
3.2 The exponential model . . . 5
3.3 The inverse power law model . . . 6
3.4 Elasticity . . . 7
3.4.1 Elastic stability . . . 8
3.4.2 The density of states (DOS) and the ratio ω∗/ω0 . . . 8
3.4.3 The bulk modulus B . . . 9
3.4.4 The shear modulus µ . . . 12
3.5 Predictions . . . 13
3.5.1 Scaling relation for µ/B for the canonical soft sphere model . . . 13
3.5.2 Scaling relation for µ/B for the ipl system . . . 16
4 Simulation details 17 5 Results 18 5.1 µ/B . . . 18
5.2 DOS . . . 21
5.3 ω∗/ω0 . . . 22
5.4 Growing length scale . . . 24
6 Discussion 25 7 Conclusion 27 8 Acknowledgements 27 A Alternative definition for the bulk modulus 28 B Microscopic expressions for µ and B 29 B.1 general derivatives . . . 29
B.2 shearmodulus . . . 31
B.3 bulk modulus . . . 31
Abstract
Some amorphous materials show an interesting property, their state of matter is something between a solid and a fluid, they are jammed. Two common models to examine the jamming transition are systems of soft or hard spheres. These models use interaction potentials that have a sharp cut-off, making it possible to describe many properties of the solid in terms of the contact number z. Potentials without such a sharp cut-off have not been considered in this context, but are of interest because descriptions in terms of the contact number now become meaningless, plus predictions of the scaling exponents can in principle be made by taking derivatives with respect to the control parameter. Two potentials are used; an inverse power law and an exponential potential. Both systems show the characteristic features of the jamming transition, among which are the vanishing of the shear to bulk modulus ratio and the emergence of an excess of low frequency modes. Using simple arguments, a mapping can be made from the control parameters of the soft sphere system to the control parameters of the exponential and ipl system, making it possible to predict the observed scaling laws. This further allows one to assign an effective contact number, for systems where the contact number is ill-defined.
1
Preface
This report is the end product of a master project which took place over a period of one academic year. It is the final step in acquiring the master of theoretical physics. The research was done under the supervision of Edan Lerner at the Institute of Theoretical Physics Amsterdam. Computations where mainly done on the SURFsara cluster computer Lisa.
2
Introduction
There exist many amorphous materials that share a curious property. Their state of matter is something between a solid and a fluid; they are jammed. Unlike a fluid, they can withstand a finite shear force before yielding, but when this yielding transition is reached, the substance can flow. Examples of such materials can be found in foams, emulsions, suspensions and granular media. The jamming transition, i.e. the transition from a fluid like state towards a more solid like state (or vice versa), can be controlled by various things. For example, for a foam it could be the gas fraction while for a suspension it could be the density of hard particles.
Although the jamming transition is a well known phenomenon, many things about it are not understood. This is mainly due to the fact that these substances are all disordered, yet have properties that depend on geometry. This complicates theoretical calculations greatly and causes effective medium theory (EMT) to be inadequate to explain all aspects of the jammed state [1]. Whenever theoretical calculations are too hard to do, one can try to make numerical simulation of a simplified model. Two very important models in the field of jamming are the packings of soft spheres and the packings of hard spheres. Both models, that will be introduced further later on, have a sharp cut-off in their interaction range. There is no clear reason to think that a model without such a sharp cut-off would not display the jamming transition, but such models are never used in this context. Furthermore, various scaling relations have been found in terms of the contact number z [1], but for a potential without a sharp cut-off the contact number loses its meaning. Hence, the aim of this research is to investigate the jamming transition for disordered solids with a potential that is continuous (in the sense that it always interacts with its neighbours) and if possible, relate the different scaling relations with each other to hopefully learn something new. Another motivation to use continuous potentials, is that derivatives with respect to the control parameter can now be taken, making it, at least in principle, possible to predict the scaling coefficients. This is not possible for the canonical systems, where not all derivatives are analytic. For example the second order derivative of the harmonic soft sphere potential would contain a jump at the cut-off radius. Any quantity that depends on this derivative would not change continuously by increasing or decreasing the pairwise distances (in which case your control parameter is the density). Derivatives with respect to density would therefore not be meaningful.
In this project, 2D solids are created of a binary system where the two types of particles interact through a pairwise potential. This is done by first letting the system equilibrate at some appropriate temperature, starting from a random starting configuration, after which the system is cooled down. To reach the zero temperature limit, a minimization algorithm is exploited to bring all net forces below some threshold value which is close enough to zero. This is then repeated for a number of system sizes. After that the control parameters is changed, minimizing the solid at each step. What that control parameter is, depends on the pairwise potential that is being used. In this project an exponential potential is investigated for which the control parameter is the density; the system is expected to unjam in the limit of zero density. Another potential that will be looked at, is an inverse power law (ipl) potential for which the control parameter is the exponent. This system is expected to unjam in the limit of an infinite exponent (the hard sphere limit).
This report is structured as follows. First there is a Theory section, wherein models are introduced plus calculations of mechanical properties and what values could be expected for them in the canonical case, followed by a section with some details of how the simulations are done, After that, the results are given with next the section Discussion, in which efforts are made to explain the results and also point out some issues in the calculations. Finally, the conclusion is given, with an overview of the most important results and some suggestions for future work.
3
Theory
Throughout this text to keep things as clear as possible, a few conventions are chosen. Particles indices are always denoted by Latin letters and always as superscripts, while components are always denoted by Greek letters and always as subscripts. Furthermore, vectors are signified by bold letters.
3.1
The canonical soft and hard sphere models
One of the most common used models in the field of jamming are the systems of soft or hard spheres. In the system of soft spheres, particles only exert a force on each other when their radii overlap and the magnitude of the force is proportional to this overlap. More formally, the potential energy function can be written as: φij = ( ij δij Ri+Rj α δij ≥ 0 0 δij < 0 (1) where δij = (Ri+Rj)−rij is the overlap (Figure 1). So Riis the size of particle i and rij= |rij|= |rj−ri|
is the distance between the two particles centers.
Figure 1: Two interacting particles
The potentials with different exponents α have various names, being the most prominent ones α = 2 and α = 5/2, which are called the harmonic and Hertzian potential respectively. To obtain a disordered solid one usually takes a 50:50 mixture of two sphere sizes where the ratio of the radii are 1:1.4. Systems of this type have a well defined jamming point, namely at the packing fraction where the particles are not overlapping but just touching. This critical point is denoted by φc and is identified with the
random close packing (RCP) fraction, although the RCP fraction is actually a dubious concept [2, 3]. Since the particles start to lose contact at the critical point, the pressure also vanishes due to all forces being zero.
A very important parameter in these systems is the contact number z, which is the average number of contacts a particle has. For these systems, z will attain the value ziso = 2d at the jamming point,
where d is the dimension and ziso is called the isostatic point [3]. To see why this is the case, one can
use a simple constrain counting argument. Since at the critical point all particles are just touching, the distances between particles are all equal to the sum of their radii. The N d positional degrees of freedom are thus subjected to N z/2 constraints (the number of pairs), which will only give solutions if z ≤ 2d. On the other hand, all particles are in equilibrium, so there are N d force balance equations with N z/2 variables, which gives z ≥ 2d. Therefore, z = 2d at the critical point.
In these systems there exist numerous scaling laws in terms of the distance towards jamming. This distance is usually expressed by the excess contact number ∆z = z − ziso or by ∆φ = φc− φ. One of
the most essential scaling laws are those for the shear modulus µ and the bulk modulus B with ∆z. In section 3.5.1 these scaling laws will be derived for the soft sphere system along with the important result that µ/B vanishes at the jamming point.
For the hard sphere model, particles only interact if the particles touch, but are not allowed to penetrate each other. The potential energy function is thus:
φij= (
∞ δij ≥ 0
This model is used to simulate granular media and is also used in the study of colloidal glasses [4]. The hard sphere model shows jamming at the same critical density as the soft sphere model, although there are subtleties [5]. The inverse power law potential, r−β, that will be used in this project will actually be the
same as the hard sphere model in the limit of the control parameter going to infinity, i.e. the exponent β.
3.2
The exponential model
As one of the new models without a sharp cut-off an exponential pairwise potential is used. The system consists of a 50:50 mixture of two types of particles of equal mass, say large (l) and small (s) particles. To make such a system, a length scale lij is introduced as well as a scale factor sij. Note that although the length scale gives rise to a certain “size” of the particles, it is not a real physical size as for the soft and hard sphere models. lij and sij just set a scale for the interaction. The pairwise potential that is
being used is then:
φ(rij) = φij= ( sije−rij/lij + (rij− r c) e −rc/lij lij − e−rc/l ij if rij ≤ r c 0 if rij > r c (3) A linear term and a constant term are added such that φij vanishes at the cutoff radius r
c as well as
its first derivative. The introduction of a cut-off radius is only to reduce calculation time and does not influence the physics. This assumption has been tested by removing the cut-off radius, that is; every particle interacts with every other particle in the system but not with itself (due to periodic boundary conditions), which gave no significant difference. The cut-off radius is set such that the second coordi-nation shell is within the interaction range. So unlike the particles in a soft sphere system, interactions are not only between nearest neighbours and more importantly, there exist no density where particles don’t interact, since the cut-off radius is always adjusted to the density.
The control parameter for the exponential system is the density ρ. As the pairwise distances increase due to a decrease in density, the interactions become softer just as for the soft sphere system. For the soft sphere system the solid will lose its stability if particles just start to lose contact, giving rise to a critical density. One could wonder if such a critical density will also be present in the exponential system. Since all derivatives change gradually with distance, it is not unlikely that the system will only unjam at zero density.
The sij and lij can take on 3 values, depending on whether the 2 particles are equal (2 possibilities,
large-large or small-small) or opposite (large-small). It turns out that the final structure of the solids is very sensitive to changes in these parameters (Figure 2).
Figure 2: Changing one of the parameters slightly causes a phase separation
How the different structures that appear depend on these parameter is an interesting question, but beyond the scope of this project. After many attempts a set of parameters was found that gave rise to
nice disordered solids (see Figure 3): lij= 1.0 i small, j small
1.2 i small, j large or vice versa 1.4 i large, j large (4) and sij = 1.0 i small, j small
1.64 i small, j large or vice versa 3.05 i large, j large
(5)
Figure 3: Example of a 2D solid created with an exponential potential, N = 1600
3.3
The inverse power law model
As a second model without a sharp cut-off an inverse power law potential (ipl) is used. Again with a system consisting of a 50:50 mixture of two types of particles of equal mass. The pairwise potential is adopted from [6] and is equal to:
φij = rij lij −β +P3 k=0c2k rij lij 2k if rij≤ r c 0 if rij> rc (6) With the coefficients given by
c2k= (−1)k+1 (6 − 2k)! ! (2k)! ! (β + 6)! ! (β − 2)! ! (β + 2k)· r −(β+2k) c . (7)
The extra terms in the potential energy function are present for the same reason as for the exponential potential; to make the potential vanish at the cut-off radius as well its first derivative. The length scales are now chosen to be:
lij = 1.0 i small, j small
1.18 i small, j large or vice versa 1.4 i large, j large
A pure inverse power law has the nice property that it is scale invariant. Still, one chooses a density such that the first coordination shell is within the cut-off radius. So just as for the exponential potential, neigbouring particles are always within the interaction range. The density is set to 0.86 and rc= 1.9.
Interesting about the inverse power law, is that in the limit of β → ∞ the potential becomes that of the hard sphere model (see Figure 4). The control parameter for this model is thus the exponent β.
Figure 4: Plot of r−β for various β. It is clear from this graph that the potential approaches the hard sphere model for increasing β
3.4
Elasticity
In this project we are always working in the zero temperature limit. This means that the solids are always sitting in some local minimum of the potential energy function (see Figure 5). When calculating the elastic properties of the solids one assumes that the deformation associated with that property will always increase the energy and that this deformation is elastic. That means that if the solid is let to relax again it will return to its original position (i.e. on the microscopic level). This is a seemingly trivial statement, but actually even small deformations can cause plastic flow [6]. It is actually this property that makes some solids so fluid like.
To describe the solids behaviour one can Taylor expand the energy around this local minimum: (9) U (r) = U r=r (0) +X l ∂U ∂rl r=r (0) ·rl− rl (0) +X k,l 1 2 rk− rk (0) ∂2U ∂rk∂rl r=r (0) rl− rl (0) + ... Where rlis the position of the lth particle and r = r
(0) are the positions of all the particles in the local
minimum (so it is a vector of vectors). All particles are in mechanical equilibrium, which means that the net force on every particle is zero. For that reason the second term in the expansion will cancel. The important part is therefore the third term in the expansion which can be rewritten as:
1 2∆ k· M kl r=r (0) · ∆l (10) Where Mkl := ∂2U ∂rk∂rl and ∆ k :=rk− rk (0) . (11)
Mkl is called the dynamical matrix and ∆k is the displacement vector from the local minimum of
the kth particle. The dynamical matrix is real and symmetric so always diagonalizable [7]. 3.4.1 Elastic stability
What makes a solid a solid is a question without a straightforward answer, especially in the context of jamming. One property it should at least have, is that it is in mechanical equilibrium, therefore ∂U/∂rl= 0 for all particles l.
A second property a solid should have, as discussed before, is that any small deformation should always increase its energy, so that once whatever is performing the deformation stops, all particles in the solid will return to its original position. Mathematically this means that all eigenvalues of the dynamical matrix are strictly positive. Because if there would exist negative or zero eigenvalues, the eigenvectors of the dynamical matrix M can be seen as displacement vectors that keep the solids energy equal or even lower it, which is precisely not what stability means. Still generally there exist symmetries in the system (e.g. translations) resulting in their associated Goldstone modes. So in equations where M−1 appears, one should actually use an eigendecomposition leaving out any zero modes.
3.4.2 The density of states (DOS) and the ratio ω∗/ω0
Important information about the mechanical properties of a solid can be extracted from the distribution of the eigenvalues of the dynamical matrix. More precisely, from the square root of the eigenvalues of M, which represents the vibrational frequencies of the solid. From Debye’s theory of solids, one expects that the density of states scales as ωd−1, where d is the dimension [8]. Systems near the jamming transition
show an excess of low frequency modes [9, 1], usually referred to as the boson peak (the ratio D(ω)/ωd−1
will show a bump at low ω). An example of this phenomenon can be seen in Figure 6, which shows the DOS for a system of soft particles as the critical packing fraction is approached.
Figure 6: The density of vibrational states for systems of 1024 spheres with a repulsive harmonic potential. From left to right with φ − φc= 10−1, 10−2, 10−3, 10−4and 10−8, where the curve with φ − φc= 10−1 is
labeled by a. ω∗ is defined by where D(ω) is half of the plateau value. The inset shows the scaling of ω∗ with δz; the excess coordination number. Adopted from [9].
The emergence of an excess of low frequency modes can clearly be seen as a plateau starts to form at small ω. Since this anomalous behavior is so characteristic for systems near the jamming transition, D(ω) will be determined for both systems investigated in this project. Different from the soft sphere system, the DOS needs to be rescaled as the control parameter is changed. The reason for this is that the local stiffness between particles, i.e. φrr, changes significantly for the different values of the control
parameter. Therefore, in the DOS the frequency is divided by the maximum frequency ω0.
In the system of soft spheres a characteristic frequency ω∗ is defined as where D(ω) is half the plateau value. This frequency scales as ω∗ ∼ δz and advocates for a diverging length scale l∗∼ 1/δz [1]. For the systems in this project an ω∗ is defined in a similar way as in [10] by
ω∗2= hδR|M∗|δRi, (12)
where |δRi is the displacement field of the particles coordinates and M∗ is the dynamical matrix without force components (i.e. φr). To see why this formula makes sense one can look at equation
10, where a displacement is related with a change of energy. Close to the jamming transition the force components will be negligibly small compared to the stiffness components (i.e. φrr). In calculating ω∗,
the force components are a cause for fluctuations in the data, since near the unjamming transition they will become negligibly small anyway, they will be omitted. Therefore M∗ is used instead of M. As displacement fields |δRi, the linear response from pushing apart two neighbouring particles is used, which is calculated by solving the equation
M∗|δRi = |di. (13)
So |di is the field where two particles are pushed apart. The obtained |δRi is then inserted in equation 12 to calculate ω∗. For both the exponential and ipl system this is done for many contacts and multiple
solids.
3.4.3 The bulk modulus B
The bulk modulus B is a measure of the systems response to decompression (or compression). It is an intensive quantity, usually defined as [8]:
B = −V ∂P ∂V
T
with V the volume and P the pressure. However, especially when doing numerics it is useful to adopt a different definition in which the derivative is taken with respect to a parameter η specifying a deformation that changes the volume (see Figure 7). Thus for all future calculations the bulk modulus B is defined as B := 1 V d2U dη2 η=0, (15)
where U is the total energy of the system.
Figure 7: Sketch of the deformation associated with the bulk modulus. The parameter η sets the amount of expansion or compression
Using
P := −∂U
∂V (16)
as the definition for the pressure, it is possible to find an expression which relates the two different definitions for B (see Appendix A).
The bulk modulus B can be expressed in terms of elastic constants [11]. These constants are just the coefficients of the multidimensional Taylor expansion of U in the variable , which is the strain tensor. That means the constants appearing in
U V = C αβ 1 αβ+ 1 2C αβνη 2 αβνη+ 1 6C αβνηκχ 3 αβνηκχ+ ... (17) Where: C1αβ= 1 V ∂U ∂αβ f, C αβνη 2 = 1 V ∂2U ∂αβ∂νη f, etc. (18)
The bold f indicates that the derivative is constrained. This is due to the fact that the derivatives are taken in the zero temperature limit, which means that the particles are in mechanical equilibrium. So the constrain on the derivative is: f = −∇U = 0. To make this calculation more explicit one can construct a simple strain tensor by defining the affine transformation of coordinates H(η) as:
H = 1 + η 0 0 0 1 + η 0 0 0 1 + η , (19)
where η is the dilatant strain. The effect of η can easily be seen by actually applying the transformation; r → H · r = (1 + η) · r. So this transformation just increases all distances by a factor (1 + η). Then the strain tensor becomes:
= 1 2 HTH − I=1 2 2η + η2 0 0 0 2η + η2 0 0 0 2η + η2 (20)
Now by Taylor expanding the energy in terms of the dilatant strain η, U = U η=0+ Uη,f η=0η + 1 2Uηη,f η=0η 2+1 6Uηηη,f η=0η 3+ ..., (21)
the elastic constants can be matched with the derivatives of U with respect to η. Where the short notation Uη,f= ∂U ∂η f and Uηη,f= ∂2U ∂η2 f,f etc. (22)
is being used. Writing out al terms and collecting all equal powers of η one finds: Uη,f η=0= V (Cxx+ Cyy+ Czz) Uηη,f
η=0= V (Cxx+ Cyy+ Czz+ Cxxxx+ Cyyyy+ Czzzz+ 2 (Cxxzz+ Cxxyy+ Cyyzz))
(23)
Note that the total derivative of U in the definition of B is the same as the constrained partial derivative of U . A completely equivalent definition of the bulk modulus is therefore:
B = 1 VUηη,f
η=0= (Cxx+ Cyy+ Czz+ Cxxxx+ Cyyyy+ Czzzz+ 2 (Cxxzz+ Cxxyy+ Cyyzz)) (24) This expression is all we need to find the values of B, by calculating all the elastic constants (how the Cijkletc. are calculated is explained in [11]). Still, from a theoretical point of view it is useful to rewrite
this expression in terms of ordinary partial derivatives. To do so, the constrained partial derivative with respect to η can be replaced by [11]
∂ ∂η f= ∂ ∂η + ∂Xl ∂η f ∂ ∂Xl, (25) where Xl
ν is the component ν of the lth particle of the additional non-affine displacement needed to keep
mechanical equilibrium. That is, the final change in coordinates with a deformation η keeping mechanical equilibrium is:
rl→ H(η) · rl+ Xl (26)
Thus the first order derivative becomes: Uη,f = ∂U ∂η + ∂Xl ∂η f ∂U ∂Xl = ∂U ∂η (27)
Since it doesn’t matter if one variates the potential with Xl or with rl, the second term cancels due to mechanical equilibrium. Taking another derivative to obtain the second order derivative:
Uηη,f= ∂2U ∂η2 + ∂Xl ∂η f ∂2U ∂Xl∂η. (28)
There is still a constrained partial derivative left, using the fact that the constrained derivative of the forces vanishes, one can get rid of this final term:
∂Fk ∂η f= ∂Fk ∂η + ∂Xl ∂η f ∂Fk ∂Xl = 0 → ∂Xl ∂η f ∂Fk ∂rl = ∂Xl ∂η f ∂2U ∂rl∂rk = ∂Xl ∂η fMkl = − ∂2U ∂η∂rk =⇒ ∂Xl ∂η f= − M −1 kl ∂2U ∂η∂rk (29)
So the final result becomes:
B = 1 V ∂2U ∂η2 − ∂2U ∂η∂rl M −1 kl ∂2U ∂η∂rk η=0 (30) If the particles interact through a pairwise potential φij this expression can be written in terms of the
positions of the particles within the solid. How this is done is explained in Appendix B. The result is: (31) B = 1 V X i<j φijqq(rij)2− X i6=l φilqqrilν × M−1 νρ kl × X i6=k φikqqrρik
3.4.4 The shear modulus µ
The calculation of the shear modulus is very similar as for the bulk modulus, only now we use a strain tensor corresponding to the transformation that characterises shearing (see Figure 8). The parameter which represents the amount of shear is given by γ and the definition of the shear modulus, µ, becomes:
µ := 1 VUγγ,f
γ=0 (32)
Figure 8: Sketch of the shear deformation and the parameter γ that sets the amount of shear The affine transformation for this type of shearing is:
H = 1 γ 0 0 1 0 0 0 1 . (33)
(Simply act with the unit vectors on H to see that this expression is correct) The strain tensor is therefore
= 1 2 0 γ 0 γ γ2 0 0 0 0 . (34)
Similar as for the bulk modulus the energy can be expanded in terms of γ: U = U γ=0 + Uγ,f γ=0 γ +1 2Uγγ,f γ=0 γ2+1 6Uγγγ,f γ=0 γ3+ ... (35) Then by also expanding the energy in terms of the components of the strain tensor as in equation (17), the derivative Uγγ,fcan be matched with the elastic constants. This yields
µ = Cyy+ Cxyxy (36)
for the shear modulus. The shear modulus can also be written in terms of ordinary partial derivatives. Since the derivation for the expression for the bulk modulus doesn’t depend on whether one uses η or γ, we can immediately write:
µ = 1 V ∂2U ∂γ2 − ∂2U ∂γ∂rl M −1 kl ∂2U ∂γ∂rk γ=0 (37) As for the bulk modulus this expression can be written in terms of the pairwise potential and the particles
coordinates, which is done in Appendix B. The final result is: µ = 1 V X i<j φijqq(r ij x)2(rijy)2 (rij)2 + X i<j φijq (r ij y)4 (rij)3 − X i6=l φilqqr il νrilxrily (ril)2 + X i6=l φilq r il xδνy+ ryilδνx ril − rνilrilxrily (ril)3 ! × M−1νρ kl × X i6=k φikqqr ik ρ r ik xr ik y (rik)2 + X i6=k φikq r ik xδρy+ riky δρx rik − rikρ rikxryik (rik)3 ! (38)
3.5
Predictions
3.5.1 Scaling relation for µ/B for the canonical soft sphere model
The scaling relations of µ and B with the distance to jamming ∆z can be derived for the soft sphere model by applying some basic linear algebra. As a start, we will go over some general facts of a matrix of the form SST. These facts can then be used to derive the scaling relations.
Any matrix S of size m × n can be factorized by using a singular value decomposition (SVD)[7]. This means it is possible to write S as:
S = U DV−1 (39)
Where D is a “diagonal” matrix and U and V orthogonal matrices. Then SST = U DV−1V DU−1= U D2U−1
STS = V DU−1U DV−1= V D2V−1 (40) D is a m × n matrix of the form
D = " M 0 0 0 # (41) Where M is a true diagonal matrix with r (= rank(S)) nonzero values. r cannot exceed the smallest of m and n.
The eigenvectors of D are of course the natural basis vectors ˆei of Rn. Calling the eigenvalues of D,
λi, the eigenvectors of SST and STS are U ˆei and V ˆei respectively, with eigenvalues λ2i. This means
that SST and STS share the same nonzero eigenvalues. Nonzero, because if λ
i= 0 it could also be that
the eigenvectors U ˆei or V ˆei are just the trivial zero vectors. But if λi 6= 0, then the eigenvectors are
nonzero and orthonormal.
To see that U ˆei is an eigenvector of SST, one can just multiply the two:
SSTU ˆei= U D2U−1U ˆei= U D2eˆi= λ2iU ˆei (42)
Applying V ˆei on S one gets:
SV ˆei= U DV−1V ˆei= λiU ˆei (43)
To say something about the number of zero modes of SST and STS one can use the rank theorem for both matrices:
rank SST+ dim Nul SST = m
rank STS + dim Nul STS = n (44) Also
rank SST = rank STS (45)
Therefore:
dim Nul SST − dim Nul STS = m − n (46)
Now getting back at the shear and bulk modulus. In bra-ket notation they are defined by: µ = 1 V ∂2U ∂γ2 − ∂2U ∂γ∂r M−1 ∂2U ∂γ∂r B = 1 V ∂2U ∂η2 − ∂2U ∂η∂r M−1 ∂2U ∂η∂r (47)
Yet again introducing a pairwise potential φ(rij) = φ(rα) = φα (α denotes a pair) one can write the shear and bulk modulus in terms of φ. Starting with the first order derivatives:
∂U ∂γ = ∂ ∂γ X α φα=X α ∂φα ∂rα ∂rα ∂γ ∂U ∂rk = X α ∂φα ∂rα ∂rα ∂rk (48)
Second order derivatives are then: ∂2U ∂γ2 = X α ∂2φα ∂rα∂rα ∂rα ∂γ 2 + * 0 X α ∂φα ∂rα ∂2rα ∂γ2 ∂2U ∂γ∂rk = X α ∂2φα ∂rα∂rα ∂rα ∂γ ∂rα ∂rk + * 0 X α ∂φα ∂rα ∂2rα ∂γ∂rk, (49)
∂2φα/∂rα∂rα is just the stiffness and can be set to 1 and near the jamming transition ∂φα/∂rα << ∂2φα/∂rα∂rα, so we can set the first order derivative to zero. Now finally, the linear algebra introduced
earlier can be used by defining S as
S = ∂r
α
∂rk, (50)
with on the rows the pairs (α) and in the columns the particles coordinates. Then the equations (49) transform to: ∂2U ∂γ2 = ∂r ∂γ I ∂r ∂γ ∂2U ∂γ∂r = ST ∂r ∂γ (51)
Here I is the identity matrix.
The dynamical matrix M can also be written in terms of S:
M = ∂ 2U ∂rl∂rk = ∂ ∂rl X α ∂φα ∂rα ∂rα ∂rk = X α > 1 ∂2φα ∂rα∂rα ∂rα ∂rl ∂rα ∂rk + * 0 X α ∂φα ∂rα ∂2rα ∂rl∂rk =⇒ M = STS (52) Thus: µ = 1 V ∂r ∂γ I ∂r ∂γ − ∂r ∂γ S(STS)−1ST ∂r ∂γ (53) And similar for B.
The dynamical matrix M has no zero eigenvalues (i.e. there are no floppy modes). In other words, dim Nul STS = 0, which implies that
dim Nul SST − 0 = m − n = N z
2 − N d = N
2∆z. (54)
As discussed before STS and SST share the same non-zero eigenvalues for which the following formulae hold:
STS|Ψωi = λ2ω|Ψωi
SST|φωi = λ2ω|φωi
S|Ψωi = λω|φωi,
(55) where |Ψωi and |φωi are both orthonormal sets.
One of the eigenvectors |φωi has an actual physical meaning, it is proportional to the contact forces.
Because ST|f i =X α ∂rα ∂rkf α=X i<j (δjk− δik)r ij rijf ij = X j in contact with k fjk= 0 =⇒ SST|f i = 0 (56)
So there exist an ω, with λω= 0, such that |φωi = |f i/phf|fi. What about all other zero modes? There
are of course more solutions to the mechanical equilibrium equation, but these contain force components that are negative, which could not be the real contact forces for this purely repulsive potential.
Since STS is a non-singular symmetric matrix one can use its spectral decomposition to write the inverse of STS and insert it in the formula for the shear modulus µ. Furthermore, by using the closure relation, the unit matrix in the representation of SST can be written as
I =X
ω
|φωihφω|. (57)
So that the shear modulus becomes µ = 1 V ∂r ∂γ X ω |φωihφω| − S X ω |ΨωihΨω| λ2 ω ! ST ! ∂r ∂γ = 1 V ∂r ∂γ X ω |φωihφω| − X λω>0 |φωihφω| ! ∂r ∂γ = 1 V ∂r ∂γ X λω=0 |φωihφω| ! ∂r ∂γ = 1 V X λω=0 φω ∂r ∂γ 2 , (58)
where the last sum is only over all zero eigenvalues of SST, for which we know there are N ∆z/2 many.
The ∂r∂γ components have no preferred sign, since there is no preferred direction in an amorphous solid (see Figure 9). So the components of
∂r ∂γ E
can essentially be seen as random numbers. The sum of N of such random numbers has an expectation value of √N (just like a random walk), a factor that is exactly cancelled by the normalization factor of 1/√N , meaning that the inner products are of order 1. Therefore µ = 1 V X ω=0 φω ∂r ∂γ 2 = ∆z 2 N V · O(1) ∼ ∆z (59)
Figure 9: The pairwise distances increase on average as frequent as they decrease for shearing For the bulk modulus we only have to change γ to η. We know that ∂r∂η = r, so now all but one inner product is again a sum of random numbers. The exception is the |φωi which is proportional to the actual
contact forces (so all elements of the sum are positive). This gives for the bulk modulus B = 1 V X ω=0 hφω|ri 2 = 1 V hf |ri2 hf |f i + 1 V X N ∆z/2−1 zero modes hφω|ri 2 . (60)
For the part involving the unphysical contact forces the previous argument applies. For the other term we can use that p ∼ f /Ld−1 and hf |ri = (p V d) to conclude that
Most important conclusion from these two scaling laws is that µ
B ∼
δz
constant + δz → 0 (62)
at the critical point. Qualitatively, this means that the solid starts to behave more and more as a fluid as the critical point is approached. Because characteristically, fluids show almost no resistance towards shearing, but are very hard to compress.
3.5.2 Scaling relation for µ/B for the ipl system
One of the motivations of this project is that the advantage of using a pairwise potential without a sharp cut-off is that, at least in principle, scaling coefficients can be calculated directly by taking derivatives with respect to the control parameter. One important dimensionless quantity is the ratio µ/B, which is a measure of how ’fluid’ the system is. For the ipl system one expect this ratio to scale like
µ B = β
θ, (63)
where θ is the critical exponent and β is the exponent of the inverse power law (i.e. the control parameter). Taking the total derivative with respect to β one gets
d dβ µ B = θβ( θ−1)= θ β µ B → θ = β B µ d dβ µ B . (64)
So by taking the derivative with respect to the control parameter one can obtain the coefficient θ. For the simulation we added terms to the power law potential to introduce a cut-off radius, for this derivation a cut-off radius is not necessary of course. Therefore the pairwise potential is just φij = rij/lij−β. For a power law potential the non-affine part of the bulk modulus is zero, because the potential is scale invariant. This can also be seen by invoking mechanical equilibrium again. For all particles the net force is zero, meaning Fnetl =X i6=l fil= X i6=l φilq r il ril = 0. (65)
In components this becomes Fnetl ν =X i6=l φilqrν r = X i6=l ril lil −β−1 rνil rillil = X i6=l ril lil −β−2 rνil (lil)2 = X i6=l φilqqrilν = 0 (66) So the non-affine part in the expression for B in equation (31) is identical zero. A further simplification of the formulas can be made by realising that φq/φqq→ 0 for β → ∞. By only considering µ/B for large
β, we can thus leave out all φq terms in the expression of the shear modulus. The ratio of the shear and
bulk modulus reduces to:
(67) µ B ≈ P i<jφ ij qq (rijx)2(rijy)2 (rij)2 − P i6=lφ il qq rilνrilxryil (ril)2 × M−1νρ kl × P i6=kφ ik qq rikρ rikxryik (rik)2 P i<jφ ij qq(rij)2
For the dynamical matrix M we have (also keeping only the second order derivative) M ≈X i6=l φilqqr il νrilρ (ril)2(δ lk− δik) =X i6=l β(β + 1) r il lil −β ril νrρil (ril)4(δ lk− δik) = β(β + 1)M 0 (68)
The total derivative with respect to β is a constrained partial derivative, because we keep mechanical equilibrium. So just as when we were dealing with derivatives with respect to the deformation parameter, we now have d dβ = ∂ ∂β f= ∂ ∂β + ∂Xl ∂β f ∂ ∂Xl = ∂ ∂β − M −1νρ kl ∂2U ∂β∂xk ν ∂ ∂xl ρ , (69) where ∂2U ∂β∂xk ν = ∂ ∂β X i6=k −β r ik lik −β rik ν (rik)2 ≈ X i6=k −β r ik lik −β rik ν (rik)2log rik lik . (70)
Since B is proportional to U the non-affine part of the derivative of the denominator is zero. Therefore: θ = β B µ d dβ µ B =β µ dµ dβ − β B ∂B ∂β (71)
The last part of this expression is then
(72) −β B ∂B ∂β = − β B β(β + 1)X i<j rij lij −β log r ij lij + (2β + 1)X i<j rij lij −β = −β1 B X i<j rij lij −β log r ij lij −(2β + 1) (β + 1)
The first term diverges for large β while the second term will converge to the value -2. Apparently the first term will cancel with a term from the part that is not considered yet. One can immediately see that the -2 will be canceled by a term that will appear in βµdµdβ where the derivative acts on the factors of β(β + 1). Then we are still left with:
θ = −β 1 B X i<j rij lij −β log r ij lij +β 2(β + 1) µ X i<j rij lij −β(rij x)2(rijy)2 (rij)4 log rij lij − 2 X i6=l ril lil −β ril νrilxrily (ril)4 log ril lil × M−10 νρ kl × X i6=k rik lik −β rik ρ r ik xr ik y (rik)4 − X i6=l ril lil −βril νr il xr il y (ril)4 × ∂M−10 ∂β νρ kl × X i6=k rik lik −βrik ρ rikxriky (rik)4 + β (β + 1)2 X i6=k rik lik −β rik ν (rik)2log rik lik × M−10 ρν kl × ∂µ ∂xl ρ (73) It not clear how a constant value can appear from this expression. On top of that, it is seemingly divergent due to the factor of β. So the objective of predicting the critical exponent for µ/B proves to be too difficult for now.
4
Simulation details
To make convincing statements about the systems with different interaction potentials, a lot of solids should be created, especially since fluctuations start to rise as the systems starts to unjam. Therefore, 1024 solids were created for each type of interaction. This was done by letting the system evolve with a certain temperature for a sufficient number of iterations, starting from a random configuration. After that, all heat is taken out of the system and an equilibrium is reached. For the minimization a conjugate gradient method was used.
For the exponential system a density ρ of 0.56 was chosen as a starting point for the calculations. Then using these configurations, the density was decreased with factors of 10−1/6, minimizing the solids at each step. This is a tedious process, but seemed necessary to make sure that the density dependency that is observed is not due to large structural changes (changes still happen of course). For densities lower than 2.6 · 10−2 quad precision was used.
For the ipl system, 1024 solids were made for an exponent of 12, where the configurations of these solids were then used to create the solids minimized for different exponents. The reason for this procedure, is that for high exponents the random starting configuration is at such a highly unfavorable position, that the time step has to be decreased dramatically to keep the system stable. It is therefore more efficient to use an already minimized solid with a lower exponent as starting position. Quad precision was used for all the calculations. For the very high exponents the extra terms added to the pure ipl were eventually set to zero, due to the coefficients c2k becoming negligibly small.
5
Results
5.1
µ/B
The ratio of µ and B for the exponential system shows a cross-over around ρ = 10−2, for high densities the ratio increases as µ/B ∼ ρ−1.2 with decreasing ρ, but for low densities the ratio vanishes as µ/B ∼ ρ0.2 (Figure 10).
Figure 10: The ratio of the shear and bulk modulus for the exponential system for N=1600 particles averaged over 1024 solids. The same was done for a system of N=3249 particles, only with less solids, but no significant difference was observed.
The shear and bulk modulus are calculated in terms of elastic constants as explained in the Theory section. These elastic constants can also be plotted separately as in Figure 11. The elastic constants themselves start to diverge even on a Log-Log scale, but do follow a power law when they are rescaled with the pressure. Luckily, when calculating µ/B the pressure terms cancel, making this rescaling useful after all. The points for CY Y and CXX are precisely 1 for all densities, which makes sense since
p = −(CXX+CY Y)/2 by definition. The other coefficients seem to follow a power law after the cross-over,
Figure 11: The elastic constants for the exponential system for N = 1600 particles, averaged over 1024 solids. The fit lines are shown in order. Due to the symmetry of the system the data points for CY Y and
CXX completely overlap, just as the CY Y Y Y and CXXXX.
For the ipl system the ratio µ/B seems to vanish as ∼ β−0.5 for β → ∞ (Figure 12). As expected no finite size effect is detected between a system of N = 1600 particles and N = 3249 particles. The elastic constants for the ipl system are shown in Figure 13.
Figure 12: The ratio µ/B for the ipl system for a system of N = 1600 particles, averaged over 1024 solids. No difference was observed with a system of N = 3249 particles. The fitted line is given as a possible asymptote.
Figure 13: The elastic constants for the ipl system of a system of N = 1600 particles. Fit functions appear in the same order as the fit lines.
5.2
DOS
For both the ipl and exponential system the density of states has been determined for multiple values of the control parameter to see if an excess of low frequency modes appears when the jamming transition is approached. This was done for two system sizes, N=1600 and N=3249, using the eigenvalues of 1024 solids, as expected no difference was observed between the two. To make a good comparison the DOS is rescaled with the maximum frequency ω0. For both system (Figure 14 and 15) the emergence of more low
frequency modes can clearly be seen as the control parameter is changed. For the ipl system exponents β = 8 and β = 512 were left out of the graphs to keep the plots from becoming to crowded, but they follow the trend that is observed.
For small ω/ω0for the exponential system there seems to appear certain peaks in the DOS. One might
think that these are due to random fluctuations or due to the binning analysis, but considering that these are the eigenvalues of 1024 solids, with N = 3249 particles this seems not to be the case. A quick attempt to match the peaks with phonon frequencies did not work, although this was problematic since the peaks are hard to locate.
Figure 14: Plots of the rescaled density of vibrational states, D(ω/ω0), for the ipl system with N=1600
particles, averaged over 1024 solids. The same thing was done for a system of N=3249 particles with no difference observed. The characteristic emergence of an excess of low frequency modes can clearly be seen.
Figure 15: Plots of the rescaled density of vibrational states, D(ω/ω0), for the exponential system with
N=1600 particles, averaged over 1024 solids. The same thing was done for a system of N=3249 particles with no difference observed. The characteristic emergence of an excess of low frequency modes can clearly be seen.
5.3
ω
∗/ω
0The characteristic frequency scale ω∗ can be calculated as described in section 3.4.2. The displacement
fields from pushing apart two neighbouring particles are used for many contacts and solids. Then from this data a distributions for the ω2∗ values can be constructed. The distributions contain large tails so the median instead of the mean is chosen as an indicator for ω2∗ . Since the average stiffnesses change by varying the control parameter for both systems the ω∗ should be rescaled. ω0 which is the rescale
frequency for the DOS serves that purpose. The data of the ipl fits very well with a power law, while the data for the exponential is very noisy. Attempts to reduce this noise failed.
Figure 16: Ratio of ω∗/ω0 for the ipl system with a system size of N = 3249. Note that ω∗/ω shows the
Figure 17: Ratio of ω∗/ω0 for the exponential system for a system of N=3249 particles. Note that the
scaling is similar to that of the ratio µ/B
5.4
Growing length scale
As the system starts to unjam, displacement fields due to a perturbation become more extended, indi-cating a growing length scale [1]. To visualise this a system of N = 10000 particles was created for the inverse power law potential. For a few exponents the resulting vector field of displacements for a dipole perturbation is plotted. Because only the decay of the displacement vectors away from the dipole is of interest, the displacement field are rescaled such that the maximum displacement for each field is the same. In Figure 18 one can clearly see that the field becomes more extended for larger exponents. Note that the dipole is in a slightly different place each time, the exponent can not be changed by large amount without changing the structure of the solid, so for each exponent a different contact pair is chosen.
Figure 18: Rescaled displacement fields caused by pushing apart two neighboring particles for different exponents β of the ipl system (N=10000). A growing length scale can clearly be seen.
6
Discussion
The aim of this project was to investigate the jamming transition for potentials without a sharp cut-off. Both the inverse power law and exponential system show the characteristic features of the jamming transition for all investigates quantities. Among these features are the emergence of an excess of low frequency modes, the vanishing of the shear to bulk modulus ratio µ/B and the vanishing frequency scale ω∗ (after correcting for the change in stiffnesses). Surprisingly, not only the inverse power law potential,
but also the exponential potential shows power law scaling relations.
The uncovered scaling laws can also be explained quantitatively, by relating the control parameters of the ipl and the exponential system to that of the soft sphere system. Since p ∼ f /hri with hri some average distance in the system and B ∼ hφrri one expect that
p B ∼
hφri
hφrri hri
This ratio will vanish as the jamming transition is approached and scales as ∼ ρ1/2 for the exponential system and ∼ β−1 for the ipl system (remember that hri ∼ ρ−1/2). For the canonical soft sphere system we know that p/B ∼ δz2[1]. Combining these relations makes it possible to map the control parameters
of the ipl and exponential system to that of the soft sphere system. So one obtains
δz ∼ β−1/2 and δz ∼ ρ1/4. (75) This implies that, although the contact number for potentials without a sharp cut-off has no meaning, one can still assign an effective (excess) contact number given a certain value of the control parameter. The shear to bulk modulus ratio scales as δz for the soft sphere system (section 3.5.1), just as the frequency scale ω∗[1], comparing that with equation 75 one sees that these results match perfectly with the ipl system and are quite close to that of the exponential system. For the exponential system one would expect µ/B ∼ ρ0.25 instead of µ/B ∼ ρ0.20.
For a good overview of all the scaling laws they are summarized in Table 1 (remember that these are measured values). System p/B µ/B ω∗ Soft spheres δz2 δz δz p/B µ/B ω∗/ω0 Exponential ρ0.5 ρ0.2 ρ0.24 ipl β−1 β−0.5 β−0.5
Table 1: Scaling relations for the different types of interaction potentials. The soft sphere scaling relations were taken from [1].
Additional to having a scaling exponent that is slightly different from what is expected, the exponential potential also shows a maximum around ρ = 10−2 for µ/B. This means that the scaling laws only hold
for very low densities. This cross-over behavior can also be seen in the DOS, were at first there is an decrease of low frequency modes, associated with an increase of µ/B (moving to lower densities), after which the exact opposite happens. What the origin of this maximum is, remains unclear. One could notice that in the equations for µ and B, terms of the form r2 e−r/l appear, which unlike the inverse
power law has a maximum as well. Still, in a fraction like µ/B one does not expect this to matter. The results for the ipl are quite robust in terms of the measurements of critical exponents as well as a theoretical explanation for them. Since ipl potentials are so widely used in model systems, it is useful to know how far away these potentials are from the unjammming transition. By comparing the known results for soft spheres [3] and the scaling relations for the ipl system more carefully, the relation between the excess contact number δz and β becomes
δz ≈ 0.5 · β−0.5. (76)
For the Lennard-Jones potential, which consist of an ipl with exponent 12 and an exponent 6, the unjamming transition is still quite far away with δz ≈ 1/10 assuming the greatest exponent of the two. The Lennard-Jones potential has many applications, for example in simulating liquids and also in the context of jamming were efforts are made to explain the boson peak [12]. Pure inverse power laws are also successfully used to predict structure factors for metals [13]. The exponents they find are of the order 10, making these metals just as the Lennard-Jones potential still far away from the unjamming transition.
Remarkably, the work of Bacher et al [14] on the quasi-universality of simple liquids seems a bit in contrast with the scaling laws found for the ipl system. According to their arguments, quasi-universality, which roughly means the existence of non-trivial similarities between systems with different potentials, is implied when both potentials belong to the so called EXP quasi-universality class. They first state that ipl’s with different exponents have similar structure and dynamics. Secondly, they argue that the quasi-universality between the hard sphere system and many other smooth potentials is explained by the fact that the ipl potential which belongs to this EXP quasi-universality class is equal to the hard sphere system in the limit of β → ∞, therefore the hard sphere potential belongs to this class. Yet, in
this project, although completely athermal, clear differences between ipl’s with different exponents are observed. It is not astonishing that people don’t find any significant changes when the exponent is only changed around the value β = 10. According to the established scaling laws, β should be varied a lot more to see any change at all, which apparently has not be done by anyone until now.
In this project only system sizes of N = 1600 and N = 3249 were used. Although no difference between the two system sizes is expected and also no difference was measured. It would be better to probe the values of µ/B over a longer range of system sizes. On top of that, the dimension could be increased as well, although in view of the results for soft sphere systems no difference should be expected. Unfortunately, for this project the computing capacities were too limited to do both of these things.
One might also comment on the way the solids in this project are created. Changing the density of the solids step by step is not very physical. Yet one would like to know that the scaling laws that are measured are not due to global structural changes, but such changes do occur when solids are made without using this procedure. Interesting, but very hard to answer, is to ask how these changes come about and how they could be quantified.
7
Conclusion
Although intuitively one expects that also for potentials without a sharp cut-off a jamming transition should be seen, it was still an open question that had to be answered. It turns out that at least for the inverse power law potential and the exponential potential the jamming transition can be seen. Predic-tions for the scaling laws can be made by relating the control parameters to the control parameters of the soft sphere system. Only the scaling laws of the exponential potential are somewhat different from what is expected, with e.g. and exponent of 0.20 instead of 0.25 for the scaling of µ/B. What the cause is for this difference is unexplained. Since the densities are very low and the shear modulus can be very sensitive to the value of the minimization criterion, it could be due to numerical error.
As a followup research one could introduce other potentials without a sharp cut-off to see whether they show the same behavior as the two potential tested in this project. Especially, to look at the ratio hφri /hφrri hri first, and maybe choose a potential that is expected to behave very different than the ipl
and exponential potential.
An interesting phenomenon arised in tuning the parameters of the exponential potential. The structure of the solids seem to be very sensitive to changes in of these parameters, only a small adjustment could produce a completely different solid. Equal ’sized’ particles seem to sometimes group together, producing a phase separation. This behavior is very interesting in itself, but could not be investigated in much detail in this project. As a possible future project, one could examine how these structures depend on the values of the parameters, the dimension of the solid, the density and the system size. Introducing a temperature in such systems one could also look at the cooling rate and how it effects the final structure of the solid.
8
Acknowledgements
I would first like to thank my master project supervisor Edan Lerner for providing this interesting subject and for all his help and support. I would also like to thank the people of the Lisa computer cluster, without their help I would never have been able to produce the necessary data.
A
Alternative definition for the bulk modulus
In many textbooks the following definition of the bulk modulus is used: B = −V ∂P
∂V
T
. (77)
For this project however, another more useful definition is implemented. Still one might want to compare the two definitions, which can be done by also writing (77) in terms of elastic constants. As a first step one can relate the volume V with the deformation parameter η. Since r → (1 + η)r it is easy to see that V (η) = (1 + η)dV0, (78)
where V0 is the initial volume and d the dimension. This equation can be inverted to
η = V
1/d− V1/d 0
V01/d
. (79)
So that the partial derivatives of η evaluated at V0 become:
∂η ∂V V 0 = V 1/d−1 V01/dd V 0 = 1 V0d ∂2η ∂V2 V 0 = 1 − d d2 V1/d−2 V01/d V 0 =1 − d d2 1 V2 0 Remembering that P = −∂U ∂V, (80)
B can now be expressed in terms of elastic constants by applying the chain rule: B = −V∂P ∂V = V ∂ ∂V ∂U ∂η ∂η ∂V = V ∂ 2U ∂η2 ∂η ∂V 2 +∂U ∂η ∂2η ∂V2 ! = 1
d2(Cxxxx+ Cyyyy+ Czzzz+ 2(Cxxyy+ Cxxzz+ Cyyzz) + (2 − d)(Cxx+ Cyy+ Czz)))
There is a clear difference between the two definitions. For example, for a 2 dimensional system the first order terms cancel, while for the other definition these terms are always there (in fact there is no dimensional dependency). What this apparent difference physically means is not immediately clear, but fortunately it is not relevant for this project.
B
Microscopic expressions for µ and B
B.1
general derivatives
Remember that for the bulk modulus we obtained B = 1 V ∂2U ∂η2 − ∂2U ∂η∂rl M −1 kl ∂2U ∂η∂rk η=0 (81) and similar for the shear modulus
µ = 1 V ∂2U ∂γ2 − ∂2U ∂γ∂rl M −1 kl ∂2U ∂γ∂rk γ=0 (82) These expressions are still very general, but in almost any practical case the energy function is given by a pairwise potential φij, making it possible to find expressions in terms of φij and the pairwise distances
only.
Equations (81) and (82) contain three distinct parts, these parts will be worked out here one at a time. The derivations will first be done for a general deformation, so a general strain tensor , then only after that the actual strain tensor will be filled in (so an belonging to the bulk or the shear modulus). When a deformation is applied the particles coordinates change, the amount of change can be calculated in terms of the strain tensor as follows:
(83) δrij = q (rij)T · HT · H · rij− rij = q (rij)T · (2 + I) · rij− rij = q 2 (rij)T · · rij+ (rij)2− rij = rij s 1 + 2 (r ij)T· · rij (rij)2 − r ij = rij ∞ X n=0 (−1)n(2n)! (1 − 2n)(n! )2(4n) 2 rijT · · rij (rij)2 !n − rij = ∞ X n=1 (−1)n(2n)! (1 − 2n)(n! )2(2n) ( rijT· · rij)n (rij)2n−1 = ∞ X n=1 (−1)n(2n)! (1 − 2n)(n! )2(2n) (rij ααβr ij β) n (rij)2n−1
Where the Taylor expansion of the square root is used.
The energy function U is not a function of the pairwise distances rij, but of the deformed pairwise
distances. To prevent any mistakes, one needs to keep a good distinction between the two, therefore qij
will serve as the deformed pairwise distances. So qij= rij+ ∞ X n=1 g(n)(r ij ααβrβij) n (rij)2n−1 , (84)
where g(n) are the coefficients of the expansion. The strain tensor could in principle be a function of multiple variables, one could e.g. shear and compress the solid at the same time. However, this complicates things greatly, therefore the deformations in this project are parametrized by only one variable, from now on denoted by γ. Then the following equalities should apply
αβ(γ) γ=0= 0 and αβ ∂n κν ∂γn γ=0= 0. (85)
Meaning simply that no deformation (γ = 0) gives no displacement and that the derivatives of should not diverge at γ = 0. These equalities allow us to through away a lot of terms that are evaluated at zero
strain. As a first step we will calculate ∂∂γ2U2, ∂U ∂γ = X i<j φijq ∂q ij ∂γ = X i<j φijq ∞ X n=1 g(n)n(r ij ααβr ij β) n−1 (rij)2n−1 rijκ ∂κχ ∂γ r ij χ ! (86) (87) ∂2U ∂γ2 = X i<j φijqq ∂q ij ∂γ 2 +X i<j φijq ∂ 2qij ∂γ2 =X i<j φijqq ∞ X n=1 g(n)n(r ij ααβr ij β) n−1 (rij)2n−1 rijκ∂κχ ∂γ r ij χ !2 +X i<j φijq ∞ X n=1 g(n)n(r ij ααβrβij) n−1 (rij)2n−1 rijκ∂ 2 κχ ∂γ2 r ij χ + ∞ X n=1 g(n)n(n − 1)(r ij ααβrβij) n−2 (rij)2n−1 rijκ∂κχ ∂γ r ij χ 2!
Evaluating this at γ = 0 almost every term of the expansion will die except a few. This result in the following : (88) ∂2U ∂γ2 γ=0= X i<j φij qq (rij)2 rijκ∂κχ ∂γ γ=0r ij χ 2 +X i<j φijq 1 (rij) rijκ∂ 2 κχ ∂γ2 r ij χ − 1 (rij)3 rijκ∂κχ ∂γ r ij χ 2! γ=0
Now we will calculate the mixed derivatives, where we will use the following equality: ∂rij ∂xl ν = ∂ ∂xl ν q (xjβ− xi β)(x j β− x i β) = rij ν rij(δ jl− δil). (89) So one gets: ∂2U ∂xl ν∂γ =X i<j φijqq r ij ν rij + ∞ X m=1 g(m)m(r ij ααβrijβ) m−1 (rij)2m−1 2νκr ij κ + ∞ X m=1 g(m)(1 − 2m)(r ij ααβrβij) m (rij)2m rij ν rij ! (δjl − δil) ∞ X n=1 g(n)n(r ij ααβrβij) n−1 (rij)2n−1 rijκχ∂κχ ∂γ r ij χ ! +X i<j φijq ∞ X n=1 g(n)n(n − 1)(r ij ααβrijβ) n−2 (rij)2n−1 rκij∂κχ ∂γ r ij χ ∂ ∂xl ν rτijτ λrijλ + ∞ X n=1 g(n)n(1 − 2n)(r ij ααβrβij) n−1 (rij)2n rijκ∂κχ ∂γ r ij χ rij ν rij(δ jl− δil) + ∞ X n=1 g(n)n(r ij ααβrijβ) n−1 (rij)2n−1 rκij∂κχ ∂γ r ij χ 2νκ ∂γr ij κ ! (90) Now taking the limit γ → 0 and collecting all non-zero components:
∂2U ∂xl ν∂γ =X i<j φijqq rij ν (rij)2 rijα∂αβ ∂γ r ij β (δjl−δil)+X i<j φijq 2r ij κ rij ∂νκ ∂γ − rνij (rij)3 rijα∂αβ ∂γ r ij β (δjl−δil) (91) For any tensor Aij which is anti symmetric in i and j the following holds:
X i<j (δjl− δil)Aij= 1 2 X i X j6=i (δjl− δil)Aij= 1 2 X i6=l Ail−1 2 X j6=l Alj=X i6=l Ail. (92)
Since rij is anti symmetric the final result becomes (93) ∂2U ∂xl ν∂γ =X i6=l φilqq ril ν (ril)2 rαil∂αβ ∂γ r il β +X i6=l φilq 2r il κ ril ∂νκ ∂γ − ril ν (ril)3 rilα∂αβ ∂γ r il β
With hindsight this way of calculating the mixed derivatives seems to be correct, having it done in another way gave wrong results; the mixed derivatives did not commute. Therefore, in Appendix C as a prove the derivative is done in the reversed order.
The final part is the dynamical matrix M. Starting with the first order derivative:
(94) ∂U ∂xl ν =X i<j φijq r ij ν rij + ∞ X n=1 g(n)n(r ij ααβr ij β) n−1 (rij)2n−1 2νκr ij κ + ∞ X n=1 g(n)(r ij ααβr ij β) n (rij)2n rijν rij ! (δjl− δil).
It is easy to see that for this derivative no contribution of the expansion will survive so we are left with: (95) ∂2U ∂xl ν∂xkρ γ =0= X i<j φijqqr ij νrρij (rij)2(δ jl− δil)(δjk− δik) +X i<j φijq δρν rij − rij νrijρ (rij)3 ! (δjl− δil)(δjk− δik)
So that the final result becomes: ∂2U ∂xl ν∂xkρ γ=0= X i6=l φilqqr il νrilρ (ril)2(δ lk − δik) +X i6=l φilq δρν ril − ril νrρil (ril)3 ! (δlk− δik) (96)
B.2
shearmodulus
Now we have found the general expression we can fill in the strain tensor belonging to an actual defor-mation. Remember that the strain tensor belonging to shearing is:
= 1 2 0 γ 0 γ γ2 0 0 0 0 . (97) Therefore: rκij ∂κχ ∂γ r ij χ γ=0 = rijxr ij y rijκ∂ 2 κχ ∂γ2 r ij χ γ=0 = (rijy)2 2r il κ ril ∂νκ ∂γ γ=0 =r il xδνy+ rilyδνx ril (98)
Then the final expression for the shear modulus becomes using (rij)2= (rijx)2+ (ryij)2:
µ = 1 V X i<j φijqq(r ij x)2(rijy)2 (rij)2 + X i<j φijq (r ij y)4 (rij)3 − X i6=l φilqqr il νrilxrily (ril)2 + X i6=l φilq r il xδνy+ ryilδνx ril − ril νrilxrily (ril)3 ! × M−1νρ kl × X i6=k φikqqr ik ρ r ik xr ik y (rik)2 + X i6=k φikq r ik xδρy+ riky δρx rik − rikρ rikxryik (rik)3 ! (99)
B.3
bulk modulus
For the bulk modulus the strain tensor is now parametrized by η instead of γ and is = 1 2 HTH − I=1 2 2η + η2 0 0 0 2η + η2 0 0 0 2η + η2 (100)
Then: rκij∂κχ ∂η r ij χ η=0 = (rij)2 rκij∂ 2 κχ ∂η2 r ij χ η=0 = (rij)2 2r il κ ril ∂νκ ∂η η=0 = r il ν ril (101)
Using these equations in the general form the bulk modulus becomes:
(102) B = 1 V X i<j φijqq(rij)2− X i6=l φilqqrilν × M −1νρ kl × X i6=k φikqqrρik
C
Mixed derivatives of U
To show that the derivative of U with respect to the coordinates and the control parameter commute, we will now take the derivative in the reversed order. The first order derivative becomes
(103) ∂U ∂xl ν =X i<j φijq rijν rij + ∞ X n=1 g(n)n(r ij ααβr ij β) n−1 (rij)2n−1 2νκr ij κ + ∞ X n=1 g(n)(r ij ααβr ij β) n (rij)2n rνij rij ! (δjl− δil).
The second order term is then ∂2U ∂γ∂xl ν =X i<j φijqq ∞ X m=1 g(m)m(r ij ααβrijβ) m−1 (rij)2m−1 rκij∂κχ ∂γ r ij χ rij ν rij + ∞ X n=1 g(n)n(r ij ααβrβij) n−1 (rij)2n−1 2νκr ij κ + ∞ X n=1 g(n)(r ij ααβrβij) n (rij)2n rij ν rij ! (δjl− δil) ! +X i<j φijq ∞ X n=1 g(n)n(r ij ααβrijβ) n−1 (rij)2n−1 2 νκ ∂γr ij κ + ∞ X n=1 g(n)n(n − 1)(r ij ααβrijβ)n−2 (rij)2n−1 2νκr ij κ ∂ ∂xl ν rijττ λr ij λ + ∞ X n=1 g(n)n(n − 1)(r ij ααβrijβ)n−2 (rij)2n rij ν rij rκij∂κχ ∂γ r ij χ ! (δjl− δil) (104) Taking again the limit γ → 0 one can see that the result is the same as (91).
References
[1] M. van Hecke, Journal of Physics: Condensed Matter 22 (2010) 033101. [2] F. Bolton and D. Weaire, Phys. Rev. Lett. 65 (1990) 3449.
[3] C.S. O’Hern et al., Phys. Rev. E 68 (2003) 011306.
[4] E. DeGiuli et al., Proceedings of the National Academy of Sciences 111 (2014) 17054. [5] V. Baranau and U. Tallarek, Soft Matter 10 (2014) 7838.
[6] H.G.E. Hentschel et al., Phys. Rev. E 83 (2011) 061101.
[7] D.C. Lay, Linear Algebra and Its Applications, 4 ed. (Pearson Education, Inc., 2012). [8] D.V. Schroeder, Thermal Physics (Addison Wesley Longman, 2000).
[9] M. Wyart, S.R. Nagel and T.A. Witten, EPL (Europhysics Letters) 72 (2005) 486. [10] L. Yan, E. DeGiuli and M. Wyart, EPL (Europhysics Letters) 114 (2016) 26003. [11] S. Karmakar, E. Lerner and I. Procaccia, Phys. Rev. E 82 (2010) 026105. [12] N. Xu et al., Phys. Rev. Lett. 98 (2007) 175502.
[13] F. Hummel et al., Phys. Rev. B 92 (2015) 174116.
